[Lme4-commits] r1773 - in pkg/lme4: R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 23 13:27:13 CEST 2012
Author: bbolker
Date: 2012-06-23 13:27:12 +0200 (Sat, 23 Jun 2012)
New Revision: 1773
Modified:
pkg/lme4/R/predict.R
pkg/lme4/R/profile.R
pkg/lme4/tests/predict.R
Log:
fixes to predict.merMod; added/modified tests
substitude [["par"]] for $par in some profile code; added globalVariables(".par") to suppress check warnings
Modified: pkg/lme4/R/predict.R
===================================================================
--- pkg/lme4/R/predict.R 2012-06-22 19:11:41 UTC (rev 1772)
+++ pkg/lme4/R/predict.R 2012-06-23 11:27:12 UTC (rev 1773)
@@ -70,8 +70,10 @@
##
ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],
- function(x) levels(droplevels(x)))
- re_x <- mapply(function(x,n) {
+ function(x) levels(droplevels(factor(x))))
+ ## FIXME: should this be unique(as.character(x)) instead?
+ ## (i.e., what is the proper way to protect against non-factors?)
+ levelfun <-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")
@@ -87,24 +89,30 @@
x <- x[rownames(x) %in% new_levels[[n]],,drop=FALSE]
}
x
- },
- re,names(re),SIMPLIFY=FALSE)
+ }
+ ## fill in/delete levels as appropriate
+ re_x <- mapply(levelfun,re,names(re),SIMPLIFY=FALSE)
## separate random effects from orig model into individual columns
- re_List <- do.call(c,lapply(re_x,as.list))
- re_names <- names(re_List)
- z_names <- mapply(paste,names(ReTrms$cnms),ReTrms$cnms,MoreArgs=list(sep="."))
+ ## re_List <- do.call(c,lapply(re_x,as.list))
+ ## re_names <- names(re_List)
+ ## names corresponding to random effects specified in 'predict' call
+ ## z_names <- mapply(paste,names(ReTrms$cnms),ReTrms$cnms,MoreArgs=list(sep="."))
+ re_new <- list()
+ if (any(!names(ReTrms$cnms) %in% names(re)))
+ stop("grouping factors specified in REform that were not present in original model")
## pick out random effects values that correspond to
## random effects incorporated in REform ...
- ## FIXME: more tests for possible things going wrong here?
- m <- match(z_names,re_names)
- if (any(is.na(m)))
- stop("random effects specified in REform that were not present in original model")
- re_new <- unlist(re_List[m])
- pred <- pred + drop(as.matrix(re_new %*% ReTrms$Zt))
+ for (i in seq_along(ReTrms$cnms)) {
+ rname <- names(ReTrms$cnms)[i]
+ if (any(!ReTrms$cnms[[rname]] %in% names(re[[rname]])))
+ stop("random effects specified in REform that were not present in original model")
+ re_new[[i]] <- re_x[[rname]][,ReTrms$cnms[[rname]]]
+ }
+ re_newvec <- unlist(lapply(re_new,t)) ## must TRANSPOSE RE matrices before unlisting
+ pred <- pred + drop(as.matrix(re_newvec %*% ReTrms$Zt))
} ## REform provided
} ## predictions with new data or new REform
- ## FIXME: would like to have an isGLMM() accessor for merMod objects?
- if (is(object at resp,"glmResp") && type=="response") {
+ if (isGLMM(object) && type=="response") {
pred <- object at resp$family$linkinv(pred)
}
return(pred)
Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R 2012-06-22 19:11:41 UTC (rev 1772)
+++ pkg/lme4/R/profile.R 2012-06-23 11:27:12 UTC (rev 1773)
@@ -3,6 +3,10 @@
##' Methods for function \code{\link{profile}} (package \pkg{stats}), here for
##' profiling (fitted) mixed effect models.
##'
+
+globalVariables(".par",package="lme4")
+
+##'
##'
##' @name profile-methods
##' @aliases profile-methods profile.merMod
@@ -800,7 +804,7 @@
sigs <- grep("^\\.sig", cn)
if (length(sigs)) {
colnames(x) <- sub("^\\.sig", ".lsig", cn)
- levels(x$.par) <- sub("^\\.sig", ".lsig", levels(x$.par))
+ levels(x[[".par"]]) <- sub("^\\.sig", ".lsig", levels(x[[".par"]]))
names(attr(x, "backward")) <-
names(attr(x, "forward")) <-
sub("^\\.sig", ".lsig", names(attr(x, "forward")))
@@ -830,7 +834,7 @@
sigs <- grep("^\\.sig", cn)
if (length(sigs)) {
colnames(x) <- sub("^\\.sig", ".sigsq", cn)
- levels(x$.par) <- sub("^\\.sig", ".sigsq", levels(x$.par))
+ levels(x[[".par"]]) <- sub("^\\.sig", ".sigsq", levels(x[[".par"]]))
names(attr(x, "backward")) <-
names(attr(x, "forward")) <-
sub("^\\.sig", ".sigsq", names(attr(x, "forward")))
@@ -914,7 +918,7 @@
onms <- names(spl) # names of original variables
vc <- onms[grep("^.sig", onms)] # variance components
ans <- subset(pr, .par %in% vc, select=c(".zeta", vc, ".par"))
- ans$.par <- factor(ans$.par) # drop unused levels
+ ans[[".par"]] <- factor(ans[[".par"]]) # drop unused levels
if (".lsig" %in% vc) ans$.lsig <- exp(ans$.lsig)
attr(ans, "forward") <- attr(ans, "backward") <- list()
for (nm in vc) {
Modified: pkg/lme4/tests/predict.R
===================================================================
--- pkg/lme4/tests/predict.R 2012-06-22 19:11:41 UTC (rev 1772)
+++ pkg/lme4/tests/predict.R 2012-06-23 11:27:12 UTC (rev 1773)
@@ -62,11 +62,35 @@
## linear model, so results should be identical patterns but smaller --
## not including intermediate days
newdata <- with(sleepstudy,expand.grid(Days=range(Days),Subject=unique(Subject)))
-p2 <- predict(fm2,newdata)
-p3 <- predict(fm2,newdata,REform=NA)
-p3 <- predict(fm2,newdata,REform=~(0+Days|Subject))
-p4 <- predict(fm2,newdata,REform=~(1|Subject))
-par(mfrow=c(2,2))
+newdata$p2 <- predict(fm2,newdata)
+newdata$p3 <- predict(fm2,newdata,REform=NA)
+newdata$p4 <- predict(fm2,newdata,REform=~(0+Days|Subject))
+newdata$p5 <- predict(fm2,newdata,REform=~(1|Subject))
+
+## reference values from an apparently-working run
+refval <- structure(list(Days = c(0, 9, 0, 9, 0, 9), Subject = structure(c(1L,
+1L, 2L, 2L, 3L, 3L), .Label = c("308", "309", "310", "330", "331",
+"332", "333", "334", "335", "337", "349", "350", "351", "352",
+"369", "370", "371", "372"), class = "factor"), p2 = c(253.663652396798,
+430.66001930835, 211.006415533628, 227.634788908917, 212.444742696829,
+257.61053840953), p3 = c(251.405104848485, 345.610678484848,
+251.405104848485, 345.610678484848, 251.405104848485, 345.610678484848
+), p4 = c(251.405104848485, 428.401471760037, 251.405104848485,
+268.033478223774, 251.405104848485, 296.570900561186), p5 = c(253.663652396798,
+347.869226033161, 211.006415533628, 305.211989169991, 212.444742696829,
+306.650316333193)), .Names = c("Days", "Subject", "p2", "p3",
+"p4", "p5"), out.attrs = structure(list(dim = structure(c(2L,
+18L), .Names = c("Days", "Subject")), dimnames = structure(list(
+ Days = c("Days=0", "Days=9"), Subject = c("Subject=308",
+ "Subject=309", "Subject=310", "Subject=330", "Subject=331",
+ "Subject=332", "Subject=333", "Subject=334", "Subject=335",
+ "Subject=337", "Subject=349", "Subject=350", "Subject=351",
+ "Subject=352", "Subject=369", "Subject=370", "Subject=371",
+ "Subject=372")), .Names = c("Days", "Subject"))), .Names = c("dim",
+"dimnames")), row.names = c(NA, 6L), class = "data.frame")
+
+stopifnot(all.equal(head(newdata),refval))
+
library(lattice)
tmpf <- function(data,...) {
data$Reaction <- predict(fm2,data,...)
@@ -77,22 +101,17 @@
tmpf(sleepstudy,REform=~(0+Days|Subject))
tmpf(sleepstudy,REform=~(1|Subject))
-## from 'Colonel Triq'
+## from 'Colonel Triq': examples using *fewer* random effect levels
+## than in original data set
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))
+predict(m, newdata=cake[-1:-18,], REform=NA)
+predict(m, newdata=cake[-1:-18,], REform=~ (1 | replicate))
+predict(m, newdata=cake[-1:-18,], REform=~ (1 | replicate), allow.new.levels=TRUE)
-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