lmer avec multiplier les données imputées

10

Comment puis-je obtenir des effets aléatoires groupés pour lmer après une imputation multiple?

J'utilise des souris pour imputer plusieurs trames de données. Et lme4 pour un modèle mixte avec interception aléatoire et pente aléatoire. La mise en commun de lmer se passe bien, sauf qu'elle ne regroupe pas les effets aléatoires. J'ai beaucoup cherché une solution sans chance. J'ai essayé le paquet mi, mais je ne vois que la sortie groupée pour l'estimation et std.error. J'ai essayé d'exporter des objets souris vers spss sans aucune chance. J'ai vu une discussion sur Zelig. J'ai pensé que cela pourrait résoudre mon problème. Je n'ai cependant pas été en mesure de comprendre comment utiliser le package avec des données imputées pour lmer.

Je sais que le package de souris ne prend en charge que la mise en commun des effets fixes. Y at-il un travail autour?

Imputation multiple:

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. 

Modèle à plusieurs niveaux:

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

Mettre à jour les résultats de lmer groupé:

> 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

Ces informations sont manquantes, ce que j'obtiens lorsque je lance lmer sans imputation multiple:

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
la source
Votre problème est-il que vous ne savez pas comment caractériser l'incertitude des ER après l'IM? Je ne sais pas quelles procédures votre code essaie de faire.
generic_user
(1 + genre | faculté): sexe comme pente aléatoire, faculté comme interception aléatoire. J'essaie d'obtenir des résultats regroupés à partir des 22 imputations pour les effets aléatoires (sexe et faculté)
Helgi Guðmundsson
Petite mise à jour. Lorsque je impute plusieurs dans SPSS et exécute un modèle mixte; SPSS regroupe uniquement les effets fixes, pas les effets aléatoires. Il en va de même pour le package mi pour R. Je commence à penser que ce n'est pas possible de le faire.
Helgi Guðmundsson
2
En réponse à Helgi: Il est statistiquement possible de le faire - Stata fournit des estimations groupées des estimations des composantes de la variance après avoir utilisé l'imputation multiple. La seule difficulté est d'obtenir les estimations et les erreurs-types des composantes de la variance, et que la mise en commun se fasse sur une échelle pour laquelle le postérieur est à peu près normal. Je crois que Stata fait la mise en commun sur l'échelle de l'écart type de log pour rendre l'approximation plus raisonnable.
Jonathan Bartlett

Réponses:

9

Vous pouvez le faire quelque peu à la main si vous profitez des lapplyfonctionnalités de R et de la structure de liste renvoyée par le Ameliapackage d'imputation multiple. Voici un exemple de script rapide.

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

Ameliaest similaire à micedonc vous pouvez simplement remplacer vos variables depuis l' miceappel ici - cet exemple provient d'un projet sur lequel je travaillais.

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

a.outest l'objet d'imputation, nous devons maintenant exécuter le modèle sur chaque ensemble de données imputé. Pour ce faire, nous utilisons la lapplyfonction dans R pour répéter une fonction sur les éléments de la liste. Cette fonction applique la fonction - qui est la spécification du modèle - à chaque ensemble de données (d) de la liste et renvoie les résultats dans une liste de modèles.

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

Nous créons maintenant un data.frame à partir de cette liste, en simulant les valeurs des effets fixes et aléatoires en utilisant les fonctions FEsim et REsim du package merTools

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

Les data.frames ci-dessus incluent des estimations distinctes pour chaque ensemble de données, maintenant nous devons les combiner en utilisant un effondrement comme un effondrement d'argument

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))

Maintenant, nous pouvons également extraire quelques statistiques sur la variance / covariance pour les effets aléatoires sur les valeurs imputées. Ici, j'ai écrit deux fonctions d'extraction simples pour ce faire.

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)))
}

Et maintenant, nous pouvons les appliquer aux modèles et les stocker en tant que vecteur:

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

Mise à jour

Les fonctions ci-dessous vous rapprocheront beaucoup plus de la sortie fournie par arm::displayen opérant sur la liste des objets lmerou glmer. Espérons que cela sera incorporé dans le merToolspaquet dans un proche avenir:

# 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
la source
1
Cette fonctionnalité - la plupart - est intégrée dans la version de développement du package merTools. Cette version sera transmise au CRAN dans la semaine à venir.
jknowles
pourriez-vous dire quelle fonction rechercher pour le faire avec le package merTools? Je n'ai rien trouvé.
smillig
1
Ce n'est pas entièrement documenté dans la version actuelle, mais consultez lmerModListet la printméthode, qui combine les résultats des listes de modèles.
jknowles
0

Vous pouvez également utiliser la fonction testEstimates après imputation à l'aide de souris, testEstimates (as.mitml.result (fm1), var.comp = T) $ var.comp

Nidhi Menon
la source