[Lme4-commits] r1768 - in pkg/lme4: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 21 09:02:40 CEST 2012
Author: bbolker
Date: 2012-06-21 09:02:39 +0200 (Thu, 21 Jun 2012)
New Revision: 1768
Modified:
pkg/lme4/R/lmer.R
pkg/lme4/R/predict.R
pkg/lme4/R/utilities.R
pkg/lme4/tests/lmer.R
pkg/lme4/tests/predict.R
Log:
fix predict with less than full set of previous levels
fix postVar of ranef when not all effects are requested
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-06-16 17:57:36 UTC (rev 1767)
+++ pkg/lme4/R/lmer.R 2012-06-21 07:02:39 UTC (rev 1768)
@@ -119,8 +119,19 @@
fr <- eval(mf, parent.frame())
# random effects and terms modules
reTrms <- mkReTrms(findbars(formula[[3]]), fr)
+ if (FALSE) {
+ ## BMB: I clearly don't know what's going on here yet.
+ ## test: lmer(angle ~ temp + recipe + (1 | replicate), data = cake)
+ ## test: sstudy9 <- subset(sleepstudy, Days == 1 | Days == 9)
+ ## m1 <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), data = sstudy9)
+ rankZ1 <- rankMatrix(bdiag(reTrms$Zt,Diagonal(ncol(reTrms$Zt))))
+ pZ1 <- nrow(reTrms$Zt)+ncol(reTrms$Zt)
+ if (rankZ1<pZ1)
+ stop(gettextf("rank of cBind(Z,1) = %d < ncol(Z)+1 = %d", rankZ1, pZ1+1))
+ }
if (any(unlist(lapply(reTrms$flist, nlevels)) >= nrow(fr)))
- stop("number of levels of each grouping factor must be less than number of obs")
+ 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
@@ -319,6 +330,15 @@
fr <- eval(mf, parent.frame())
# random-effects module
reTrms <- mkReTrms(findbars(formula[[3]]), fr)
+ if ((maxlevels <- max(unlist(lapply(reTrms$flist, nlevels)))) > nrow(fr))
+ stop("number of levels of each grouping factor must be",
+ "greater than or equal to number of obs")
+ ## FIXME: adjust test for families with estimated scale;
+ ## useSc is not defined yet/not defined properly?
+ ## if (useSc && maxlevels == nrow(fr))
+ ## stop("number of levels of each grouping factor must be",
+ ## "greater than number of obs")
+
## fixed-effects model matrix X - remove random parts from formula:
form <- formula
form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
@@ -1102,7 +1122,7 @@
if (postVar) {
sigsqr <- sigma(object)^2
vv <- .Call(merPredDcondVar, object at pp$ptr(), as.environment(rePos$new(object)))
- for (i in seq_along(ans))
+ for (i in names(ans)) ## seq_along(ans))
attr(ans[[i]], "postVar") <- vv[[i]] * sigsqr
}
if (drop)
Modified: pkg/lme4/R/predict.R
===================================================================
--- pkg/lme4/R/predict.R 2012-06-16 17:57:36 UTC (rev 1767)
+++ pkg/lme4/R/predict.R 2012-06-21 07:02:39 UTC (rev 1768)
@@ -69,8 +69,10 @@
re <- ranef(object)
##
ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
- new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],levels)
+ new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],
+ function(x) levels(droplevels(x)))
re_x <- mapply(function(x,n) {
+ ## find and deal with new levels
if (any(!new_levels[[n]] %in% rownames(x))) {
if (!allow.new.levels) stop("new levels detected in newdata")
## create an all-zero data frame corresponding to the new set of levels ...
@@ -80,6 +82,10 @@
newx[rownames(x),] <- x
x <- newx
}
+ ## find and deal with missing old levels
+ if (any(!rownames(x) %in% new_levels[[n]])) {
+ x <- x[rownames(x) %in% new_levels[[n]],,drop=FALSE]
+ }
x
},
re,names(re),SIMPLIFY=FALSE)
Modified: pkg/lme4/R/utilities.R
===================================================================
--- pkg/lme4/R/utilities.R 2012-06-16 17:57:36 UTC (rev 1767)
+++ pkg/lme4/R/utilities.R 2012-06-21 07:02:39 UTC (rev 1768)
@@ -506,6 +506,7 @@
nth=length(pp$theta), q=nrow(pp$Zt),
nAGQ=rho$nAGQ,
compDev=rho$compDev,
+ ## FIXME: need to consider GLM families with estimated scale parameter!
useSc=(rcl != "glmResp"),
reTrms=length(reTrms$cnms),
spFe=0L,
Modified: pkg/lme4/tests/lmer.R
===================================================================
--- pkg/lme4/tests/lmer.R 2012-06-16 17:57:36 UTC (rev 1767)
+++ pkg/lme4/tests/lmer.R 2012-06-21 07:02:39 UTC (rev 1768)
@@ -20,8 +20,16 @@
sstudy9 <- subset(sleepstudy, Days == 1 | Days == 9)
try({## This "did work" in lme4.0 and nlme -- FIXME ??
m1 <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), data = sstudy9)
+ if (FALSE) {
+ image(reTrms$Zt)
+ testRankZt <- function(object) {
+ rankMatrix(cBind(reTrms$Zt,1))<nrow(reTrms$Zt)
+ }
+ }
## -> Error in ptr() : Downdated VtV is not positive definite
## FIXME?(2): More helpful error message, or rank-checking diagnostics?
+ ## how would we check the rank here?
+
print(sm1 <- summary(m1))
fm1 <- fitted(m1)
})
Modified: pkg/lme4/tests/predict.R
===================================================================
--- pkg/lme4/tests/predict.R 2012-06-16 17:57:36 UTC (rev 1767)
+++ pkg/lme4/tests/predict.R 2012-06-21 07:02:39 UTC (rev 1768)
@@ -77,5 +77,22 @@
tmpf(sleepstudy,REform=~(0+Days|Subject))
tmpf(sleepstudy,REform=~(1|Subject))
+## from 'Colonel Triq'
+m <- lmer(angle ~ temp + recipe + (1 | replicate), data=cake)
+summary(m)
+# replicate 1 only appears in rows 1:18.
+rownames(cake[cake$replicate==1,])
+## Prediction using "new" data frame
+## that includes no replicate 1 examples will fail upon prediction.
+## Including at least one row with replicate=1 will allow predictions.
+## REform parameter is based on syntax found in example(predict.merMod)
+
+predict(m, newdata=cake[-1:-17,], REform=~ (1 | replicate)) # works
+predict(m, newdata=cake[-1:-18,], REform=NA) # works
+predict(m, newdata=cake[-1:-18,], REform=~ (1 | replicate)) # doesn't work
+predict(m, newdata=cake[-1:-18,], REform=~ (1 | replicate), allow.new.levels=TRUE) # doesn't work
+
+
+
More information about the Lme4-commits
mailing list