[WinBUGS] Hierarchical Analysis in WinBUGS14 "Undefined Real Result"

[rneuenhoff]

Hello, I'm new here so I hope that I am posting in the correct forum. I am building a stock-recruitment hierarchical model for Atlantic cod in WinBUGS14. In this case, I'm estimating several hyperparameters across several different and exchangeable stocks.

BEFORE I can do that, I must standardize across stocks because recruit age is different for each stock according to the following:

Rp[i] = Rs*exp(-(maxA-sAge)*M)

where Rp is the predicted (number of recruits)
Rs is the actual observed number of recruits for a given stock, s
maxA is the maximum age recruit that appears for ALL stocks
sAge is the recruitment age for stock s
M is the mortality rate for recruits

In this case, I am ONLY interested in looking at how much M varies to determine if the standardization is appropriate to carry through to my stock-recruit model. I have a datafile and initial values for 2 chains. stockstart indicates the starting value of stock s, and stockfinish indicates the finishing value of stock s. I had to set it up this way so that all of my data is in one file. So when you see stockstart[Group[i]]:stockfinish[Group[i]], it simply denotes the start and finish values of stock s with the designation Group[1....12] for 12 groups. CodSR indicates the sum total of all values across 12 stocks (in this case I have 423 recruit values for 12 Groups[i]'s). NGroup indicates the total number of Groups (I'm using 12 for my model). When I run the script, I am getting an "Undefined Real Result" Trap, but I'm not sure what the problem is. Here is my code:.

#S-R Standardization

model
{

for (i in 1:NGroup)
{
M[i]~dnorm(Mo,Mprec)I(,0.000000001)
tau[i]~dnorm(mutau,tauprec)I(,0.000000001)
}

for (i in 1:CodSR)
{ #loop over all Ro-S observations for Ncod
Rs[i]~dlnorm(Rp[Group[i]],tau[Group[i]]) #evaluate prob. density of data
Rp[i]<-Rs[stockstart[Group[i]]:stockfinish[Group[i]]]*exp(-(maxA-sAge[stockstart[Group[i]]:stockfinish[Group[i]]])*M[Group[i]])
}
Mo~dlnorm(0.3,0.1)
Mprec~dgamma(0.1,0.01)
mutau~dgamma(0.01,0.01)
tauprec~dgamma(0.01,0.01)
maxA<-4

Mmu<-mean(M[]) #mean mortality of stocks
sdM<-sd(M[])
cvM<-sdM/Mmu #CV for mortality across stocks
}

I have attempted changing my starting values, truncating my distributions, providing less information for my prior distributions (gamma distributions), and keeping all my distributions lognormal. Nothing has worked so I'm pretty stuck at this point. Can anyone help?

[benarcha]

At first sight it appears you have a prior problem with M values getting negative and going pretty high. Isn't M[i]~dnorm(Mo,Mprec)I(,0.000000001) supposed to be M[i]~dnorm(Mo,Mprec)I(0.000000001,).

Note that even with this correction you will have strange values for M, going super high. therefore I would try finding a more approriate specification for M priors at first. Why not constrain variance or try an uniform distribution.

Ps : I am working on hierarchical S-R models too, not on cod tho, and would be interested in knowing how it develops for you. I am confronter to standardisation issues as I am in a multispecific context