Marginalization - dimension mismatch for 'dHMM' distribution (nimble ecology)

Sarah Durham

unread,

Jul 14, 2022, 5:31:01 PM7/14/22

to nimble-users

Hello,

I am attempting to marginalize my multi-state capture recapture model in order to speed up computation time, but I'm running into an error in model$checkBasics() regarding dimensionality. It reads:

Error in model$checkBasics() :
Dimension of distribution argument(s) 'value' does not match required dimension(s) for the distribution 'dHMM'. Necessary dimension(s) are 1. You may need to ensure that you have explicit vectors and not one-row or one-column matrices.

I can't seem to figure out what the issue is. I've tried feeding in the data different ways and I keep getting the same error. I can't seem to figure out which specific argument(s) are the dimension distribution argument(s). Any suggestions would be greatly appreciated! The model code is pasted below along with the data simulation.

Cheers!

##############################################

library(nimble)
library(nimbleEcology)

# Generation of simulated data
# Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals
phiA <- 0.69 # Apparent survival probability for juveniles
phiB <- 0.80 # Apparent survival probability for sub-adults
phiC <- 0.90 # Apparent survival probability for adults
pA <- 0.10 # Recapture probability for juveniles
pB <- 0.20 # Recapture probability for sub-adults
pC <- 0.30 # Recapture probability for adults

n.occasions <- 25
n.states <- 9
n.obs <- 9
marked <- matrix(NA, ncol = n.states, nrow = n.occasions)
marked[,1] <- rep(300, n.occasions)
marked[,2] <- rep(0, n.occasions)
marked[,3] <- rep(0, n.occasions)
marked[,4] <- rep(0, n.occasions)
marked[,5] <- rep(0, n.occasions)
marked[,6] <- rep(0, n.occasions)
marked[,7] <- rep(0, n.occasions)
marked[,8] <- rep(0, n.occasions)
marked[,9] <- rep(0, n.occasions)

# Define matrices with survival, transition and recapture probabilities
# These are 4-dimensional matrices, with
# Dimension 1: state of departure
# Dimension 2: state of arrival
# Dimension 3: individual
# Dimension 4: time
# 1. State process matrix
totrel <- sum(marked)*(n.occasions-1)
PSI.STATE <- array(NA, dim=c(n.states, n.states, totrel, n.occasions-1))
for (i in 1:totrel){
for (t in 1:(n.occasions-1)){
PSI.STATE[,,i,t] <- matrix(c(
0,phiA,0,0,0,0,0,0,1-phiA,
0,0,phiA,0,0,0,0,0,1-phiA,
0,0,0,phiA,0,0,0,0,1-phiA,
0,0,0,0,phiB,0,0,0,1-phiB,
0,0,0,0,0,phiB,0,0,1-phiB,
0,0,0,0,0,0,phiB,0,1-phiB,
0,0,0,0,0,0,0,phiB,1-phiB,
0,0,0,0,0,0,0,phiC,1-phiC,
0,0,0,0,0,0,0,0,1), nrow = n.states, byrow = TRUE)
} #t
} #i

# 2.Observation process matrix
PSI.OBS <- array(NA, dim=c(n.states, n.obs, totrel, n.occasions-1))
for (i in 1:totrel){
for (t in 1:(n.occasions-1)){
PSI.OBS[,,i,t] <- matrix(c(
1,0,0,0,0,0,0,0,0,
0,pA,0,0,0,0,0,0,1-pA,
0,0,pA,0,0,0,0,0,1-pA,
0,0,0,pB,0,0,0,0,1-pB,
0,0,0,0,pB,0,0,0,1-pB,
0,0,0,0,0,pB,0,0,1-pB,
0,0,0,0,0,0,pB,0,1-pB,
0,0,0,0,0,0,0,pC,1-pC,
0,0,0,0,0,0,0,0,1), nrow = n.states, byrow = TRUE)
} #t
} #i

