R/utils.R

Defines functions summarise_varx_posthoc summarise_tvarstar_posthoc summarise_placebo_top3_posthoc summarise_te_top3_by_type_posthoc summarise_te_top3_posthoc summarise_te_top3_posthoc .read_te_all summarise_hurdle_top3_posthoc .write_sheet add_qsig .first_pvalue .scalar1_chr .scalar1

Documented in add_qsig .first_pvalue .read_te_all .scalar1 .scalar1_chr summarise_hurdle_top3_posthoc summarise_placebo_top3_posthoc summarise_te_top3_by_type_posthoc summarise_te_top3_posthoc summarise_tvarstar_posthoc summarise_varx_posthoc .write_sheet

#' Coerce to numeric scalar safely
#'
#' Helper to coerce an object to numeric and return a single scalar
#' if possible. If coercion fails, the result is non-finite, or the
#' length is not exactly 1, a default value is returned.
#'
#' @param x Object to be coerced to numeric.
#' @param default Numeric scalar returned when a valid scalar cannot
#'   be extracted (default is \code{NA_real_}).
#'
#' @return A numeric scalar or \code{default} if extraction fails.
#'
#' @keywords internal

.scalar1 <- function(x, default = NA_real_) {
  xnum <- suppressWarnings(as.numeric(x))
  if (length(xnum) == 1 && is.finite(xnum)) xnum else default
}

#' Coerce to character scalar safely
#'
#' Helper to coerce an object to character and return a single scalar
#' if possible. If the input has length not equal to 1 or is \code{NA},
#' a default value is returned.
#'
#' @param x Object to be coerced to character.
#' @param default Character scalar returned when a valid scalar cannot
#'   be extracted (default is \code{NA_character_}).
#'
#' @return A character scalar or \code{default} if extraction fails.
#'
#' @keywords internal

.scalar1_chr <- function(x, default = NA_character_) {
  if (length(x) == 1 && !is.na(x)) as.character(x) else default
}

#' Extract a p-value from nested test objects
#'
#' Attempts to extract a single p-value from a variety of test result
#' objects, including nested lists produced by functions in the
#' \pkg{vars} package and related diagnostics.
#'
#' @param x An object potentially containing a p-value, such as:
#'   \itemize{
#'     \item A list with element \code{p.value}.
#'     \item A list with nested elements like \code{LMh}, \code{LMFh},
#'           \code{pt.mul}, \code{jb.mul}, \code{arch.mul}, \code{arch.uni}.
#'     \item A numeric value that can be interpreted as a p-value.
#'   }
#'
#' @details
#' The function recursively explores nested list components and attempts
#' to find a scalar p-value. Special handling is included for structures
#' like \code{jb.mul$JB}. If nothing suitable is found, \code{NA_real_}
#' is returned.
#'
#' @return A numeric scalar with the first p-value found, or
#'   \code{NA_real_} if no p-value can be extracted.
#'
#' @keywords internal

.first_pvalue <- function(x) {
  if (inherits(x, "try-error") || is.null(x)) return(NA_real_)
  if (is.list(x)) {
    if (!is.null(x$p.value) && length(x$p.value) == 1) return(.scalar1(x$p.value))
    if (!is.null(x$LMh))    { pv <- .first_pvalue(x$LMh);    if (!is.na(pv)) return(pv) }
    if (!is.null(x$LMFh))   { pv <- .first_pvalue(x$LMFh);   if (!is.na(pv)) return(pv) }
    if (!is.null(x$pt.mul)) { pv <- .first_pvalue(x$pt.mul); if (!is.na(pv)) return(pv) }
    if (!is.null(x$jb.mul)) {
      if (!is.null(x$jb.mul$JB)) {
        pv <- .first_pvalue(x$jb.mul$JB); if (!is.na(pv)) return(pv)
      }
      for (el in x$jb.mul) { pv <- .first_pvalue(el); if (!is.na(pv)) return(pv) }
    }
    if (!is.null(x$arch.mul)) { pv <- .first_pvalue(x$arch.mul); if (!is.na(pv)) return(pv) }
    if (!is.null(x$arch.uni)) for (el in x$arch.uni) {
      pv <- .first_pvalue(el); if (!is.na(pv)) return(pv) }
    for (el in x) { pv <- try(.first_pvalue(el), silent = TRUE); if (!inherits(pv,"try-error") && !is.na(pv)) return(pv) }
    return(NA_real_)
  }
  .scalar1(x)
}

