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

Question

我正在按模拟 ID 拆分我的数据集，并同时对每个子集应用运行jags 函数。

现在，每个模拟包含 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)
})

Answer 1

您有多种选择

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

您可以向 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, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
25L), SubjectID = 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, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L
), Dose = c(3, 4, 4, 3, 3, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 1, 1, 3, 4, 4, 3, 3, 2, 1, 1, 1, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1), Toxicity = c(0, 0, 1, 0, 
1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 
0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 
1, 0, 1, 1), Trial = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 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, .y) {
  d <- list(XHat=.x$Dose, Toxicity=.x$Toxicity)
  run.jags(modelString, data=d)
}

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

我个人倾向于选项 2。

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

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

Update number of observations in each simulation (using libraries: runjags, parallel)

simulation

parallel-processing

r

mcmc

runjags