一位审稿人要求我评估逻辑回归结果的稳健性,因为估计可能会因结果中的类别不平衡而产生偏差。
为了具体化,我运行了三个不同的模型,结果都有大约 25% 的概率(我什至不会说这是不平衡的)。目的是分析结果与一些协变量之间的关联,而不是预测未来的观察结果。
我看到了相关的问题,但没有一个帮助我为这个论点制定答案。我应该反驳它,说结果不是真的不平衡吗?或者,正如@whuber 所写,指出
逻辑回归往往效果很好,即使结果不平衡,也会给出合理接近正确参数的值。
任何提示或参考表示赞赏。
谢谢!
一位审稿人要求我评估逻辑回归结果的稳健性,因为估计可能会因结果中的类别不平衡而产生偏差。
为了具体化,我运行了三个不同的模型,结果都有大约 25% 的概率(我什至不会说这是不平衡的)。目的是分析结果与一些协变量之间的关联,而不是预测未来的观察结果。
我看到了相关的问题,但没有一个帮助我为这个论点制定答案。我应该反驳它,说结果不是真的不平衡吗?或者,正如@whuber 所写,指出
逻辑回归往往效果很好,即使结果不平衡,也会给出合理接近正确参数的值。
任何提示或参考表示赞赏。
谢谢!
让我们来了解一下。
首先,平衡数据集会发生什么?
这是一个包含观测值的数据集的散点图,其中为零,其余为 1。在它上面,我绘制了基础(“真实”)概率和符合逻辑回归的概率。这两张图非常吻合,表明逻辑回归在这种情况下做得很好。
为了更好地理解它,我保留了相同的值,但随机次。在这个逻辑回归中, )的估计。这是这些估计的散点图。
中间的红色三角形绘制了真实系数 椭圆是该点云的二阶近似:一个用于包含大约一半的点,另一个用于包含大约 95% 的点。他们这样做表明他们给出了这样一个数据集的任何给定估计可能有多不确定性的可靠指示:截距可能偏离大约(外椭圆的宽度),斜率可能偏离大约(外椭圆的高度)。这些是误差范围。
不平衡的数据集会发生什么?
这是一个类似的设置,但一个类中只有的分数(给予或拿几分,取决于进行这些观察所涉及的随机性):
(确实很小:它告诉我们期望在一个类中只看到左右的值,而在另一个类中看到另外。)
现在的拟合明显偏离了真实的图表——但是逻辑回归的这个证据是不是“稳健”的?同样,我们可以通过多次重复生成随机数据和估计拟合的过程来找出答案。这是估计值的散点图。
总的来说,估计值接近 从这个意义上说,逻辑回归看起来“稳健”。(我保持与以前相同的斜率并调整截距以降低观察的比率。)
误差范围发生了变化:(某种程度上)描述截距不确定性的外椭圆的半轴已从增长到超过,而另一个半轴从缩小到大约整个画面已经倾斜。
椭圆不再像以前那样很好地描述点云:现在,逻辑回归有一些趋势来估计极负的斜率和截距。倾斜表明估计之间存在明显的相关性:低(负)截距往往与低负斜率相关(通过预测值来补偿小截距。 )但是这种相关性是可以预料的:这只要数据的平均值不靠近垂直轴,它看起来就像普通的最小二乘回归。
这些实验表明了什么?
对于这种大小(或更大)的数据集,至少:
逻辑回归往往效果很好,即使结果不平衡,也会给出合理接近正确参数的值。
参数估计(逻辑回归的常规输出)之间相关性的二阶描述并不能完全捕捉到估计可能同时与事实相距甚远的可能性。
一个元结论
您可以通过在根据已知现实模型生成的数据上重复运行它并跟踪对您很重要的输出来评估任何过程(例如逻辑回归)的“稳健性”(或更一般地说,显着的统计属性)。
这是R产生这些数字的代码。对于前两个数字,第一行更改为p <- 50/100. 删除set.seed调用以生成其他随机示例。
尝试像这样的模拟(扩展到更多解释变量)可能会让您相信标准经验法则的效用:
让较小类别中的观察数量指导模型的复杂性。
而在普通最小二乘回归中,我们可能对每个解释变量有 10 个观察值(总计)感到满意,而对于逻辑回归,我们希望每个解释变量在较小的类别中有 10 个观察值。
p <- 5/100 # Proportion of one class
n <- 200 # Dataset size
x <- seq(-1, 1, length.out=n) # Explanatory variable
beta <- -3 # Slope
#
# Find an intercept that yields `p` as the true proportion for these `x`.
#
logistic <- function(z) 1 - 1/(1 + exp(z))
alpha <- uniroot(function(a) mean(logistic(a + beta*x)) - p, c(-5,5))$root
#
# Create and plot a dataset with an expected value of `p`.
#
set.seed(17)
y <- rbinom(n, 1, logistic(alpha + beta*x))
plot(range(x), c(-0.015, 1.015), type="n", bty="n", xlab="x", ylab="y",
main="Data with True (Solid) and Fitted (Dashed) Probabilities")
curve(logistic(alpha + beta*x), add=TRUE, col="Gray", lwd=2)
rug(x[y==0], side=1, col="Red")
rug(x[y==1], side=3, col="Red")
points(x, y, pch=21, bg="#00000020")
#
# Fit a logistic model.
#
X <- data.frame(x=x, y=y)
fit <- glm(y ~ x, data=X, family="binomial")
summary(fit)
#
# Sketch the fit.
#
b <- coefficients(fit)
curve(logistic(b[1] + b[2]*x), add=TRUE, col="Black", lty=3, lwd=2)
#
# Evaluate the robustness of the fit.
#
sim <- replicate(500, {
X$y.new <- with(X, rbinom(n, 1, logistic(alpha + beta*x)))
coefficients(glm(y.new ~ x, data=X, family="binomial"))
})
plot(t(sim), main="Estimated Coefficients", ylab="")
mtext(expression(hat(beta)), side=2, line=2.5, las=2, cex=1.25)
points(alpha, beta, pch=24, bg="#ff0000c0", cex=1.6)
#
# Plot second moment ellipses.
#
V <- cov(t(sim))
obj <- eigen(V)
a <- seq(0, 2*pi, length.out=361)
for (level in c(.50, .95)) {
lambda <- sqrt(obj$values) * sqrt(qchisq(level, 2))
st <- obj$vectors %*% (rbind(cos(a), sin(a)) * lambda) + c(alpha, beta)
polygon(st[1,], st[2,], col="#ffff0010")
}