[Lme4-commits] r1583 - in pkg/lme4Eigen: . R tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 8 23:40:21 CET 2012
Author: bbolker
Date: 2012-02-08 23:40:20 +0100 (Wed, 08 Feb 2012)
New Revision: 1583
Added:
pkg/lme4Eigen/R/predict.R
pkg/lme4Eigen/tests/predict.R
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/tests/simulate.R
Log:
add predict.merMod
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-02-08 21:21:15 UTC (rev 1582)
+++ pkg/lme4Eigen/DESCRIPTION 2012-02-08 22:40:20 UTC (rev 1583)
@@ -47,3 +47,4 @@
'sparsegrid.R'
'utilities.R'
'optimizer.R'
+ 'predict.R'
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-02-08 21:21:15 UTC (rev 1582)
+++ pkg/lme4Eigen/NAMESPACE 2012-02-08 22:40:20 UTC (rev 1583)
@@ -113,6 +113,7 @@
S3method(plot,coef.mer)
S3method(plot,lmList.confint)
S3method(plot,ranef.mer)
+S3method(predict,merMod)
S3method(print,merMod)
S3method(print,summary.mer)
S3method(print,VarCorr.merMod)
Added: pkg/lme4Eigen/R/predict.R
===================================================================
--- pkg/lme4Eigen/R/predict.R (rev 0)
+++ pkg/lme4Eigen/R/predict.R 2012-02-08 22:40:20 UTC (rev 1583)
@@ -0,0 +1,77 @@
+##' @param newdata data frame for which to evaluate predictions
+##' @param REform formula for random effects to include. If NULL, include all random effects; if NA, include no random effects
+##' @param allow.new.levels (logical) if FALSE, then any new levels detected in \code{newdata} will trigger an error; if TRUE, then the prediction will use the unconditional (population-level) values for data with previously unobserved levels
+##' @note explain why there is no option for computing standard errors of predictions!
+##' @note offsets not yet handled
+
+## FIXME: take some of the stuff from tests/predict.R and incorporate it as examples
+
+predict.merMod <- function(object, newdata=NULL, REform=NULL,
+ terms=NULL, type=c("link","response"),
+ allow.new.levels=FALSE, ...) {
+ ## FIXME: appropriate names for result vector?
+ if (any(getME(object,"offset")!=0)) stop("offsets not handled yet") ## FIXME
+ type <- match.arg(type)
+ if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
+ X_orig <- getME(object, "X")
+ ## FIXME/WARNING: how do we do this in an eval-safe way???
+ form_orig <- eval(object at call$formula,parent.frame())
+ if (is.null(newdata) && is.null(REform)) {
+ if (is(object at resp,"lmerResp")) return(fitted(object))
+ return(switch(type,response=object at resp$mu, ## fitted(object),
+ link=object at resp$eta)) ## fixme: getME() ?
+ } else {
+ 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"))
+ }
+ pred <- drop(X %*% fixef(object))
+ if (is.null(REform)) {
+ REform <- form_orig[-2]
+ }
+ ## FIXME: ??? can't apply is.na() to a 'language' object?
+ ## what's the appropriate test?
+ if (is.language(REform)) {
+ re <- ranef(object)
+ ##
+ ReTrms <- mkReTrms(findbars(REform[[2]]),newdata)
+ new_levels <- lapply(newdata[unique(sort(names(ReTrms$cnms)))],levels)
+ re_x <- mapply(function(x,n) {
+ 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 ...
+ newx <- as.data.frame(matrix(0,nrow=length(new_levels[[n]]),ncol=ncol(x),
+ dimnames=list(new_levels[[n]],names(x))))
+ ## then paste in the matching RE values from the original fit/set of levels
+ newx[rownames(x),] <- x
+ x <- newx
+ }
+ x
+ },
+ 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="."))
+ ## 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))
+ } ## 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") {
+ pred <- object at resp$family$linkinv(pred)
+ }
+ return(pred)
+}
+
+
Added: pkg/lme4Eigen/tests/predict.R
===================================================================
--- pkg/lme4Eigen/tests/predict.R (rev 0)
+++ pkg/lme4Eigen/tests/predict.R 2012-02-08 22:40:20 UTC (rev 1583)
@@ -0,0 +1,81 @@
+library(lme4Eigen)
+gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ data = cbpp, family = binomial)
+
+## fitted values
+p0 <- predict(gm1)
+## fitted values, unconditional (level-0)
+p1 <- predict(gm1,REform=NA)
+stopifnot(length(unique(p1))==length(unique(cbpp$period)))
+
+matplot(cbind(p0,p1),col=1:2,type="b")
+newdata <- with(cbpp,expand.grid(period=unique(period),herd=unique(herd)))
+## new data, all RE
+p2 <- predict(gm1,newdata)
+## new data, level-0
+p3 <- predict(gm1,newdata,REform=NA)
+## explicitly specify RE
+p4 <- predict(gm1,newdata,REform=~(1|herd))
+stopifnot(all.equal(p2,p4))
+
+
+p5 <- predict(gm1,type="response")
+stopifnot(all.equal(p5,plogis(p0)))
+
+matplot(cbind(p2,p3),col=1:2,type="b")
+
+## effects of new RE levels
+newdata2 <- rbind(newdata,
+ data.frame(period=as.character(1:4),herd=rep("new",4)))
+try(predict(gm1,newdata2))
+p6 <- predict(gm1,newdata2,allow.new.levels=TRUE)
+stopifnot(all.equal(p2,p6[1:length(p2)])) ## original values should match
+## last 4 values should match unconditional values
+stopifnot(all(tail(p6,4)==predict(gm1,newdata=data.frame(period=factor(1:4)),REform=NA)))
+
+## multi-group model
+fm1 <- lmer(diameter ~ (1|plate) + (1|sample), Penicillin)
+## fitted values
+p0 <- predict(fm1)
+## fitted values, unconditional (level-0)
+p1 <- predict(fm1,REform=NA)
+
+matplot(cbind(p0,p1),col=1:2,type="b")
+newdata <- with(Penicillin,expand.grid(plate=unique(plate),sample=unique(sample)))
+## new data, all RE
+p2 <- predict(fm1,newdata)
+## new data, level-0
+p3 <- predict(fm1,newdata,REform=NA)
+## explicitly specify RE
+p4 <- predict(fm1,newdata,REform=~(1|plate)+(~1|sample))
+p4B <- predict(fm1,newdata,REform=~(1|sample)+(~1|plate))
+stopifnot(all.equal(p2,p4,p4B))
+
+p5 <- predict(fm1,newdata,REform=~(1|sample))
+p6 <- predict(fm1,newdata,REform=~(1|plate))
+
+matplot(cbind(p2,p3,p5,p6),type="b",lty=1,pch=16)
+
+fm2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+p0 <- predict(fm2)
+p1 <- predict(fm2,REform=NA)
+## 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))
+library(lattice)
+tmpf <- function(data,...) {
+ data$Reaction <- predict(fm2,data,...)
+ xyplot(Reaction~Days,group=Subject,data=data,type="l")
+}
+tmpf(sleepstudy)
+tmpf(sleepstudy,REform=NA)
+tmpf(sleepstudy,REform=~(0+Days|Subject))
+tmpf(sleepstudy,REform=~(1|Subject))
+
+
+
Modified: pkg/lme4Eigen/tests/simulate.R
===================================================================
--- pkg/lme4Eigen/tests/simulate.R 2012-02-08 21:21:15 UTC (rev 1582)
+++ pkg/lme4Eigen/tests/simulate.R 2012-02-08 22:40:20 UTC (rev 1583)
@@ -1,5 +1,5 @@
library(lme4Eigen)
-debug(simulate.merMod)
+## debug(simulate.merMod)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
More information about the Lme4-commits
mailing list