Binary data example

Introduction

Population adjustment methods such as matching-adjusted indirect comparison (MAIC) are increasingly used to compare marginal treatment effects when there are cross-trial differences in effect modifiers and limited patient-level data. MAIC is based on propensity score weighting, which is sensitive to poor covariate overlap and cannot extrapolate beyond the observed covariate space. Current outcome regression-based alternatives can extrapolate but target a conditional treatment effect that is incompatible in the indirect comparison. When adjusting for covariates, one must integrate or average the conditional estimate over the relevant population to recover a compatible marginal treatment effect.

We propose a marginalization method based on parametric G-computation that can be easily applied where the outcome regression is a generalized linear model or a Cox model. The approach views the covariate adjustment regression as a nuisance model and separates its estimation from the evaluation of the marginal treatment effect of interest. The method can accommodate a Bayesian statistical framework, which naturally integrates the analysis into a probabilistic framework. A simulation study provides proof-of-principle and benchmarks the method’s performance against MAIC and the conventional outcome regression. Parametric G-computation achieves more precise and more accurate estimates than MAIC, particularly when covariate overlap is poor, and yields unbiased marginal treatment effect estimates under no failures of assumptions. Furthermore, the marginalized regression-adjusted estimates provide greater precision and accuracy than the conditional estimates produced by the conventional outcome regression.

General problem

Consider one AB trial, for which the company has IPD, and one AC trial, for which only published aggregate data are available. We wish to estimate a comparison of the effects of treatments B and C on an appropriate scale in some target population P, denoted by the parameter $d_{BC(P)}$ . We make use of bracketed subscripts to denote a specific population. Within the AB population there are parameters $\mu_{A( AB)}$ , $\mu_{B(AB)}$ and $\mu_{C(AB)}$ representing the expected outcome on each treatment (including parameters for treatments not studied in the AB trial, e.g. treatment C). The AB trial provides estimators $\bar{Y}_{A(AB)}$ and $\bar{Y}_{B(AB)}$ of $\mu_{A( AB)}$ , $\mu_{B(AB)}$ , respectively, which are the summary outcomes. It is the same situation for the AC trial.

For a suitable scale, for example a log-odds ratio, or risk difference, we form estimators $\Delta_{AB(AB)}$ and $\Delta_{AC(AC)}$ of the trial level (or marginal) relative treatment effects. This is always as a difference so, for example, for the risk ratio this is on the log scale.

$\Delta_{AB(AB)} = g(\bar{Y}_{B{(AB)}}) - g(\bar{Y}_{A{(AB)}})$

Example analysis

First, let us load necessary packages.

library(boot)      # non-parametric bootstrap in MAIC and ML G-computation
library(copula)    # simulating BC covariates from Gaussian copula
library(rstanarm)  # fit outcome regression, draw outcomes in Bayesian G-computation
library(outstandR)
library(simcovariates)

Data

We consider binary outcomes using the log-odds ratio as the measure of effect. The binary outcome may be response to treatment or the occurrence of an adverse event. For trials AC and BC, outcome $y_n$ for subject $n$ is simulated from a Bernoulli distribution with probabilities of success generated from logistic regression.

For the BC trial, the individual-level covariates and outcomes are aggregated to obtain summaries. The continuous covariates are summarized as means and standard deviations, which would be available to the analyst in the published study in a table of baseline characteristics in the RCT publication. The binary outcomes are summarized in an overall event table. Typically, the published study only provides aggregate information to the analyst.

Use the gen_data() function available with the outstandR package.

library(dplyr)
library(MASS)

N <- 200
allocation <- 2/3      # active treatment vs. placebo allocation ratio (2:1)
b_trt <- log(0.17)     # conditional effect of active treatment vs. common comparator
b_X <- -log(0.5)       # conditional effect of each prognostic variable
b_EM <- -log(0.67)     # conditional interaction effect of each effect modifier
meanX_AC <- c(0.45, 0.45)       # mean of normally-distributed covariate in AC trial
meanX_BC <- c(0.6, 0.6)         # mean of each normally-distributed covariate in BC
meanX_EM_AC <- c(0.45, 0.45)    # mean of normally-distributed EM covariate in AC trial
meanX_EM_BC <- c(0.6, 0.6)      # mean of each normally-distributed EM covariate in BC
sdX <- c(0.4, 0.4)     # standard deviation of each covariate (same for AC and BC)
sdX_EM <- c(0.4, 0.4)  # standard deviation of each EM covariate
corX <- 0.2            # covariate correlation coefficient  
b_0 <- -0.6            # baseline intercept coefficient  ##TODO: fixed value

