brms 中的 Kfold CV

Question

我正在尝试使用 kfold CV 作为使用 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 CV 来评估模型。我可以使用这样的东西：

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

但我的数据中的响应高度不平衡（~18% = 1，~82% = 0），我的阅读表明我需要使用分层 kfold cv 来解释这一点。如果我使用：

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

数据按照我预期的方式拆分，但我不确定从这里开始进行 CV 流程的最佳方式是什么。我看到了这个 post 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 是响应而不是模型中的预测变量。在这种情况下，我的团队会是什么？我 missing/misinterpreting 在这里吗？我假设这里有一个非常简单的解决方案，我忽略了但感谢任何想法。

Answer 1

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

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

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

kfold1$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

brms 中的 Kfold CV

Kfold CV in brms

r

cross-validation

k-fold

brms