[R] gamm in mgcv random effect significance

[Simon Wood]

I would be very nervous about relying on an anova call here. It will 
attempt a generalized likelihood ratio test, but gamm is using penalized 
quasi likelihood and there is really no likelihood here (even without 
the problem that if there was a likelihood the null hypothesis would 
still be on the edge of the feasible parameter space making the GLRT 
problematic). The best hope might be to model the random effect of xc 
using a term s(xc,bs="re") in the model formula (xc will need to be a 
factor for this), and then use summary on the gam part of the fitted 
model object to assess significance. If you do this you'll need to 
include the grouping factor explicitly in corAR1 (at present it's picked 
up from the random effect, so is nested in xc).
i.e
   g2 <- gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random =
list(xc=~1),correlation=corAR1())
becomes something like...
   xf <- factor(xc)
   g2 <- gamm(y ~ s(xc) +z+ int + s(xf,bs="re"),family=binomial, 
weights=trial,
             correlation=corAR1(form=~1|xf))
   summary(g2$gam)

... I'm also a bit nervous about xc entering as an iid random effect and 
the argument of a smooth, however - does that model structure really 
make sense?

best,
Simon

On 11/06/13 18:08, William Shadish wrote:
> Gavin et al.,
>
> Thanks so much for the help. Unfortunately, the command
>
> > anova(g1$lme, g2$lme)
>
> gives "Error in eval(expr, envir, enclos) : object 'fixed' not found
>
> and for bam (which is the one that can use a known ar1 term), the 
> error is
>
> > AR1 parameter rho unused with generalized model
>
> Apparently it cannot run for binomial distributions, and presumably 
> also Poisson.
>
> I did find a Frequently Asked Questions for package mgcv page that said
>
> "How can I compare gamm models? In the identity link normal errors 
> case, then AIC and hypotheis testing based methods are fine. Otherwise 
> it is best to work out a strategy based on the summary.gam"
>
> So putting all this together, I take it that my binomial example will 
> not support a direct model comparison to test the significance of the 
> random effects. I'm guessing the best strategy based on the 
> summary.gam is probably just to compare fit indices like Log Likelihoods.
>
> If anyone has any other suggestions, though, please do let me know.
>
> Thanks so much.
>
> Will Shadish
>
> On 6/7/2013 3:02 PM, Gavin Simpson wrote:
>> 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"
>>>
>>>
>>
>