本文介绍了计算成功/失败总数的二项式GLMM的brier分数的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
当响应变量为成功/失败的cbind时，如何计算具有二项式误差的广义线性混合效应模型的brier分数?
我已经尝试过DescTools和scoring软件包，但是我在理解BrierScore和brierscore功能的帮助文件时遇到了麻烦...
这是到目前为止我对一些示例数据所做的工作.
#Data
df <- structure(list(Dose = c(20, 14, 14, 20, 0, 0, 14, 14, 14, 16,
10, 20, 20, 20, 16, 10, 10, 0, 16, 16, 16, 10, 0, 12, 10, 12,
12, 0, 0, 20, 12, 16, 10, 12, 12, 0, 14, 14, 16, 0, 14, 20, 16,
20, 14, 12, 12, 20, 20, 0, 0, 14, 12, 10, 10, 20, 16, 16, 14,
10, 10, 10, 20, 16, 10, 0, 12, 12, 0, 12, 16, 14, 16, 14, 0,
0, 12, 20, 0, 12, 14, 14, 0, 0, 20, 20, 20, 14, 14, 10, 10, 20,
16, 16, 0, 12, 10, 10, 10, 16, 16, 12, 20, 10, 12, 12, 16, 14,
0, 16, 20, 12, 14, 10, 10, 0, 0, 12, 12, 10, 10, 0, 0, 0, 14,
12, 12, 20, 20, 14, 14, 14, 12, 20, 20, 20, 16, 16, 14, 10, 10,
16, 16, 16), Success = c(100, 91, 87, 100, 0, 0, 91, 96, 89,
96, 82, 99, 99, 99, 92, 59, 45, 0, 100, 95, 100, 83, 8, 82, 63,
98, 74, 9, 0, 99, 78, 98, 53, 96, 52, 0, 62, 85, 98, 4, 89, 99,
99, 97, 82, 80, 91, 99, 96, 0, 0, 95, 80, 68, 74, 100, 97, 93,
87, 34, 32, 47, 99, 96, 86, 15, 93, 86, 0, 77, 89, 80, 98, 96,
31, 0, 61, 100, 0, 84, 88, 97, 0, 0, 99, 100, 100, 92, 88, 46,
51, 99, 97, 100, 0, 93, 61, 91, 57, 76, 95, 50, 98, 16, 87, 93,
87, 88, 11, 92, 98, 60, 96, 0, 64, 72, 0, 74, 77, 0, 62, 0, 0,
0, 84, 47, 69, 98, 100, 61, 90, 79, 11, 100, 98, 100, 98, 93,
91, 58, 58, 93, 78, 69), Failure = c(0, 9, 13, 0, 100, 100, 9,
4, 11, 4, 18, 1, 1, 1, 8, 41, 55, 100, 0, 5, 0, 17, 92, 18, 37,
2, 26, 91, 100, 1, 22, 2, 47, 4, 48, 100, 38, 15, 2, 96, 11,
1, 1, 3, 18, 20, 9, 1, 4, 100, 100, 5, 20, 32, 26, 0, 3, 7, 13,
66, 68, 53, 1, 4, 14, 85, 7, 14, 100, 23, 11, 20, 2, 4, 69, 100,
39, 0, 100, 16, 12, 3, 100, 100, 1, 0, 0, 8, 12, 54, 49, 1, 3,
0, 100, 7, 39, 9, 43, 24, 5, 50, 2, 84, 13, 7, 13, 12, 89, 8,
2, 40, 4, 100, 36, 28, 100, 26, 23, 100, 38, 100, 100, 100, 16,
53, 31, 2, 0, 39, 10, 21, 89, 0, 2, 0, 2, 7, 9, 42, 42, 7, 22,
31), Rep = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L)), row.names = 433:576, class = "data.frame")
#Model fitting
fit <- lme4::glmer(cbind(Success, Failure) ~ Dose + (1|factor(Rep)), family = binomial, data = df)
#Brier scores
DescTools::BrierScore(fit)
#Returns following error
#Error in mean(resp * (1 - pred)^2 + (1 - resp) * pred^2)
scoring::brierscore(fit)
#Returns following error
#argument "pred" is missing, with no default
是否可以为混合效应模型计算刺骨分数?如果可以，怎么办?
解决方案
我想您可以简单地摆脱cbind结构，获取长数据格式，并使用DescTools :: BrierScore()的预测/响应接口. /p>
类似的事情可能会起作用:
library(DescTools)
tt <- aggregate(df[, c("Success", "Failure")], by=as.list(df[, c("Dose", "Rep")]), FUN = sum)
d.frm <- Untable(reshape(tt, idvar=c("Dose","Rep"), varying = c("Success", "Failure"),
direction = "long",
timevar="Success", v.names="Freq", times=c(1, 0)))
fit <- lme4::glmer(Success ~ Dose + (1|Rep), family = binomial, data = d.frm)
BrierScore(pred=predict(fit, type="response"), resp = d.frm$Success)
How do I calculate brier scores for generalised linear mixed-effects model with binomial errors when my response variable is a cbind of success/failures?
I have tried the DescTools and scoring packages but I am having trouble understanding the help files for the BrierScore and brierscore functions...
Here is what I have done so far with some sample data.
#Data
df <- structure(list(Dose = c(20, 14, 14, 20, 0, 0, 14, 14, 14, 16,
10, 20, 20, 20, 16, 10, 10, 0, 16, 16, 16, 10, 0, 12, 10, 12,
12, 0, 0, 20, 12, 16, 10, 12, 12, 0, 14, 14, 16, 0, 14, 20, 16,
20, 14, 12, 12, 20, 20, 0, 0, 14, 12, 10, 10, 20, 16, 16, 14,
10, 10, 10, 20, 16, 10, 0, 12, 12, 0, 12, 16, 14, 16, 14, 0,
0, 12, 20, 0, 12, 14, 14, 0, 0, 20, 20, 20, 14, 14, 10, 10, 20,
16, 16, 0, 12, 10, 10, 10, 16, 16, 12, 20, 10, 12, 12, 16, 14,
0, 16, 20, 12, 14, 10, 10, 0, 0, 12, 12, 10, 10, 0, 0, 0, 14,
12, 12, 20, 20, 14, 14, 14, 12, 20, 20, 20, 16, 16, 14, 10, 10,
16, 16, 16), Success = c(100, 91, 87, 100, 0, 0, 91, 96, 89,
96, 82, 99, 99, 99, 92, 59, 45, 0, 100, 95, 100, 83, 8, 82, 63,
98, 74, 9, 0, 99, 78, 98, 53, 96, 52, 0, 62, 85, 98, 4, 89, 99,
99, 97, 82, 80, 91, 99, 96, 0, 0, 95, 80, 68, 74, 100, 97, 93,
87, 34, 32, 47, 99, 96, 86, 15, 93, 86, 0, 77, 89, 80, 98, 96,
31, 0, 61, 100, 0, 84, 88, 97, 0, 0, 99, 100, 100, 92, 88, 46,
51, 99, 97, 100, 0, 93, 61, 91, 57, 76, 95, 50, 98, 16, 87, 93,
87, 88, 11, 92, 98, 60, 96, 0, 64, 72, 0, 74, 77, 0, 62, 0, 0,
0, 84, 47, 69, 98, 100, 61, 90, 79, 11, 100, 98, 100, 98, 93,
91, 58, 58, 93, 78, 69), Failure = c(0, 9, 13, 0, 100, 100, 9,
4, 11, 4, 18, 1, 1, 1, 8, 41, 55, 100, 0, 5, 0, 17, 92, 18, 37,
2, 26, 91, 100, 1, 22, 2, 47, 4, 48, 100, 38, 15, 2, 96, 11,
1, 1, 3, 18, 20, 9, 1, 4, 100, 100, 5, 20, 32, 26, 0, 3, 7, 13,
66, 68, 53, 1, 4, 14, 85, 7, 14, 100, 23, 11, 20, 2, 4, 69, 100,
39, 0, 100, 16, 12, 3, 100, 100, 1, 0, 0, 8, 12, 54, 49, 1, 3,
0, 100, 7, 39, 9, 43, 24, 5, 50, 2, 84, 13, 7, 13, 12, 89, 8,
2, 40, 4, 100, 36, 28, 100, 26, 23, 100, 38, 100, 100, 100, 16,
53, 31, 2, 0, 39, 10, 21, 89, 0, 2, 0, 2, 7, 9, 42, 42, 7, 22,
31), Rep = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L)), row.names = 433:576, class = "data.frame")
#Model fitting
fit <- lme4::glmer(cbind(Success, Failure) ~ Dose + (1|factor(Rep)), family = binomial, data = df)
#Brier scores
DescTools::BrierScore(fit)
#Returns following error
#Error in mean(resp * (1 - pred)^2 + (1 - resp) * pred^2)
scoring::brierscore(fit)
#Returns following error
#argument "pred" is missing, with no default
Is it possible to calculate brier scores for a mixed effects model? If so, how?
解决方案
I suppose you can simply get rid of the cbind construction, get a long data format and use the predict/response interface of DescTools::BrierScore().
Something like this might work:
library(DescTools)
tt <- aggregate(df[, c("Success", "Failure")], by=as.list(df[, c("Dose", "Rep")]), FUN = sum)
d.frm <- Untable(reshape(tt, idvar=c("Dose","Rep"), varying = c("Success", "Failure"),
direction = "long",
timevar="Success", v.names="Freq", times=c(1, 0)))
fit <- lme4::glmer(Success ~ Dose + (1|Rep), family = binomial, data = d.frm)
BrierScore(pred=predict(fit, type="response"), resp = d.frm$Success)
这篇关于计算成功/失败总数的二项式GLMM的brier分数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！