AC.IPD <- gen_data(N, b_trt, b_X, b_EM, b_0,
                   meanX_AC, sdX, 
                   meanX_EM_AC, sdX_EM, 
                   corX, allocation,
                   family = binomial("logit"))

Similarly, for the aggregate data but with the additional summarise step.

BC.IPD <- gen_data(N, b_trt, b_X, b_EM, b_0,
                   meanX_BC, sdX, 
                   meanX_EM_BC, sdX_EM, 
                   corX, allocation,
                   family = binomial("logit"))

cov.X <- BC.IPD %>%
  summarise(across(starts_with("X"),
                   list(mean = mean, sd = sd),
                   .names = "{fn}.{col}"))

out.B <- dplyr::filter(BC.IPD, trt == 1) %>%
  summarise(y.B.sum = sum(y),
            y.B.bar = mean(y),
            y.B.sd = sd(y),
            N.B = n())

out.C <- dplyr::filter(BC.IPD, trt == 0) %>%
  summarise(y.C.sum = sum(y),
            y.C.bar = mean(y),
            y.C.sd = sd(y),
            N.C = n())

BC.ALD <- cbind.data.frame(cov.X, out.C, out.B)

The true logistic outcome model which we use to simulate the data will be:

$\text{logit}(p_{it}) = \beta_0 + \beta_X (X_3 + X_4) + [\beta_{trt} + \beta_{EM} (X_1 + X_2)] \; \text{I}(trt \neq C)$

So, the log OR for the treatment effect of interest, i.e. AC or BC, is $\beta_{trt}$ (b_trt).

This general format of data sets consist of the following.

`AC.IPD`: Individual patient data

X*: patient measurements
trt: treatment ID (integer)
y: (logical) indicator of whether event was observed

`BC.ALD`: Aggregate-level data

mean.X*: mean patient measurement
sd.X*: standard deviation of patient measurement
y.*.sum: total number of events
y.*.bar: proportion of events
N.*: total number of individuals

Note that the wildcard * here is usually an integer from 1 or the trial identifier B, C.

Our data look like the following.

head(AC.IPD)
#>           X1         X2         X3         X4 trt y
#> 1 0.42066874  1.0957407 0.37118099  1.3291540   1 1
#> 2 0.51227329  0.9079984 0.30560144 -0.1842358   1 0
#> 3 1.00347612  0.8136847 1.25054503  1.1986315   1 0
#> 4 0.05641685  0.5318140 0.54799796  0.6694090   1 0
#> 5 0.42997845  0.5274304 0.09140568  0.1222184   1 0
#> 6 0.06916082 -0.3125090 0.99150666 -0.1102693   1 0

There are 4 correlated continuous covariates generated per subject, simulated from a multivariate normal distribution. Treatment trt 1 corresponds to new treatment A, and 0 is standard of care or status quo C. The ITC is ‘anchored’ via C, the common treatment.

BC.ALD
#>    mean.X1     sd.X1  mean.X2     sd.X2   mean.X3     sd.X3  mean.X4     sd.X4
#> 1 0.568822 0.3972807 0.623277 0.3850408 0.5787528 0.3906499 0.565223 0.3789089
#>   y.C.sum   y.C.bar    y.C.sd N.C y.B.sum   y.B.bar    y.B.sd N.B
#> 1      34 0.5074627 0.5037175  67      34 0.2556391 0.4378691 133

In this case, we have 4 covariate mean and standard deviation values; and the event total, average and sample size for each treatment B, and C.

Output statistics

We will implement for MAIC, STC, and G-computation methods to obtain the marginal treatment effect and marginal variance. The definition by which these are calculated depends on the type of data and outcome scale. For our current example of binary data and log-odds ratio the marginal treatment effect is

$\log\left( \frac{n_B/(N_B-n_B)}{n_C/(N_B-n_{B})} \right) = \log(n_B n_{\bar{C}}) - \log(n_C n_{\bar{B}})$

and marginal variance is

