如何将大量的Brm保存为不同的R对象或Rds / Rda文件？

编程问答 2022-06-05

问题描述

我正在尝试使用brms::brm()替换模型中一个参数的先验来拟合贝叶斯模型 purrr::map2()一个接一个（我对该参数有63个优先级）。由于我上一个问题的answer，我可以list2env()“理论上”将每个拟合模型另存为全局环境（即我的工作空间）中的不同对象。通过list2env()，将在所有63次迭代完成后保存对象（brmsfit）。但是，当我运行整个代码时，总会收到一条错误消息，说Error in scan(con,nlines = 1,sep = ",",quiet = TRUE) : Could not allocate memory (0 Mb) in C function 'R_AllocStringBuffer'，并且brmsfit对象没有存储在全局环境中，尽管模型拟合似乎完成了63次，最多我以前的存在。

因此，我想将每个brmsfit对象保存为R对象，并将其保存为.Rds或.Rda文件，在迭代完成后立即保存避免这种内存问题。但是，我应该怎么做才能实现呢？

MWE

请注意，以下命令只是一个“示意图”示例，其中包含我要尝试的公开数据。它可以工作，而不会遇到我上面提到的问题，因为在此示例中，brm()中的数据和模型更简单，并且先验数量比我要解决的要少得多。 

library(lme4)
library(tidyverse)
library(magrittr)
library(brms)
library(cmdstanr)
library(rstan)

## Parallelize the chains using all the cores:
options(mc.cores = parallel::detectCores())

prior_test <- data.frame(
  sd = c(0.01,0.01,0.01),mean = c(50,-50,0)
) %>% 
  mutate(
    id = row_number()
  )

list2env(
  purrr::map2(
    prior_test$sd,prior_test$mean,function(psd,pm){
          gc()
          gc()
          cbfm1 <- brm(
            Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),data = sleepstudy,family = "normal",prior =
              c(
                prior(normal(0,1),class = b,coef = Intercept),set_prior(
                  paste0(
                    "normal(",pm,psd,")"
                  ),class = "b",coef = "Days"
                ),prior(normal(0,class = sd),prior(lkj(2),class = cor)    
              ),save_pars = save_pars(all = TRUE),backend = "cmdstanr"
          )
        }
  ) %>%
    setNames(paste0('model',prior_test$id)),.GlobalEnv
)

带有追溯信息的错误消息

Error in scan(con,quiet = TRUE) : 
  Could not allocate memory (0 Mb) in C function 'R_AllocStringBuffer' 
15.
scan(con,quiet = TRUE) 
14.
rstan::read_stan_csv(out$output_files()) 
13.
.fit_model(model,...) 
12.
.fun(model = .x1,sdata = .x2,algorithm = .x3,backend = .x4,iter = .x5,warmup = .x6,thin = .x7,chains = .x8,cores = .x9,threads = .x10,inits = .x11,exclude = .x12,control = .x13,future = .x14,seed = .x15,silent = .x16) 
11.
eval(expr,envir,...) 
10.
eval(expr,...) 
9.
eval2(call,envir = args,enclos = parent.frame()) 
8.
do_call(fit_model,fit_args) 
7.
brm(voice ~ 0 + Intercept + agent + patient + REGION1 + REGION2 + 
    agent:patient + agent:REGION1 + agent:REGION2 + patient:REGION1 + 
    patient:REGION2 + agent:patient:REGION1 + agent:patient:REGION2 + 
    (0 + Intercept + agent + patient + agent:patient | subj) +  ... 
6.
.f(.x[[i]],.y[[i]],...) 
5.
map2(R1data$prior_sd,R1data$prior_mean,pm) {
    gc()
    gc()
    brm(voice ~ 0 + Intercept + agent + patient + REGION1 + REGION2 +  ... 
4.
eval(lhs,parent,parent) 
3.
eval(lhs,parent) 
2.
map2(R1data$prior_sd,pm) {
    gc()
    gc()
    brm(voice ~ 0 + Intercept + agent + patient + REGION1 + REGION2 +  ... 
1.
list2env(map2(R1data$prior_sd,pm) {
    gc()
    gc() ...

解决方法

我认为您可以使用类似以下的内容（这是一个更简单的版本，它避免了数据的来源，变量名，模型规格等）


dir.create("models")

# This function generates a file path,runs the `brm` function and saves the object
# returned by `brm()` as a .rds file
my_fun = function(prior_sd,prior_mean) {
  file_name = paste0("models/model_sd_",prior_sd,"_mean_",prior_mean,".rds")
  brms_object = brm(...) # pass formula,prior,data,options,etc.
  saveRDS(brms_object,file_name) # you can control compression too if you want
}

sds = c(0.01,0.01,0.01)
means = c(50,-50,0)
purrr::walk2(sds,means,my_fun)

您必须确保要使用的数据框存在于全局环境中

,

我对the accepted answer中的函数做了一些修改，以便（1）我们可以保留brmsfit 使用assign(...,envir = .GlobalEnv)和（2）在全局环境中创建对象，我们可以使用-来避免对象名称，当我们使用负值作为先验均值时，R解释为减法运算符。

data <- data.frame(
  prior_sd = prior_sd <- c(0.01,0.01),prior_mean = prior_mean <- c(50,0),id = case_when(
    prior_mean == 0 ~ paste0("model_mean_","_sd_",prior_sd),prior_mean < 0 ~ paste0("model_mean_m",abs(prior_mean),prior_mean > 0 ~ paste0("model_mean_p",prior_sd)
    )
)

computation = function(prior_mean,prior_sd) {
  file_name = case_when(
    prior_mean == 0 ~ paste0("path-to-your-directory/model_mean_",".rds"),prior_mean < 0 ~ paste0("path-to-your-directory/model_mean_m",prior_mean > 0 ~ paste0("path-to-your-directory/model_mean_p",".rds")
  )
  # pass formula,etc.
  brms_object = brm(
    Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),data = sleepstudy,family = "normal",prior =
      c(
        prior(normal(0,1),class = b,coef = Intercept),set_prior(
          paste0(
            "normal(",",")"
          ),class = "b",coef = "Days"
        ),prior(normal(0,class = sd),prior(lkj(2),class = cor)    
      ),save_pars = save_pars(all = TRUE),backend = "cmdstanr"
  )
  # save the computation results with object name
  assign(
    case_when(
      prior_mean == 0 ~ paste0("model_mean_",prior_mean <  0 ~ paste0("model_mean_m",prior_mean >  0 ~ paste0("model_mean_p",prior_sd)
    ),brms_object,envir = .GlobalEnv)
  # you can control compression too if you want
  saveRDS(brms_object,file_name) 
}

purrr::walk2(data$prior_mean,data$prior_sd,computation)
out-of-memorypurrrrrrstan