Code used in "Density-structured dynamic species distribution models for the top shells Osilinus lineatus and Gibbula umbilicalis"

Nova Mieszkowska
Gregg Milligan
Mike Burrows
Rob Freckleton
Matt Spencer

This directory contains the most important bits of code used in the paper.

Updated 2 June 2016: fixed the bug reported below, and should now work with recent versions of JAGS (4.2.) and rjags (4.6).

Tested with R version 2.11.1 on 32-bit linux, and with R version 2.14.1, JAGS 3.1.0 on 64-bit Ubuntu 12.04.2. We have had a report (thanks to Remi Vergnon) that with JAGS 3.3.0 on Mac OSX 10.7.4, running the example script gives the error message

> Error in setParameters(init.values[[i]], i) : Error in node alpha[1]
> Attempt to overwrite value of observed node

Using JAGS 3.1.0 instead apparently fixes this.


Requires the packages rjags, nnet (only for ML fitting, not used in the paper), mgcv, and lattice, plus their dependencies. Note that rjags depends on the jags program, which is not part of R, but is available at

http://mcmc-jags.sourceforge.net/

(definitely available for Linux and Mac. Not sure about Windows).

RUNNING THE EXAMPLE

In R,

source("example.script.R")

should do the following:

-generate a small simulated data set (by calling simulate.densitystructure.R)
-fit a baseline-category logit model by MCMC (the functions that do this are in chainfuncsmulti.R, and the model is specified in marclim.logistic.multi.bug). Note that in this test code, the MCMC is not run for long enough to converge. Also, because we simulate only a small sample, you might get errors further down connected with smoothing in the model diagnostics (if this happens, just run the simulation again).
-print out some MCMC diagnostics
-display matrices of observed and predicted transition counts (might reveal gross errors e.g. using the wrong data. In other cases, even a poorly specified model will give the right predicted counts)
-display some graphical model diagnostics using plotP.R, as in the supporting information for the paper. For each transition probability, plot the marginal relationship between the transition probability and an explanatory variable, with marginal dotplots of the indicator variable for occurrence of that transition. A binomial GAM smoother is used to approximate the marginal relationship between the observations and the explanatory variable, although it often oversmooths.
-plot contours of the relationships between probability of presence and the two explanatory variables, and normalized entropy and the two explanatory variables. This is done by a call to example.envelope.R, which first uses construct.P.R to calculate fitted transition probability matrices on a grid of values of the explanatory variables, and then uses stat.dist.R to calculate summary statistics, including probability of presence and entropy. Finally, envelope.plot.R is used to produce the graphs.

Matthew Spencer, 19/09/2012

m.spencer@liverpool.ac.uk
