绘制逻辑回归预测概率的置信区间

机器算法验证 r 物流 置信区间
2022-01-17 06:11:26

好的,我有一个逻辑回归,并使用该predict()函数根据我的估计绘制了一条概率曲线。

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

这很好,但我很好奇绘制概率的置信区间。我试过plot.ci()但没有运气。谁能指出我完成这项工作的一些方法,最好是使用car包或基础 R。

2个回答

您使用的代码使用该glm函数估计逻辑回归模型。你没有包括数据,所以我只是弥补一些。

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

逻辑回归模型对二元响应变量和在这种情况下的一个连续预测变量之间的关系进行建模。结果是一个对数转换的概率,作为与预测变量的线性关系。在您的情况下,结果是对应于赌博中赢或不赢的二元响应,它是通过赌注的价值来预测的。的系数mod1以对数赔率(难以解释)给出,根据:

logit(p)=log(p(1p))=β0+β1x1

要将记录的赔率转换为概率,我们可以将上述转换为

p=exp(β0+β1x1)(1+exp(β0+β1x1))

您可以使用此信息来设置绘图。首先,您需要一个预测变量范围:

plotdat <- data.frame(bid=(0:1000))

然后使用predict,您可以根据您的模型获得预测

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

请注意,拟合值也可以通过

mod1$fitted

通过指定se.fit=TRUE,您还可以获得与每个拟合值相关的标准误差。结果data.frame是一个包含以下组件的矩阵:拟合预测 ( fit)、估计的标准误差 ( ) 和一个标量,该标量给出了用于计算标准误差 ( )se.fit的离散度平方根。residual.scale在二项式 logit 的情况下,该值将为 1(您可以通过输入 来查看preddat$residual.scaleR如果您想查看到目前为止所计算的示例,您可以键入head(data.frame(preddat)).

下一步是设置情节。我喜欢先用参数设置一个空白绘图区域:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

现在您可以看到知道如何计算拟合概率的重要之处。您可以根据上面的第二个公式绘制与拟合概率相对应的线。使用 ,preddat data.frame您可以将拟合值转换为概率,并使用它来根据预测变量的值绘制一条线。

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

最后,回答您的问题,可以通过计算拟合值的概率+/- 1.96乘以标准误差来将置信区间添加到图中:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

结果图(来自随机生成的数据)应如下所示:

在此处输入图像描述

为方便起见,以下是一大块中的所有代码:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(注意:这是一个经过大量编辑的答案,试图使其与 stats.stackexchange 更相关。)

这是@smillig 解决方案的修改。我在这里使用 tidyverse 工具,也使用linkinv作为 GLM 模型对象一部分的函数mod1这样,您不必手动反转逻辑函数,并且无论您适合哪种特定 GLM,这种方法都将起作用。

library(tidyverse)
library(magrittr)


set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 


# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))