$\frac{1}{n_C} + \frac{1}{n_{\bar{C}}} + \frac{1}{n_B} + \frac{1}{n_{\bar{B}}}$ where $n_B, n_C$ are the number of events in each arm and $\bar{C}$ is the compliment of $C$ so e.g. $n_{\bar{C}} = N_C - n_c$ . Other scales will be discussed below.

Model fitting in R

The outstandR package has been written to be easy to use and essential consists of a single function, outstandR(). This can be used to run all of the different types of model, which we will call strategies. The first two arguments of outstandR() are the individual and aggregate-level data, respectively.

A strategy argument of outstandR takes functions called strategy_*(), where the wildcard * is replaced by the name of the particular method required, e.g. strategy_maic() for MAIC. Each specific example is provided below.

MAIC

Using the individual level data for AC firstly we perform non-parametric bootstrap of the maic.boot function with R = 1000 replicates. This function fits treatment coefficient for the marginal effect for A vs C. The returned value is an object of class boot from the {boot} package. We then calculate the bootstrap mean and variance in the wrapper function maic_boot_stats.

Defining $X_1,X2$ as effect modifiers and $X_3,X_4$ as prognostic variables then the formula used in this model is

$y = X_3 + X_4 + \beta_t X_1 + \beta_t X_2$

which corresponds to the following R formula object passed as an argument to the strategy function.

lin_form <- as.formula("y ~ X3 + X4 + trt*X1 + trt*X2")

outstandR_maic <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_maic(formula = lin_form,
                                     family = binomial(link = "logit")))

The returned object is of class outstandR.

str(outstandR_maic)
#> List of 3
#>  $ contrasts:List of 3
#>   ..$ AB: num -1
#>   ..$ AC: num -2.1
#>   ..$ BC: num -1.1
#>  $ variances:List of 3
#>   ..$ AB: num 0.259
#>   ..$ AC: num 0.16
#>   ..$ BC: num 0.0992
#>  $ CI       :List of 3
#>   ..$ AB: num [1:2] -1.99888 -0.00345
#>   ..$ AC: num [1:2] -2.88 -1.32
#>   ..$ BC: num [1:2] -1.716 -0.481
#>  - attr(*, "CI")= num 0.95
#>  - attr(*, "scale")= chr "log_odds"
#>  - attr(*, "class")= chr [1:2] "outstandR" "list"

We see that this is a list object with 3 parts, each containing statistics between each pair of treatments. These are the mean contrasts, variances and confidence intervals (CI), respectively. The default CI is for 95% but can be altered in outstandR with the CI argument.

A print method is available for outstandR objects for more human-readable output

outstandR_maic
#> Object of class 'outstandR' 
#> Scale: log_odds 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB            -1.00    0.259       -2.00   -0.00345
#> 2 AC            -2.10    0.160       -2.88   -1.32   
#> 3 BC            -1.10    0.0992      -1.72   -0.481

Outcome scale

If we do not explicitly specify the outcome scale, the default is that used to fit to the data in the regression model. As we saw, in this case, the default is log-odds ratio corresponding to the "logit" link function for binary data. However, we can change this to some other scale which may be more appropriate for a particular analysis. So far implemented in the package, the links and their corresponding relative treatment effect scales are as follows:

Data Type	Model	Scale	Argument
Binary	`logit`	Log-odds ratio	`log_odds`
Count	`log`	Log-risk ratio	`log_relative_risk`
Continuous	`mean`	Mean difference	`risk_difference`

The full list of possible transformed treatment effect scales are: log-odds ratio, log-risk ratio, mean difference, risk difference, hazard ratio, hazard difference.

For binary data the marginal treatment effect and variance are

Log-risk ratio

Treatment effect is $\log(n_B/N_B) - \log(n_A/N_A)$ and variance $\frac{1}{n_B} - \frac{1}{N_B} + \frac{1}{n_A} - \frac{1}{N_A}$

Risk difference

Treatment effect is $\frac{n_B}{N_B} - \frac{n_A}{N_A}$ and variance $\frac{n_B}{N_B} \left( 1 - \frac{n_B}{N_B} \right) + \frac{n_A}{N_A} \left( 1 - \frac{n_A}{N_A} \right)$

To change the outcome scale, we can pass the scale argument in the outstandR() function. For example, to change the scale to risk difference, we can use the following code.