# Define function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, marked, unobservable = NA){
# Unobservable: number of state that is unobservable
n.occasions <- dim(PSI.STATE)[4] + 1
CH <- CH.TRUE <- matrix(NA, ncol = n.occasions, nrow = sum(marked))
# Define a vector with the occasion of marking
mark.occ <- matrix(0, ncol = dim(PSI.STATE)[1], nrow = sum(marked))
g <- colSums(marked)
for (s in 1:dim(PSI.STATE)[1]){
if (g[s]==0) next # To avoid error message if nothing to replace
mark.occ[(cumsum(g[1:s])-g[s]+1)[s]:cumsum(g[1:s])[s],s] <-
rep(1:n.occasions, marked[1:n.occasions,s])
} #s
for (i in 1:sum(marked)){
for (s in 1:dim(PSI.STATE)[1]){
if (mark.occ[i,s]==0) next
first <- mark.occ[i,s]
CH[i,first] <- s
CH.TRUE[i,first] <- s
} #s
for (t in (first+1):n.occasions){
# Multinomial trials for state transitions
if (first==n.occasions) next
state <- which(rmultinom(1, 1, PSI.STATE[CH.TRUE[i,t-1],,i,t-1])==1)
CH.TRUE[i,t] <- state
# Multinomial trials for observation process
event <- which(rmultinom(1, 1, PSI.OBS[CH.TRUE[i,t],,i,t-1])==1)
CH[i,t] <- event
} #t
} #i
# Replace the NA and the highest state number (dead) in the file by 0
CH[is.na(CH)] <- 0
CH[CH==dim(PSI.STATE)[1]] <- 0
CH[CH==unobservable] <- 0
id <- numeric(0)
for (i in 1:dim(CH)[1]){
z <- min(which(CH[i,]!=0))
ifelse(z==dim(CH)[2], id <- c(id,i), id <- c(id))
}
return(list(CH=CH[-id,], CH.TRUE=CH.TRUE[-id,]))
# CH: capture histories to be used
# CH.TRUE: capture histories with perfect observation
}

# Execute simulation function
sim <- simul.ms(PSI.STATE, PSI.OBS, marked)
CH <- sim$CH

# Compute vector with occasion of first capture
get.first <- function(x) min(which(x!=0))
first <- apply(CH, 1, get.first)

# Recode CH matrix: note, a 0 is not allowed!
rCH <- CH # Recoded CH
rCH[rCH==0] <- 9

rCH_fun <- function(ch, f){
z <- ch
for (i in 1:nrow(ch)) {
for (j in 1:ncol(ch)) {
if (j < (f[i])){z[i,j] <- NA}

}
}
z <- as.matrix(z)
return(z)
}

y=as.matrix(rCH_fun(rCH,first))

