Re: [R] Extracting Phi from gls/lme

From: John Logsdon <j.logsdon_at_quantex-research.com>
Date: Fri 14 Jul 2006 - 00:03:23 EST


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@quantex-research.com              a.einstein@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@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@quantex-research.com a.einstein@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@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@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@biostat.ku.dk) FAX: (+45) 35327907
>



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Jul 14 00:12:39 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 14 Jul 2006 - 02:14:15 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.