University of Waterloo Department of Statistics and Actuarial Science Assessing the Health of Richibucto Estuary with the Latent Health Factor Index

The ability to quantitatively assess ecological health is of great interest to those tasked with monitoring and conserving ecosystems. For decades, biomonitoring research and policies have relied on multimetric health indices of various forms. Although indices are numbers, many are constructed based on qualitative procedures, thus limiting the quantitative rigor of the practical interpretations of such indices. The statistical modeling approach to construct the latent health factor index (LHFI) was recently developed. With ecological data that otherwise are used to construct conventional multimetric indices, the LHFI framework expresses such data in a rigorous quantitative model, integrating qualitative features of ecosystem health and preconceived ecological relationships among such features. This hierarchical modeling approach allows unified statistical inference of health for observed sites (along with prediction of health for partially observed sites, if desired) and of the relevance of ecological drivers, all accompanied by formal uncertainty statements from a single, integrated analysis. Thus far, the LHFI approach has been demonstrated and validated in a freshwater context. We adapt this approach to modeling estuarine health, and illustrate it on the previously unassessed system in Richibucto in New Brunswick, Canada, where active oyster farming is a potential stressor through its effects on sediment properties. Field data correspond to health metrics that constitute the popular AZTI marine biotic index and the infaunal trophic index, as well as abiotic predictors preconceived to influence biota. Our paper is the first to construct a scientifically sensible model that rigorously identifies the collective explanatory capacity of salinity, distance downstream, channel depth, and silt–clay content–all regarded a priori as qualitatively important abiotic drivers–towards site health in the Richibucto ecosystem. This suggests the potential effectiveness of the LHFI approach for assessing not only freshwater systems but aquatic ecosystems in general.