multisite.marginalized <- nimbleCode({
# ---------------------------------
# Parameters:
# phiA: survival probability for juveniles
# phiB: survival probability for sub-adults
# phiC: survival probability for adults
# pA: recapture probability of juvenile
# pB: recapture probability of sub-adults
# pC: recapture probability of adults
# ---------------------------------
# States (S):
# 1 chick
# 2 age 1 bird
# 3 age 2 bird
# 4 age 3 bird
# 5 age 4 bird
# 6 age 5 bird
# 7 age 6 bird
# 8 age 7+ bird
# 9 Dead

# Observations (O):
# 1 seen as chick
# 2 seen as age 1 bird
# 3 seen as age 2 bird
# 4 seen as age 3 bird
# 5 seen as age 4 bird
# 6 seen as age 5 bird
# 7 seen as age 6 bird
# 8 seen as age 7+ bird
# 9 not seen
# ---------------------------------
# Priors and constraints
phiA ~ dunif(0, 1)
phiB ~ dunif(0, 1)
phiC ~ dunif(0, 1)
pA ~ dunif(0, 1)
pB ~ dunif(0, 1)
pC ~ dunif(0, 1)

# Define state-transition and observation matrices
# Define probabilities of state S(t+1) given S(t)

gamma[1,1] <- 0
gamma[1,2] <- phiA
gamma[1,3] <- 0
gamma[1,4] <- 0
gamma[1,5] <- 0
gamma[1,6] <- 0
gamma[1,7] <- 0
gamma[1,8] <- 0
gamma[1,9] <- 1-phiA

gamma[2,1] <- 0
gamma[2,2] <- 0
gamma[2,3] <- phiA
gamma[2,4] <- 0
gamma[2,5] <- 0
gamma[2,6] <- 0
gamma[2,7] <- 0
gamma[2,8] <- 0
gamma[2,9] <- 1-phiA

gamma[3,1] <- 0
gamma[3,2] <- 0
gamma[3,3] <- 0
gamma[3,4] <- phiA
gamma[3,5] <- 0
gamma[3,6] <- 0
gamma[3,7] <- 0
gamma[3,8] <- 0
gamma[3,9] <- 1-phiA

gamma[4,1] <- 0
gamma[4,2] <- 0
gamma[4,3] <- 0
gamma[4,4] <- 0
gamma[4,5] <- phiB
gamma[4,6] <- 0
gamma[4,7] <- 0
gamma[4,8] <- 0
gamma[4,9] <- 1-phiB

gamma[5,1] <- 0
gamma[5,2] <- 0
gamma[5,3] <- 0
gamma[5,4] <- 0
gamma[5,5] <- 0
gamma[5,6] <- phiB
gamma[5,7] <- 0
gamma[5,8] <- 0
gamma[5,9] <- 1-phiB

gamma[6,1] <- 0
gamma[6,2] <- 0
gamma[6,3] <- 0
gamma[6,4] <- 0
gamma[6,5] <- 0
gamma[6,6] <- 0
gamma[6,7] <- phiB
gamma[6,8] <- 0
gamma[6,9] <- 1-phiB

gamma[7,1] <- 0
gamma[7,2] <- 0
gamma[7,3] <- 0
gamma[7,4] <- 0
gamma[7,5] <- 0
gamma[7,6] <- 0
gamma[7,7] <- 0
gamma[7,8] <- phiB
gamma[7,9] <- 1-phiB

gamma[8,1] <- 0
gamma[8,2] <- 0
gamma[8,3] <- 0
gamma[8,4] <- 0
gamma[8,5] <- 0
gamma[8,6] <- 0
gamma[8,7] <- 0
gamma[8,8] <- phiC
gamma[8,9] <- 1-phiC

gamma[9,1] <- 0
gamma[9,2] <- 0
gamma[9,3] <- 0
gamma[9,4] <- 0
gamma[9,5] <- 0
gamma[9,6] <- 0
gamma[9,7] <- 0
gamma[9,8] <- 0
gamma[9,9] <- 1

# Define probabilities of O(t) given S(t)
omega[1,1] <- 1
omega[1,2] <- 0
omega[1,3] <- 0
omega[1,4] <- 0
omega[1,5] <- 0
omega[1,6] <- 0
omega[1,7] <- 0
omega[1,8] <- 0
omega[1,9] <- 0

omega[2,1] <- 0
omega[2,2] <- pA
omega[2,3] <- 0
omega[2,4] <- 0
omega[2,5] <- 0
omega[2,6] <- 0
omega[2,7] <- 0
omega[2,8] <- 0
omega[2,9] <- 1-pA

omega[3,1] <- 0
omega[3,2] <- 0
omega[3,3] <- pA
omega[3,4] <- 0
omega[3,5] <- 0
omega[3,6] <- 0
omega[3,7] <- 0
omega[3,8] <- 0
omega[3,9] <- 1-pA

omega[4,1] <- 0
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- pB
omega[4,5] <- 0
omega[4,6] <- 0
omega[4,7] <- 0
omega[4,8] <- 0
omega[4,9] <- 1-pB

omega[5,1] <- 0
omega[5,2] <- 0
omega[5,3] <- 0
omega[5,4] <- 0
omega[5,5] <- pB
omega[5,6] <- 0
omega[5,7] <- 0
omega[5,8] <- 0
omega[5,9] <- 1-pB

omega[6,1] <- 0
omega[6,2] <- 0
omega[6,3] <- 0
omega[6,4] <- 0
omega[6,5] <- 0
omega[6,6] <- pB
omega[6,7] <- 0
omega[6,8] <- 0
omega[6,9] <- 1-pB

omega[7,1] <- 0
omega[7,2] <- 0
omega[7,3] <- 0
omega[7,4] <- 0
omega[7,5] <- 0
omega[7,6] <- 0
omega[7,7] <- pB
omega[7,8] <- 0
omega[7,9] <- 1-pB

omega[8,1] <- 0
omega[8,2] <- 0
omega[8,3] <- 0
omega[8,4] <- 0
omega[8,5] <- 0
omega[8,6] <- 0
omega[8,7] <- 0
omega[8,8] <- pC
omega[8,9] <- 1-pC

omega[9,1] <- 0
omega[9,2] <- 0
omega[9,3] <- 0
omega[9,4] <- 0
omega[9,5] <- 0
omega[9,6] <- 0
omega[9,7] <- 0
omega[9,8] <- 0
omega[9,9] <- 1

#Initial state probabilities

delta[1] <- 1 # Pr(alive in state 1 at t = 1) = 1
delta[2] <- 0 # Pr(alive in state 2 at t = 1) = 0
delta[3] <- 0 # Pr(alive in state 3 at t = 1) = 0
delta[4] <- 0 # Pr(alive in state 4 at t = 1) = 0
delta[5] <- 0 # Pr(alive in state 5 at t = 1) = 0
delta[6] <- 0 # Pr(alive in state 6 at t = 1) = 0
delta[7] <- 0 # Pr(alive in state 7 at t = 1) = 0
delta[8] <- 0 # Pr(alive in state 8 at t = 1) = 0
delta[9] <- 0 # Pr(alive in state 9 at t = 1) = 0

# initial state probs
for(i in 1:N) {
init[i, 1:9] <- delta[1:9] # First state propagation
}
# likelihood
for (i in 1:N){
y[i,(first[i]+1):K] ~ dHMM(init = init[i,1:9], # count data from first[i] + 1
probObs = omega[1:9,1:9], # observation matrix
probTrans = gamma[1:9,1:9], # transition matrix
len = K - first[i], # nb of occasions
checkRowSums = 0) # do not check whether elements in a row sum tp 1
}
})

