need help with predicted values from lmer with newdata

Hi Paul,

I think this may be a bug in predict.merMod. ?It is not terribly
elegant, but I suggest changing line 102 from:

re_new <- unlist(re_List[m])

to

re_new <- as.vector(do.call("rbind", re_List[m]))

mypredict <- lme4:::predict.merMod

I believe the issue with the current is that the random effects are
unlisted which essentially 'stacks' the vector, so you have the random
effects one at a time by effect, not by level (i.e., all intercepts,
then all slopes). ?Later this is post multiplied by a matrix that is
sorted by level, the dimensions are correct so the matrices multiply
fine, but the results are not what I believe was intended.

I cannot see a way to add a patch to R-forge or submit a bug report so
I cced Doug on this.

Cheers,

Josh

On Thu, Jun 21, 2012 at 8:28 PM, Paul Johnson <pauljohn32 at gmail.com> wrote:

Greetings.

I'm getting some crazy results out of predict on an lmer model. ?I've
gone through
this lots of times, and feel certain I'm missing some obvious mistake.
?If you browse
through and get to the bottom, there's some output that shows my predict output
is just wrong, but I can't see why.

I will send the first person with the right answer a new KU t-shirt in
the size of
your choice. ?That's better than "thanks in advance". Isn't it?

## Paul Johnson <pauljohn at ku.edu>
## 2012-06-21

## I want to compare the "separate ols regressions on clusters"
## results with the mixed model estimates. I've run into
## trouble because the output of predict with a newdata object
## seems to be completely wrong. ?I demonstrate "manual" calculations
## to try to convince you I understand what's going on (but
## peculiarity of predict output may make things appear otherwise).

## I pasted in some output at the bottom in case you just want to
## see where I think there is a problem.

## My system:
## > sessionInfo()
## R version 2.15.0 (2012-03-30)
## Platform: x86_64-pc-linux-gnu (64-bit)

## locale:
## ?[1] LC_CTYPE=en_US.UTF-8 ? ? ? LC_NUMERIC=C
## ?[3] LC_TIME=en_US.UTF-8 ? ? ? ?LC_COLLATE=en_US.UTF-8
## ?[5] LC_MONETARY=en_US.UTF-8 ? ?LC_MESSAGES=en_US.UTF-8
## ?[7] LC_PAPER=C ? ? ? ? ? ? ? ? LC_NAME=C
## ?[9] LC_ADDRESS=C ? ? ? ? ? ? ? LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

## attached base packages:
## [1] stats ? ? graphics ?grDevices utils ? ? datasets ?methods ? base

## other attached packages:
## [1] lme4_0.999902344-0 Matrix_1.0-6 ? ? ? lattice_0.20-6 ? ? MASS_7.3-18
## loaded via a namespace (and not attached):
## [1] compiler_2.15.0 grid_2.15.0 ? ? minqa_1.2.1 ? ? nlme_3.1-104
## [5] splines_2.15.0 ?tools_2.15.0


## Topic: create multi-level data and use regression and lmer to
## estimate it. ?This creates 3 different dependent variables,
## y1: ordinary "homoskedastic" regression (no grouping effect)
## y2: clustered "random intercepts"
## y3: clustered "random intercepts" and "random slope" for 1 predictor (x1)

library(MASS)
set.seed(1234)

## Step 1. Create a Data Frame.

## M respondents, N observations for each one.
## In a repeated measures context, this is called "longitudinal data".
## In cross sectional approach, this is M groups (classrooms) with
## several people in each)

M <- 10
N <- 30

## get M unique gray colors, will use later for plotting
grays <- gray.colors(M)

## Standard deviation of error term at individual level
STDE <- 48

## STDEM: standard deviation of clustered intercepts.
### In a longitudinal "repeated measures" exercise, this is an
## individual-level effect. In a cross section, this is a random
## classroom effect.
STDEM <- 30

## STEx1: standard deviation of slopes across cluster units
STDEx1 <- 5

## The true effects of b0, b1, b2, b3 in
## y = b0 + b1*x1 + b2*x2 + b3*x3
bslopes <- c(0.2, 15.4, -0.2, 3)

## Now generate the data frame with x1, x2 and x3.
## Let's suppose the predictors are Multivariate Normal, with
## means (100,200,150), standard deviations (10, 20, 30), and
## intercorrelations of 0.4. ?These can, of course, be adjusted
## for interesting variations.

## Mind, an "indicator" that says with which cluster an observation belongs.
Mind <- 1:M %x% rep(1,N)

