从逻辑回归模型拟合中获取预测值(Y=1 或 0)

机器算法验证 r 广义线性模型 物流
2022-02-08 08:22:54

假设我有一个类对象glm(对应于逻辑回归模型),我想将predict.glm使用参数给出的预测概率type="response"转换为二元响应,即Y=1或者Y=0. 在 R 中执行此操作的最快和最规范的方法是什么?

虽然我再次知道predict.glm,但我不知道截止值到底在哪里P(Yi=1|X^i)生活——我想这是我的主要绊脚石。

1个回答

一旦你有了预测的概率,你就可以自己决定使用什么阈值了。您可以选择阈值来优化灵敏度、特异性或在应用程序上下文中最重要的任何度量(一些附加信息将有助于获得更具体的答案)。您可能想查看 ROC 曲线和其他与最佳分类相关的度量。

编辑:为了稍微澄清这个答案,我将举一个例子。真正的答案是最佳截止值取决于分类器的哪些属性在应用程序的上下文中是重要的。Yi成为观察的真实值i, 和Y^i是预测的类。一些常见的绩效衡量标准是

(1) 灵敏度:P(Y^i=1|Yi=1)- 正确识别为“1”的比例。

(2) 特异性:P(Y^i=0|Yi=0)- 正确识别为“0”的比例

(3)(正确)分类率:P(Yi=Y^i)- 正确预测的比例。

(1) 也称为 True Positive Rate,(2) 也称为 True Negative Rate。

例如,如果您的分类器旨在评估具有相对安全治愈方法的严重疾病的诊断测试,则敏感性远比特异性重要。在另一种情况下,如果疾病相对较轻并且治疗有风险,那么特异性控制将更重要。对于一般的分类问题,联合优化灵敏度和规范被认为是“好的” - 例如,您可以使用最小化它们与点的欧几里得距离的分类器(1,1)

δ=[P(Yi=1|Y^i=1)1]2+[P(Yi=0|Y^i=0)1]2

δ可以以另一种方式加权或修改,以反映更合理的距离测量(1,1)在应用程序的上下文中 - 出于说明目的,此处任意选择了 (1,1) 的欧几里德距离。在任何情况下,所有这四种措施都可能是最合适的,具体取决于应用程序。

下面是一个模拟示例,使用逻辑回归模型的预测进行分类。改变截止值以查看在这三个度量中的每一个下哪个截止值给出了“最佳”分类器。在此示例中,数据来自具有三个预测变量的逻辑回归模型(请参见下图的 R 代码)。正如您从这个示例中看到的,“最佳”截止值取决于这些措施中哪一个最重要——这完全取决于应用程序。

编辑2: P(Yi=1|Y^i=1)P(Yi=0|Y^i=0),阳性预测值和阴性预测值(注意这些与敏感性和特异性不同)也可能是有用的性能度量。

在此处输入图像描述

# data y simulated from a logistic regression model 
# with with three predictors, n=10000
x = matrix(rnorm(30000),10000,3)
lp = 0 + x[,1] - 1.42*x[2] + .67*x[,3] + 1.1*x[,1]*x[,2] - 1.5*x[,1]*x[,3] +2.2*x[,2]*x[,3] + x[,1]*x[,2]*x[,3]
p = 1/(1+exp(-lp))
y = runif(10000)<p

# fit a logistic regression model
mod = glm(y~x[,1]*x[,2]*x[,3],family="binomial")

# using a cutoff of cut, calculate sensitivity, specificity, and classification rate
perf = function(cut, mod, y)
{
   yhat = (mod$fit>cut)
   w = which(y==1)
   sensitivity = mean( yhat[w] == 1 ) 
   specificity = mean( yhat[-w] == 0 ) 
   c.rate = mean( y==yhat ) 
   d = cbind(sensitivity,specificity)-c(1,1)
   d = sqrt( d[1]^2 + d[2]^2 ) 
   out = t(as.matrix(c(sensitivity, specificity, c.rate,d)))
   colnames(out) = c("sensitivity", "specificity", "c.rate", "distance")
   return(out)
}

s = seq(.01,.99,length=1000)
OUT = matrix(0,1000,4)
for(i in 1:1000) OUT[i,]=perf(s[i],mod,y)
plot(s,OUT[,1],xlab="Cutoff",ylab="Value",cex.lab=1.5,cex.axis=1.5,ylim=c(0,1),type="l",lwd=2,axes=FALSE,col=2)
axis(1,seq(0,1,length=5),seq(0,1,length=5),cex.lab=1.5)
axis(2,seq(0,1,length=5),seq(0,1,length=5),cex.lab=1.5)
lines(s,OUT[,2],col="darkgreen",lwd=2)
lines(s,OUT[,3],col=4,lwd=2)
lines(s,OUT[,4],col="darkred",lwd=2)
box()
legend(0,.25,col=c(2,"darkgreen",4,"darkred"),lwd=c(2,2,2,2),c("Sensitivity","Specificity","Classification Rate","Distance"))