[R] mgcv, testing gamm vs lme, which degrees of freedom?
Carlo Fezzi
c.fezzi at uea.ac.uk
Wed Jun 16 21:33:24 CEST 2010
Dear all,
I am using the "mgcv" package by Simon Wood to estimate an additive mixed
model in which I assume normal distribution for the residuals. I would
like to test this model vs a standard parametric mixed model, such as the
ones which are possible to estimate with "lme".
Since the smoothing splines can be written as random effects, is it
correct to use an (approximate) likelihood ratio test for this? If so,
which is the correct number of degrees of freedom? Sometime the function
LogLik() seems to provide strange results regarding the number of degrees
of freedom (df) for the gam, for instance in the example I copied below
the df for the "gamm" are equal to the ones for the "lme", but the
summary(model.gam) seems to indicate a much higher edf for the gamm.
I would be very grateful to anybody who could point out a solution,
Best wishes,
Carlo
Example below:
----
rm(list = ls())
library(mgcv)
library(nlme)
set.seed(123)
x <- runif(100,1,10) # regressor
b0 <- rep(rnorm(10,mean=1,sd=2),each=10) # random intercept
id <- rep(1:10, each=10) # identifier
y <- b0 + x - 0.1 * x^3 + rnorm(100,0,1) # dependent variable
f1 <- lme(y ~ x + I(x^2), random = list(id=~1) , method="ML" ) # lme model
f2 <- gamm(y ~ s(x), random = list(id=~1), method="ML" ) # gamm
## same number of "df" according to logLik:
logLik(f1)
logLik(f2$lme)
## much higher edf according to summary:
summary(f2$gam)
-----------
More information about the R-help
mailing list