Failed to create a loop over different samples with CJS model with individual heterogeneity

48 views
Skip to first unread message

blanca sarzo

unread,
May 28, 2021, 7:16:10 AM5/28/21
to nimble-users
Good afternoon,

I am trying to run a "simple" CJS model with individual heterogeneity with NIMBLE. The code works fine with one sample, however, when I include a loop at the beginning of the code to run it for several different samples a get different errors or even R crashes.
I am new on NIMBLE so I am not sure what is happening or whereas I have to indicate that loop differently. Any help would be appreciated. Thanks.
Below is the code:

code1 <- nimbleCode( {
 
  for (i in 1:nind){
    for (t in f[i]:(n.occasions-1)){
      logit(phi[i,t]) <- alpha[x[i,t]] + beta[t] + epsilon[i]
      p[i,t] <- alpha.p[x[i,t]]
    }
  }
 
  # Priors
  for (i in 1:nind){
    epsilon[i] ~ dnorm(0,tau)
  }
 
  ## To avoid the lack of identifiability
  ## we set alpha[1]=0
 
  for (t in 1:(n.occasions-1)){
    beta[t] ~ dnorm(mu, tau.beta)
  }
  tau.beta <- 1/(sigma.s*sigma.s)
  sigma.s ~ dunif(0,10)  
  mu ~ dnorm(0,0.1)
 
  alpha[1] <- 0

  for (u in 2:4){
    alpha[u] ~ dnorm(0,0.25)
  }
 
  for (u in 1:4){
    alpha.p[u] ~ dbeta(1,1)
  }
 
  sigma ~ dunif(0, 2)
  tau <- 1/(sigma*sigma)
 
  # Likelihood
  for (i in 1:nind){
    
    # Define latent state at first capture
    z[i,f[i]] ~ dbern(1)
    
    for (t in (f[i]+1):n.occasions){
      z[i,t] ~ dbern(mu1[i,t])
      mu1[i,t] <- phi[i,t-1]*z[i,t-1]
      
      # Observation process
      y[i,t] ~ dbern(mu2[i,t])
      mu2[i,t] <- p[i,t-1]*z[i,t]
    }
  }
}
)

################################################################
## Data, constants, initials
################################################################

set.seed(1)

## Functions
last.seen <- function(x) max(which(x!=0))
get.first <- function(x) min(which(x!=0))

for (s in 1:3){ ## Here is the problem!!!!

sub1 <- readRDS(paste0("subsample.real_", s, ".RDS"))

CH <- as.matrix(sub1[,-c(1,13:14)])

## Ringing vector
f <- apply(CH, 1, get.first)
l <- apply(CH, 1, last.seen)

## Create matrices X indicating age classes

x <- matrix(NA, ncol=dim(CH)[2]-1, nrow=dim(CH)[1])
vec.age <- c(1,2, 3, rep(4,7))
f.inverse <- dim(CH)[2]- f
for(i in 1:dim(CH)[1]){
  x[i,f[i]:dim(x)[2]] <- vec.age[1:f.inverse[i]]
}

## Known values of the latent variable Z (1's, 0's, NAs)
known.state.cjs <- function(ch){
  state <- ch
  for (i in 1:dim(ch)[1]){
    n1 <- min(which(ch[i,]==1))         
    n2 <- max(which(ch[i,]==1))         
    state[i, n1:n2] <- 1             
  }
  state[state==0] <- NA                 
  return(state)
}

## Function  to create initial values of Z

cjs.init.z <- function(ch,f,l){
  inz <- ch
  for(i in 1:dim(ch)[1]){
    if(sum(ch[i,])==1) next
    n2 <- max(which(ch[i,]==1))
    inz[i, f[i]:n2] <- 1
  }
  for(i in 1: dim(ch)[1]){
    inz[i, 1:f[i]] <- NA
  }
 
  # To randomly add in time of death   
 
  for(i in 1: dim(ch)[1]){
    if (l[i] != dim(ch)[2]) {
      death <- sample(c(l[i]:dim(ch)[2]),1)
      if (death != l[i]){
        for (j in (l[i]+1):death){
          inz[i,j] <- 1
        }
      }
    }
  }
 
  return(inz)
}

z <- known.state.cjs(CH)
z2 <- cjs.init.z(CH, f, l)
nind <- dim(CH)[1]
n.occasions <- dim(CH)[2]
y <- CH

inepsilon1 <- rnorm(nind,0,1)
inbeta <- rep(0,(n.occasions-1))

jags.data <- list(y=y,z=z)

constants <- list(nind=nind,n.occasions=n.occasions,x=x,f=f)

### Initial values
inits <- list(sigma=1, alpha=c(NA,2,1,1),alpha.p=c(0.02,0.2,0.4,0.4),z=z2,epsilon=inepsilon1,beta=inbeta)

### Parameters
parameters <- c("sigma", "alpha", "alpha.p","beta","sigma.s","mu")  

NIMBLE_mod <- nimbleModel(code1, data=jags.data,
                          constants = constants,
                          inits = inits)

NIMBLE_real <- nimbleMCMC(code1, data=jags.data,constants = constants,
                          inits = inits, niter = 35000, nburnin = 1500,
                          thin = 1, nchains=1, monitors=parameters)

saveRDS(NIMBLE_real, paste0("NIMBLE_real", s, ".RDS"))

}


