[R-sig-eco] logistic and exponential model fitting

Lara R. Appleby 04 Lara.R.Appleby.04 at Alum.Dartmouth.ORG
Sat Jan 28 21:42:56 CET 2012


I reared approximately 300 new ant queens until their pupae were about to hatch.  
I then collected all pupae, shuffled them and dealt them out to the queens in  
groups of 1, 2, 5, 7, or 9. These are the five treatments in my experiment. For  
example, approximately 70 queens received 2 pupae. I then counted the adult workers  
present in each colony every week for the next 21 weeks. So basically, my response  
variable is colony size.

I want to fit a function to each treatment's data. In other words, I want to model  
the growth rate of colonies in each treatment separately. I have started by fitting  
all the data to a logistic function and an exponential function and comparing  
the fits using AIC.

The exponential function fits just fine; when I start at r = 0.1, it yields a  
best fit r of 0.1339. I can even calculate the 95% confidnece intervals of r as  
[0.1338,0.1340]. It is also pretty simple to find the best fit r in the logistc  
function. If I start fitting at r=0.1 and K=10000, I get a best fit r of 0.1339.  
I can even calculate the 95% confidence intervals for r as [0.1339,0.1341]. The  
best fit r value seems to be pretty stable to changes in the starting r value.  


Problem 1: There is something strange with the logistic function's K.
The best fit K, seems to depend on the starting value. When the starting value  
is 10000, best fit K is 9999.75. In fact, the best fit K is just slightly below  
the starting value no matter how high of a K I start at. And I cannot calculate  
the confidence intervals of K. Instead, R produces this warning message: "In .local(fitted,  
...) : stepsize effectively zero/flat profile (K)".

Leaving the problematic logistic function alone for a minute, let's try using  
the exponential function to answer my remaining biological questions. One thing  
I want to know is whether the parameter (r) values are different for exponential  
functions fit to different subsets of the data. There were five treatments in  
my experiment. Colonies in different treatments differed in the number of pupae  
they received at the beginning of the experiment. I want to fit an exponential  
curve to each data set separately and see if the confidence intervals of the r  
values for the functions overlap.

Problem 2: The exponential function's r behaves strangely when different functions  
are fit to different treatments.
If I use the same intial value (N0) for all the treatments, calculations are smooth  
only for treatments 2 and 9. For the other treatments, I get tons of error messages.  
For example, I can find the best fit r for treatments 1,5, and 7, but when I try  
to find the confidence intervals for r for these treatments, R produces this message:  


"Profiling...
Profiling has found a better solution, so original fit had not converged:
(new deviance=3.232e+05, old deviance=4.009e+05, diff=-7.767e+04)
Returning better fit ...
returning better fit

Call:
mle2(minuslogl = function (r)
{
     count <- exponential(f$week, r)
     -sum(dnbinom(f$workers, size = dim(f)[1], prob = exp(count)/(1 +
         exp(count)), log = TRUE))
}, start = list(r = 0.0928385773286232), fixed = list(r = 0.0927786619191297),  

     skip.hessian = TRUE, lower = -Inf, upper = Inf, control = list())

Coefficients:
        r.r
0.09277866

Log-likelihood: -161616.4
Error in if (object at details$convergence > 0) cat("\nWarning: optimization did  
not converge (code ",  :
   argument is of length zero"

Anyways, it may be more appropriate to change the initial value with treatment.  
For example, the 2 pupae treatment should be fit to a function with an initial  
value of 2, not 1. If I change the initial values for each treatment, some of  
the best fit r's are negative (e.g. -0.056 for treatment 9) and some are positive  
(e.g. 0.044 for treatment 2) and the confidence intervals for r for treatment  
1 cannot be calculated. Instead R produces this error:

"Profiling...
Error in spline(x = pv, y = pro[, 1]) :
   NA/NaN/Inf in foreign function call (arg 4)"



More information about the R-sig-ecology mailing list