Ottorino-Luca Pantani
2009-Aug-17 08:28 UTC
[R] [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? > > KevinActually 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 ?
Douglas Bates
2009-Aug-17 13:08 UTC
[R] [Fwd: Re: R code to reproduce (while studying) Bates & Watts 1988]]]
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. >