[R] lme4 and lmeSplines
Rod Ball
rod.ball at ensisjv.com
Sun Aug 13 23:59:53 CEST 2006
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
>
More information about the R-help
mailing list