[R-sig-ME] ar (1) with data and output inline

[Ben Bolker]

  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