lmer con datos imputados multiplicados

10

¿Cómo puedo obtener efectos aleatorios agrupados para lmer después de una imputación múltiple?

Estoy usando ratones para imputar múltiples un marco de datos. Y lme4 para un modelo mixto con intercepción aleatoria y pendiente aleatoria. Agrupar lmer funciona bien, excepto que no agrupa los efectos aleatorios. He buscado mucho una solución sin suerte. Probé el paquete mi, sin embargo, solo veo resultados agrupados para la estimación y std.error. Intenté exportar objetos de ratones a spss sin suerte. Vi un poco de discusión sobre Zelig. Pensé que eso podría resolver mi problema. Sin embargo, no pude descubrir cómo usar el paquete con datos imputados para lmer.

Sé que el paquete de ratones solo admite la combinación de los efectos fijos. ¿Hay alguna solución?

Imputación múltiple:

library(mice)
Data <- subset(Data0, select=c(id, faculty, gender, age, age_sqr, occupation, degree, private_sector, overtime, wage))
ini <- mice(Data, maxit=0, pri=F) #get predictor matrix
pred <- ini$pred
    pred[,"id"] <- 0 #don't use id as predictor
    meth <- ini$meth
meth[c("id", "faculty", "gender", "age", "age_sqr", "occupation", "degree", "private_sector", "overtime", "wage")] <- "" #don't impute these variables, use only as predictors.
imp <- mice(Data, m=22, maxit=10, printFlag=TRUE, pred=pred, meth=meth) #impute Data with 22 imputations and 10 iterations. 

Modelo multinivel:

library(lme4)
    fm1 <- with(imp, lmer(log(wage) ~ gender + age + age_sqr + occupation + degree + private_sector + overtime + (1+gender|faculty))) #my multilevel model
    summary(est <- pool(fm1)) #pool my results

Actualizar resultados de lmer agrupados:

> summary(est <- pool(fm1))
                                est           se            t       df     Pr(>|t|)         lo 95         hi 95 nmis       fmi    lambda
(Intercept)   7,635148e+00 0,1749178710 43,649905006 212,5553 0,000000e+00  7,2903525425  7,9799443672   NA 0,2632782 0,2563786
Gender        -1,094186e-01 0,0286629154 -3,817427078 117,1059 2,171066e-04 -0,1661834550 -0,0526537238   NA 0,3846276 0,3742069
Occupation1   1,125022e-01 0,0250082538  4,498601518 157,6557 1,320753e-05  0,0631077322  0,1618966049   NA 0,3207350 0,3121722
Occupation2   2,753089e-02 0,0176032487  1,563966385 215,6197 1,192919e-01 -0,0071655902  0,0622273689   NA 0,2606725 0,2538465
Occupation3   1,881908e-04 0,0221992053  0,008477365 235,3705 9,932433e-01 -0,0435463305  0,0439227120   NA 0,2449795 0,2385910
Age           1,131147e-02 0,0087366178  1,294719230 187,0021 1,970135e-01 -0,0059235288  0,0285464629    0 0,2871640 0,2795807
Age_sqr       -7,790476e-05 0,0001033263 -0,753968159 185,4630 4,518245e-01 -0,0002817508  0,0001259413    0 0,2887420 0,2811131
Overtime      -2,376501e-03 0,0004065466 -5,845581504 243,3563 1,614693e-08 -0,0031773002 -0,0015757019    9 0,2391179 0,2328903
Private_sector  8,322438e-02 0,0203047665  4,098760934 371,9971 5,102752e-05  0,0432978716  0,1231508962   NA 0,1688478 0,1643912

Falta esta información, que obtengo cuando ejecuto lmer sin imputación múltiple:

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Faculty  (Intercept) 0,008383 0,09156      
          Genderfemale0,002240 0,04732  1,00
 Residual             0,041845 0,20456      
