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

Rolf Turner r.turner at auckland.ac.nz
Fri Mar 23 23:30:17 CET 2018


On 24/03/18 10:52, Rolf Turner wrote:
> 
> 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.

Nope!  ***NOT*** bingo!!!  Louise just pointed out to me that what I got 
by applying as.vector(), i.e. 0.3533897, is not the same as the number 
one gets by just typing fit.gls$modelStruct$corStruct, which is 
0.1748786 --- which is what Andrzej's approach gives (and which would
appear to be the right answer).  Bit of a duhhh on my part.

Setting unconstrained=FALSE as in Andrzej's approach seems to be 
crucial.  I missed that bit ....

But it still mystifies me why as.vector() doesn't give the same numeric 
value as is obtained from just printing the object in question.  Also
looking at the code of nlme:::coef.corStruct I would have thought that
this function would simply throw an error when unconstrained is set to 
FALSE.  (Just type "nlme:::coef.corStruct" and you'll see what I mean.)

If I set foo <- nlme:::coef.corStruct and then do

     foo(fit.gls$modelStruct$corStruct)

I do indeed get an error.  WTF?

There would appear to be no way for the results to be different (look
at the code) --- but different they are.

But that sort of nonsense is to be expected when S4 classes and methods 
are involved in any way.

cheers,

Rolf

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276



More information about the R-sig-mixed-models mailing list