#' Add BH-adjusted q-values and significance stars
#'
#' Adds Benjamini-Hochberg adjusted q-values and a simple significance
#' code column based on p-values contained in a data frame.
#'
#' @param df A data frame containing at least a numeric column
#'   \code{p_value}. If \code{df} is \code{NULL} or has zero rows,
#'   it is returned unchanged.
#'
#' @details
#' The function:
#' \itemize{
#'   \item Computes \code{q_value} using \code{p.adjust(method = "BH")}.
#'   \item Creates a \code{sig} column with significance codes:
#'     \itemize{
#'       \item \code{"***"} for \code{q_value <= 0.001}
#'       \item \code{"**"}  for \code{0.001 < q_value <= 0.01}
#'       \item \code{"*"}   for \code{0.01 < q_value <= 0.05}
#'       \item \code{""}    otherwise
#'     }
#' }
#'
#' @return The input data frame with added columns \code{q_value} and
#'   \code{sig}. If \code{df} is \code{NULL} or empty, it is returned as is.
#'
#' @export

add_qsig <- function(df) {
  if (is.null(df) || !nrow(df)) return(df)
  df$q_value <- p.adjust(as.numeric(df$p_value), method = "BH")
  df$sig <- cut(df$q_value,
                breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
                labels = c("***","**","*",""))
  df
}

#' Safely write a data frame to an Excel worksheet
#'
#' Helper to write a data frame into an \pkg{openxlsx} workbook as a
#' table, replacing any existing sheet with the same name, applying
#' basic formatting, and handling empty data frames gracefully.
#'
#' @param wb An \code{openxlsx} workbook object.
#' @param sheet_name Character scalar; name of the worksheet to create
#'   or replace.
#' @param df A data frame to write. If \code{df} is \code{NULL} or has
#'   zero rows or columns, a placeholder data frame with a single column
#'   \code{info = "Sin datos"} is written instead.
#'
#' @details
#' The function:
#' \itemize{
#'   \item Removes the sheet \code{sheet_name} if it already exists.
#'   \item Adds a new worksheet with that name.
#'   \item Writes \code{df} as a data table and freezes the first row.
#'   \item Sets column widths to \code{"auto"}.
#'   \item Applies a bold style to the header row.
#' }
#'
#' @return Invisibly returns \code{NULL}. The workbook \code{wb} is
#'   modified in place.
#'
#' @keywords internal

.write_sheet <- function(wb, sheet_name, df) {
  if (sheet_name %in% openxlsx::sheets(wb)) openxlsx::removeWorksheet(wb, sheet_name)
  openxlsx::addWorksheet(wb, sheet_name)
  if (is.null(df) || !nrow(df) || !ncol(df)) {
    df <- data.frame(info = "Sin datos", stringsAsFactors = FALSE)
  }
  openxlsx::writeDataTable(wb, sheet_name, df)
  openxlsx::freezePane(wb, sheet_name, firstActiveRow = 2)
  openxlsx::setColWidths(wb, sheet_name, cols = 1:ncol(df), widths = "auto")
  hdr <- openxlsx::createStyle(textDecoration = "bold")
  openxlsx::addStyle(wb, sheet_name, hdr, rows = 1, cols = 1:ncol(df), gridExpand = TRUE)
}

#' Summarise top-3 Hurdle-NB models across control combos
#'
#' Extracts and summarises the top three Hurdle-NB specifications (by
#' estimated ELPD) from BMA selection tables, either taken from an
#' in-memory list of results or read from CSV files on disk.
#'
#' @param bma_per_combo Optional named list of BMA results by control
#'   combination, where each element contains a component \code{$table}
#'   with columns such as \code{elpd}, \code{elpd_se}, \code{weight},
#'   \code{k}, \code{hs_tau0}, \code{hs_slab_scale}, \code{hs_slab_df},
#'   etc.
#' @param dir_csv Character scalar; directory where BMA weight CSV files
#'   \code{"bma_weights_specC_ctrl*.csv"} are stored if
#'   \code{bma_per_combo} is \code{NULL} or empty.
#'
#' @details
#' If \code{bma_per_combo} is provided and non-empty, the function uses
#' its \code{$table} components. Otherwise, it scans \code{dir_csv} for
#' BMA weight files matching the pattern
#' \code{"bma_weights_specC_ctrl*.csv"} and reads them.
#'
#' All valid rows are combined, ordered by decreasing \code{elpd}, and
#' the top three models are retained. For each, a human-readable
#' configuration string summarising \code{k}, the horseshoe hyperparameters
#' and the control combo is constructed.
#'
#' @return A data frame with up to three rows and columns:
#'   \itemize{
#'     \item \code{model}: constant string \code{"Hurdle-NB"}.
#'     \item \code{config}: textual description of the specification.
#'     \item \code{elpd}, \code{elpd_se}, \code{weight}:
#'           selection metrics from the BMA table.
#'     \item \code{k}, \code{hs_tau0}, \code{hs_slab_scale},
#'           \code{hs_slab_df}, \code{combo}: numeric tuning parameters
#'           and control-combo tag.
#'   }
#'   If no valid tables are found, a single-row data frame with \code{NA}
#'   entries is returned.
#'
#' @export

summarise_hurdle_top3_posthoc <- function(bma_per_combo, dir_csv) {
  tabs <- list()
  if (!is.null(bma_per_combo) && length(bma_per_combo)) {
    for (tag in names(bma_per_combo)) {
      tb <- try(bma_per_combo[[tag]]$table, silent = TRUE)
      if (!inherits(tb, "try-error") && !is.null(tb) && nrow(tb)) {
        tb$combo <- tag; tabs[[tag]] <- tb
      }
    }
  } else {
    files <- list.files(dir_csv, pattern = "^bma_weights_specC_ctrl.*\\.csv$", full.names = TRUE)
    for (fp in files) {
      tb <- try(suppressWarnings(readr::read_csv(fp, show_col_types = FALSE)), silent = TRUE)
      if (!inherits(tb, "try-error") && nrow(tb)) {
        tag <- sub("^.*ctrl(.*)\\.csv$", "\\1", basename(fp)); tb$combo <- tag; tabs[[fp]] <- tb
      }
    }
  }
  if (!length(tabs)) {
    return(data.frame(model="Hurdle-NB", config=NA_character_, elpd=NA_real_, elpd_se=NA_real_, weight=NA_real_,
                      k=NA_real_, hs_tau0=NA_real_, hs_slab_scale=NA_real_, hs_slab_df=NA_real_, combo=NA_character_,
                      stringsAsFactors = FALSE))
  }
  all <- dplyr::bind_rows(tabs)
  all <- dplyr::arrange(all, dplyr::desc(elpd))
  top3 <- utils::head(all, 3)
  conf <- function(r) {
    sprintf("spec=C, k=%s, hs=(tau0=%.3f, slab=%.2f, df=%.0f), controls=%s",
            .scalar1_chr(r$k,"?"), .scalar1(r$hs_tau0), .scalar1(r$hs_slab_scale),
            .scalar1(r$hs_slab_df), .scalar1_chr(r$combo,"None"))
  }
  data.frame(
    model         = rep("Hurdle-NB", nrow(top3)),
    config        = vapply(seq_len(nrow(top3)), function(i) conf(top3[i,]), character(1)),
    elpd          = as.numeric(top3$elpd),
    elpd_se       = as.numeric(top3$elpd_se),
    weight        = as.numeric(top3$weight),
    k             = as.integer(top3$k),
    hs_tau0       = as.numeric(top3$hs_tau0),
    hs_slab_scale = as.numeric(top3$hs_slab_scale),
    hs_slab_df    = as.numeric(top3$hs_slab_df),
    combo         = as.character(top3$combo),
    stringsAsFactors = FALSE
  )
}

#' Read transfer entropy results from CSV files
#'
#' Helper to load transfer entropy results from a combined CSV file or
#' from separate per-type files (counts, rates, binary) if the combined
#' file is not available.
#'
#' @param dir_csv Character scalar; directory where the transfer entropy
#'   CSV files are stored.
#'
#* @details
#' The function first looks for \code{"transfer_entropy.csv"}. If present,
#' it is read and returned. Otherwise it attempts to read:
#' \itemize{
#'   \item \code{"transfer_entropy_counts.csv"}
#'   \item \code{"transfer_entropy_rates.csv"}
#'   \item \code{"transfer_entropy_binary.csv"}
#' }
#' and combines them into a single data frame with a \code{type} column.
#'
#' @return A data frame with transfer entropy results (potentially
#'   combining several types) or \code{NULL} if no files are found.
#'
#' @keywords internal

.read_te_all <- function(dir_csv) {
  fp_all <- file.path(dir_csv, "transfer_entropy.csv")
  if (file.exists(fp_all)) {
    tb <- suppressWarnings(readr::read_csv(fp_all, show_col_types = FALSE))
    return(tb)
  }
  gl <- list(
    counts = "transfer_entropy_counts.csv",
    rates  = "transfer_entropy_rates.csv",
    binary = "transfer_entropy_binary.csv"
  )
  tabs <- list()
  for (nm in names(gl)) {
    fp <- file.path(dir_csv, gl[[nm]])
    if (file.exists(fp)) {
      tb <- suppressWarnings(readr::read_csv(fp, show_col_types = FALSE))
      if (nrow(tb)) { tb$type <- nm; tabs[[nm]] <- tb }
    }
  }
  if (!length(tabs)) return(NULL)
  dplyr::bind_rows(tabs)
}

summarise_te_top3_posthoc <- function(te_tab, dir_csv) {
  if (is.null(te_tab) || !nrow(te_tab)) te_tab <- .read_te_all(dir_csv)
  if (is.null(te_tab) || !nrow(te_tab)) {
    return(data.frame(model="TransferEntropy", config=NA_character_, stat=NA_real_, p_value=NA_real_,
                      stringsAsFactors = FALSE))
  }
  long <- dplyr::bind_rows(
    dplyr::transmute(te_tab, dir="I->C", lag=lag, stat=TE_ItoC, p=p_ItoC, type = .data[["type"]]),
    dplyr::transmute(te_tab, dir="C->I", lag=lag, stat=TE_CtoI, p=p_CtoI, type = .data[["type"]])
  )
  long$p <- as.numeric(long$p)
  long <- long[order(long$p, long$lag), ]
  if (!nrow(long)) {
    return(data.frame(model="TransferEntropy", config=NA_character_, stat=NA_real_, p_value=NA_real_,
                      stringsAsFactors = FALSE))
  }
  top3 <- utils::head(long, 3)
  data.frame(
    model   = rep("TransferEntropy", nrow(top3)),
    config  = sprintf("dir=%s, lag=%d%s",
                      top3$dir, top3$lag,
                      ifelse(is.na(top3$type), "", paste0(", type=", top3$type))),
    stat    = as.numeric(top3$stat),
    p_value = as.numeric(top3$p),
    stringsAsFactors = FALSE
  )
}