Existing Methods to Quantify Ecosystem Health
Assessment of the "health" of an ecosystem is often of great importance to those interested in the monitoring and conservation of ecosystems.Health is a complex concept often involving many diverse factors, and therefore is not straightforward to quantify.A popular method to estimate the health of an ecosystem is through one or more multimetric indices, each of which is a scalar number collapsed from several indicator variables of health, or metrics.Ecosystem health metrics are frequently measures of faunal abundance and diversity.For aquatic ecosystems, these biotic metrics often focus on benthic populations -organisms living on or in the sediment at the bottom of a body of water -since they are useful indicators of underlying health conditions (Bilyard, 1987;Dauer, 1993).For example, the AZTI2 marine biotic index (AMBI) (Borja et al., 2000) is a quantitative measure of health for an estuarine ecosystem based on the sample counts of categorised benthos.Its popularity is evident from its use across the globe, including Africa (Bazairi et al., 2005), Asia (Cai et al., 2012), Europe (Medeiros et al., 2012), North America (Teixeira et al., 2012), and South America (Muniz et al., 2012).
AMBI and other common multimetric indices, e.g.infaunal trophic index (ITI) (Word, 1980), estuarine biotic integrity index (Deegan et al., 1997), benthic response index (Smith et al., 2001), benthic quality index (Rosenberg et al., 2004), infaunal quality index (Mackie, 2009), have the main appeal that they are conceptually simple and thus easily interpretable.They also contain a high amount of biological content from subject-matter scientists being involved at all stages of the design of the index.On the other hand, many of these stages during index construction can involve a non-trivial amount of arbitrariness.Consequently, rigorous evaluation of index reliability and other quantitative aspects is difficult with conventional indices: for example, detecting relationships between health and environmental or impact-related covariates such as water depth or urbanisation; and formally assessing the uncertainty in these estimates of health.Recent multi-step approaches towards addressing such concerns (e.g.Smith et al., 2001;Johnston et al., 2009) do not address propagation of uncertainty from one step to another, thereby resulting in inference that is less reliable than that from an integrated statistical methodology.Chiu & Guttorp (2006) proposed the SHIPSL approach, a statistically enhanced multimetric index construction scheme that improves various quantitative aspects of conventional indices (Dobbie & Dail, to appear), although it and others share unresolved issues such as non-transferability in space or time, and the need for follow-up analyses to determine its relationship with non-faunal (abiotic) variables in method evaluation or policy-making contexts.
Recently Chiu et al. (2011) devised the latent health factor index (LHFI), a novel statistical model-based ecological index aimed to retain the advantages of conventional multimetric indices while addressing some of their shortcomings.The LHFI methodology involves a multi-level analysis of covariance generalised linear mixed-effects (regression) model (e.g.Hoff, 2009), or ANOCOVA GLMM: instead of being treated as measures of health, metrics are regarded as indicators of underlying health conditions; thus, these indicators are regressed as response variables upon a latent health quantity (latent since it is not directly observable) which is site-specific, forming the main level of the regression; health in turn can be regressed upon available covariates, such as environmental (e.g.salinity, silt-clay content) and impact related (e.g.urbanisation) variables, forming the optional sublevel in the model hierarchy.
With data on metrics and covariates, latent health can be estimated as a scalar, so that interpretability is retained; the estimated quantity is the value of the index.Additionally, the effect of the covariates on health can be estimated in a single integrated statistical framework.Importantly, statistical modelling is what directly produces the health index under the integrated LHFI framework, as opposed to being employed merely to select relevant metrics before index construction (e.g. in Deegan et al., 1997) or to evaluate the resulting index (e.g. in Borja et al., 2009).Thus, the LHFI is much more rigorous than conventional indices, as its definition utilises universal modelling practices for the definition of the index; its hierarchical modelling framework also allows for comprehensive statistical inference without the need for sequential analyses through which the propagation of uncertainty is lost from one analysis to the next.As well, the LHFI model, once specified for a set of existing sites, allows for cost effective yet rigorous interpolation of health for a new site: prediction can be accomplished simply with information on covariates at this new site, thus bypassing the expensive benthic taxonomic laboratory procedures that are required to gather the metric data as required by conventional indices.These desirable properties are gained without sacrificing biological integrity which can be embedded through subject-matter expertise in the identification of useful metrics and covariates for constructing the LHFI.It is also straightforward to use the LHFI framework to handle data that have certain spatial and/or temporal features, thus resolving the non-transferability issue of previous indices.
Recently, Schliep and Hoeting (2012) integrated formal point-referenced spatial modelling Banerjee et al. (2004) with LHFI principles to model the hierarchical relationship among four levels of quantities: five-point ordinal health metrics (scaled from "poor" to "excellent"), latent continuous quantities that determine the ordinal metrics, latent health, and drivers of health that pertain to geographical and environmental characteristics.They illustrate the type of unified statistical inference that can be drawn from such an LHFI-based approach for assessing biotic integrity of river basins in Colorado, USA.In contrast, Chiu et al. (2011) and Wu (unpublished) directly model the quantitative health indicators based upon which ordinal metrics such as those of Schliep and Hoeting (2012) are defined.This avoids loss of information due to the mapping of quantitative health metrics to a coarse ordinal scale.
More information on popular indices, the LHFI, and their properties are given by Chiu et al. (2011) and Wu (unpublished).