outstandR_maic_lrr <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_maic(formula = lin_form,
                                     family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_maic_lrr
#> Object of class 'outstandR' 
#> Scale: log_relative_risk 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.458    0.0870      -1.04      0.120
#> 2 AC           -1.14     0.0507      -1.59     -0.703
#> 3 BC           -0.686    0.0364      -1.06     -0.312

Simulated treatment comparison

Simulated treatment comparison (STC) is the conventional outcome regression method. It involves fitting a regression model of outcome on treatment and covariates to the IPD. IPD effect modifiers are centred at the mean BC values.

$g(\mu_n) = \beta_0 + \beta_X (\boldsymbol{x}_n - \boldsymbol{\theta}) + \boldsymbol{\beta_{EM}} (\beta_{trt} + (\boldsymbol{x_n^{EM}} - \boldsymbol{\theta^{EM}}) ) \; \mbox{I}(trt \neq C)$

where $\beta_0$ is the intercept, $\beta_1$ are the covariate coefficients, $\beta_z$ and $\beta_2$ are the effect modifier coefficients, $z_n$ are the indicator variables of effect alternative treatment. $g(\cdot)$ is the link function e.g. $\log$ .

As already mentioned, running the STC analysis is almost identical to the previous analysis but we now use the strategy_stc() strategy function instead and a formula with centered covariates.

$y = X_3 + X_4 + \beta_t(X_1 - \bar{X_1}) + \beta_t(X_2 - \bar{X_2})$

However, outstandR() knows how to handle this so we can simply pass the same (uncentred) formula as before.

outstandR_stc <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_stc(formula = lin_form,
                                    family = binomial(link = "logit")))
outstandR_stc
#> Object of class 'outstandR' 
#> Scale: log_odds 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB            0.366    1.85        -2.30      3.03 
#> 2 AC           -0.733    1.75        -3.33      1.86 
#> 3 BC           -1.10     0.0992      -1.72     -0.481

Change the outcome scale

outstandR_stc_lrr <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_stc(formula = lin_form,
                                     family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_stc_lrr
#> Object of class 'outstandR' 
#> Scale: log_relative_risk 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB            0.212    0.699       -1.43      1.85 
#> 2 AC           -0.474    0.662       -2.07      1.12 
#> 3 BC           -0.686    0.0364      -1.06     -0.312

For the last two approaches, we perform G-computation firstly with a frequentist MLE and then a Bayesian approach.

Parametric G-computation with maximum-likelihood estimation

G-computation marginalizes the conditional estimates by separating the regression modelling from the estimation of the marginal treatment effect for A versus C. First, a regression model of the observed outcome $y$ on the covariates $x$ and treatment $z$ is fitted to the AC IPD:

$g(\mu_n) = \beta_0 + \boldsymbol{x}_n \boldsymbol{\beta_X} + (\beta_z + \boldsymbol{x_n^{EM}} \boldsymbol{\beta_{EM}}) \; \mbox{I}(trt \neq C)$

In the context of G-computation, this regression model is often called the “Q-model.” Having fitted the Q-model, the regression coefficients are treated as nuisance parameters. The parameters are applied to the simulated covariates $x*$ to predict hypothetical outcomes for each subject under both possible treatments. Namely, a pair of predicted outcomes, also called potential outcomes, under A and under C, is generated for each subject.

By plugging treatment C into the regression fit for every simulated observation, we predict the marginal outcome mean in the hypothetical scenario in which all units are under treatment C:

$\hat{\mu}_0 = \int_{x^*} g^{-1} (\hat{\beta}_0 + x^* \hat{\beta}_1 ) p(x^*) \; \text{d}x^*$

To estimate the marginal or population-average treatment effect for A versus C in the linear predictor scale, one back-transforms to this scale the average predictions, taken over all subjects on the natural outcome scale, and calculates the difference between the average linear predictions:

$\hat{\Delta}^{(2)}_{10} = g(\hat{\mu}_1) - g(\hat{\mu}_0)$

outstandR_gcomp_ml <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_gcomp_ml(formula = lin_form,
                                         family = binomial(link = "logit")))
outstandR_gcomp_ml
#> Object of class 'outstandR' 
#> Scale: log_odds 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.966    0.242       -1.93   -0.00274
#> 2 AC           -2.06     0.142       -2.80   -1.33   
#> 3 BC           -1.10     0.0992      -1.72   -0.481

Change the outcome scale

outstandR_gcomp_ml_lrr <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_gcomp_ml(formula = lin_form,
                                         family = binomial(link = "logit")),
            scale = "log_relative_risk")

