[R] R code to reproduce (while studying) Bates & Watts 1988
Ottorino-Luca Pantani
ottorino-luca.pantani at unifi.it
Thu Aug 13 13:11:56 CEST 2009
Hi R users,
I'm here trying to understand correlated residuals in nonlinear estimation.
I'm reading/studying the book Bates, D. M. and D. G. Watts, (1988),
/Nonlinear regression analysis and its applications/, Wiley, NY. pages
92-94, trying to reproduce the figures and to find out the code in R to
perform the necessary calculations.
I also consulted Pinheiro and Bates, but without success
Here below you'll find are my efforts.
I'm in trouble at plotting the lag plot (fig 3.6 b) and in fitting the
new model with autocorrelated residuals (most probably misusing
"correlation = corAR1()" in updating a nls model).
Could someone be so kind to help me ?
Thanks a lot
Ottorino
"df.Chloride"<-
structure(.Data = list(time = c(2.45, 2.55, 2.65, 2.75, 2.85, 2.95,
3.05, 3.15,
3.25, 3.35, 3.45, 3.55, 3.65, 3.75, 3.85, 3.95, 4.05, 4.15, 4.25, 4.35,
4.45, 4.55, 4.65, 4.75, 4.85, 4.95, 5.05, 5.15, 5.25, 5.35, 5.45, 5.55,
5.65, 5.75, 5.85, 5.95, 6.05, 6.15, 6.25, 6.35, 6.45, 6.55, 6.65, 6.75,
6.85, 6.95, 7.05, 7.15, 7.25, 7.35, 7.45, 7.55, 7.65, 7.75), conc = c(
17.3, 17.6, 17.9, 18.3, 18.5, 18.9, 19, 19.3, 19.8, 19.9, 20.2, 20.5,
20.6, 21.1, 21.5, 21.9, 22, 22.3, 22.6, 22.8, 23, 23.2, 23.4, 23.7,
24, 24.2, 24.5, 25, 25.4, 25.5, 25.9, 25.9, 26.3, 26.2, 26.5, 26.5,
26.6, 27, 27, 27, 27, 27.3, 27.8, 28.1, 28.1, 28.1, 28.4, 28.6, 29,
29.2, 29.3, 29.4, 29.4, 29.4)), row.names = c("1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17",
"18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29",
"30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41",
"42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53",
"54"), class = "data.frame", reference = "A1.9, p. 274")
##Figure 3.5, pag 92
plot(df.Chloride$time, df.Chloride$conc, pch = 20)
## fitting the model
nls.1 <- nls(conc ~ t1*(1-t2*exp(-k*time)),
data = df.Chloride,
start = list(
t1 = 35,
t2 = 0.91,
k = 0.22))
##Figure 3.56a, pag 93
plot(nls.1, pch = 20)
##Figure 3.56b, pag 93
???????????not the foggiest idea
##Figure 3.57, pag 94
acf(resid(nls.1), xlim = c(1, 15), lab = c(5, 4, 7))
##Try to fit a model with autocorrelated residues (Problems !!)
nls.2 <- update(nls.1, corr = corAR1(0.67))
--
Ottorino-Luca Pantani, Università di Firenze
Dip. Scienza del Suolo e Nutrizione della Pianta
P.zle Cascine 28 50144 Firenze Italia
Tel 39 055 3288 202 (348 lab) Fax 39 055 333 273
OLPantani at unifi.it http://www4.unifi.it/dssnp/
More information about the R-help
mailing list