Previous LHFIs for the Richibucto Estuary
The LHFI modelling methodology was first demonstrated and validated on freshwater ecosystems by Chiu et al. (2011).Following this, Wu (unpublished) applied the methodology successfully to an estuarine ecosystem, utilising the dataset from Lu et al. (2008) on the previously unassessed Richibucto estuary in the Canadian province of New Brunswick.The data were collected in the estuary at 18 sites (Fig. 1).Sites 2, 4-7, and 14 were closest to active oyster farms.Oyster farming activity is perceived to impact site health through its direct influence on sediment properties, although Lu et al. (2008) report that different biotic indicators show different types of association with proximity to oyster farms.For example, macrobenthic faunal abundance is medium for 5 of these 6 sites but high for Site 2; Shannon's diversity is relatively even among all 18 sites with a slight upward trend as the site moves away from the upper channel instead of from an oyster farm.This lack of obvious association remains despite these authors' efforts to consider certain individual species as separate health indicators.This motivated us to develop LHFI models for Richibucto based on indicator metrics (Table 1) used to construct the AMBI and ITI, two popular estuarine ecosystem health indices.Specifically, AMBI and ITI metrics are better tailored to estuarine ecosystems than the generic indicators of abundance, richness, and diversity; and they are more comprehensive than indicators based on single species.However, biotic health indicators alone do not explicitly reveal the collective impact on overall health from abiotic variables, including sediment properties (which, in this case, can be affected by oyster farming activity), oceanographic properties (salinity, water temperature, etc.), and their interactions.The work by Wu (unpublished) is the precursor to our current paper.Wu (unpublished) developed two sets of LHFI models, the first using only metrics from AMBI (denoted by LHFI-A), and the second using metrics from both AMBI and ITI (denoted by LHFI-A-I).These models were considered in a Bayesian statistical framework and implemented using Markov chain Monte Carlo (MCMC) techniques.Both sets of models were successful in that they were able to make reasonable distinctions between health levels at different sites, while allowing rigorous assessment of reliability.As well, the resulting health estimates were justifiable by subject-matter expertise.

