[R-sig-ME] [FORGED] Re: Question about gls

[Rolf Turner]

On 24/03/18 08:58, Andrzej Galecki wrote:

> Hi Louise,
> 
> One possible way ...
> 
> mSt <- fit.gls$modelStruct
> cSt <- mSt$corStruct
> coef(cSt, unconstrained = FALSE)
> 
> #       Phi
> # 0.1748786

Not exactly obvious, but, is it?  How on earth would one ever figure 
this out, except by appealing to r-sig-mixed-models?

Looking at str(fit.gls) one sees that there is a (list) component 
"modelStruct" and that this component in turn has a component 
"corStruct".  If one is an optimist, one might try

     fit.gls$modelStruct$corStruct

which (mirabile dictu!) gives:

Correlation structure of class corAR1 representing
       Phi
0.1748786

Again if one is an optimist, one might try stripping away all the 
unwanted baggage by doing

     as.vector(fit.gls$modelStruct$corStruct)

giving:

[1] 0.3533897

Bingo.

Still it's all pretty obscure.  And typical of the obfuscation that 
results from using S4 classes and methods.

cheers,

Rolf

> 
> On Fri, Mar 23, 2018 at 3:27 PM, Louise Ryan <Louise.M.Ryan at uts.edu.au>
> wrote:
> 
>> Hi there,
>>
>> I am using your function gls in R to fit some AR models (it is part of a
>> larger thing I am trying to do). I can see from the printed output that the
>> function is estimating the AR parameter.  But I CANNOT for the life of me
>> figure out how to extract this parameter so I can use it for something else.
>>
>> Here’s a bit of code that simulates some data and runs the model.
>>
>> #####
>> # Set true parameter values:
>> beta0True <- 0.3
>> beta1True <- 1.6
>> sigmaTrue <- 0.5
>> rhoTrue <- 0.2
>>
>> # Generate data:
>> set.seed(1)
>> n <- 500
>> x <- seq(0,1,length=n)
>> epsilon <- rep(NA,n)
>> epsilon[1] <- sigmaTrue*rnorm(1)
>> for (i in 2:n) {epsilon[i] <- rhoTrue*epsilon[i-1] + sigmaTrue*rnorm(1)}
>> y <- beta0True + beta1True*x + epsilon
>>
>> fit.gls <- gls(y~x,  correlation=corAR1())
>> print(summary(fit.gls))
>> ###########################
>>
>> For my little simulation, the parameter Phi is being estimated as .1748786.
>>
>> I have looked at attributes(fit.gls) and just cannot figure out where this
>> estimate is stored!
>>
>> I’d appreciate your help!