Causal inference: a simulated clinical trial

1 Data simulation

First, let’s create a simulated dataset of a clinical trial, in which we want to evaluate the treatment effect of a new drug on a target outcome. A sample of 20 patients have been equally assigned into the Treatment group (n=10) and Placebo group (n=10). The target outcome is daily measured for each patient during 30 days before the treatment and 30 days during the treatment.

library(tidyverse)

df = tibble(Subject = NA, Time = NA, Outcome = NA)

for (i in 1:10){
  set.seed(i+100)
  
  sim_df = mgcv::gamSim(eg=1,n=30,dist="normal",scale=1,verbose=F)
  
  temp_df = tibble(Subject = paste("Patient_",i,sep = ""),
                   Time = c(1:60),
                   Outcome = c(sim_df$x1*sim_df$f1,
                               sim_df$x2*sim_df$f2)*10)
  
  df = bind_rows(df,temp_df)%>%na.omit()
  rm(sim_df,temp_df,i)
}

df2 = tibble(Subject = NA, Time = NA, Outcome = NA)

for (i in 1:10){
  set.seed(i+200)
  
  sim_df = mgcv::gamSim(eg=1,n=60,dist="normal",scale=1,verbose=F)
  
  temp_df = tibble(Subject = paste("Patient_",i,sep = ""),
                   Time = c(1:60),
                   Outcome = (sim_df$x2*sim_df$f2)*10)
  
  df2 = bind_rows(df2,temp_df)%>%na.omit()
  rm(sim_df,temp_df,i)
}

For each group, we have 10 data sequences of the target outcome, indexed from 1 to 60. In the treatment group, the changing patterns look like this:

df %>% mutate(Phase = Time<31)%>%
  ggplot(aes(x=Time,y=Outcome))+
  geom_rect(xmin=0,xmax=30,ymin=0,ymax=70,fill="red",alpha=0.005,color=NA)+
  geom_rect(xmin=31,xmax=60,ymin=0,ymax=70,fill="green",alpha=0.005,color=NA)+
  geom_line(aes(col=Phase),
            size=0.8)+
  geom_point(aes(col=Phase),
             size=2)+
  facet_wrap(~ Subject, ncol=2)+
  geom_vline(xintercept = 31)+
  scale_colour_manual(values = c("green4","red"))+
  theme_bw(8)

In contrast, the Placebo group did not show any change in outcome during the last 30 days

df2 %>% mutate(Phase = Time<31)%>%
  ggplot(aes(x=Time,y=Outcome))+
  geom_rect(xmin=0,xmax=30,ymin=0,ymax=70,fill="red",alpha=0.005,color=NA)+
  geom_rect(xmin=31,xmax=60,ymin=0,ymax=70,fill="green",alpha=0.005,color=NA)+
  geom_line(aes(col=Phase),
            size=1)+
  geom_point(aes(col=Phase),
             size=2)+
  facet_wrap(~ Subject, ncol=2)+
  geom_vline(xintercept = 31)+
  scale_colour_manual(values = c("green4","red"))+
  theme_bw()

2 Application of the CausalImpact algorithm

We will apply the Causal Impact algorithm on each outcome sequences in two datasets, then compile the results into a meta-data

library(CausalImpact)

model_func = function(data=x,
                      pre=c(1,30),
                      post=c(31,60),...){
  cimp = CausalImpact(data[c(1:60),3],c(1,30), c(31,60), ...)%>%
    .$summary%>%
    mutate(Type=rownames(.))
  return(cimp)
}

df %>% split(.$Subject)%>%
  map_df(~model_func(data=., 
                     alpha=0.01)) -> meta_df

meta_df$Patient = rep(paste("Treat",1:10),each=2)
meta_df$Group = "Treatment"

df2 %>% split(.$Subject)%>% 
  map_df(~model_func(data=., 
                     alpha=0.01)) -> meta_df2

meta_df2$Patient = rep(paste("Placebo_",1:10),each=2)
meta_df2$Group = "Placebo"

comb_df = bind_rows(meta_df,meta_df2)

Now, we can evaluate the absolute treatment effects between two groups

At individual level:

