使用 pglm 的比例与二进制响应答案

【问题标题】：Proportion vs. binary response with pglm使用 pglm 的比例与二进制响应
【发布时间】：2020-10-22 02:09:08
【问题描述】：

我正在处理包含多年对学校观察的面板数据。我的 DV 是考试通过者的比例，但不是正态分布的，DV 的许多观察值 > 0.8。因此，使用plm()（来自包plm）的面板线性模型是不合适的，所以我尝试使用treat the DV as a binary response and use logistic regression 和pglm()（来自包pglm）。我统计了应试者和通过者的数量。

我已确定我需要对这些数据使用固定效应（单位内）估计，因为我对学校内考试通过率的平均变化感兴趣。我有太多的观察结果无法发布完整的数据集，但这里是错误消息的一个可重复的小示例：

id <- c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4)
year <- rep(c(2017, 2018, 2019), 4)
proportion <- c(.67, .77, .79, .88, .89, .85, .79, .81, .79, .87, .75, .74)
X1 <- c(.05, .041, .037, .015, .012, .021, .081, .055, .062, .034, .031, .022)
X2 <- c(145, 146, 145, 155, 154, 154, 150, 152, 156, 148, 150, 151)
takers <- c(50, 62, 55, 112, 101, 119, 44, 45, 48, 66, 69, 60)
passers <- c(34, 48, 43, 99, 90, 101, 35, 36, 38, 57, 52, 44)
fails <- takers - passers

data <- as.data.frame(cbind(id, year, proportion, X1, X2, takers, passers, fails))

pglm::pglm(cbind(passers, fails) ~ X1 + X2, index = c("id", "year"), model = "within",  family = binomial(link = "logit"), data = data)
#> Error in `.rowNamesDF<-`(x, value = value): duplicate 'row.names' are not allowed

^{由reprex package (v0.3.0) 于 2020 年 10 月 21 日创建}

我没有遇到运行常规 logit 的问题：

glm(cbind(passers, fails) ~ X1 + X2,family = binomial(link = "logit"), data = data)

而且我也熟悉将 DV 视为二进制方法的替代方法，即使用 beta 回归的 betareg() 包]2，但我不明白为什么要使用 betareg 的固定效果()。我也可以使用 glmer() 并设置随机截距 (1|id) 运行此代码，但考虑到我的研究问题，随机效应方法在理论上没有意义，而且 Hausman 测试表明我无论如何都需要固定效应。

我对错误消息的解释是行名以某种方式重复；我通过将所有行名设置为 NULL 来确保不是这种情况，但这并没有解决问题：

row.names(data) <- NULL

我在这个问题上也提到了看似相似的问题such as this，但我确保在 id-year 配对中没有重复。

因此，非常感谢您对找出错误所在的任何帮助。当然，也欢迎对方法发表评论。

【问题讨论】：

你想用公式中的cbind(passers, fails)实现什么？你看到这个命令的结果了吗？那“通过考试的比例”怎么算？虽然glm 似乎支持这样的规范（参见?glm），但pglm 似乎不支持它。
嗨 Helix，只运行该命令会返回一个 12x2 的路人矩阵并失败。在那里，我只是试图让 pglm 函数根据观察到的那所学校的通过概率来预测考试通过。我自己不习惯这种方法，但在我的帖子链接的来源中阅读了它。如果 pglm 不支持这一点（感谢您提供的信息），您能推荐一种更好的方法吗？

标签： r error-handling plm

【解决方案1】：

有关重复行名称的错误消息有点误导，因为pglm 无法处理特定输入glm 可以处理指定比例的两列矩阵（代码中的cbind(passers, fails)）。 glm 在各种输入可能性方面更加灵活，请参阅 ?glm。

pglm 只能处理一个二进制因变量作为公式左侧的输入。因此，您希望将数据降低到“个人级别”（这里是使用 glm http://www.simonqueenborough.info/R/statistics/lessons/Binomial_Data.html;http://pages.stat.wisc.edu/~mchung/teaching/MIA/reading/GLM.logistic.Rpackage.pdf 对个人结果（二元响应）和群体结果（比例）的主题进行更好的讨论）。

这里有一些代码可以为您提供数据转换，以复制您使用glm 和pglm 估计的模型。了解如何使用考试总数 takers（passers 和 fails）将数据从组级别（proportion）降低到个人结果（二进制 response）。

# glm - your reference
summary(mod1 <- glm(cbind(passers, fails) ~ X1 + X2, family = binomial(link = "logit"), data = data))
# glm - same with weights
data$prop <- data$passers / data$takers
summary(mod2 <- glm(prop ~ X1 + X2, family = binomial(link = "logit"), data = data, weights = takers))

# construct data suitable for pglm
df2 <- df[rep(seq_along(data$takers), data$takers), ]
df2$ID <- paste(df2$id, unlist(lapply(df$takers, seq_len)), sep = '')
vec <- numeric()
for (i in 1:nrow(data)) {
    vec  <- c(vec, (c(rep(1, data$passers[i]), rep(0, data$fails[i]))))
}
df2$resp <- vec  
pdf2 <- pdata.frame(df2, index = "id")

# same with pglm
summary(mod3 <- pglm(resp ~ X1 + X2, family = binomial(link = "logit"), data = pdf2, model = "pooling"))

如果您要估计除"pooling" 模型之外的任何其他指标，您将需要构建一个不同的索引（否则您会得到错误的结果，我假设）——您可能没有相关信息（个人- pdf2/df2 中所有行的时间组合。

【讨论】：