1
$\begingroup$

I am using generalized additive models in R with an correlated structure on the residuals, however I am unable to extract the actual coefficients of the AR(p) structure even though they are provide in the output. Here is an example of an output and I am interested in being able to extract the Phi value of 0.6100947.

>GAM_sim1<-(gamm(y_sim~s(x_sim, bs="cs", k=40), correlation = corARMA(p=length(AR))))
>GAM_sim1$lme

Linear mixed-effects model fit by maximum likelihood
  Data: strip.offset(mf) 
  Log-likelihood: -250.7398
  Fixed: y ~ X - 1 
        X 
-7.442257 

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
               Xr1        Xr2        Xr3        Xr4        Xr5        Xr6        Xr7        Xr8
StdDev: 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312
               Xr9       Xr10       Xr11       Xr12       Xr13       Xr14       Xr15       Xr16
StdDev: 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312
              Xr17       Xr18       Xr19       Xr20       Xr21       Xr22       Xr23       Xr24
StdDev: 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312
              Xr25       Xr26       Xr27       Xr28       Xr29       Xr30       Xr31       Xr32
StdDev: 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312
              Xr33       Xr34       Xr35       Xr36       Xr37       Xr38       Xr39 Residual
StdDev: 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 0.01667312 1.557552

Correlation Structure: AR(1)
 Formula: ~1 | g 
 Parameter estimate(s):
      Phi 
0.6100947 
Number of Observations: 151
Number of Groups: 1 

The summary above can be shorten even more by using,

>GAM_sim1$lme$modelStruct$corStruct
Correlation structure of class corAR1 representing
      Phi 
0.6100947

And now this is where the fun happens as when I try to extract this value, I get a completely different value and no further information is given.

>GAM_sim1$lme$modelStruct$corStruct[1]
[1] 1.418144

I haven't been able to figure out where the 1.4 comes from and even if I expand this out to auto-correlated errors of an AR(2) then the same situation happens when I try to extract the two Phi values.

> GAM_sim2<-(gamm(y_sim~s(x_sim, bs="cs", k=40), correlation = corARMA(p=2)))
> GAM_sim2$lme$modelStruct$corStruct
Correlation structure of class corARMA representing
     Phi1      Phi2 
0.5257324 0.1708431 
> GAM_sim2$lme$modelStruct$corStruct[1:2]
[1] 1.496342 0.345070

Does anyone know how to extract the AR values from the GAMM function in the mgcv package as I haven't been able to get a proper answer in the documentation?

Thank you

$\endgroup$

1 Answer 1

1
$\begingroup$

Try coef(GAM_sim2$lme$modelStruct$corStruct, unconstrained = FALSE).

I found this by looking through the {nlme} functions nlme:::print.corStruct and nlme:::coef.corAR1. The latter function is

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
}

Not sure exactly sure why the transform is doing that. It resembles a hyperbolic tangent function.

$\endgroup$
2
  • $\begingroup$ Thank you Alex, I'll check this out to see if this resolves the issue and see if it also works an AR(2) and higher order. $\endgroup$ Feb 2, 2023 at 0:23
  • $\begingroup$ I checked it out this morning and this resolved the issue. Here is a copy of the results when setting unconstrained either as "TRUE" or "FALSE." By setting this to be FALSE I'm able to extract the correct AR values. > coef(GAM_sim2$lme$modelStruct$corStruct, unconstrained = TRUE) [1] 1.496342 0.345070 > coef(GAM_sim2$lme$modelStruct$corStruct, unconstrained = FALSE) Phi1 Phi2 0.5257324 0.1708431 Thank you very much, I struggled for about 2 hours trying to get a work around to work. $\endgroup$ Feb 2, 2023 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.