[R] Question on CAR appendix on NLS
Ajay Shah
ajayshah at mayin.org
Wed Apr 21 18:12:51 CEST 2004
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)))
More information about the R-help
mailing list