从聚合二项式回归生成预测

使用伯努利结果评估模型准确性相当容易，但我不确定如何从聚合二项式回归中生成有意义的预测。

以这个例子为例。我们想要模拟numCouns客户在 12 周内参加的药物咨询会议的数量（变量）：（1）他们在开始治疗之前定期使用大麻durationRegUse的年数（变量）和（2）数量他们平均每天使用的大麻克数（变量gms）。每个客户最多可以参加六次咨询会议。

这是数据

df <- data.frame(durationRegUse = c(19, 9, 13, 19, 10, 13, 2, 14, 11, 12, 7, 6, 3, 18, 17, 9, 9, 10, 0, 20, 4, 4, 8, 5, 4, 19, 25, 10, 27, 1, 10, 25, 8, 24, 8, 18, 15, 10, 6, 14, 16, 13, 4, 4, 5, 17, 13, 21, 8, 7, 10, 17, 13, 12, 28, 38, 23, 19, 36, 3, 14, 14, 22, 11, 26, 17, 4, 8, 25, 35, 14, 28, 32, 29, 22, 21, 2, 23, 35, 34, 31, 34, 15, 14, 26, 6, 3, 25, 24, 31, 31, 27, 30, 14.5, 12, 9, 3, 13, 5, 6, 23, 21, 27, 7, 36, 19, 22, 15, 11, 17, 11, 26, 21, 15),
                 gms = c(3.5, 2, 0.5, 10, 3, 3, 4, 4, 2, 2, 2, 2, 2, 2, 1, 1.75, 4, 1.75, 0.33, 5, 2.5, 1.25, 1, 0.5, 3, 2, 5, 3, 3, 0.571, 1, 0.5, 2, 4, 2.5, 1.25, 1.5, 1, 2.5, 2, 1, 2, 1.5, 2, 0.2, 1, 1, 2, 14, 2, 3.5, 3, 2, 1.75, 2, 0.55, 1, 2, 6, 0.5, 0.5, 0.5, 3, 1, 2.75, 4.5, 3, 3, 3, 2, 2, 1, 2.5, 1.75, 1, 1.5, 2, 0.7, 7, 0.5, 2, 1.2, 0.4, 3, 0.8, 1.3, 1.2, 2, 1.5, 3, 2, 2, 4, 3, 1, 6, 1, 0.5, 1.5, 2.5, 1, 2.5, 1.5, 1, 1.5, 2.5, 1.5, 2.5, 10, 1.5, 1.5, 0.5, 5, 1.5),
                 numCouns = c(6, 1, 2, 6, 0, 6, 0, 0, 2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 2, 5, 6, 0, 0, 6, 0, 6, 3, 6, 0, 0, 0, 4, 5, 0, 0, 4, 0, 4, 3, 0, 1, 2, 6, 4, 2, 4, 3, 1, 0, 2, 2, 5, 2, 0, 1, 3, 0, 3, 2, 1, 6, 0, 0, 1, 0, 1, 2, 0, 0, 5, 1, 1, 1, 5, 3, 5, 6, 6, 5, 3, 6, 2, 4, 3, 4, 6, 1, 0, 6, 4, 3, 3, 1, 5, 0, 1, 1, 6, 6, 6, 3, 3, 2, 0, 0, 5, 1, 6, 3, 0, 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()函数，选择"response"类型作为我们从伯努利响应量表（即 [0,1]）中得出 1的预测每次试验概率。

predBin <- predict(aggBinMod, type = "response")
predBin

# (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

因此，按照这个逻辑，为了从我们的聚合二项式回归模型中为每个客户生成会话数的预测，我们应该能够简单地将该值乘以我们希望预测的试验数，在我们的案例 6 中。所以为了生成我们将运行的预测

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 

从那里可以直接通过均方误差评估模型准确性

error <- predBin6 - df$numCouns
mse <- mean(error^2)
mse

# output
# [1] 4.871892

所以我的问题是这是从聚合二项式回归生成预测的正确方法吗？