2.1 Modeling pieces of evidence at face value

Consider a meta-analysis of RCTs with binary outcomes, where \(y_{1,i}^{(AD)}\) denotes the number of events in the control group of study \(i\) (\(i=1, \ldots, N\)) arising from \(n_{1,i}^{(AD)}\) subjects, and \(y_{2,i}^{(AD)}\) and \(n_{2,i}^{(AD)}\) denote the equivalent quantities in the treatment group.

In addition, we have RWD of patients who have been treated with routine medical care, similarly to the RCT’s control group. We assume that in the RWD the patients’ outcome is the same as in the RCTs. We denote by \(y^{(IPD)}_{1, N+1, j}\) (for \(j=1, \ldots, M\)) the individual participant outcome variable. Moreover, the RWD provides information from several baseline covariates \(x_{1,j}, x_{2,j}, \ldots, x_{p,j}\).

We model the outcome variables \(y_{1,i}^{(AD)}\), \(y_{2,i}^{(AD)}\) and \(y^{(IPD)}_{1, N+1, j}\) with the following Binomial distributions: \[\begin{eqnarray} y_{1,i}^{(AD)} &\sim& \text{Binomial}(p_{1,i}^{(AD)}, n_{1,i}^{(AD)}), \tag{2.1}\\ y_{2,i}^{(AD)} &\sim& \text{Binomial}(p_{2,i}^{(AD)}, n_{2,i}^{(AD)}), \tag{2.2}\\ y^{(IPD)}_{1, N+1, j} &\sim& \text{Binomial}(p^{(IPD)}_{1, N+1, j}, 1), \tag{2.3} \end{eqnarray}\] where \(p_{1,i}^{(AD)}\) and \(p_{2,i}^{(AD)}\) are the event rates for each treatment group and \(p^{(IPD)}_{1, N+1, j} = \text{Pr}\left(y^{(IPD)}_{1, N+1, j} = 1\right)\) the event rate of the IPD-RWD.

Although, in this section we make emphasis in the HMR with dichotomous outcomes, the model can be easily extended e.g. for continues or time to event outcomes.

2.2 Random effects and models for biases

The HMR is built upon the combination of several submodels represented by conditional probability distributions. These submodels break down the diversity of data characteristics and the quality of the evidence we aim to combine. We model these aspects of the between-study variability with two correlated random effects:

\[ \theta_{1, i} = \text{logit}(p_{1,i}) \tag{2.4} \] and \[ \theta_{2, i} = \text{logit}(p_{2,i}) - \text{logit}(p_{1,i}), \tag{2.5} \]

where \(\text{logit}(p) = \ln(p/(1-p))\). The \(\text{logit}(\cdot)\) link is chosen arbitrarily to map probabilities from the (0, 1) interval into the real line. Of course, this is not the only link function that we can apply in HMR, the complementary log-log is an alternative that can be used to model asymmetry.

The random effect \(\theta_{1, i}\) represents an external validity bias and \(\theta_{2, i}\) represents the relative treatment effect. The idea of including the random effect \(\theta_{1,i}\) is to summarize the number of patients’ characteristics and study design features that may influence the treatment effect \(\theta_{2,i}\), but are not directly observed. For this reason we have called \(\theta_{1,i}\) the baseline risk effect of study \(i\). Independently of the study design, every piece of evidence may suffer from quality performance. Therefore, HMR includes a submodel of study quality bias represented by the random weight \(q_i\), where \(q_i \sim \chi^2(\nu)\). The idea is to penalize studies with low quality or high risk of bias. In addition, by using a \(\chi^2(\nu)\) distribution for the random weights \(q_i\) we implicitly defines a marginal t-distribution for the random effects \(\theta_1\) and \(\theta_2\). This bivariate t-distribution makes the model resistant against outlying results.

The HMR explicitly models the RWD intrinsic bias by a random component \(\phi\), which has a mean \(\mu_{\phi}\) and a variance \(\sigma^2_{\phi} = \sigma_1^2/q_{N+1}\). This intrinsic bias could be the result of a multiplicity of uncontrolled factors such as participant selection bias, errors in measurement of outcomes or information bias, intensity bias in the applications of an intervention and so on.

