[R-sig-ME] extract phi parameter from lme model

[Ben Bolker]

On 13-12-03 12:27 PM, Rasmussen, Paul W - DNR wrote:
> Hi,
> 
> Can anyone tell me how to extract an AR(1) parameter from an lme
> model?
> 
> Here is a summary of the model:
> 
>> summary(WAElme.modAR1)
> Linear mixed-effects model fit by maximum likelihood Data:
> WAElmeData AIC       BIC   logLik -103.6577 -83.43384 57.82883
> 
> Random effects: Formula: ~1 | WBIC (Intercept)  Residual StdDev:
> 0.08148838 0.2578487
> 
> Correlation Structure: ARMA(1,0) Formula: ~year | WBIC Parameter
> estimate(s): Phi1 0.7262637 Fixed effects: log10pe ~ year + area 
> Value  Std.Error  DF  t-value p-value (Intercept)  3.468126
> 0.08024473 203 43.21937  0.0000 year        -0.009656 0.00491928 203
> -1.96292  0.0510 area         0.276664 0.04673097   9  5.92035
> 0.0002 Correlation: (Intr) year year -0.823 area -0.034  0.065
> 
> Standardized Within-Group Residuals: Min           Q1          Med
> Q3          Max -3.182366526 -0.526080732 -0.002701109  0.573640884
> 2.406525772
> 
> Number of Observations: 215 Number of Groups: 11
> 
> I can print the value of the phi parameter as follows:
> 
>> WAElme.modAR1$modelStruct$corStruct
> Correlation structure of class corARMA representing Phi1 0.7262637
> 
> but I am not able to use this in further calculations.  Can anyone
> tell me how to extract this as a scalar for use in further
> calculations?
> 
> Thank you, Paul Rasmussen Wisconsin DNR
> 

A reproducible example is always nice ...

set.seed(101)
d <- data.frame(y=rnorm(100),x=rnorm(100),f=factor(rep(1:10,10)))
library(nlme)
m <- lme(y~x,random=~1|f,correlation=corAR1(),data=d)

cS
str(cS <- m$modelStruct$corStruct)
## ... something complicated
coef(cS)  ## not same as printed value
?coef.corStruct
coef(cS,unconstrained=FALSE)