[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