每次模拟中更新观察的数量(使用库: runjags，并行)

Question

我通过模拟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)
})

Answer 1

你有几个选择

您可以动态构造模型字符串。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论文中使用的一个示例。为了重现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在线文档中使用的约定相匹配。