Issues with a latent variable model

[Drew Tyre]

Hi all,

I'm running a model with a large latent variable matrix representing
population size in each year for several size classes.

var lambda[T], N[T,A];
model {
    # Constant p, constant S model
    # main matrix remains in age by year form
    # N[] is the matrix of latent population size
    # C[] is the matrix of losses by age and year (observed)
    # T is the maximium years, A is maximum age class 
    # X[] is a vector of total observed mortality for initial cohort

    # initial year handled differently 
    for (a in 1:A){
        C[1,a] ~ dbin(1/(1+exp(-p0)),trunc(N[1,a]))
    } # end loop over A
    lambda[1] <- sum((N[1,2:A] - C[1,2:A])*exp(f0))

    for (t in 2:T){
        N[t,1] ~ dpois(lambda[t-1]) # recruitment
        for (a in 2:(A-1)){
            N[t,a] ~ dbin(1/(1+exp(-S0)),trunc(N[t-1,a-1]))
            C[t,a] ~ dbin(1/(1+exp(-p0)),trunc(N[t,a]))
        }   # end loop over A
        lambda[t] <- sum((N[t,2:A] - C[t,2:A])*exp(f0))
        N[t,A] ~ dbin(1/(1+exp(-S0)),trunc(N[t-1,A-1]+N[t-1,A]))
        C[t,A] ~ dbin(1/(1+exp(-p0)),trunc(N[t,A]))

    }   # end loop over T

    # priors on N for all ages at initial time? 
    for (a in 1:A){
        N[1,a] ~ dunif(X[a],1000)   
    }
    # Other priors
    p0 ~ dnorm(-1.15268,13.7) # prior drawn from Cooley etal 2009, Hunting mortality rate for heavily hunted area
    S0 ~ dnorm(1.516347,13.7) # prior drawn from Cooley etal 2009, 1-natural mortality rate for heavily hunted area
    f0 ~ dnorm(0,1)
}

I run the model using rjags and the following script

load("testdata.Rdata")
# Run basic model
if (!require(rjags)) stop("Install the rjags package first.")

data.list = list(C = C, T = T, A = A, X = X)
inits.list = list(list(N=matrix(rep(initN-50,times=T),byrow=TRUE,ncol=A),p0=-1,S0=1,f0=0),
                list(N=matrix(rep(initN,times=T),byrow=TRUE,ncol=A),p0=-2,S0=2,f0=1),
                list(N=matrix(rep(initN+50,times=T),byrow=TRUE,ncol=A),p0=0,S0=0,f0=-1))

# setup to run 3 chains
model.MpS <- jags.model("tyre-jagsmodel2-2012-5-24.r",data=data.list,inits=inits.list,n.chains=3)

# this works
MpS.samples = coda.samples(model.MpS,c("N","p0","S0","f0"),1000) # looks good

# but I get warnings and errors on 
summary(MpS.samples)
# and there are missing values in the time series SE column

# takes a while - go get coffee
MpS.samples = coda.samples(model.MpS,c("N","p0","S0","f0"),100000)
# this crashes everything (including Rgui) at the end of the run

I run into a few issues when I run the model, especially if I do a very long
run. Short runs give all sorts of warnings and errors when I look at the
summary. Long runs (100000 iterations) are needed because this model doesn't
mix well, but this version crashes Rgui at the end of long runs.

One strange thing - when I use list.samplers(model.MpS) I find that one entry
in the matrix N has no sampler associated with it - N. This node doesn't seem
to update after I call jags.model(...) - it still has the initial value
associated with it. However it does seem to update after calling
coda.samples(...) .

Slaps on the wrist, comments, insights all welcome. Thank you.

The file 'testdata.Rdata' (0 KB) is available for download at

< http://dropbox.unl.edu/uploads/20120613/a29f8360a6188769/testdata.Rdata >

for the next 14 days.

It will be removed after Wednesday, June 13, 2012.

[Drew Tyre]

The error message when Rgui crashes is from the Microsoft C++ library, and
says something like "the application has requested a shutdown in an unusual
way".

Very helpful!

[Drew Tyre]

OK, so I'm glad no one has responded, because I'm a victim of my own
ignorance. The crash at the end is a memory problem - if I aggresively thin
the samples it goes away. I have to do that anyway, so no harm done. The issue
with the missing samplers has to do with me goofing and not including part of
the observed data for the first age class.

Now that I've fixed that, I have other issues, but JAGS is now asking me to
submit a bug report, so I'll go do that!