Accurately estimating correlations between demographic parameters: A comment on Deane et al. (2023)

Abstract Estimating correlations among demographic parameters is an important method in population ecology. A recent paper by Deane et al. (Ecology and Evolution 13:e9847, 2023) attempted to explore the effects of different priors for covariance matrices on inference when using mark‐recovery data. Unfortunately, Deane et al. (2023) made a mistake when parameterizing some of their models. Rather than exploring the effects of different priors, they examined the effects of the use of incorrect equations on inference. In this manuscript, we clearly describe the mistake in Deane et al. (2023). We then demonstrate the use of an alternative and appropriate method and reach different conclusions regarding the effects of priors on inference. Consistent with other recent literature, informative inverse Wishart priors can lead to flawed inference, while vague priors on covariance matrix components have little impact when sample sizes are adequate.

Models that estimate correlated random effects from multivariate normal distributions can be useful for these types of analyses and allow for inference regarding processes such as synchronous responses to environmental stressors (Riecke et al., 2019), the effects of human harvest on survival rates of wildlife populations (Ergon et al., 2018), or latent variation in individual quality (Fay et al., 2022).
When these analyses are performed using Bayesian statistical approaches, prior choice can affect posterior estimates and subsequent inference about the strength of correlations and resulting parameter estimates.Recent papers have sought to explore the consequences of prior choice on estimates of correlations at both the population-and individual-level using capture-recapture data (Fay et al., 2022;Riecke et al., 2019), and concluded that commonly used priors for covariance matrices, such as the conjugate inverse-Wishart distribution, may lead to erroneous inference.Deane et al. (2023) explored a similar line of inquiry using markrecovery data.Specifically, they attempted to examine the effects of three prior structures for covariance matrices, as well as different sample sizes (i.e., number of releases), on posterior estimates of the correlation between survival and band-recovery rates.This approach may allow for inference regarding the effects of hunting on survival in harvested populations.Two of the three priors ('Wishart' and 'Uniform') used in Deane et al. for capture-recovery data are identical to priors described in Riecke et al. (2019) for capture-recapture data.Deane et al. (2023) also described a third prior, which they referred to as the 'Gamma' prior.
Unfortunately, incorrect equations related to the 'Gamma' prior were used in both the text and code of Deane et al. (2023) that led to uninterpretable parameter estimates and associated erroneous inference.
In this Letter to the Editor, we first demonstrate basic errors in equations in Deane et al. (2023) so that others can avoid similar mistakes (Section 2).We then demonstrate an accurate and effective method for placing gamma priors on the standard deviations of a covariance matrix for bivariate normal distributions (Section 3).Third, we discuss some of the errors in interpretation in Deane et al. (2023) after obtaining parameter estimates from correctly parameterized models (Section 4).Finally, we note that many different effective approaches for placing priors on covariance matrices exist, highlight several key areas of agreement with Deane et al. (2023), and describe best practices for avoiding similar errors in the future (Section 5).Deane et al. (2023) estimated correlations (ρ) between logittransformed survival (S) and band-recovery (f) probabilities.Given a correctly specified covariance matrix (Σ) equal to,

| ERROR S IN IMPLEMENTATION
the correct formulation of a precision matrix (i.e., the inverse of a covariance matrix) is, Precision matrices are useful for the calculation of partial correlations, and are used by some MCMC samplers (e.g., JAGS; Plummer, 2003) to parameterize multivariate normal distributions.
In equation ( 7) of their manuscript, Deane et al. (2023) specified the inverse of a covariance matrix (i.e., the precision matrix) as, This differs from Equation (2) and is incorrect (Lauritzen, 1996).
The Authors indicated in the text that this parameterization (equation 7 in Deane et al., 2023) is 'like the default priors specified for multivariate normal distributions in Program MARK (White & Burnham, 1999)'.This statement is inaccurate.The gamma priors used in This is also incorrect (Lauritzen, 1996).To be clear, Deane et al. (2023) assumed that they were placing a gamma prior on the precisions = 1 2 of survival and band-recovery probabilities.In both the text and the code, they assigned this prior to a complex function including all the parameters in the covariance matrix.For example, the authors assign this prior to Unfortunately, this parameterization is incorrect.The disparity between the correct formulation of a precision matrix and the two incorrect equations used by Deane et al. (2023) led to errors in the results of their manuscript.Furthermore, it led to an erroneous conclusion that the use of gamma priors for precisions resulted in estimates of correlations that are more extreme (i.e., further from 0) than truth.In Section 3, we demonstrate that this conclusion was inaccurate.Unfortunately, given the use of incorrect equations in their analyses, we suggest that the estimates and conclusions from Deane et al. (2023) regarding the 'Gamma' prior they described are unreliable.In the next section of this letter, we demonstrate how to appropriately apply gamma priors to standard deviations of a covariance matrix.

