Code to accompany Allen et al, among-site variability in the stochastic dynamics of East African coral reefs.

The following instructions assume that you are using Linux (tested on a system with 4 3.2 GHz Intel Xeon cores, 16 GiB RAM, 64-bit Ubuntu 12.04, Stan 2.7.0, R 3.2.3). They may also work in OSX.

1. Create the required directory structure. The R script and shell script included in the zip file assume this structure exists when opening data files and producing output. The Stan code will of course work fine without setting up this structure.

mkdir coral
mkdir coral/stanresults
mkdir coral/project
mv distribution.zip coral/project
cd coral/project
mkdir figures
mkdir figures/obs_truestates_sep
unzip distribution.zip

After unzipping, coral/project should contain the following subdirectories:

Data
STAN
R

2. Install Stan, following the instructions at http://mc-stan.org/
You will need both the command-line interface (cmdstan) and the rstan R package. We used cmdstan to run the Monte Carlo simulation, because unlike rstan, it does not hold the whole of the simulation output in memory. We used rstan only for producing input files in stan-readable form, and extracting the results.

3. Install any required R packages. As well as rstan, you will need coda, RColorBrewer, grDevices, compositions, MASS, numDeriv, loo, mvoutlier, MVN, moments, mvtnorm, latex2exp, ncf, hdrcde, and the packages on which they depend.

4. Compile the Stan code. 

In a terminal, cd to the directory containing your cmdstan installation and type

make [path to STAN directory]/coralss_t

5. Run the Stan code. 

In a terminal, cd to the STAN directory
num_samples and num_warmup in run4chains_t_estimatedf.sh determine the size of the sample from the posterior distribution and the length of the warmup period before sampling. With these variables set to the values used in the manuscript (both 5000), sampling will take about 2 hours on 4 Intel Xeon 3.2GHz cores, with 16GiB RAM. If you want to test that everything is working first, reduce these numbers to something small, such as 50. Then type

./run4chains_t_estimatedf.sh

At the start of the run, you will see messages like this:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception thrown at line 47:
stan::math::multi_student_t: LDLT_Factor of scale parameter is not positive definite. last conditional variance is -nan.
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

These can be ignored. Note that after the last chain has finished, the shell script doesn't automatically return you to the command prompt (pressing enter brings up the prompt).

6. Process the results using R.

cd to the R directory
open the script coralsummary.R in an editor 
set which outputs you want using the logical variables on lines 23-56. The following options (which are the current settings) will produce a basic set of figures similar to those in the main text of the manuscript:

makestan <- FALSE #create a .Rdump file (name dumpname) from a data frame (ddatas) including reefid, year, transect and yt (two transformed composition components)
uset <- TRUE #t rather than normal errors
savesmall <- TRUE #save a smaller .RData file containing everything extracted from standata
loadfull <- TRUE #if FALSE, load from smaller .RData file instead of full csv files: doesn't work with printsummary, dotraceplots, savesmall, dopostmeans
dotraceplots <- FALSE #MCMC trace plots
dolowcoral <- TRUE #calculate probability of low coral cover
printsummary <- TRUE #summary of MCMC diagnostics
doeigen <- FALSE #stuff related to eigenvalues and volume ratio
doreefmeans <- TRUE #calculate and plot reef means
dopostmeans <- TRUE #calculate and plot posterior means of true states
doresids <- FALSE #residuals: calculate and plot
docovplots <- TRUE #covariance matrix ellipses
dosst <- FALSE #SST calculations
doparms <- TRUE #plot parameters
docorrs <- FALSE #correlations between single true state and stationary mean
plotlowcoral <- TRUE #plot derivatives of probability of low coral cover, and probability of low coral against observed mean coral cover
makepairsderiv <- FALSE #pairwise scatter plots of derivatives
checkderivs <- FALSE #check whether derivatives satisfy conditions on B, Sigma, Zeta
makeanalyticalderivs <- FALSE #save analytical derivatives of prob low coral cover wrt B
plotanalyticalderivs <- FALSE #plot and analyze analytical derivatives of prob low coral cover wrt B (loaded from file)
savereefeffects <- FALSE #save reef effects and lp__
doparamcheck <- FALSE #check derivatives are sane by making moderately large changes in parameters
doepsilon <- FALSE #save epsilons from a few representative iterations
diagnostics <- FALSE #diagnostic plots for random parts
simulatedata <- FALSE #generate simulated data sets using posterior mean estimates from real data: only implemented for t version
comparemodels <- FALSE #comparison of normal and t models based on leave-one-out cross-validation and WAIC
checksimulated <- FALSE #compare posterior distributions of parameters to true values from simulated data
checkmanysimulated <- FALSE #compare posterior distributions to true values in multiple simulated data sets
readdumpfiles <- FALSE #read all the raw dump files from multiple simulated data sets (if false, load from file)
plotderivintuition <- FALSE #plot to get intuition about derivatives of mustar wrt B
checkPlow <- FALSE #check that probability of low coral calculations are correct, by simulation
plotilrintuition <- FALSE #plot to get intuition about ilr transformation
plowcoralrandom <- FALSE #report P(coral <= 0.1) for a randomly-chosen reef
plotsplinecorrel <- FALSE #plot spline correlograms for P(coral <= 0.1)
docountpatch <- FALSE #count number of patch and fringing transects
plotalphavsamplesize <- FALSE #plot generalized variance of alpha against sample size
simulatelowcoralvsmean <- FALSE #simulate relationship between sample mean coral cover and probability of coral cover <= 0.1

Once the options are set, start R and in the R console type

source("coralsummary.R")

The script will take a few minutes to run.

7. Inspect the outputs. In the figures directory:

The directory obs_truestates_sep contains time series plots for each site, with posterior means from the model, as in Figure 1. For example, Mombasa1_tback.eps is for site Mombasa1.

The file parmst.eps should be similar to Figure 2 (although the dots representing the distribution of parameter estimates will be very faint if the size of the Monte Carlo sample was small).

The file amongreef_t.eps should be similar to Figure 3.

The file covellipsesback_t.eps should be similar to Figure 4.

The file problowcoralobs_0.1t.eps should be similar to Figure 5.

The file derivs_0.1t.eps should be similar to Figure 6.
