0

I've been trying to analyze the results of multiple imputations in a multi-level model in R. I know that my lme4 model works. I know that my imputation is working. But using Zelig OR using mitools::MIcombine() both throw similar errors about S4 objects.

library(lme4)
library(mitools)

data(smi)
# this model works fine
model.f.glm <- function(df) {
  glm(drinkreg~wave*sex,family=binomial(), data=df)
}
# but this one will not
model.f.glmer <- function(df) {
  glmer(drinkreg~wave*sex + (1|id),family=binomial(), data=df)
}

# check the base model first:
model.f.glm(smi$imputations[[1]])
model.f.glmer(smi$imputations[[1]])

# use the mitools methods:
models.glm <- lapply(smi$imputations, model.f.glm)
models.glmer <- lapply(smi$imputations, model.f.glmer)

summary(MIcombine(models.glm))
# error:
summary(MIcombine(models.glmer))
# Error in cbar + results[[i]] : non-numeric argument to binary operator

# get MIcombine for them:
betas.glm <- MIextract(models.glm, fun = coef)
vars.glm  <- MIextract(models.glm, fun = vcov)
summary(MIcombine(betas.glm, vars.glm))
betas.glmer <- MIextract(models.glmer, fun = fixef)
vars.glmer  <- MIextract(models.glmer, fun = vcov)
# error:
summary(MIcombine(betas.glmer, vars.glmer))
# Error in diag(vcov(object)) : 
    # no method for coercing this S4 class to a vector

With Zelig and Amelia:

library(Zelig)
library(ZeligMultilevel)
library(Amelia)

data(africa)
a.out <- amelia(x = africa, cs = "country", ts = "year", logs = "gdp_pc")
# glm:
z.out <- zelig(gdp_pc ~ trade + civlib, model = "ls", data = a.out)
summary(z.out)
# glmer:
z.out <- zelig(gdp_pc ~ trade + civlib + (1 | year), model = "ls.mixed", data = a.out)
# ERROR:
summary(z.out)
#Error in diag(vcovlist[[i]]) : 
  # no method for coercing this S4 class to a vector

I've looked over this previous post, where @Roman Luštrik put a bunch of effort into fixing the summary about four years ago -- but it seems to be an issue again. Any thoughts?

JDB
  • 110
  • 6
  • My _guess_ thus far is that something changed in the `lme4` or `Matrix` packages. I find the `Zelig` documentation to be very confusing, but the `mitools::MIcombine()` function seems to have trouble in part because it's called `is.matrix()` on a dsyMatrix class (which returns FALSE despite being a matrix). The work @roman-luštrik did on the previous post mentioned above seems to be outdated now, so I think these packages just aren't up to date? – JDB Jun 30 '17 at 23:22
  • 1
    Looks like `ZeligMultilevel` has been changed to "abandoned" on github; the `merTools` package has some functions that solved this. – JDB Jul 03 '17 at 21:54

1 Answers1

0

To expand on JDB's suggestion, try

vars.glmer  <- MIextract(models.glmer, fun = function(mod) as.matrix(vcov(mod)))
TLM
  • 23
  • 1
  • 7