I've been tasked with re-running someone else's model with updated data (I did the data update and formatting). I understand how the model works but didn't write it, and it's quite long and detailed. On top of that, I exclusively use JAGS and this is my first time wading into WinBUGS, so the interface and error messages are new to me. The model takes two days to compile (and it does compile), but then I get a Trap error that says Index Out of Range followed by a few dozen lines that are unintelligible to me. I've gone over the entire model, for-loops, and data structure and I can't see where there are any indexing issues. And, frustratingly, unlike JAGS, WinBUGS does not appear to tell you which index is out of range.
I'd greatly appreciate any assistance decoding the Trap error message, in case there's more information in there than I realize to specify the location of the index mismatch. I'm going to post the model below, along with the trap message, but I don't see the point in posting the actual data as no one is going to spend 2 days waiting for the model to compile in order to reproduce the error. Instead I'll just post the data structure so the dimensions of each data object are obvious.
A few considerations: I'm using the exact model that was previously used (literally copied and pasted), so I know that it ought to work. The previous data files are in rectangular format, whereas the ones I've created and exported are in S format. My understanding is that WinBUGS should be able to handle both, but am wondering if that might be part of the problem. I'm running the model directly from WinBUGS, loading model and data files manually. I tried doing it from R with r2winbugs previously but was getting memory errors and crashes.
Model:
model{
for(i in 1:n+nzeros){ # loop through the observed and unobserved areas
for(j in 1:J){ #loop through the J fisheries areas
for(t in 1:T){ # loop through the years
y[i,j,t] ~ dpois(mu.y[i,j,t])
mu.y[i,j,t] <- lambda[i,j,t]*x[i,t]
log(lambda[i,j,t]) <- mu.lambda[j]+ theta[z[i],j] + e.lambda[j,t]
}}}
for(i in 1:n+nzeros){
x[i,1] ~ dbern(gamma[i,1])
recruitable[i,1]<-1
for(t in 2:T){
x[i,t] ~ dbern(mu.x[i,t])
mu.x[i,t] <- survived[i,t] + gamma[i,t]*recruitable[i,t]
recruitable[i,t] <- recruitable[i,t-1]*(1-x[i,t-1])
survived[i,t] <- x[i,t-1]*phi[i,t-1]
}}
#priors
for(c in 1:C){
theta[c,1:J] ~ dmnorm(mu.th[], tau.th[,])
for(j in 1:J){
lcapture.clus[c,j]<-mu.lambda[j]+theta[c,j]
pcapture.clus[c,j]<-1-exp(-exp(lcapture.clus[c,j]))
}
}
for(j in 1:J){
mu.th[j] <-0
mu.lambda[j]~dnorm(0, 0.001)
pcapture.overall[j]<-1-exp(-exp(mu.lambda[j]))
for (t in 1:T){
e.lambda[j,t] ~ dnorm(0, tau.e)
lambda.area[j,t]<-mu.lambda[j]+e.lambda[j,t]
p.area[j,t]<-1-exp(-exp(lambda.area[j,t]))
#n.area[j,t]~dbin(p.area[j,t],N.area[j,t])
N.area[j,t]<-n.area[j,t]/p.area[j,t]
#N.area[j,t] ~ dunif(n.area[j,t],500)
#N.area.round[j,t]<-round(N.area[j,t])
}}
tau.e <- 1/(sd.e*sd.e)
sd.e ~ dunif(0,10)
tau.th[1:J, 1:J] ~ dwish(B[,], v)
sigma2[1:J, 1:J] <- inverse(tau.th[1:J,1:J])
v<-J
for(j1 in 1:J){
B[j1,j1]<-1
for(j2 in 1:j1-1){
B[j1,j2]<-0
}
for(j2 in j1+1:J){
B[j1,j2]<-0
}
for(j2 in 1:J){
pd[j1,j2] <- step(sigma2[j1,j2])
}}
for(i in 1:n+nzeros){
for(t in 1:T-1){
logit(phi[i,t]) <- logit(mu.phi) + e.phi[z.cut[i],t]
}
gamma[i,1]~dunif(0,1)
for(t in 2:T){
logit(gamma[i,t]) <-logit(mu.gamma) + e.gamma[z.cut[i],t]
}}
for(c in 1:C){
for (t in 1:T-1){
e.phi[c,t] ~ dnorm(0, tau.phi)
prob.phi[c,t]<-step(e.phi[c,t])
}
for(t in 2:T){
e.gamma[c,t] ~ dnorm(0, tau.gamma)
}}
mu.gamma ~ dunif(0,1)
tau.gamma <- 1/(sd.gamma*sd.gamma)
sd.gamma ~ dunif(0,10)
mu.phi ~ dunif(0,1)
tau.phi <- 1/(sd.phi*sd.phi)
sd.phi ~ dunif(0,10)
# Clustering
for(i in 1:n+nzeros){
z[i] ~ dcat(w[])
z.cut[i]<-cut(z[i])
for(c in 1:C){
post.prob[i,c] <- equals(z[i],c)}}
for(c in 1:C){
w[c] <- ww[c]/sum(ww[])
r[c] ~ dbeta(1,a)
NC[c] <- sum(post.prob[,c])
cl[c] <- step(NC[c]-1)}
ww[1] <- r[1]
for(c in 2:C){
ww[c] <- r[c]*(1-r[c-1])*ww[c-1]/r[c-1]}
CLUSTERS <- sum(cl[])
a<- 1 # controls prior smoothing
#derived parameters
for(i in 1:n+nzeros){
for(c in 1:C){
recruit[i,c,1]<-x[i,1]*equals(z[i],c)
for(t in 2:T){
recruit[i,c,t]<-(1-x[i,t-1])*(x[i,t])*equals(z[i],c)
death[i,c,t] <- x[i,t-1]*(1-x[i,t])*equals(z[i],c)
alive[i,c,t]<-x[i,t]*equals(z[i],c)
}}}
for(c in 1:C){
for(t in 2:T){
R[c,t]<-sum(recruit[1:n+nzeros,c,t])
D[c,t]<-sum(death[1:n+nzeros,c,t])
N[c,t]<-sum(alive[1:n+nzeros,c,t])
Rpc[c,t]<-R[c,t]/N[c,t]
recruitment.clus[c,t]<- exp(logit(mu.gamma) + e.gamma[c,t])/(1+exp(logit(mu.gamma) + e.gamma[c,t]))
}
for(t in 1:T-1){
survival.clus[c,t] <- exp(logit(mu.phi) + e.phi[c,t])/(1+exp(logit(mu.phi) + e.phi[c,t]))
}}
for(t in 2:T){
NN[t]<-sum(N[1:C, t])
DD[t]<-sum(D[1:C, t])
RR[t]<-sum(R[1:C, t])
RRpc[t]<-RR[t]/NN[t]
}
#partitioned p values
for (i in 1:n){
for(j in 1:J){
for (t in 1:T){
y.new[i,j,t]~dpois(mu.y[i,j,t])
# for predictive goodness of fit
p0[i,j,t] <- equals(y[i,j,t],0)
p0new[i,j,t]<-equals(y.new[i,j,t],0)
count[i,j,t]<-y[i,j,t]
countnew[i,j,t]<-y.new[i,j,t]
for (c in 1:C){
d.p0[i,j,t,c]<-p0[i,j,t]*equals(z[i],c)
d.p0new[i,j,t,c]<-p0new[i,j,t]*equals(z[i],c)
d.count[i,j,t,c]<-count[i,j,t]*equals(z[i],c)
d.countnew[i,j,t,c]<-countnew[i,j,t]*equals(z[i],c)
}}}}
for(j in 1:J){
for(t in 1:T){
for(c in 1:C){
fit1.p0[j,t,c] <- sum(d.p0[1:n,j,t,c])
fit.new1.p0[j,t,c] <- sum(d.p0new[1:n,j,t,c])
fit1.count[j,t,c] <- sum(d.count[1:n,j,t,c])
fit.new1.count[j,t,c] <- sum(d.countnew[1:n,j,t,c])
}}}
for(t in 1:T){
for(c in 1:C){
fit2.p0[t,c] <- sum(fit1.p0[1:J,t,c])
fit.new2.p0[t,c] <- sum(fit.new1.p0[1:J,t,c])
fit2.count[t,c] <- sum(fit1.count[1:J,t,c])
fit.new2.count[t,c] <- sum(fit.new1.count[1:J,t,c])
}}
for(c in 1:C){
DF.p0[c] <- sum(fit2.p0[,c])
DF.new.p0[c] <- sum(fit.new2.p0[,c])
pvalue.p0[c] <- step(DF.new.p0[c] - DF.p0[c])
DF.count[c] <- sum(fit2.count[,c])
DF.new.count[c] <- sum(fit.new2.count[,c])
pvalue.count[c] <- step(DF.new.count[c] - DF.count[c])
}
DFoverall.p0<-sum(DF.p0[1:C])
DFoverall.new.p0<-sum(DF.new.p0[1:C])
pvalueoverall.p0<-step(DFoverall.new.p0-DFoverall.p0)
DFoverall.count<-sum(DF.count[1:C])
DFoverall.new.count<-sum(DF.new.count[1:C])
pvalueoverall.count<-step(DFoverall.new.count-DFoverall.count)
}
Input data (I did all of the data formatting in R, so I compile this as a list and then export it with r2winbugs in S/R format, so presumably it should work in WinBUGS):
bugs.data <- list(y = y,
x = x.data,
J = dim(y)[2],
T = dim(y)[3],
n = dim(y)[1]-100,
nzeros = 100,
C = 10,
n.area = n.area)
Data Structure (and I confirmed these manually in the actual text file):
> dim(y)
[1] 727 7 41
> dim(x.data)
[1] 727 41
> dim(n.area)
[1] 7 41
> J
[1] 7
> T
[1] 41
> n
[1] 627
And finally the Trap error:
index out of range
HostWindows.AppendInt [000002B7H]
.d ARRAY 12 OF CHAR "080'069'784'"
.i INTEGER 12
.j INTEGER 1971081876
.len INTEGER 1971172856
.n INTEGER 1
.s ARRAY 256 OF CHAR " Allocated Memory: " ...
.useSeparators BOOLEAN TRUE
HostWindows.UpdateInfo [000048B8H]
.res INTEGER 2291852
.sstr ARRAY 256 OF SHORTCHAR "l˜"" ...
.str ARRAY 256 OF CHAR " Allocated Memory: " ...
HostWindows.Idle [00004A73H]
.focus BOOLEAN FALSE
.tick Controllers.TickMsg Fields
.w HostWindows.Window NIL
HostMenus.TimerTick [00003422H]
.lParam INTEGER 0
.ops Controllers.PollOpsMsg Fields
.wParam INTEGER 1
.wnd INTEGER 459558
Kernel.Try [00003A61H]
.a INTEGER 459558
.b INTEGER 1
.c INTEGER 0
.h PROCEDURE HostMenus.TimerTick
HostMenus.ApplWinHandler [00003841H]
.Proc PROCEDURE NIL
.hit BOOLEAN Undefined117
.lParam INTEGER 0
.message INTEGER 275
.res INTEGER 2292416
.s ARRAY 256 OF SHORTCHAR 9AX, 4X, 0CX ...
.w INTEGER 0
.wParam INTEGER 1
.wnd INTEGER 459558
<system> (pc=757CC4E6H, fp=0022FB3CH)
<system> (pc=757CC5E6H, fp=0022FBB4H)
<system> (pc=757CCC18H, fp=0022FC14H)
<system> (pc=757C2E40H, fp=0022FC24H)
HostMenus.Loop [00003BDEH]
.done BOOLEAN FALSE
.f SET {0..5}
.n INTEGER 10
.res INTEGER 0
.w HostWindows.Window NIL
Kernel.Start [00002B8CH]
.code PROCEDURE HostMenus.Loop
Thanks a ton in advance for your help! -Josh
[winbugs] index out of range? Reading your comment on the possibility of porting to JAGS, I also wonder if using OpenBUGS might be explored. – Heeler