Effects of marine reserves in the context of spatial and temporal variation : an analysis using Bayesian zero-inflated mixed models

Evaluating the effects of marine reserves on exploited species can be challenging because they occur within a context of natural spatial and temporal variation at many scales. For rigorous inferences to be made, such evaluations require monitoring programmes that are replicated at appropriate scales. We analysed monitoring data of snapper Pagrus auratus (Sparidae) in northeastern New Zealand, comprised of counts from baited-underwater-video surveys from inside and outside 3 marine reserves. Surveys were replicated at many levels, including areas inside and outside of marine reserves at 3 locations in 2 seasons, over a period of up to 14 yr, in an unbalanced design. The Bayesian modelling approach allowed the use of some familiar aspects of ANOVA, including mixed models of fixed and random effects, hierarchically nested structures, and variance decomposition, while allowing for overdispersion and excess zeros in the counts. Model selection and estimates of variance components revealed that protection by marine reserves was by far the strongest measured source of variation for relative densities of legal-sized snapper. The size of the effect varied across years among the 3 reserves, with relative densities between 7 and 20 times greater in reserves than in nearby areas. Other than the reserve effect, the temporal factors of season and year were generally more important than the spatial factors at explaining variation in counts. In particular, overall relative densities were ~2 to 3 times greater in autumn than in spring for legal-sized snapper, although the seasonal effect was also variable among locations and years. We consider that the Bayesian generalised linear mixed modelling approach, as used here, provides an extremely useful and flexible tool for estimating the effects of management actions and comparing them directly with other sources of spatial and temporal variation in natural systems.