| Data simulation and analysis
Here we demonstrate how to specify gamma priors for the standard deviations of demographic parameters estimated using a bivariate normal distribution.For each simulation (n = 500) we generated mean hunting = − 1 and natural = − 1 mortality hazard rates.We then simulated correlated temporal variation ( ) between the two mortality hazard rates given a multivariate normal distribution with mean zero and an appropriate covariance matrix constructed using values of σ and ρ that varied among simulations, We simulated correlated hunting h and natural h mortality hazard rates (Ergon et al., 2018) and transformed those rates to survival (S), natural mortality (η), hunting mortality (κ) and bandrecovery (f) probabilities assuming constant crippling loss (c = .2) and band reporting (b = .5)probabilities to correspond to classic ringrecovery models (Brownie et al., 1985), We then simulated mark-recovery m-arrays (M), where we first specified cell probabilities (P) for each cell of the m-array, and then simulated outcomes (i.e., band recoveries) as a multinomial trial, (2) (5) Ω 11 ∼ gamma(1.001,0.001).

(6)
∼ gamma(20, 80), ∼ uniform( − 0.9, 0), Like Deane et al.'s (2023) 'intensive' monitoring scenario, we released 10k individuals per year for 36 years.We analyzed each simulated dataset (Supplementary Code) using the data generating model and two different ('Uniform' and 'Gamma') priors on the standard deviations to recover parameter estimates from the simulated data.In the first model, which follows the 'Uniform' parameterization described in Riecke et al. (2019), we applied vague uniform priors to the standard deviations of the mortality hazard rates, ∼ uniform(0, 5) .In the second model, we applied vague gamma priors to the standard deviations of the mortality hazard rates, ∼ gamma(1, 1).We specified a uniform prior for the correlation between mortality hazard rates, ∼ uniform(−1, 1) and weakly informative priors for the means of the mortality hazard rates, ∼ normal(−1, 1) for both the 'Uniform' and 'Gamma' approaches.We fixed the values of crippling loss probability, c = .2

| Computational details
We conducted all analyses in R (R Core Team, 2023) and Stan (Stan Development Team, 2024a) using the rstan package (Stan Development Team, 2024b).For each simulation (n = 500) we sampled four MCMC chains for 25k iterations, discarding the first 10k iterations and retaining every 5th sample.In the following figures we report posterior distribution medians from each iteration.We excluded simulations (n = 0) in which either model type did not achieve >1000 effective samples or � R < 1.01 for the correlation (ρ) parameters.This led to a total of 500 simulations that were used to calculate summary statistics.

