更新每次模拟中的观察次数使用库：runjags、parallel

编程问答 2022-04-25

问题描述

我正在按模拟 ID 拆分数据集，并同时将 runjags 函数应用于每个子集。

现在，每个模拟都包含 1000 个观察值。我知道有时观察的数量会有所不同，因为我会删除符合某些标准的行。我不知道会丢弃多少观察值，但我可以使用 groupobs <- fulldata %>% count(SimulID,sort=TRUE).

来计算

有没有办法在每次模拟运行期间更改 N=1000。这意味着必须在每次运行模拟时重写 tempModel.txt 文件。

谢谢。

#Subset data by SimulID
subsetdata <- split(fulldata,as.factor(fulldata$SimulID))
#Count obs within each group
groupobs <- fulldata %>% count(SimulID,sort=TRUE)

modelString <- "
  model{
#Model specification
   for (i in 1:1000) {
      y[i]~dnorm(muy[i],Inv_sig2_e)
      muy[i]<-b0+b1*x1[i]+b2*x2[i]
   }
#priors
   b0~dnorm(0,1.0E-6)
   b1~dnorm(0,1.0E-6)
   b2~dnorm(0,1.0E-6)
   Inv_sig2_e~dgamma(1.0E-3,1.0E-3)
#parameter transformation
   Sig2_e<-1/Inv_sig2_e
  }
"

writeLines(modelString,"tempModel.txt")

output_models <- lapply(subsetdata,function(x){
  model_data = x
  initsList1 <- list(b0=1,b1=1,b2=1,Inv_sig2_e=1)
  initsList2 <- list(b0=1,b1=2,b2=3,Inv_sig2_e=1)
  initsList3 <- list(b0=2,b1=3,b2=4,Inv_sig2_e=1)





 runJagsOut <- run.jags(method = "parallel",model = "tempModel.txt",# NOTE: theta and omega are vectors:
                         monitor = c( "b0","b1","b2","Sig2_e"),data = model_data,inits = list(initsList1,initsList2,initsList3),# NOTE: Let JAGS initialize.
                         n.chains = 3,# NOTE: Not only 1 chain.
                         adapt = 500,burnin = 2500,sample = 2500,thin = 1,summarise = FALSE,plots = FALSE)
})

解决方法

您有多种选择

您可以即时构建模型字符串。 [model 的 run.jags 参数可以包含字符串而不是文件名，因此无需写入文件然后再次读取。]

您可以向 data 列表（代码中的 x）添加一个包含观察次数的元素，

x[["groupobs"]] <- fulldata %>% count(SimulID,sort=TRUE)

并在您的 model_string 中引用：

for (i in 1:groupobs)

您可以实时计算观察次数：

for (i in 1:length(y))

在您的 model_string 中。

编辑 为了回应 OP 的评论，这里是我上面三个建议中每一个的实现。 OP 的代码不可重现，因为他们没有提供他们的数据，所以我将重新分析 O'Quigley 等人在他们的 1990 CRM paper 中使用的一个例子。为了重现 OP 的分组分析，我将复制数据并简单地分析两次。

输入数据：

dput(observedData)
structure(list(Cohort = c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L,14L,15L,16L,17L,18L,19L,20L,21L,22L,23L,24L,25L,1L,25L),SubjectID = c(1L,25L
),Dose = c(3,4,3,2,1,1),Toxicity = c(0,Trial = c(1,2)),row.names = c(NA,-50L),class = c("tbl_df","tbl","data.frame"))

我发现 tidyverse 的 group_map 函数提供的代码比 lapply 更紧凑、更易于理解，所以我会使用它。

library(tidyverse)
library(runjags)

选项 1：将观察计数粘贴到模型字符串中。

modelString <- 
"model { 
  #Prior 
  a ~ dexp(1) 
  #Likelihood 
  for (i in 1:n) { 
    Toxicity[i] ~ dbern(((tanh(XHat[i]) + 1)/2)**a) 
  } 
} 
#monitor# a"

fit1 <- function(.x,.y) {
  modelString <- paste0(
    "model { 
      #Prior 
      a ~ dexp(1) 
      #Likelihood 
      for (i in 1:",.x %>% nrow(),") { 
        Toxicity[i] ~ dbern(((tanh(XHat[i]) + 1)/2)**a) 
      } 
    } 
    #monitor# a")
  d <- list(XHat=.x$Dose,Toxicity=.x$Toxicity)
  run.jags(modelString,data=d)
}

observedData %>% group_by(Trial) %>% group_map(fit1)

选项 2：将观察计数作为 data

的元素传递 
modelString <- 
  "model { 
    #Prior 
    a ~ dexp(1) 
    #Likelihood 
    for (i in 1:n) { 
      Toxicity[i] ~ dbern(((tanh(XHat[i]) + 1)/2)**a) 
    } 
   } 
   #monitor# a"

fit2 <- function(.x,.y) {
  d <- list(XHat=.x$Dose,Toxicity=.x$Toxicity,n=.x %>% nrow())
  run.jags(modelString,data=d)
}

observedData %>% group_by(Trial) %>% group_map(fit2)

选项 3：让 JAGS 计算观察次数

modelString <- 
  "model { 
    #Prior 
    a ~ dexp(1) 
    #Likelihood 
    for (i in 1:length(Toxicity)) { 
      Toxicity[i] ~ dbern(((tanh(XHat[i]) + 1)/2)**a) 
    } 
  } 
  #monitor# a"

fit3 <- function(.x,data=d)
}

observedData %>% group_by(Trial) %>% group_map(fit3)

我的个人偏好是选项 2。

我使用 .x 和 .y 作为三个 fitX 函数的参数名称，以匹配 group_map 的在线文档中使用的约定。

mcmcparallel-processingrrrunjagssimulation