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

[Douglas Bates]

On Mon, Aug 17, 2009 at 3:28 AM, Ottorino-Luca
Pantani<ottorino-luca.pantani at unifi.it> wrote:
> 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 have been without Internet access for a couple of weeks (in the US
AT&T is now competing with the cable companies and has succeeded in
emulating the cable TV companies' terrible service) and I missed the
beginning of this discussion.  Is there a reason that you have not
used the example in the NRAIA package to obtain that model fit?

Try

install.packages("NRAIA")
library(NRAIA)
example(Chloride)

That gets you the data and the model fit without the AR1 correlation.
I guess that I didn't put in the general optimization in that example
yet.

> 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 ?
>
> ______________________________________________
> R-help at 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.
>