R help - Aug 2009 - R code to reproduce (while studying) Bates & Watts 1988

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/

library(nlme)
m2 <- gnls(conc  ~ t1*(1-t2*exp(-k*time)),
           data  = df.Chloride,
           start =  list(
             t1 = 35,
             t2 = 0.91,
             k = 0.22))
summary(m2)
plot(m2)
lag.plot(resid(m2), do.lines=FALSE)
acf(resid(m2))

m3 <- update(m2, corr=corAR1(.67))
summary(m3)
plot(m3)
lag.plot(resid(m3), do.lines=FALSE)
acf(resid(m3))

The residual plots for model m3 still show structure, unlike in Bates &
Watts, so maybe this is not the correct model?

Kevin



On Thu, Aug 13, 2009 at 6:11 AM, Ottorino-Luca Pantani <
ottorino-luca.pantani@unifi.it> wrote:
> 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@unifi.it
> http://www4.unifi.it/dssnp/
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]