Kevin Gross, Anthony R. Ives, and Erik V. Nordheim. 2005. Estimating fluctuating vital rates from time-series data: a case study of aphid biocontrol. Ecology 86:740–752.

Appendix E. Markov-chain Monte Carlo. See Appendix H for references cited.

Markov-chain Monte Carlo (MCMC) is a numerical technique used to estimate probability distributions (Metropolis et al. 1953, Hastings 1970; Geman and Geman 1984, Gilks et al. 1996). MCMC uses simple transition rules to define a Markov chain whose stationary distribution equals the target distribution to be estimated. In a Bayesian context, the target distribution is the joint posterior distribution of the model parameters. The Markov-chain produces sequences of parameter values, and these values can be aggregated and treated as a pseudo-random sample from the posterior distribution. 

The challenge of implementing MCMC is to design transition rules so that the chain moves (or mixes) through the posterior distribution efficiently. This ensures that the MCMC random sample obtained reflects the shape of the posterior, and not the movement tendencies of the chain algorithm itself. For this model, the primary obstacle to an efficient MCMC algorithm arose from the fact that any change in _t changed the values of all subsequent hidden state variables X_t+1, X_t+2, …, X_T. Consequently, values of _t closer to the beginning of a cutting cycle did not mix as freely as values closer to the end of the cutting cycle. We overcame this problem by designing an algorithm that proposed new values for early _t more frequently than for later _t. 

The updating algorithm is a single-component Metropolis updating scheme, meaning that parameters were updated individually, and all proposal distributions were symmetric (Metropolis et al. 1953). At each iteration of the Markov chain, we proposed updates for all parameters by looping through the parameter set according to the following steps:

Update

Update (5X)

Repeat for each cutting cycle {

Update initial total aphid density

Update initial juvenile fraction

Update initial juvenile parasitism

Update initial adult parasitism

Update initial total mummy density

Update initial fraction hyperparasitized

Loop forward from k=1 to k=n-1 {

Loop backward from j=k to j=1 {

Update

Update

Update

Update

Update

Update

}

Update X1 (as above)

}

}

Update kA (10X)

Update kM (10X)

Update s (10X)

Update tc (10X)

Update W (10X)

Update W (10X)

Update H (10X)

Update H (10X)

For parameters with prior support on the positive real line (k_A and k_M only), proposal distributions were normals folded at 0. For all other parameters, proposal distributions were doubly folded normals. All proposal distributions were centered at the current value of the parameter. For the parameters in  and (i.e., those that did not depend on cutting cycle), proposal distributions were selected to achieve rejection rates between 50% and 85% (Roberts et al. 1994, Gelman et al. 1996). For X₁'s and _t's, proposal distributions were chosen so that rejection rates were in the same range, although for lower rejection rates were occasionally tolerated (instead of using customized proposal distributions for each parameter). Variance parameters for the proposal distributions are given in table E1.

We ran 15 separate chains, and aggregated the output to form the MCMC sample. Each chain was run for 12,000 iterations, with the first 2000 iterations discarded and every 100th iteration saved thereafter.  Thus, the final MCMC sample consisted of 1500 sets of parameter values. Chains were initiated from parameter sets randomly selected from the prior, with the exception that the initial values of _j were the same for each _j and for each cutting cycle, and the initial values of X₁ were the same for each cutting cycle. Each chain took 1.5 – 2.5 hours to run, using the computing resource pool available at UW-Madison through condor(http://www.cs.wisc.edu/condor/). 

Convergence was determined by calculating estimated potential scale reduction factors () (Gelman 1996) for each of the 484 individual estimands in the posterior.  \epsr\ is the square root of the ratio of the estimated posterior variance from all chains pooled together to the estimated variance from the chains considered separately. As the separate chains converge to their stationary distribution,  approaches 1 from above; Gelman suggests running chains until all
's are less than “1.1 or 1.2”. Of the 484 's calculated here, the maximum was 1.02, indicating that the chains had converged to the posterior distribution. We verified that the burn-in was sufficiently long by comparing the average log posterior of the first 20 recorded iterations from each chain (12.127393) to the average log posterior from the final 50 recorded iterations for each chain (12.127386).

Table E1. Values of ² for proposal distributions.

Parameter	2	Parameter	2
	0.025		0.1
	0.02		0.5
Aphid density	0.25		1.5
Juvenile fraction	0.04	kA	2.5
Juvenile parasitism	0.75	kM	5
Adult parasitism	0.1	s	0.25
Mummy density	0.2	tc	0.5
Mummy hyperparasitism	0.5	W	0.025
	0.05	W	0.05
	0.05	H	0.05
	0.2	H	0.05