Number of obs: 698, groups:  Faculty, 17
Helgi Guðmundsson
fuente
¿Su problema es que no sabe cómo caracterizar la incertidumbre RE después de un IM? No entiendo qué procedimientos intenta hacer su código.
generic_user
(1 + género | facultad): género como pendiente aleatoria, facultad como intercepción aleatoria. Estoy tratando de obtener resultados agrupados de las 22 imputaciones para los efectos aleatorios (género y facultad)
Helgi Guðmundsson
Pequeña actualización Cuando imputo múltiples en SPSS y ejecuto un modelo mixto; SPSS solo agrupa los efectos fijos, no los efectos aleatorios. Lo mismo ocurre con el paquete mi para R. Estoy empezando a pensar que esto no es posible.
Helgi Guðmundsson
2
En respuesta a Helgi: es estadísticamente posible hacerlo: Stata proporciona estimaciones agrupadas de las estimaciones de los componentes de la varianza después de usar la imputación múltiple. La única dificultad es obtener las estimaciones y los errores estándar de los componentes de la varianza, y que la agrupación debe realizarse en una escala para la cual el posterior es aproximadamente normal. Creo que Stata hace la agrupación en la escala de desviación estándar de registro para hacer la aproximación más razonable.
Jonathan Bartlett

Respuestas:

9

Puede hacerlo un poco a mano si aprovecha la lapplyfuncionalidad de R y la estructura de lista que devuelve el Ameliapaquete de imputación múltiple. Aquí hay un script de ejemplo rápido.

library(Amelia)
library(lme4)
library(merTools)
library(plyr) # for collapsing estimates

Ameliaes similar a micelo que puede sustituir sus variables desde la micellamada aquí; este ejemplo es de un proyecto en el que estaba trabajando.

 a.out <- amelia(dat[sub1, varIndex], idvars = "SCH_ID", 
            noms = varIndex[!varIndex %in% c("SCH_ID", "math12")], 
            m = 10)

