I'm getting some odd results computing log-likelihoods with gls using splines with increasing degrees of freedom -- the deviance *increases* substantially with increasing df. (Since spline models with increasing df aren't nested, it need not decline monotonically but I would expect it to have a decreasing trend!)
I may just be confused, but I *think* the issue is somewhere within the log-likelihood calculation for gls. The example below compares a series of spline fits to a simulated exponentially decreasing data set with some heteroscedasticity; a more detailed example comparing fits of lm() and gls() to a simpler simulation [linear regression] shows similar, although less extreme, behavior, and is posted at http://www.zoo.ufl.edu/bolker/splinelik.Rnw (or http://www.zoo.ufl.edu/bolker/splinelik.pdf ). Apologies for boneheadedness, but I have taken this as far as I can for now ... ------------- ## simulate "data" (exponential decrease w/ heterosced.) set.seed(1001) n=1000 x = sort(runif(n)) grow_det = exp(-2*x) grow_var = 0.1*grow_det^2 y = rnorm(n,mean=grow_det,sd=sqrt(grow_var)) dat = data.frame(x=x,y=y) ## nlme likes to have data= specified ## fit true model library(nlme) g1 = gnls(y~a*exp(-b*x), start=list(a=1,b=2), weights=varPower(form=~fitted(.)), data=dat) expdev = -2*logLik(g1) ## Fitting the true model recovers ## the true parameters nicely: coef(g1) ## Fit a series of splines: sfit <- function(d) { form <- bquote(y~ns(x,df=.(d))) gls(eval(form), weights=varPower(form=~fitted(.)), data=dat) } spline_df = c(3:15,20,25,40) spline_list = lapply(as.list(spline_df),sfit) ## calculate deviances: spline_dev = -2*sapply(spline_list,logLik) plot(spline_df,spline_dev,type="b", ylim=range(c(expdev,spline_dev))) abline(h=expdev,col=2) ______________________________________________ R-help@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.
