While health researchers may have a simple clinical question, e.g. What is the relative treatment effect of an intervention in a population of interest?, the structure of evidence available to answer it may be very complex. These pieces of evidence could be a combination of published results with different grades of quality, observational studies, or real world data (RWD) in general. Consequentially, statistical models for combining disparate pieces of evidence are of a practical importance.
In this section, we will present recent advances in statistical methods for crossdesign synthesis and crossdata synthesis. These methods have been developed to combine studies with different statistical designs, i.e RCTs with RWD, and with different data types, i.e. aggregated data (AD) and individual participant data (IPD).
These crosssynthesis approaches constitute ongoing research, and the aim is to broaden the view of metaanalysis in this complex area. The crucial difference between crosssynthesis approaches and other popular approaches to metaanalysis (e.g. randomeffects models) is the explicit modeling of internal and external validity bias. Ignoring these types of biases in metaanalysis may lead to misleading conclusions.
In some situations, we may have IPD from RWD, but relevant results about treatment comparisons from RCTs are only available as a metaanalysis of AD. Crossdata synthesis combines IPD from RWD with AD from RCTs. The hope is that by combining multiple types of studies and data we can potentially gain new insights from RCTs’ results, which cannot be seen using only a metaanalysis of RCTs.
The following situation illustrates the aim of the Hierarchical MetaRegression (HMR) model presented in this section: A metaanalysis of RCTs showed efficacy of a new intervention compared to routine clinical practice. Although this kind of metaanalysis is at the top in the hierarchy of clinicalbased evidence to prove efficacy, we may be interested in extrapolating RCTs’ results to subgroups of patients that may not meet inclusion criteria (e.g. they may have several comorbidities), but may benefit from this new intervention.
This question about effectiveness in subgroups of patients was called the empty cell problem by Droitcour, Silberman, and Chelimsky (1993), where the authors introduced the term crossdesign synthesis. In the empty cell problem we want to assess effectiveness in a subgroup of patients, but those patients do not meet inclusion criteria in RCTs. An schematic representation of the empty cell problem is given in the following table:
Study Type



Meet inclusion criteria

RCTs

OS