INTRODUCTION
The exploitation of marine species by humans has caused the depletion of many stocks of fishes worldwide (Pauly et al. 2005, Worm & Branch 2012. N otake marine reserves, designated areas in which all harvesting and damaging of marine life is prohibited (Lubchenco et al. 2003), are increasingly being used as part of the effort to ameliorate this trend. If sufficiently enforced, marine reserves have been shown to increase the size and abundance of exploited species within their borders (Mosqueira et al. 2000, Micheli et al. 2004, Claudet et al. 2010. This, in turn, may produce secondary ecologi-cal effects, such as enhancing populations of exploited species beyond the boundaries of the reserve (Roberts & Polunin 1991, Dugan & Davis 1993, Stoner & Ray 1996, Bohnsack 1998 or facilitating changes in habitat through trophic cascades (Babcock et al. 1999, Shears et al. 2008, Leleu et al. 2012. The value of marine reserves is primarily as a means to manage and protect exploited or endangered species in a particular area, which may then produce broader benefits in terms of increased biodiversity and ecosystem function.
For marine reserves to be used effectively as a management tool, it is critical to be able to estimate and predict their effects. Studies that monitor the abundance of exploited species in existing marine reserves are an essential source of information on which to base such predictions. Accurately quantifying the effects of marine reserves on exploited species can be challenging, however. Data from such studies, often in the form of counts, can be overdispersed or contain excess zeros (Smith et al. 2012), requiring statistical models to be based on nonstandard distributions. Furthermore, marine ecosystems exhibit considerable variation at several temporal and spatial scales (Underwood et al. 2000). Hierarchical sampling regimes that span these scales of variation are therefore necessary to obtain rigorous estimates of the effects of reserves (Andrew & Mapstone 1987, García-Charton & Ruzafa 1999, García-Charton et al. 2000, Willis et al. 2003b. For example, if the abundance of an organism varies from year to year, then a study that spans several years will enable far more accurate estimates of long-term effects as well as provide information on inter-annual variation. The extent to which abundance varies in time and space at different scales is interesting in itself and provides a context of the underlying 'natural' variation with which to compare any measured effect of marine reserves. While some authors have stressed the need to make such comparisons (García-Charton & Ruzafa 1999, García-Charton et al. 2000, 2004, appropriate statistical methods for directly comparing sources of variation in studies of reserve effects have not been explicitly specified. Here, we analyse a long-term, spatially replicated monitoring dataset of counts of snapper Pagrus auratus (Sparidae) from areas inside and outside each of 3 marine reserves in northeastern N ew Zealand. The analysis used a Bayesian approach outlined by Gelman (2005) for ANOVA, which was extended here to more complex zero-inflated generalised linear mixed models (GLMMs). This approach easily incorporated the unbalanced hierarchical structure of the study design in combination with nonstandard error distributions to allow for overdispersion and excess zeros. The primary aim was to estimate the effects of marine reserve protection on counts of snapper while simultaneously accounting for other sources of variation at various spatial and temporal scales. We then compared the estimated reserve effects with other sources of variation in the study design using variance components. The consistency of reserve effects in time and space was also evaluated by estimating interactions between the reserve effect and other factors.

Background and sampling design
Snapper is an important coastal species in temperate northeastern New Zealand, supporting the country's largest inshore commercial and recreational fisheries (Maunder & Starr 2001). Stocks of snapper in this region (fisheries management area known as SN A1) are believed to be slowly rebuilding since being heavily exploited and re duced below the maximum sustainable yield in the latter half of the 20th century (Ministry for Primary Industries 2013). Snapper is also ecologically important, with strong evidence that its predation of sea urchins Evechinus chloroticus can contribute to a trophic cascade that allows the restoration of kelp Ecklonia radiata forests within marine reserves in some contexts (Babcock et al. 1999, Shears et al. 2008. There is also some evidence that small crypto-benthic fishes may be affected by large densities of snapper in a marine reserve (Willis & Anderson 2003).
An ongoing monitoring programme of the relative density of snapper in areas inside and outside (adjacent to) marine reserves at 3 locations in the northeastern bioregion of N ew Zealand, namely Leigh, Tawharanui, and Hahei (Table 1, Fig. 1), began in 1997. These 3 locations have broadly similar habitat and environmental conditions (Shears et al. 2008). Refer to Willis et al. (2003a) for a description and analysis of the first 3 yr of data and Drake (2006) for preliminary Bayesian modelling of the data from Leigh only. The programme used a baitedunderwater-video (BUV) sampling method (Willis & Babcock 2000), which was developed following reports that snapper were differentially attracted to divers within reserves compared to outside reserves, thereby introducing bias into the usual underwatervisual-survey method (Cole 1994). The data are in the form of counts, taken as the maximum number of snapper seen in any 1 frame of a 30 min long underwater video deployment ('MaxN '). This is assumed here to be a measure of the relative density of snapper. Snapper were divided into those below ('sublegal') and above ('legal') the recreational minimum legal size of 27 cm fork length (scheduled to increase to 30 cm in April 2014), and these 2 size classes were modelled separately.
At each of the 3 locations, the coastline was divided into several areas, some falling inside and some falling outside the marine reserves (Fig. 1). Note that at each location, the areas falling outside of the reserve occurred in both directions along the coastline to avoid spatial confounding of areas with reserve effects. Monitoring surveys began in 1997 and were carried out twice per year in each of 2 seasons, spring (primarily September to December) and au tumn (primarily March to June), but surveys were not repeated consistently at all locations after the autumn survey of 1999, yielding an overall sampling design that is highly unbalanced in its cell structure ( Table 2). The 2 seasons were included in the monitoring design because some individuals of this species undergo a seasonal inshore migration which causes inshore densities to increase during summer months and subsequently decline during winter months (Willis et al. 2003a, Willis & Millar 2005. At the time of each survey at a particular location, n = 3 to 6 (usually 4) replicate BUV deployments were done at haphazardly chosen positions within each area. A total of n = 1045 deployments were included in the models described here (Table 2). This sampling design yielded 5 factors: Reserve (fixed with 2 levels, inside and outside), Season (fixed with 2 levels, autumn and spring), Year (random with up to 12 levels), Location (fixed with 3 levels), and Area (random, nested in Reserve × Location, with up to 6 levels per combination of reserve and location; see Fig. 1). We chose to treat Location as a fixed effect because the focus was to estimate the effects for these particular reserves rather than for reserves in general.

Candidate models and model selection
Counts of sublegal and legal snapper from the monitoring programme were analysed using Bayesian zero-inflated GLMMs. Following Gelman (2005), we used the Bayesian approach to model variation associated with the effects of marine reserves, as well as seasons, locations, areas, and years, in a structured ANOVA framework. The term ANOVA is used here to refer specifically to the structuring of the coefficients into 'batches', so that the levels of each categorical factor are grouped together. Variance components are estimated for each batch of coefficients to compare the relative importance of the terms in the model. In our models, Gelman's (2005) framework was extended to include error structures considerably more complex than those of traditional Gaussian AN OVA, as required to account for overdispersion and zero inflation.
Models were implemented using Markov chain Monte Carlo (MCMC) methodology with the software OpenBUGS (Lunn et al. 2009), called from within R (R Development Core Team 2013) by the R2OpenBUGS library (Sturtz et al. 2005). Each model was run with 3 chains, each having a length of 100 000 iterations, from which a burn-in of 50 000 was discarded. The chains were thinned at a rate of 1 in 5, resulting in a total of 30 000 values being kept for each model. Convergence was checked using the Brooks-Gelman-Rubin statistic (Gelman & Rubin 1992, Brooks & Gelman 1998. The full 5-factor experimental design, including all of the potential interactions among factors, was highly complex, having a total of 19 terms ( Table 3). Because of the number of missing cells in the sampling design (Table 2), leading to non-identifiability, heuristics were used first to identify appropriate candidate predictor terms for model selection. Specifically, candidate terms for model selection did not include interactions higher than third order and also did not include the third-order interaction that did not involve Reserves (i.e. Season × Year × Location; see Table 3). A formal model-selection procedure was then used to choose the most favourable model out of hundreds of remaining available candidate models for each of the sublegal-and legal-sized snapper datasets.
The candidate models differed in 2 key respects: the structure of the distribution of errors and the factors (including interactions) that were included. The base distributions considered for the errors were the zero-inflated Poisson (ZIP) and the zero-inflated negative binomial (ZINB). Both the ZIP and the ZINB had the parameters λ, the mean of the Poisson distribution conditional on the absence of excess zero, and π, the probability of the occurrence of an excess zero. In addition, the ZIN B had the parameter δ, which allowed for aggregation (overdispersion) in the counts. For either of these error distributions, the overall mean is given by The conditional mean, λ, was modelled as a linear predictor of candidate terms with a log-link function. Zero inflation, π, was incorporated using 1 of 4 alternative types of models (Smith et al. 2012): (1) no zero inflation (π = 0); (2) constant zero inflation (π = α, where α is a single constant parameter to be estimated); (3) zero inflation linked to the conditional mean (Liu & Chan 2011, Smith et al. 2012; and (4) zero inflation modelled as a separate linear predictor of the candidate terms with a logitlink function. A computing cluster with multiple processors allowed us to conduct a thorough search for the best combination of terms (including 2-and 3-way interactions; see Table 3) for modelling λ, and also π in the case of zero inflation by way of a separate model (type 4 above). The general approach began by fitting the most complex model with the full set of candidate terms (as listed in 48 terms removed, and then the process was repeated. Third-order interactions were removed prior to second-order interactions, in a logical sequence, and no models included interaction terms involving main effects that were not also included in the candidate model. This approach for selecting appropriate terms was done separately for both types of error distributions and for all 4 types of zeroinflated models for estimating π. Model selection was based on the deviance information criterion (DIC; Spiegelhalter et al. 2002, Millar 2009), using half the variance of the posterior deviance for estimating the effective number of parameters, p D (Gelman et al. 2004). Some models were excluded because of very high variance in the posterior distributions of some parameters, which was probably caused by poor identifiability (Omlin & Reichert 1999). The final models chosen, from those that remained, were those with the fewest parameters within 2 units of the lowest DIC score.

Structure of the selected models
For both sublegal-and legal-sized snapper, the count (y) in replicate m in Year l, Area k (nested in Reserve × Location), Location j, Season i, and Reserve status h was best modelled using the negative binomial distribution (2) For sublegal-sized snapper, the linear predictor for the conditional mean was (3) and excess zeros required the use of a separate linear predictor (type 4), namely For legal-sized snapper, the linear predictor for λ was (5) and the model for the zero inflation parameter was Source of variation Abbreviation df Fixed or random Selected for sublegal (S) or legal (L) models Not included as a candidate for model selection, based on preliminary heuristics Table 3. Sources of variation for the full ANOVA model based on all factors in the study design. The abbreviation for each term, as shown, was used to indicate the model parameters associated with that term in the generalised linear mixed models, given in Eqs. (3−5) in the text. Terms that were chosen to be included in the final models of relative densities of legal or sublegal snapper (above or below the recreational minimum 27 cm fork length, respectively), obtained using model selection on the basis of the deviance information criterion, are also provided count, using the parameters γ 0 and γ 1 which were estimated (Type 3). In Eqs. (3−5), β 0 is an overall fitted mean, and the subsequent abbreviations correspond to parameters for individual terms in the model, as indicated in Table 3. Within each factor, the coefficients were centred on zero (see Eq. 9 below), so that estimates of mean counts of interest (e.g. the overall mean within reserves) could be constructed based on the above equations, where the values of λ and π are obtained by adding the appropriate estimates of coefficients to the global mean and back-transforming through the above equations. For example, the overall mean inside reserves μ R for legal snapper was calculated as (7) where the conditional mean count within reserves is λ R = e β 0 + R reserve .

Parameterisation and prior distributions for model terms
Let A be a factor represented by a vector of coefficients β = (β 1 , …, β ᐉ ), where ᐉ is the number of levels in A. If the factor A was fixed, coefficients β 1 to β ᐉ−1 were each given prior distributions β~ Norm(0,100) A sum-to-zero constraint was used for fixed factors, such that 1 coefficient was set to (9) For interactions between fixed and random factors, this constraint was also used for the fixed factor within each level of the random factor. Components of variation for fixed factors were defined as (10) We refer to these as 'variance components' in the following, although for fixed factors these are, strictly speaking, not variances but sums of squared fixed effects divided by the appropriate degrees of freedom.
If A was a random factor, the coefficients were given prior distributions (11) where σ 2 A is common to all coefficients and represents the variance component for factor A. The square roots of variance components for random factors were given standard half-Cauchy priors (Gelman 2006). The dispersion parameter for models with the ZINB distribution was given the prior distribution of δ~ Gamma(0.0001, 0.0001) For type 3 zero-inflated models, the parameters γ 0 and γ 1 were both given the prior distribution of Code for fitting the selected models for both sublegal-and legal-sized snapper in R and OpenBUGS is provided as a Supplement (available at www.intres.com/articles/suppl/m499p203_supp/).

Spatial factors: effect of reserve status and variation among locations and areas
For legal snapper, reserve status was by far the greatest source of variation (Fig. 2). After controlling for variation among locations, areas, seasons, and years, the overall reserve effect (i.e. the ratio of mean MaxN counts in reserve vs. non-reserve areas) was estimated to be 13.4 (see Table 4 for credible intervals). However, the reserve effect differed substantially among locations, as evidenced by inclusion of the Reserve × Location interaction term in the model, with the greatest effect observed at Leigh (effect size of 19.3), followed by Hahei (16.0) and then Tawharanui (7.8). Estimated mean MaxN counts per BUV deployment (mean relative densities) in non-reserve areas were around 0.4 for both Leigh and Tawharanui and 0.2 for Hahei. In protected areas, Leigh had by far the greatest mean relative density at 7.5, compared to ~3 in Tawharanui and Hahei. In contrast, for sublegal snapper, reserve status was not included in the chosen model at all. Instead, the 2 spatial factors (Location and Area) were most important for determining the occurrence of excess zeros in sublegal snapper, with the smaller scale of areas being most important (Fig. 2). For predicting the conditional mean counts of sublegal snapper, temporal factors were most important, especially year. Densities of sublegal fish varied among locations, however, with Leigh and Tawharanui supporting densities ~1.5 times that of Hahei (Table 4).
At the finer spatial scale of areas, mean relative densities of sublegal snapper were similar inside and outside reserves (Fig. 3). Credible intervals around the estimated means of areas were too large to make strong conclusions about fine-scale spatial patterns. However, there was potentially a gradient of increas- ing density from the westernmost area (Area 1) to Cape Rodney (Areas 9 and 10). At Tawharanui, there was little variation in estimated mean densities of sublegal snapper among areas. At Hahei, the highest estimated mean density of sublegal snapper was from the central area of the reserve. For legal snapper, relative densities were consistently very low outside of the reserves at all locations, and there was no apparent trend with proximity to the reserve. The greatest densities of legal snapper were found in the central areas of the reserve at Leigh (Areas 5 and 6), with densities declining steeply towards the eastern and western boundaries of the reserve. Within the reserves at Tawharanui and Hahei, however, there were gradients of increasing density from east to west and west to east, respectively. There did not appear to be any consistent relationship between the densities of sublegal and legal snapper among areas, except perhaps in non-reserve areas at Leigh, where similar spatial patterns were apparent for the 2 size classes.
Regarding causal inferences, we note that the mar ine reserves in the present study were established long before this monitoring programme began, so a before-after-control-impact (BACI-) type design (Underwood 1991), which would have provided stronger evidence for the causal effects of the establishment of the marine reserves, was not possible. Thus, the estimated reserve effects might be due to differences in the habitat or environment between the existing designated reserve and non-reserve areas. However, our general conclusion that the differences observed were caused by the absence of fishing in the reserves is supported by the fact that strong effects for legal-sized snapper oc curred in all 3 marine reserves, there was no spatial pseudo -209 Fig. 2. Variance components plot (Gelman 2005) showing the variation associated with each term in the chosen models, expressed as the estimate of the standard deviation (σ) among levels, for predicting the mean relative density of legal or sublegal (above or below the recreational minimum 27 cm fork length, respectively) snapper. For the sublegal snapper, separate linear predictors were used to model the probability of an excess zero (π) and the conditional mean of the counts (λ), so a separate panel is used for each. Point estimates (means of posterior distributions) are represented by vertical lines, with 50 and 95% credible intervals for the means as thick and thin horizontal lines, respectively  Table 4. Point estimates (mean of the posterior distribution, represented by the set of values given by Markov chain Monte Carlo, MCMC) and 95% credible intervals (0.025 and 0.975 quantiles of the posterior distribution) of the mean relative densities for either sublegal or legal snapper in reserve and non-reserve areas at each of 3 locations. Reserve and non-reserve densities for sublegal snapper were pooled because there was no reserve effect in the model. Estimates of the ratio of reserve to non-reserve densities are also provided for legal snapper as an index of the 'reserve effect'. Point estimates for the ratios were obtained by first calculating the ratios for each MCMC iteration, taking the natural log of the ratios, calculating the mean, and then back-transforming replication of the areas sampled at any of these locations, and no such reserve effects were observed for sublegal-sized snapper.

Temporal factors: seasonal effects and variation among years
Mean counts of sublegal-and legal-sized snapper were greater in autumn than in spring (Table 5). The effect was strong for legal snapper and was the second most important source of variation (Fig. 2), with an estimated seasonal effect size (ratio of densities in autumn vs. spring) of 2.6. The overall effect was less convincing for sublegal snapper, with an estimated effect size of 1.8 and a 95% credible interval that included 1. The models selected for both size classes included interaction terms, indicating that the effect of season differed among years and locations (Table 3). For both size classes, the seasonal effect was greater for Hahei than at other locations, driven by relatively low densities in spring. There was little evidence of a strong seasonal effect on sublegal fish at Tawharanui, but Fig. 3. Fine-scale spatial patterns in the estimated mean relative density (MaxN count from 30-min baited-unterwater-video) of sublegal (triangles) and legal (circles) snapper in areas within 3 locations (see Fig. 1 Table 4) of the mean relative densities for either legal or sublegal snapper in each of 2 seasons at each of 3 locations. Estimates for ratios of seasonal effects were obtained as described for reserve effects in the caption for Table 4 there was only 1 yr in which this location was surveyed in both seasons. Annual variation was also important for both size classes (Fig. 2), with mean relative densities varying substantially among the 12 yr of the study (Fig. 4). The model for sublegal snapper included an interaction between location and year, suggesting that different inter-annual patterns were observed among locations. However, at all locations, the largest densities of sublegal snapper were observed over the period from 1999 to 2001. Inter-annual patterns were consistent among locations for legal snapper, as reflected by the absence of a Location × Year interaction in the model for this size class. Autumnal densities of legal snapper within the reserves at Leigh and Hahei appeared to peak in 2003 and decline thereafter (Fig. 4). The most recent survey in 2010 at Hahei recorded the lowest autumnal density yet recorded at any location. At Tawharanui, densities did not appear to vary substantially for the 4 yr in which surveys were done.

Effects of protection by marine reserves on snapper
Marine reserve protection was by far the most important determinant of the relative density of legal snapper, with the estimated component of variation associated with the reserve effect being much greater than any of the other spatial or temporal factors (Fig. 2). When averaged across all other factors, the relative density of legal snapper was estimated to be 13 times greater inside reserves than outside reserves (Table 4), a result similar to the value of 14 times greater, which was reported from an analysis of the first 3 yr of this monitoring program (Willis et al. 2003a).
Our results indicated large differences in the effects of reserves on legal snapper among locations, which were not reported by Willis et al. (2003a). Differences in reserve effects have also been observed in a recent study of another species (rock lobster Jasus edwardsii) in a set of reserves which included the 3 studied here (Freeman et al. 2012). The largest effect for legal snapper was observed at Leigh, with densities estimated to be nearly 20 times greater within the reserve than outside the reserve, while Tawharanui and Hahei had effect sizes of 8 and 16, respectively. This range in effect sizes also compares favourably with those estimated from the Poor Knights Islands Marine Reserve, located offshore within the same bioregion, where densities of legalsized snapper were estimated to be 22 and 11 times greater than those at 2 comparable non-reserve locations (Denny et al. 2004). The densities of 211 Fig. 4. Inter-annual and seasonal patterns in the estimated mean relative density (MaxN count from 30-min baited-unterwater-video) of sublegal (triangles) and legal (circles) snapper at 3 locations. Open and closed symbols represent the point estimates (means of posterior distributions) for spring and autumn, respectively. Error bars are 95% credible intervals for the means. For legal snapper, estimates for within the reserves only are shown because too few snapper were observed outside the reserves to show any interpretable patterns. Note the different scales on the y-axes for sublegal (left) and legal (right) panels legal snapper outside the reserves at Leigh and Tawharanui were roughly the same (Table 4), which is not surprising given their close proximity and similar environmental conditions. However, the density within the reserve at Leigh was over twice that of the reserves at Tawharanui and Hahei, which were similar to one another. Several potential factors might explain differences in the measured effect sizes of marine reserves placed in different locations. First, theory suggests that the size of a reserve is an important factor determining the extent of the recovery of populations within it (e.g. Kramer & Chapman 1999). Yet, results from recent meta-analyses examining the relationship between the size of the reserve and its effects on populations have been mixed: a positive relationship was evident in some studies (Claudet et al. 2008, Stewart et al. 2009) but not in others (Halpern 2003). In general, the effects of reserve size must be considered in light of the home-range dynamics of a species, as this will influence the proportion of fish that will move into adjacent fished areas (Kramer & Chapman 1999, Moffitt et al. 2009). The spatial dynamics of snapper are complex, as this species shows considerable variation in movement patterns among individuals. Tagging studies of snapper in this region have shown that some snapper make seasonal inshore-offshore migrations, travelling up to tens or even hundreds of kilometres, while others are resident on reefs and move only hundreds of metres (Crossland 1976, Willis et al. 2001, Parsons et al. 2003, Egli & Babcock 2004). Summertime onshore migration of fish that subsequently become resident on inshore reefs is thought to be an important mechanism responsible for increases in densities of adult snapper within reserves (Willis et al. 2001, 2003a, Denny et al. 2004, Willis & Millar 2005. The position of reserves with respect to patterns of onshore migration in this species is therefore likely to be an important factor in determining their success. Patterns of settlement of larvae near reserves could also potentially influence densities of sublegal-and legal-sized snapper in reserves, but post-settlement processes such as mortality and dispersal are likely to moderate the influence of larval supply on adult densities (Freeman et al. 2012). The colonisation of reserves by seasonal migrants potentially allows the number of resident adult snapper to accumulate more rapidly than would be possible through the progression of juvenile fish to adulthood (as documented at the nearby Poor Knights Islands Marine Reserve by Denny et al. 2004), provided the reserve is large enough to protect them once they are resident. Drawing on knowledge from tagging studies of snapper near Leigh and Tawharanui, a recent simulation study concluded that both the Leigh and Tawharanui reserves were of insufficient size to restore densities to unfished levels (Babcock et al. 2012). Thus, the size of the reserves may well be an important factor contributing to the differential reserve effects shown here. Densities of sublegal snapper, and legal snapper outside reserves, were similar at Leigh and Tawharanui (Table 4). Yet, densities of legal snapper were much greater in the larger reserve at Leigh. Legal snapper in the Hahei reserve had a mean density similar to that of the Tawharanui reserve, despite this location having much lower densities of sublegal snapper and legal snapper outside the reserve, consistent with the general southward decrease in the abundance of this species. This could be because the reserve at Hahei is much larger than the other two and, perhaps more importantly (Freeman et al. 2012), has more than twice the offshore extent (Table 1). These patterns are consistent with the hypothesis that the size of the reserves plays a key role in producing the observed variation in their effects.
Second, several recent meta-analyses have demonstrated a positive relationship between the duration of protection and the effect size of marine reserves (Micheli et al. 2004, Claudet et al. 2008, Molloy et al. 2009). If reserve age was an important factor in this study, a trend of increasing density over time inside these reserves would be expected, yet no such trend was present in these data (Fig. 4). While the reserve at Leigh, the oldest of the 3 reserves studied here (Table 1), showed the greatest effect size, it is only 4 yr older than the one at Tawharanui, which showed a markedly smaller effect. One might also expect greater densities in the reserve at Tawharanui than in the one at Hahei, being 11 yr older (Table 1), but they were in fact very similar (Table 4). Thus, differences in the ages of the reserves do not appear to contribute substantially to the differential effect sizes seen here.
A third group of variables potentially responsible for differential effects of marine reserves includes differences in environmental conditions and habitat at these locations (García-Charton et al. 2004, Huntington et al. 2010. Environmental conditions such as water clarity, sedimentation levels, wind fetch, and wave exposure are broadly similar among the locations studied here (Appendix B in Shears et al. 2008). The speed of the local current could potentially influence the distance covered by the bait plume, and thus the number of fish drawn to a BUV, but we consider it unlikely that general current regimes varied significantly among the locations studied. Differences in habitat among locations are more likely to have contributed to the differences observed in this study. Seasonal inshore migrants will presumably be more likely to remain as residents on reefs that are of sufficient size and quality. Moreover, less favourable habitat is expected to support lower densities, as fish would be required to move over greater areas to satisfy their nutritional needs, therefore putting them at greater risk of moving outside of the reserve and into fished areas. The reserves at which the strongest effects were observed, Leigh and Hahei, contain more extensive reefs than the reserve at Tawharanui and include features such as islands (providing shelter and shallow zones) and vertical reef walls. Indeed, the largest densities at Leigh were observed in the central areas of the reserve where these features are located, although larger densities of targeted fish at the centre of a reserve are expected in any event because of the increased risk of them exiting the reserve in areas nearer its borders (Kramer & Chapman 1999).
Finally, differential levels of fishing effort at these locations may also contribute to the differential reserve effects in many ways. Fishing effort in nearby non-reserve areas is likely to be similar among these locations, all of which are very popular with recreational fishers. It has been suggested that more illegal fishing may occur within the reserve at Tawharanui than within the one at Leigh (Babcock et al. 2012). Poor enforcement is thought to be a major issue potentially compromising the effectiveness of marine reserves in many regions of the world (Guidetti et al. 2008). Thus, a lack of compliance to the no-take status may contribute to the relatively modest estimated effect of the reserve at Tawharanui. Commercial fishing, which occurs primarily offshore, might also potentially moderate the numbers of fish available to make the seasonal inshore migration.
We suggest that variation in the estimated effects of these 3 reserves is likely caused by a combination of factors, including size, habitat, degree of compliance with their no-take status, and patterns of inshore migration. Environmental planners need to consider these factors carefully when planning future marine reserves. Perhaps the most important point is that variation in the effects of reserves exists and should be expected, even within the same geographic region. The sources of such variation in snapper clearly require further study.

Temporal and spatial variation in snapper
Other than the reserve effect, temporal factors (season and year) were generally more important than the other spatial factors for predicting relative densities of snapper in this study. In particular, the seasonal effect was strong (Fig. 2), with counts in autumn ~2 to 3 times greater than those in spring (Table 5, Figs. 3 & 4). Seasonal changes in inshore snapper numbers have been documented in many other studies of this species in this region and are thought to be a result of inshore migration for spawning (Francis 1995, Millar et al. 1997, Millar & Willis 1999, Willis et al. 2003a, Willis & Millar 2005. This explanation is consistent with a stronger seasonal effect for legal than sublegal snapper, as found here, because fewer sublegal fish will be reproductively active. The seasonal effect was variable among years and locations, supporting the results of Francis (1995). The effect was notably absent from Tawharanui for sublegal snapper and was strongest at Hahei for both size classes. Although Willis & Millar (2005) found that the seasonal effect for legal snapper was different inside versus outside the marine reserves, no such interaction was apparent in the present analysis. This is due to differences in the structure of the statistical models: Willis & Millar (2005) used an additive identity-link function as opposed to the log-link model presented here. Thus, an interaction may exist on an additive scale but not on a multiplicative scale.
For sublegal snapper, the effects of the spatial factors on the overall density were difficult to interpret because they were split between separate predictors for the excess zeros and the counts, an unfortunate property of this type of zero-inflated model (Smith et al. 2012). However, the pattern of excess zeros was apparently driven by spatial rather than temporal factors and at the finer spatial scale of individual areas in particular (Fig. 2). This indicates that some areas are consistently more likely than others to give counts of zero for sublegal snapper, perhaps because of spatial variation in the suitability of habitat or environmental conditions among areas and locations (Francis 1995, Ross et al. 2007. Inter-annual variation in both size classes was relatively large (Fig. 2), which is consistent with studies showing highly variable recruitment in this species related to temperature (Francis 1993) or prevailing wind patterns (Zeldis et al. 2005). There were peaks in the relative densities of sublegal snapper in 1999 to 2001 and of legal snapper in around 2003. Considering the growth curve for this species ), this may correspond to a strong recruitment pulse observed in the mid-1990s (Maunder & Starr 2001), which then boosted densities of legal fish in reserves in the early 2000s. In years subsequent to 2003, a trend was observed that suggests that snapper densities declined inside reserves. Although these inter-annual patterns may reflect region-wide temporal changes in snapper populations, they might also to some extent be caused by changes in the personnel conducting the monitoring from year to year. N onetheless, it is clear that any attempts to understand temporal trends and make accurate estimates of the effects of reserves or seasons require that reserves be monitored consistently over several years.

Concluding remarks
Here, we demonstrated the use of Bayesian zeroinflated GLMMs for simultaneously quantifying the effects of marine reserves and variation associated with a number of spatial and temporal factors, including 3 locations divided into 26 areas, 2 seasons, and multiple years, in an unbalanced design. The Bayesian approach easily accommodated the hierarchical sampling designs and mixture of fixed and random effects and their interactions in an ANOVAtype analysis while also incorporating various nonstandard error distributions to account for overdispersion and excess zeros, which are a common issue in ecology (see also Smith et al. 2012). Using the output from the MCMC, it was straightforward to estimate effect sizes of interest while accounting for the other factors. The results obtained by our models were generally consistent with those published earlier for this species, with the distinction that interaction terms were also apparent in our models, indicating important variation in the effects of reserves in time and space and at a variety of scales. Rigorous estimates of (and credible intervals for) components of variation attributable to different sources of variation, expressed as the estimated standard deviation among the levels of each factor (Fig. 2), were a particularly useful output from our analysis. Following Gelman (2005), components of variation were calculated for both fixed and random factors so that the relative contribution of all factors and their interactions could be directly compared. This allowed us to ascertain the most important factors for explaining variation in counts of snapper, which complemented the estimation of the effects of interest. The results herein have a wide range of potential benefits, in -cluding greater understanding of the interplay between the effects of management and spatial and temporal ecological patterns, the provision of valuable data for stochastic simulation models of ecosystems, and enabling more accurate predictions for future reserves.
While classical approaches to estimating effect sizes and components of variation in mixed models have been used for many years in ecological studies (Lewis 1978, Underwood & Chapman 1996, Underwood 1997, Anderson & Millar 2004), many authors have noted advantages of the Bayesian approach over its classical counterparts (Ellison 1996, 2004, Clark 2005, Cressie et al. 2009. We refer readers to the recent work of Bolker et al. (2009Bolker et al. ( , 2013, for general comparisons and guidelines for a range of methods for fitting GLMMs, and to Link et al. (2002), for a more directed discussion of the advantages of MCMC and the Bayesian approach. The present study highlights a particular advantage of contemporary Bayesian software (e.g. OpenBUGS) in that it provides modellers with the flexibility to develop new and innovative model structures, such as the linked zero-inflated model used here (Smith et al. 2012). We note that elements of the dataset used here made it particularly well suited to modelling with Bayesian MCMC, such as the highly unbalanced design, the presence of multiple fixed and random effects, and the need for nonstandard error distributions to account for overdispersion and excess zeros. Simultaneously incorporating all these features in a single model using any other approach would be very challenging. Yet, such complexities are common in monitoring data and should not be overlooked. More generally, we consider that our approach provides a useful and flexible framework for placing the effects of management actions, such as protection by marine reserves, into a broader context of natural underlying variation in biological systems.