[R] Question on CAR appendix on NLS

[John Fox]

Dear Alay,

I'm leaving town this morning for several days, so I don't have time to
check through your code, but I did rerun the examples from the appendix (see
below), and all three approaches produce identical parameter correlations
(using R 1.9.0 under Win XP). Perhaps you made an error in entering a
derivative.

I hope this helps,
 John

> library(car)
> data(US.pop)
> attach(US.pop)
> time <- 0:20
> pop.mod <- nls(population ~ beta1/(1 + exp(beta2 + beta3*time)), 
+     start=list(beta1=350, beta2=4.5, beta3=-0.3),
+     trace=T)
13007.48 :  350.0   4.5  -0.3 
609.5727 :  351.8074862   3.8405002  -0.2270578 
365.4396 :  383.7045367   3.9911148  -0.2276690 
356.4056 :  389.1350277   3.9897242  -0.2265769 
356.4001 :  389.1462893   3.9903758  -0.2266276 
356.4001 :  389.1665304   3.9903412  -0.2266193 
356.4001 :  389.1655126   3.9903457  -0.2266199 
> summary(pop.mod)

Formula: population ~ beta1/(1 + exp(beta2 + beta3 * time))

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 389.16551   30.81197   12.63 2.20e-10 ***
beta2   3.99035    0.07032   56.74  < 2e-16 ***
beta3  -0.22662    0.01086  -20.87 4.60e-14 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 4.45 on 18 degrees of freedom

Correlation of Parameter Estimates:
        beta1   beta2
beta2 -0.1662        
beta3  0.9145 -0.5407

> 
> 
>  model <- function(beta1, beta2, beta3, time){
+     model <- beta1/(1 + exp(beta2 + beta3*time))
+     term <- exp(beta2 + beta3*time)
+     gradient <- cbind((1 + term)^-1, # in proper order
+         -beta1*(1 + term)^-2 * term,
+         -beta1*(1 + term)^-2 * term * time)
+     attr(model, "gradient") <- gradient
+     model
+     }
>     
>  summary(nls(population ~ model(beta1, beta2, beta3, time),
+     start=list(beta1=350, beta2=4.5, beta3=-0.3)))

Formula: population ~ model(beta1, beta2, beta3, time)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 389.16551   30.81196   12.63 2.20e-10 ***
beta2   3.99035    0.07032   56.74  < 2e-16 ***
beta3  -0.22662    0.01086  -20.87 4.60e-14 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 4.45 on 18 degrees of freedom

Correlation of Parameter Estimates:
        beta1   beta2
beta2 -0.1662        
beta3  0.9145 -0.5407

> 
>  model <- deriv(~ beta1/(1 + exp(beta2 + beta3*time)), # rhs of model
+     c("beta1", "beta2", "beta3"), # parameter names
+     function(beta1, beta2, beta3, time){} # arguments for result
+     )
>  summary(nls(population ~ model(beta1, beta2, beta3, time),
+  start=list(beta1=350, beta2=4.5, beta3=-0.3)))

Formula: population ~ model(beta1, beta2, beta3, time)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 389.16551   30.81196   12.63 2.20e-10 ***
beta2   3.99035    0.07032   56.74  < 2e-16 ***
beta3  -0.22662    0.01086  -20.87 4.60e-14 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 4.45 on 18 degrees of freedom

Correlation of Parameter Estimates:
        beta1   beta2
beta2 -0.1662        
beta3  0.9145 -0.5407

> 
>  

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Ajay Shah
> Sent: Wednesday, April 21, 2004 11:13 AM
> To: r-help
> Subject: [R] Question on CAR appendix on NLS
> 
> The PDF file on the web, which is an appendix on nonlinear 
> regression associated with the CAR book, is very nice.
> 
> When I ran through the code presented there, I found 
> something odd. The code does a certain model in 3 ways: 
> Vanilla NLS (using numerical differentation), Analytical 
> derivatives (where the user supplies the derivatives) and 
> analytical derivatives (using automatic differentiation). The 
> three results agree, except for the correlation of parameter 
> estimates :
> 
>           beta1   beta2
>   beta2 -0.1662                      Numerical derivatives
>   beta3  0.9145 -0.5407
> 
>   beta2 -0.7950                      Analytical derivatives
>   beta3  0.9145 -0.9661
> 
>   beta2 -0.1662                      Automatic differentiation
>   beta3  0.9145 -0.5407
> 
> Is this just a glitch of a small sample, or should I worry? 
> My source file (which should be the same as John Fox's file; 
> I typed it in while reading the PDF file, and made minor 
> changes) is attached.
> 
> -- 
> Ajay Shah                                                   Consultant
> ajayshah at mayin.org                      Department of Economic Affairs
> http://www.mayin.org/ajayshah           Ministry of Finance, New Delhi
>