[R] gamm in mgcv random effect significance

[Gavin Simpson]

On Fri, 2013-06-07 at 13:12 -0700, William Shadish wrote:
> Dear R-helpers,
> 
> I'd like to understand how to test the statistical significance of a 
> random effect in gamm. I am using gamm because I want to test a model 
> with an AR(1) error structure, and it is my understanding neither gam 
> nor gamm4 will do the latter.

gamm4() can't yes and out of the box mgcv::gam can't either but
see ?magic for an example of correlated errors and how the fits can be
manipulated to take the AR(1) (or any structure really as far as I can
tell) into account.

You might like to look at mgcv::bam() which allows an known AR(1) term
but do check that it does what you think; with a random effect spline
I'm not at all certain that it will nest the AR(1) in the random effect
level.

<snip />
> Consider, for example, two models, both with AR(1) but one allowing a 
> random effect on xc:
> 
> g1 <- gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, 
> correlation=corAR1())
> g2 <- gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random = 
> list(xc=~1),correlation=corAR1())

Shouldn't you specify how the AR(1) is nested in the hierarchy here,
i.e. AR(1) within xc? maybe I'm not following your data structure
correctly.

> I include the output for g1 and g2 below, but the question is how to 
> test the significance of the random effect on xc. I considered a test 
> comparing the Log-Likelihoods, but have no idea what the degrees of 
> freedom would be given that s(xc) is smoothed. I also tried:
> 
> anova(g1$gam, g2$gam)

gamm() fits via the lme() function of package nlme. To do what you want,
you need the anova() method for objects of class "lme", e.g.

anova(g1$lme, g2$lme)

Then I think you should check if the fits were done via REML and also be
aware of the issue of testing wether a variance term is 0.

> that did not seem to return anything useful for this question.
> 
> A related question is how to test the significance of adding a second 
> random effect to a model that already has a random effect, such as:
> 
> g3 <- gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = 
> list(Case=~1, z=~1),correlation=corAR1())
> g4 <- gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = 
> list(Case=~1, z=~1, int=~1),correlation=corAR1())

Again, I think you need anova() on the $lme components.

HTH

G

> Any help would be appreciated.
> 
> Thanks.
> 
> Will Shadish
> ********************************************
> g1
> $lme
> Linear mixed-effects model fit by maximum likelihood
>    Data: data
>    Log-likelihood: -437.696
>    Fixed: fixed
> X(Intercept)           Xz         Xint    Xs(xc)Fx1
>     0.6738466   -2.5688317    0.0137415   -0.1801294
> 
> Random effects:
>   Formula: ~Xr - 1 | g
>   Structure: pdIdnot
>                   Xr1          Xr2          Xr3 Xr4          
> Xr5          Xr6          Xr7          Xr8 Residual
> StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 
> 0.0004377781 0.0004377781 0.0004377781 1.693177
> 
> Correlation Structure: AR(1)
>   Formula: ~1 | g
>   Parameter estimate(s):
>        Phi
> 0.3110725
> Variance function:
>   Structure: fixed weights
>   Formula: ~invwt
> Number of Observations: 264
> Number of Groups: 1
> 
> $gam
> 
> Family: binomial
> Link function: logit
> 
> Formula:
> y ~ s(xc) + z + int
> 
> Estimated degrees of freedom:
> 1  total = 4
> 
> attr(,"class")
> [1] "gamm" "list"
> ****************************
>  > g2
> $lme
> Linear mixed-effects model fit by maximum likelihood
>    Data: data
>    Log-likelihood: -443.9495
>    Fixed: fixed
> X(Intercept)           Xz         Xint    Xs(xc)Fx1
>   0.720018143 -2.562155820  0.003457463 -0.045821030
> 
> Random effects:
>   Formula: ~Xr - 1 | g
>   Structure: pdIdnot
>                   Xr1          Xr2          Xr3 Xr4          
> Xr5          Xr6          Xr7          Xr8
> StdDev: 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 
> 7.056078e-06 7.056078e-06 7.056078e-06
> 
>   Formula: ~1 | xc %in% g
>           (Intercept) Residual
> StdDev: 6.277279e-05 1.683007
> 
> Correlation Structure: AR(1)
>   Formula: ~1 | g/xc
>   Parameter estimate(s):
>        Phi
> 0.1809409
> Variance function:
>   Structure: fixed weights
>   Formula: ~invwt
> Number of Observations: 264
> Number of Groups:
>          g xc %in% g
>          1        34
> 
> $gam
> 
> Family: binomial
> Link function: logit
> 
> Formula:
> y ~ s(xc) + z + int
> 
> Estimated degrees of freedom:
> 1  total = 4
> 
> attr(,"class")
> [1] "gamm" "list"
> 
> 

-- 
Gavin Simpson, PhD                          [t] +1 306 337 8863
Adjunct Professor, Department of Biology    [f] +1 306 337 2410
Institute of Environmental Change & Society [e] gavin.simpson at uregina.ca
523 Research and Innovation Centre          [tw] @ucfagls
University of Regina
Regina, SK S4S 0A2, Canada