The HMR model is described by the following system of distributions:

\[\begin{eqnarray} \theta_{1, i}|q_i &\sim& \text{Normal}\left(\mu_1, \sigma_1^2/q_i \right),\quad (i = 1, 2, \ldots, N) \tag{2.6} \\ \theta_{2, i} | \theta_{1, i}, q_i &\sim& \text{Normal} \left(\mu_2 - \rho \frac{\sigma_2}{\sigma_1}\, \mu_1 + \rho \frac{\sigma_2}{\sigma_1} (\theta_{1, i}-\mu_1), (1-\rho^2)\sigma_2^2/q_i \right), \tag{2.7} \\ \phi &\sim& \text{Normal}(\mu_{\phi}, \sigma^2_{\phi}), \tag{2.8} \\ \theta_{N+1,1} + \phi|q_{N+1} &\sim& \text{Normal}(\mu_1 + \mu_{\phi} + \beta_1 x_1 + \ldots + \beta_p x_p, \sigma_1^2/q_{N+1}), \tag{2.9}\\ q_i &\sim& \chi^2(\nu), \quad (i = 1, 2, \ldots, N+1) \tag{2.10}. \end{eqnarray}\]

The first two equations (2.6) - (2.7) describe a submodel for the variability between RCTs, which is based on a bivariate scale mixture distribution. Equation (2.7) is the distribution of the baseline risk \(\theta_{1, i}\), which has a baseline mean \(\mu_1\) and variance \(\sigma_1^2/q_i\) relative to the study’s quality weight \(q_i\).

Equation (2.6) represents the distribution of the relative treatment effect \(\theta_{2, i}\) conditional to the baseline risk \(\theta_{1, i}\). This distribution is used to extrapolate a treatment effect for a particular value of \(\theta_{1}\). Clearly, the mean of this distribution changes as a linear function which depends on the study’s baseline risk. Equation (2.7) can also be extended by including further observed patient or design characteristics that may influence the relative treatment effect.

This regression line has the slope \(\sigma_2/\sigma_1 \rho\) and is centered at \(\mu_1\). The intercept is \(\mu_2 - \rho \frac{\sigma_2}{\sigma_1}\, \mu_1\), which can be interpreted as the treatment effect adjusted by the effect of the baseline risk. If the posterior distribution of \(\rho\) is concentrated around zero, the treatment effect is summarized by the posterior of \(\mu_2\).

Equations (2.8) - (2.9) represent a bias model for the RWD, which can be constructed as follows: The RWD has a baseline risk with a mean

\[ E(\theta_{1, N+1} + \phi) = \mu_1 + \mu_{\phi}. \]

The individual participant characteristics \(x_1, x_2, \ldots, x_p\) are used to model risk factors, and to reduce the influence of the patient selection bias in the average bias \(\mu_{\phi}\):

\[ E(\theta_{1, N+1} + \phi |x_1,\ldots, x_p ) = \mu_1 + \mu_{\phi} + \beta_1 x_{1,j} + \ldots + \beta_p x_{p,j} \quad (j = 1, 2,\ldots, M). \quad \tag{2.11} \]

In addition to adjusting for the mean observational bias \(\mu_{\phi}\), equation (2.9) connects the baseline effect of the RCTs \(\mu_1\) with the effect of the IPD. In a similar way, the parameter \(\sigma_1^2\) accounts for the baseline variability across study types. In this context, we call the parameters \(\mu_1\) and \(\sigma_1^2\) shared parameters to highlight that they are estimated across different data and study types.

2.3 Assessing effectiveness in subgroups of patients

The HMR model is a specially designed Bayesian approach to assess the effectiveness \(\theta_2\) in a subgroup characterized by the risk factor \(x_k\) (\(k=1,\ldots,p\)).

