HUMPFIT - Least squares estimation of no-interaction model for biomass species richness data

I suggested in Journal of Ecology 84, 293-295 (1996) that a simple `no-interaction' model could be used to produce a humped-back species richness - biomass curve. I presented a set of GLIM macros that can be used for maximum likelihood fitting of the model, together with Fisher's link function to translate number of units ('plants') into number of species. These macros may be difficult to use, especially without access to GLIM. Therefore I wrote a small program humpfit.exe to perform least squares estimation of the model.

Program humpfit fits the following model to species richness (spno) biomass (mass) data:

	x = mass/hump		if (mass <= hump)	
	x = (hump/mass)**2	otherwise	
	fitted = alpha*log(1+b*x/alpha)

The model has the following parameters

The program uses Levenberg-Marquardt method as implemented by Numerical Recipes. for estimating the parameters. The program can be invoked with command

	humpfit < infile >outfile		
	humpfit < infile 

Where infile is the name of the input file. If outfile is not given, output will be sent directly to the screen. Input file must have only pairs of values of biomass species-number and no other text or numbers. Admittedly this is clumsy, but hopefully you can live with that (or write yourself better input routines).

Iteration process and its parametrization

The program uses iterative process to find the least squares estimates of the parameters. The program estimates log(alpha) instead of alpha, since this gives much more stable iterations, better convergence, and prevents alpha getting negative values which is should not have in theory. The program starts with the following default values:

The final results may differ when these initial values are changed. The starting values can be changed with giving new values on the command line in the order hump, alpha, b. For instance: humpfit < infile.dat 500 2 14 (the order of '< infile' and starting values can be changed.) If you want to fix a parameter value, give it as a negative number. Corresponding positive number will be used, but the parameter is not changed in iterations (and so its s.e. is 0). Convergence cannot be guaranteed in non-linear fitting, and so it may be advisable to vary input values to see if the results are stable.

Interpretation of the output

The standard errors of the parameter values are computed with the Numerical Recipes funtion mrqmin, together with the parameters. After parameters, the programs gives basic goodness of fit statistics. R-squared is estimated from the residual sum of squares and sum of squares assuming constant average.

Distributing, copying and developing of the program

The program uses Numerical Recipes routines for Levenberg-Marquardt method. This zip file comes with the program code written by me (humpfit3.c), but the least squares routines cannot be distributed since this would be against Numerical Recipes copyright. However, I have interpreted their web pages to give permission to distribute executable (I rely on their correcting me if I'm wrong). The C program written by me is in public domain and it can be used to build improved versions of the routines. In that case I would appreciate receiving these improved versions and a notification if there are plans to develop the program.

Does the program work?

In some test data sets the program seems to give fairly stable and reliable results. Often the results are very close to the Poisson likelihood model fitted with GLIM. However, convergence cannot be guaranteed in non-linear fitting. If the guesses for initial values are bad, the program may converge to a bad solution (even with negative R²) or end up in an error due to numerical problems (SIGFPE = floating point error). In these cases you can try with other starting values. Multiple starting values can be used to check the stability and optimality of the results as well. The program is not tested thoroughly, and I cannot guarantee that it works. The program is distributed in hope it is useful, but no warranty is given. I'm interested in all reports on the performance of the program: bad and good.

Some examples on fitted models produces with humpfit can be found through my web pages.


24/10/97 Jari Oksanen