What May Influence Richibucto's Health?
AMBI metrics, which for the most part pertain to organic enrichment, were the main focus when Wu (unpublished) constructed the Richibucto LHFIs.This was Table 1: Metrics, based on definitions of the AMBI and ITI, that are used here and by Wu (unpublished) to construct LHFIs for the Richibucto estuary.
Preconceived association AMBI abundance 3 metric w/ health 1 species (including specialist carnivores and some deposit-feeding + tubicolous polychaetes) very sensitive to organic enrichment and present under unpolluted conditions 2 species (including suspension feeders, less selective carnivores and ± 4 scavengers) indifferent to enrichment, always present in low densities with non-significant variations with time 3 species tolerant to excess organic matter enrichment (including − surface deposit-feeding species, e.g.tubicolous spionids) 4 second-order oppotunistic species; mainly small-sized polychaetes: − subsurface deposit-feeders, e.g.cirratulids 5 first-order opportunistic species: deposit-feeders, which proliferate − in reduced sediments ITI abundance metric 5 1 suspension feeders: feed on detritus from the water column and + usually lack sediment grains in their stomach contents 2 interface / surface detrital feeders: obtain the same types of food as + suspension feeders but usually from the upper 0.5 cm of the sediment 3 deposit feeders: invertebrates (including carnivores); generally feed ± from the top few cm of the sediment and feed on encrusted mineral aggregates, deposit particles or biological remains 4 specialised environment feeders: mobile burrowers that feed on − deposited organic material; all adapted to live in highly anaerobic sediment because prior to Wu's analyses, it had been perceived that benthic fauna in Richibucto were related to organic enrichment, as well as freshwater input (salinity gradient), variability of particle size and topography (channel and water depth) (Lu et al., 2008).These non-enrichment aspects of the estuary were observed as covariates alongside benthic fauna; see Lu et al. (2008) for details.
As such, in addition to assessing the health of the Richibucto system, the above LHFI models were used to investigate which and how covariates may influence health as reflected by biotic metrics.As discussed by Chiu et al. (2011) and Wu (unpublished), a thorough understanding of the relationship between covariates and health is key to rigorous yet cost effective interpolation of site health, and could prove to be an enormous asset to scientists.To this end, each of the above LHFI-A and LHFI-A-I models had been implemented with a different combination of covariates.Several covariates and interactions were found to have a statistically significant relationship with health.6For LHFI-A, a model with site-associated covariates log(SC), log(depth), salinity, and interaction log(SC)×log(depth), and another model with a single covariate distance downstream (DD), were the two best-fitting models among those investigated.For LHFI-A-I, there were three best models: two corresponded to the same sets of covariates as those for LHFI-A, and another model with covariates log(depth), log(SC) and their interaction.SC denotes the fraction of silt-clay (grains of size <63 µm) out of the sediment pooled from all (2 to 3) grab sample replicates.Depth is the distance (m) from the water surface to the estuary bed at the location of the site from which grab samples were obtained.Salinity (parts per thousand) was measured based on one in situ water sample obtained at the site.To determine DD (km), first a straight line was extended from the western-most site (Site 1) to the eastern-most site (Site 18); DD of any site is defined as the distance between the site's perpendicular projection onto the straight line and Site 1 (Fig. 1).Other site-specific covariates also considered but found to be insignificant or confounded with others were water temperature ( o C), time of year (September or October), median grain size of sediment, sorting (a unitless measure of variability of grain size) and organic content (%).However, attempts by Wu (unpublished) to include covariates from various best-fitting models together in a single LHFI-A or LHFI-A-I model were unsatisfactory; in such combined models, DD remained highly significant, while all other covariates and their interactions were no longer significant at a reasonable credible level.Indeed, salinity and DD are highly correlated (Table 2), and it was unsurprising that the two are not simultaneously significant.However, no strong correlation exist among log(depth), log(SC), and DD (Table 2), and so why did DD "trump" all others in a combined model, despite non-DD covariates being significant when DD was absent?This question had yet to be addressed.As well, while relationships between health and covariates were quite strong for the LHFI-A models, they were less clear for LHFI-A-I models (significance at credible levels ≈60-85% in the best-fitting LHFI-A-I models as opposed to >90% for LHFI-A).This indicated that the extra data from ITI metrics weakened the overall relationship between health and covariates.
One possible explanation for this phenomenon is that the LHFI construct was appropriate for describing health using AMBI metrics and the available covariates, but ITI metrics have weak ecological relevance to Richibucto.This is plausi-ble from a qualitative perspective, in light of our prior beliefs as stated in the start of Section 2.1.To quantitatively address this, one may determine additional covariates that can be more appropriately paired with ITI metrics, then model these alongside the original covariates.However, this would require further field activities, and we do not pursue it in this article due to cost constraints.
On the other hand, it is also possible that a) the ITI data were too noisy for Wu's specific LHFI models to detect any patterns in their relationship with the covariates, or that b) the data were not too noisy, but these LHFI models were inadequate for revealing a clear relationship among latent health and covariates.Scenario a) is certainly conceivable given the type of study at hand, in which data often involve substantial measurement error.If this were the case, there is unlikely much room for the proposed models to be improved upon within the same modelling framework, given that they are already quite ecologically informative.This leaves us with b) to consider; indeed, Wu's models were merely preliminary models under the general LHFI construct.This could also explain the "trumping" phenomenon in the simpler LHFI models.Specifically, distance likely contained much less measurement error than the other covariates, it being easier to measure with precision than the environmental covariates which are intrinsincally more variable in nature.With an LHFI model that was perhaps too simplistic, an effect on health from distance could therefore manifest itself more clearly than effects from other covariates even if all of them are equally important qualitatively.
Existing data can be used immediately to address b), but they require an improved quantitative framework.In light of our prior beliefs, and the fact that the environmental covariates contained specialised information that distance did not, a more sophisticated LHFI modelling framework yet may prove to be helpful in providing a common thread through health, DD and the other ecologically relevant covariates.
Thus, in the rest of this article, we focus on addressing the "trumping" phenonmenon in the context of b).To do so, we wish to determine if implementing either or both of the following will allow us to properly identify the nature of the relationship among latent health and the available covariates: 1. Introduce a covariance structure for the metric effects, instead of independence which was assumed for the preliminary models to reduce computational burden.
2. Introduce additional level(s) to the regression hierarchy based upon the known associations between the available covariates.
Since these steps pertain to different parts of the LHFI model, below we treat each as a stand-alone investigation.Moreover, we consider LHFI-A models only: introducing extra model complexity for LHFI-A-I models can be impractical for proper inference (e.g. via MCMC), given that AMBI and ITI metrics are dependent according to their definitions, so that extra model parameters are required to account for this.(For example, see Wu (unpublished) for the extra parameters involved even when such dependence was only informally accounted for.)

