Question

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 
Was it helpful?

Solution

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!

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top