2

以下是我将使用的模型和示例数据。我需要设置先验以生成数字的数据中有一些 NA，但这种方式可能会导致一些错误。我想知道我是否可以让 JAGS 跳过 NA，就像拥有一个具有不同行和列的矩阵一样。

NA 在ex_expectancy和中ex_shock。

# data
# 3 subjects * 14 trials
ex_expectancy <- structure(list(`1` = c(9L, 5L, 1L), `2` = c(5L, 6L, 1L), `3` = c(2L, 7L, 4L), `4` = c(3L, 6L, 2L), `5` = c(9L, 6L, 4L), `6` = c(9L, 7L, 1L), `7` = c(3L, 5L, 5L), `8` = c(8L, 5L, 1L), `9` = c(10L, 5L, NA), `10` = c(9L, NA, NA), `11` = c(2L, NA, NA), `12` = c(3L,NA, NA), `13` = c(3L, NA, NA), `14` = c(4L, NA, NA)), row.names = c(NA,-3L), class = c("data.table", "data.frame"))
ex_shock <- structure(list(`1` = c(0L, 1L, 1L), `2` = c(0L, 1L, 0L), `3` = c(1L, 0L, 1L), `4` = c(1L, 0L, 1L), `5` = c(0L, 1L, 1L), `6` = c(0L, 0L, 0L), `7` = c(1L, 0L, 1L), `8` = c(1L, 1L, 1L), `9` = c(0L,1L, NA), `10` = c(1L, NA, NA), `11` = c(1L, NA, NA), `12` = c(0L, NA, NA), `13` = c(1L, NA, NA), `14` = c(0L, NA, NA)), row.names = c(NA,-3L), class = c("data.table", "data.frame"))

v <- matrix(NA, nrow=3,ncol=14)
v[,1] <- 0 # first v is 0

dlist <- list(
  Nsubjects = 3,
  Ntrials = 14,
  expectancy = ex_expectancy, 
  shock = ex_shock,
  v=v
)

myinits <-  list(list(
  alpha = runif (3,0,1))) # 3 subjects
parameters <- c('alpha','v','predk','scale','c','tau')

# model
RW <- function(){
  for (i in 1:Nsubjects)
  {
    for (j in 2:Ntrials) # for each trial
    {
      expectancy [i,j] ~ dnorm (scale [i] * v[i,j] + c[i],tau[i,j])
      # posteiror predictive
      predk [i,j] ~ dnorm (scale [i] * v[i,j] + c[i],tau[i,j])
      
      
      pe [i,j-1] <- shock [i,j-1] - v [i,j-1]
      v [i,j] <- v [i,j-1] + alpha [i]  * pe [i,j-1]
    }
  }
  # priors
  for (i in 1: Nsubjects){
    alpha [i] ~ dunif (0,1)
    scale [i] ~ dunif (0,10)
    c[i] ~ dunif (0,5)
    
    for (j in 1:Ntrials){
      sigma[i,j] ~ dunif (0,5)
      tau [i,j] <- 1/pow(sigma [i,j],2)
    }}
}


samples <- jags(dlist, inits=myinits, parameters,
                model.file = RW,
                n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T)
4

1 回答 1

1

因此，这里的解决方法比我的标准嵌套索引要容易一些，因为您总是在矩阵的右侧丢失数据（即，一旦数据为 NA，则该列的其余部分为 NA）。因此，无需在循环中应用嵌套索引，您只需将其应用于第二个for循环（我在runjags这里使用它，因为这是我最熟悉的）。

# data
ex_expectancy <- structure(list(`1` = c(9L, 5L, 1L), `2` = c(5L, 6L, 1L), `3` = c(2L, 7L, 4L), `4` = c(3L, 6L, 2L), `5` = c(9L, 6L, 4L), `6` = c(9L, 7L, 1L), `7` = c(3L, 5L, 5L), `8` = c(8L, 5L, 1L), `9` = c(10L, 5L, NA), `10` = c(9L, NA, NA), `11` = c(2L, NA, NA), `12` = c(3L,NA, NA), `13` = c(3L, NA, NA), `14` = c(4L, NA, NA)), row.names = c(NA,-3L), class = c("data.table", "data.frame"))
ex_shock <- structure(list(`1` = c(0L, 1L, 1L), `2` = c(0L, 1L, 0L), `3` = c(1L, 0L, 1L), `4` = c(1L, 0L, 1L), `5` = c(0L, 1L, 1L), `6` = c(0L, 0L, 0L), `7` = c(1L, 0L, 1L), `8` = c(1L, 1L, 1L), `9` = c(0L,1L, NA), `10` = c(1L, NA, NA), `11` = c(1L, NA, NA), `12` = c(0L, NA, NA), `13` = c(1L, NA, NA), `14` = c(0L, NA, NA)), row.names = c(NA,-3L), class = c("data.table", "data.frame"))


v <- matrix(NA, nrow=3,ncol=14)
v[,1] <- 0

dlist <- list(
  NSubjects = 3,
  Ntrials = 14 - rowSums(is.na(ex_shock)),
  maxTrials = 14,
  expectancy = as.matrix(ex_expectancy), 
  shock = as.matrix(ex_shock),
  v = v
)

