[R] lme4 and lmeSplines

[Rod Ball]

The source of the problem is that lmeSplines requires a set of iid random 
effects corresponding to a transformed smoothing spline basis (`pdIdent' 
structure in nlme) i.e. pdIdent(~Zs -1),  where Zs is the transformed set of 
spline basis functions. Kevin's model does not work because each of his 
random effects has a separately estimated variance parameter, equivalent to a 
`pdDiag' structure in nlme. lmeSplines really does require nlme.

There is currently no way in lme4 to get the equivalent of a pdIdent 
structure, i.e. set of iid random effects with common variance parameter (see 
mesage from Doug Bates below), unless these correspond to a grouping factor 
as in (~1|row). There are a number of other nlme models e.g. genetics models, 
that cannot be fitted in lme4, examples are in my message to Doug Bates 
below.

Regards,
Rod Ball
-- 
Dr Roderick D. Ball, 
Statistician, 
Ensis Wood Quality - Linking Quality to Value
Ensis - The Joint Forces of CSIRO and SCION
Te Papa Tipu Innovation Park, 49 Sala Street
P.B. 3020, Rotorua, New Zealand
Phone: +64 7 343 5777
Facsimile: +64 7 348 0952
email: rod.ball at ensisjv.com
url: www.ensisjv.com

Re: pdIdent in lme4
Date: Mon, 7 Aug 2006 11:25:31 -0500
From: "Douglas Bates" <bates at stat.wisc.edu>
To: "Rod Ball" <rod.ball at ensisjv.com>
Cc: bates at wisc.edu, "Kevin Wright" <kwright68 at gmail.com>

Most of the time what would have been created using the pdIdent
structure in lme can be created from a model matrix in lmer.  However,
I don't know if this will apply to the spline basis functions.  I
imagine that what you want is to have a random effects for each of the
basis functions at each observation and then to externally impose the
structure of pdIdent on the variance-covariance matrix.  I'm afraid I
don't see a way of doing this in the current lmer structures.
However, I am documenting the structures much more extensively that
was done for the lme structures and also factoring the computation
into distinct steps.  Perhaps when I am done it will be possible for
someone to add in such a capability.

On 8/7/06, Rod Ball <rod.ball at ensisjv.com> wrote:
> Dear Doug,
>
> I have just being trying your lme4 package. The information I've read is 
that
> lme4 (lmer() function) is more general than lme. However I can't see how to
> do the equivalent of pdIdent(~Z -1) where Z is a matrix, or any pdIdent
> structure for that matter. Cf below or the help from my lmeSplines package.
> Type
> > library(lmeSplines)
> > ?smspline
> to see the example.
>
> I've had a query about how to do this (see below). The attempt below fails
> because adding multiple terms for each column of the Z-matrix gives a pdDiag
> structure, not a pdIdent structure, as needed for a smoothing spline. I 
can't
> see how to get a pdIdent structure at all, or how to use it. I did find a
> pdCompSymm class listed in lme4 on the web, but no indication how to use it,
> nor does it seem to be in my copy (R 2.3.1, just upgraded).
>
> > library(help=lme4)
>
>                 Information on package 'lme4'
>
> Description:
>
> Package:       lme4
> Version:       0.995-2
> Date:          2006-01-17
> Title:         Linear mixed-effects models using S4 classes
> Author:        Douglas Bates <bates at stat.wisc.edu> and Deepayan Sarkar
>                <deepayan at stat.wisc.edu>
> Maintainer:    Douglas Bates <bates at stat.wisc.edu>
> Description:   Fit linear and generalized linear mixed-effects models.
> Depends:       methods, R(>= 2.2.0), Matrix(>= 0.995-2), lattice
> Imports:       graphics, stats
> Suggests:      mlmRev
> SaveImage:     no
> LazyLoad:      yes
> License:       GPL version 2 or later
> Packaged:      Thu Jan 19 12:39:40 2006; bates
> Built:         R 2.3.1; ; 2006-08-02 20:15:09; unix
>
> Index:
>
> gsummary                Summarize a data frame by group
> lmList                  List of lm Objects with a Common Model
> lmList-class            Class "lmList"
> pooledSD                Extract pooled standard deviation
>
> Further information is available in the following vignettes in
> directory '/home/rod/R/library/lme4/doc':
>
> Implementation: Implementation Details (source, pdf)
>
> Another example is the heart rate data from your Rnews article. Suppose I 
want
> to examine if Drug effects and rates vary randomly with patients. How can I
> include for example Patient=pdIdent(~Drug -1)? Are there an lmer()
> equivalents for the following?
>
> fm4 <- update(fm1, random = list(Patient = ~ Time, Patient= pdIdent(~Drug
> -1)))
> or
> fm4a <- update(fm1, random=list(Patient = ~ Time, Patient= pdCompSymm(~Drug
> -1))
> or
> fm5 <- update(fm1, random=list(Patient = ~ Time, Patient= pdIdent(~Drug -1),
> Patient= pdIdent(~ Drug:Time-1)))
>
> Incidentally, I discovered also, when trying these examples, the intervals
> function does not seem to work in lme4, and the lattice package is loaded
> twice.
>
> > data("HR",package="SASmixed")
> > library(lme4)
> Loading required package: Matrix
> Loading required package: lattice
> Loading required package: lattice
> > fm1 <- lmer(HR ~ baseHR + Time*Drug + (1|Patient), data=HR)
> > intervals(fm1)
> Error: could not find function "intervals"
>
> Is lme4 intended to be a replacement for nlme, or only an alternative for
> certain models?
>
>
> Thanks,
> Rod Ball
> --
> Dr Roderick D. Ball,
> Statistician,
> Ensis Wood Quality - Linking Quality to Value
> Ensis - The Joint Forces of CSIRO and SCION
> Te Papa Tipu Innovation Park, 49 Sala Street
> P.B. 3020, Rotorua, New Zealand
> Phone: +64 7 343 5777
> Facsimile: +64 7 348 0952
> email: rod.ball at ensisjv.com
> url: www.ensisjv.com
>
>
> lme4 and lmeSplines
> Date: Wed, 2 Aug 2006 14:25:20 -0500
> From: "Kevin Wright" <kwright68 at gmail.com>
> To: r-help at stat.math.ethz.ch
>
> 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.
>
> Kevin Wright
>
>
>
>



On Thu, 10 Aug 2006 02:35, Spencer Graves wrote:
>    At least at one time, to use lmer after lme, you had to quit R
> before loading 'lme4' because of substantive conflicts between 'nlme'
> and 'lme4' / 'Matrix'.  That may be the source of your problem.  I can
> think of two possible ways to get around this:
>
> 	  1.  I might try quitting R and restarting after your use of
> 'lmeSplines' but before 'lmer'.
>
> 	  2.  If that failed, I might try making local copies of everything I
> needed from 'lmeSplines' and using them with 'lmer'.  If that failed,
> I'd use 'debug' to walk through the code (or local copies of whatever I
> needed, possibly with name changes) until I figured it out.
>
> 	  I know this is not what you wanted to hear, but I hope it helps.
> 	  Spencer Graves
>