From: Gavin Simpson <gavin.simpson_at_ucl.ac.uk>

Date: Wed, 19 Nov 2008 17:27:47 +0000

Date: Wed, 19 Nov 2008 17:27:47 +0000

On Wed, 2008-11-19 at 09:33 -0500, David M Warner wrote:

> Greetings all

*> The help file for GAMM in mgcv indicates that the log likelihood for a
**> GAMM reported using
**>
**> summary(my.gamm$lme) (as an example) is not correct.
*

It is slightly more specific than that. It mentions that glmmPQL was used, and in your examples below, it isn't used as your GAMM's are Gaussian.

Having said that, the log likelihood for the Gaussian GAMM will be based on whatever you set argument 'method' to. It defaults to "ML" which is full Maximum Likelihood and IIRC is the correct form for checking the fixed effects in a likelihood ratio test. REML produces unbiased estimates of the variances, say of random effects, whereas ML doesn't.

If you are fitting Gaussian GAMMs then the book by Pinhiero and Bates that supports the nlme package (where the lme function comes from) will be very useful as the GAMM is just an LME. Simon Woods book on GAMs is perhaps the best place to find out about what gamm() is doing currently.

**HTH
**
G

*>
*

> However, in a past R-help post (included below), there is some indication

*> that the likelihood ratio test in anova.lme(mygamm$lme, mygamm1$lme) is
**> valid.
**>
**> How can I tell if anova.lme results are meaningful (are AIC, BIC, and
**> logLik estimates accurate)?
**>
**> The data include hydroacoustic estimates of fish biomass (lbloat) in 1,000
**> meter long intervals (elementary sampling units) from multiple transects
**> (each 20-30 km long, tranf) in two different lakes and three different
**> years.
**>
**> bloat.gamm1 <- gamm(lbloat ~ s(depth), correlation=corSpher(c(30000,
**> 0.01),form = ~ x+y|tranf, nugget=TRUE), data=fish3)
**>
**> bloat.gamm2 <- gamm(lbloat ~ lakef + s(depth),
**> correlation=corSpher(c(30000, 0.01),form = ~ x+y|tranf, nugget=TRUE),
**> data=fish3)
**>
**> bloat.gamm3 <- gamm(lbloat ~ lakef + s(depth, by=lakef),
**> correlation=corSpher(c(30000, 0.01),form = ~ x+y|tranf, nugget=TRUE),
**> data=fish3)
**>
**> > anova.lme(bloat.gamm1$lme, bloat.gamm2$lme, bloat.gamm3$lme)
**> Model df AIC BIC logLik Test L.Ratio
**> p-value
**> bloat.gamm1$lme 1 6 7916.315 7950.702 -3952.158
**> bloat.gamm2$lme 2 7 7902.718 7942.835 -3944.359 1 vs 2 15.597489
**> 0.0001
**> bloat.gamm3$lme 3 9 7910.987 7962.567 -3946.494 2 vs 3 4.269119
**> 0.1183
**>
**> Thanks
**> Dave
**>
**>
**>
**>
**>
**>
**> Hi R user,
**>
**> I am using the gamm() function of the mgcv-package. Now I would like to
**> decide on the random effects to include in the model. Within a GAMM
**> framework, is it allowed to compare the following two models
**>
**> inv_1<-gamm(y~te(sat,inv),data=daten_final, random=list(proband=~1))
**>
**> inv_2<-gamm(y~te(sat,inv),data=daten_final, random=list(proband=~sat))
**>
**>
**> with a likelihood ratio test for a traditional GLMM, like this:
**>
**> anova(inv_1$lme, inv_2$lme)
**>
**> The output is as follows:
**>
**> Model df AIC BIC logLik Test L.Ratio p-value
**> inv_2$lme 1 10 21495.90 21557.59 -10737.95
**> inv_1$lme 2 8 23211.12 23260.46 -11597.56 1 vs 2 1719.214 <.0001
**>
**>
**> Or is this not in tune with the automatic smoothing parameter selection
**> (i.e. it is not exactly the same for both model alternatives)?
**>
**>
**> If not, how can I decide on the selection of random effects?
**>
**>
**>
**>
**> This comparison is just as valid as it is for a regular linear mixed
**> model,
**> which is all that the GAMM is in this case --- the smoothing parameters
**> are
**> just variance components in your example.
**>
**> In general you have to be a bit careful with generalized likelihood ratio
**> tests involving variance components, of course, since the null hypothesis
**>
**> often involves restricting some variance parameters to the edge of their
**> possible range, which rather messes up the Taylor expansion about the null
**>
**> parameter values that underpins the large sample distributional results
**> used
**> in the glrt. Your example does involve such a problematic comparison, but
**> the
**> result is so clear cut here that there is not really any doubt that inv_2
**> is
**> better in this case (I wonder if inv_1 even passes basic model checking?).
**>
**> See Pinheiro and Bates (2000) for more info.
**>
**> hope this is some use,
**> Simon
**>
**>
**>
**>
**>
**>
**>
**>
**>
**>
**>
**>
**>
**> David Warner
**> Research Fishery Biologist
**> USGS Great Lakes Science Center
**> 1451 Green Road
**> Ann Arbor MI 48105
**> 734.214.9392
**> [[alternative HTML version deleted]]
**>
**> ______________________________________________
**> R-help_at_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.
*

-- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr. Gavin Simpson [t] +44 (0)20 7679 0522 ECRC, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% ______________________________________________ R-help_at_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.Received on Wed 19 Nov 2008 - 17:36:45 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Wed 19 Nov 2008 - 19:30:50 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*