[R] gamm in mgcv random effect significance

[Gavin Simpson]

On Tue, 2013-06-11 at 10:08 -0700, 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

This is with mgcv:::gamm yes? Strange - did you load nlme first? I think
mgcv now imports from the nlme package so to use its functions/methods
explicitly you need to load nlme - before loading mgcv also loaded nlme,
but it no longer does so.

That *should* work.

HTH

G

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

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