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)))
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 >
Ajay,
I beleve there is an error in your hand-calculated derivatives.
model = function(beta1, beta2, beta3, time) {
m = beta1/(1+exp(beta2+beta3*time))
term = exp(beta2 + beta3*time)
gradient = cbind((1+term)^-2, <-------- the exponent should be -1
-beta1*(1+term)^-2 * term,
-beta1*(1+term)^-2 * term * time)
attr(m, 'gradient') <- gradient
m
}
Cheers,
Andy
__________________________________
Andy Jaworski
518-1-01
Process Laboratory
3M Corporate Research Laboratory
-----
E-mail: apjaworski at mmm.com
Tel: (651) 733-6092
Fax: (651) 736-3122
Ajay Shah
<ajayshah at mayin.o
rg> To
Sent by: r-help <r-help at
stat.math.ethz.ch>
r-help-bounces at st cc
at.math.ethz.ch
Subject
[R] Question on CAR appendix on NLS
04/21/2004 11:12
AM
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
[attachment "nls.R" deleted by Andrzej P. Jaworski/US-Corporate/3M/US]
______________________________________________
R-help at stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html