Modelling Covariance of Metric Effects for LHFI-A
Without additional data, estimating an unstructured covariance matrix may be impractical due to the considerable increase in new parameters in the model.Instead, assuming a structured covariance matrix would be preferred.The form of the LHFI model, for instance, might inherently imply a certain covariance structure.Thus, we examine what covariance structure might be appropriate for the metric effects under previous LHFI-A models, and would then implement them to determine if the extra model complexity can bring about more detectable relationships between latent health and covariates.

Extending the LHFI-A Models for Richibucto
We follow the notation of Wu (unpublished).AMBI metrics, denoted by Y in the LHFI framework, are abundances of five disjoint taxonomic groups.Due to the different preconceived directions of their association with health (Table 1), we split the metrics into two groups: s="−" for Metrics 3-5 (negatively related to health), and s="+" for the remaining metrics.In the LHFI model, each member of Group s is modelled as a multinomial random variable.The link function for the GLMM is a generalised logit for s="+", and an inverted generalised logit for s="−", so that large metric values for s="+" and "−" reflect good/neutral and poor health, respectively.
More precisely, let Y i×j(s)×k , written Y ijks for simplicity, denote the value of the jth metric (nested within the sth metric group) for the kth replicate grab sample at the ith site.Let N ik denote the cardinality (total number of benthic organisms) of the kth replicate sample at the ith site; and p ijs denote the unknown probability of an organism from the ith site belonging to the j(s)th taxonomic group.Thus, we have multinomial distributions Next, let H i denote the latent health of the ith site; and θ s and β j(s) respectively denote the metric group effect and individual metric effect (both unknown) in the regression model.Then, the linear predictor in the LHFI framework is For ( 5), we model θ s as a fixed effect and take θ + =0 (as is customary when considering one of the categories as baseline) to ensure model identifiability, and we model site health and metric effects as random.Specifically, the latent regression of H i is where x i is the vector of a given combination of the aforementioned covariates,7 α 0 and α are the unknown coefficients of the corresponding latent regression, and ε i is the normally distributed regression error with unknown variance σ 2 H .Note that there is overlap and thus dependency between the two multinomials of ( 1) and (2).Wu (unpublished) explained that this dependency is crudely accounted for by θ s ; similarly, the mean-zero β j(s) crudely accounts for the dependency among νs within group s.Thus, independent β j(s) s were assumed.More rigorously, we now replace independence of metric effects by where , Σ is the unknown covariance matrix for β, and "MVN" denotes the multivariate normal distribution.Thus, Σ + (2×2), Σ − (3×3), and Σ ± (2×3) denote the covariance matrices for metric groups positively and negatively related to health, and their cross-covariance matrix, respectively.
We use relatively diffuse8 distributions (with the same parametrisation as Wu (unpublished)) as priors for α 0 , α, θ − (univariate Gaussian with mean 0 and variance 100) and σ 2 H (inverse-Gamma with unit shape and scale).To complete the Bayesian modelling hierarchy, we must specify the structure for Σ, as we now explore.

