25 views

Skip to first unread message

Jun 2, 2021, 11:00:31 AM6/2/21

to nimble-users

I am interested in using dCJS_vv for quickly simulating data for teaching exercises. The following code runs and seems to give reasonable answers for the dipper data set. But, when I try to use it to simulate, I get an error about the dimension for phi not being equal to len-1. I understand that only 6 phi's are involved here, so that makes sense. But, when I changed the code to only have phi have 6 values instead of 7, then the code wouldn't work for the actual dipper data. Any help would be greatly appreciated. Thanks!

library(nimble)

library(nimbleEcology)

library(MCMCvis)

# 'dipper_example_dir' is set to point at directory with dipper.csv

dipper <- read.csv(file.path(dipper_example_dir, "dipper.csv"))

y <- as.matrix(dipper[ , 1:7])

first <- apply(y, 1, function(x) min(which(x !=0)))

y <- y[ first != 7, ]

# re-apply first to only apply to those remaining

first <- apply(y, 1, function(x) min(which(x !=0)))

head(y)

T = ncol(y)

n.ch <- nrow(y)

first <- apply(y, 1, function(x) min(which(x !=0)))

y <- y[ first != 7, ]

# re-apply first to only apply to those remaining

first <- apply(y, 1, function(x) min(which(x !=0)))

head(y)

T = ncol(y)

n.ch <- nrow(y)

code <- nimbleCode({

for (t in 1:T){

phi[t] ~ dunif(0, 1)

p[t] ~ dunif(0, 1)

}

phi[t] ~ dunif(0, 1)

p[t] ~ dunif(0, 1)

}

for( i in 1:n.ch ){

y[i, first[i]:T] ~ dCJS_vv(probSurvive = phi[first[i]:T],

probCapture = p[first[i]:T],

len = T-first[i]+1)

}

})

dipper_inits <- function() list(phi = rep(0.1, 7),

p = rep(0.5, 7))

dipper_constants <- list(n.ch = n.ch,

T = T,

first = first)

dipper_data <- list(y = y)

samples <- nimbleMCMC(code,

constants = dipper_constants,

data = dipper_data,

dipper_inits(), samplesAsCodaMCMC = FALSE,

nchains = 3, nburnin = 1000,

niter = 6000, WAIC = FALSE)

MCMCsummary(samples)

Jun 2, 2021, 11:05:06 AM6/2/21

to nimble-users

here is the code that generates errors when I try to simulate

# Try to simulate

simModel <- nimbleModel(code = code,

constants = list(n.ch = n.ch,

T = T,

first = first),

inits = list(phi = c(0.6,0.9,0.75,0.6,0.8,0.7,1),

p = c(1,0.4,0.35,0.6,0.3,0.5,0.4)))

simModel$calculate('y')

nodesToSim <- simModel$getDependencies(c("phi", "p"),

self = F, downstream = T)

nodesToSim

simModel$simulate(nodesToSim)

simModel$y

simModel <- nimbleModel(code = code,

constants = list(n.ch = n.ch,

T = T,

first = first),

inits = list(phi = c(0.6,0.9,0.75,0.6,0.8,0.7,1),

p = c(1,0.4,0.35,0.6,0.3,0.5,0.4)))

simModel$calculate('y')

nodesToSim <- simModel$getDependencies(c("phi", "p"),

self = F, downstream = T)

nodesToSim

simModel$simulate(nodesToSim)

simModel$y

Jun 2, 2021, 12:10:16 PM6/2/21

to Jay Rotella, nimble-users

Hi Jay,

Thanks for sharing this issue. Let me know if anything in the below explanation is unclear.

The code as written won't work. The key thing is that, as you've noticed, "phi" needs to be length 6, not 7, while p needs to be length 7. In addition to changing the dimensions of your initial values, you need to change your code to conform to that dimension. Specifically, (1) you can't use the indexing "phi[first[i]:T]" any more because phi is smaller than length T, and (2) you need to change the length of the for-loop for the prior on phi (it'll need to go in a separate loop from p).

You'll notice 2 printed errors when the model is built without data but these can be safely ignored.

This was a helpful question because I didn't realize before now that the probability density calculation dCJS_vv is more lenient to this constraint than the simulation equivalent. I'll open an issue on Github; we may want to make dCJS_vv more strict.

Ben

--

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/1f658575-4d2e-4c7c-b7c6-232e60156ed7n%40googlegroups.com.

Jun 2, 2021, 12:33:09 PM6/2/21

to Rotella, Jay, nimble-users

Ah, that's illustrative.

I bet you have individuals whose first captures were at time T-1 (i.e. first[i] == 6). This will cause a problem because NIMBLE wants the probSurvive argument to be a vector, but in that case it'd be coerced to a scalar.

Unfortunately I don't have a good solution to offer here as NIMBLE would need a totally separate dCJS function to accept a scalar probSurvive. Here are two suggestions:

(1) If this is just an exercise, you could drop all individuals with first capture at time T-1, with "y <- y[ first < 6, ]"

(2) You could separate the data for these individuals and make a second for-loop in your model code that uses dCJS_sv instead of dCJS_vv, which would allow you to specify a scalar detection probability. Like:

for( i in 1:n.ch.vec ){

y_onecase[i, first[i]:T] ~ dCJS_vv(probSurvive = phi[first[i]:(T-1)],

probCapture = p[first[i]:T],

len = T-first[i]+1)

}