| Results
Both the 'Uniform' and 'Gamma' parameterizations we describe in Section 3 adequately recovered simulated correlation parameter values (Figure 1).Parameter coverage, or the proportion of the data generating parameter values that fall within the 95% credible interval of the correlation posterior distribution, was very slightly less than acceptable for both the gamma (.936) and uniform (.936) model parameterizations.The average distance between the median of the posterior distribution for the correlation parameter and the values used to generate simulated estimates were .023for models using gamma priors, and .021for uniform priors (Figure 1).If we compare parameter estimates to the actual simulated correlations (i.e., the correlations of the simulated points), coverage improved for both the gamma (.998) and uniform (.998) parameterizations, and the distance between the median of the posterior distributions and the true correlation of the simulated points was minimal for both gamma (−.0003) and uniform (−.002) models.In other words, the models accurately estimated the correlation of the observed data, but the (8) F I G U R E 1 Scatterplots of the median estimated values of correlations between hunting and natural mortality hazard rates from models using uniform (left) and gamma (center) priors for mortality hazard rate standard deviations regressed against the data generating values of ρ, as well as the estimates from the two models regressed against each other (right).Note that these two models produce nearly identical estimates of correlations between mortality hazard rates (right), as there is little difference between the priors for standard deviations, no difference between the priors for correlations, and no dependency between the correlation and standard deviation parameters (Figure 2) for either model.correlation of the observed data may differ slightly from the true underlying correlation due to the duration of the study (T = 36).In summary, both approaches recovered the simulated parameter values and produced nearly identical results (Figure 1).Deane et al. (2023) incorrectly interpreted the meaning of disparate results from their three priors.Specifically, they concluded that 'using Gamma priors will lead to overestimating the magnitude of negative correlation and Wishart priors likely underestimated the magnitude of negative correlation between survival and recovery'.

| ERROR S IN INTERPRE TATI ON AND INFEREN CE
Later in the discussion, they concluded that 'variable prior influence among different priors is an expected outcome of Bayesian estimation'.This is not entirely correct.Our simulations reveal quite clearly that the Uniform and Gamma priors we describe here produced nearly identical estimates (Figure 1).When different vague priors are used with adequate data and implemented appropriately, we should expect posterior estimates to be robust.Critically, major changes in posterior distributions as a result of minor changes in prior distributions should serve as a warning sign for ecologists writing complex models with robust datasets.It is only when informative priors (e.g., the inverse Wishart prior), very 'weak' data, or simply incorrect equations (Deane et al., 2023) are used that we should expect major disparities in inference among priors.
The conjugate inverse-Wishart prior has been repeatedly shown to be problematic in these types of analyses (Alvarez et al., 2014;Fay et al., 2022;Link & Barker, 2005;Riecke et al., 2019).While its use often leads to underestimation of correlations in capturereencounter analyses due to an inherent dependency between the standard deviations and the correlation parameter(s), the shape of the logit-link, and the ranges of typical demographic parameter estimates, this problem can cause either under-or overestimation of correlations (Figure 2).In Section 3, we demonstrated that uniform and gamma priors for standard deviations of mortality hazard rates resulted in nearly identical parameter estimates for correlations (Figure 1).Our results here and elsewhere (Riecke et al., 2019) suggest that these two approaches, which place vague priors on the standard deviations of the demographic parameters and do not induce dependencies between the standard deviations and the correlations (Figure 2), recover accurate estimates of simulated correlations when sample sizes are adequate (Riecke et al., 2019).
As our simulations and basic theory indicate, the use of different vague priors in Bayesian analyses should not induce major variation in posterior estimates if the data are sufficient to estimate the parameters of interest.We suggest that Deane et al.'s (2023) conclusions are a consequence of misunderstandings of the equations and relevant literature underlying the simulations and analyses presented in their manuscript.We recommend that interested readers refer to other work (Alvarez et al., 2014;Ergon et al., 2018;Fay et al., 2022).Finally, we note that many expansions of these model types have been developed relatively recently, and that continued advances in Bayesian software development and statistics have led to additional useful parameterizations.

