提问人：Agata 提问时间：9/8/2023 最后编辑：Agata 更新时间：9/10/2023 访问量：88
predict（） 函数生成的值与 GLMER 中的手动计算不同
predict() function produces different values than hand calculation in a glmer
问：
我正在尝试从 glmer 模型中获取数据源（编码为 0 或 1，对于源 A 和源 B）的预测概率。 使用示例数据：
set.seed(123)
n<-7052
Df <- data.frame(
source = sample(c(0, 1), n, replace = TRUE,
prob = c(0.719, 0.221)),
Response.number = sample(1:20, n, replace = TRUE),
Item.number = sample(1:40, n, replace = TRUE),
Ps.number = sample(1:40, n, replace = TRUE)
)
Model1 <- glmer(source ~ (1|Response.number/Item.number) +
(1|Ps.number),
data=Df, family = binomial,
glmerControl(optimizer="bobyqa"))
根据 https://sebastiansauer.github.io/convert_logit2prob/，手计算产生的预测概率与以下函数相同：(exp(b)/(1+(exp(b))
probability <- predict(Model1, type="response")
mean(probability)
我尝试了多种类型的练习数据，这通常有效（在上面的示例中，它是 0.23199）。但是，当我使用实际数据时，我从预测函数 （0.59） 获得的值与手动 （0.57） 的值略有不同。我知道这不是很多，但是当我使用任何其他数据时，不会发生这种差异。
head(Df_real)
source Response.number Item.number Ps.number
0 1 1 1
0 2 1 1
1 3 1 1
1 4 1 1
0 5 1 1
0 6 1 1
0 1 2 1
0 2 2 1
1 3 2 1
1 4 2 1
0 5 2 1
0 6 2 1
0 1 1 2
0 2 1 2
1 3 1 2
1 4 1 2
0 5 1 2
0 6 1 2
0 1 2 2
0 2 2 2
1 3 2 2
1 4 2 2
0 5 2 2
等。
数据是嵌套的，也就是说，每个响应值的参与者数量大致相同，每个项目值的响应数量相同，依此类推。这可能是差异的根源吗？如果是这样，如何处理？功能合适吗？predict()
答：
3赞
Allan Cameron
9/8/2023
#1
当您在 中运行时，它使用原始数据中存在的变量（包括随机效应）来估计概率，因此您不会返回与通过运行固定效应系数获得的单个值相同的值向量。predictglmerpredictexp(b)/(1 + exp(b))
为了看到这一点，让我们尝试将随机效应变量的小数据框传递给以下参数：newdatapredict
predict(Model1, newdata = data.frame(Item.number = 1,
Response.number = c(1, 2),
Ps.number = 1), type = 'response')
#> 1 2
#> 0.2261900 0.2405297
由于模型中没有任何固定效应，因此总体概率（考虑随机效应）为：
b <- fixef(Model1)
exp(b)/(1 + exp(b))
#> (Intercept)
#> 0.2319048
正如 Ben Bolker 在评论中指出的那样，由于 glmms 中使用的偏差调整，这与数据中的原始比例不同。他还指出，我们可以从使用 中删除随机效应，这将为您提供与转换后的截距相同的值：predictre.form = NA
mean(predict(Model1, type= 'response', re.form = NA)) == plogis(fixef(Model1))
#> (Intercept)
#> TRUE
所以这实际上取决于你想预测什么，即你是否希望考虑随机变量。如果这样做，你可以使用 ，否则你可以从固定效果中手工计算或使用内部predictre.form = NApredict
顺便说一句，基本 R 函数可能是将对数赔率转换为概率的最简单方法，它在这里显然有效 - 我们可以看到 using 等价于plogistype = "response"plogis(predict(Model1, type = "link"))
all(
plogis(predict(Model1, type = "link")) == predict(Model1, type = "response")
)
#> [1] TRUE
手动计算是可以的，尽管你会得到非常小的浮点差异：
b <- predict(Model1, type = "link")
hist(exp(b)/(1 + exp(b)) - predict(Model1, type = 'response'))
因此，从模型中手动计算总体概率的明智方法是
plogis(fixef(Model1))
#> (Intercept)
#> 0.2319048
评论
0赞
Ben Bolker
9/8/2023
这是一个很好的答案，但有一个非常重要的区别。如果我们运行 GLM（无随机效应），则数据的平均值将与反向变换的截距相同：。但是，对于 GLMM 来说，情况并非如此;例如，请参阅 cran.r-project.org/web/packages/emmeans/vignettes/...set.seed(101); x <- rbinom(100, size = 1, prob = 0.2); g <- glm(x ~ 1, family = binomial); all.equal(unname(plogis(coef(g))), mean(x))
1赞
Ben Bolker
9/8/2023
此外，如果你想忽略随机效果，你可以使用predict(model, re.form = NA)
0赞
Allan Cameron
9/8/2023
@BenBolker 我认为这基本上是OP正在寻找的答案。 与mean(predict(Model1, type= 'response', re.form = NA))plogis(fixef(Model1))
0赞
Agata
9/8/2023
谢谢！“re.form = NA”部分是 R 和我的手计算之间的区别 - 但我认为我不想忽略随机效应，它们是我在模型中的唯一效应。那么，即使它与手计算不匹配并将差异归因于随机效应，也可以使用 R 中的值吗？
0赞
Allan Cameron
9/8/2023
@Agata有什么区别呢？由 给出的概率输出的平均值 ？predict
0赞
Ben Bolker
9/10/2023
#2
我认为在比较数据均值和预测均值时，您可能会遗漏一个重要的点。@AllanCameron的评论与（这是詹森的不等式）不同。plogis(mean(predict(model)))mean(plogis(predict(model)))
library(lme4)
library(emmeans)
set.seed(123)
n <- 7052
Df <- data.frame(
Response.number = sample(1:20, n, replace = TRUE),
Item.number = sample(1:40, n, replace = TRUE),
Ps.number = sample(1:40, n, replace = TRUE)
)
Df$source <- simulate(~(1|Response.number/Item.number) + (1|Ps.number),
family = binomial,
newdata = Df,
newparams = list(beta = qlogis(0.7), theta = c(1, 1, 1)))[[1]]
fit <- glmer(source ~(1|Response.number/Item.number) + (1|Ps.number),
family = binomial,
data = Df)
mean(Df$source) ## 0.6498866
(p1 <- predict(fit, newdata = data.frame(dummy = 1), re.form = NA)) ## 0.9280772
plogis(p1) ## 0.716685
(p2 <- predict(fit, newdata = data.frame(dummy = 1),
re.form = NA, type = "response")) ## 0.716685
emmeans(fit, ~ 1)
emmeans(fit, ~ 1)
## 1 emmean SE df asymp.LCL asymp.UCL
## overall 0.928 0.263 Inf 0.412 1.44
emmeans(fit, ~ 1, type = "response")
## 1 prob SE df asymp.LCL asymp.UCL
## overall 0.717 0.0535 Inf 0.602 0.809
来自 emmeans 小插图
vars <- sapply(VarCorr(fit), c)
total.SD <- sqrt(sum(vars^2))
emmeans(fit, ~ 1, type = "response", bias.adj = TRUE,
sigma = total.SD)
## 1 prob SE df asymp.LCL asymp.UCL
## overall 0.614 0.0398 Inf 0.545 0.698
偏差校正并不精确（它使用增量方法近似），所以这不太正确，但它更接近。
这稍微好一点：
library(logitnorm)
momentsLogitnorm(mu = fixef(fit), sigma = total.SD)
## mean var
## 0.65790176 0.06472473
艺术
mean(predict(fit, type = "response")) ## 0.6500409
评论
plogis(x)exp(x)/(1+exp(x))xxpredictb