r - 使用pglm的比例与二进制响应

Question

我正在处理包含对学校的几年观察的面板数据。我的DV是考试合格者的一部分，但不是正态分布的，许多DV观察值均大于0.8。因此，使用plm()(来自plm包)的面板线性模型是不合适的，所以我尝试使用pglm()(来自pglm包)来treat the DV as a binary response and use logistic regression。我统计了考生和合格者的人数。
我确定我需要对这些数据使用固定效果(单位内)估算，因为我对学校内部考试通过率的平均变化感兴趣。我有太多的观察结果无法发布完整的数据集，但这是错误消息的一个可重现的小示例:

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的方法”的替代方法，即betareg()程序包，该程序包使用beta回归] 2，但是我看不出为什么要对betareg()使用固定效果。我也可以使用glmer()并设置一个随机截距(1 | id)来运行此代码，但是考虑到我的研究问题，随机效果方法在理论上没有意义，并且Hausman测试表明我仍然需要固定效果。
我对错误消息的解释是，以某种方式重复了行名；我通过将所有行名设置为NULL来确保不是这种情况，但这不能解决问题:

row.names(data) <- NULL

在这个问题上，我也提到了看似相似的问题 such as this，但我确保在id-year配对中没有重复。
因此，对于找出错误原因的任何帮助将不胜感激。当然，也欢迎对方法论发表意见。

Answer 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)将数据从组级别(resp ortion)带到单个结果(二进制prop onse)。

# 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中)。