views:

956

answers:

3

I have fit a regression using lme4 thanks to a previous answer. Now that I have a regression fit for each state I'd like to use lattice to plot QQ plots for each state. I would also like to plot error plots for each state in a lattice format. How do I make a lattice plot using the results of a lme4 regression?

Below is a simple sample (yeah, I like a good alliteration) using two states. I would like to make a two panel lattice made from the object fits.

library(lme4)
d <- data.frame(state=rep(c('NY', 'CA'), c(10, 10)), year=rep(1:10, 2), response=c(rnorm(10), rnorm(10)))
fits <- lmList(response ~ year | state, data=d)
+2  A: 

I am not sure you can get this into lattice easily. What you have in fits is a an S4 object containing a .Data slot with a list of standard lm objects:

R> class(fits)
[1] "lmList"
attr(,"package")
[1] "lme4"
R> class([email protected])
[1] "list"
R> class([email protected][[1]])
[1] "lm"
R> op <- par(mfrow=c(2,4))
R> invisible(lapply([email protected], plot))

This last plot simply plots you the standard 2x2 plot for lm objects twice, once for each element of the list of fitted objects. Use the which argument to plot to select subsets of these or for other regression diagnostics.

If you actually want lattice plots of predicted vs actual, you may have to program this.

Dirk Eddelbuettel
+4  A: 

Instead of using lmList, I'd recommend the more general plyr package.

library(plyr)

d <- data.frame(
 state = rep(c('NY', 'CA'), c(10, 10)), 
 year = rep(1:10, 2), 
 response = c(rnorm(10), rnorm(10))
)

# Create a list of models
# dlply = data frame -> list
models <- dlply(d, ~ state, function(df) { 
  lm(response ~ year, data = df)
})

# Extract the coefficients in a useful form
# ldply = list -> data frame
ldply(models, coef)

# We can get the predictions in a similar way, but we need
# to cast to a data frame so the numbers come out as rows,
# not columns.
predictions <- ldply(models, as.data.frame(predict))

predictions is a regular R data frame and so is easy to plot.

hadley
well you just sold me on learning the plyr package. :) I had gone with lme4 because I could grok what it was doing easier than what plyr does. Now I see an advantage to plyr. Thanks.
JD Long
I'm having a heck of a time joining my predicted values back to the original data set. In the example above NY is first in the input data.frame and CA is first in the output.The field [year] is not retained in the data.frame [predictions]. This is probably obvious, but how do I join these back?
JD Long
instead of `as.data.frame(predict)`, use something like function(model) { transform(d, pred = predict(model)) }
hadley
I've gotten so much good use out of this method this week. Thanks for your help!
JD Long
A: 

I've had some trouble with lme4::lmList. For example, summary doesn't seem to work. So you might run into some problems because of that.

So even though I use lmer, instead of lme, I've been explicitly calling nlme::lmList instead. Then summary etc. will work.

Dean Eckles