The idea is straightforward: For a baseline risk estimate \(\theta_1^k\) representing the subgroup defined by \(x_k\), we apply the regression model (2.7) to estimate \(\theta_2^k\): \[ \theta_2^k = \alpha_0 + \alpha_1 \, (\theta_1^k - \mu_1), \] where \(\alpha_0 = \mu_2 - \rho \frac{\sigma_2}{\sigma_1}\) and \(\alpha_1 =\rho \frac{\sigma_2}{\sigma_1}\). Hence, in principle, the posterior distribution of \(\theta_2^k\) is used to infer effectiveness of the subgroup \(k\).

The crucial problem of this method is the uncertainty of the location parameter \(\theta_1^k\), which involves a partial identified bias correction parameter \(\mu_{\phi}\) between the AD-RCTs and the IPD-RWD. We propose the following location for the baseline risk characterized by \(x_k\):

\[ \theta_1^{k}(B) = \mu_1 + (1-B)\,\mu_{\phi} + \beta_k\, x_k. \tag{2.12} \]

Hence, the effectiveness in a subgroup \(k\) is a function of the amount of bias correction \(B\) given by

\[ \theta_2^k(B) = \alpha_0 + \alpha_1 \, (\theta_1^k(B) - \mu_1). \tag{2.13} \]

It is worth highlighting the following remarks for the location of the baseline risk \(\theta_1^k(B)\):

• \(B = 0\) the parameter \(\mu_1+\mu_{\phi}\) corresponds to the intercept of the logistic regression using the IPD-RWD and \(\theta_1^{k}(0) = \mu_1 + \mu_{\phi} + \beta_k \, x_k\). We proposed this approach in Verde et al. (2016), but it has the drawback that no active bias correction is involved by using \(\mu_{\phi}\). Therefore, if we expect that a subgroup of frail patients has a baseline risk lower than the mean \(\mu_1\), this implies that \(|\beta_k| > |\mu_{\phi}|\), which does not necessarily happen for a IPD-RWD.

• Similarly, for \(B = 1\) the location \(\theta_1^{k}(0)\) centers the IPD-RWD at the mean of RCTs \(\mu_1\), which could not reflect the possible baseline risk of a subgroup of frail patients.

• We proposed \(B = 2\) in Verde (2019) and we showed by simulations that this location \(\theta_1^{k}(2)\) is a good bias correction candidate.

• In Verde et al. (2016) and Verde (2019), we plugged into (2.12) the posterior means of \(\mu_1\), \(\mu_{\phi}\) and \(\beta_k\). This was a pragmatic approach, which clearly underestimates the variability of \(\theta_1^k(B)\).

In general, we cannot be certain about the location of the baseline risk \(\theta_1^{k}(B)\). Therefore, we proceed by giving a prior distribution to \(B\): \[ B \sim \text{Uniform} (B.lower, B.upper), \] where the range of the Uniform distribution gives the amount of bias correction. In the example below, we use \(B.lower = 0\) and \(B.upper = 3\), which goes from \(\mu_1+\mu_{\phi}\) to \(\mu_1-2\,\mu_{\phi}\).

Finally, the join posterior distribution of \(\theta_1^{k}(B)\) and \(\theta_2^{k}(B)\) is presented to assess effectiveness of the subgroup \(k\). We will illustrate this procedure in action in the case study of this section.

2.4 Model diagnostic and conflict of evidence

In the HMR model we assess conflict of evidence between the RCTs and the RWD by the posterior distribution of \(\mu_{\phi}\), where a posterior concentrated at zero with high probability indicates no conflict between study types. Given that \(\mu_{\phi}\) is partially identifiable, a prior-to-posterior sensitivity analysis of \(\mu_{\phi}\) to check if the data contains information to estimate its posterior mean.

A priori, all studies included have a mean of \(E(q_i) = \nu\). We can expect that low quality studies could have a posterior for \(q_i\) substantially less than \(\nu\). For this reason, the studies’ weights \(q_i\) can be interpreted as an adjustment of studies’ internal validity bias.

In the HMR framework we define \(w_i = \nu /q_i\), where large values of \(w_i\) should correspond to outlying studies. In our applications we display the forest plot of the posterior distributions of \(w\). This plot can be used to detect studies that are outliers, i.e. studies that could be in conflict with the other studies.

These model diagnostics are implemented in the function diagnostic() in jarbes. The outcome is a plot showing the prior-to-posterior effect in \(\mu_{\phi}\), and the forest plot for the \(w\)’s posteriors.

