gamm in mgcv random effect significance

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"