R/placebo.R

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

#' Temporal Placebo Test via Time-Index Permutations
#'
#' Implements a temporal placebo test for the bivariate hurdle model by
#' randomly permuting the time ordering of \code{DT}, re-estimating the
#' model on each permuted dataset, and comparing the PSIS-LOO ELPD of the
#' original fit against the permuted fits.
#'
#' @param DT A \code{data.table} (or \code{data.frame}) containing the data
#'   used by \code{fit_one()}.
#' @param spec Character scalar; model specification (e.g. \code{"A"},
#'   \code{"B"}, \code{"C"}, \code{"D"}) passed to \code{fit_one()}.
#' @param k Integer; lag order passed to \code{fit_one()}.
#' @param controls Character vector of control variable names passed to
#'   \code{fit_one()}.
#' @param n_perm Integer; number of temporal permutations (placebo
#'   datasets) to run.
#' @param seed Integer; base random seed used for reproducibility of the
#'   original fit and the permutations.
#' @param dir_csv Character scalar; directory path to save the summary CSV.
#'   If \code{NULL} (default), the CSV is not saved to disk.
#'
#' @details
#' The function:
#' \itemize{
#'   \item Fits the model on the original \code{DT} via \code{fit_one()},
#'     extracts \code{"log_lik_joint"} and computes PSIS-LOO (with
#'     \code{moment_match = TRUE}).
#'   \item For each of \code{n_perm} iterations, permutes the row order of
#'     \code{DT}, refits the model on the permuted data, recomputes
#'     PSIS-LOO, and stores the permuted ELPD.
#'   \item Reports, for each permutation, the original ELPD, the permuted
#'     ELPD, and their difference (\code{elpd_orig - elpd_perm}).
#' }
#'
#' This procedure evaluates whether the temporal structure captured by the
#' model is informative: if the model is exploiting genuine time dependence,
#' the original ELPD should typically be higher than that of the permuted
#' (time-scrambled) datasets.
#'
#' The function assumes that \code{fit_one()} is available in the search path.
#'
#' @return A \code{data.frame} with one row per permutation and columns:
#' \itemize{
#'   \item \code{perm}: permutation index (1, \dots, \code{n_perm}).
#'   \item \code{elpd_orig}: ELPD of the original (non-permuted) fit.
#'   \item \code{elpd_perm}: ELPD of the model fit on the permuted data.
#'   \item \code{diff}: difference \code{elpd_orig - elpd_perm}.
#' }
#'
#' @examples
#' \donttest{
#' # This example fits Stan models and therefore runs only in an
#' # interactive session with 'cmdstanr' and a working CmdStan.
#' if (interactive() && requireNamespace("cmdstanr", quietly = TRUE)) {
#'   N <- 50
#'   DT <- data.frame(
#'     time = seq_len(N),
#'     y = rpois(N, lambda = 4),
#'     X1 = rnorm(N),
#'     X2 = rnorm(N)
#'   )
#'
#'   out_placebo <- placebo_temporal(DT, spec = "C", k = 1,
#'                                   controls = c("X1", "X2"),
#'                                   n_perm = 2, seed = 999,
#'                                   dir_csv = tempdir())
#'   head(out_placebo)
#' }
#' }
#'
#' @export
placebo_temporal <- function(DT, spec="C", k=2, controls=character(0),
                             n_perm=10, seed=999, dir_csv=NULL) {
  
  # Ensure DT is a data.table to support DT[idx] syntax
  if (!requireNamespace("data.table", quietly = TRUE)) {
    stop("Package 'data.table' is required.")
  }
  DT <- data.table::as.data.table(DT)
  
  set.seed(seed)
  
  # Fit Original
  fit_orig <- fit_one(DT, k=k, spec=spec, controls=controls,
                      iter_warmup=900, iter_sampling=1200, chains=2, seed=seed)
  
  loglik_orig <- fit_orig$fit$draws("log_lik_joint", format="draws_matrix")
  
  # Handle parallel workers safely
  n_cores <- max(1L, if (requireNamespace("future", quietly = TRUE)) future::nbrOfWorkers() else 1L)
  loo_orig <- loo::loo(loglik_orig, moment_match = TRUE, cores = n_cores)
  
  res <- data.frame(perm=integer(0), elpd_orig=numeric(0), elpd_perm=numeric(0), diff=numeric(0))
  
  # Initialize progress bar
  if (requireNamespace("progressr", quietly = TRUE)) {
    progressr::with_progress({
      p <- progressr::progressor(steps = n_perm)
      
      for (i in 1:n_perm) {
        idx <- sample(nrow(DT))
        DTp <- data.table::copy(DT)[idx] # Permute rows
        
        fit_p <- fit_one(DTp, k=k, spec=spec, controls=controls,
                         iter_warmup=700, iter_sampling=1000, chains=2, seed=seed+i)
        
        llp <- fit_p$fit$draws("log_lik_joint", format="draws_matrix")
        loo_p <- loo::loo(llp, moment_match=TRUE, cores = n_cores)
        
        diff_val <- loo_orig$estimates["elpd_loo","Estimate"] - loo_p$estimates["elpd_loo","Estimate"]
        
        res <- rbind(res, data.frame(perm=i,
                                     elpd_orig=loo_orig$estimates["elpd_loo","Estimate"],
                                     elpd_perm=loo_p$estimates["elpd_loo","Estimate"],
                                     diff=diff_val))
        
        p(message = sprintf("Placebo temporal %d/%d", i, n_perm))
      }
    })
  } else {
    # Fallback without progressr
    for (i in 1:n_perm) {
      idx <- sample(nrow(DT))
      DTp <- data.table::copy(DT)[idx]
      
      fit_p <- fit_one(DTp, k=k, spec=spec, controls=controls,
                       iter_warmup=700, iter_sampling=1000, chains=2, seed=seed+i)
      
      llp <- fit_p$fit$draws("log_lik_joint", format="draws_matrix")
      loo_p <- loo::loo(llp, moment_match=TRUE, cores = n_cores)
      
      diff_val <- loo_orig$estimates["elpd_loo","Estimate"] - loo_p$estimates["elpd_loo","Estimate"]
      
      res <- rbind(res, data.frame(perm=i,
                                   elpd_orig=loo_orig$estimates["elpd_loo","Estimate"],
                                   elpd_perm=loo_p$estimates["elpd_loo","Estimate"],
                                   diff=diff_val))
      message(sprintf("Placebo temporal %d/%d", i, n_perm))
    }
  }
  
  # Save only if dir_csv is provided
  if (!is.null(dir_csv)) {
    if (!dir.exists(dir_csv)) {
      dir.create(dir_csv, recursive = TRUE)
    }
    if (requireNamespace("readr", quietly = TRUE)) {
      readr::write_csv(res, file.path(dir_csv, "placebo_temporal.csv"))
    } else {
      utils::write.csv(res, file.path(dir_csv, "placebo_temporal.csv"), row.names = FALSE)
    }
  }
  
  res
}