0
votes

I have a simple Bayesian hierachical model (sir model) which should be easy to run. The problem is that after successfully loading the model and data, I get the following error when I try to compile the model, "array index is greater than array upper bound for muc." It seems like this should be an easy fix, but I've repeatedly checked the indexing, and the data, and cannot find a problem.

Any suggestions would be greatly appreciated!

#model
model{
for (i in 1:M){
rem[i,1]<-0
susc[i,1]<-susint[i]
muc[i,1]<-susc[i,1]
cpos[i,1]~dpois(muc[i,1])
}

for (i in 1:M){
for (j in 2: T){

rem[i,j]<-betaR*cpos[i,j]
susc[i,j]<-susc[i,j-1]-cpos[i,j-1]-rem[i,j-1]
cpos[i,j]~dpois(muc[i,j])
log(muc[i,j])<-bet0+log(susc[i,j]+0.001)+log(cpos[i,j-1]+0.001)+b1[i]
}

muct1[i]<-muc[i,1]
remt1[i]<-rem[i,1]
sust1[i]<-susc[i,1]
}
for (j in 1:T){

mucrich[j]<-muc[40,j]
mucchar[j]<-muc[10,j]
muchor[j]<-muc[26,j]
mucbea[j]<-muc[7,j]}

b1[1:14] ~ car.normal(adj[], weights[], num[], tau.b1)
for(k in 1:sumNumNeigh) {
  weights[k] <- 1 }

bet0~dflat()
tau.b1~dgamma(0.01,0.01)
#betaR~dgamma(0.01,1.0)
betaR<-0.001
}

#data

list( M=14,T=1,cpos=structure(.Data= c(
129, 843, 933, 1160, 1387, 2374, 506, 337, 1700, 779, 730, 1354, 3443, 1821),    
.Dim=c(14,1)),susint = c(
236200, 1958100, 1546800, 2351300, 1695900, 5004600, 995600, 753500, 3312400, 1513100, 1094300, 1595000, 3219200,2452800),
num = c(
1, 3, 2, 5, 1, 4, 4,2, 3, 6, 2, 3, 1, 1),  
adj = c(
2,          
4, 3, 1,    
4, 2,       
8, 7, 5, 3, 2, 
8,       
4, 7, 10, 9, 
11, 12, 8, 7,  
12, 10,     
11, 10, 7,  
6, 5, 4, 8, 10, 12, 
7, 5,   
7, 6, 4,    
13, 
15),

sumNumNeigh=38,     #sumNumNeigh=total of neighbors

)

inits

list(bet0=0.01,tau.b1=0.1,b1=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0)
list(bet0=0.03,tau.b1=0.2,b1=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0)
1

1 Answers

0
votes

Your data specifies M=14 and T=1 which implies that muc has 14 rows and 1 column. You then try to access the non-existent 40th row of muc with mucrich[j]<-muc[40,j] and get the error you quote.

I guess the mistake is really with the data - in particular I doubt you meant T=1 (or otherwise your model has a lot of unnecessary code in it!).