myinits <-  list(list(
  alpha = runif (3,0,1)))
parameters <- c('alpha','v','predk','scale','c','tau')


{sink("model.txt")
cat("
model{
    for (i in 1:NSubjects){
      for (j in 2:Ntrials[i]){
        expectancy[i,j] ~ dnorm (scale[i] * v[i,j] + c[i],tau[i,j])
        # posteiror predictive
        predk[i,j] ~ dnorm (scale [i] * v[i,j] + c[i],tau[i,j])
        pe[i,j-1] <- shock[i,j-1] - v[i,j-1]
        v[i,j] <- v[i,j-1] + alpha[i]  * pe[i,j-1]
      }
    }
    # priors
    for (i in 1: NSubjects){
      alpha[i] ~ dunif (0,1)
      scale[i] ~ dunif (0,10)
      c[i] ~ dunif (0,5)
      
      for (j in 1:maxTrials){
        sigma[i,j] ~ dunif (0,5)
        tau[i,j] <- 1/pow(sigma [i,j],2)
      }}
  }"
  ,fill = TRUE)
}
sink()



library(runjags)


samples <- run.jags("model.txt", monitor = parameters, data = dlist,
                    n.chains = 2,sample = 10000, burnin = 5000,
                    thin = 1)

基本上Ntrials变成一个长度向量NSubjects。通过应用微小的更改，模型将编译并运行。但是，这并不能解决模型的任何潜在拟合问题。由于我不确定您实际适合什么，因此我不知道模型是否按指定正确。查看 mcmc 的输出，它看起来好像仍在发生一些奇怪的事情（predk和tau的某些部分NA）。

library(coda)
 my_mcmc <- as.matrix(as.mcmc.list(samples))

round(my_mcmc[1,],2)
  alpha[1]    alpha[2]    alpha[3]      v[1,1]      v[2,1]      v[3,1] 
       0.23        0.13        0.75        0.48        0.32        0.12 
     v[1,2]      v[2,2]      v[3,2]      v[1,3]      v[2,3]      v[3,3] 
       0.23        0.09        0.05        0.08        0.29        7.60 
     v[1,4]      v[2,4]      v[3,4]      v[1,5]      v[2,5]      v[3,5] 
       0.05        0.11        0.10        0.15        0.62        0.00 
     v[1,6]      v[2,6]      v[3,6]      v[1,7]      v[2,7]      v[3,7] 
       0.00        0.00        0.00        0.13        0.75        0.00 
     v[1,8]      v[2,8]      v[3,8]      v[1,9]      v[2,9]     v[1,10] 
       0.24        0.19        0.23        0.21        0.80        0.41 
    v[1,11]     v[1,12]     v[1,13]     v[1,14]  predk[1,2]  predk[2,2] 
       0.18        0.95        0.32        0.29          NA          NA 
 predk[3,2]  predk[1,3]  predk[2,3]  predk[3,3]  predk[1,4]  predk[2,4] 
         NA        5.30       -0.48        1.83        2.01        7.30 
 predk[3,4]  predk[1,5]  predk[2,5]  predk[3,5]  predk[1,6]  predk[2,6] 
       0.57        2.77        6.49        2.37        2.82        5.76 
 predk[3,6]  predk[1,7]  predk[2,7]  predk[3,7]  predk[1,8]  predk[2,8] 
       6.23       -4.78        7.10        5.12        3.10        0.95 
 predk[3,8]  predk[1,9]  predk[2,9] predk[1,10] predk[1,11] predk[1,12] 
      -0.34       -0.31       10.04        4.13        2.60        9.53 
predk[1,13] predk[1,14]    scale[1]    scale[2]    scale[3]        c[1] 
         NA       10.83          NA          NA        1.48        2.46 
       c[2]        c[3]    tau[1,1]    tau[2,1]    tau[3,1]    tau[1,2] 
       4.75        1.36          NA          NA       10.40          NA 
   tau[2,2]    tau[3,2]    tau[1,3]    tau[2,3]    tau[3,3]    tau[1,4] 
         NA        2.61          NA          NA       -2.68          NA 
   tau[2,4]    tau[3,4]    tau[1,5]    tau[2,5]    tau[3,5]    tau[1,6] 
         NA        1.35        8.25        1.19        0.14        0.11 
   tau[2,6]    tau[3,6]    tau[1,7]    tau[2,7]    tau[3,7]    tau[1,8] 
       0.08        0.10        0.09        0.22        0.76        4.85 
   tau[2,8]    tau[3,8]    tau[1,9]    tau[2,9]    tau[3,9]   tau[1,10] 
       0.70       16.36        1.66        1.59        0.05        3.98 
  tau[2,10]   tau[3,10]   tau[1,11]   tau[2,11]   tau[3,11]   tau[1,12] 
       0.07        0.05       53.40        0.08        9.30        0.08 
  tau[2,12]   tau[3,12]   tau[1,13]   tau[2,13]   tau[3,13]   tau[1,14] 
       0.18        0.10        0.12        0.87        0.08        0.09 
  tau[2,14]   tau[3,14] 
       5.40        0.04 ```
于 2021-03-15T15:19:55.493 回答