Hello,
I’m interested in fitting some amalgamation of a multi-scale/staggered entry/dynamic occupancy models as a hidden Markov Model (following a previous suggestion in the thread here). For the time being at least, the general structure I’m pursuing has a complete data likelihood that looks a little something like this, where i indexes a site, t indexes a temporal occasions, and c indexes a replicate detectors operating during i and t:
z_i~Bernoulli(\psi_i)
logit(\psi_i) = s(space) +\bf\beta X_i
a_{it}~Bernoulli(\theta_{it})
logit(\theta_{it}) = \bf\delta s(time) X_i+\bf\lambda X_{it}+\gamma a_{it-1}
y_{itc}~Bernoulli(z_ia_{it}p_c)
logit(p_c) = \bf\alpha X_c
So, the general idea is that y can = 1 if z_i and a_it =1, and that a_it is a first-order (hidden) Markov model. The smoothers are implemented in the name of efficiency/dimension reduction (although tbh, the program as I’ve written is slow almost certainly because I’ve just hacked away at getting it running). To the extent that I’ve been able to interrogate the problem, treating availability for detection as a first order markov process seems to be more computationally efficient (perhaps because it constrains the state transitions) than a model in which transitions in z_it are of interest.
At any rate, the marginal approach to HMM and the forward algorithm is an exciting but very new way for me to think about these problems. While there are some excellent resources on HMM and dynamic in the manual, in the BPA translations on github, and vignettes focusing on HMM for animal movement models, I’m a little curious about a few issues/questions I haven’t seen fully addressed elsewhere. My interpretation of the transition and emission sort of matrices looks a little like so:
for (c in 1:ndetectors){
//state, detector # , observation
po[1, c, 1]=1;
po[1, c, 2]=0;
po[2, c, 1]=1;
po[2, c, 2]=0;
po[3, c, 1]=1-p[c];
po[3, c, 2]=p[c];
}
for (i in 1:ncells){
// ps1[i, 1]=1-psi[i];
//ps1[i, 1]=psi[i]*(1-inv_logit(logittheta[i, 1]));
//ps1[3, i]=psi[i]*inv_logit(logittheta[i, 1]);
for (t in 1:364){
ps[1, i, t, 1]=1; //unoccupied remains unoccupied
ps[1, i, t, 2]=0;
ps[1, i, t, 3]=0;
ps[2, i, t, 1]=0;
ps[2, i, t, 2]=1-inv_logit(logittheta[i, t+1]); //occupied
ps[2, i, t, 3]=inv_logit(logittheta[i, t+1]);
ps[3, i, t, 1]=0;
ps[3, i, t, 2]=1-inv_logit(logittheta[i, t+1]+a_AL);
ps[3, i, t, 3]=inv_logit(logittheta[i, t+1]+a_AL);
}
}
}
instead of
model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}
