Introduction

For survival outcomes, MECs can be used as a control to compare the effect of a control treatment against the effect of an observed experimental treatment on overall survival in patients.

We will implement MECs using the ‘psc’ package which specifically allows for the comparison of patient groups treated with an experimental treatment against a counterfactual model (CFM).

Load necessary packages:

library(psc)
library(survival)
library(ggpubr)

An example

In this example, we will look at survival outcomes in patients with pancreatic cancer from the ESPAC (European Study for Pancreatic Cancer)-4 trial and compare how the survival of patients differs with two different treatments: GEM Vs. GEMCAP…

Load the model that will act as the control treatment (CFM). This model is a flexible parametric model.

load("flsm.R")

Now that we have our control treatment, monotherapy gemcitabine (GEM), we will load the ESPAC-4 dataset. ESPAC-4 consists of patients that have been treated with the experimental treatment, adjuvant gemcitabine and capecitabine, GEMCAP.

Our aim is to compare the patients treated with GEMCAP in the ESPAC-4 dataset against the GEM model to determine which of the treatments is more effective.

load("espac4gemcap.R")

Fit counterfactual model :)

We can turn the flexible parametric model into a CFM using the pscCFM() function:

cfm <- pscCFM(flsm)

It is possible to visualise and summarise the data that was used to fit this model with the built in functions within pscCFM():

cfm$datasumm$summ_Table

ggarrange(plotlist = cfm$datavis, ncol = 2)

$`1`

$`2`

attr(,"class")
[1] "list"      "ggarrange"

Make comparison!

The CFM can be compared against the ESPAC-4 data cohort. The comparison will be carried out using the pscfcit() function.

psc <- pscfit(cfm, espac4_gemcap)

The plot below visualises the effect of each treatment on ESPAC-4 patients. The pink line represents the model’s predicted survival estimates (if the ESPAC-4 patients had been treated with GEM) and the orange line represents the data cohort’s observed survival estimates.

plot(psc)

The summary of the model can be obtained using the generic summary() function to show the posterior density. The posterior estimates of the deviance information criterion (DIC) and the efficacy parameter, β, are calculated and shown.

β is a measurement of the distance between the observed data and the model estimate. The 2.5 and 97.5% quantiles are also given. The odds ratio comparing the experimental treatment to the control can be calculated as exp(β), in this case exp(-0.466452).

summary(psc)

Summary: 
 
362 observations selected from the data cohort for comparison 
CFM of type flexsurvreg identified  
linear predictor succesfully obtained with median: 
 trt: 0.812
Average expected response: 
 trt: 28.337
Average observed response: 31.67 

Counterfactual Model (CFM): 
A model of class 'flexsurvreg' 
 Fit with 5 internal knots

Formula: 
Surv(time, cen) ~ LymphN + ResecM + Diff_Status + PostOpCA199
<environment: 0x0000016a1fa155b0>

Call:
 CFM model + beta

Coefficients:
      median     2.5%       97.5%      Pr(x<0)    Pr(x>0)  
beta    -0.4627    -0.5988    -0.3378     1.0000     0.0000
DIC   1320.7933  1297.5917  1355.0851         NA         NA

We can extract the posterior distribution of the ‘psc’ and can view its autocorrelation and trace

posterior <- psc$posterior
head(posterior[1:3,])

     gamma0   gamma1   gamma2    gamma3    gamma4     gamma5     gamma6
1 -11.38080 3.835982 1.291162 -1.317656 1.1190182 -0.8876277  0.3372484
2 -13.30464 4.898178 4.394070 -5.271692 1.6488526  0.7120773 -0.6273199
3 -10.92023 3.903896 2.945046 -2.635881 0.5773368 -0.8796604  0.8259487
    LymphN1     ResecM1 Diff_Status1 Diff_Status2 PostOpCA199       beta
1 0.4876152  0.18053219   -0.4160534   -0.5897823   0.2671471 -0.4552224
2 0.6315325 -0.08449907   -0.5292838   -0.2732867   0.2633989 -0.4552224
3 0.7238781  0.08044458   -0.3553484   -0.5314689   0.2202701 -0.4552224
       DIC
1       NA
2 1310.674
3 1314.744

# autocorrelation
acf(posterior$beta)

# trace
plot(posterior$beta, type = 's')