2.5 Priors for hyper-parameters and regression coefficients

We complete the model specification by giving the following set of independent weakly informative priors for the hyperparameters \(\mu_1, \mu_2, \mu_{\phi}, \sigma_1, \sigma_2\), \(\rho\) and \(\nu\). For mean and standard deviation parameters we use:

\[ \mu_1, \mu_2, \mu_{\phi} \sim \text{Logistic}(0, 1) \] and \[ \sigma_1 \sim \text{Uniform}(0, 5), \quad \sigma_2 \sim \text{Uniform}(0, 5). \]

The priors for the means \(\mu_1, \mu_2\), and \(\mu_{\phi}\) are non-informative in the sense that in the probability scale, by transforming each mean parameter using \(f(x)= \exp(x)/(1+\exp(x))\), they result in uniform distributions between [0, 1].

We chose to use the uniform distributions for the priors of \(\sigma_1\) and \(\sigma_2\) to weakly constrain the values of these parameters. These priors give equal probabilities within a range between 0 to 5, which includes the case of “non-heterogeneity between studies’ results” to “impossible to combine results”.

The correlation parameter \(\rho\) is transformed by using the Fisher transformation, \[ z = \text{logit}\left( \frac{\rho+1}{2} \right) \] and a Normal prior is used for \(z \sim \text{Normal}(0, 1.25)\). This prior is also weakly informative and it centers the prior of \(\rho\) scale around 0 with uniform probabilities between [-1, 1].

Using independent priors that constrain \(\sigma_1 >0\), \(\sigma_2>0\) and \(|\rho|<1\) guarantees that in each Markov Chain Monte Carlo (MCMC) iteration the variance-covariance matrix of the random effects \(\theta_1\) and \(\theta_2\) is positive definite.

In order to have a marginal finite variance of \(\theta_1\) and \(\theta_2\), the degrees of freedom parameter must be \(\nu>2\). We recommend to use \(\nu = 4\), which strongly penalizes studies with low quality weights \(q_i\) (outlying results) while avoiding infinite marginal variances of \(\theta_1\) and \(\theta_2\). In this setup, \(q_i \sim \chi^2(4)\) with a prior mean of \(q_i\) equal to 4, but the probability to have \(q_i < 4\) is \(59.4\%\), which gives a large chance to observe outliers.

If we combine, i.e. \(N >20\) studies, an alternative approach to give a prior for \(\nu\) is the following: The \(\nu\) parameter is transformed by \[ U = 1/\nu \] and we give a uniform distribution for \(U\): \[ U \sim \text{Uniform}(a, b). \]

For example, candidate values for \(a\) and \(b\) are \(a = 1 / 10\) and \(b = 1 / 3\), which allows to explore random-effects distributions that go from a t-distribution with 3 to 10 degrees of freedom. This prior is designed to favor long-tailed distributions and to explore conflict of evidence in meta-analysis.

Priors for regression coefficients act as variable selection procedures, in the HMR we consider the following hierarchical prior for the regression coefficients \(\beta_1 \ldots \beta_K\): \[ \beta_k \sim \text{Normal}(0, \sigma_{\beta}^2), \quad \sigma_{\beta} \sim \text{Uniform}(0, 5). \tag{2.14} \] This prior models the coefficients as exchangeable with a prior mean equal to zero and unknown variance \(\sigma_{\beta}^2\). This prior produces a shrinkage effect toward 0 and it penalizes the inclusion of uninteresting covariates. In general, we have to scale the covariates or use binary covariates in order to make the assumption of exchangeability of (2.14) reasonable.

2.6 Statistical computations

The posterior distributions of the HMR are calculated with MCMC. We implemented these computations in the R package jarbes (Just a rather Bayesian evidence synthesis) (Verde (2024)). The function hmr() implements the statistical analysis of the HMR, performs model diagnostics, sensitivity analysis of prior-to-posterior, and visualization of results.

The statistical analysis of this section is based on two parallel MCMC simulations with 10,000 iterations for each chain and taks the first 1000 iterations as burn-in period.