my.data <- list(y = y)

my.constants <- list(first = first,
K = ncol(y),
N = nrow(y))

initial.values <- function(){list(phiA = runif(1, 0, 1),
phiB = runif(1, 0, 1),
phiC = runif(1, 0, 1),
pA = runif(1, 0, 1),
pB = runif(1, 0, 1),
pC = runif(1, 0, 1))}

parameters.to.save <- c("phiA", "phiB", "phiC","pA", "pB", "pC")

n.iter <- 1
n.burnin <-0
n.chains <- 3

system.time(multisite.marginalized.out <- nimbleMCMC(code = multisite.marginalized,
constants = my.constants,
data = my.data,
inits = initial.values(),
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains))

Daniel Turek

unread,

Jul 15, 2022, 7:57:59 AM7/15/22

to Sarah Durham, nimble-users

Sarah, that's great that you're using nimbleEcology, and I hope it helps out with your analysis.

The short answer is that the "value" argument to the dHMM distribution (and any other distribution in nimble model code) is the node appearing on the left-hand-side of the stochastic "~" declaration. So in your case:

y[i,(first[i]+1):K] ~ dHMM(init = init[i,1:9], # count data from first[i] + 1
probObs = omega[1:9,1:9], # observation matrix
probTrans = gamma[1:9,1:9], # transition matrix
len = K - first[i], # nb of occasions
checkRowSums = 0) # do not check whether elements in a row sum tp 1

the "value" argument to dHMM is "y[i,(first[i]+1):K]".

