[R] Question on CAR appendix on NLS

[Ajay Shah]

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
-------------- next part --------------
# John Fox has a book "An R and S+ companion to applied regression"
# (abbreviated CAR).
#
# An appendix associated with this book, titled
#   "Nonlinear regression and NLS"
# is up on the web, and I strongly recommend that you go read it.
#
# All the data and code of this program is from there.

# First take some data - from the CAR book --
library(car)
data(US.pop)
attach(US.pop)
plot(year, population, type="l", col="blue")

# So you see, we have a time-series of the US population. We want to
# fit a nonlinear model to it.

library(nls)                            # The nonlinear regression library.
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)
# Note that you just write in the formula that you want to fit,
# and supply starting values. "trace=T" makes him show iterations go by.

summary(pop.mod)
# Add in predicted values into the plot
lines(year, fitted.values(pop.mod), lwd=3, col="red")
z <- locator(1)

# Look at residuals
plot(year, residuals(pop.mod), type="b")
abline(h=0, lty=2)
z <- locator(1)

# Using analytical derivatives:
model = function(beta1, beta2, beta3, time) {
   m = beta1/(1+exp(beta2+beta3*time))
   term = exp(beta2 + beta3*time)
   gradient = cbind((1+term)^-2,
     -beta1*(1+term)^-2 * term,
     -beta1*(1+term)^-2 * term * time)
   attr(m, 'gradient') <- gradient
   m
 }

summary(nls(population ~ model(beta1, beta2, beta3, time),
            start=list(beta1=350, beta2=4.5, beta3=-0.3)))

# Using analytical derivatives, using automatic differentiation (!!!):
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)))