问题描述
以下是我将使用的模型和示例数据。数据中有一些 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,`3` = c(2L,7L,4L),`4` = c(3L,2L),`5` = c(9L,`6` = c(9L,`7` = c(3L,5L),`8` = c(8L,`9` = c(10L,NA),`10` = c(9L,NA,`11` = c(2L,`12` = c(3L,`13` = c(3L,`14` = c(4L,NA)),row.names = c(NA,-3L),class = c("data.table","data.frame"))
ex_shock <- structure(list(`1` = c(0L,1L,`2` = c(0L,0L),`3` = c(1L,0L,`4` = c(1L,`5` = c(0L,`6` = c(0L,`7` = c(1L,`8` = c(1L,`9` = c(0L,`10` = c(1L,`11` = c(1L,`12` = c(0L,`13` = c(1L,`14` = c(0L,"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,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])
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,`3` = c(2L,7L,4L),`4` = c(3L,2L),`5` = c(9L,`6` = c(9L,`7` = c(3L,5L),`8` = c(8L,`9` = c(10L,NA),`10` = c(9L,NA,`11` = c(2L,`12` = c(3L,`13` = c(3L,`14` = c(4L,NA)),row.names = c(NA,-3L),class = c("data.table","data.frame"))
ex_shock <- structure(list(`1` = c(0L,1L,`2` = c(0L,0L),`3` = c(1L,0L,`4` = c(1L,`5` = c(0L,`6` = c(0L,`7` = c(1L,`8` = c(1L,`9` = c(0L,`10` = c(1L,`11` = c(1L,`12` = c(0L,`13` = c(1L,`14` = c(0L,"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,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])
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 ```