[R-sig-ME] inconsistent gamm likelihoods

[Alan Swanson]

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