Perry de Valpine

unread,
May 28, 2021, 8:49:17 PM5/28/21
to blanca sarzo, nimble-users
Hi Blanca,
Thanks for posting this.
Can you please share the error messages and/or data to reproduce what you're seeing?
-Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a02e25d2-7757-4ca3-81df-84ec935ee97cn%40googlegroups.com.

Blanca Sarzo

unread,
May 30, 2021, 12:56:32 PM5/30/21
to Perry de Valpine, nimble-users

Hi Perry,

This is the error I have and attached two toy examples of my databases. What I do not understand is that with only one sample it works perfectly well, however, if I include the loop is when I get errors.

Thanks again,

Blanca

Error message:

defining model...
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ...
checking model sizes and dimensions... This model is not fully initialized. This is not an error. To see which variables are not initialized, use model$initializeInfo(). For more information on model initialization, see help(modelInitialization).
checking model calculations...
NAs were detected in model variables: sigma.s, logProb_sigma.s, mu, logProb_mu, tau.beta, lifted_d1_over_sqrt_oPtau_dot_beta_cP, logProb_beta.
model building finished.
compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors.
Además: Warning messages:
1: Ignoring non-NA values in inits for data nodes: z.
2: In system2(SHLIBcmd, SHLIBargs, stdout = logFile, stderr = errorFile) :
  system call failed: Cannot allocate memory
3: In system2(SHLIBcmd, SHLIBargs, stdout = logFile, stderr = errorFile) :
  error in running command


El 29/5/21 a las 2:49, Perry de Valpine escribió:
-- 
Blanca Sarzo Carles
PhD Statistics
Foundation for the Promotion of Health and Biomedical Research of Valencia Region (FISABIO)
sub_1.RDS
sub_2.RDS

Perry de Valpine

unread,
Jun 23, 2021, 1:51:26 PM6/23/21
to Blanca Sarzo, nimble-users
Hi Blanca,

Sorry this fell through the cracks.  We've tried to keep up our support efforts on nimble-users amid summer travels, but it came to my attention that your issue was left dangling.  Did you find a solution or at least a work-around?  When I ran your code, replacing paste0("subsample.real_", s, ".RDS") with "sub_1.RDS", one of the files you sent, I got to the point that creation of f and l generated errors, putting some values of Inf or -Inf in the result, and I was not sure they are well formed.  Is that correct?

On studying your code, IIUC the issue is having a for-loop that contains nimbleMCMC.  Usually people are asking about for loops in their model code, so my thoughts had started in the wrong direction.  A side comment is that you do not need to do nimbleModel before nimbleMCMC in your code.

Crashes are sometimes easy but sometimes hard to reproduce.  Can you let me know if this is an ongoing issue and if so can you guide me on reproducing it?  I wasn't sure if have Inf and -Inf values in an indexing vector would make sense.

Thanks.
-Perry 

blanca sarzo

unread,
Jun 25, 2021, 3:29:36 AM6/25/21
to nimble-users
Hi again Perry,

Sorry, but it workeed with 2 samples, but with more than 2 R crashes. The "real" size of the subsamples is a matrix with dimensions 5789 x 11, also, the number of iterations I am using is 35,000 with 1,500 of burn-in... Besides, in the code I sent you there is a bug and that is why we got -Inf values, below the corrected code:
CH <- as.matrix(sub1)

The code again works perfectly well for 1 sample and at most 2 at a time, but when I insert the loop everything crashes...

Thanks again for your help,

Blanca

Perry de Valpine

unread,
Jul 3, 2021, 7:24:07 PM7/3/21
to blanca sarzo, nimble-users
Hi Blanca,

I'm not able to reproduce this error.

A possible workaround would be to launch each iteration via a system call (using system() or system2()) to Rscript, which will run a separate R session.  That might avoid some of the kinds of issues that in general could end up with a crash, but I'm guessing here.  You can give the Rscript call a code file to run and other arguments.  I hope that helps.  I'm not sure what we can do to pin down your error more.  Possibly if you provide all the data files, we would then be running exactly what you are.  When I tried it, I was repeatedly loading the same one of the two data files you sent before.  It's hard to imagine what difference it would make, but when we have mysterious bugs it can help to start from doing *exactly* what is generating the bug.  However, if the Rscript method works, I would go with that for now.

-Perry


Reply all
Reply to author
Forward
0 new messages