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.simp...@uregina.ca 523 Research and Innovation Centre [tw] @ucfagls University of Regina Regina, SK S4S 0A2, Canada ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