#' Summarise top-3 transfer entropy results (global)
#'
#' Produces a compact summary of the three most statistically significant
#' transfer entropy estimates across directions and lags, optionally
#' combining information from counts, rates, and binary specifications.
#'
#' @param te_tab Optional data frame with transfer entropy results,
#'   containing at least columns \code{lag}, \code{TE_ItoC},
#'   \code{TE_CtoI}, \code{p_ItoC}, \code{p_CtoI}, and optionally
#'   \code{type}. If \code{NULL} or empty, the function attempts to read
#'   the data from CSV files via the internal helper \code{.read_te_all()}.
#' @param dir_csv Character scalar; directory where the transfer entropy
#'   CSV files are stored (used when \code{te_tab} is missing).
#'
#' @details
#' The function reshapes \code{te_tab} into a long format with directions
#' \code{"I->C"} and \code{"C->I"}, orders by p-value (ascending) and
#' lag, and keeps the three rows with the smallest p-values.
#'
#' @return A data frame with up to three rows and columns:
#'   \itemize{
#'     \item \code{model}: constant string \code{"TransferEntropy"}.
#'     \item \code{config}: textual description of direction, lag, and,
#'           if available, type (counts, rates, binary).
#'     \item \code{stat}: transfer entropy estimate.
#'     \item \code{p_value}: associated p-value.
#'   }
#'   If no results are available, a single-row data frame with \code{NA}
#'   entries is returned.
#'
#' @export

summarise_te_top3_posthoc <- function(te_tab, dir_csv) {
  if (is.null(te_tab) || !nrow(te_tab)) te_tab <- .read_te_all(dir_csv)
  if (is.null(te_tab) || !nrow(te_tab)) {
    return(data.frame(model="TransferEntropy", config=NA_character_, stat=NA_real_, p_value=NA_real_,
                      stringsAsFactors = FALSE))
  }
  long <- dplyr::bind_rows(
    dplyr::transmute(te_tab, dir="I->C", lag=lag, stat=TE_ItoC, p=p_ItoC, type = .data[["type"]]),
    dplyr::transmute(te_tab, dir="C->I", lag=lag, stat=TE_CtoI, p=p_CtoI, type = .data[["type"]])
  )
  long$p <- as.numeric(long$p)
  long <- long[order(long$p, long$lag), ]
  if (!nrow(long)) {
    return(data.frame(model="TransferEntropy", config=NA_character_, stat=NA_real_, p_value=NA_real_,
                      stringsAsFactors = FALSE))
  }
  top3 <- utils::head(long, 3)
  data.frame(
    model   = rep("TransferEntropy", nrow(top3)),
    config  = sprintf("dir=%s, lag=%d%s",
                      top3$dir, top3$lag,
                      ifelse(is.na(top3$type), "", paste0(", type=", top3$type))),
    stat    = as.numeric(top3$stat),
    p_value = as.numeric(top3$p),
    stringsAsFactors = FALSE
  )
}

#' Summarise top-3 transfer entropy results by type
#'
#' Produces a list of small tables with the three most significant
#' transfer entropy estimates for each data type (counts, rates, binary)
#' separately.
#'
#' @param te_tab Optional data frame with transfer entropy results,
#'   including a \code{type} column and at least \code{lag},
#'   \code{TE_ItoC}, \code{TE_CtoI}, \code{p_ItoC}, \code{p_CtoI}.
#'   If \code{NULL} or empty, the function attempts to read the data
#'   from CSV files via \code{.read_te_all()}.
#' @param dir_csv Character scalar; directory where the transfer entropy
#'   CSV files are stored (used when \code{te_tab} is missing).
#'
#' @details
#' For each type in \code{c("counts", "rates", "binary")}, the function
#' ranks all direction-lag combinations by p-value and retains the top
#' three. Types with no valid rows remain \code{NULL} in the output list.
#'
#' @return A named list with up to three elements:
#'   \itemize{
#'     \item \code{$counts}, \code{$rates}, \code{$binary}: each is a
#'           data frame with columns \code{model}, \code{type},
#'           \code{config} (direction and lag), \code{stat}, and
#'           \code{p_value}, or \code{NULL} if no results for that type.
#'   }
#'
#' @export

summarise_te_top3_by_type_posthoc <- function(te_tab, dir_csv) {
  if (is.null(te_tab) || !nrow(te_tab)) te_tab <- .read_te_all(dir_csv)
  out <- list(counts = NULL, rates = NULL, binary = NULL)
  if (is.null(te_tab) || !nrow(te_tab) || !"type" %in% names(te_tab)) return(out)
  types <- intersect(unique(te_tab$type), c("counts","rates","binary"))
  for (tp in types) {
    sub <- te_tab[te_tab$type == tp, , drop = FALSE]
    long <- dplyr::bind_rows(
      dplyr::transmute(sub, dir="I->C", lag=lag, stat=TE_ItoC, p=p_ItoC),
      dplyr::transmute(sub, dir="C->I", lag=lag, stat=TE_CtoI, p=p_CtoI)
    )
    long$p <- as.numeric(long$p)
    long <- long[order(long$p, long$lag), ]
    top3 <- utils::head(long, 3)
    if (nrow(top3)) {
      out[[tp]] <- data.frame(
        model   = "TransferEntropy",
        type    = tp,
        config  = sprintf("dir=%s, lag=%d", top3$dir, top3$lag),
        stat    = as.numeric(top3$stat),
        p_value = as.numeric(top3$p),
        stringsAsFactors = FALSE
      )
    }
  }
  out
}

#' Summarise top-3 temporal placebo results
#'
#' Summarises the three strongest temporal placebo results (based on
#' the difference between original and permuted ELPD) from a temporal
#' permutation test.
#'
#' @param placebo_tab Optional data frame with placebo results, typically
#'   containing columns \code{perm}, \code{elpd_orig}, \code{elpd_perm},
#'   and \code{diff}. If \code{NULL} or empty, the function attempts to
#'   read \code{"placebo_temporal.csv"} from \code{dir_csv}.
#' @param dir_csv Character scalar; directory where the placebo CSV file
#'   is stored.
#'
#' @details
#' The table is ordered by decreasing \code{diff} (ELPD gain of the
#' original fit over the permuted fit), and the top three permutations
#' are retained.
#'
#' @return A data frame with up to three rows and columns:
#'   \itemize{
#'     \item \code{model}: constant string \code{"PlaceboTemporal"}.
#'     \item \code{config}: text of the form \code{"perm=<id>"}.
#'     \item \code{elpd_orig}, \code{elpd_perm}, \code{diff}: original
#'           ELPD, permuted ELPD, and their difference.
#'   }
#'   If no data are available, a single-row data frame with \code{NA}
#'   entries is returned.
#'
#' @export

summarise_placebo_top3_posthoc <- function(placebo_tab, dir_csv) {
  if (is.null(placebo_tab) || !nrow(placebo_tab)) {
    fp <- file.path(dir_csv, "placebo_temporal.csv")
    if (file.exists(fp)) placebo_tab <- suppressWarnings(readr::read_csv(fp, show_col_types = FALSE))
  }
  if (is.null(placebo_tab) || !nrow(placebo_tab)) {
    return(data.frame(model="PlaceboTemporal", config=NA_character_, elpd_orig=NA_real_,
                      elpd_perm=NA_real_, diff=NA_real_, stringsAsFactors = FALSE))
  }
  placebo_tab <- dplyr::arrange(placebo_tab, dplyr::desc(diff))
  top3 <- utils::head(placebo_tab, 3)
  data.frame(
    model     = rep("PlaceboTemporal", nrow(top3)),
    config    = sprintf("perm=%d", as.integer(top3$perm)),
    elpd_orig = as.numeric(top3$elpd_orig),
    elpd_perm = as.numeric(top3$elpd_perm),
    diff      = as.numeric(top3$diff),
    stringsAsFactors = FALSE
  )
}

#' Summarise nonlinear time-series models (TVAR and LSTAR)
#'
#' Produces a small summary table for nonlinear time-series models
#' such as TVAR and LSTAR, focusing on model status and AIC.
#'
#' @param tsdyn_res A list of model objects, typically with elements
#'   \code{$TVAR}, \code{$LSTAR_I}, \code{$LSTAR_C}, as returned by
#'   a fitting routine based on the \pkg{tsDyn} package.
#'
#' @details
#' For each of the three models (TVAR, LSTAR for I, LSTAR for C),
#' the function extracts:
#' \itemize{
#'   \item A textual status (class names of the object).
#'   \item The AIC, if \code{stats::AIC()} can be computed.
#' }
#' If \code{tsdyn_res} is \code{NULL}, default rows with \code{NA}
#' values are returned.
#'
#' @return A data frame with one row per model and columns:
#'   \itemize{
#'     \item \code{model}: "TVAR", "LSTAR_I", "LSTAR_C".
#'     \item \code{status}: model class string or \code{NA}.
#'     \item \code{aic}: numeric AIC value or \code{NA}.
#'   }
#'
#' @export

summarise_tvarstar_posthoc <- function(tsdyn_res) {
  if (is.null(tsdyn_res)) {
    return(data.frame(model=c("TVAR","LSTAR_I","LSTAR_C"),
                      status=NA_character_, aic=NA_real_, stringsAsFactors = FALSE))
  }
  get_aic <- function(obj) tryCatch(stats::AIC(obj), error = function(e) NA_real_)
  row <- function(name, obj) data.frame(
    model  = name,
    status = .scalar1_chr(if (is.null(obj)) NA_character_ else paste(class(obj), collapse="+")),
    aic    = .scalar1(get_aic(obj)),
    stringsAsFactors = FALSE
  )
  dplyr::bind_rows(
    row("TVAR",    tsdyn_res$TVAR),
    row("LSTAR_I", tsdyn_res$LSTAR_I),
    row("LSTAR_C", tsdyn_res$LSTAR_C)
  )
}

#' Summarise VARX model fit and diagnostics
#'
#' Produces a compact summary of a VARX model, including information
#' about lag order, exogenous variables, information criteria, and
#' selected diagnostic p-values.
#'
#' @param varx_res A list returned by \code{\link{run_varx}()}, typically
#'   containing elements \code{$fit}, \code{$serial}, \code{$normal},
#'   and \code{$arch}.
#'
#' @details
#' The function extracts:
#' \itemize{
#'   \item Lag order \code{p} from \code{fit$p}, if available.
#'   \item AIC and BIC via \code{stats::AIC()} and \code{stats::BIC()}.
#'   \item P-values from serial correlation, normality, and ARCH tests
#'         using the helper \code{.first_pvalue()}.
#' }
#'
#' If \code{varx_res} or \code{varx_res$fit} is \code{NULL}, a default
#' row with \code{NA} values is returned.
#'
#' @return A data frame with one row and columns:
#'   \itemize{
#'     \item \code{model}: constant string \code{"VARX"}.
#'     \item \code{config}: textual description of the lag order and
#'           exogenous variables.
#'     \item \code{AIC}, \code{BIC}: information criteria.
#'     \item \code{p_serial}, \code{p_normal}, \code{p_arch}: p-values
#'           from diagnostic tests.
#'   }
#'
#' @export

summarise_varx_posthoc <- function(varx_res) {
  if (is.null(varx_res) || is.null(varx_res$fit)) {
    return(data.frame(model="VARX", config=NA_character_, AIC=NA_real_, BIC=NA_real_,
                      p_serial=NA_real_, p_normal=NA_real_, p_arch=NA_real_, stringsAsFactors = FALSE))
  }
  fit <- varx_res$fit
  p_order <- tryCatch(fit$p, error = function(e) NA_integer_)
  exog_vec <- c("EconCycle","PopDensity","Epidemics","Climate","War","t_norm")
  cfg <- sprintf("p=%s, exogen=%s",
                 ifelse(is.na(p_order), "?", as.character(p_order)),
                 paste(exog_vec, collapse = "+"))
  data.frame(
    model    = "VARX",
    config   = cfg,
    AIC      = .scalar1(tryCatch(stats::AIC(fit), error = function(e) NA_real_)),
    BIC      = .scalar1(tryCatch(stats::BIC(fit), error = function(e) NA_real_)),
    p_serial = .first_pvalue(varx_res$serial),
    p_normal = .first_pvalue(varx_res$normal),
    p_arch   = .first_pvalue(varx_res$arch),
    stringsAsFactors = FALSE
  )
}

Try the bivarhr package in your browser

Any scripts or data that you put into this service are public.

bivarhr documentation built on July 7, 2026, 1:06 a.m.