[R] logLik calculation in gls (nlme)

[Ben Bolker]

   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)