[R] lme4 and lmeSplines

[Douglas Bates]

On 8/2/06, Kevin Wright <kwright68 at gmail.com> wrote:
> I'm trying to use the lmeSplines package together with lme4.
>
> Below is (1) an example of lmeSplines together with nlme (2) an
> attempt to use lmeSplines with lme4 (3) then a comparison of the
> random effects from the two different methods.
>
> (1)
>
> require(lmeSplines)
> data(smSplineEx1)
> dat <- smSplineEx1
> dat.lo <- loess(y~time, data=dat)
> plot(dat.lo)
> dat$all <- rep(1,nrow(dat))
> times20 <- seq(1,100,length=20)
> Zt20 <- smspline(times20)
> dat$Zt20 <- approx.Z(Zt20, times20, dat$time)
> fit1.20 <- lme(y~time, data=dat, random=list(all=pdIdent(~Zt20-1)))
> # Loess model
> dat.lo <- loess(y~time, data=dat)
> plot(dat.lo)
> # Spline model
> with(dat, lines(fitted(fit1.20)~time, col="red"))
> # Save random effects for later
> ranef.nlme <- unlist(ranef(fit1.20))
>
> (2) Now an attempt to use lme4:
>
> library(lmeSplines)
> detach(package:nlme)
> library(lme4)
> data(smSplineEx1)
> # Use 20 spline in lme4
> dat <- smSplineEx1
> times20 <- seq(1,100,length=20)
> Zt20 <- smspline(times20)
> dat <- cbind(dat, approx.Z(Zt20, times20, dat$time))
> names(dat)[4:21] <- paste("Zt",names(dat)[4:21],sep="")
> dat$all <- rep(1, nrow(dat))
> fit1.20 <- lmer(y~time
>              +(-1+Zt1|all)+(-1+Zt2|all)+(-1+Zt3|all)+(-1+Zt4|all)+(-1+Zt5|all)+(-1+Zt6|all)
>              +(-1+Zt7|all)+(-1+Zt8|all)+(-1+Zt9|all)+(-1+Zt10|all)+(-1+Zt11|all)+(-1+Zt12|all)
>              +(-1+Zt13|all)+(-1+Zt14|all)+(-1+Zt15|all)+(-1+Zt16|all)+(-1+Zt17|all)+(-1+Zt18|all),
>              data=dat)
> #summary(fit1)
> # Plot the data and loess fit
> dat.lo <- loess(y~time, data=dat)
> plot(dat.lo)
> # Fitting with splines
> with(dat, lines(fitted(fit1.20)~time, col="red"))
> ranef.lme4 <- unlist(ranef(fit1.20))
>
> (3) Compare nlme lme4 random effects
>
> plot(ranef.nlme~ranef.lme4)
>
> The plot of fitted values from lme4 is visually appealing, but the
> random effects from lme4 are peculiar--three are non-zero and the rest
> are essentially zero.
>
> Any help in getting lme4 + lmeSplines working would be appreciated.
> It is not unlikely that I have the lmer syntax wrong.

It is not surprising that you get different answers from lme and lme4
because you are fitting different models.  The variances of the random
effects for the spline basis in the lme fit are constrained to be
equal.  In the lmer fit they are not constrained to be equal.  It is
interesting that you get all but three of the variances essentially
zero.  That means that there are only three active components in your
spline basis, out of 20, for the fit.

I exchanged some mail off-list with Rod Ball about the definition of
the random effects needed for lmeSplines and we concluded that the
current capabilities in lmer are not sufficiently flexible to use it
for lmeSplines.  However, the sources for lmer are freely available
and any enterprising programmer who would like to use the components
for a more flexible model is welcome to do so.

The point is that the tools are available in lmer to represent a
mixed-effects model and to evaluate the log-likelihood or restricted
log-likelihood from such a model very efficiently.  To optimize a
model such as that being used in the lme fit one needs to go from the
parameters to the model representation.  This is the part that would
need to be written and after that you could hook into existing code.

> Kevin Wright
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>