Lmer dengan data yang dimasukkan multiply

Bagaimana saya bisa mendapatkan efek acak dikumpulkan untuk lmer setelah beberapa kali imputasi?

Saya menggunakan mouse untuk menghubungkan beberapa frame data. Dan lme4 untuk model campuran dengan intersep acak dan kemiringan acak. Pooling lmer baik-baik saja, kecuali bahwa itu tidak mengumpulkan efek acak. Saya telah mencari banyak solusi tanpa keberuntungan. Saya mencoba paket mi, namun saya hanya melihat output gabungan untuk perkiraan dan std.error. Saya sudah mencoba mengekspor objek tikus ke spss tanpa hasil. Saya melihat beberapa diskusi tentang Zelig. Saya pikir itu bisa menyelesaikan masalah saya. Namun saya tidak dapat menemukan cara menggunakan paket dengan data yang diperhitungkan untuk lmer.

Saya tahu paket mouse hanya mendukung penyatuan efek tetap. Apakah ada pekerjaan?

Beberapa imputasi:

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. 

Model bertingkat:

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

Perbarui Hasil dari pooled lmer:

> 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

Informasi ini tidak ada, yang saya dapatkan ketika menjalankan lmer tanpa banyak imputasi:

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

Anda dapat melakukan ini sedikit dengan tangan jika dengan memanfaatkan lapplyfungsi dalam R dan struktur daftar yang dikembalikan oleh Ameliapaket imputasi berganda. Berikut ini contoh skrip cepat.

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

Ameliamirip dengan micesehingga Anda bisa mengganti variabel Anda dari micepanggilan di sini - contoh ini dari proyek yang saya kerjakan.

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

a.outadalah objek imputasi, sekarang kita perlu menjalankan model pada setiap dataset yang diimputasi. Untuk melakukan ini, kami menggunakan lapplyfungsi dalam R untuk mengulangi fungsi lebih dari elemen daftar. Fungsi ini menerapkan fungsi - yang merupakan spesifikasi model - untuk setiap dataset (d) dalam daftar dan mengembalikan hasilnya dalam daftar model.

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

Sekarang kita membuat data.frame dari daftar itu, dengan mensimulasikan nilai efek tetap dan acak menggunakan fungsi FEsim dan REsim dari paket merTools

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

Data.frame di atas termasuk perkiraan terpisah untuk setiap dataset, sekarang kita perlu menggabungkan mereka menggunakan runtuh seperti keruntuhan argumen

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

Sekarang kita juga dapat mengekstraksi beberapa statistik pada varians / kovarians untuk efek acak di seluruh nilai imputasi. Di sini saya telah menulis dua fungsi ekstraktor sederhana untuk melakukan ini.

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

Dan sekarang kita dapat menerapkannya pada model dan menyimpannya sebagai vektor:

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

Memperbarui

Fungsi-fungsi di bawah ini akan membuat Anda lebih dekat dengan output yang disediakan arm::displayoleh operasi pada daftar lmeratau glmerobjek. Semoga ini akan dimasukkan ke dalam merToolspaket dalam waktu dekat:

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

Anda juga dapat menggunakan fungsi testEstimates setelah imputasi menggunakan mouse, testEstimates (as.mitml.result (fm1), var.comp = T) $ var.comp