[R] [Fwd: Re: R code to reproduce (while studying) Bates & Watts 1988]]]
bates at stat.wisc.edu
Mon Aug 17 15:08:58 CEST 2009
On Mon, Aug 17, 2009 at 3:28 AM, Ottorino-Luca
Pantani<ottorino-luca.pantani at unifi.it> wrote:
> Kevin Wright wrote:
>> 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.
>> lag.plot(resid(m2), do.lines=FALSE)
>> m3 <- update(m2, corr=corAR1(.67))
>> lag.plot(resid(m3), do.lines=FALSE)
>> The residual plots for model m3 still show structure, unlike in Bates &
>> Watts, so maybe this is not the correct model?
> 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?
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
> 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),
> 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
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help