Kfold在brms中的应用

Question

我正在尝试使用kfold作为评估模型运行使用brms的手段，我觉得我错过了什么。作为一个可重复的例子，我的数据被构造成一个二进制响应(0，1)，这取决于个体的长度。下面是一些生成和绘制数据的代码，类似于我正在使用的代码：

library(brms)
library(tidyverse)
library(loo)

length <- seq(0, 100, by = 1)
n_fish_per_length <- 10

a0 <- -48
a1 <- 2
a2 <- -0.02

prob <- plogis(a0 + a1 * length + a2 * length^2)

plot(length, prob , type = 'l')

sim_data <-
  expand_grid(fish_id = seq_len(n_fish_per_length),
              length = length) %>%
  mutate(prob_use =  plogis(a0 + a1 * length + a2 * length^2)) %>%
  mutate(is_carp = rbinom(n = n(), size = 1, prob= prob_use))

ggplot(sim_data, aes(x = length, y = is_carp)) +
  geom_jitter(width = 0, height = 0.05) +
  geom_smooth(method = "glm", formula = y ~ x + I(x^2),
              method.args = list(family = binomial(link = "logit")))

然后我使用brms运行我的模型。

Bayes_Model_Binary <- brm(formula = is_carp ~ length + I(length^2),  
                          data=sim_data, 
                          family = bernoulli(link = "logit"),
                          warmup = 2500, 
                          iter = 5000, 
                          chains = 4, 
                          inits= "0", 
                          cores=4,
                          seed = 123)

summary(Bayes_Model_Binary)

我想用kfold来评价这个模型。我可以用这样的东西：

kfold(Bayes_Model_Binary, K = 10, chains = 1, save_fits = T)

但是我的数据中的反应是高度不平衡的(~18% = 1，~82% = 0)，而我的读数表明我需要使用分层的k折叠cv来解释这一点。如果我用：

sim_data$fold <- kfold_split_stratified(K = 10, x = sim_data$is_carp)

数据是按我预期的方式分开的，但我不知道从这里开始进行简历过程的最佳方式是什么。我看到了这个帖子https://mc-stan.org/loo/articles/loo2-elpd.html，但我不知道如何修改它来处理brmsfit对象。或者，我似乎应该能够使用：

kfold(Bayes_Model_Binary, K = 10, folds = 'stratified', group = sim_data$is_carp)

但这会造成一个错误。可能是因为is_carp是模型中的响应而不是预测器。在这种情况下，我的小组将是什么？我是不是漏掉了/误解了什么？我假设这里有一个非常简单的解决方案，我忽略了它，但欣赏任何想法。

Answer 1

经过一些额外的挖掘和学习如何访问分析中的每个折叠的信息之后，我能够确定数据的结构(响应中的0和1s的比例)是使用k折叠()函数中的默认设置来维护的。为此，我使用了以下代码。

首先，将kfold分析保存为一个对象。

kfold1 <- kfold(Bayes_Model_Binary, K = 10, save_fits = T)

kfold1 1$fits是模型拟合结果的列表，以及测试数据集(省略)中用于每个折叠的观测结果的列表。

根据这些信息，我创建了一个循环，在每个训练数据集中输出观察值的比例，其中is_carp =1(也可以对每个测试数据集这样做)使用以下代码。

for(i in 1:10){
    print(length(which(sim_data$is_carp[-kfold1$fits[i, ]$omitted] == 1)) / 
           nrow(sim_data[-kfold1$fits[i, ]$omitted, ]))
}

[1] 0.1859186
[1] 0.1925193
[1] 0.1991199
[1] 0.1914191
[1] 0.1881188
[1] 0.1848185
[1] 0.1936194
[1] 0.1980198
[1] 0.190319
[1] 0.1870187

然后很容易将这些比例与原始数据集中is_carp =1的观测值的比例进行比较。

length(which(sim_data$is_carp == 1)) / nrow(sim_data)

[1] 0.1910891