What nimble is trying to tell you is that the dHMM distribution expects a dimension=1 node for the "value", which means a vector. Generally, "y[i,(first[i]+1):K]" will be a vector, *except* when first[i] = K for some value(s) of i, in which case that reduces to a single value, or a scalar, or a dimension=0 quantity. Unfortunately, at this point, nimble cannot distinguish between dimension=0 scalars, and vectors of length 1, which is causing the error you're seeing.

What you'll have to do to get around this, since capture histories of length=1 generally do not contribute anything to inference, is to remove the capture histories for which first[i]=K from the dataset - remove the individuals that were *first* sighted on the final sampling occasion K.

This generally won't affect inference, unless you were also doing inference on the initial state probabilities "init[1:9]", which it appears you are not, since those appear to be hard-coded as c(1,0,0,0,....) in your model, meaning you condition on the initial state of capture being state 1, and also individuals being observed in that first time period. So removing these individuals won't affect the inferences from your model, it will just require a little bit of data manipulation, and changing the value of N.

Does this make sense? Please let us know if this gets you running.

Cheers,

Daniel

--
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/2146978f-9022-4728-9f9f-adf04ce45b7en%40googlegroups.com.

Sarah Durham

unread,

Jul 15, 2022, 9:56:55 AM7/15/22

to nimble-users

Hi Daniel,

Thanks so much for this! That definitely makes sense but when I add the code in there to do that I still get the same error message. When I look back at the raw capture histories I simulated there weren't any individuals that were simulated to have been first captured at the last occasion, so it seems like that shouldn't be the issue but I still can't seem to figure out where the issue is coming from.

Cheers!

Sarah

Daniel Turek

unread,

Jul 15, 2022, 10:32:44 AM7/15/22

to Sarah Durham, nimble-users

I answered a little quickly, but briefly noting a correction: this would be a problem for your model when (first[i]+1) = K.

Generally, for dHMM distribution, I think you might want the declaration to be y[i, first[i] : K], but double-check the documentation for it. Sorry, I can't look in more detail right now.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/75e2ab02-9e2f-4bfb-882d-3b9deb6ad0b1n%40googlegroups.com.

Sarah Durham

unread,

Jul 18, 2022, 11:41:13 AM7/18/22

to nimble-users

Okay, so I've been working on this and I think I figured out a workaround for my particular problem. I ran the model above with my actual data collected in the field and it worked. The only difference between the simulated and actual data is that we aren't analyzing individuals less than 5 years ago so they are removed from the study set. Sure enough, when I did this with the simulated data the marginalized model worked just fine. Not really sure why this is given that the example I was following didn't remove any individuals but those first marked on the last occasion, but at least it's a workaround for me for now. Thanks for your help!

Cheers!

Ben Goldstein

unread,

Jul 18, 2022, 4:40:33 PM7/18/22

to Sarah Durham, nimble-users

Hi Sarah--

I think Daniel's response is exactly correct. dHMM doesn't condition on first capture, meaning that we expect you to have at least two observations of each individual. Since a one-observation individual has a simpler likelihood you could probably handle those separately; that's another workaround that could work for you. Or, you could include the first instance of each individual and provide y[i, first[i] : K] as Daniel suggested. If you want to condition on first capture I think you could do that pretty easily by hard-coding your detection matrix for the first event for each individual. Let me know if this is unclear and if you want any support going in one of these directions!

Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/937a150a-7e70-4163-ba02-fa5b809afdcbn%40googlegroups.com.

Sarah Durham

unread,

Jul 25, 2022, 11:42:23 AM7/25/22

to nimble-users

Hi Ben!

Thanks so much for this! I'll be sure to reach out if I have any issues. As a sidenote, when marginalzing, am I unable to calculate WAIC? This seems to be the case since we aren't monitoring the latent state, but I just wanted to ask in case there was a workaround.

Cheers!

Sarah

Perry de Valpine

unread,

Jul 25, 2022, 11:57:43 AM7/25/22

to Sarah Durham, nimble-users

Hi Sarah,