Yes  ADData  IPDRWD 
No  Empty Cell  IPDRWD 
The HMR model aims to give an answer to the empty cell problem by combining aggregated data results from RCTs showing efficacy of an intervention and RWD containing IPD of patients that may benefit from this intervention (Verde et al. (2016), Verde (2019)).
On one hand, the idea behind the HMR is to consider the IPDRWD as representing an observational control group, where each subgroup of interest has its own baseline risk. On the other hand, the key regression aspect in the HMR is the potential relationship between baseline risk and efficacy between the ADRCTs. In this way, given a baseline risk value of a subgroup of patients, we can estimate its corresponding effectiveness. Therefore, the regression model acts as an external validity bias correction for the metaanalysis of RCTs when we are interested in a RWD subgroup.
Consider a metaanalysis 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 IPDRWD.
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.
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 betweenstudy 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/(1p))\). 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 loglog 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 tdistribution for the random effects \(\theta_1\) and \(\theta_2\). This bivariate tdistribution 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} + \phiq_{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.
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 ADRCTs and the IPDRWD. We propose the following location for the baseline risk characterized by \(x_k\):
\[ \theta_1^{k}(B) = \mu_1 + (1B)\,\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 IPDRWD 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 IPDRWD.
Similarly, for \(B = 1\) the location \(\theta_1^{k}(0)\) centers the IPDRWD 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_12\,\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.
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 priortoposterior 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 priortoposterior effect in \(\mu_{\phi}\), and the forest plot for the \(w\)’s posteriors.
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 noninformative 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 “nonheterogeneity 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 variancecovariance 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 randomeffects distributions that go from a tdistribution with 3 to 10 degrees of freedom. This prior is designed to favor longtailed distributions and to explore conflict of evidence in metaanalysis.
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.
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 priortoposterior, 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 burnin period.
In this example we illustrate the HMR in action by performing a crossevidence and crossdata synthesis. We fit a HMR to a metaanalysis of AD and a cohort of patients from an IPDRWD. In the supplementary material of this chapter we provide the R script to run the model and to display the results. The data sets of this example are available in the R package jarbes.
The source data of the AD metaanalysis is a large systematic review of RCTs from the Clinical Practice at NICE (UK and others) (2011). This systematic review investigated the efficacy of adjunctive treatments in managing diabetic foot problems compared with routine medical care only. The resulting posterior mean odds ratio was 1.98 with a 95% posterior interval of 1.42, 2.74, which showed a protective effect of adding adjunctive treatments. The aim of applying the HMR is to quantify if this effect can be generalized to subgroups of patients in an IPDRWD that only received routine medical care. Of course, this statistical analysis only makes sense if the clinical context is adequate.
The cohort of IPDRWD is a set of 260 patients. Some of these patients have similar characteristics to those represented in the RCTs metaanalysis. However, other patients have comorbidities, where one or more risk factors prevent them to participate in the RCTs due to ethical reasons. For example, 118 patients have severe ulcer lesions with a Wagner score 3 and 4, and 77 patients suffer from both, severe ulcer lesions and peripheral arterial disease (PAD). The question is: Can we generalize the benefit observed in the RCTs to the subgroups of patients?
We fitted a HMR model to these data sources by using the function hmr() in jarbes. We start our analysis by visualizing the conflict of evidence between the different types of data and study types. The left panel of Figure 3.1 presents the posterior distribution of \(\mu_{\phi}\), which is the mean bias of the IPDRWD compared to the ADRCTs control groups. The posterior distribution has a substantial probability mass on the right of zero, which indicates that in average the IPDRWD patients present a better prognoses than the ADRCTs control groups. That means that taking the IPDRWD results at face value would be misleading if we aim to combine them with a metaanalysis of ADRCTs.
The right panel of Figure 3.1 presents the posterior distribution of the weights \(w_{i}\) for each study included in the HMR. These posteriors are summarized using a forest plot, where posterior intervals substantially greater than one indicate outliers. One important aspect of the HMR is that those outliers are automatically downweighted in the analysis.
#' hmr
#'
set.seed(2022)
hmr.model.1 = hmr(AD, # Data frame for Aggregated
two.by.two = FALSE, # If TRUE indicates that the trial results are with names: yc, nc, yt, nt
dataIPD = IPD, # Data frame of the IPD
re = "sm", # Random effects model: "sm" scale mixtures
link = "logit", # Link function of the random effects
sd.mu.1 = 1, # Scale parameter for the prior of mu.1
sd.mu.2 = 1, # Scale parameter for the prior of mu.2
sd.mu.phi = 1, # Scale parameter for the prior of mu.phi
sigma.1.upper = 5, # Upper bound of the prior of sigma.1
sigma.2.upper = 5, # Upper bound of the prior of sigma.2
sigma.beta.upper = 5, # Upper bound of the prior of sigma.beta
sd.Fisher.rho = 1.25, # Scale parameter for the prior of rho under Fisher's transformation
df.estimate = TRUE, # If TRUE the degrees of freedom are estimated
df.lower = 3, # Lower bound of the df's prior
df.upper = 10, # Upper bound of the df's prior
nr.chains = 2, # Number of MCMC chains
nr.iterations = 10000, # Total number of iterations
nr.adapt = 1000, # Number of iteration discarded for burnin period
nr.thin = 1) # Thinning rate
summary(hmr.model.1, digits= 2)
#> Model specification:
#> Random effects: sm
#> Link function: logit
#> Split weights: FALSE
#> Estimate degrees of freedom: TRUE
#>
#> Fixed effects:
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> mu.1 0.56 0.21 0.96 0.71 0.57 0.43 0.14 1 1800
#> mu.2 0.66 0.14 0.40 0.57 0.66 0.76 0.95 1 1400
#> intercept 0.59 0.15 0.30 0.48 0.58 0.68 0.90 1 890
#> slope 0.14 0.14 0.41 0.23 0.14 0.05 0.15 1 970
#> Odds.pool 1.82 0.28 1.35 1.62 1.78 1.98 2.47 1 840
#> P_control.pool 0.36 0.05 0.28 0.33 0.36 0.40 0.47 1 1900
#>
#> Random effects:
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> sigma.1 1.04 0.20 0.70 0.90 1.02 1.16 1.47 1 730
#> sigma.2 0.53 0.14 0.29 0.43 0.52 0.62 0.86 1 6200
#> rho 0.27 0.26 0.71 0.46 0.29 0.10 0.28 1 580
#> df 4.07 1.20 3.02 3.26 3.65 4.44 7.65 1 1100
#>
#> Bias parameters:
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> mu.phi 0.79 1.09 1.35 0.09 0.78 1.48 2.96 1 5000
#>
#> 
#> Idividual Participant Data Predictors:
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> sigma.beta 0.64 0.20 0.35 0.51 0.62 0.74 1.08 1 13000
#> beta.PAD 0.69 0.29 1.29 0.88 0.68 0.49 0.14 1 18000
#> beta.neuropathy 0.13 0.35 0.55 0.10 0.13 0.36 0.82 1 18000
#> beta.first.ever.lesion 0.10 0.28 0.64 0.28 0.10 0.09 0.44 1 18000
#> beta.no.continuous.care 0.27 0.28 0.26 0.09 0.27 0.45 0.82 1 8000
#> beta.male 0.20 0.30 0.79 0.40 0.20 0.00 0.38 1 18000
#> beta.diab.typ2 0.88 0.44 0.08 0.57 0.86 1.17 1.80 1 18000
#> beta.insulin 0.10 0.30 0.70 0.29 0.09 0.11 0.49 1 18000
#> beta.HOCHD 0.05 0.31 0.67 0.25 0.04 0.16 0.57 1 6100
#> beta.HOS 0.12 0.31 0.47 0.09 0.11 0.32 0.73 1 18000
#> beta.CRF 0.24 0.32 0.87 0.45 0.24 0.03 0.37 1 18000
#> beta.dialysis 0.47 0.51 1.55 0.78 0.44 0.12 0.46 1 18000
#> beta.DNOAP 0.02 0.38 0.76 0.28 0.02 0.23 0.72 1 18000
#> beta.smoking.ever 0.08 0.29 0.49 0.11 0.08 0.27 0.65 1 3900
#> beta.diabdur 0.01 0.02 0.02 0.00 0.01 0.02 0.05 1 18000
#> beta.wagner.class 1.41 0.31 2.03 1.62 1.41 1.20 0.83 1 18000
#>
#> 
#> MCMC setup (fit using jags): 2 chains, each with 10000 iterations (first 1000 discarded)
#> DIC: 708.16
#> pD: 84.87
# Diagnostic plot
# Analysis of conflict of evidence ................
diagnostic(hmr.model.1, study.names = AD$Study,
title.plot.mu.phi = "Prior to Posterior\nSensitivity Analysis",
title.plot.weights = "Outlier Detection",
post.p.value.cut = 0.1,
lwd.forest = 1, shape.forest = 1,
size.forest = 0.4)
Figure 3.2 displays the results of the submodel corresponding to the IPDRWD that received only medical routine care. The posteriors of the regression coefficients \(\beta_k\) are summarized in a forest plot. This submodel detects risk factors that can reduce the chance of getting healed. We see that the group of patients with a Wagner score greater than 2 have substantially less chance of getting healed compared to the group with lower scores. This can also be observed in the group of patients with PAD.
Interestingly, these subgroups of patients that have lower chances of getting healed are underrepresented in the RCTs populations. Therefore, by combining IPDRWD with ADRCT we can learn new insights about these patients that cannot be learned neither from AD nor from IPD alone.
attach.jags(hmr.model.1, overwrite = TRUE)
# Figure: Posterior distribution of the regression coefficients IPD
# Forest plot for the 95% posterior intervals of the regression coefficients
# Variable names
var.names = names(IPD[1])
var.names = c("PAD", "neuropathy", "first.ever.lesion", "no.continuous.care",
"male", "diab.typ2", "insulin", "HOCHD", "HOS", "CRF", "dialysis",
"DNOAP", "smoking.ever", "diabdur", "Wagner.4")
# Coefficient names
v = paste("beta", names(IPD[1]), sep = ".")
mcmc.x.2 = as.mcmc.rjags(hmr.model.1)
greek.names = paste(paste("beta[",1:15, sep=""),"]", sep="")
par.names = paste(paste("beta.IPD[",1:15, sep=""),"]", sep="")
caterplot(mcmc.x.2,
parms = par.names,
col = "black", lty = 1,
labels = greek.names,
greek = T,
labels.loc="axis", cex =0.7,
style = "plain",reorder = F, denstrip = F)
caterplot(mcmc.x.2,
parms = par.names,
col = "black", lty = 1,
labels = var.names,
labels.loc="above",
greek = F,
cex =0.7,
style = "plain",reorder = F, denstrip = F,
add = T)
abline(v=0, lty = 2, lwd = 2)
The results of the submodel that correlates the baseline risk with the relative treatment effect are presented in Figure 3.3, where results are displayed in the logit scale. This HMR’s submodel is used to extrapolate treatment effects to subgroups of patients. The posterior median and the 95% posterior intervals show that increasing the healing rate tends to decrease the relative treatment effect. In other words, healthier patients benefit less from the adjunctive therapy.
The model is centered at 0.565, corresponding to the posterior mean of \(\mu_1\), the RCTs’ baseline risk. To the right of \(\mu_1\) we have the posterior mean of the IPDRWD \(\mu_1 +\mu_{\phi}\), which has a posterior mean of 0.222. This shows an important bias captured by the introduction of \(\mu_{\phi}\) in the model.
# Generalization of treatment effects
plot(hmr.model.1,
x.lim = c(5, 3),
y.lim = c(2, 6),
x.lab = "Event rate of The Control Group (logit scale)",
y.lab = "No improvement = Effectiveness > Improvement",
title.plot = "HMR: Effectiveness Against Baseline Risk",
Study.Types = c("ADRCTs", "IPDRWD")
)
Figure 3.4 presents the posterior effectiveness contours of \((\theta_1(B), \theta_2(B))\) for the subgroups of patients not included in the RCTs and with low chances of getting healed. On the left panel we have the resulting contour for patients with PAD and on the right panel for patients with Wagner score 3 and 4.
The horizontal axis displays the uncertainty in the location of the baseline risk \(\theta_1(B)\) of these subgroups. This uncertainty resulted from the posterior variability of \(\mu_1\), \(\mu_{\phi}\), \(\beta_k\) and the amount of bias correction \(B\). We can see that for both subgroups the posterior effectiveness \(\theta_2(B)\) is above the horizontal line of no effectiveness for the full range of \(\theta_1(B)\). These results indicates that these subgroup of patients may benefit from this new intervention.
# PDA
p.PDA = effect(hmr.model.1,
title.plot = "Subgroup with PDA",
k = 1, # Regression coefficient
x.lim = c(7, 2.5),
y.lim = c(.5, 2.5),
y.lab = "Effectiveness (logit scale)",
x.lab = "Baseline risk (logit scale)",
kde2d.n= 150, S = 15000,
color.line = "blue",
color.hist = "lightblue",
font.size.title = 8)
# Wanger
p.Wagner = effect(hmr.model.1,
title.plot = "Subgroup with Wagner Score > 2",
k = 15, # Regression coefficient
x.lim = c(7, 2.5),
y.lim = c(.5, 2.5),
y.lab = "Effectiveness (logit scale)",
x.lab = "Baseline risk (logit scale)",
kde2d.n= 150, S = 15000,
color.line = "red",
color.hist = "#DE1A1A", #medium red
display.probability = FALSE,
line.no.effect = 0,
font.size.title = 8)
grid.arrange(p.PDA, p.Wagner, ncol = 2, nrow = 1)
The effect plot can be presented in the OR against probability scales with the following:
# PDA
p.PDA = effect(hmr.model.1,
title.plot = "Subgroup with PDA",
k = 1, # Regression coefficient
x.lim = c(0, 1),
y.lim = c(0, 5),
y.lab = "Effectiveness (Odds Ratio)",
x.lab = "Baseline risk (probability)",
kde2d.n= 150, S = 15000,
color.line = "red",
color.hist = "#DE1A1A", #medium red
display.probability = TRUE,
line.no.effect = 1,
font.size.title = 8)
# Wanger
p.Wagner = effect(hmr.model.1,
title.plot = "Subgroup with Wagner Score > 2",
k = 15, # Regression coefficient
x.lim = c(0, 1),
y.lim = c(0, 5),
y.lab = "Effectiveness (Odds Ratio)",
x.lab = "Baseline risk (probability)",
kde2d.n= 150, S = 15000,
color.line = "blue",
color.hist = "lightblue",
display.probability = TRUE,
line.no.effect = 1,
font.size.title = 8)
grid.arrange(p.PDA, p.Wagner, ncol = 2, nrow = 1)