4
votes

I'm trying to compute confidence intervals for groups with a common model using lmList from the lme4 package. It works fine for normal linear model, but fails when the dependent variable is dichotomous. For example, this works fine:

d <- data.frame(
  g = sample(c("A","B","C","D","E"), 250, replace=TRUE),
  y1 = runif(250, max=100),
  y2 = sample(c(0,1), 250, replace=TRUE)
)

library(lme4)

fm1 <- lmList(y1 ~ 1 | g, data=d)

I can extract the coefficients using coef(fm1) and the confidence interval for the coefficients using confint(fm1). Then I run a model with the dichotomous outcome:

fm2 <- lmList(y2 ~ 1 | g, data=d, family=binomial)

I can still get the coefficients using coef(fm2), but when I try to get the confidence intervals, I get an error:

> confint(fm2)
Waiting for profiling to be done...
Waiting for profiling to be done...
Error in val[, , i] <- eval(mCall) : incorrect number of subscripts

I was originally going to post this to stats.stackexchange because I thought it might be something I don't understand about the confidence intervals of GLMs, but then I figured out that I can still get the confidence intervals using

by(d, d$g, function(x) confint(glm(y2 ~ 1, data=x, family=binomial))) 

Is there some way to do this using lmList?

> sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252  
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                  
[5] LC_TIME=German_Germany.1252   

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

other attached packages:
[1] lme4_0.999375-42 Matrix_1.0-6     lattice_0.20-6 

loaded via a namespace (and not attached):
[1] grid_2.15.0   MASS_7.3-17   nlme_3.1-103  stats4_2.15.0 tools_2.15.0 
1
I think this is a bug in the "!missing(family)" branch of lmList, or its processing by confint.lmList. Note that when you try "family=gaussian", it also produces the error in confint. There is nothing wrong in using by() or ddply in this case.Dieter Menne
Ah ha, ok. Nice catch. I'll just go with the by() method then. Thanks for the comment.smillig
@BenBolker Should we file a bug report on this minor issue, or let Douglas focus on more important things such?Dieter Menne

1 Answers

2
votes

It is indeed a bug, I'm working on a fix. It turns out to be fairly straightforward: the basic problem is that confint.lm returns a matrix while confint.glm returns a data frame.

In general I strongly recommend posting all but the most trivial bug fixes to the bug tracker at http://r-forge.r-project.org/tracker/?atid=298&group_id=60&func=browse (I complain about R-core's attitude toward bug reports, so for consistency I have to follow my own rules ...) Thanks for the report!