R 中的 xtlogit - logit 固定效果答案

【问题标题】：xtlogit in R - logit fixed effectsR 中的 xtlogit - logit 固定效果
【发布时间】：2018-07-31 07:25:57
【问题描述】：

我试图弄清楚如何在 R 中执行固定效应 logit 回归（类似于 Stata 的 xtlogit 命令）。我阅读了几个包，例如“pglm”或“bife”，但无法让我的模型运行。

我的数据保存为数据框，如下所示：

ID   Time  Y  X
1    2000  1  0
1    2001  0  1
1    2002  1  1
...
1    2016  1  0
...
n

基本上我想运行固定效应 logit 回归：

 y_jt = beta*x_jt + mu_j + pi_t + epsilon_jt

其中 j 是 ID，t 是时间，mu 是 ID 固定效应，pi 是时间固定效应，epsilon 是误差项。

我愿意使用任何软件包。我从“bife”开始，但不知道如何设置 ID 和时间固定效果。我试过了：

 mod.no <- bife(y ~ x | ID, data = panel)

我是否可能需要像 Stata 的“xtset”命令一样将数据设置为面板？

非常感谢！

编辑

我想在 R 中复制的 Stata 命令是：

xi: xtlogit Y X i.Time, fe

【问题讨论】：

标签： r binary panel

【解决方案1】：

总的来说，我认为这里的策略是执行以下操作：

创建一个包含预测变量的个体水平平均值的变量。使用dplyr 最容易做到这一点：

data <- data $>$ group_by(ID) %>% mutate(X_mean = mean(X))

这里的神奇之处在于group_by 函数，它使mean 运算计算组均值而不是全局均值。

使用lme4 将 logit 模型估计为多级模型。以下是我指定模型的方式：

glmer(Y ~ X + X_mean + Time + (1 | ID), family = binomial)

“固定”和“随机”这两个术语在面板数据、多级建模和其他一些文献之间确实混为一谈，所以我不太清楚您如何概念化“时间的固定效应”。这个模型给你的是X 的固定效应，因为X 的系数将代表X 的主体内效应。我将Time 作为一个预测变量，它将年份视为一个额外的预测变量，其解释取决于它是连续的还是分类的。有些人会将其视为“随机”效应（如随机斜率或某些文献中的“增长曲线”）。你可以这样做：

glmer(Y ~ X + X_mean + Time + (Time | ID), family = binomial)

它估计每个人的不同时间影响。

第一个模型中的(1 | ID) 和第二个模型中的(Time | ID) 告诉lme4 分组变量是什么，在您的情况下是ID。您在第一个模型中获得ID 的随机截距，在第二个模型中获得Time 的随机截距和随机斜率。对您的第一篇文章的另一种解释是，您也希望随机截取Time，在这种情况下，您可以执行以下操作：

glmer(Y ~ X + X_mean + (1 | Time) + (1 | ID), family = binomial)

或者，或者，如果波浪很少，您可以通过将Time 包含为预测变量并将该变量作为输入数据中的一个因素来到达同一个地方。如果有很多时间点会使输出变得笨拙。

受来自 Stata 的 xt 套件的启发，我一直在开发一个包来自动化其中的一些内容，尽管此时我的包受到了更多限制。它叫panelr，目前必须从GitHub下载。更多信息在这里：https://github.com/jacob-long/panelr

在这种情况下，使用panelr，您的情况将如下所示：

library(panelr)
pdata <- panel_data(data, id = ID, wave = Time)
model <- wbm(Y ~ X, data = pdata, use.wave = TRUE, family = binomial)

panelr 所做的一切都是自动化我上面解释过的事情。您可以使用 model = "within" 参数删除单个平均变量，而不会影响对 X 的主体内效应的估计。

panelr 此时距离 CRAN 提交可能还有几周的时间，因为有些事情需要记录，有一些边缘情况会意外中断，我希望在时间处理上更加灵活。

【讨论】：

非常感谢。我会试试看。这对我来说似乎相当复杂。有没有办法用 pglm 或 bife 包做到这一点？

【解决方案2】：

也许尝试将 Time 变量强制转换为一个因子，就像您在 Stata 中所做的那样，并将其包含在 bife 公式中：

panel$Time <– as.factor(panel$Time)

mod.no <- bife(y ~ x + Time | ID, data = panel, bias_corr = "ana")

注意：您可能需要在末尾使用bias_corr = "ana" 来纠正附带的参数偏差。

或者，使用 survival 包中的 clogit 函数也应该可以工作：

mod.no <– clogit(y ~ x + strata(Time) + strata(ID), data = panel

strata() 选项表示固定效果。

【讨论】：