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