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
inconsistent gamm likelihoods
1 message · Alan Swanson