State-space models are a special type of hierarchical model designed for dealing with time series of counts or other data. With state-space models, we decompose the time series into 2 components:
A process component that describes how the actual state of the system (e.g., population size, age distribution, spatial location) varies through time, usually as a function of initial conditions (e.g., initial abundance) and transition parameters (such as survival, reproduction, and movement rates).
An observation component that relates the system states and transition parameters to field observations. For example, Nt might be abundance at time t and Ct sample counts of animals.
Kéry and Schaub provide a series of examples in Chapter 5 to illustrate state-space models. We will take their real data example on House Martins in Switzerland, using the implementation in the JAGS language via jagUI. I re-formatted Kéry and Schaub's code slightly to make it clearer that N is the system state (unknown) and hm.counts are the vector of count data. The data are transformed to the log scale in R outside the JAGS model and read in on the log scale. In the model, the state space transition is formulated on the log scale via
logN[t+1] <- logN[t] + r[t]
which is equivalent to the growth model
N[t+1]<-N[t]exp(r[t]).
Population growth is allowed to be a random effect by
r[t] ~ dnorm(mean.r, tau.proc).
Finally, the log-scale counts are related to the abundance state by
y[t] ~ dnorm(logN[t], tau.obs)
where
y[t]<-log(C[t]).
This formulation allows for separate estimates of the variation due to the process model (tau.proc) from that due to sampling error (tau.obs), but otherwise assumes that the counts provide unbiased inference (E[y[t]]= logN[t]). Priors must be placed on initial abundance logN[1], mean.r, and tau.proc.
logN[1] ~ dnorm(5.6, 0.01) # Prior for initial population size
mean.r ~ dnorm(0, 0.001) # Prior for mean growth rate
sigma.proc ~ dunif(0, 10) # Prior for sd of state process
sigma2.proc <- pow(sigma.proc, 2)
tau.proc <- pow(sigma.proc, -2)
sigma.obs ~ dunif(0, 10) # Prior for sd of observation process
sigma2.obs <- pow(sigma.obs, 2)
tau.obs <- pow(sigma.obs, -2)
The code produces posterior inference on r, N, tau.proc and tau.obs and plots of the observed and "estimated" (predicted from the state-space model) and abundance. Note that predictions go base the data; here we are basically extending the model to predicting the future population trajectory assuming that the parameters estimated from the data hold remain stationary. You should not be surprised to see that the 95% credibility intervals (gray area of graph) becomes very wide as we get farther into the future.
You should also look over the examples in Chapter 5 that involve more general assumptions about the observation process. For example, in section 5.3 there is an example involving systematic bias (underestimation) in counts, where counts are assumed to come from binomial process. In this example, the counts provide a reliable index (so proportional to N) but are incapable of providing unbiased estimates of N unless additional information on the detection process is available. This requires either prior information on the detection process (e.g., an informative prior from another study) or combining data structures, something considered in the next lab.
Integrated Population Models (Chapter 11 Kéry and Schaub 2011 )
Often we are faced with a situation where we have several different kinds of data that relate to a population process, and we would like to integrate them in a coherent framework. For instance, it is fairly typical to have population counts from one source, capture-recapture data for survival estimates, and estimates or indices of reproductive success from another source (such as a nest survey, or harvest data). We could get independent estimates of abundance (N), survival (S), and reproductive rates (f), and then piece these together into a demographic model, but that approach is inefficient (and sometimes difficult because of the need to estimate variance for the population predictions) compared to integrated population models (IPMs), an extension of the state-space idea (Chapter 5) that we've already considered.
We can illustrate the IPM approach with the simulated data set in Section 11.3, which I've implemented in JAGS rather than BUGS (although the syntax is almost identical). The steps involved to build the IPM for this data are pretty straightforward and can be summarized as follows:
Build a state-space model for change in N through time. For this example the model can be constructed as a very simple matrix model (see page 351) involving numbers of juveniles (N1) and adults (Nad) at each year after the first year. Because the model provides no information for year 1, N1[1] and Nad[1] must be initialized using prior distributions.
Specify a likelihood for the population counts in each year.
Note that the counts are simply total numbers of animals observed in year t, which in turn is implicitly the sum of surviving adults and juveniles from the previous year. Because these values are not separately observable (i.e., we cannot distinguish 1-year-old from older birds in the counts) N1 and Nad are known as latent variables (or latent states).
Also note that the likelihood in this example assumes that the counts are distributed as Poisson variables with mean N1 +Nad , i.e., they are unbiased for estimating total abundance. To relax this assumption we would need additional data in order to estimate detection probabilities.
Specify a likelihood for the capture-recapture data. Since these are conventional CMR data, we will use CJS multinomial likelihoods previously introduced. This will give posterior values for survival St, which is in turn used in the state-space model to project abundance growth through time.
Specify a likelihood for the reproductive success data. In this example, Rt broods are surveyed each year and Jt nestlings counted. We assume that the broods are a random sample from the population, and model Jt using a Poisson distribution with mean Rtft, where ft is the productivity parameter that is in turn used in the state-space model to project abundance growth through time.
The IPM is implemented in R using jagsUI. The code produces posterior values for abundance, survival, and reproductive rates, as well as a graph comparing population counts with the posterior predictions for total abundance (K & S call these 'estimates', I think 'predictions' is a more accurate term). I have added code to compute posterior predictions for lambda[i] from the saved traces (sim.list) for Ntot[i], i=1,..10. As mentioned by Kéry and Schaub, computation of lambda from within the JAGS code sometimes fails because of weird initial values for Ntot[i] (e.g., zero) and this approach essentially moves that calculation outside of JAGS after convergence.