| CON CLUS ION
The continued expansion of Bayesian software and the careful and  et al., 2014;Riecke et al., 2019).The other two approaches do not induce the same dependency.
In this example, this leads to accurate estimation of correlations among mortality hazard rates when using vague uniform or gamma priors for standard deviations of mortality hazard rates.
priors (Lewandowski et al., 2009) as well as various modifications of the Cholesky decomposition approach (Chen & Dunson, 2003;Dunson, 2008;Fay et al., 2022).When using bivariate normal distributions, researchers can simply assign hyperpriors to each component of the covariance matrix (Riecke et al., 2019).While the approaches that involve inverse Wishart distributions (Huang & Wand, 2013;O'Malley & Zaslavsky, 2008) retain some of the same challenges of the inverse Wishart, all of the approaches listed above have advantages over the conjugate inverse Wishart prior.
Further, they can be implemented using a wide array of Bayesian software such as JAGS (Plummer, 2003), NIMBLE (de Valpine et al., 2017), or Stan (Gelman et al., 2015), and many are readily extendable to covariance matrices larger than 2 × 2. We note that these approaches are complex (Alvarez et al., 2014;Tokuda et al., 2011) and we strongly recommend that ecologists and wildlife managers interested in exploring these concepts use one of the many well-documented existing approaches developed by formally trained statisticians.
Although we have been critical in our letter, we do find areas of agreement with Deane et al. (2023).Like other related manuscripts (Fay et al., 2022;Riecke et al., 2019) We again acknowledge that these techniques are complicated.
As such, we strongly recommend that researchers that are unfamiliar with these approaches consult the literature and confer with one or more statisticians if possible as they implement these model types.Failure to do so can lead to inappropriate inference.As a final note, we strongly encourage data and code sharing as part of the publication process (Jenkins et al., 2023).Much of our response would not have been possible if the code from Deane et al. (2023) was not publicly available.
Program MARK are assigned to standard deviations, not precisions.More importantly, the equations used by Program MARK to specify covariance matrices are valid (G.C. White, personal communication).Deane et al. (2023) used a second incorrect formulation of this equation in the code they wrote to analyze their simulated data as well as the real mallard data, where they formulated the covariance (Σ) and precision (Ω) matrices as,

,
and band reporting probability, b = .5,for both models.Readers should note that our simulation study differs slightly from Deane et al. (2023) in that we explore the effects of different numbers of releases on inference.Furthermore, we use Stan (Stan Development Team, 2024b) rather than JAGS to take advantage of more efficient sampling.Thus, the goal of our simulation was simply to that vague gamma and uniform priors on the standard deviations of mortality hazard rates or probabilities should not have major impacts on inference if researchers use appropriate equations.
thoughtful work of statisticians has led to a plethora of options for placing priors on correlated random effects drawn from multivariate normal distributions.These include but are not limited to; scaled inverse Wishart distributions (O'Malley & Zaslavsky, 2008), hierarchical half-t priors (Huang & Wand, 2013), the separation strategy proposed by Barnard et al. (2000), Lewandoski-Kurowicka-Joe F I G U R E 2 Scatterplots of the natural log of standard deviations (σ) drawn from inverse Wishart (left), Uniform (center), and Gamma (right) priors in which lighter blue colors indicate increased absolute values of the correlation parameter (ρ).Note that the inverse Wishart induces a dependency between the standard deviations and the correlation which can lead to under-or over-estimation of the correlation as a function of the data and the prior (Alvarez ,Deane et al. (2023) demonstrated challenges associated with use of the inverse Wishart prior when estimating correlations among parameters.Similar to other recent work(Fay et al., 2022;Riecke et al., 2019),Deane et al. (2023) demonstrated a need for assessing variation in the annual number of marked individuals for estimating correlations.We note that study duration may also be an important metric to consider when assessing posterior distributions of correlations, and that simply estimating correlations without considering underlying causal processes may lead to flawed inference(Riecke, Lohman, et al., 2022; Riecke, Sedinger, et al., 2022).Perhaps most importantly, we agree with Deane et al. (2023) that power analyses, or simulation-based model validation, are an important component of assessing quality of inference.Deane et al. (2023) wrote that, 'studies should demonstrate through simulation or power analysis that the ecological questions being assessed can be answered with the data and statistical methods employed' We fully agree and add that when complex problems are addressed, practitioners should take additional care to check their equations and code to ensure that the model parameterizations they have chosen are sensible, and that equations are accurate.