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

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 ?

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.
>