JAGS：使用单元格矩阵作为排除 NA 的数据输入
JAGS: Use cell matrix as data input for excluding NA
以下是我将使用的模型和示例数据。数据中有一些 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)
所以这里的解决方法比我的标准嵌套索引更容易一点，因为您总是在矩阵的右侧丢失数据（即，一旦数据为 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 ```
以下是我将使用的模型和示例数据。数据中有一些 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)
所以这里的解决方法比我的标准嵌套索引更容易一点，因为您总是在矩阵的右侧丢失数据（即，一旦数据为 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 ```