outstandR_gcomp_ml_lrr
#> Object of class 'outstandR' 
#> Scale: log_relative_risk 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.463    0.0869      -1.04      0.115
#> 2 AC           -1.15     0.0505      -1.59     -0.708
#> 3 BC           -0.686    0.0364      -1.06     -0.312

Bayesian G-computation with MCMC

The difference between Bayesian G-computation and its maximum-likelihood counterpart is in the estimated distribution of the predicted outcomes. The Bayesian approach also marginalizes, integrates or standardizes over the joint posterior distribution of the conditional nuisance parameters of the outcome regression, as well as the joint covariate distribution.

Draw a vector of size $N^*$ of predicted outcomes $y^*_z$ under each set intervention $z^* \in \{0, 1\}$ from its posterior predictive distribution under the specific treatment. This is defined as $p(y^*_{z^*} \mid \mathcal{D}_{AC}) = \int_{\beta} p(y^*_{z^*} \mid \beta) p(\beta \mid \mathcal{D}_{AC}) d\beta$ where $p(\beta \mid \mathcal{D}_{AC})$ is the posterior distribution of the outcome regression coefficients $\beta$ , which encode the predictor-outcome relationships observed in the AC trial IPD. This is given by:

$p(y^*_{^z*} \mid \mathcal{D}_{AC}) = \int_{x^*} p(y^* \mid z^*, x^*, \mathcal{D}_{AC}) p(x^* \mid \mathcal{D}_{AC})\; \text{d}x^*$

$= \int_{x^*} \int_{\beta} p(y^* \mid z^*, x^*, \beta) p(x^* \mid \beta) p(\beta \mid \mathcal{D}_{AC})\; d\beta \; \text{d}x^*$

In practice, the integrals above can be approximated numerically, using full Bayesian estimation via Markov chain Monte Carlo (MCMC) sampling.

The average, variance and interval estimates of the marginal treatment effect can be derived empirically from draws of the posterior density.

We can draw a vector of size $N^*$ of predicted outcomes $y^*_z$ under each set intervention $z^*$ from its posterior predictive distribution under the specific treatment.

outstandR_gcomp_stan <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_gcomp_stan(formula = lin_form,
                                           family = binomial(link = "logit")))
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 5.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.54 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.199 seconds (Warm-up)
#> Chain 1:                0.21 seconds (Sampling)
#> Chain 1:                0.409 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.191 seconds (Warm-up)
#> Chain 2:                0.197 seconds (Sampling)
#> Chain 2:                0.388 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 1.2e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.183 seconds (Warm-up)
#> Chain 3:                0.186 seconds (Sampling)
#> Chain 3:                0.369 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 1.3e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.2 seconds (Warm-up)
#> Chain 4:                0.191 seconds (Sampling)
#> Chain 4:                0.391 seconds (Total)
#> Chain 4:
outstandR_gcomp_stan
#> Object of class 'outstandR' 
#> Scale: log_odds 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.980    0.234       -1.93    -0.0315
#> 2 AC           -2.08     0.135       -2.80    -1.36  
#> 3 BC           -1.10     0.0992      -1.72    -0.481

Change the outcome scale

outstandR_gcomp_stan_lrr <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_gcomp_stan(formula = lin_form,
                                           family = binomial(link = "logit")),
            scale = "log_relative_risk")
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 1.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.198 seconds (Warm-up)
#> Chain 1:                0.207 seconds (Sampling)
#> Chain 1:                0.405 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.2e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.188 seconds (Warm-up)
#> Chain 2:                0.196 seconds (Sampling)
#> Chain 2:                0.384 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 1.3e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.186 seconds (Warm-up)
#> Chain 3:                0.182 seconds (Sampling)
#> Chain 3:                0.368 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 1.3e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.2 seconds (Warm-up)
#> Chain 4:                0.186 seconds (Sampling)
#> Chain 4:                0.386 seconds (Total)
#> Chain 4:

outstandR_gcomp_stan_lrr
#> Object of class 'outstandR' 
#> Scale: log_relative_risk 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.457    0.0829      -1.02      0.107
#> 2 AC           -1.14     0.0465      -1.57     -0.720
#> 3 BC           -0.686    0.0364      -1.06     -0.312

Multiple imputation marginalisation

ref

$equation here$

