I am currently using the CausalImpact package for some research and in this context I need to know and be able to explain, how the posterior tail-area probability is calculated in order to reproduce that value for validation purposes. Does anyone know, how that value can be reproduced given the data and the estimation series provided by the model? Thanks in advance!
Asked
Active
Viewed 362 times
0
-
Have you read [the associated research paper](https://storage.googleapis.com/pub-tools-public-publication-data/pdf/41854.pdf)? – merv Feb 14 '19 at 02:57
-
I did. However, if you can tell me where exactly they talk about the tail-area probability, I'd appreciate it! In section 2.3 they talk about inference, yet I couldn't find (or just did not understand it when I saw it) where and how the p-value is calculated – Tiko Feb 14 '19 at 10:06
1 Answers
0
I've never used this library, but skimming through the code, it appears that they compute the quantiles (alpha/2
, 1-alpha/2
) of the samples from the posterior predictive distribution.
From the relevant section of code (Apache v2.0 License)
ComputeCumulativePredictions <- function(y.samples, point.pred, y,
post.period.begin, alpha = 0.05) {
# Computes summary statistics for the cumulative posterior predictions over
# the unobserved data points in the post-intervention period.
#
# Args:
# y.samples: Matrix of simulated response trajectories, as returned
# by \code{ComputeResponseTrajectories()}.
# point.pred: Data frame of point predictions, as returned by
# \code{ComputePointPredictions()}.
# y: Actual observed response, from the beginning of the
# pre-period to the end of the observed period.
# post.period.begin: Index of the first data point of the post-period.
# alpha: The resulting coverage of the posterior intervals will
# be \code{1 - alpha}.
#
# Returns:
# data frame with 3 columns:
# cum.pred: posterior predictive expectation
# cum.pred.lower: lower limit of a \code{(1 - alpha)*100}% interval
# cum.pred.upper: upper limit
... # [Computing the cum.pred.mean]
prob.lower <- alpha / 2 # e.g., 0.025 when alpha = 0.05
prob.upper <- 1 - alpha / 2 # e.g., 0.975 when alpha = 0.05
cum.pred.lower.post <- as.numeric(t(apply(y.samples.cum.post, 2, quantile,
prob.lower)))
cum.pred.upper.post <- as.numeric(t(apply(y.samples.cum.post, 2, quantile,
prob.upper)))
cum.pred.lower <- c(cum.pred.lower.pre, cum.pred.lower.post)
cum.pred.upper <- c(cum.pred.upper.pre, cum.pred.upper.post)
# Put cumulative prediction together
cum.pred <- data.frame(cum.pred = cum.pred.mean,
cum.pred.lower, cum.pred.upper)
return(cum.pred)
}

merv
- 67,214
- 13
- 180
- 245