2
$\begingroup$

I've tough luck with the use of nls() in R for the following model

$$N_e = N_o\{1-exp[\frac{(d+bN_o)(T_h N_e - T)}{(1+c N_o)}]\}$$

where $b>0$, $c\geq 0$, $T_h>0$, and $T=72$.

This code

T <- 72
NLS.Fit3 <- nls(Ne~No*(1-exp((d+b*No)*(Th*Ne-T)/(1+c*No))), data = Data,
            start = list(d = 0.01, b = 0.01, Th = 0.01, c = 0.01),
            control = nls.control(maxiter=50, tol=1e-05, minFactor=1/1024))

gives the following error message:

Error in nls(Ne ~ No * (1 - exp((d + b * No) * (Th * Ne - T)/(1 + c * : singular gradient

And the following

NLS.Fit31 <- nls(Ne~No*(1-exp((d+b*No)*(Th*Ne-T)/(1+c*No))), data = Data,
            start = list(d = 0.01, b = 0.01, Th = 0.01, c = 0.01),
            control = nls.control(maxiter=50, tol=1e-05, minFactor=1/1024),
            algorithm = "port", lower=c(0, 0, 0, 0))
summary(NLS.Fit31)

code converges but provides the wrong results (drastically different from PROC NLIN)

Formula: Ne ~ No * (1 - exp((d + b * No) * (Th * Ne - T)/(1 + c * No)))

Parameters:
Estimate Std. Error t value Pr(>|t|)  
d  0.008325   0.003488   2.387   0.0192 *
b  0.000000   0.000064   0.000   1.0000  
Th 0.000000   0.614220   0.000   1.0000  
c  0.020670   0.034439   0.600   0.5500  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.631 on 85 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

I'd prefer to do this in R rather than SAS and how the constrains can be placed on only few paramters. Any help in this regard will be highly appreciated. Thanks

Data is here:

No  Ne
5   0
5   1
5   1
5   2
5   2
5   2
5   2
5   3
7   0
7   0
7   1
7   1
7   2
7   2
7   2
7   3
10  1
10  1
10  2
10  2
10  3
10  3
10  3
10  4
10  7
15  1
15  1
15  3
15  3
15  4
15  5
15  5
15  5
20  3
20  4
20  7
20  7
20  8
20  8
20  9
20  11
25  4
25  5
25  6
25  7
25  9
25  9
25  13
25  14
30  5
30  8
30  10
30  11
30  11
30  12
30  14
30  20
45  4
45  7
45  8
45  10
45  11
45  14
45  15
45  19
60  9
60  14
60  14
60  16
60  18
60  21
60  24
60  26
80  7
80  11
80  12
80  15
80  17
80  12
80  21
80  23
100 7
100 8
100 10
100 11
100 15
100 24
100 26
100 33
$\endgroup$
  • $\begingroup$ The notation is strange: the only independent variable is $N_o$ and you are trying to model $N_e$ in terms of itself! What really are you trying to do here? $\endgroup$ – whuber Apr 23 '11 at 4:05
  • $\begingroup$ @Whuber: I'd appreciate if you have a look on this link esapubs.org/Archive/ecol/E088/187/appendix-B.htm. (B.1). This notation also confused me first. $\endgroup$ – MYaseen208 Apr 23 '11 at 4:16
  • $\begingroup$ I don't know if this is the case, but everytime I've used nls and I've got those error it was due to problems in the start parameters... I'll leave a proper answer to the real statisticians ;) $\endgroup$ – nico Apr 23 '11 at 6:36
  • $\begingroup$ This is a numerics problem -- in both cases you have no true convergence. Can you post the data you're fitting to, and/or at least SAS output? And if this $N_e$ is the same on the both sides, you should solve the eq for it first. $\endgroup$ – user88 Apr 23 '11 at 8:04
  • 1
    $\begingroup$ Here is the answer from @Ben Bolker. stackoverflow.com/questions/6701051/… $\endgroup$ – MYaseen208 Sep 22 '11 at 8:09
2
$\begingroup$

I now think you have a problem with data. Why?

First of all, let's get rid of $e^x$, solving it out and taking $\log$ of both sides

$$\log\left(\frac{N_0-N_e}{N_0}\right)=\frac{(d+bN_0)(T_h N_e -T)}{1+cN_0}.$$

Now, this log should be a linear function of $N_e$ for fixed $N_0$. As I understand, $T-T_h N_e $ is a time available to the predator to predate and this $-\frac{d+bN_0}{1+cN_0}$ (let's call it $\lambda$) is the frequency of attacks, which is told to be somehow related to the number of prey. Putting this $\lambda$ into the model, we have

$$\log\left(\frac{N_0-N_e}{N_0}\right)=\lambda T- \lambda T_h N_e,$$

linear function of $N_e$. Is this in data? It seems so: enter image description here

Thus, we fit lines, get line parameters

       lamT       -lamTh   No
1  0.05753546 -0.29739723   5
2  0.01408135 -0.18090769   7
3  0.13772614 -0.18005256  10
4  0.02065162 -0.08429763  15
5  0.09416886 -0.07751886  20
6  0.11148562 -0.06477383  25
7  0.19165200 -0.06134858  30
8  0.04328295 -0.03025566  45
9  0.06706140 -0.02399182  60
10 0.02236404 -0.01553757  80
11 0.01947130 -0.01247877 100

and try to calculate the $\lambda$ and $T_h$; but something is wrong at this point:

enter image description here

$-\lambda T_h$ is a clearly $(\alpha N_0+\beta)^{-1}$ (what corresponds to $b=0$), but $T/T_h$ is a total mess; thus my guess is that $T$ was not equal to 72 for all samples and this is the main origin of your problems.

$\endgroup$
  • 1
    $\begingroup$ @mbq: The LHS must be $$\log\left(\frac{N_0-N_e}{N_0}\right)$$ $\endgroup$ – MYaseen208 Apr 24 '11 at 6:29
  • $\begingroup$ Thanks all, good point with $N_0$. My answer still doesn't solve the problem, but I hope it will be more useful now. $\endgroup$ – user88 Apr 24 '11 at 18:06
  • $\begingroup$ @mbq: The following SAS code was used by the author for this dataset: PROC NLIN DATA=NOTONEC2; PARMS BHAT= 0.001 0.01 0.1 CHAT= 0.001 0.01 0.1 DHAT= 0 THHAT=3.0; BOUNDS BHAT>0,CHAT>=0, THHAT>0; T=72; X=NE; A=(DHAT+BHATN0)/(1+CHATN0); C1=EXP(-AT); C2=ATHHAT; H=N0*C1*EXP(C2*X)+X-N0; ITER=0; DO WHILE(ABS(H)>0.0001 AND ITER<50); X=X-H/(N0*C1*C2*EXP(C2*X)+1); H=N0*C1*EXP(C2*X)+X-N0; ITER=ITER+1; END; MODEL NE=X; Run; $\endgroup$ – MYaseen208 Apr 24 '11 at 19:44
  • $\begingroup$ On the revised graph, these are not straight lines, but curves of the form $y=\log(1-\frac{x}{k})$. They may look fairly straight, but some are obviously concave. $\endgroup$ – Henry Apr 24 '11 at 20:12
  • $\begingroup$ @MYaseen208 I have no access to SAS so I can't check, however I'm pretty sure that this model does not fit this data and thus the SAS solution is either unstable or degenerated (as your R one). $\endgroup$ – user88 Apr 24 '11 at 20:18
1
$\begingroup$

I think this is is not a good question: first having $N_e$ on both sides of the model, and second I suspect it is ill conditioned.

Here are my explorations. Let's use R to try to minimise the sum of squares of the difference between the right hand side and $N_e$ by setting up two functions with param[1]=b, param[2]=c, param[3]=d and param[4]=T

RHS <- function(par) { 
        Data$No * (1 - exp( (par[3] + par[1]*Data$No) * (par[4]*Data$Ne - 72)
                           / (1 + par[2]*Data$No) ) )
     }
sumsq <- function(par) { 
     sum( ( Data$Ne - RHS(par) )^2 )
     }

We now need to start from some initial estimate, and this one seems to work quite well

param <- c(2,3,4,11)

We can then use optim. We could set up the constraints, either with optim or with constrOptim, but with luck we can find a solution which meets the constraints anyway. We may need to run optim several times until the estimates of the parameters stop changing. For example

> o <- optim(param, sumsq) 
> (param <- o$par)
[1] -0.200  3.825  4.825 11.825
> o <- optim(param, sumsq) 
> (param <- o$par)
[1]  6.081155e-05  1.821490e+00 -2.820120e-02  9.624836e+01
> o <- optim(param, sumsq) 
> (param <- o$par)
[1]  5.278305e-05  3.244559e+00 -2.819326e-02  1.631129e+02
> o <- optim(param, sumsq) 
> (param <- o$par)
[1]  5.278179e-05  3.244649e+00 -2.819326e-02  1.631168e+02
> o <- optim(param, sumsq) 
> (param <- o$par)
[1]  5.278179e-05  3.244649e+00 -2.819326e-02  1.631168e+02

and looking at the details we can see the sum of squares of the differences have been reduced to under 30, which is not bad given there are 89 data points

> o
$par
[1]  5.278179e-05  3.244649e+00 -2.819326e-02  1.631168e+02

$value
[1] 29.79642

$counts
function gradient 
     195       NA 

$convergence
[1] 0

$message
NULL 

Using these parameters we can plot the right hand side of the model against the left hand side with a line to show where $y=x$

plot(RHS(param) ~ Data$Ne)
abline(0,1)

model against actual

and that does not look too bad, though perhaps it should curve a bit lower for high $N_e$. The biggest differences come with the points (No=30,Ne=20) and (No=100,Ne=33).

The problem comes when starting with other initial values for the parameters in the optimisation: some fail completely, some fail the constraints, and most worryingly some settle down on good but different values. For example starting with

param <- c(2,3,4,9)

the repeated optim calls settle down on the parameters

[1]  8.291305e-06  9.521925e-01 -4.292615e-03  3.153295e+02

which are very different to the values found earlier. The sum of squares of the differences are still under 34 and the graph looks almost identical, but if the parameters are supposed to have physical meaning then this is not a good way to estimate them; $Th$ has almost doubled from about 163 to 315.

So the parameters produced can not be trusted.

$\endgroup$
  • $\begingroup$ The given model might be wrong but numerous articles in prestigious journals has been published based on this model. And this is one of the commonly used models in ecology. $\endgroup$ – MYaseen208 Apr 25 '11 at 11:58
  • $\begingroup$ I am not suggesting that the model is wrong as such, but that it is unhelpful in that it is difficult to get accurate estimates of the parameters: different parameters would produce similar results, and you cannot easily get estimates of $N_e$ from $N_o$ since $N_e$ is an input to the calculation of itself. $\endgroup$ – Henry Apr 25 '11 at 12:20
0
$\begingroup$

@Ben Bolker answer the same question on Stack Overflow here. I'm reproducing his answer here:

@Ben Bolker answer

I haven't done it exactly the same way as you, but I think the answer should be OK.

Set up starting conditions

s0 <- list(b=0.1,c=0.1,d=0,Th=3)
X <- read.table("rogersdat.txt",header=TRUE)

Function for computed expected number eaten:

predfun <- function(b,c,d,Th,Te,N0,debug=FALSE) {
  a <- (d+b*N0)/(1+c*N0)
  r <- N0 - (1/(a*Th))*lambertW(a*Th*N0*exp(a*(Th*N0-Te)))
  if (debug) cat(mean(a),b,c,d,Th,mean(r),"\n")
  r
}

Plot the data (just to make sure):

with(X,plot(Ne~N0))
## check starting value
lines(1:100,with(c(s0,X),predfun(b,c,d,Th,Te=72,N0=1:100))) 

Fit:

library(emdbook)
n1 <- nls(Ne~predfun(b,c,d,Th,Te=72,N0),data=X,
    lower=c(1e-6,1e-6,-Inf,1e-6),algorithm="port",
    start=list(b=0.1,c=0.1,d=0.002,Th=3))
summary(n1)

Formula: Ne ~ predfun(b, c, d, Th, Te = 72, N0)

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
b  0.0004155  0.0008745   0.475  0.63591   
c  0.0000010  0.0657237   0.000  0.99999   
d  0.0008318  0.0067374   0.123  0.90203   
Th 4.0639937  1.4665102   2.771  0.00686 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 4.516 on 85 degrees of freedom

Algorithm "port", convergence message: relative convergence (4) 

confint.default(n1)

          2.5 %      97.5 %
b  -0.001298523 0.002129547
c  -0.128815083 0.128817083
d  -0.012373153 0.014036799
Th  1.189686504 6.938300956

Estimates match pretty well, confidence intervals don't (these types of confidence intervals are a little dicey on the boundary anyway ...)

I would actually suggest that maximum likelihood estimation with binomial errors would be a little better.

pdat <- data.frame(N0=1:100,Ne=predict(n1,newdata=data.frame(N0=1:100)))
with(pdat,lines(N0,Ne,col=2))
library(ggplot2)
ggplot(X,aes(x=N0,y=Ne))+stat_sum(aes(size=factor(..n..)),alpha=0.5)+theme_bw()+
  geom_line(data=pdat,colour="red")
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.