outstandR_mim <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_mim(formula = lin_form,
                                    family = binomial(link = "logit")))
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 1.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 1001 / 4000 [ 25%]  (Sampling)
#> Chain 1: Iteration: 1400 / 4000 [ 35%]  (Sampling)
#> Chain 1: Iteration: 1800 / 4000 [ 45%]  (Sampling)
#> Chain 1: Iteration: 2200 / 4000 [ 55%]  (Sampling)
#> Chain 1: Iteration: 2600 / 4000 [ 65%]  (Sampling)
#> Chain 1: Iteration: 3000 / 4000 [ 75%]  (Sampling)
#> Chain 1: Iteration: 3400 / 4000 [ 85%]  (Sampling)
#> Chain 1: Iteration: 3800 / 4000 [ 95%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.196 seconds (Warm-up)
#> Chain 1:                0.588 seconds (Sampling)
#> Chain 1:                0.784 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1001 / 4000 [ 25%]  (Sampling)
#> Chain 2: Iteration: 1400 / 4000 [ 35%]  (Sampling)
#> Chain 2: Iteration: 1800 / 4000 [ 45%]  (Sampling)
#> Chain 2: Iteration: 2200 / 4000 [ 55%]  (Sampling)
#> Chain 2: Iteration: 2600 / 4000 [ 65%]  (Sampling)
#> Chain 2: Iteration: 3000 / 4000 [ 75%]  (Sampling)
#> Chain 2: Iteration: 3400 / 4000 [ 85%]  (Sampling)
#> Chain 2: Iteration: 3800 / 4000 [ 95%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.19 seconds (Warm-up)
#> Chain 2:                0.568 seconds (Sampling)
#> Chain 2:                0.758 seconds (Total)
#> Chain 2:
outstandR_mim
#> Object of class 'outstandR' 
#> Scale: log_odds 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.960    0.230       -1.90    -0.0206
#> 2 AC           -2.06     0.130       -2.77    -1.35  
#> 3 BC           -1.10     0.0992      -1.72    -0.481

Change the outcome scale

outstandR_mim_lrr <-
  outstandR(AC.IPD, BC.ALD,
            strategy = strategy_mim(formula = lin_form,
                                    family = binomial(link = "logit")),
            scale = "log_relative_risk")
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 1.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 1001 / 4000 [ 25%]  (Sampling)
#> Chain 1: Iteration: 1400 / 4000 [ 35%]  (Sampling)
#> Chain 1: Iteration: 1800 / 4000 [ 45%]  (Sampling)
#> Chain 1: Iteration: 2200 / 4000 [ 55%]  (Sampling)
#> Chain 1: Iteration: 2600 / 4000 [ 65%]  (Sampling)
#> Chain 1: Iteration: 3000 / 4000 [ 75%]  (Sampling)
#> Chain 1: Iteration: 3400 / 4000 [ 85%]  (Sampling)
#> Chain 1: Iteration: 3800 / 4000 [ 95%]  (Sampling)
#> Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.2 seconds (Warm-up)
#> Chain 1:                0.6 seconds (Sampling)
#> Chain 1:                0.8 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1.2e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
#> Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
#> Chain 2: Iteration: 1001 / 4000 [ 25%]  (Sampling)
#> Chain 2: Iteration: 1400 / 4000 [ 35%]  (Sampling)
#> Chain 2: Iteration: 1800 / 4000 [ 45%]  (Sampling)
#> Chain 2: Iteration: 2200 / 4000 [ 55%]  (Sampling)
#> Chain 2: Iteration: 2600 / 4000 [ 65%]  (Sampling)
#> Chain 2: Iteration: 3000 / 4000 [ 75%]  (Sampling)
#> Chain 2: Iteration: 3400 / 4000 [ 85%]  (Sampling)
#> Chain 2: Iteration: 3800 / 4000 [ 95%]  (Sampling)
#> Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.19 seconds (Warm-up)
#> Chain 2:                0.592 seconds (Sampling)
#> Chain 2:                0.782 seconds (Total)
#> Chain 2:

outstandR_mim_lrr
#> Object of class 'outstandR' 
#> Scale: log_relative_risk 
#> Common treatment: C 
#> Individual patient data study: AC 
#> Aggregate level data study: BC 
#> Confidence interval level: 0.95 
#> 
#> # A tibble: 3 × 5
#>   Treatments Estimate Std.Error lower.0.95 upper.0.95
#>   <chr>         <dbl>     <dbl>      <dbl>      <dbl>
#> 1 AB           -0.445    0.0734     -0.976     0.0858
#> 2 AC           -1.13     0.0370     -1.51     -0.754 
#> 3 BC           -0.686    0.0364     -1.06     -0.312

Model comparison

$AC$ effect in $BC$ population

The true $AC$ effect on the log OR scale in the $BC$ (aggregate trial data) population is $\beta_{trt}^{AC} + \beta_{EM} (\bar{X}^{AC}_1 + \bar{X}_2^{AC})$ . Calculated by

d_AC_true <- b_trt + b_EM * (BC.ALD$mean.X1 + BC.ALD$mean.X2)

The naive approach is to just convert directly from one population to another, ignoring the imbalance in effect modifiers.

d_AC_naive <- 
  AC.IPD |> 
  group_by(trt) |> 
  summarise(y_sum = sum(y), y_bar = mean(y), n = n()) |> 
  tidyr::pivot_wider(names_from = trt,
                     values_from = c(y_sum, y_bar, n)) |> 
  mutate(d_AC =
           log(y_bar_1/(1-y_bar_1)) - log(y_bar_0/(1-y_bar_0)),
         var_AC =
           1/(n_0-y_sum_0) + 1/y_sum_0 + 1/(n_1-y_sum_1) + 1/y_sum_1)

$AB$ effect in $BC$ population

This is the indirect effect. The true $AB$ effect in the $BC$ population is $\beta_{trt}^{AC} - \beta_{trt}^{BC}$ .

Following the simulation study in Remiro et al (2020) these cancel out and the true effect is zero.

The naive comparison calculating $AB$ effect in the $BC$ population is

d_BC <-
  with(BC.ALD, log(y.B.bar/(1-y.B.bar)) - log(y.C.bar/(1-y.C.bar)))

d_AB_naive <- d_AC_naive$d_AC - d_BC

var.d.BC <- with(BC.ALD, 1/y.B.sum + 1/(N.B - y.B.sum) + 1/y.C.sum + 1/(N.C - y.C.sum))
var.d.AB.naive <- d_AC_naive$var_AC + var.d.BC

Of course, the $BC$ contrast is calculated directly.

Results

We will combine all outputs for a log-odds ratio table of all contrasts and methods in the $BC$ population.

res_tab <- 
  data.frame(
    d_true = c(d_AB_true, d_AC_true, d_BC),
    d_naive = c(d_AB_naive, d_AC_naive$d_AC, d_BC),
    `MAIC` = unlist(outstandR_maic$contrasts),
    `STC` = unlist(outstandR_stc$contrasts),
    `Gcomp ML` = unlist(outstandR_gcomp_ml$contrasts),
    `Gcomp Bayes` = unlist(outstandR_gcomp_stan$contrasts),
    `MIM` = unlist(outstandR_mim$contrasts))

res_tab_var <- 
  data.frame(
    d_true = c(NA, NA, NA),
    d_naive = c(var.d.AB.naive, d_AC_naive$var_AC, var.d.BC),
    `MAIC` = unlist(outstandR_maic$variances),
    `STC` = unlist(outstandR_stc$variances),
    `Gcomp ML` = unlist(outstandR_gcomp_ml$variances),
    `Gcomp Bayes` = unlist(outstandR_gcomp_stan$variances),
    `MIM` = unlist(outstandR_mim$variances))

knitr::kable(res_tab)

	d_true	d_naive	MAIC	STC	Gcomp.ML	Gcomp.Bayes	MIM
AB	0.000000	-0.9477894	-1.001163	0.3656568	-0.9660545	-0.9795259	-0.959866
AC	-1.294548	-2.0464017	-2.099775	-0.7329555	-2.0646668	-2.0781382	-2.058478
BC	-1.098612	-1.0986123	-1.098612	-1.0986123	-1.0986123	-1.0986123	-1.098612

d_AB_true_lrr <- 0
d_AC_true_lrr <- log(plogis(d_A_true) / plogis(d_C_true))

On log relative risk scale the table is