comb_df%>%filter(Type=="Average")%>%
  ggplot(aes(x=Patient, col=Group))+
  geom_errorbar(aes(ymin= AbsEffect.lower, 
                    ymax= AbsEffect.upper), 
                width=0.2, size=1)+
  geom_point(size=3,aes(y=AbsEffect))+
  geom_hline(yintercept = 0)+
  theme_bw()+
  theme(axis.text.x = element_text(angle =45,hjust=1,vjust=1))

or for the pooled sample:

comb_df%>%filter(Type=="Average")%>%
  ggplot(aes(x=Group, 
             y=AbsEffect,
             fill=Group))+
  geom_boxplot()+
  geom_hline(yintercept = 0)+
  theme_bw()+
  coord_flip()+
  theme(axis.text.x = element_text(angle =45,
                                   hjust=1,
                                   vjust=1))

What did happen ?

For each patient, the CausalImpact algorithm fitted a Bayesian structural timeseries model using the first 30 data points, then used this model to forecast the tendency of the outcome during the next 30 days, if the patients would have not received any treatment. It’s called a “hypothetical baseline”, or an approximation of the counterfactual effect which cannot be seen in reality. A causal inference was performed by measuring the pointwise differences between the true observed values and the hypothetical baseline. A statistical null hypothesis testing was done to verify whether the Causal impact is not null.

df%>%filter(Subject=="Patient_1")%>% .[,3]%>%
  CausalImpact(., c(1,30), c(31,60)) -> mod

summary(mod)

## Posterior inference {CausalImpact}
## 
##                          Average          Cumulative    
## Actual                   9.5              284.5         
## Prediction (s.d.)        19 (4.8)         562 (142.8)   
## 95% CI                   [9.3, 29]        [279.2, 861]  
##                                                         
## Absolute effect (s.d.)   -9.3 (4.8)       -277.7 (142.8)
## 95% CI                   [-19, 0.18]      [-576, 5.39]  
##                                                         
## Relative effect (s.d.)   -49% (25%)       -49% (25%)    
## 95% CI                   [-102%, 0.96%]   [-102%, 0.96%]
## 
## Posterior tail-area probability p:   0.029
## Posterior prob. of a causal effect:  97.1%
## 
## For more details, type: summary(impact, "report")

summary(mod, "report")

## Analysis report {CausalImpact}
## 
## 
## During the post-intervention period, the response variable had an average value of approx. 9.48. In the absence of an intervention, we would have expected an average response of 18.74. The 95% interval of this counterfactual prediction is [9.31, 28.68]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is -9.26 with a 95% interval of [-19.20, 0.18]. For a discussion of the significance of this effect, see below.
## 
## Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of 284.54. Had the intervention not taken place, we would have expected a sum of 562.23. The 95% interval of this prediction is [279.15, 860.52].
## 
## The above results are given in terms of absolute numbers. In relative terms, the response variable showed a decrease of-49%. The 95% interval of this percentage is [-102%, +1%].
## 
## This means that, although it may look as though the intervention has exerted a negative effect on the response variable when considering the intervention period as a whole, this effect is not statistically significant, and so cannot be meaningfully interpreted. The apparent effect could be the result of random fluctuations that are unrelated to the intervention. This is often the case when the intervention period is very long and includes much of the time when the effect has already worn off. It can also be the case when the intervention period is too short to distinguish the signal from the noise. Finally, failing to find a significant effect can happen when there are not enough control variables or when these variables do not correlate well with the response variable during the learning period.
## 
## The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.029). This means the causal effect can be considered statistically significant.

plot(mod)

## Warning: Removed 60 rows containing missing values (geom_path).

## Warning: Removed 120 rows containing missing values (geom_path).

3 Discussion

In our days, daily monitoring of a respiratory outcome is an easy task, with the presence of various types of portable and home-based devices such as spirometer, FOT, pulse oximeter or FENO analyzers. Such longitudinal measurement generates sequential data and we need to extract the information from this kind of data.

The CausalImpact algorithm and sequential data based causal inference present many advantages, including:

1. Possibility of measuring both treatment effect and counterfacts

2. Possibility of making causal inference at both population and individual levels, thus compatible to the clinical studies in personalized medicine,

3. The Bayesian timeserie model is able to include other confounding factors