Message-ID: <ff98e7c6-0b02-fa91-6e47-26773eb743e4@auckland.ac.nz>
Date: 2018-03-23T21:52:34Z
From: Rolf Turner
Subject: [FORGED] Re: Question about gls
In-Reply-To: <CA+XOvOTtW8L_yVpteRdsVg0iPS1v902jwYAf1TiR95H4Pt=azA@mail.gmail.com>
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!