повторить Stata mi Impute регрессии в R или проанализировать вмененный набор данных из Stata

Question

Я пытаюсь повторить анализ Stata (в R), который использует множественное вменение. Шаг mi принимает форму:

mi impute regress y x1, add(5) rseed(123)
mi estimate:reg y trt

Я знаю из Stata docs , что эта команда вменяется с использованием линейной регрессии. Параметр add(5) указывает 5 импутаций.

Было бы полезно, если бы кто-то

1. знал, как повторить это в R, или
2. знает, как я может анализировать вмененный набор данных Stata в R (например, lm(y ~ trt, data=df_mi) ... с вмененными данными)

Для репрезента я создал набор игрушечных данных в R, выполнил команду impute в Stata 12 и экспортировал вмененный набор данных. Начиная с конца следующего фрагмента кода, я покажу, как получить вмененные данные игрушки.

В исходном наборе данных игрушки было 100 единиц измерения, а 20 отсутствовало на y. Вмененный набор данных из Stata имеет 200 obs: исходные 100 плюс 5 наборов вмененных значений для 20 пропусков.

# simulate original data with missing
# set.seed(1)
# df <- data.frame(id = seq(1:100),
#                  trt = c(rep(0, 50), rep(1, 50)),
#                  y = c(sample(0:2, 20, replace = TRUE),
#                        sample(2:5, 20, replace = TRUE),
#                        rep(NA, 10),
#                        sample(5:7, 20, replace = TRUE),
#                        sample(7:9, 20, replace = TRUE),
#                        rep(NA, 10)),
#                  x1 = c(rep(0, 20), 
#                         rep(1, 20), 
#                         rep(0, 5),
#                         rep(1, 5),
#                         rep(0, 20), 
#                         rep(1, 20), 
#                         rep(0, 5),
#                         rep(1, 5)))

# in stata: mi impute regress y x1, add(5) rseed(34908) noisily

# export imputed data from stata
# exported data has 200 obs
#   100 original, including 20 obs with missing
#   5 x 20 imputed obs for missing = 100 obs

# import data from gist
# library("devtools")
# install_github("leeper/rio")

  library("rio")
  url <- "https://gist.githubusercontent.com/ericpgreen/a1a640184db13a070753288782f76d06/raw/0c61a32cdc8c1eca1ecf0d9c3e1d2f502f0cd991/stata_mi_example.R"
  df_mi <- rio::import(url)

Вывод Stata из mi estimate:reg y trt

. mi estimate:reg y trt

Multiple-imputation estimates                     Imputations     =          5
Linear regression                                 Number of obs   =        100
                                                  Average RVI     =     1.0286
                                                  Largest FMI     =     0.3865
                                                  Complete DF     =         98
DF adjustment:   Small sample                     DF:     min     =      21.82
                                                          avg     =      28.02
                                                          max     =      34.23
Model F test:       Equal FMI                     F(   1,   21.8) =      63.66
Within VCE type:          OLS                     Prob > F        =     0.0000

------------------------------------------------------------------------------
          y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         trt |   4.169307   .5225426     7.98   0.000     3.085088    5.253526
       _cons |   3.472899   .3445782    10.08   0.000     2.772806    4.172992
------------------------------------------------------------------------------

Обновление :

Вот моя попытка (1) с {mice}. Он точно не копирует результаты Stata, но я не уверен, насколько близко я должен ожидать результатов для репликации.

library(mice)
set.seed(1)
df <- data.frame(id = seq(1:100),
                 trt = c(rep(0, 50), rep(1, 50)),
                 y = c(sample(0:2, 20, replace = TRUE),
                       sample(2:5, 20, replace = TRUE),
                       rep(NA, 10),
                       sample(5:7, 20, replace = TRUE),
                       sample(7:9, 20, replace = TRUE),
                       rep(NA, 10)),
                 x1 = c(rep(0, 20),
                        rep(1, 20),
                        rep(0, 5),
                        rep(1, 5),
                        rep(0, 20),
                        rep(1, 20),
                        rep(0, 5),
                        rep(1, 5)))

df_ <- mice(df, m = 5, seed = 123, method = "norm.predict")
pred <- df_$pred
pred[, "id"] <- 0
pred[, "trt"] <- 0
pred[, "y"] <- 0

df_mi <- mice(df, m = 5, seed = 123, method = "norm.predict", 
              pred = pred)

fit <- with(df_mi, lm(y ~ trt))
pool.fit <- pool(fit)
summary(pool.fit)

#          term estimate std.error statistic       df p.value
# 1 (Intercept)     2.46 0.2302793  10.68268 96.04978       0
# 2         trt     4.08 0.3256642  12.52824 96.04978       0