R/synthetic-control.R

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

#' Synthetic Control via BSTS (CausalImpact)
#'
#' Builds a simple synthetic-control-style analysis using
#' \pkg{CausalImpact}/BSTS for either \code{I} or \code{C} as the outcome,
#' with treatment defined endogenously by a high level of a chosen
#' control variable.
#'
#' @param DT A \code{data.frame} or \code{data.table} containing at least:
#'   \itemize{
#'     \item \code{I}, \code{C}: outcome candidates (counts or rates).
#'     \item \code{EconCycle}, \code{PopDensity}, \code{Epidemics},
#'           \code{Climate}, \code{War}, \code{t_norm}: predictors used to
#'           build the synthetic control.
#'     \item The column named in \code{control_var}, used to define
#'           the treated period.
#'   }
#' @param outcome Character; which outcome series to use as the response,
#'   one of \code{"I"} or \code{"C"}.
#' @param control_var Character scalar; name of a column in \code{DT} whose
#'   high values define the treated period (e.g., intensity of some
#'   intervention or shock proxy).
#' @param seed Integer; random seed for reproducibility of the BSTS fit.
#' @param dir_csv Character scalar or \code{NULL}; directory where the
#'   effect summary CSV is written. If \code{NULL} (default), nothing is
#'   written to disk.
#'
#' @details
#' This function requires the \pkg{CausalImpact} package (listed under
#' \code{Suggests}); an informative error is raised at call time if it is
#' not installed.
#'
#' The function:
#' \enumerate{
#'   \item Selects the outcome series \code{y <- DT[[outcome]]}.
#'   \item Builds the predictor matrix from
#'         \code{EconCycle}, \code{PopDensity}, \code{Epidemics},
#'         \code{Climate}, \code{War}, and \code{t_norm}.
#'   \item Uses \code{control_var} to define a treated period as
#'         observations where \code{control_var} is in the top third
#'         (\code{>= 2/3} quantile). If fewer than 5 treated observations
#'         are found, the function returns \code{NULL}.
#'   \item Sets the intervention start time \code{t0} as one period before
#'         the first treated index (with a minimum of 10 observations in
#'         the pre-period). The pre- and post-intervention windows are:
#'         \code{pre.period = c(1, t0)} and
#'         \code{post.period = c(t0 + 1, length(y))}.
#'   \item Calls \code{CausalImpact::CausalImpact()} on the combined
#'         \code{cbind(y, preds)} matrix, with \code{model.args = list(nseasons = 1)}.
#' }
#'
#' From the resulting \code{impact} object, the function extracts the
#' average absolute and relative effects from \code{impact$summary} and
#' stores them in a small summary table with two rows:
#' \code{"abs_effect_mean"} and \code{"rel_effect_mean"}.
#'
#' When \code{dir_csv} is supplied, a CSV file named
#' \code{"causalimpact_<control_var>_on_<outcome>.csv"} is written to that
#' directory. If \code{CausalImpact()} fails, the function returns
#' \code{NULL}.
#'
#' @return On success, a list with components:
#' \itemize{
#'   \item \code{impact}: the full \code{CausalImpact} object.
#'   \item \code{summary}: a \code{data.frame} with the mean absolute and
#'         relative effects.
#' }
#' If the treated period is too short or the model fit fails, the function
#' returns \code{NULL}.
#'
#' @examples
#' \donttest{
#' # This example runs only when 'CausalImpact' is installed.
#' if (requireNamespace("CausalImpact", quietly = TRUE)) {
#'   DT <- data.frame(
#'     I = rpois(30, lambda = 10),
#'     C = rpois(30, lambda = 8),
#'     EconCycle = rnorm(30),
#'     PopDensity = rnorm(30),
#'     Epidemics = rnorm(30),
#'     Climate = rnorm(30),
#'     War = rnorm(30),
#'     t_norm = seq(-1, 1, length.out = 30)
#'   )
#'
#'   res_I <- run_synth_bsts(DT, outcome = "I", control_var = "War", seed = 123)
#'   if (!is.null(res_I)) {
#'     print(res_I$summary)
#'   }
#' }
#' }
#'
#' @export

run_synth_bsts <- function(DT, outcome = c("I", "C"), control_var, seed = 123,
                           dir_csv = NULL) {
  if (!requireNamespace("CausalImpact", quietly = TRUE)) {
    stop("Package 'CausalImpact' is required for run_synth_bsts(). Please install it.",
         call. = FALSE)
  }
  outcome <- match.arg(outcome)
  set.seed(seed)
  y <- DT[[outcome]]
  preds <- DT %>% 
    dplyr::select(EconCycle, PopDensity, Epidemics, Climate, War, t_norm) %>% 
    as.matrix()
  ctrl_vals <- DT[[control_var]]
  qh <- quantile(ctrl_vals, 2/3, na.rm=TRUE)
  treated_idx <- which(ctrl_vals >= qh)
  if (length(treated_idx) < 5) return(NULL)
  t0 <- min(treated_idx) - 1L
  if (t0 < 10) t0 <- 10
  pre.period <- c(1, t0)
  post.period <- c(t0+1, length(y))
  dat <- cbind(y, preds)
  out <- try({
    impact <- CausalImpact::CausalImpact(dat, pre.period, post.period, model.args = list(nseasons=1))
    summary <- data.frame(stat=c("abs_effect_mean","rel_effect_mean"),
                          value=c(impact$summary$AbsEffect["Average"], impact$summary$RelEffect["Average"]))
    if (!is.null(dir_csv)) {
      if (!dir.exists(dir_csv)) dir.create(dir_csv, recursive = TRUE)
      readr::write_csv(summary, file.path(dir_csv, sprintf("causalimpact_%s_on_%s.csv", control_var, outcome)))
    }
    list(impact=impact, summary=summary)
  }, silent=TRUE)
  if (inherits(out, "try-error")) return(NULL) else out
}