{'fraction_non_alphanumeric': 0.040291866880227796, 'fraction_numerical': 0.018935753692827906, 'mean_word_length': 4.631589496893165, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantification of empirical determinacy: the impact of likelihood weighting on posterior location and spread in Bayesian meta-analysis estimated with JAGS and INLA

The popular Bayesian meta-analysis expressed by Bayesian normal-normal hierarchical model (NNHM) synthesizes knowledge from several studies and is highly relevant in practice. Moreover, NNHM is the simplest Bayesian hierarchical model (BHM), which illustrates problems typical in more complex BHMs. Until now, it has been unclear to what extent the data determines the marginal posterior distributions of the parameters in NNHM. To address this issue we computed the second derivative of the Bhattacharyya coefficient with respect to the weighted likelihood, defined the total empirical determinacy (TED), the proportion of the empirical determinacy of location to TED (pEDL), and the proportion of the empirical determinacy of spread to TED (pEDS). We implemented this method in the R package \texttt{ed4bhm} and considered two case studies and one simulation study. We quantified TED, pEDL and pEDS under different modeling conditions such as model parametrization, the primary outcome, and the prior. This clarified to what extent the location and spread of the marginal posterior distributions of the parameters are determined by the data. Although these investigations focused on Bayesian NNHM, the method proposed is applicable more generally to complex BHMs.


Introduction
Bayesian meta-analysis, which synthesizes the evidence from several already published studies, is an indispensable tool for evidence-based medicine. Usually, Bayesian meta-analysis is based on a Bayesian normal-normal hierarchical model (NNHM) [Bender et al., 2018, Friede et al., 2017a,b, Röver, 2020. This is the simplest Bayesian hierarchical model (BHM), but presents problems typical in more complex BHMs .
If different values of the parameters in a likelihood map to a single data model, then these values are indeterminate, or in other words non-identified, and are not estimable from the likelihood alone [Gustafson, 2015, Lewbel, 2019, Poirier, 1998]. The Bayesian approach is known to mitigate this problem, because a likelihood combined with the proper priors for the parameters that are not likelihood identifiable still yields posterior estimates [Eberly and Carlin, 2000, Gelfand and Sahu, 1999, Gustafson, 2015, Kadane, 1975. This raises a legitimate question: to what extent are Bayesian posterior parameter estimates determined by the data?
In classical statistical models, non-identified parameters in the likelihood are known to cause problems [Lele, 2010, Lele et al., 2007, Lewbel, 2019, Skrondal and Rabe-Hesketh, 2004, Sólymos, 2010. One method that detects non-identified parameters in classical statistical models is data cloning [Lele, 2010, Lele et al., 2007, Sólymos, 2010. Data cloning increases the total sample size by artificially replicating the data and focuses on the standard errors of the parameters. If the standard error of a parameter does not decrease at a rate of 1/ √ n with increasing sample size n, such a parameter is deemed to be non-identified [Gustafson, 2015, Lele, 2010, Lele et al., 2007, Sólymos, 2010.
In contrast, non-identification of parameters is not a problem for complex BHMs, because proper priors lead to proper posteriors [Eberly and Carlin, 2000, Gelfand and Sahu, 1999, Gustafson, 2015. However, it is unknown to what extent the posterior estimates of such parameters are determined by the data. In complex BHMs, it can happen that the posterior is not determined by the data at all but is completely determined by the prior. Because the detection of non-identified parameters in BHMs is very challenging, Carlin and Louis [1996] and Eberly and Carlin [2000] suggest that one should avoid BHMs for which the identifiability issues have not been clarified. For this reason, a formal method applicable to Bayesian NNHM and BHMs is needed that can quantify the empirical determinacy of the parameters, thus answering the question of to what extent the posterior estimates are determined by the data.
Recently, Roos et al. [2021] considered Bayesian NNHM and proposed a two-dimensional method for sensitivity assessment based on the second derivative of the Bhattacharyya coefficient (BC). In one dimension, they used a factor f to scale the within-study standard deviation provided by the data. This scaling perturbed the total number of individual patient data provided for Bayesian meta-analysis, which is related to data cloning , Sólymos, 2010. Alternatively, perturbation of the total number of observations could be expressed by a weighted likelihood with w = 1/f 2 [Bissiri et al., 2016, Holmes andWalker, 2017]. An open question remains whether the applicability of the scaling proposed by Roos et al. [2021] could be generalized to weighting other likelihoods, or extended to general-purpose Bayesian estimation techniques such as the Bayesian numerical approximation INLA [Rue et al., 2009], and Bayesian MCMC sampling with JAGS [Plummer, 2016].
We expect that likelihood weighting will affect both the location and the spread of the marginal posterior distributions. Because the method proposed by Roos et al. [2021] quantifies the total impact of the scaling on the posteriors without focusing on location and spread, a refined method is needed. Ideally, such a refined method should quantify what proportion of the total impact of the data on the posterior affects the location and what proportion affects the spread.
To provide a refinement, we consider the second derivative of the BC induced by the weighted likelihood, which quantifies the total empirical determinacy (TED) for each parameter. Based on the properties of both the BC and the rules for the computation of the second derivative, we split TED into two parts: one for location (EDL) and one for the spread (EDS). This enables the quantification of pEDL, the proportion of EDL to TED, and pEDS, the proportion of EDS to TED. This method is implemented in the R package ed4bhm and works with JAGS and INLA. The method proposed quantifies how much weighted data impact posterior parameter estimates and whether location or spread is more affected by the weighting.
To demonstrate how the proposed method works in applications, we consider two case studies and one simulation study. We use INLA and JAGS for the estimation. We consider both centered and non-centered parametrizations for Bayesian NNHM and demonstrate how TED, pEDL, and pEDS depend on the amount of pooling induced by the heterogeneity prior. Moreover, we consider both the Bayesian NNHM applied to log(OR) and a logit model. To demonstrate that the proposed method is useful in applications, for clarifying the empirical determinacy of parameters, we consider three variants of the NNHM combined with different prior assumptions in a simulation study.
In Section 2, we review the theoretical background behind our method and its implementation through R package ed4bhm. In Section 3, we review the case studies and the models assumed. Then, in Section 4 we present the results which demonstrate the performance of the measures TED, pEDL and pEDS in applications. We conclude with a discussion in Section 5.

Methods
In this section, we review the Bayesian NNHM, present the idea of weighting the likelihood, and define the TED measure, which quantifies the total empirical determinacy of posterior parameters. Moreover, we define the proportion of empirical determinacy for posterior location (pEDL) and spread (pEDS). In addition, we specify the implementation of TED, pEDL and pEDS in INLA and JAGS. Finally, a short description of the R package ed4bhm is provided.

Likelihood weighting
Likelihood weighting is applicable to BHMs with various primary outcomes. A base model is the model without weighting (w = 1). Posterior distribution in the base model can be obtained as where π(y | η, θ) is the likelihood, π(θ | η) is the latent field, π(η) is the assumed prior distribution, and η is the set of all the parameters in the model except the random effects θ.
For posteriors of precision parameters, we apply this approximation to log-transformed posteriors. This approach is justified by the invariance of BC.

Quantification of the total empirical determinacy
Following Roos et al. [2021], we define the TED(ψ) measure as the negative second derivative of the Bhattacharyya coefficient (BC) between the base and weighted posterior distributions for each parameter of the model. Here, where π w (ψ | y) is the marginal posterior distribution from the model with weighted likelihood with ψ ∈ {µ, τ, θ 1 , . . . , θ k }. By (H 2 ) = (1 − BC) = −BC , TED evaluates the second derivative of the squared Hellinger distance H 2 . Note that, w = 1 denotes the base model with the original likelihood, which induces the base marginal posterior density π 1 (ψ | y). We quantify the total empirical determinacy (TED) of the marginal posterior of the parameter. TED can be used to compare the empirical determinacy between different parameters ψ and φ by computing a ratio TED(ψ)/TED(φ) [Roos et al., 2021]. Note that TED in Equation (6) measures the impact of likelihood weighting on the whole marginal posterior distribution. We refine this approach and consider the impact of likelihood weighting on location (L) and spread (S) of the marginal posterior distribution in the next section.

Empirical determinacy of location and spread
In this section, we show how to quantify the impact of likelihood weighting on posterior location (L) and spread (S). This gives rise to empirical determinacy of location (EDL) and spread (EDS).
Given the BC between two normal densities in Equation (5) .
where (7) facilitates computation of the proportion of the empirical determinacy of the location (pEDL), which is the proportion of EDL to TED, and the empirical determinacy of the spread (pEDS), which is the proportion of the EDS to TED (Equation (8)) where ψ ∈ {µ, τ, θ 1 , . . . , θ k }.

Implementation of TED in INLA
In order to implement TED(ψ) in INLA, we use the INLA object from the base model (w = 1). The function inla.rerun fits an INLA model with weighted likelihood. We use this function for w = 1−δ and w = 1+δ. Moreover, we extract the summary statistics (mean and standard deviation) of marginal posterior distributions directly from INLA objects. For log-precision, we extract the estimates on internal scale provided in the INLA object.

Implementation of TED in JAGS
For JAGS, we derive a general formula for weighting the likelihood in BHMs fit by MCMC sampling. The posterior distributions for the base and weighted models are given by Equations (2) and (3). From Equation (2) we obtain We then plug in the formula (9) into (3) to get .
The MCMC samples ψ (m) from π 1 (ψ | y) and dev(ψ) values can be extracted from the JAGS object provided by JAGS. To estimate the mean and standard deviation of π w (ψ | y), we compute and where c m = (π(y | ψ)) w−1 [Held and Sabanés Bové, 2020]. For log-precisions, we take the logarithmic transformation of the MCMC sample for the precision parameter and apply Equations (11) and (12) on that transformed sample.

Relative latent model complexity
In the Bayesian NNHM, the interplay between the within-study standard deviation σ i values provided by the data from k studies and the between-study standard deviation τ can be characterized by the effective number of parameters in the model [Spiegelhalter et al., 2002]. Like Roos et al. [2021], we consider the standardized ratio, the relative latent model complexity (RLMC) RLMC defined in Equation (13) attains values between 0 and 1 and expresses the amount of pooling induced by a heterogeneity prior [Gelman and Hill, 2007]. In Section 4.2, a grid of scale parameter values (4.1, 10.4, 18, 31.2, 78.4) for the heterogeneity prior τ induced by the grid of RLMC values fixed at 0.05, 0.25, 0.5, 0.75, 0.95. Whereas HN(4.1), which corresponds to RLMC=0.05, assigns much probability mass to τ values close to 0 and induces much pooling, HN(78.4), which corresponds to RLMC=0.95, assigns less probability mass to values of τ close to 0 and induces little pooling.

Bayesian computation and convergence diagnostics
The MCMC simulations performed in this paper are based on four parallel chains, with a length of 4 × 10 5 iterations. In each chain, we removed the first half of iterations as a burn-in period and from the remaining 2 × 10 5 iterations we kept every 100th iteration in each of the four chains. Our choice of these parameters was guided by raftery.diag for Model C1 (Section 4.4) and supported by .
To assess the convergence to a stationary distribution, we applied Convergence Diagnostics and Output Analysis implemented in the package coda [Plummer et al., 2006]. Moreover, we implemented the rank plots proposed by , which are histograms of posterior draws ranked over all chains and plotted separately for each chain. Nearly uniform rank plots for each chain indicate good mixing of chains. In addition, we used the function n.eff from the package stableGR [Vats and Knudson, 2021], which calculates the effective sample size using the lugsail variance estimators and determines whether Markov chains have converged.

The R package ed4bhm
The open-access R package ed4bhm Empirical determinacy for Bayesian hierarchical models

Data and models
In this section, we review data and models for two case studies. Moreover, we review the design of a simulation study.

Eight schools
To quantify the effect of a coaching program on SAT-V (Scholastic Aptitude Test-Verbal) scores in eight high schools (Table 3 on page 18 of the Supplementary Material), data from a randomized study was pre-analyzed and used for a Bayesian meta-analysis. The data from these eight schools has been used to study the performance of the Bayesian NNHM and to demonstrate typical issues which arise for BHMs [Gelman and Hill, 2007, Gelman et al., 2014.
We consider two parametrizations of the Bayesian NNHM: the centered and non-centered parametrizations . The model with centered parametrization is defined as where j = 1, . . . , 8. This parametrization is used for both INLA and JAGS. On the other hand, the model with the non-centered parametrization [Gelman et al., 2014 reads for j = 1, . . . , 8, which we implemented in JAGS.

Respiratory tract infections
In this section, we review the meta-analysis data set including 22 randomized, controlled clinical trials on the prevention of respiratory tract infections (RTI) by selective decontamination of the digestive tract in intensive care unit patients [Bodnar et al., 2017] that are presented in Table 4 on page 19 of the Supplementary Material.
For the RTI data set, we consider two different models. First, a Bayesian NNHM model and second, a binomial model with logistic transformation where x is a stacked vector attaining value 1 for the treatment and 0 for the control group, z = (f.t, f.c), n = (n.t, n.c), where f.t, f.c, n.t and f.c are vectors of length 22 (see in Table 4 of the Supplementary Material). Originally, the weakly informative HC(1) prior assumption for τ was suggested by Bodnar et al. [2017]. Ott et al. [2021] suggested a prior N(0, 16) for µ and used the HC(1) prior for τ . We assume priors for α and β similar to the prior for µ.

Simulation study
Our simulation study extends the original simulation suggested by Sólymos [2010]. We simulate a sample of random observations of size n = 50 under NNHM with parameters µ = 2.5, σ = 0.2, τ = 0.5. See Section 7 of the Supplementary Material for further details. For the analysis of the simulated data, we use three types of models A, B, and C. Model A is a normal model that does not assume random effects. Models B and C assume a Bayesian NNHM with random effects. Model B assumes that the within-study standard deviation is known and assigns a prior to the between-study standard deviation τ . This defines the usual Bayesian NNHM. In contrast, model C assigns priors to both the within-study (σ) and between-study (τ ) standard deviations. Note that the parameters σ and τ in model C are non-identified, because only the sum of the within-study and between-study variances is identified by the likelihood [Bayarri and Berger, 2004, Gelman et al., 2014, Lele, 2010, Sólymos, 2010.

Sections 3 and 4 of the Supplementary Material consider two motivating examples.
Both examples use a normal likelihood with identified parameters. One example considers the posterior of the mean. The other example focuses on the posterior of the precision. Because the parameters in the likelihoods of both motivating examples are identified, the rate of the decrease of the sample size is close to 1/ √ w, which matches the result provided by data cloning [Lele, 2010, Lele et al., 2007, Sólymos, 2010. In both motivating examples, the estimates of the pEDS are close to 1. This indicates that likelihood weighting affects mostly the spread of the parameters. In fact, both examples motivate that the 1/ √ w rate used by cloning translates into the properties of the pEDS measure.

Eight schools example
We analyze the data from the eight schools (Table 3 on Table 1 on page 11. For JAGS, both the centered (JAGSc, Equation (14)) and non-centered (JAGSnc, Equation (15) This issue is addressed by the estimates of both the total empirical determinacy (TED) and the within-parameter empirical determinacy (EDL and EDS). For example, the values of TED(µ) are larger than those of TED(log(τ −2 )). This indicates that the data impacts the posterior of µ more than the posterior of log(τ −2 ). pEDL and pEDS provide further details and demonstrate that for both µ and log(τ −2 ) the location of the marginal posterior distribution is more impacted by the data than its spread. Whereas the estimates of pEDS(µ) provided by INLA, JAGSc and JAGSnc are close to 0.18, those of pEDS(log(τ −2 )) differ depending on the centered (JAGSc and INLA) and non-centered (JAGSnc) parametrization. This result indicates that the parametrization used for MCMC sampling may affect the impact of the data on the parameter estimates. Figure 1 puts the results of Table 1 on page 11 in a broader context and shows a transition phase plot for the data of the eight schools (Table 3 on Figure 1 demonstrates that TED is always larger for µ (red) than for log(τ −2 ) (blue) across both the grid of RLMC values and the estimation techniques (INLA, JAGSc, and JAGSnc). Independently of the amount of pooling induced by the heterogeneity prior, the posterior of µ is more determined by the data than is the posterior of log(τ −2 ).
The proportion of the estimates of the within-parameter empirical determinacy (pEDS) for the scale of µ (middle panel) and log(τ −2 ) (bottom panel) of Figure 1  with ESS of MCMC samples for the data provided in Table 3 Table  3  78.4 and for three estimation techniques (INLA, JAGSc, JAGSnc), the values of pEDS(µ) remain at a low level of at most 20%. This indicates that the data determine the spread of the posterior of µ less than its location. In contrast, the estimates of pEDS(log(τ −2 )) highly depend on the amount of pooling induced by the heterogeneity prior. The estimates of pEDS(log(τ −2 )) increase from 10% to 90% across the values of RLMC. This means that heterogeneity priors that induce little pooling lead to posteriors of log(τ −2 ) with spreads more determined by the data than for heterogeneity priors that induce much pooling. Although the estimates of pEDS(µ) are similar for the three estimation techniques (INLA, JAGSc, JAGSnc), the values of pEDS(log(τ −2 )) depend on the estimation technique.  (16) with a normal primary outcome and a Bayesian logit model defined in Equation (17) Table 3 on page 15 provides estimates of the empirical determinacy for the simulation study described in Section 3.3, which considers three types of NNHM models (A, B and C), which are fit by INLA and JAGS. The results for JAGS are based on MCMC samples with ESS reported in Table 4 on page 16. For model A, the estimates of TED(µ) and TED(log(γ −2 )) in Table 3 on page 15 are similar. Moreover, pEDS(µ) and pEDS(log(γ −2 )) are high. For example, pEDS(log(γ −2 )) ≥ 0.86 demonstrates that the spread of the posterior of log(γ −2 ) is highly impacted by the data.

Simulation: empirical determinacy
For models B1, B2 and B3, the values of TED(log(τ −2 )) are larger than the values of TED(µ) in Table 3 on page 15 and indicate that the posterior of log(τ −2 ) is more impacted by the data than is the posterior of µ. Similarly to the phase transition in the example of the eight schools discussed in Section 4.2, the values of TED(log(τ −2 )) depend on the heterogeneity prior. Moreover, the large estimates of pEDS(µ) demonstrate that a large proportion of the change in the posterior distribution of µ is due to the change in spread rather than in location. In contrast, the low estimates of pEDS(log(τ −2 )) indicate that the posterior spread is not much determined by the data.  Table 2: Models for RTI: NNHM (above) and logit (below). Mean, sd, 95% CrI, median and TED, EDL, EDS, pEDL and pEDS for marginal posterior distributions for the parameters µ and log(τ −2 ) for NNHM, and α, β and log(τ −2 ) for logit models defined in Equations (16) and (17) Table 4: ESS for the JAGS samples for the models A, B.1, B.2, B.3, C.1, C.2 and C.3 for the simulation study described in Section 3.3. For the MCMC sampling the total number of iterations used is 4 × 10 5 , the burn-in is 2 × 10 5 , thinning = 100.
Models of type C assume priors on the between-study standard deviation and on the within-study standard deviation. The estimation of models C1, C2 and C3 is very challenging. For example, the ESS of MCMC samples for model C1 provided by JAGS (Table 4 on page 16) are very small. Although the parameters σ and τ are non-identified in models C2 and C3, ESS is reasonably high. Table 3 on page 15 shows that TED of the posteriors of log(σ −2 ) and log(τ −2 ) is larger than that of the parameter µ. Again, the values of TED(log(σ −2 )) and TED(log(τ −2 )) depend on heterogeneity priors. For C1 and C2, the estimates of pEDS(µ) differ between INLA and JAGS attaining large values for INLA and low values for JAGS. This indicates that numerical Bayesian approximation (INLA) and MCMC sampling (JAGS) can react differently to models with a non-identified likelihood. For model C1, pEDS(log(σ −2 )) and pEDS(log(τ −2 )) values shown by INLA and JAGS disagree. This is due to the small values of the ESS value of the MCMC samples provided by JAGS. In contrast, these estimates of pEDS agree well for models C2 and C3. The difference between the values of pEDS(log(σ −2 )) and pEDS(log(τ −2 )) for models C2 and C3 indicates that the prior assumptions impact the values of pEDS.

Discussion
We considered two case studies and one simulation study. For the well-known eight-school example we applied Bayesian NNHM, considering INLA analysis and both centered and non-centered parametrizations for JAGS. Moreover, we provided a transition phase plot for the estimates of TED, pEDL and pEDS across a grid of heterogeneity prior scale parameters, which govern the amount of pooling induced by the heterogeneity prior. This provided insights into how TED, pEDL and pEDS change depending on the properties of the heterogeneity prior. For the RTI data set, we used both the Bayesian NNHM applied to log(OR) and a logit model providing insights into how TED, pEDL and pEDS change depending on the primary outcome. To challenge the method proposed, we relaxed the assumption that the within-study standard deviation is known and assumed priors for both the within-study and between-study standard deviation. This provided an insight into how TED, pEDL and pEDS perform when the underlying model has two non-identified parameters for both INLA and JAGS. The proposed method provided novel insight and proved useful in clarifying the empirical determinacy of the parameters in the Bayesian NNHM.
The analysis of two simple motivating examples, normal mean and normal precision, translated the results provided by data cloning, which are based on the 1/ √ n criterion, to the unified pEDS measure. They showed that for identified parameters of non-hierarchical likelihoods the spread of the posterior is mainly affected by the likelihood weighting and leads to pEDS estimates close to 1. We prefer the use of pEDS, because the application of the 1/ √ n criterion to BHMs is not well justified [Jiang, 2017, Lewbel, 2019. This method considerably refines and generalizes the original method proposed by Roos et al. [2021]. We proposed a unified method that quantifies what proportion of the total impact of likelihood weighting on the posterior is due to the change in the location (pEDL) and what proportion is due to the change in the spread (pEDS). This was achieved by matching posteriors moments with those of a normal distribution for the computation of the BC. Note that the normal distribution encapsulates the information about its location and spread in two unrelated parameters. This property proved useful for the definition of pEDL and pEDS. Because likelihood weighting affects both the location and the spread of the marginal posterior distributions, the proposed method is better suited for the quantification of empirical determinacy than other methods that are focused only on the total impact [Kallioinen et al., 2021, Roos et al., 2021.
We successfully applied the proposed method to a non-normal likelihood. This shows that the proposed method can also be applied to other BHMs with different primary outcomes. Moreover, we implemented this method in general-purpose Bayesian estimation programs, such as INLA and JAGS. For JAGS, we developed and applied a technique for the fast and efficient computation of likelihood weighting, which dispenses with re-runs of MCMC chains. All these refinements and generalizations enable future extensions of the method proposed to complex BHMs and to other general-purpose Bayesian estimation programs, such as Stan [Stan Development Team, 2016]. The estimates of the empirical determinacy will help researchers understand to what extent posterior estimates are determined by the data in complex BHMs.
Currently, the weighting approach focuses on weights that are very close to 1. This is very useful to study the impact of the total number of patients, which is induced by withinstudy sample size changes, on the posteriors of parameters in Bayesian NNHM. However, the data cloning approach, which changes the total number of studies included for the Bayesian meta-analysis, indicates that integer weights that are larger than 1 could also be very useful in assessing the impact of the data on BHMs. More work is necessary to extend the proposed method to cope with integer weights.
The proposed method is based on estimates of the mean and standard deviation for location and spread. These posterior descriptive statistics are provided by default by generalpurpose software for Bayesian computation. Moreover, they are used by other modern approaches to quantify the impact of priors on posteriors [Reimherr et al., 2021]. However, for stability reasons,  recommend the use of location and spread estimates based on quantiles. In addition, Vehtari et al. [2017] recommend Pareto smoothed importance sampling (PSIS) to regularize importance weights, which flow into the computation of the location and spread of the weighted posteriors. Therefore, anchoring the method proposed on quantiles and the use of PSIS are two additional goals that need to be addressed by our future research. This paper focuses mainly on Bayesian NNHM. However, the application of the method proposed to other complex BHMs would provide a further insight into the empirical determinacy of posterior parameter estimates in other applications. The openly accessible R package ed4bhm (https://github.com/hunansona/ed4bhm) conveniently facilitates the application of the proposed method in other settings.