从 R 中的二项式 logit 输出计算累积概率

Question

我有一个拟合的二项式 logit 模型，并且想要计算经历事件的累积概率

例如，如果我有一个基于连续距离范围（0-8.5 公里）预测和结果的拟合模型，我可能想找出距离

我有如下估计概率和相关距离的向量

dat <- structure(list(km = c(0, 0.447368421052632, 0.894736842105263, 
1.34210526315789, 1.78947368421053, 2.23684210526316, 2.68421052631579, 
3.13157894736842, 3.57894736842105, 4.02631578947368, 4.47368421052632, 
4.92105263157895, 5.36842105263158, 5.81578947368421, 6.26315789473684, 
6.71052631578947, 7.15789473684211, 7.60526315789474, 8.05263157894737, 
8.5), prob = c(0.99010519543441, 0.985413663823809, 0.97854588563623, 
0.968547716962174, 0.954108659036907, 0.933496091194704, 0.904551377544634, 
0.864833064332603, 0.81202174997839, 0.744668375529677, 0.663191827576796, 
0.570704402277059, 0.47300143764816, 0.377323442817887, 0.290336664745317, 
0.216433162546689, 0.157174982015906, 0.111825887625402, 0.0783449309507567, 
0.054275681518511)), .Names = c("km", "prob"), row.names = c(NA, 
-20L), class = "data.frame")

我最终想说的是“x 距离内的 x% 的观察预计会经历一个事件”。这是正确的做法吗？

还有一种简单的方法可以计算在哪个距离（从 0 到任意）包含 50% 的累积概率。

谢谢，蒂姆

Answer 1

可能有一些方法可以从您的模型中提取它，但如果您是从头开始，我会尝试将您的数据拟合到一个分布中，然后提取您的相关数据点。

首先定义一个误差函数：

rmse <- function(x,y) sqrt(sum((x-y)^2)/length(x)) # or some other error fxn

现在假设您的数据看起来像伽马分布，那么试试：

gdf <- function(x, d=dat$km) pgamma(d,shape=x[1], scale=x[2])

因此，您要优化的函数将是数据的误差函数和拟合分布：

error_fxn <- function(x) rmse(rev(dat$prob),gdf(x)) # rev data to make ascending

现在优化此函数以获取您感兴趣的分布参数：

rr <- optim(c(1,1),error_fxn)

让我们看看合身有多好（还可以……）；

rr
# $par
# [1] 3.108392 1.112584

# $value
# [1] 0.0333369

# $counts
# function gradient 
     119       NA 

# $convergence
# [1] 0

# $message
# NULL

或图形化：

with(dat,plot(km,prob,xlim=c(10,0)))
with(dat,lines(rev(km),pgamma(km,shape=rr$par[1], scale=rr$par[2]),col='red'))

查看 CDF 的值：

kms <- seq(0,8.5,0.5)
data.frame(dist = kms, cdf = pgamma(kms,shape=rr$par[1], scale=rr$par[2]))

#       dist         cdf
#    1   0.0 0.000000000
#    2   0.5 0.008634055
#    3   1.0 0.053615340
#    4   1.5 0.137291689
#    5   2.0 0.245961242
#    6   2.5 0.363956061
#    7   3.0 0.479070721
#    8   3.5 0.583659363
#    9   4.0 0.673982194
#    10  4.5 0.749075757
#    11  5.0 0.809691054
#    12  5.5 0.857478086
#    13  6.0 0.894431622
#    14  6.5 0.922551998
#    15  7.0 0.943661710
#    16  7.5 0.959325076
#    17  8.0 0.970830577
#    18  8.5 0.979207658

为了回答你的最后一个问题，距离为 CDF 的 50%：

qgamma(0.5,shape=rr$par[1], scale=rr$par[2])
# [1] 3.095395