[R] Extracting Phi from gls/lme

John Logsdon j.logsdon at quantex-research.com
Thu Jul 13 14:29:47 CEST 2006


Peter

This looks very promising:

x<-mod.lme$model$Struct$corStruct
Correlation structure of class corAR1 representing
       Phi
-0.1996813

which is the value I want.  

Yippee (save the bricks)

But:

coef(x,unconstrained=FALSE)
[1] -0.4048011

and any attempt to coerce x into a scalar always returns -0.404...

This is not an obvious transformation of the -0.1996813 I think.

Looking at str(x) returns the first line:

Classes 'corAR1', 'corStruct' atomic [1:1] -0.405

Not a -0.199 ... in sight in the attributes various.

How does summary.lme/gls do it?

Best wishes

John

John Logsdon                               "Try to make things as simple
Quantex Research Ltd, Manchester UK         as possible but not simpler"
j.logsdon at quantex-research.com              a.einstein at relativity.org
+44(0)161 445 4951/G:+44(0)7717758675       www.quantex-research.com


On 13 Jul 2006, Peter Dalgaard wrote:

> John Logsdon <j.logsdon at quantex-research.com> writes:
> 
> > I am trying to extract into a scalar the value of Phi from the printed
> > output of gls or lme using corAR1 correlation.  ie I want the estimate of
> > the autocorrelation.  I can't see how to do this and haven't seen it
> > anywhere in str(model.lme).
> > 
> > I can get all the other information - fixed and random effects etc.
> > 
> > Is there an obvious way so that I can save the brick wall some damage?
> 
> Just be soft in the head...
> 
> Seriously, I think the recipe is to drill down into model.lme until
> you find the corAR1 class object,  like this, I think.
> 
> x <- model.lme$modelStruct$corStruct
> coef(x,unconstrained=FALSE). 
> 
> -- 
>    O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
>   c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
>  (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
> ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907
>



More information about the R-help mailing list