[R-sig-ME] inconsistent gamm likelihoods
Alan Swanson
as146608 at umconnect.umt.edu
Wed Jun 23 01:16:28 CEST 2010
I have some soil measurements along transects and I'd like to test for
a treatment effect. GAMM models with a random intercept fit the data
well, although there is some residual autocorrelation. My approach is
to use a likelihood ratio test to compare GAMMs with no treatment effect
to GAMMs that allows for a separate fit by treatment.
My problem is that I'm getting very different results when I use
penalized cubic regression splines versus cr splines with fixed degrees
of freedom. For some response variables, the results are as I would
expect: likelihood values are about the same when effective degrees of
freedom are similar. For other response variables, penalized cr splines
give inferior likelihoods when the treatment effect is included, so
likelihood ratio test statistic is much smaller (4 vs 30), even though
the fits look very similar.
Does anybody know why penalized cr splines would give a significantly
different likelihood than a very similar fit with fixed df? Sample code
(same syntax but I am unable to reproduce the inconsistency) and session
info are pasted below.
Many thanks in advance,
Alan Swanson
#simulate data#
d_vec <- -25:25; n <- length(d_vec)
simdata <-
data.frame(treatment=factor(rep(1:3,rep(n*3,3))),dist=rep(d_vec,9),rep=factor(rep(1:9,n)))
t1_mean <- .002*d_vec^2+.00005*d_vec^3
t2_mean <- .002*d_vec^2+.00007*d_vec^3+1
t3_mean <- .002*d_vec^2+.00009*d_vec^3+2
plot(t3_mean~d_vec,type="l");lines(t2_mean~d_vec,lty=2);lines(t1_mean~d_vec,lty=3)
simdata$response <-
c(rep(t1_mean,3),rep(t2_mean,3),rep(t3_mean,3))+rep(rnorm(9,0,10),rep(n,9))+rnorm(nrow(simdata),0,20)
#fit gamm models with fixed df#
require(mgcv)
m0_fixed<-gamm(response~s(dist, fx = T, k = 4, bs =
"cr"),random=list(rep=~1),data=simdata)
m1_fixed<-gamm(response~treatment + s(dist,by=treatment,fx = T, k = 4,
bs = "cr"),random=list(rep=~1),data=simdata)
#fit gamm models with penalized cubic regression splines.
m0_pen<-gamm(response~s(dist, fx = F, k = -1, bs =
"cr"),random=list(rep=~1),data=simdata)
m1_pen<-gamm(response~treatment + s(dist,by=treatment,fx = F, k = -1, bs
= "cr"),random=list(rep=~1),data=simdata)
#likelihood ratio tests
print(anova(m0_fixed$lme,m1_fixed$lme))
print(anova(m0_pen$lme,m1_pen$lme))
> sessionInfo()
R version 2.11.1 (2010-05-31)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] mgcv_1.6-2 nlme_3.1-96 RWinEdt_1.8-2
loaded via a namespace (and not attached):
[1] grid_2.11.1 lattice_0.18-8 Matrix_0.999375-40
More information about the R-sig-mixed-models
mailing list