问题描述
使用伯努利结果评估模型准确性相当容易,但我不确定如何从聚合二项式回归中生成有意义的预测。
以这个例子为例。我们希望对客户在 12 周内参加的药物咨询次数(变量 numCouns
)进行建模,其基础是:(1) 在开始治疗之前他们定期使用大麻的年数(变量 {{ 1}}) 和 (2) 他们平均每天使用的大麻克数(变量 durationRegUse
)。每个客户最多可以参加六次咨询会议。
这是数据
gms
要将其建模为聚合二项式回归,我们需要创建一个覆盖变量(最大会话数)。
df <- data.frame(durationRegUse = c(19,9,13,19,10,2,14,11,12,7,6,3,18,17,20,4,8,5,25,27,1,24,15,16,21,28,38,23,36,22,26,35,32,29,34,31,30,14.5,15),gms = c(3.5,0.5,1.75,0.33,2.5,1.25,0.571,1.5,0.2,3.5,0.55,2.75,4.5,0.7,1.2,0.4,0.8,1.3,1.5),numCouns = c(6,0))
现在我们可以创建聚合二项式回归模型
df$coverage <- 6
这是输出
aggBinMod <- glm(
formula = cbind(numCouns,coverage - numCouns) ~ durationRegUse + gms,data = df,family = binomial(link = "logit"))
现在是我不确定的部分:如何生成用于评估模型准确性的预测。现在,据我所知,如果我们使用 summary(aggBinMod)
#output
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.157570 0.183116 -6.322 2.59e-10 ***
# durationRegUse 0.035975 0.008455 4.255 2.09e-05 ***
# gms 0.075838 0.039273 1.931 0.0535 .
函数,选择 predict()
作为我们得到预测的每次试验概率从伯努利响应量表中抽取 1 的概率(即 [0,1])。
"response"
predBin <- predict(aggBinMod,type = "response")
predBin
因此,按照这个逻辑,为了从我们的聚合二项式回归模型中为每个客户生成会话数的预测,我们应该能够简单地将此值乘以我们希望预测的试验数,在我们的案例 6. 所以为了生成我们将运行的预测
# (predicted bernoulli probability for first 16 participants)
# 1 2 3 4 5 6 7 8
# 0.4480346 0.3357882 0.3425441 0.5706073 0.3611657 0.3864206 0.3138308 0.4132440
# 9 10 11 12 13 14 15 16
# 0.3520203 0.3602692 0.3199350 0.3121589 0.2894678 0.4113600 0.3845787 0.3315728
从那里可以直接通过均方误差评估模型准确性
predBin6 <- predict(aggBinMod,type = "response")*6
predBin6
# predicted number of sessions,out of a possible 6),for first 18 clients
# 1 2 3 4 5 6 7 8 9
# 2.688208 2.014729 2.055265 3.423644 2.166994 2.318524 1.882985 2.479464 2.112122
# 10 11 12 13 14 15 16 17 18
# 2.161615 1.919610 1.872954 1.736807 2.468160 2.307472 1.989437 2.222478 2.037563
解决方法
或多或少,是的。
我建议不要硬编码每次观察有 6 次试验(在某些应用中,试验次数因观察而异),而是建议
predBin6 <- predict(aggBinMod,type = "response")*weights(aggBinMod)
(在您的情况下应该给出相同的答案)。
我还要说 MSE 是合理的,但不一定是二项式模型预测准确性的最佳衡量标准(它没有考虑方差对均值的依赖性)。 (我没有特别的替代建议,但偏差 (deviance(aggBinMod)
) 或类似的东西可能是合适的。)