means <- c(100, 200, 150)
sds <- c(10, 20, 30)
rho <- 0.4
corr.mat <- matrix(c(1, rho, rho, rho, 1, rho, rho, rho, 1), nrow = 3)
sigma <- diag(sds) %*% corr.mat %*% diag(sds)
x.mat <- mvrnorm(n = N * M, mu = means, Sigma = sigma)
dimnames(x.mat)[[2]] <- c("x1", "x2", "x3")
## Add an intercept on the front of x.mat
x.mat <- cbind("(Intercept)"= 1, x.mat)


## Create a dependent variable that has no clustering effects. ?This
## is an "ordinary" regression with slopes and random error designated
## above random noise
y1 <- x.mat %*% bslopes + rnorm(M*N, m=0, s= STDE)
dat <- data.frame(id=1:(N*M), Mind, y1, as.data.frame(x.mat))
rm(Mind, y1, x.mat) ## cleanup workspace

## Layer on additive group level error in y2
## Add a group-level intercept.
## A vector of M disturbances is created, and each is
## added to all N case within its cluster
dat$y2 <- dat$y1 + ?rnorm(M, 0, STDEM) %x% rep(1, N)

## In y3, add in random slope effect
## are changing the slope for x1 from b1 to N(b1, STDEx1^2).
## Same as newb1 ~ b1 + N(0, STDEx1) * x1
dat$y3 <- dat$y2 + ? (rnorm(M, 0, STDEx1) %x% rep(1, N)) * dat$x1

## In the lme4 package, there is an "easy" tool to run M separate lm regressions
library(lme4)

##100 separate regressions, with 3 predictors
m3list <- lmList(y3 ~ x1 + x2 + x3 | Mind, data=dat, pool = FALSE)


## Now I see a pattern. I'll use the unique range on x1 for each cluster
## I think that makes the lines look more "in" the data.
plot(y3 ~ x1, data=dat, col=grays[Mind], main = "lm on clusters")
for( i in seq_along(m3list)){
? ?m3mf <- model.frame(m3list[[i]]) #data set for group i
? ?x1range <- range(m3mf$x1) ## use group-specific ranges this time
? ?pgroup <- predict( m3list[[i]], newdata = data.frame(x1=x1range,
x2=mean(m3mf$x2), x3=mean(m3mf$x3)))
? ?lines(x1range, pgroup, col=grays[i])
}

## Hmm. That data looks like it might be handy for comparison later on.
## Build a data frame of it. I've never thought of doing this work
## cluster-by-cluster before.

m3newdat <- lapply(m3list, function(x) {
? ?m3mf <- model.frame(x)
? ?ndf = data.frame(x1=range(m3mf$x1), x2=mean(m3mf$x2), x3=mean(m3mf$x3))
? ?ndf$m3pred <- predict(x, newdata = ndf)
? ?ndf} )
## Smash the list of data frames together
m3newdat <- do.call("rbind", m3newdat)
## Interesting!
m3newdat[1:20, ]
## Better add a variable Mind. This looks stupid, but works.
m3newdat$Mind <- as.integer( do.call("rbind",
strsplit(row.names(m3newdat), split="\\.", perl=T))[ ,1] )


## Draw new graphs on a new device, so we can compare
dev.new()

## The "right" model allows the random effects to be correlated.
mm4 <- lmer( y3 ~ x1 + x2 + x3 + (x1 | Mind), data=dat)
summary(mm4)



mm4ranef <- ranef(mm4, postVar = TRUE) ## a ranef.mer object, a list
that includes one data frame


### What do the predicted values mean from lmer?
### here they are. But what do they include?
head(mm4pred <- predict(mm4))
tail(mm4pred <- predict(mm4))


## I can't tell what those numbers are. Can't be sure from ?predict.merMod
## So I'm going to go through the step-by-step
## process of calculating predicted values. These should follow a
formula like this
## for model 4. j is the "cluster" index.
##
## ?j'th intercept ? ? ? ? ? ? ? ?j'th slope
## ?(b0 + raneff(intercept, j)) + (b1 + raneff(x1,j))*x1 + b2 x2 + b3 x3

mm4b <- coef(summary(mm4))[ ,1]
mm4vnames <- names(mm4b)
mm4mm <- model.matrix(mm4) ## predictors

## Doing this bit by bit. First
## b0 + raneff(intercept, j)
mm4inteffect <- mm4b["(Intercept)"] + ?mm4ranef[[1]][dat$Mind, 1]
mm4x1effect <- mm4mm[ , c("x1")] * (mm4b["x1"] + mm4ranef[[1]][dat$Mind, 2])
mm4pred2 <- mm4inteffect + mm4x1effect + ?mm4mm[ ,c("x2","x3") ] %*%
mm4b[c("x2","x3")]
head(mm4pred2)
tail(mm4pred2)

## Aha! Those exactly match predict for mm4. So I understand what "predict" does




