[R] [Fwd: Re: R code to reproduce (while studying) Bates & Watts 1988]]]

Ottorino-Luca Pantani ottorino-luca.pantani at unifi.it
Mon Aug 17 10:28:05 CEST 2009


Kevin Wright wrote:

> library(nlme)
> m2 <- gnls(conc  ~ t1*(1-t2*exp(-k*time)),
>            data  = df.Chloride,
>            start =  list(
>              t1 = 35,
>              t2 = 0.91,
>              k = 0.22))
So my error was to use nls instead that gnls. Thanks a lot, Kevin.
> 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
Actually is not even the same model described in the book.

There the model is told to have the following parameters;
t1  37.58
t2    0.849
k     0.178
Phi   0.69,

while the one obtained with R has the following
t1  38.98
t2   0.825
k     0.158
Phi  0.682.

I run this code, but without success. I obtain again the m3 model.

m4 <- gnls(conc  ~ t1*(1-t2*exp(-k*time)),
        data  = df.Chloride,
        start =  list(
          t1 = 37.58,
          t2 = 0.849,
          k = 0.178),
        corr=corAR1(.69))

I cannot understand why

anova(m2,  m3)
Model df       AIC       BIC   logLik   Test  L.Ratio p-value
m2     1  4 -20.09053 -12.13459 14.04526                       
m3     2  5 -51.08761 -41.14269 30.54380 1 vs 2 32.99708  <.0001
indicate a susbstantial improvement in the model, but the plot of 
residual is quite the same (slight differences)

plot(m2,  pch  = 20);x11();plot(m3,  pch  = 20)

Any other idea ?




More information about the R-help mailing list