r语言 - 分层狄利克雷回归(jags)中的随机截距

Question

我有以下数据结构：

• y：3列，是多年来观察到的死亡人数比例。
• x1：国内生产总值 - 与每年相关的连续变量
• x2：与死亡相关的年龄

这里是型号规格：

型

这里模拟数据：

library(tidyverse)
Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)
# X2 AGES
n.ages <- length(Ages)
# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)
# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)

Year.tot <- rep(Year,length(unique(Ages)))
DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),
C1=C_Porp[,1],
C2=C_Porp[,2],
C3=C_Porp[,3]) 
# X1 GDP OVER YEARS
GDP <- rnorm(n.year,50,15)
GDP <- as.data.frame(cbind(GDP,Year))
# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)

锯齿模型

library(runjags)
dirlichet.model = 


"model {
#setup priors for each species
for(j in 1:N.c){
m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
}

#implement dirlichet
for(i in 1:N){

y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 

for(j in 1:N.c){

a0[i,j] ~ dnorm(mu[i,j],0.001) 

log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])

}

# Effect of site:
for (s in 1:S){

alpha[s] ~ dnorm(0,0.001)

}}

}
"
jags.data <- list(y = D[,c(3,4,5)],mat = D[,5], Ages=D[,1],N = nrow(D[,c(3,4,5)]), 
N.c = ncol(D[,c(3,4,5)]),S=(length(unique(D[,1]))))
jags.out <- run.jags(dirlichet.model,
data=jags.data,
adapt = 500,
burnin = 5000,
sample = 10000,
n.chains=5,
monitor=c('m0','m1',"mu"))
out <- summary(jags.out)

我收到此错误：

Errore: The following error occured when compiling and adapting the model using rjags:
Error in rjags::jags.model(model, data = dataenv, n.chains = length(runjags.object$end.state),  : 
RUNTIME ERROR:
Compilation error on line 12.
Index out of range taking subset of  alpha

对型号规格有什么建议吗？

Answer 1

还有其他几个小错误。不过，这段代码现在应该可以工作了。 问题是：

1. Age被定义为从 0 到 50 x 5 的序列，但随后您使用Age来索引alpha。 我认为您真正想要的是每个值的不同 alphaAge. 我通过在数据中执行以下操作来解决此问题：Ages=match(D[,1], unique(D[,1])),
2. alpha正在被重新定义，因为N上的循环同时包围了N.c和S上的循环。 通过在循环开始之前关闭N循环S，可以解决问题。
3. GDP在数据中被定义为mat，所以我用GDP = D[,5]替换了mat = D[,5]。

完成这些更改后，模型将运行。

library(tidyverse)
Year <-2000:2010 
n.year <- length(Year)
Ages <- rep(seq(0,50,5),n.year)
# X2 AGES
n.ages <- length(Ages)
# Y Causes of Death 
Cause1 <- rpois(n.ages,60)
Cause2 <- rpois(n.ages,20)
Cause3 <- rpois(n.ages,90)
# Get proportion
C <- cbind(Cause1,Cause2,Cause3)
C_Porp <- C/rowSums(C)

Year.tot <- rep(Year,length(unique(Ages)))
DATA <- data.frame(Year=Year.tot,Ages=sort(Ages),
C1=C_Porp[,1],
C2=C_Porp[,2],
C3=C_Porp[,3]) 
# X1 GDP OVER YEARS
GDP <- rnorm(n.year,50,15)
GDP <- as.data.frame(cbind(GDP,Year))
# TOT
D <- as.matrix(merge(DATA,GDP))[,-1]
head(D)

library(runjags)
dirlichet.model = 


"model {
#setup priors for each species
for(j in 1:N.c){
m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
}

#implement dirlichet
for(i in 1:N){

y[i,1:N.c] ~ ddirch(a0[i,1:N.c]) 

for(j in 1:N.c){

a0[i,j] ~ dnorm(mu[i,j],0.001) 

log(mu[i,j]) <- (m0[j]+alpha[Ages[i]] + m1[j] * GDP[i])

}
}
# Effect of site:
for (s in 1:S){

alpha[s] ~ dnorm(0,0.001)

}

}
"
jags.data <- list(y = D[,c(2,3,4)],
GDP = D[,5], 
Ages=match(D[,1], unique(D[,1])),
N = nrow(D[,c(3,4,5)]), 
N.c = ncol(D[,c(3,4,5)]),
S=(length(unique(D[,1]))))
jags.out <- run.jags(dirlichet.model,
data=jags.data,
adapt = 500,
burnin = 5000,
sample = 10000,
n.chains=5,
monitor=c('m0','m1',"mu"))
out <- summary(jags.out)