[R] Extracting Phi from gls/lme

[Peter Dalgaard]

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