[R-sig-ME] ar (1) with data and output inline
Ben Bolker
bbolker at gmail.com
Wed Jan 26 15:50:48 CET 2011
Hmm. When I look at the autocorrelation functions of the residuals
from the first model (by doing plot(ACF(mod),alpha=0.05) makes it look
like R's estimate of the autocorrelation parameter (0.8) is better than
SAS's (0.2). Furthermore, it looks like something a little more
complicated than AR(1) is going on in the residuals (the ACF goes
significantly negative, which shouldn't happen with AR(1) -- an
ARMA(2,1) model does a little better as measured by AIC (with various
caveats about what the 'effective sample size' is here), but it still
doesn't seem to totally capture the short-term autocorrelation.
You said you believed the SAS fit more on biological (or whatever, I
don't remember) grounds. On what basis ... ? (Not trying to be
difficult, just trying to understand what's going on here.)
If I had more time to play with this I would augment the data with a
"time within run" variable and
xyplot(residuals(mod)~run.time|run,type="l")
to get a better handle -- are there deterministic (smooth) temporal
patterns that you're trying to model with an AR structure?
Ben Bolker
library(nlme)
mod<-lme(y ~ a + b + c, random=~1|name/run, na.action=na.omit,data=eg)
plot(ACF(mod),alpha=0.05)
## However, when trying to introduce a correlation structure in SAS by
altering the 2nd random statement to
## random run(name)/type=ar(1), differences then appear. I've tried
coding this in R:
mod2<-lme(y ~ a + b + c, random=~1|name/run,
correlation=corAR1(form=~1|name/run), na.action=na.omit,
data=eg)
plot(ACF(mod2,resType="normalized"),alpha=0.05)
mod3<-lme(y ~ a + b + c, random=~1|name/run,
correlation=corARMA(p=2,q=1,form=~1|name/run), na.action=na.omit,
data=eg)
AIC(mod,mod2,mod3)
plot(ACF(mod3,resType="normalized"),alpha=0.05)
On 11-01-26 06:16 AM, Paul Chatfield wrote:
> Hopefully this should get through clearly now with data and output - any thoughts of why ar(1) is not agreeing with SAS?
>
>
> Paul
>
> Data accessible from below:
>
> eg<-read.csv("http://www.personal.rdg.ac.uk/~sns97aal/eg.csv", header=T, na=".")
> eg$run<-as.factor(eg$run)
>
> Output below:
> SAS output for first model (no autocorrelation):
>
> Covariance Parameter
> Estimates
>
> Cov Parm Estimate
>
> name 0.003483
> run(name) 0.01600
> Residual 0.0288
>
> Solution for Fixed Effects
>
> Standard
> Effect Estimate Error DF t Value Pr > |t|
>
> Intercept 1.7264 0.03213 4 53.73 <.0001
> a -0.00137 0.004333 1114 -0.32 0.7515
> b 0.6129 0.2150 1114 2.85 0.0045
> c -0.08898 0.01455 1114 -6.11 <.0001
>
>
> R output which agrees
>
> Random effects:
> Formula: ~1 | name
> (Intercept)
> StdDev: 0.05900847
>
> Formula: ~1 | run %in% name
> (Intercept) Residual
> StdDev: 0.1265071 0.1697051
>
> Value Std.Error DF t-value p-value
> (Intercept) 1.726390 0.03212781 1114 53.73507 0.0000
> cct1.1 -0.001371 0.00433359 1114 -0.31636 0.7518
> k21.1 0.612844 0.21502840 1114 2.85006 0.0045
> pup01.1 -0.088977 0.01455320 1114 -6.11391 0.0000
>
> SAS output for 2nd model (with autocorrelation):
>
> Covariance Parameter
> Estimates
>
> Cov Parm Estimate
>
> name 0.002904
> Variance 0.01671
> AR(1) 0.2038
> Residual 0.02880
>
> Solution for Fixed Effects
>
> Standard
> Effect Estimate Error DF t Value Pr > |t|
>
> Intercept 1.7274 0.03279 4 52.67 <.0001
> a -0.00134 0.004366 1114 -0.31 0.7583
> b 0.6034 0.2154 1114 2.80 0.0052
> c -0.09075 0.01566 1114 -5.79 <.0001
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list