[R] Extracting Phi from gls/lme

[John Logsdon]

Peter

Yes spelling would help!  I transcribed this from one screen to another
and still the penny didn't drop.  And the transformation works too.

getAnywhere() is yet another useful function that I hadn't heard about...

Many thanks for saving my house.

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:
> 
> > 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...
> 
> Works for me. Is there a rogue coef() around?? Or did you misspell
> "unconstrained" and not tell us?
> 
> > coef(x,uncostrained=FALSE)
> [1] 0.1171201
> > coef(x,unconstrained=FALSE)
>        Phi
> 0.05849318
>  
> 
>  
> > This is not an obvious transformation of the -0.1996813 I think.
> 
> Obvious is in the eyes of the beholder, but:
> 
> > aux <- exp(-0.4048011)
> > (aux - 1) / (aux + 1)
> [1] -0.1996813
> 
> > 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?
> 
> I don't think they do anything, but their print methods call
> print.corStruct which calls coef.corAR1. And
> 
> > getAnywhere(coef.corAR1)
> A single object matching coef.corAR1 was found
> It was found in the following places
>   registered S3 method for coef from namespace nlme
>   namespace:nlme
> with value
> 
> function (object, unconstrained = TRUE, ...)
> {
>     if (unconstrained) {
>         if (attr(object, "fixed")) {
>             return(numeric(0))
>         }
>         else {
>             return(as.vector(object))
>         }
>     }
>     aux <- exp(as.vector(object))
>     aux <- c((aux - 1)/(aux + 1))
>     names(aux) <- "Phi"
>     aux
> }
> <environment: namespace:nlme>
> > aux <- exp(-0.4048011
> 
> 
>  
> > 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
> > > 
> > 
> > 
> 
> -- 
>    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
>