14 views

Skip to first unread message

May 28, 2021, 7:16:10 AMMay 28

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"))

}

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"))

}

May 28, 2021, 8:49:17 PMMay 28

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.

May 30, 2021, 12:56:32 PMMay 30

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)

Jun 23, 2021, 1:51:26 PMJun 23

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

Jun 25, 2021, 3:29:36 AMJun 25

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

Jul 3, 2021, 7:24:07 PMJul 3

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

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/ef4fb38c-756b-4a50-b529-bf2ddcb25d0fn%40googlegroups.com.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu