R/dbn.R

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

#' Fit a Two-Slice Dynamic Bayesian Network (DBN) for I, C, and Regime
#'
#' Constructs and estimates a simple two-slice Dynamic Bayesian Network
#' (DBN) over discretized versions of \code{I}, \code{C}, and \code{Regime}
#' using \pkg{bnlearn}. The network includes current and lag-1 nodes for
#' each variable, with structural constraints enforcing the DBN topology.
#'
#' @param DT A \code{data.frame} or \code{data.table} containing at least:
#'   \itemize{
#'     \item \code{I_cat}, \code{C_cat}: discretized (e.g., tercile) versions
#'       of \code{I} and \code{C}.
#'     \item \code{Regime}: categorical regime indicator.
#'   }
#'   The function internally renames these to \code{Ic}, \code{Cc}, and
#'   \code{R}, constructs their lag-1 counterparts, and drops rows with
#'   missing lags.
#' @param dir_csv Character scalar or \code{NULL}; directory where the
#'   preprocessed data (\code{"dbn_data.csv"}) is written. If \code{NULL}
#'   (default), no CSV is written.
#' @param dir_out Character scalar or \code{NULL}; directory where the
#'   fitted objects (\code{"dbn_fit.rds"}) are saved. If \code{NULL}
#'   (default), nothing is saved to disk.
#' @param dir_figs Character scalar or \code{NULL}; directory where the
#'   DAG plot (\code{"dbn_graph.png"}) is written, if \pkg{Rgraphviz} is
#'   available. If \code{NULL} (default), no figure is produced.
#'
#' @details
#' This function requires the \pkg{bnlearn} package (listed under
#' \code{Suggests}); an informative error is raised at call time if it is
#' not installed.
#' The DBN is defined on the nodes
#' \code{Ic}, \code{Cc}, \code{R}, \code{Ic_l1}, \code{Cc_l1}, \code{R_l1}.
#' A blacklist is used to forbid arrows from current to lagged nodes, while
#' a whitelist ensures arrows from lagged to current nodes:
#' \itemize{
#'   \item Blacklist: \code{Ic → Ic_l1}, \code{Cc → Cc_l1}, \code{R → R_l1}.
#'   \item Whitelist: \code{Ic_l1 → Ic}, \code{Cc_l1 → Cc}, \code{R_l1 → R}.
#' }
#'
#' The structure is learned via hill-climbing (\code{bnlearn::hc()}) with
#' BDe score (\code{score = "bde"}) and imaginary sample size \code{iss = 10}.
#' Parameters are then estimated via \code{bnlearn::bn.fit()} using Bayesian
#' estimation with the same \code{iss}.
#'
#' If \pkg{Rgraphviz} is available and \code{dir_figs} is supplied, a graph
#' of the learned DAG is produced and saved as \code{"dbn_graph.png"} in
#' that directory. When \code{dir_csv} is supplied, the preprocessed data
#' used to fit the DBN are written to \code{"dbn_data.csv"}; when
#' \code{dir_out} is supplied, the fitted objects are saved as
#' \code{"dbn_fit.rds"}.
#'
#' @return A list with components:
#' \itemize{
#'   \item \code{dag}: the learned Bayesian network structure
#'         (\code{bnlearn} \code{"bn"} object).
#'   \item \code{fit}: the fitted DBN (\code{"bn.fit"} object).
#'   \item \code{data}: the processed data frame (\code{Ic}, \code{Cc},
#'         \code{R}, and their lag-1 versions) used to learn/fit the DBN.
#' }
#'
#' @examples
#' \donttest{
#' # This example runs only when 'bnlearn' is installed.
#' if (requireNamespace("bnlearn", quietly = TRUE)) {
#'   DT <- data.frame(
#'     I_cat  = factor(sample(c("Low", "Medium", "High"), 100, replace = TRUE)),
#'     C_cat  = factor(sample(c("Low", "Medium", "High"), 100, replace = TRUE)),
#'     Regime = factor(sample(c("Growth", "Crisis"), 100, replace = TRUE))
#'   )
#'
#'   dbn_res <- run_dbn(DT)
#'   print(dbn_res$dag)
#' }
#' }
#'
#' @export

run_dbn <- function(DT, dir_csv = NULL, dir_out = NULL, dir_figs = NULL) {
  if (!requireNamespace("bnlearn", quietly = TRUE)) {
    stop("Package 'bnlearn' is required for run_dbn(). Please install it.",
         call. = FALSE)
  }
  df <- DT %>%
    mutate(Ic = I_cat, Cc = C_cat, R = Regime) %>%
    mutate(Ic_l1 = dplyr::lag(Ic,1), Cc_l1 = dplyr::lag(Cc,1), R_l1 = dplyr::lag(R,1)) %>%
    filter(!is.na(Ic_l1), !is.na(Cc_l1), !is.na(R_l1)) %>%
    dplyr::select(Ic, Cc, R, Ic_l1, Cc_l1, R_l1)

  bl <- data.frame(from = c("Ic","Cc","R"), to = c("Ic_l1","Cc_l1","R_l1"))
  wl <- data.frame(from = c("Ic_l1","Cc_l1","R_l1"), to = c("Ic","Cc","R"))

  dag <- bnlearn::hc(df, whitelist = wl, blacklist = bl, score = "bde", iss = 10)
  fit <- bnlearn::bn.fit(dag, data=df, method = "bayes", iss = 10)

  if (!is.null(dir_figs) && requireNamespace("Rgraphviz", quietly = TRUE)) {
    if (!dir.exists(dir_figs)) dir.create(dir_figs, recursive = TRUE)
    bnlearn::graphviz.plot(dag, shape = "ellipse")
    grDevices::dev.copy(png, filename = file.path(dir_figs, "dbn_graph.png"),
                        width = 900, height = 700)
    grDevices::dev.off()
  }

  if (!is.null(dir_csv)) {
    if (!dir.exists(dir_csv)) dir.create(dir_csv, recursive = TRUE)
    readr::write_csv(df, file.path(dir_csv, "dbn_data.csv"))
  }
  if (!is.null(dir_out)) {
    if (!dir.exists(dir_out)) dir.create(dir_out, recursive = TRUE)
    saveRDS(list(dag = dag, fit = fit), file.path(dir_out, "dbn_fit.rds"))
  }
  list(dag = dag, fit = fit, data = df)
}