0

我正在尝试使用 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 流程的最佳方法是什么。我看到这篇文章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 是模型中的响应而不是预测变量。在这种情况下，我的小组会是什么？我在这里遗漏/误解了某些东西吗？我假设这里有一个非常简单的解决方案，我可以忽略但感谢任何想法。

4

1 回答 1

1

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

首先，将 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
于 2022-01-20T15:41:04.997 回答