Dependence Structures of βs and νs
The most general form of Σ is to take it as fully unstructured, and thus it can be assumed to have an inverse-Wishart (IW) distribution where "IW d " denotes the inverse of a d × d Wishart matrix with d degrees of freedom and scale matrix equal to the identity, which is a relatively diffuse prior for a d × d unstructured covariance matrix.Then, one can take advantage of existing MCMC software such as OpenBUGS (http://www.openbugs.info)for straightforward implementation of the LHFI model, although in our experience, non-trivial hierarchical centring is essential to improve MCMC mixing (see Chiu et al., 2011;Wu, unpublished, for details).However, one concern for assuming ( 9) is that with a small dataset from 18 sites each with only 2 to 3 replicate grab samples, an unstructured Σ may be only weakly identifiable.9This issue was encountered by Chiu et al. (2011) for a freshwater benthic dataset that also involved 18 sites with 3 replicates per site, although there were nine metrics altogether.To address this concern, we instead consider a structure for Σ which involves fewer unknown parameters.
To this end, let us momentarily consider a frequentist's viewpoint.Without loss of generality, we may drop the subscript s from p and ν in (3)-( 5) since the value of j determines s.Then, given site i, let Σ ν denote the 5×5 covariance matrix whose (j, j )th element is Σ ν jj =Cov(ν ij , ν ij ) which does not depend on i.In the frequentist context, one can show that β has covariance matrix where J dd is a d × d matrix of 1s.Furthermore, we partition Σ ν into "+", "−" and "±" blocks accordingly.Thus, as an alternative to ( 9), one may wish to consider the structure of Σ ν while limiting the complexity of the resulting structure ultimately assumed for Σ.For this, we assume that the dependence between the vectors of ν + s and ν − s is adequately addressed by θ s in (5), and consequently Σ ν ± is a matrix of 0s.However, we allow Σ ν + and Σ ν − to be unstructured.The form of (10) in light of the above consideration suggests that a reasonable structure for Σ may be Our choice of distributions and hyperparameters in ( 11) and ( 12) yields relatively diffuse priors; it also allows Σ to be decomposed as the sum of the constant matrix ςJ 55 and a random block-diagonal matrix whose blocks are unstructured, thus giving Σ a general form that mimics that in (10) while keeping Σ positive definite.Altogether, our extended model comprises ( 3)-( 8) and ( 11)-( 12).Note that the structure of ( 11)-( 12) for Σ corresponds purely to (8); it will not be the covariance structure in the posterior inference for β.

Results of Implementation
For the investigation, we focus on a single LHFI-A model involving the covariates DD, log(depth), log(SC) and log(depth)×log(SC).Bayesian parameter estimates Figure 2: Brooks-Gelman-Rubin diagnostics (Spiegelhalter et al., 2011) for the MCMC samples of each matrix element of a block diagonal Σ. "Sig.b[i,j]"denotes Σ ij , i.e. the (i, j)th element of Σ. Convergence is suggested by a red curve approaching 1, and green and blue curves approaching the same constant.and credible intervals then are compared to the LHFI-A models from Wu (unpublished) with the same set of covariates.Though, instead of fitting ( 11)-( 12) which would require an implementation outside of OpenBUGS, we implemented a block diagonal Σ (i.e.ς≡0) as a limiting case.Inference for Richibucto latent health as a whole10 from this limiting case was very similar to that from a diagonal Σ as assumed by Wu (unpublished), suggesting reasonable robustness of the health inference to Σ. Additionally, the concern of weak identifiability associated with a non-diagonal Σ proves to be somewhat irrelevant for these data, as our two independently generated MCMC chains mixed readily after a burn-in of around 20,000 iterations: parameters of the non-diagonal Σ required this longer burnin (Fig. 2), although all other model parameters each required a burn-in of only 1,00011 (Fig. 3).
Finally, we observe that the significance of the covariates was also virtually unaffected by assuming a Σ that is more complex than that for Wu (unpublished).

An Additional Level in the Latent Regression
The investigation of Section 3 suggested that the extra complexity from a nontrivial dependence structure among metric effects does not help to clarify the re-Figure 3: Trace plots of MCMC iterations (thinned by 100) from the posterior for fitting (3)-( 9); each covariate in ( 6) has been centred (see Footnote 7).Counter-clockwise from top-left: regression coefficient of the centred interaction log(depth) × log(SC), intercept α 0 , random effect β 2 , fixed effect θ − , latent health H 18 , and standard deviation σ H . Trace plots for all other non-Σ model parameters show similar patterns that suggest convergence after a burn-in of merely 1,000.lationships between latent health and covariates given the current Richibuto data.In this section, we revert to the naive independence assumption for β but introduce extra model complexity through an additional level in the regression of latent health on covariates.
Specifically, although the strong correlation between salinity and DD reflects ecological reasoning for coastal sea waters entering an estuary, it is the only clear empirical relationship detected among the available covariates.Therefore, instead of considering salinity and DD to be complementary covariates, we now take salinity as a response of DD, and in turn, latent health as a response of salinity and the remaining covariates from Section 3.3 (Fig. 4).
Then, model statements of this section include (1)-( 8), and additionally, where Σ = σ 2 β I in (8), I is the identity matrix, and x i in (6) denotes the vector of centred covariates for site i including salinity x sal,i (and possibly other covariates) but excluding DD x DD,i .Hence, ( 6) and ( 13) can be collapsed as follows: where α −sal is α with α sal removed, and similarly for x −sal,i .Thus, (15) regards salinity as an implicit covariate, so that when latent health is explicitly regressed on x −sal,i and DD, the implicit covariate decomposes the total error variation into Hence, a smaller ratio σ 2 H /(α 2 sal σ 2 δ + σ 2 H ) reflects a higher contribution from the implicit covariate towards explaining the total error variation of the latent health regression.
For the assessment and prediction of latent health, the model parameters of main interest include H i , σ H , and α, while the rest of the parameters are regarded as nuisance (Wu, unpublished).
It is evident from the 95% CIs in Table 3 that, when DD is considered a driver of salinity, the two are simultaneously relevant to explaining latent health.In fact,   4), ( 5), ( 1), (2), and (3) of Table 3, respectively.Lu et al. (2008) partition Richibucto sites into six groups according to their benthic community composition: red (lower channel), violet (estuarine mouth), green, blue (lower shallow), turqoise (upper shallow), and pink (upper channel).
Model (3) suggest that the interaction log(depth)×log(SC) is an additional credible driver of latent health, complementing the explanatory capacity of salinity-on-DD.This is the first time that a scientifically sensible model has been successfully constructed to rigorously identify the collective explanatory capacity of salinity, DD, depth, and SC -all regarded a priori as qualitatively important -towards site health in the Richibucto ecosystem.
For Models (1)-(3), the posterior mean for the ratio σ 2 H /(α 2 sal σ 2 δ + σ 2 H ) ranges from around 0.55 to 0.65, respectively; corresponding 95% CIs suggest a smallest ratio for Model (1).Thus, despite the high credibility of the correlation between α sal and α DD in Model (2) and of the influence on health from (the interaction between) depth and SC in Model (3), the least complex Model (1) provides slightly clearer evidence for the explanatory capacity of the salinity-on-DD structure.In terms of the model's predictive power, the least and most complex among the three models share the same deviance information criterion (DIC) (Spiegelhalter et al., 2002) which is slightly smaller (better) than that of Model (2).This predictive power corresponds to observed AMBI metrics (not latent health) being predicted by the model.To assess the model's predictive ability for latent health, one could conduct a simulation study in which unobservable H i values are generated then estimated, although such an approach for hierarchical models has its shortcomings (Marshall & Spiegelhalter, 2007) or requires intensive computations (Dey et al., 1998) that may be impractical.
Instead, we compare CIs for H i among models; in Fig. 5, they appear nearly identical across all Models (1)-( 5), i.e. the inference for health is essentially equally credible across various models.Within model, the ranking of sites according to 0.4 0.5 0.6 0.7 0.8 0.9 their LHFI scores and associated CIs do not coincide with the clustering by Lu et al. (2008) based on similarity in benthic community composition (highly correlated with site location).This indicates that the LHFI approach does not merely represent community composition or site location; instead, it rigorously and comprehensively models biotic indicators, abiotic drivers, the abstract notion of health, and the relationship among them.Note that the health CIs from the 18 sites mutually overlap, suggesting that the small dataset does not allow us to distinguish sites at a high credible level based on health; this was also the case for those models by Wu (unpublished), all with single-level covariates.Despite (i) suboptimal distinguishability and (ii) weaker predictive power for AMBI metrics compared to the single-covariate Models (4)-( 5), our two-level-covariate Models (1)-(3) clearly resolved the earlier counterintuitive phenomenon of covariates not being simul-taneously significant.Indeed, (i) is an improvement over conventional methods in quantitative rigour due to the integrated manner from which our uncertainty estimates are obtained.Moreover, (ii) is of secondary concern when the response of key interest is H i instead of the metrics.Technical note: The LHFI framework is built on the fundamental principles of analysis-of-covariance, so that one can only interpret H i values in a relative sense.However, Chiu et al. (2011) explain that including in the study any site that is qualitatively pre-identified as very healthy or very unhealthy would facilitate the interpretation of the magnitude of H i for an individual i.This is slightly different from the approach of L ópez & Fennessy (2002) who include sites that span the spectrum of individual covariates.Finally, aside from nearly identifical CIs for H i across models, Fig. 6 indicates that the five models perform equally well with respect to the uncertainty (width of CIs) of various nuisance parameters, but with one exception: two-level-covariate models have clearly less uncertainty in their inference for σ β (Fig. 6(f)).As this parameter directly contributes to the uncertainty in the linear predictor ν of AMBI metrics, Fig. 6(f) indicates that having two levels of covariates lead to more reliable inference for the model as a whole.

