### R-INLA code to run the analyses outlined In Perkins et al. (2020) Monitoring the resilience of a ### no take marine reserve to a range extending species using benthic imagery. PLOS One. ### Please cite this paper when using or modifying the below code ### For any questions, please contact Nick Perkins: Nicholas.Perkins@utas.edu.au #Set working directory to where data is stored - change this to your local directory setwd("XXXXXXXXXXXXX") rm(list = ls()) #clear workspace library(INLA) #load the INLA library. Please ensure this package is installed. #Read in data file. Note that this data file is available in the same repository as the code #The data has been subsetted to every fifth image to avoid overlapping images B <- read.csv("All_barrens_fifthimg.csv", stringsAsFactors = FALSE) #Remove data where there are NAs for rugosity (vrm_3). These were areas where there was no underlying #bathymetric mapping. B <- B[-which(is.na(B$vrm_3)), ] #Re-scale year column B$year <- B$year - 2010 #year is ordered so that 2011 is year '1' #Rename columns to match the terms used in the paper #Rename "actual depth" to "depth" to match terms used in paper names(B)[5] <- "depth" #Rename vrm_3 to rugosity to match terms used in paper names(B)[11] <- "rugosity" #Rename "MPA" to "NTR" to match terms used in paper names(B)[10] <- "NTR" #Logit transform rugosity B$rugosity <- log(B$rugosity/(1-B$rugosity)) #Scale covariates B$depth <- scale(B$depth) B$depth2 <- (B$depth)^2 B$rugosity <- scale(B$rugosity) ###INLA NON-SPATIAL MODEL### #This is model 1 described in the paper. There are no spatial or temporal random effects inla.formula <- barrens ~ NTR*year + rugosity + depth + depth2 mod.inla.1 <- inla(inla.formula, family = "binomial", data = B, control.fixed = list (prec = 0.01)) summary(mod.inla.1) ###SPATIAL MODEL### #This is model 2 in the paper. Spatial random effects are introduced. #Create a 'mesh' for Governor Island and controls #Non-convex hull mesh bnd <- inla.nonconvex.hull(cbind(B$utm.x, B$utm.y), convex = -0.01, resolution = c(24, 105)) mesh.b <- inla.mesh.2d(boundary = bnd, loc = cbind(B$utm.x, B$utm.y), max.edge = c(100,2000), cutoff = 5) plot(mesh.b) #visualisation of the mesh #Parameters for modelling #Spatial parameters specified using pc priors. For more details see Appendix B in the paper. spde <- inla.spde2.pcmatern( mesh=mesh.b, alpha=2, ### mesh and smoothness parameter -- alpha=3/2 is exponential prior.range=c(10, 0.1), ### P(practic.range<10)=0.1 prior.sigma=c(1, 0.02)) ### P(sigma>1)=0.02 #Prior for the temporal autoregressive parameter, see also spatio_temporal_prior_specs a.phi <- exp(0.95) b.phi <- 0.45 b.phi.prec <- 1/b.phi^2 h.spec <- list(theta = list(prior = 'normal', param = c(a.phi, b.phi.prec))) #Create index #'i' defines the index for the vertices of the mesh for calculating the spatial random effect #'that is defined by the spde model i.index <- inla.spde.make.index(name = "i", n.spde = spde$n.spde) #Create a sparse weight matrix 'A' by identifying the data locations in the mesh and organising the corresponding values of the basis functions A <- inla.spde.make.A(mesh = mesh.b, loc = cbind( B$utm.x, B$utm.y)) #Define the model matrix form <- ~ -1 + NTR*year + rugosity + depth + depth2 #formula for fixed effects - used for both spatial and spatio-temporal models X <- as.data.frame(model.matrix(form, data = B)) names(X)[6] <- "NTR_year" N <- nrow(B) #Create data stack stack.1 <- inla.stack(data=list(z=B$barrens, ntrial=B$ntrial), A=list(1, 1, A), effects=list( intercept = rep(1, N), X = X, i = i.index )) #Spatial model formula formula.inla.sp <- z ~ -1 + intercept + NTR*year + rugosity + depth + depth2 + f(i, model = spde) #Fit the spatial model inla.sp.res <- inla(formula.inla.sp, data=inla.stack.data(stack.1), control.predictor=list(compute = TRUE, A=inla.stack.A(stack.1)), control.inla=list(int.strategy='eb'), control.compute=list(dic = TRUE, cpo = TRUE), control.fixed = list(expand.factor.strategy='model.matrix', prec = 0.01), family = "binomial", Ntrials = ntrial) summary(inla.sp.res) ###SPATIO-TEMPORAL MODEL### #This is the full model (Model 3 in the paper) #Create a one dimensional 'mesh' for time mesh.t <- inla.mesh.1d(unique(B$year)) #Create a new index that also includes the temporal mesh #'i' defines the index for the vertices of the mesh for calculating the spatial random effect that is defined by the spde model #n.group is the number of time points in the data i.index2 <- inla.spde.make.index(name = "i", n.spde = spde$n.spde, n.group = length(unique(B$year))) #Create a sparse weight matrix 'A' by identifying the data locations in the mesh and organising the corresponding values of the basis functions #Also include the temporal projection by specifying the time points and temporal mesh A2 <- inla.spde.make.A(mesh = mesh.b, loc = cbind( B$utm.x, B$utm.y), group = B$year, group.mesh = mesh.t, n.group = length(unique(B$year))) #Combine data into an inla.stack object stack.2 <- inla.stack(data=list(z=B$barrens, ntrial=B$ntrial), A=list(1, 1, A2), effects=list( intercept = rep(1, N), X = X, i = i.index2 )) #Model specification with main effect for year and interaction effect for NTR and year, a main effect for VRM (rugosity measure) #and with a smoothed effect for depth. Default intercpet removed and renamed as "intercept". Year fitted as an ar1 process inla.spt.form <- z ~ -1 + intercept + NTR*year + rugosity + depth + depth2 + f(i, model = spde, group = i.group, control.group = list(model = 'ar1', hyper = h.spec)) #Fit spatio-temporal model - this may take some time!! inla.spt.res <- inla(inla.spt.form, data=inla.stack.data(stack.2), control.predictor=list(compute = TRUE, A=inla.stack.A(stack.2)), control.inla=list(int.strategy='eb'), control.compute=list(dic = TRUE, cpo = TRUE, config = TRUE), control.fixed = list(expand.factor.strategy='model.matrix', prec = 0.01), family = "binomial", Ntrials = ntrial) #Examine output summary(inla.spt.res) plot(inla.spt.res) #plot posterior summaries ###CODE TO TAKE POSTERIOR SAMPLE DRAWS AND PLOT FIGURES 3 AND 4 IN THE PAPER### #Create list of coefficients to retain to.keep <- list(1,1,1,1,1,1,1) names(to.keep) <- c("intercept", "MPA", "year", "rugosity", "depth", "depth2", "MPA:year") #Take 5000 posterior samples sim <- inla.posterior.sample(5000, result = inla.spt.res, selection = to.keep) #5000 posterior samples #Calculate probabilities over depths 2-40 metres sim.mpa.yr <- 6 #year 2016 sim.yr <- 6 #year 2016 in the MPA sim.rug <- mean(B$rugosity[which(B$MPA == 1)]) #mean value for rugosity in the MPA sim.depth <- scale(seq(2, 40, by = 1)) sim.depth2 <- (sim.depth)^2 sim.res.depth <- sapply(1:39, function (x) sapply(sim, function (s) s$latent[ "intercept:1",] + s$latent["MPA:1", ]*1 + s$latent[ "year:1",]*sim.yr + s$latent[ "rugosity:1",]*sim.rug + s$latent[ "depth:1",]*sim.depth[x] + s$latent[ "depth2:1",]*sim.depth2[x] + s$latent[ "MPA:year:1",]*sim.mpa.yr )) #Calculate means over 5000 samples sim.res.final.depth.means <- colMeans(inv.logit(sim.res.depth)) #Calculate credible intervals lo.95.depth <- apply(inv.logit(sim.res.depth), 2, quantile, 0.025) up.95.depth <- apply(inv.logit(sim.res.depth), 2, quantile, 0.975) #Bind ci's together for plotting sim.res.final.depth.cis <- rbind(lo.95.depth, up.95.depth) #Calculate probabilities over range of rugosity values from 0 to 0.75 sim.rug <- scale(seq(0, 0.75, by = 0.05)) sim.depth.sec <- mean(B$depth[which(B$MPA == 1)]) #mean value for depth in the MPA sim.depth2.sec <- mean(B$depth2[which(B$MPA == 1)]) #mean value for depth2 in the MPA sim.res.rug <- sapply(1:16, function (x) sapply(sim, function (s) s$latent[ "intercept:1",] + s$latent["MPA:1", ]*1 + s$latent[ "year:1",]*sim.yr + s$latent[ "rugosity:1",]*sim.rug[x] + s$latent[ "depth:1",]*sim.depth.sec + s$latent[ "depth2:1",]*sim.depth2.sec + s$latent[ "MPA:year:1",]*sim.mpa.yr )) #Calculate means over 5000 samples sim.res.final.rug.means <- inv.logit(colMeans(sim.res.rug)) #Calculate credible intervals lo.95.rug <- apply(inv.logit(sim.res.rug), 2, quantile, 0.025) up.95.rug <- apply(inv.logit(sim.res.rug), 2, quantile, 0.975) #Bind ci's together for plotting sim.res.final.rug.cis <- rbind(lo.95.rug, up.95.rug) #Examine prob. barrens through time for both MPA and reference sim.yr.2 <- sort(unique(B$year)) sim.mpa.yr.2 <- sort(unique(B$year)) sim.rug.mpa <- mean(B$rugosity[which(B$MPA == 1)]) sim.rug.ref <- mean(B$rugosity[which(B$MPA == 0)]) sim.depth.mpa <- mean(B$depth[which(B$MPA == 1)]) sim.depth2.mpa <- mean(B$depth2[which(B$MPA == 1)]) sim.depth.ref <- mean(B$depth[which(B$MPA == 0)]) sim.depth2.ref <- mean(B$depth2[which(B$MPA == 0)]) #MPA sim.res.mpa <- sapply(1:4, function (x) sapply(sim, function (s) s$latent[ "intercept:1",] + s$latent["MPA:1", ]*1 + s$latent[ "year:1",]*sim.yr.2[x] + s$latent[ "rugosity:1",]*sim.rug.mpa + s$latent[ "depth:1",]*sim.depth.mpa + s$latent[ "depth2:1",]*sim.depth2.mpa + s$latent[ "MPA:year:1",]*sim.mpa.yr.2[x] )) sim.res.final.mpa.means <- inv.logit(colMeans(sim.res.mpa)) lo.95.mpa <- apply(inv.logit(sim.res.mpa), 2, quantile, 0.025) up.95.mpa <- apply(inv.logit(sim.res.mpa), 2, quantile, 0.975) sim.res.final.mpa.cis <- rbind(lo.95.mpa, up.95.mpa) #REF sim.res.ref <- sapply(1:4, function (x) sapply(sim, function (s) s$latent[ "intercept:1",] + s$latent["MPA:1", ]*0 + s$latent[ "year:1",]*sim.yr[x] + s$latent[ "rugosity:1",]*sim.rug.ref + s$latent[ "depth:1",]*sim.depth.ref + s$latent[ "depth2:1",]*sim.depth2.ref)) sim.res.final.ref.means <- inv.logit(colMeans(sim.res.ref)) lo.95.ref <- apply(inv.logit(sim.res.ref), 2, quantile, 0.025) up.95.ref <- apply(inv.logit(sim.res.ref), 2, quantile, 0.975) sim.res.final.ref.cis <- rbind(lo.95.ref, up.95.ref) #Plots #The below code will produce the plots Fig. 3 and Fig. 4 in the paper #Plots will be saved in the working directory png("DepthRugProbs.png", height = 3200, width = 2500, res = 300) par(mfrow = c(2,1)) par(mar = c(6.1, 5.1, 4.1, 2.1)) #Depth depth.orig <- t(apply(sim.depth, 1, function(r)r*attr(sim.depth,'scaled:scale') + attr(sim.depth, 'scaled:center'))) plot(depth.orig, sim.res.final.depth.means, type = 'l', xlab = "Depth", ylab = "Probability of barrens presence", cex.lab = 1.4, cex.axis = 1.1, ylim = c(min(sim.res.final.depth.cis[1,]), max(sim.res.final.depth.cis[2,])) ) polygon(c(depth.orig, rev(depth.orig)) ,c(sim.res.final.depth.cis[1,], rev(sim.res.final.depth.cis[2,])), col = adjustcolor("grey", alpha.f = 0.3), border=NA) mtext("A", side = 3, line = 0.5, adj = 0, cex = 3) #Rugosity rug.orig <- t(apply(sim.rug, 1, function(r)r*attr(sim.rug,'scaled:scale') + attr(sim.rug, 'scaled:center'))) plot(rug.orig, sim.res.final.rug.means, type = 'l', xlab = "Rugosity", ylab = "Probability of barrens presence", cex.lab = 1.4, cex.axis = 1.1, ylim = c(min(sim.res.final.rug.cis[1,]), max(sim.res.final.rug.cis[2,])) ) polygon(c(rug.orig, rev(rug.orig)) ,c(sim.res.final.rug.cis[1,], rev(sim.res.final.rug.cis[2,])), col = adjustcolor("grey", alpha.f = 0.3), border=NA) mtext("B", side = 3, line = 0.5, adj = 0, cex = 3) dev.off() #MPA/REF through time plot png("MPA_REF_time.png", height = 3000, width = 4000, res = 300) par(mar = c(5.1, 4.9, 4.1, 2.1)) yrs <- c(2011, 2013, 2014, 2016) plot(yrs, sim.res.final.ref.means, type = 'l', xlab = "Year", ylim = c(1e-06, to = 0.02), ylab = "Probability of barrens presence", lty = 1, log = 'y', cex.lab = 2.2, cex.axis = 1.6) polygon(c(yrs, rev(yrs)) ,c(sim.res.final.ref.cis[1,], rev(sim.res.final.ref.cis[2,])), col = adjustcolor("grey", alpha.f = 0.3), border=NA) lines(yrs, sim.res.final.mpa.means, ylim = c(1e-06, to = 0.02), lty = 2) polygon(c(yrs, rev(yrs)) ,c(sim.res.final.mpa.cis[1,], rev(sim.res.final.mpa.cis[2,])), col = adjustcolor("grey", alpha.f = 0.3), border=NA) legend("topleft", legend = c("outside NTR", "NTR"), lty = c(1,2), bty = 'n', cex = 2) dev.off()