[Lme4-commits] r1783 - in pkg/lme4: R tests

[noreply at r-forge.r-project.org]

Author: bbolker
Date: 2012-06-27 19:01:45 +0200 (Wed, 27 Jun 2012)
New Revision: 1783

Added:
   pkg/lme4/tests/predict_basis.R
Modified:
   pkg/lme4/R/lmer.R
   pkg/lme4/R/predict.R
   pkg/lme4/tests/offset.R
Log:

  update predict, terms to allow prediction with data-dependent bases
(spline, orthogonal polynomial, etc.); add relevant tests

  tweaks to offset tests (remove some #boom! comments that are no longer valid)



Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-06-27 09:48:31 UTC (rev 1782)
+++ pkg/lme4/R/lmer.R	2012-06-27 17:01:45 UTC (rev 1783)
@@ -139,9 +139,13 @@
         stop("number of levels of each grouping factor must be ",
              "less than number of obs")
     ## fixed-effects model matrix X - remove random effects from formula:
-    form <- formula
-    form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
-    X <- model.matrix(form, fr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet
+    fixedform <- formula
+    fixedform[[3]] <- if(is.null(nb <- nobars(fixedform[[3]]))) 1 else nb
+    mf$formula <- fixedform
+    ## re-evaluate model frame to extract predvars component
+    fixedfr <- eval(mf, parent.frame())
+    attr(attr(fr,"terms"),"predvars.fixed") <- attr(attr(fixedfr,"terms"),"predvars")
+    X <- model.matrix(fixedform, fixedfr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet
     p <- ncol(X)
     if ((qrX <- qr(X))$rank < p)
         stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
@@ -1444,8 +1448,11 @@
 ##' @S3method terms merMod
 terms.merMod <- function(x, fixed.only=TRUE, ...) {
   if (fixed.only) {
-      terms.formula(formula(x,fixed.only=TRUE))
-  } else attr(x at frame,"terms")
+      tt <- terms.formula(formula(x,fixed.only=TRUE))
+      attr(tt,"predvars") <- attr(attr(x at frame,"terms"),"predvars.fixed")
+      tt
+  }
+  else attr(x at frame,"terms")
 }
 
 ##' @importFrom stats update

Modified: pkg/lme4/R/predict.R
===================================================================
--- pkg/lme4/R/predict.R	2012-06-27 09:48:31 UTC (rev 1782)
+++ pkg/lme4/R/predict.R	2012-06-27 17:01:45 UTC (rev 1783)
@@ -15,7 +15,6 @@
 ##' @param ... optional additional parameters.  None are used at present.
 ##' @return a numeric vector of predicted values
 ##' @note explain why there is no option for computing standard errors of predictions!
-##' @note offsets not yet handled
 ##' @examples
 ##' (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 |herd), cbpp, binomial))
 ##' str(p0 <- predict(gm1))            # fitted values
@@ -43,10 +42,10 @@
         if (is.null(newdata)) {
             X <- X_orig
         } else {
-            form <- form_orig
-            form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
-            RHS <- form[-2]
-            X <- model.matrix(RHS, newdata, contrasts.arg=attr(X_orig,"contrasts"))
+            RHS <- formula(object,fixed.only=TRUE)[-2]
+            Terms <- terms(object,fixed.only=TRUE)
+            X <- model.matrix(RHS, model.frame(delete.response(Terms), newdata),
+                              contrasts.arg=attr(X_orig,"contrasts"))
         }
         pred <- drop(X %*% fixef(object))
         ## modified from predict.glm ...

Modified: pkg/lme4/tests/offset.R
===================================================================
--- pkg/lme4/tests/offset.R	2012-06-27 09:48:31 UTC (rev 1782)
+++ pkg/lme4/tests/offset.R	2012-06-27 17:01:45 UTC (rev 1783)
@@ -22,7 +22,7 @@
 
 p0 <- predict(fm1)
 p1 <- predict(fm1,newdata=d)
-p2 <- predict(fm1off,newdata=d) ## boom
+p2 <- predict(fm1off,newdata=d)
 stopifnot(all.equal(p0,unname(p1),unname(p2)))
 
 
@@ -31,14 +31,14 @@
 d$y <- rpois(nrow(d),d$mu)
 
 gm1 <- glmer(y~x+(1|block),data=d,family=poisson)
-gm1off <- lmer(y~x+(1|block)+offset(3*x),data=d,family=poisson)
+gm1off <- glmer(y~x+(1|block)+offset(3*x),data=d,family=poisson)
 
 ## check equality
 stopifnot(all.equal(fixef(gm1)[2]-3,fixef(gm1off)[2]))
 
 p0 <- predict(gm1)
 p1 <- predict(gm1,newdata=d)
-p2 <- predict(gm1off,newdata=d) ## boom
+p2 <- predict(gm1off,newdata=d)
 stopifnot(all.equal(p0,unname(p1),unname(p2)))
 
 ## FIXME: should also test simulations

Added: pkg/lme4/tests/predict_basis.R
===================================================================
--- pkg/lme4/tests/predict_basis.R	                        (rev 0)
+++ pkg/lme4/tests/predict_basis.R	2012-06-27 17:01:45 UTC (rev 1783)
@@ -0,0 +1,39 @@
+## test for models containing data-defined bases
+
+## ?makepredictcall
+## ?model.frame
+## ????
+
+data(sleepstudy,package="lme4")
+library(splines)
+
+## lm0 <- lm(Reaction~ns(Days,2),sleepstudy)
+## attr(terms(lm0),"predvars")
+## library(nlme)
+## lme1 <- lme(Reaction~ns(Days,2),random=~1|Subject,sleepstudy)
+## attr(terms(lme1),"predvars")  ## no!
+## attr(lme1$terms,"predvars")   ## yes
+## detach("package:nlme")
+
+library(lme4)
+library(testthat)
+fm1 <- lmer(Reaction ~ ns(Days,2) + (1|Subject), sleepstudy)
+fm2 <- lmer(Reaction ~ poly(Days,2) + (1|Subject), sleepstudy)
+fm3 <- lmer(Reaction ~ poly(Days,2,raw=TRUE) + (1|Subject), sleepstudy)
+
+
+newdat0 <- data.frame(Days=unique(sleepstudy$Days))
+newdat <- data.frame(Days=5:12)
+tmpf <- function(fit) {
+    with(sleepstudy,plot(Reaction~Days,xlim=c(0,12)))
+    with(sleepstudy,points(Days,predict(fit),col=2))
+    with(newdat0,lines(Days,predict(fit,REform=NA,newdata=newdat0),col=4))
+    with(newdat,lines(Days,predict(fit,REform=NA,newdata=newdat),col=5))
+}
+
+## pictures
+tmpf(fm1)
+tmpf(fm2)
+tmpf(fm3)
+stopifnot(all.equal(predict(fm2,newdat,REform=NA),
+                    predict(fm3,newdat,REform=NA)))