Conclusions
Unlike conventional multimetric health indices, the integrated LHFI approach yields health scores, asseses the influence of health drivers, and provides their associated uncertainty, all in a single, unified analysis for a given model.
LHFI models by Wu (unpublished) with single-level covariates were satisfactory as preliminary models for the 18 Richibucto sites, but lacking the important ability to rigorously identify relationships between health and abiotic drivers that are deemed ecologically important for the Richibucto estuarine system.Wu (unpublished) proposed, without implementation, two ways to address this issue: (a) to introduce a covariance structure on the random metric effects, and (b) to introduce additional layers of regression given preconceived relationships among the covariates.In this paper, we implemented (a) and (b) with AMBI biotic metrics only, but the approach would be applicable in principle to combining AMBI and ITI biotic metrics.Though, with merely 18 sites in Richibucto, combined AMBI-ITI models by Wu (unpublished) suggest that ITI metrics potentially weaken any signal in the health-covariate relationship.
In this paper, we have found (a) to be ineffective.Various block diagonal covariance structures were considered, the most general of which was reported in this paper.However, none substantially affected the health inference nor improved credible levels of the covariates.
On the other hand, our efforts on (b) proved to be well rewarded.An additional layer of covariates based upon the empirical relationship between salinity and distance downstream allowed the model to identify the simultaneous significance of distance and those abiotic covariates that were previously shown by Wu (unpublished) to be significant only when distance was excluded.Moreover, we also show that model inference is more reliable overall when compared to singlelevel-covariate models.Thus, our two-level-covariate modelling framework more comprehensively exploits the ecological relationship among health, biotic metrics, and abiotic covariates, and it yields less uncertainty in model inference.We implemented three variants of the two-level-covariate model: (a) salinity-on-distance alone, with a priori independent regression coefficients and metric effects, (b) same as (a) but assuming bivariate regression coefficients, and (c) same as (a) but including channel depth and silt-clay content (both on the log scale), as well as their interaction.Overall, (a)-(c) are almost equal in statistical performance, with slightly better predictive power of biotic metrics by (a) and (c).Finally, (a) corresponds to marginally stronger evidence for the two-level structure between salinity and distance to influence site health.
Although field data, especially biotic data, are costly to collect and process for the use in quantitative assessment of ecosystem health, our work has shedded light on one practical concern: more than 18 sites and/or more precise measurements on abiotic covariates are needed in order for the LHFI framework to rigorously distinguish Richibucto sites according to AMBI and/or ITI metrics as indicators of ecosystem health.This point will be considered as part of future health monitoring and conservation efforts for the Richibucto estuarine system.

Figure 1 :
Figure 1: Map of Richibucto estuary with the 18 monitored sites (labelled 'G' followed by the site number); the straight lines illustrate the calculation of distance downstream (DD) for Sites 3 and 5.

Figure 4 :
Figure 4: Graphical depiction of regressing salinity on DD as an additional level in the hierarchical latent health model, while, as usual, health is the response of salinity and other covariates, and AMBI metrics are responses of latent health.

Table 2 :
Sample correlation coefficients among covariates for latent health of

Table 3 :
Selected summary statistics of posterior draws.Boldfaced CI limits suggest that the corresponding parameter differs from 0 with high credibility.from the 99% CI (not shown) for both α sal and α DD in each of Models (1)-(3), suggesting very high credibility for the two-level structure.The CIs from