0

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

1 回答 1

0

你有几个选择

您可以即时构建模型字符串。[ to 的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。

我使用.xand.y作为三个fitX函数的参数名称，以匹配group_map.

于 2021-06-08T15:29:51.467 回答