R（Amelia、zelig、lme4）の多重帰属データセットのマルチレベル回帰モデル

Question

私は多重帰属データ（Ameliaで作成）にマルチレベルモデルを実行しようとしています。

Error in object[[1]]$result$call : 
$ operator not defined for this S4 class 

私はOLS回帰を実行する場合、それは動作します：

サンプルは、このコードは、次のエラーコードを生成グループとクラスタ化されたサンプル= 24に基づいて、N = 150

library("ZeligMultilevel") 
ML.model.0 <- zelig(dv~1 + tag(1|group), model="ls.mixed", 
data=a.out$imputations) 
summary(ML.model.0) 

されます

model.0 <- zelig(dv~1, model="ls", data=a.out$imputations) 
m.0 <- coef(summary(model.0)) 
print(m.0, digits = 2) 

     Value Std. Error t-stat p-value 
[1,] 45  0.34 130 2.6e-285 

の作業例を提供しています。

require(Zelig) 
require(Amelia) 
require(ZeligMultilevel) 

data(freetrade) 
length(freetrade$country) #grouping variable 

#Imputation of missing data 

a.out <- amelia(freetrade, m=5, ts="year", cs="country") 

# Models: (1) OLS; (2) multi-level 

model.0 <- zelig(polity~1, model="ls", data=a.out$imputations) 
m.0 <- coef(summary(model.0)) 
print(m.0, digits = 2) 

ML.model.0 <- zelig(polity~1 + tag(1|country), model="ls.mixed", data=a.out$imputations) 
summary(ML.model.0) 

私の問題は、ZeligがAmeliaのmiクラスとどのようにインターフェイスするかで問題になると思います。したがって、私は別のRパッケージ、lme4に目を向ける。

require(lme4) 
write.amelia(obj=a.out, file.stem="inmi", format="csv", na="NA") 
diff <-list(5) # a list to store each model, 5 is the number of the imputed datasets 

for (i in 1:5) { 
file.name <- paste("inmi", 5 ,".csv",sep="") 
data.to.use <- read.csv(file.name) 
diff[[5]] <- lmer(polity ~ 1 + (1 | country), 
data = data.to.use)} 
diff 

結果は以下の通りです：私はdiff[[4]]でdiff[[5]]を交換する際に

[[1]] 
[1] 5 

[[2]] 
NULL 

[[3]] 
NULL 

[[4]] 
NULL 

[[5]] 
Linear mixed model fit by REML 
Formula: polity ~ 1 + (1 | country) 
    Data: data.to.use 
    AIC BIC logLik deviance REMLdev 
1006 1015 -499.9  1002 999.9 
Random effects: 
Groups Name  Variance Std.Dev. 
country (Intercept) 14.609 3.8222 
Residual    17.839 4.2236 
Number of obs: 171, groups: country, 9 

Fixed effects: 
      Estimate Std. Error t value 
(Intercept) 2.878  1.314 2.19 

結果は同じまま、diff[[3]]などそれでも、私は、これは実際に組み合わせたデータセットの結果であるかどうかを疑問に思ってまたは1つの帰属データセットに対して実行されます。何かご意見は？ありがとう！

Answer 1

このオブジェクトのサマリー機能を変更しました（ソースを取得し、./R/summary.Rファイルを開きました）。コードフローを作成するためにいくつかの中括弧を追加し、getcoefをcoefに変更しました。これはこの特定のケースではうまくいくはずですが、一般的であるかどうかはわかりません。機能getcoefはスロットcoef3を検索し、これは見たことがありません。多分@BenBolkerはここに目を投げることができますか？これは結果がどのようなものかは保証できませんが、出力は当然です。おそらく、将来のバージョンでこれを修正するために、パッケージ作成者に連絡することができます。

要約（ML.model.0）

Model: ls.mixed 
    Number of multiply imputed data sets: 5 

Combined results: 

Call: 
zelig(formula = polity ~ 1 + tag(1 | country), model = "ls.mixed", 
    data = a.out$imputations) 

Coefficients: 
     Value Std. Error t-stat p-value 
[1,] 2.902863 1.311427 2.213515 0.02686218 

For combined results from datasets i to j, use summary(x, subset = i:j). 
For separate results, use print(summary(x), subset = i:j). 

が変更された機能：

summary.MI <- function (object, subset = NULL, ...) { 
    if (length(object) == 0) { 
    stop('Invalid input for "subset"') 
    } else { 
    if (length(object) == 1) { 
     return(summary(object[[1]])) 
    } 
    } 

    # Roman: This function isn't fecthing coefficients robustly. Something goes wrong. Contact package author. 
    getcoef <- function(obj) { 
    # S4 
    if (!isS4(obj)) { 
     coef(obj) 
    } else { 
     if ("coef3" %in% slotNames(obj)) { 
     obj@coef3 
     } else { 
     obj@coef 
     } 
    } 
    } 

    # 
    res <- list() 

    # Get indices 
    subset <- if (is.null(subset)) { 
     1:length(object) 
    } else { 
     c(subset) 
    } 

    # Compute the summary of all objects 
    for (k in subset) { 
     res[[k]] <- summary(object[[k]]) 
    } 


    # Answer 
    ans <- list(
     zelig = object[[1]]$name, 
     call = object[[1]]$result@call, 
     all = res 
    ) 

    # 
    coef1 <- se1 <- NULL 

    # 
    for (k in subset) { 
#  tmp <- getcoef(res[[k]]) # Roman: I changed this to coef, not 100% sure if the output is the same 
     tmp <- coef(res[[k]]) 
     coef1 <- cbind(coef1, tmp[, 1]) 
     se1 <- cbind(se1, tmp[, 2]) 
    } 

    rows <- nrow(coef1) 
    Q <- apply(coef1, 1, mean) 
    U <- apply(se1^2, 1, mean) 
    B <- apply((coef1-Q)^2, 1, sum)/(length(subset)-1) 
    var <- U+(1+1/length(subset))*B 
    nu <- (length(subset)-1)*(1+U/((1+1/length(subset))*B))^2 

    coef.table <- matrix(NA, nrow = rows, ncol = 4) 
    dimnames(coef.table) <- list(rownames(coef1), 
           c("Value", "Std. Error", "t-stat", "p-value")) 
    coef.table[,1] <- Q 
    coef.table[,2] <- sqrt(var) 
    coef.table[,3] <- Q/sqrt(var) 
    coef.table[,4] <- pt(abs(Q/sqrt(var)), df=nu, lower.tail=F)*2 
    ans$coefficients <- coef.table 
    ans$cov.scaled <- ans$cov.unscaled <- NULL 

    for (i in 1:length(ans)) { 
     if (is.numeric(ans[[i]]) && !names(ans)[i] %in% c("coefficients")) { 
     tmp <- NULL 
     for (j in subset) { 
      r <- res[[j]] 
      tmp <- cbind(tmp, r[[pmatch(names(ans)[i], names(res[[j]]))]]) 
     } 
     ans[[i]] <- apply(tmp, 1, mean) 
     } 
    } 

    class(ans) <- "summaryMI" 
    ans 
    }