a.outes el objeto de imputación, ahora necesitamos ejecutar el modelo en cada conjunto de datos imputado. Para hacer esto, usamos la lapplyfunción en R para repetir una función sobre los elementos de la lista. Esta función aplica la función, que es la especificación del modelo, a cada conjunto de datos (d) en la lista y devuelve los resultados en una lista de modelos.

 mods <- lapply(a.out$imputations,
           function(d) lmer((log(wage) ~ gender + age + age_sqr + 
            occupation + degree + private_sector + overtime + 
             (1+gender|faculty), data = d)

Ahora creamos un data.frame a partir de esa lista, simulando los valores de los efectos fijos y aleatorios usando las funciones FEsim y REsim del paquete merTools

imputeFEs <- ldply(mods, FEsim, nsims = 1000)
imputeREs <- ldply(mods, REsim, nsims = 1000)

Los data.frames anteriores incluyen estimaciones separadas para cada conjunto de datos, ahora necesitamos combinarlos usando un colapso como colapso de argumento

imputeREs <- ddply(imputeREs, .(X1, X2), summarize, mean = mean(mean), 
               median = mean(median), sd = mean(sd), 
               level = level[1])

imputeFEs <- ddply(imputeFEs, .(var), summarize, meanEff = mean(meanEff), 
               medEff = mean(medEff), sdEff = mean(sdEff))

Ahora también podemos extraer algunas estadísticas sobre la varianza / covarianza de los efectos aleatorios en los valores imputados. Aquí he escrito dos funciones extractoras simples para hacer esto.

REsdExtract <- function(model){
  out <- unlist(lapply(VarCorr(model), attr, "stddev"))
  return(out)
}

REcorrExtract <- function(model){
  out <- unlist(lapply(VarCorr(model), attr, "corre"))
  return(min(unique(out)))
}

Y ahora podemos aplicarlos a los modelos y almacenarlos como un vector:

modStats <- cbind(ldply(mods, REsdExtract), ldply(mods, REcorrExtract))

Actualizar

Las funciones a continuación te llevará mucho más cerca de la salida proporcionada por arm::displayoperando en la lista de lmero glmerobjetos. Esperemos que esto se incorpore al merToolspaquete en un futuro cercano:

# Functions to extract standard deviation of random effects from model
REsdExtract <- function(model){
  out <- unlist(lapply(VarCorr(model), attr, "stddev"))
  return(out)
}

#slope intercept correlation from model
REcorrExtract <- function(model){
  out <- unlist(lapply(VarCorr(model), attr, "corre"))
  return(min(unique(out)))
}

modelRandEffStats <- function(modList){
  SDs <- ldply(modList, REsdExtract)
  corrs <- ldply(modList, REcorrExtract)
  tmp <- cbind(SDs, corrs)
  names(tmp) <- c("Imp", "Int", "Slope", "id", "Corr")
  out <- data.frame(IntSD_mean = mean(tmp$Int), 
                        SlopeSD_mean = mean(tmp$Slope), 
                    Corr_mean = mean(tmp$Corr), 
                        IntSD_sd = sd(tmp$Int),
                    SlopeSD_sd = sd(tmp$Slope), 
                        Corr_sd = sd(tmp$Corr))
  return(out)
}

modelFixedEff <- function(modList){
  require(broom)
  fixEst <- ldply(modList, tidy, effects = "fixed")
  # Collapse
  out <- ddply(fixEst, .(term), summarize,
               estimate = mean(estimate), 
               std.error = mean(std.error))
  out$statistic <- out$estimate / out$std.error
  return(out)
}

print.merModList <- function(modList, digits = 3){
  len <- length(modList)
  form <- modList[[1]]@call
  print(form)
  cat("\nFixed Effects:\n")
  fedat <- modelFixedEff(modList)
  dimnames(fedat)[[1]] <- fedat$term
  pfround(fedat[-1, -1], digits)
  cat("\nError Terms Random Effect Std. Devs\n")
  cat("and covariances:\n")
  cat("\n")
  ngrps <- length(VarCorr(modmathG[[1]]))
  errorList <- vector(mode = 'list', length = ngrps)
  corrList <- vector(mode = 'list', length = ngrps)
  for(i in 1:ngrps){
    subList <- lapply(modList, function(x) VarCorr(x)[[i]])
    subList <- apply(simplify2array(subList), 1:2, mean)
    errorList[[i]] <- subList
    subList <- lapply(modList, function(x) attr(VarCorr(x)[[i]], "corre"))
    subList <- min(unique(apply(simplify2array(subList), 1:2, function(x) mean(x))))
    corrList[[i]] <- subList
  }
  errorList <- lapply(errorList, function(x) {
    diag(x) <- sqrt(diag(x))
    return(x)
    })

  lapply(errorList, pfround, digits)
  cat("\nError Term Correlations:\n")
  lapply(corrList, pfround, digits)
  residError <- mean(unlist(lapply(modList, function(x) attr(VarCorr(x), "sc"))))
  cat("\nResidual Error =", fround(residError,
                                             digits), "\n")
  cat("\n---Groups\n")
  ngrps <- lapply(modList[[1]]@flist, function(x) length(levels(x)))
  modn <- getME(modList[[1]], "devcomp")$dims["n"]
  cat(sprintf("number of obs: %d, groups: ", modn))
  cat(paste(paste(names(ngrps), ngrps, sep = ", "),
            collapse = "; "))
  cat("\n")
  cat("\nModel Fit Stats")
  mAIC <- mean(unlist(lapply(modList, AIC)))
  cat(sprintf("\nAIC = %g", round(mAIC, 1)))
  moDsigma.hat <- mean(unlist(lapply(modmathG, sigma)))
  cat("\nOverdispersion parameter =", fround(moDsigma.hat,
                                             digits), "\n")
}
jknowles
fuente
1
Esta funcionalidad, la mayor parte, está integrada en la versión de desarrollo del paquete merTools. Esta versión se lanzará a CRAN en la próxima semana.
jknowles
¿podría decir qué función buscar para hacer esto con el paquete merTools? No pude encontrar nada.
smillig
1
No está completamente documentado en la versión actual, pero echa un vistazo lmerModListal printmétodo que combina los resultados de las listas de modelos.
jknowles
0

También puede usar la función testEstimates después de la imputación con ratones, testEstimates (as.mitml.result (fm1), var.comp = T) $ var.comp

Nidhi Menon
fuente