with \(y\left(s_{i},\ t\right)\ \) the categorisation of barren presence for an image located at \(s_{i}\)in year \(t\), space and time-varying covariate matrix \(X\left(s_{i},\ t\right)\) and separable spatio-temporal random effects \(z\left(s,t\right)\). The temporal covariance function is a stationary autoregressive process, such that \(-1<\ \phi<1\), where\(z\left(s,\ t\right)=\ \phi z\left(s,\ t-1\right)+\ \omega(s,\ t)\)with \(z(s,\ 1)\) drawn from the stationary distribution.
The spatial random effects in year \(t\) have a stationary spatial covariance function with spatial correlation function\(H(s-\ s^{{}^{\prime}},\rho)\) for a site \(s\) and \(s^{{}^{\prime}}\). The spatial range is defined by Lindgren et al. (2011) as the distance at which the spatial correlation drops to close to 0.1. The internal parameterisation of the range \(\rho\) and spatial variance \(\sigma^{2}\) used by INLA is given in Appendix A.
The model covariates for which β coefficients were estimated included fixed effects for the NTR (binary, whether an image was located inside the NTR or not), year, an interaction term between NTR and year, depth, depth-squared (to capture expected non-linear effects in deep and shallow images), and rugosity (see above). Depth was included as it has been previously found to be an important predictor of barrens presence (Johnson et al., 2005, Perkins et al., 2015).
Model M 2 is obtained fromM 3 by omitting the temporal dependence component (i.e., setting the temporal dependence parameter to zero). ModelM 1 is then obtained from modelM 2 by also omitting the spatial dependence component (i.e., setting the spatial variance parameter to zero). The hypothesised models are considered equally likely a priori and compared on the basis of the marginal likelihoods given the observed data. The marginal likelihood (or evidence) of each model is proportional to the posterior model probabilities given the a priori equal model weights.
We used the Integrated Nested Laplace Approximation approach (INLA; see Rue et al., 2009) for Bayesian spatial and spatio-temporal modelling. All statistical analyses were conducted within the R statistical computing package (R Core Team, 2019).
As the values of rugosity were right-skewed (i.e. mostly small values, with a few larger values), to avoid leverage issues a logit transformation was applied to the raw rugosity values. All physical model covariates (i.e. depth, depth-squared and logit-transformed rugosity) were centred by their mean and scaled by their standard deviation.
We examined both the percentage change in odds ratios and predicted changes in the probability of urchin barren presence to interpret the relationships with respect to model covariates. For the percentage change in odds ratios given by an increase of one unit in thei th covariate, we used the formula\(\operatorname{(exp}\left(\beta_{i}\right)-1)*100\). We examined the influence of covariate effects by plotting the expected probability of barrens presence over the range of values for that covariate in our sample space while holding all other covariates at their mean and excluding spatial random effects. That is, we plot the marginal relationships between the probability of barrens presence for both depth and rugosity. This was accomplished by taking 5000 joint posterior draws of the unknown β coefficients from the fitted model. We chose to examine the influence of each covariate within the NTR in 2016 (i.e. the last year surveyed).