Casual Impact Notes
The follorwing is a collection of notes based on Googles reseach and practical application of Causal Inference. This demonstration is conducted in an R Notebook using the R package CausalImpact created by Google for this purpose.
Installing the package
To install CausalImpact, type the following commands into an R session:
install.packages("devtools")
library(devtools)
devtools::install_github("google/CausalImpact")Note: You can now use
install.packages("CausalImpact")library(bsts) # required for Causal Impact
library(CausalImpact)Once installed, the package can be loaded in a given R session using:
library(CausalImpact)Creating an example dataset
To illustrate how the package works, we create a simple toy dataset. It consists of a response variable y and a predictor x1. Note that in practice, we’d strive for including many more predictor variables and let the model choose an appropriate subset. The example data has 100 observations. We create an intervention effect by lifting the response variable by 10 units after timepoint 71.
set.seed(1)
x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 100)
y <- 1.2 * x1 + rnorm(100)
y[71:100] <- y[71:100] + 10
data <- cbind(y, x1)We now have a simple matrix with 100 rows and two columns:
dim(data)[1] 100 2[1] 100 2
head(data) y x1
[1,] 105.2950 88.21513
[2,] 105.8943 88.48415
[3,] 106.6209 87.87684
[4,] 106.1572 86.77954
[5,] 101.2812 84.62243
[6,] 101.4484 84.60650We can visualize the generated data using:
matplot(data, type = "l")
Running an analysis
To estimate a causal effect, we begin by specifying which period in the data should be used for training the model (pre-intervention period) and which period for computing a counterfactual prediction (post-intervention period).
pre.period <- c(1, 70)
post.period <- c(71, 100)This says that time points 1 … 70 will be used for training, and time points 71 … 100 will be used for computing predictions. Alternatively, we could specify the periods in terms of dates or time points; see Section 5 for an example.
To perform inference, we run the analysis using:
impact <- CausalImpact(data, pre.period, post.period)This instructs the package to assemble a structural time-series model, perform posterior inference, and compute estimates of the causal effect. The return value is a CausalImpact object.
Plotting the results
The easiest way of visualizing the results is to use the plot() function that is part of the package:
plot(impact)
By default, the plot contains three panels. The first panel shows the data and a counterfactual prediction for the post-treatment period. The second panel shows the difference between observed data and counterfactual predictions. This is the pointwise causal effect, as estimated by the model. The third panel adds up the pointwise contributions from the second panel, resulting in a plot of the cumulative effect of the intervention.
Remember, once again, that all of the above inferences depend critically on the assumption that the covariates were not themselves affected by the intervention. The model also assumes that the relationship between covariates and treated time series, as established during the pre-period, remains stable throughout the post-period.
- Working with dates and times
It is often more natural to feed a time-series object into CausalImpact() rather than a data frame. For example, we might create a data variable as follows:
time.points <- seq.Date(as.Date("2014-01-01"), by = 1, length.out = 100)
data <- zoo(cbind(y, x1), time.points)
head(data) y x1
2014-01-01 105.2950 88.21513
2014-01-02 105.8943 88.48415
2014-01-03 106.6209 87.87684
2014-01-04 106.1572 86.77954
2014-01-05 101.2812 84.62243
2014-01-06 101.4484 84.60650We can now specify the pre-period and the post-period in terms of time points rather than indices:
pre.period <- as.Date(c("2014-01-01", "2014-03-11"))
post.period <- as.Date(c("2014-03-12", "2014-04-10"))As a result, the x-axis of the plot shows time points instead of indices:
impact <- CausalImpact(data, pre.period, post.period)
plot(impact)
Printing a summary table
To obtain a numerical summary of the analysis, we use:
summary(impact)Posterior inference {CausalImpact}
Average Cumulative
Actual 117 3511
Prediction (s.d.) 107 (0.37) 3196 (11.02)
95% CI [106, 107] [3174, 3217]
Absolute effect (s.d.) 11 (0.37) 315 (11.02)
95% CI [9.8, 11] [294.4, 338]
Relative effect (s.d.) 9.9% (0.34%) 9.9% (0.34%)
95% CI [9.2%, 11%] [9.2%, 11%]
Posterior tail-area probability p: 0.001
Posterior prob. of a causal effect: 99.9%
For more details, type: summary(impact, "report")The Average column talks about the average (across time) during the post-intervention period (in the example: time points 71 through 100). The Cumulative column sums up individual time points, which is a useful perspective if the response variable represents a flow quantity (such as queries, clicks, visits, installs, sales, or revenue) rather than a stock quantity (such as number of users or stock price).
In the example, the estimated average causal effect of treatment was 11 (rounded to a whole number; for full precision see impact$summary). This is because we observed an average value of 99 but would have expected an average value of only 89. The 95% posterior interval of the average effect is [9.8, 11]. Since this excludes 0, we (correctly) conclude that the intervention had a causal effect on the response variable. Since we generated the data ourselves, we know that we injected a true effect of 10, and so the model accurately recovered ground truth. One reason for this is that we ensured, by design, that the covariate x1 was not itself affected by the intervention. In practice, we must always reason whether this assumption is justified.
For additional guidance about the correct interpretation of the summary table, the package provides a verbal interpretation, which we can print using:
summary(impact, "report")Analysis report {CausalImpact}
During the post-intervention period, the response variable had an average value of approx. 117.05. By contrast, in the absence of an intervention, we would have expected an average response of 106.54. The 95% interval of this counterfactual prediction is [105.80, 107.24]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 10.51 with a 95% interval of [9.81, 11.25]. 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 3.51K. By contrast, had the intervention not taken place, we would have expected a sum of 3.20K. The 95% interval of this prediction is [3.17K, 3.22K].
The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +10%. The 95% interval of this percentage is [+9%, +11%].
This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (10.51) to the original goal of the underlying intervention.
The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.001). This means the causal effect can be considered statistically significant. The individual numbers in the table, at full precision, can be accessed using:
impact$summarySee below for tips on how to use these commands with knitr / R Markdown.
Adjusting the model
So far, we’ve simply let the package decide how to construct a time-series model for the available data. However, there are several options that allow us to gain a little more control over this process. These options are passed into model.args as individual list elements, for example:
impact <- CausalImpact(..., model.args = list(niter = 5000, nseasons = 7))Available options
niterNumber of MCMC samples to draw. More samples lead to more accurate inferences. Defaults to 1000.
standardize.dataWhether to standardize all columns of the data before fitting the model. This is equivalent to an empirical Bayes approach to setting the priors. It ensures that results are invariant to linear transformations of the data. Defaults to TRUE.
prior.level.sdPrior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.
nseasonsPeriod of the seasonal components. In order to include a seasonal component, set this to a whole number greater than 1. For example, if the data represent daily observations, use 7 for a day-of-week component. This interface currently only supports up to one seasonal component. To specify multiple seasonal components, use bsts to specify the model directly, then pass the fitted model in as bsts.model. Defaults to 1, which means no seasonal component is used.
season.durationDuration of each season, i.e., number of data points each season spans. Defaults to 1. For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
dynamic.regressionWhether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
Using a custom model
Instead of using the default model constructed by the CausalImpact package, we can use the bsts package to specify our own model. This provides the greatest degree of flexibility.
Before constructing a custom model, we set the observed data in the post-treatment period to NA, reflecting the fact that the counterfactual response is unobserved after the intervention. We keep a copy of the actual observed response in the variable post.period.response.
post.period <- c(71, 100)
post.period.response <- y[post.period[1] : post.period[2]]
y[post.period[1] : post.period[2]] <- NAWe next set up and estimate a time-series model using the bsts package. Here is a simple example:
ss <- AddLocalLevel(list(), y)
bsts.model <- bsts(y ~ x1, ss, niter = 1000)=-=-=-=-= Iteration 0 Mon May 8 14:40:52 2017 =-=-=-=-=
=-=-=-=-= Iteration 100 Mon May 8 14:40:52 2017 =-=-=-=-=
=-=-=-=-= Iteration 200 Mon May 8 14:40:52 2017 =-=-=-=-=
=-=-=-=-= Iteration 300 Mon May 8 14:40:52 2017 =-=-=-=-=
=-=-=-=-= Iteration 400 Mon May 8 14:40:53 2017 =-=-=-=-=
=-=-=-=-= Iteration 500 Mon May 8 14:40:53 2017 =-=-=-=-=
=-=-=-=-= Iteration 600 Mon May 8 14:40:53 2017 =-=-=-=-=
=-=-=-=-= Iteration 700 Mon May 8 14:40:53 2017 =-=-=-=-=
=-=-=-=-= Iteration 800 Mon May 8 14:40:53 2017 =-=-=-=-=
=-=-=-=-= Iteration 900 Mon May 8 14:40:53 2017 =-=-=-=-=Finally, we call CausalImpact Instead of providing input data, we simply pass in the fitted model object (bsts.model). We also need to provide the actual observed response. This is needed so that the package can compute the difference between predicted response (stored in bsts.model) and actual observed response (stored in post.period.response).
impact <- CausalImpact(bsts.model = bsts.model,
post.period.response = post.period.response)The results can be inspected in the usual way:
plot(impact)
summary(impact)Posterior inference {CausalImpact}
Average Cumulative
Actual 117 3511
Prediction (s.d.) 106 (0.49) 3190 (14.59)
95% CI [105, 107] [3162, 3218]
Absolute effect (s.d.) 11 (0.49) 321 (14.59)
95% CI [9.8, 12] [293.4, 350]
Relative effect (s.d.) 10% (0.46%) 10% (0.46%)
95% CI [9.2%, 11%] [9.2%, 11%]
Posterior tail-area probability p: 0.00105
Posterior prob. of a causal effect: 99.8954%
For more details, type: summary(impact, "report")summary(impact, "report")Analysis report {CausalImpact}
During the post-intervention period, the response variable had an average value of approx. 117.05. By contrast, in the absence of an intervention, we would have expected an average response of 106.34. The 95% interval of this counterfactual prediction is [105.39, 107.27]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 10.71 with a 95% interval of [9.78, 11.65]. 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 3.51K. By contrast, had the intervention not taken place, we would have expected a sum of 3.19K. The 95% interval of this prediction is [3.16K, 3.22K].
The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +10%. The 95% interval of this percentage is [+9%, +11%].
This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (10.71) to the original goal of the underlying intervention.
The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.001). This means the causal effect can be considered statistically significant.