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

[Berwin A Turlach]

G'day Rolf,

On Sat, 24 Mar 2018 11:30:17 +1300
Rolf Turner <r.turner at auckland.ac.nz> wrote:

> 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
> [...] 0.1748786 [...] Bit of a duhhh on my part.

:-)

> 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.  

Because, as you said, you are printing the object. Essentially,
as.vector() just returns the value stored in the object will all
attributes stripped.  Printing on the other hand might call some
function, and here it does:

R> class(fit.gls$modelStruct$corStruct)
[1] "corAR1"    "corStruct"
R> getAnywhere("print.corAR1")
no object named ‘print.corAR1’ was found
R> getAnywhere("print.corStruct")
A single object matching ‘print.corStruct’ was found
It was found in the following places
  registered S3 method for print from namespace nlme
  namespace:nlme
with value
[...]

> 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?

What error did you get?  I do not get an error:

R> foo <- nlme:::coef.corStruct
R> foo(fit.gls$modelStruct$corStruct)
[1] 0.3533897

Or did you mean

R> foo(fit.gls$modelStruct$corStruct, unconstrained=FALSE)
Error in foo(fit.gls$modelStruct$corStruct, unconstrained = FALSE) : 
  do not know how to obtain parameters of “corAR1” object
>
 There would appear to be no way for the results to be different (look
> at the code) --- but different they are.

Well, you are missing that 

R> class(fit.gls$modelStruct$corStruct)
[1] "corAR1"    "corStruct"

And S3 will try to dispatch initially on the first class, then the
second class, and so forth...

And coef.corAR1 exists in the nlme package:

R> nlme:::coef.corAR1
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
}
<bytecode: 0x2935e60>
<environment: namespace:nlme>

From this code it should now be clear how one can change from the
unconstrained parameter formulation used when fitting the model to the
parameter one is actually interested in for an AR1 error model:

R> (tmp <- as.numeric( fit.gls$modelStruct$corStruct))
[1] 0.3533897
R> (tmp <- exp(tmp))
[1] 1.423886
R> (tmp-1)/(tmp+1)
[1] 0.1748786
 
Cheers,

	Berwin