Great question. There are different levels of AIC. If you calculate WAIC using a marginal probability distribution (i.e. dHMM) as a likelihood, you would get a corresponding marginal WAIC. The WAIC would be based on predictions of capture histories from parameters, not conditioned on discrete latent states within the capture series, if that makes sense. In many situations, the marginal version would be more useful in terms of selecting on the level of the model of interest. If you want to read more, there are some references in section 7.7 of our user manual. (Those refer to calculating marginal WAIC by simulation when the likelihood used is not marginalized, but the articles might still be helpful on the distinction involved.)

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4a122d5d-5d0a-4566-887d-dc64e05406a2n%40googlegroups.com.

Sarah Durham

unread,

Jul 27, 2022, 12:16:06 PM7/27/22

to nimble-users

Hi Perry,

Thanks so much for the information! Would it not make sense to calculate WAIC post hoc here (like example below)? That seems like the simplest solution but I see that the manual says it only provides conidtional WAIC without grouping of data nodes. I'm unclear if I need to group the nodes since I'm using the marginalized distribution.

Cheers!

##################################################################################################################################################################

library(nimble)
library(nimbleEcology)

# Multistate Capture-mark-recapture model

y <- as.matrix(CH) # Recoded CH
y[y==0] <- 9

# Remove individuals marked in last 5 years of study
# B/c 1) birds haven't had chance to recruit (same as we did with ATPU data) and 2) causes issues with nimbleEcology
mask <- which(first!=ncol(y) & first!=ncol(y)-1 & first!=ncol(y)-2 & first!=ncol(y)-3& first!=ncol(y)-4 & first!=ncol(y)-5)
y <- y[mask, ] # keep only these
first <- first[mask]

for(i in 1:N) {
init[i, 1:9] <- gamma[ y[i, first[i] ] , 1:9 ]# First state propagation

}
# likelihood
for (i in 1:N){
y[i,(first[i]+1):K] ~ dHMM(init = init[i,1:9], # count data from first[i] + 1
probObs = omega[1:9,1:9], # observation matrix
probTrans = gamma[1:9,1:9], # transition matrix
len = K - first[i], # nb of occasions
checkRowSums = 0) # do not check whether elements in a row sum tp 1
}
})

my.data <- list(y = y)

my.constants <- list(first = first,
K = ncol(y),
N = nrow(y))

initial.values <- function(){list(phiA = runif(1, 0, 1),
phiB = runif(1, 0, 1),
phiC = runif(1, 0, 1),
pA = runif(1, 0, 1),
pB = runif(1, 0, 1),
pC = runif(1, 0, 1))}

parameters.to.save <- c("phiA", "phiB", "phiC","pA", "pB", "pC")

n.iter <- 100000
n.burnin <-10000
n.thin <- 10
n.chains <- 3

myModel <- nimbleModel(code = multisite.marginalized,

constants = my.constants,
data = my.data,

inits = initial.values())

CmyModel <- compileNimble(myModel)

myMCMC <- buildMCMC(CmyModel, monitors = parameters.to.save)

CmyMCMC <- compileNimble(myMCMC)

results <- runMCMC(CmyMCMC, niter = n.iter, nburnin = n.burnin, nchains =n.chains, thin=n.thin)

calculateWAIC(CmyMCMC)

Perry de Valpine

unread,

Jul 27, 2022, 12:41:09 PM7/27/22

to Sarah Durham, nimble-users

Hi Sarah,

Yes, I think doing it post-hoc is reasonable here. Each capture history will be treated as something to predict. The main advantage of the online mode of WAIC is performance: faster speed and reduced memory use (and potentially not monitoring all nodes that would be needed for post-hoc WAIC). If those issues aren't obstacles, there is no reason not to do it post-hoc when conditional WAIC without groupings are what you want.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/92d0e0f4-0e8f-4e99-97c1-ae8a64ebb6c5n%40googlegroups.com.

Sarah Durham

unread,

Jul 27, 2022, 1:04:48 PM7/27/22

to nimble-users

Perfect, thanks Perry!

Reply all

Reply to author

Forward