specifying different data for first iteration only?
Rebecca Taylor
Hi all,
I am working with an age-structured Jolly-Seber model (superpopulation
parameterization with data augmentation.) Ages are observed without error for
real individuals. The time 1 population is comprised of individuals of many ages.
Entries in time 2 and beyond are by birth only. For augmented individuals, I
would like to initialize all ages at “NA” so the ages can simply be determined by
the code.
That works fine if survival and recapture rates are not age-specific. However, if rates are age-specific, then the initialization fails because there is an “NA” being used as an index. For example, any individual in time t has a survival rate of phi[t, NA].
I am able to work around this by assigning all augmented individuals an age in time 1. But that means I need many groups of augmented individuals—e.g., if there are 10 ages, there need to be 11 groups of augmented individuals (the extra group is for not yet born). And all groups of augmented individuals have to be pretty large, so I end up with a much larger augmented population and run time than I would like.
Is there a way I can do something different on the first iteration to get around the initialization problem? Perhaps assign everyone a non-age-dependent survival rate on the first iteration only? Or assign everyone an initial age in the first time period for the first iteration, but then let those initial ages be overwritten in subsequent iterations?
I hope this makes sense. I’ll copy code below. Thank you!!
Rebecca
code <- nimbleCode({
psi ~ dbeta(1,1) # psi = inclusion parameter (is individual real or fake)
for (t in 1:K){
for(aa in 1:max.age){ # age 0 means not yet born
p[t,aa+1] ~ dbeta(1,1)# p[t,aa+1] = cap rate for time t, age aa; age 0 cap rate = 0
}
}
for (t in 1:(K-1)){
for(aa in 1:(max.age-1)){
phi[t,aa+1] ~ dbeta(1,1) # phi[t,aa+1] = surv rate for time t, age aa; age 0 surv rate = 0
}
}
beta[1:K] ~ ddirch(b[1:K]) # betas = unconditional entry probabilities
eta[1] <- beta[1] # etas = corresponding conditional probabilities
for(k in 2:K){eta[k] <- beta[k]/(1-sum(beta[1:(k-1)]))}
# piAGE = age distribution of animals alive in the population at time 1
# piAGEuncond has an additional "age 0" which means "not yet born".
# Entries after time 1 are by birth only
piAGE[1:max.age] ~ ddirch(a[1:max.age])
piAGEuncond[1:(max.age+1)] <- c( (1-eta[1]), eta[1]*piAGE[1:max.age] )
# Likelihoods
for (i in 1:M){ # M = size of augmented population
w[i] ~ dbern(psi)# w[i] = indicator if individual i is real
# agePlusOne are time 1 age data + 1
agePlusOne[i] ~ dcat(piAGEuncond[1:(max.age+1)])
age[i,1] <- (agePlusOne[i]-1)
# state process
u[i,1] <- step(age[i,1]-.1) # u[i,j] = indicator for alive
z[i,1] <- u[i,1]*w[i]
# Observation process
yprob_t1[i,1] <- z[i,1]*p[1,age[i,1]+1]
y[i,1] ~ dbern(yprob_t1[i,1] )
# avail[i,j] = indicator if not yet recruited from superpop to pop
avail[i,1] <- 1 - u[i,1]
for (t in 2:K){
# State process
uprob[i,t] <- u[i,t-1]*phi[t-1,age[i,t-1]+1] + avail[i,t-1]*eta[t]
u[i,t] ~ dbern(uprob[i,t])
z[i,t] <- u[i,t]*w[i]
# Age process: grows older by one year after recruitment
age[i,t] <- age[i,t-1] + max(u[i,1:t])
# Observation process
yprob[i,t] <- z[i,t]*p[t,age[i,t]+1]
y[i,t] ~ dbern(yprob[i,t])
# avail[i,j] = indicator if not yet recruited from superpop to pop
avail[i,t] <- 1- max(u[i,1:t])
} #t
} #i
})#END model
Daniel Gibson
# State Process
######################################
for (i in 1:M){ # Loop through real + psuedo individuals
w[i] ~ dbern(psi) # Latent inclusion state
z[i,1] ~ dbern(eta[1]) # Entry into first occasion
for (t in 2:n.occasions) {
q[i,t-1] <- 1 - z[i,t-1] # Have not entered yet
mu1[i,t] <- phi[i,t-1] * z[i,t-1] + eta[t] * prod(q[i,1:(t-1)]) # Presence = survived since last + entered since last
z[i,t] ~ dbern(mu1[i,t])
}
}
######################################
# Observation Process
######################################
for (i in 1:n.ch){ # Loop through real individuals
y[i,1] ~ dbern(z[i,1] * p[1] * w[i]) # Latent state * detection * latent inclusion
for (t in 2:n.occasions){
y[i,t] ~ dbern(z[i,t] * p[t] * w[i]) # Latent state * detection * latent inclusion
}
}
######################################
# Chandler trick
######################################
for(i in (n.ch+1):M) { # Loop through psuedo individuals
pzo[i,1] <- 1 - (p[1] * z[i,1] * w[i]) # probability of not being detected on first occasion
for (t in 2:n.occasions){
pzo[i,t] <- 1 - (p[t] * z[i,t] * w[i]) # probability of not being detected on each other occasion
}
zeros[i] ~ dbern(1-prod(pzo[i,])) # probability of at being detected at least once (zeros = vector of 0's for each augmented individual)
}
--
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/3c0504a2-ceb1-4ccd-85d8-255642b7e38en%40googlegroups.com.
Taylor, Rebecca L
Sent: Saturday, October 15, 2022 3:37 PM
To: Taylor, Rebecca L <rebecc...@usgs.gov>
Cc: nimble-users <nimble...@googlegroups.com>
Subject: [EXTERNAL] Re: specifying different data for first iteration only?
|
This email has been received from outside of DOI - Use caution before clicking on links, opening attachments, or responding. |
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nimble-users/YFi4aBAQkcA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAD5tc9xK3LFa07SPcz77%2BFsXCU_4N9S8o6SwOE2cVrL_%2Bwe-FQ%40mail.gmail.com.
Perry de Valpine
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/BY5PR09MB5154E83FB631C9E921763558D4269%40BY5PR09MB5154.namprd09.prod.outlook.com.
Taylor, Rebecca L
Hi Perry,
Thank you for your response and for all you do to help us.
The error was indeed on initialization, and the first line to cause the error was
yprob_t1[i,1] <- z[i,1]*p[1,age[i,1]+1] because age[i,1]equaled “NA” for the augmented individuals, making the index evaluate to [1,NA+1]. At the time I posted my question, the code only worked if time 1 age data were specified numerically for all individuals, including the augmented individuals.
I previously tried providing numerical inits for all of the entries where I wanted to have “NA” data values, but it did not work. However, since you suggested it again, I tried it again, and this time it did work—which means something else was wrong the last time I tried it. I am taking the next online course y’all are offering, so hopefully that will help me become less reliant on this forum.
I see what you are saying about inefficient sampling for this kind of data augmentation—thank you for taking the extra time to mention this! On my to-do list is combining w with the betas, which I think will alleviate the particular problem you are referencing (although please correct me if I am missing something). Next is marginalizing over individuals so I am looping through unique capture histories, not unique individuals. After that is deciding if I need to abandon data augmentation. As you can tell, however, I am taking baby steps….but if there are any written resources you can direct me to for making this more efficient, that would be outstanding.
Thank you again.
Rebecca
Sent: Wednesday, October 19, 2022 7:55 AM
Subject: Re: [EXTERNAL] Re: specifying different data for first iteration only?