Extracting Phi from gls/lme
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