使用带有 Amelia 数据集的 Zelig“sim”函数来获得跨 R 中的插补数据集汇集的估计值

Question

我正在使用 Amelia 的多重插补数据集，然后希望 Zelig 从回归模型中计算预测值。 Zelig's documentation 声明“当绘制感兴趣的数量时，例如预期值和预测值以及一阶差分，它们会正确地汇集在来自每个 m 个估算数据集的那些中”。这是真的，但我也想获得跨每个估算数据集汇集的估计值作为“sim”命令的输出。

这里是复制 Zelig 网站上的说明并生成相同输出的示例代码：

library("Amelia")
data(africa)
a.out <- amelia(x = africa, m=5, cs = "country", ts = "year", logs = "gdp_pc")
z.out <- zelig(gdp_pc ~ trade + civlib, model = "ls", data = a.out)
summary(z.out)

然后，当“trade”设置为 50 和 100 时，我使用“setx”来估计 DV (gdp_pc) 的预测值。

x.out <- setx (z.out, trade = c(50,100))
x.out
range:
  (Intercept) trade civlib
1           1    50  0.289
2           1   100  0.289

Next step: Use 'sim' method

如果我随后使用“sim”和“plot”，R 会生成一个带有我要求的估计值的图：

s.out <- sim (z.out, x = x.out)
plot(s.out)

但是，我希望打印输出预测值及其标准误差和不同置信区间的值根据鲁宾规则汇集所有估算数据集。这不是“摘要”命令似乎正在做的事情：

summary(s.out)
[1] 50


 sim range :
 -----
ev
     mean     sd      50%     2.5%   97.5%
1 844.843 30.567 845.1218 791.8107 908.658
pv
         mean       sd      50%     2.5%    97.5%
[1,] 857.6479 372.9689 852.9239 157.7842 1553.552

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 836.2505 36.72892 833.3876 770.7931 908.7371
pv
         mean      sd      50%     2.5%    97.5%
[1,] 821.3542 359.461 790.5742 204.7687 1483.275

 sim range :
 -----
ev
     mean       sd      50%     2.5%    97.5%
1 837.307 34.99979 839.4895 765.0043 896.1513
pv
         mean       sd      50%     2.5%    97.5%
[1,] 831.6275 347.4005 844.0667 120.8968 1526.509

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 838.1396 33.49521 837.6317 776.3413 901.4235
pv
         mean       sd      50%     2.5%    97.5%
[1,] 866.5946 364.2909 830.9851 263.8757 1594.664

 sim range :
 -----
ev
     mean       sd      50%     2.5%    97.5%
1 842.784 35.18827 843.5563 779.9052 914.5869
pv
         mean       sd      50%     2.5%    97.5%
[1,] 834.7425 350.5647 834.0003 228.0261 1527.293


[1] 100


 sim range :
 -----
ev
      mean       sd      50%    2.5%    97.5%
1 1743.969 54.06692 1742.795 1627.39 1840.744
pv
        mean       sd      50%     2.5%    97.5%
[1,] 1700.53 350.1268 1718.504 1047.998 2322.216

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 1748.554 58.46152 1755.443 1634.345 1854.652
pv
         mean       sd      50%     2.5%    97.5%
[1,] 1734.831 340.8356 1734.907 1071.973 2347.156

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 1741.014 63.86164 1741.492 1615.497 1863.306
pv
         mean       sd      50%   2.5%    97.5%
[1,] 1759.305 329.6513 1746.153 1172.5 2435.067

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 1738.422 64.75221 1738.474 1615.078 1854.675
pv
         mean       sd      50%     2.5%    97.5%
[1,] 1728.152 386.8327 1761.047 849.7188 2395.825

 sim range :
 -----
ev
      mean       sd      50%     2.5%    97.5%
1 1746.575 53.02558 1744.919 1638.602 1848.114
pv
         mean       sd      50%    2.5%    97.5%
[1,] 1710.864 342.1865 1702.769 1050.85 2288.021

在这里，我得到了每个插补数据集的所有值，而不是所有多重插补数据集的汇总值。在提供预测估计的汇总统计数据以及基于它们绘制图表时，是否有办法让 Zelig 将 Rubin 规则应用于多重插补数据集？

注意：我需要的应用程序需要 负二项式回归，而不是线性回归，作为 Zelig 中使用的模型。我已使用此示例复制了 Zelig 开发人员提供的示例。

非常感谢您的帮助，祝您有美好的一天！

Answer 1

在这种情况下，您不需要使用鲁宾规则，因为不确定性是根据模拟中的方差计算得出的。 Zelig 没有为你平均这些，我有点惊讶，但你可以自己做，没有太多困难：

qi.out <- zelig_qi_to_df(s.out)

lapply(split(qi.out, qi.out["trade"]),
       function(x) c(trade = unique(x$trade),
                     mean = mean(x$expected_value),
                     sd = sd(x$expected_value),
                     median = median(x$expected_value),
                     quantile(x$expected_value, probs = c(0.5, 0.025, 0.975))))

lapply(split(qi.out, qi.out["trade"]),
       function(x) c(trade = unique(x$trade),
                     mean = mean(x$predicted_value),
                     sd = sd(x$predicted_value),
                     median = median(x$predicted_value),
                     quantile(x$predicted_value, probs = c(0.5, 0.025, 0.975))))