res_tab_lrr <- 
  data.frame(
    # d_true = c(d_AB_true_lrr, d_AC_true_lrr, d_BC_lrr),
    # d_naive = c(d_AB_naive_lrr, d_AC_naive_lrr$d_AC, d_BC_lrr),
    `MAIC` = unlist(outstandR_maic_lrr$contrasts),
    `STC` = unlist(outstandR_stc_lrr$contrasts),
    `Gcomp ML` = unlist(outstandR_gcomp_ml_lrr$contrasts),
    `Gcomp Bayes` = unlist(outstandR_gcomp_stan_lrr$contrasts),
    `MIM` = unlist(outstandR_mim_lrr$contrasts))

res_tab_var_lrr <- 
  data.frame(
    # d_true = c(NA, NA, NA),
    # d_naive = c(var.d.AB.naive, d_AC_naive$var_AC, var.d.BC),
    `MAIC` = unlist(outstandR_maic_lrr$variances),
    `STC` = unlist(outstandR_stc_lrr$variances),
    `Gcomp ML` = unlist(outstandR_gcomp_ml_lrr$variances),
    `Gcomp Bayes` = unlist(outstandR_gcomp_stan_lrr$variances),
    `MIM` = unlist(outstandR_mim_lrr$variances))

knitr::kable(res_tab_lrr)

	MAIC	STC	Gcomp.ML	Gcomp.Bayes	MIM
AB	-0.4582824	0.2117284	-0.4626977	-0.4574343	-0.4450257
AC	-1.1439389	-0.4739281	-1.1483542	-1.1430909	-1.1306822
BC	-0.6856565	-0.6856565	-0.6856565	-0.6856565	-0.6856565

Plots

The same output in forest plots is

library(ggplot2)

var_dat <- 
  t(res_tab_var) |> 
  as.data.frame() |> 
  tibble::rownames_to_column("type") |>
  reshape2::melt(variable.name = "Comparison",
                 value.name = "var")

plotdat <- 
  t(res_tab) |> 
  as.data.frame() |> 
  tibble::rownames_to_column("type") |>
  reshape2::melt(variable.name = "Comparison",
                 value.name = "Estimate") |> 
  mutate(id = 1:n(),
         type = as.factor(type)) |> 
  merge(var_dat) |> 
  mutate(lo = Estimate + qnorm(0.025) * sqrt(var),
         hi = Estimate + qnorm(0.975) * sqrt(var))

ggplot(aes(x = Estimate, y = id, col = type), data = plotdat) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(aes(y = id, yend = id, x = lo, xend = hi), na.rm = TRUE) +
  xlab("Estimate (Log OR)") +
  facet_grid(Comparison~., switch = "y", scales = "free_y", space = "free_y") +
  scale_y_reverse(name = "Comparison in BC population",
                  breaks = NULL, expand = c(0, 0.6))

library(ggplot2)

var_dat <- 
  t(res_tab_var_lrr) |> 
  as.data.frame() |> 
  tibble::rownames_to_column("type") |>
  reshape2::melt(variable.name = "Comparison",
                 value.name = "var")

plotdat <- 
  t(res_tab_lrr) |> 
  as.data.frame() |> 
  tibble::rownames_to_column("type") |>
  reshape2::melt(variable.name = "Comparison",
                 value.name = "Estimate") |> 
  mutate(id = 1:n(),
         type = as.factor(type)) |> 
  merge(var_dat) |> 
  mutate(lo = Estimate + qnorm(0.025) * sqrt(var),
         hi = Estimate + qnorm(0.975) * sqrt(var))

ggplot(aes(x = Estimate, y = id, col = type), data = plotdat) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(aes(y = id, yend = id, x = lo, xend = hi), na.rm = TRUE) +
  xlab("Estimate (Log RR)") +
  facet_grid(Comparison~., switch = "y", scales = "free_y", space = "free_y") +
  scale_y_reverse(name = "Comparison in BC population",
                  breaks = NULL, expand = c(0, 0.6))

Introduction

General problem

Example analysis

Data

AC.IPD: Individual patient data

BC.ALD: Aggregate-level data

Output statistics

Model fitting in R

MAIC

Outcome scale

Simulated treatment comparison

Parametric G-computation with maximum-likelihood estimation

Bayesian G-computation with MCMC

Multiple imputation marginalisation

Model comparison

ACAC effect in BCBC population

ABAB effect in BCBC population

Results

Plots

`AC.IPD`: Individual patient data

`BC.ALD`: Aggregate-level data

$AC$ effect in $BC$ population

$AB$ effect in $BC$ population