Causal Impact using Bayesian Structural Time-Series Models

Suppose we like to investigate if an intervention in some process or experiment made an impact,causal impact which relies on bayesian structural models provide the tool to accomplish this task. We have a series which starts in january with record of the number of vehicles which are tested at a vehicle testing site. Suppose the recording of these vehicles in a database was not accurately recorded until September 25,2017. We would like to investigate whether the intervention in September made any significant impact on forecasting future numbers vehicles which will be tested.

We load the required packages at this step.

pacman::p_load(CausalImpact,bsts,rio,imputeTS)

Load the data at this step.

Data<- rio::import("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/actvsfcdatajuly311.xlsx",sheet="Sheet4")

Data1<-Data%>%filter(Date<"2017-09-25")
Data2<-Data%>%filter(Date>"2017-09-25")

Impute the missing values.

Data_ts=Data1%>%tk_zoo()
plotNA.gapsize(Data_ts)

statsNA(Data_ts)

[1] "Length of time series:"
[1] 91
[1] "-------------------------"
[1] "Number of Missing Values:"
[1] 4
[1] "-------------------------"
[1] "Percentage of Missing Values:"
[1] "4.4%"
[1] "-------------------------"
[1] "Stats for Bins"
[1] "  Bin 1 (23 values from 1 to 23) :      1 NAs (4.35%)"
[1] "  Bin 2 (23 values from 24 to 46) :      0 NAs (0%)"
[1] "  Bin 3 (23 values from 47 to 69) :      3 NAs (13%)"
[1] "  Bin 4 (22 values from 70 to 91) :      0 NAs (0%)"
[1] "-------------------------"
[1] "Longest NA gap (series of consecutive NAs)"
[1] "2 in a row"
[1] "-------------------------"
[1] "Most frequent gap size (series of consecutive NA series)"
[1] "1 NA in a row (occuring 2 times)"
[1] "-------------------------"
[1] "Gap size accounting for most NAs"
[1] "2 NA in a row (occuring 1 times, making up for overall 2 NAs)"
[1] "-------------------------"
[1] "Overview NA series"
[1] "  1 NA in a row: 2 times"
[1] "  2 NA in a row: 1 times"

plotNA.distributionBar(Data_ts)

plotNA.distribution(Data_ts)

Data$Actual[which(is.na(Data$Actual))]

[1] NA NA NA NA

#
#Uses Kalman Smoothing on structural time series models (or on the state space representation of an arima model) for imputation
Data_ts=na.kalman(Data_ts)
#plotNA.imputations(Data_ts)

Data1_ts=Data1%>%tk_zoo()
Data2_ts=Data2%>%tk_zoo()
Data_ts<-rbind(Data1_ts,Data2_ts)
Data1_ts=na.kalman(Data1_ts)

# fitting weekly data. number of seasons is 52 
ss <- AddSemilocalLinearTrend(list(), Data1_ts)
ss <- AddSeasonal(ss, Data1_ts, nseasons = 52)
modeld <- bsts(Data1_ts, state.specification = ss, niter = 1000)

=-=-=-=-= Iteration 0 Tue Jan  9 15:26:17 2018 =-=-=-=-=
=-=-=-=-= Iteration 100 Tue Jan  9 15:26:18 2018 =-=-=-=-=
=-=-=-=-= Iteration 200 Tue Jan  9 15:26:18 2018 =-=-=-=-=
=-=-=-=-= Iteration 300 Tue Jan  9 15:26:19 2018 =-=-=-=-=
=-=-=-=-= Iteration 400 Tue Jan  9 15:26:19 2018 =-=-=-=-=
=-=-=-=-= Iteration 500 Tue Jan  9 15:26:20 2018 =-=-=-=-=
=-=-=-=-= Iteration 600 Tue Jan  9 15:26:20 2018 =-=-=-=-=
=-=-=-=-= Iteration 700 Tue Jan  9 15:26:21 2018 =-=-=-=-=
=-=-=-=-= Iteration 800 Tue Jan  9 15:26:21 2018 =-=-=-=-=
=-=-=-=-= Iteration 900 Tue Jan  9 15:26:22 2018 =-=-=-=-=

burn <- SuggestBurn(0.1, modeld)
predd <- predict(modeld, horizon = 15, burn = burn, quantiles = c(.025, .975))

pred1 <- predict(modeld, horizon = 15,quantiles = c(.025, .975))
plot(predd, plot.original = 156)

pre.period <- c(1, 91)
post.period <- c(92, 106)
impact <- CausalImpact(Data$Actual,pre.period,
                       post.period)

# Print and plot results
summary(impact)

Posterior inference {CausalImpact}

                         Average      Cumulative
Actual                   52           775       
Prediction (s.d.)        30 (1.2)     446 (18.4)
95% CI                   [27, 32]     [409, 482]
                                                
Absolute effect (s.d.)   22 (1.2)     329 (18.4)
95% CI                   [20, 24]     [293, 366]
                                                
Relative effect (s.d.)   74% (4.1%)   74% (4.1%)
95% CI                   [66%, 82%]   [66%, 82%]

Posterior tail-area probability p:   0.001
Posterior prob. of a causal effect:  99.8999%

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. 51.67. By contrast, in the absence of an intervention, we would have expected an average response of 29.70. The 95% interval of this counterfactual prediction is [27.29, 32.13]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 21.96 with a 95% interval of [19.53, 24.38]. 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 775.00. By contrast, had the intervention not taken place, we would have expected a sum of 445.57. The 95% interval of this prediction is [409.34, 482.00].

The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +74%. The 95% interval of this percentage is [+66%, +82%].

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 (21.96) 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. 

plot(impact)

plot(impact, "original")

#plot(impact$model$bsts.model, "coefficients")

# Example 3
  #
  # For full flexibility, specify a custom model and pass the
  # fitted model to CausalImpact(). To run this example, run
  # the code for Example 1 first.
Data=Data%>%drop_na()
pre.period <- c(1, 87)
post.period <- c(88, 102)
y=Data$Actual
x1=Data$Date
  post.period.response <- y[post.period[1] : post.period[2]]
  y[post.period[1] : post.period[2]] <- NA
  ss <- AddLocalLevel(list(), y)
  bsts.model <- bsts(y ~ x1, ss, niter = 1000)

=-=-=-=-= Iteration 0 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 100 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 200 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 300 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 400 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 500 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 600 Tue Jan  9 15:26:25 2018 =-=-=-=-=
=-=-=-=-= Iteration 700 Tue Jan  9 15:26:26 2018 =-=-=-=-=
=-=-=-=-= Iteration 800 Tue Jan  9 15:26:26 2018 =-=-=-=-=
=-=-=-=-= Iteration 900 Tue Jan  9 15:26:26 2018 =-=-=-=-=

  impact <- CausalImpact(bsts.model = bsts.model,
                         post.period.response = post.period.response)
plot(impact)