## Now, I want to run predict for some particular values of x1, x2, x3.

## Let's try running m3newdat through the predict function for mm4.
m3newdat$mm4.3 <- predict(mm4, newdata = m3newdat)
m3newdat


## Disaster. Predicted values completely out of whack. Lots of
## lines fall off the bottom of the graph.

plot(y3 ~ x1, data=dat, col=grays[Mind], main="lmer mixed model predictions")
by(m3newdat, m3newdat$Mind, function(x) { lines(x$x1, x$mm4.3,
col=grays[x$Mind])})


## They don't match the manual calculations I perform here
mm4b <- coef(summary(mm4))[ ,1]
mm4vnames <- names(mm4b)
mm4mm <- as.matrix(m3newdat)

## ?b0 + raneff(intercept, j)
mm4inteffect <- mm4b["(Intercept)"] + rep(mm4ranef[[1]][, 1], each=2)
mm4x1effect <- mm4mm[ , c("x1")] * rep(mm4b["x1"] + mm4ranef[[1]][, 2], each=2)
mm4manualpred2 <- mm4inteffect + mm4x1effect + ?mm4mm[ ,c("x2","x3") ]
%*% ?mm4b[c("x2","x3")]

mm4mm <- cbind(as.data.frame(mm4mm), "mm4manualpred" = mm4manualpred2)
mm4mm

## Here's what I see when I do that. m3pred: the lmList predictions.
## mm4manualpred: my "manual" calculations of predicted values for lmer mm4
## mm4.3 output from predict for same "newdat" used with mm4manualpred.

## > mm4mm
## ? ? ? ? ? ? x1 ? ? ? x2 ? ? ? x3 ? m3pred Mind ? ? ?mm4.3 mm4manualpred
## 1.1 ? 88.51074 196.0615 141.3859 1221.477 ? ?1 ?-434.7602 ? ? ?1239.152
## 1.2 ?110.88258 196.0615 141.3859 1483.957 ? ?1 ?-630.0463 ? ? ?1466.961
## 2.1 ? 80.26167 194.1252 133.7265 1405.672 ? ?2 ?-614.3256 ? ? ?1388.751
## 2.2 ?117.08533 194.1252 133.7265 1845.499 ? ?2 -1051.4105 ? ? ?1866.505
## 3.1 ? 80.84091 202.4517 154.3976 1302.205 ? ?3 ?3527.8154 ? ? ?1284.706
## 3.2 ?120.33775 202.4517 154.3976 1692.819 ? ?3 ?5037.6392 ? ? ?1717.210
## 4.1 ? 68.96320 204.5145 151.6662 1321.927 ? ?4 ?5097.5261 ? ? ?1379.251
## 4.2 ?120.17029 204.5145 151.6662 2137.589 ? ?4 ?8559.2362 ? ? ?2108.628
## 5.1 ? 84.85202 200.7799 153.1595 1821.135 ? ?5 ?4736.2240 ? ? ?1813.415
## 5.2 ?122.58730 200.7799 153.1595 2418.381 ? ?5 ?6641.2026 ? ? ?2426.855
## 6.1 ? 84.60214 202.3324 153.2583 2015.183 ? ?6 ?1521.6831 ? ? ?2001.700
## 6.2 ?128.13418 202.3324 153.2583 2778.684 ? ?6 ?2086.4726 ? ? ?2799.053
## 7.1 ? 74.73548 198.5761 155.8984 1780.196 ? ?7 ?1497.6761 ? ? ?1802.381
## 7.2 ?118.97636 198.5761 155.8984 2622.083 ? ?7 ?2127.8282 ? ? ?2603.626
## 8.1 ? 83.40882 200.4991 152.0515 2342.043 ? ?8 ?1954.4592 ? ? ?2323.477
## 8.2 ?118.73997 200.4991 152.0515 3088.686 ? ?8 ?2601.6007 ? ? ?3105.182
## 9.1 ? 82.48090 205.4736 152.2774 1979.852 ? ?9 ?2253.5650 ? ? ?1950.388
## 9.2 ?117.48279 205.4736 152.2774 2554.581 ? ?9 ?3027.9842 ? ? ?2587.319
## 10.1 ?84.77764 205.9689 152.3325 2256.434 ? 10 ?2236.3346 ? ? ?2268.620
## 10.2 119.66417 205.9689 152.3325 3022.053 ? 10 ?2980.1202 ? ? ?3012.406


--
Paul E. Johnson
Professor, Political Science ? ?Assoc. Director
1541 Lilac Lane, Room 504 ? ? Center for Research Methods
University of Kansas ? ? ? ? ? ? ? University of Kansas
http://pj.freefaculty.org ? ? ? ? ? ?http://quant.ku.edu

--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/