for( i in 1:n.ch.scal ){

y_othercase[i, first[i]:T] ~ dCJS_sv(probSurvive = phi[first[i]],

probCapture = p[first[i]:T],

len = T-first[i]+1)

}

In this example you'd provide data = list(y_onecase, y_othercase) where y_othercase contains the detection histories for individuals with first[i] == 6, and you'd need two indexing constants.

On Wed, Jun 2, 2021 at 9:25 AM Rotella, Jay <rot...@montana.edu> wrote:

I did try that one earlier and again just now (code below) … I get the following error

Dimension of distribution argument(s) 'probSurvive' does not match required dimension(s) for the distribution 'dCJS_vv'. Necessary dimension(s) are 1. You may need to ensure that you have explicit vectors and not one-row or one-column matrices.

Best!Jay

code <- nimbleCode({

for (t in 1:(T-1)){

y[i, first[i]:T] ~ dCJS_vv(probSurvive = phi[first[i]:(T-1)],

probCapture = p[first[i]:T],len = T-first[i]+1)}})

dipper_inits <- function() list(phi = rep(0.1, 6),

samples <- nimbleMCMC(code, dipper_constants, dipper_data,

dipper_inits(), samplesAsCodaMCMC = FALSE,nchains = 3, nburnin = 1000,niter = 6000, WAIC = FALSE)

On Jun 2, 2021, at 10:21 AM, Ben Goldstein <ben.go...@berkeley.edu> wrote:

probSurvive = phi[first[i]:(T-1)]?

The prior piece becomes

for (t in 1:(T-1)){

phi[t] ~ dunif(0, 1)

}

Jun 2, 2021, 3:23:17 PM6/2/21

to nimble-users

Thanks, Ben - that helped greatly. In case it's of use to anyone, here's a version that seems to work.

Jay

y <- as.matrix(dipper[ , 1:7])

first <- apply(y, 1, function(x) min(which(x !=0)))

y <- y[ first != 7, ]

# re-apply first to only apply to those remaining

first <- apply(y, 1, function(x) min(which(x !=0)))

head(y)

T = ncol(y)

n.ch <- nrow(y)

first <- apply(y, 1, function(x) min(which(x !=0)))

y <- y[ first != 7, ]

# re-apply first to only apply to those remaining

first <- apply(y, 1, function(x) min(which(x !=0)))

head(y)

T = ncol(y)

n.ch <- nrow(y)

# Individuals released on occasion T-1 will only have a single phi[t]

# and require a scalar for their phi values rather than a vector.

# Thus, need separate data and loops for individuals first releasted

# on occasion T-1 and for those released earlier

y_onecase <- y[ first == (T-1), ]

n.ch.scal <- nrow(y_onecase)

first.scal <- T-1

y_multi <- y[ first < (T-1), ]

n.ch.vec <- nrow(y_multi)

first.vec <- first[ first < (T-1)]

# and require a scalar for their phi values rather than a vector.

# Thus, need separate data and loops for individuals first releasted

# on occasion T-1 and for those released earlier

y_onecase <- y[ first == (T-1), ]

n.ch.scal <- nrow(y_onecase)

first.scal <- T-1

y_multi <- y[ first < (T-1), ]

n.ch.vec <- nrow(y_multi)

first.vec <- first[ first < (T-1)]

code <- nimbleCode({

for (t in 1:(T-1)){

phi[t] ~ dunif(0, 1)

}

for (t in 1:T){

p[t] ~ dunif(0, 1)

}

for( i in 1:n.ch.vec ){

y_multi[i, first.vec[i]:T] ~ dCJS_vv(probSurvive = phi[first.vec[i]:(T-1)],

probCapture = p[first.vec[i]:T],

len = T-first.vec[i]+1)

y_multi[i, first.vec[i]:T] ~ dCJS_vv(probSurvive = phi[first.vec[i]:(T-1)],

probCapture = p[first.vec[i]:T],

len = T-first.vec[i]+1)

}

for( i in 1:n.ch.scal ){

y_onecase[i, first.scal:T] ~ dCJS_sv(probSurvive = phi[first.scal],

probCapture = p[first.scal:T],

len = T-first.scal+1)

}

})

dipper_inits <- function() list(phi = runif((T-1), 0, 1),

p = runif(T, 0, 1))

dipper_constants <- list(n.ch.vec = n.ch.vec,

n.ch.scal = n.ch.scal,

first.vec = first.vec,

first.scal = first.scal,

T = T)

dipper_data <- list(y_onecase = y_onecase,

y_multi = y_multi)

probCapture = p[first.scal:T],

len = T-first.scal+1)

}

})

dipper_inits <- function() list(phi = runif((T-1), 0, 1),

p = runif(T, 0, 1))

dipper_constants <- list(n.ch.vec = n.ch.vec,

n.ch.scal = n.ch.scal,

first.vec = first.vec,

first.scal = first.scal,

T = T)

dipper_data <- list(y_onecase = y_onecase,

y_multi = y_multi)

samples <- nimbleMCMC(code, dipper_constants, dipper_data,

dipper_inits(),

nchains = 3,

nburnin = 1000,

nburnin = 1000,

niter = 16000,

#thin = 10,

WAIC = FALSE,

samplesAsCodaMCMC = FALSE)

#thin = 10,

WAIC = FALSE,

samplesAsCodaMCMC = FALSE)

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu