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
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 Kevin Wright 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. > > 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.
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. >