# Q: R - CausalImpact | order of control time series [closed]

My colleagues and I have been using Google's R-package CausalImpact for a while but recently discovered something that we can't quite explain. Depending on the order of including the different control time series into the variable "data" CausalImpact generates different outputs.

Example:

Using a slightly different version of the example on CausalImpact's Github we first create 2 control time series and a response variable.

library(CausalImpact)
set.seed(1)
x1 <- 100 + arima.sim(model = list(ar = 0.9), n = 100)
x2 <- 99  + arima.sim(model = list(ar = 0.8), n = 100)
y <- 1.2 * x1 + rnorm(100)
y[71:100] <- y[71:100] + 10


Define pre- and post-period.

pre.period <- c(1, 70)
post.period <- c(71, 100)


And finally

# order 1
data <- cbind(y, x1, x2)
set.seed(1)
impact <- CausalImpact(data, pre.period, post.period)
summary(impact)

# order 2
data <- cbind(y, x2, x1)
set.seed(1)
impact <- CausalImpact(data, pre.period, post.period)
summary(impact)


These generate (in this case slightly) different outputs, even though (to our understanding) it shouldn't. Does CausalImpact weigh the "first" control time series different to the "second"? Any help would be much appreciated!

## closed as off-topic by gung♦, mdewey, Michael Chernick, kjetil b halvorsen, FirebugJun 13 '17 at 20:10

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – gung, mdewey, Michael Chernick, kjetil b halvorsen, Firebug
If this question can be reworded to fit the rules in the help center, please edit the question.

• it might be more appropriate to file an issue directly on GitHub... – Antoine Jun 8 '17 at 8:26
• This appears to be a question about how software works, not about statistics per se. As such it is off topic here. If you have a question about time series analysis, please edit to clarify. – gung Jun 12 '17 at 20:14

Even though you've set the seed to be the same for both runs, changing the order of the predictors means that the random draws for each predictor within BSTS are going to be different between the runs.

Note that when Causal Impact (CI) calls BSTS to generate the model, it resets the random seed to 1. So it doesn't matter what the seed is going into the CI call.

This seed reset prevents you from seeing the random variation that should be appearing in the results from each call of Causal Impact.

By default, CI calls BSTS with a parameter of 1000 iterations. If you crank up the iterations in the call to BSTS, the results converge. Half million iterations may be overkill, but does the trick:

# order 1
data1 <- cbind(y, x1, x2)
set.seed(1)
impact1 <- CausalImpact(data1, pre.period, post.period, model.args = list(niter = 500000))
summary(impact1)

# order 2
data2 <- cbind(y, x2, x1)
set.seed(1)
impact2 <- CausalImpact(data2, pre.period, post.period,  model.args = list(niter = 500000))
summary(impact2)


As an aside, you can see the variation in BSTS results directly if you build the BSTS model externally and look at the variation in the estimated coefficient.

results_data1 <- vector(mode = "numeric", length = 0)
results_data2 <- vector(mode = "numeric", length = 0)
for(i in 1:1000) {

external.model1 <- bsts(y ~ .,
data = data1,
state.specification = ss,
expected.model.size = 3,
expected.r2 = 0.9,
prior.df = 50,
niter = 1000,
ping = 0,
model.options = BstsOptions(save.prediction.errors      = TRUE)
)
results_data1  <- append(results_data1,     mean(external.model1$coefficients[,2]) , after = length(results_data1)) external.model2 <- bsts(y ~ ., data = data2, state.specification = ss, expected.model.size = 3, expected.r2 = 0.9, prior.df = 50, niter = 1000, ping = 0, model.options = BstsOptions(save.prediction.errors = TRUE) ) results_data2 <- append(results_data2, mean(external.model2$coefficients[,3]) , after = length(results_data2))

print(i)
}

summary(results_data1)


| Min. | 1st Qu.| Median | Mean | 3rd Qu. | Max. |
-0.0496 | 0.3055 | 0.3156 | 0.3140 | 0.3260 | 0.3591 |

 summary(results_data2)


| Min. | 1st Qu.| Median | Mean | 3rd Qu. | Max. |
0.0260 | 0.3052 | 0.3161 | 0.3146 | 0.3260 | 0.3698

results_dif <- results_data1 - results_data2
summary(results_dif)


| Min. | 1st Qu.| Median | Mean | 3rd Qu. | Max. |
|-0.3633710 | -0.0143417 | -0.0009655 | -0.0006387 | 0.0137142 | 0.3039858

Over a 1000 calls of a 1000 iterations, you can see the sizeable variation in the coefficient estimated from any given 1000 iteration call, and a small average difference remains.

t.test(results_data1,results_data2)


Welch Two Sample t-test
data: results_data1 and results_data2
t = -0.62233, df = 1975.4, p-value = 0.5338
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.002651600 0.001374129
sample estimates:
mean of x mean of y
0.3139709 0.3146096

But the differences between the runs with X1 first and the runs with X2 first are not significant.