R/balqual.R

In vecmatch: Generalized Propensity Score Estimation and Matching for Multiple Groups

#' @title Evaluate Matching Quality
#'
#' @description The `balqual()` function evaluates the balance quality of a
#'   dataset after matching, comparing it to the original unbalanced dataset. It
#'   computes various summary statistics and provides an easy interpretation
#'   using user-specified cutoff values.
#' @param matched_data An object of class `matched`, generated by the
#'   [match_gps()] function. This object is essential for the `balqual()`
#'   function as it contains the final data.frame and attributes required to
#'   compute the quality coefficients.
#' @param formula A valid R formula used to compute generalized propensity
#'   scores during the first step of the vector matching algorithm in
#'   [estimate_gps()]. This formula must match the one used in `estimate_gps()`.
#' @param type A character vector specifying the quality metrics to calculate,
#'   provided in a vector created by the `c()`. The values fall into two
#'   groups. `smd`, `r` and `var_ratio` are balance metrics summarizing
#'   *pairwise* comparisons between the treatment levels. `desc_full` and
#'   `desc_reduced` instead describe the marginal distribution of each
#'   covariate within every treatment level, which makes it possible to
#'   directly compare the covariate distributions of the unmatched and the
#'   matched dataset. Note that they describe the covariates as they are given
#'   in the `formula`, whereas the balance metrics are computed on the model
#'   matrix: a factor is described once, by its levels, rather than once per
#'   dummy-coded column, and interaction terms are not described at all, since
#'   they are model terms rather than covariates. The balance metrics still
#'   cover them. Neither descriptive metric is part of the default, and both
#'   have to be requested explicitly, either on their own or together with the
#'   balance metrics, e.g. `type = c("smd", "desc_reduced")`. As `desc_reduced`
#'   is a subset of `desc_full`, the two are mutually exclusive. Possible
#'   values include:
#' * `smd` - Calculates standardized mean differences (SMD) between groups,
#'   defined as the difference in means divided by the standard deviation of the
#'   treatment group (Rubin, 2001).
#' * `r` - Computes Pearson's r coefficient using the Z statistic from the
#'   U-Mann-Whitney test.
#' * `var_ratio` - Measures the dispersion differences between groups,
#'   calculated as the ratio of the larger variance to the smaller one.
#' * `desc_full` - Computes standard descriptive statistics for every balancing
#'   variable, separately for each treatment level, before and after matching.
#'   Numeric covariates are summarized by the number of observations (`N`),
#'   mean, standard deviation, minimum, first quartile, median, third quartile,
#'   maximum, skewness and excess kurtosis, where skewness and excess kurtosis
#'   are the classical moment (type 1) estimators. Categorical covariates
#'   (factor, character and logical) are instead cross-tabulated, reporting the
#'   count and the percentage of every level within each treatment level, as
#'   means and quantiles are not meaningful for them. The two are printed as
#'   two separate tables.
#' * `desc_reduced` - As `desc_full`, but the summary of the numeric covariates
#'   is restricted to the four most commonly reported statistics: minimum,
#'   mean, median and maximum. The cross-tabulation of the categorical
#'   covariates is unaffected. Useful when the full table is too wide to read
#'   comfortably in the console.
#' @param statistic A character vector specifying the type of statistics used to
#'   summarize the quality metrics. Since quality metrics are calculated for all
#'   pairwise comparisons between treatment levels, they need to be aggregated
#'   for the entire dataset.
#'   - `max`: Returns the maximum values of the statistics defined in the `type`
#'   argument (as suggested by Lopez and Gutman, 2017).
#'   - `mean`: Returns the corresponding averages.
#'
#'   To compute both, provide both names using the `c()` function. This argument
#'   is ignored by the `desc_full` and `desc_reduced` metrics, which are not
#'   based on pairwise comparisons and are therefore never aggregated across
#'   them.
#' @param cutoffs A numeric vector with the same length as the number of
#'   balance metrics specified in the `type` argument. Defines the cutoffs for
#'   each corresponding metric, below which the dataset is considered balanced.
#'   If `NULL`, the default cutoffs are used: 0.1 for `smd` and `r`, and 2 for
#'   `var_ratio`. The descriptive metrics have no cutoff and are not counted
#'   here, so `type = c("smd", "desc_full")` still requires a single cutoff
#'   value.
#' @param round A single non-negative integer specifying the number of
#'   decimal places to round the output to.
#'
#' @return If assigned to a name, returns a list of summary statistics of class
#'   `quality` containing:
#'  * `quality_mean` - A data frame with the mean values of the statistics
#'   specified in the `type` argument for all balancing variables used in
#'   `formula`.
#'  * `quality_max` - A data frame with the maximal values of the statistics
#'   specified in the `type` argument for all balancing variables used in
#'   `formula`.
#'  * `quality_desc` - A data frame with the descriptive statistics of every
#'   balancing variable, laid out like the balance tables above: the statistics
#'   are the rows and the two matching stages are the `Before` and `After`
#'   columns. It holds the numeric and the categorical covariates in a single
#'   table, with the columns `Variable`, `Group`, `Type`, `Statistic`,
#'   `Before`, `After`, `Percent_Before` and `Percent_After`. `Type` is either
#'   `"continuous"` or `"categorical"`. For a continuous covariate, `Statistic`
#'   names the reported statistic - `N`, `Mean`, `SD`, `Min`, `Q1`, `Median`,
#'   `Q3`, `Max`, `Skewness` and `Kurtosis` for `desc_full`, or `N`, `Min`,
#'   `Mean`, `Median` and `Max` for `desc_reduced` - and the two `Percent`
#'   columns are `NA`. For a categorical covariate, `Statistic` names a level
#'   of that covariate, `Before` and `After` hold its counts, and the two
#'   `Percent` columns hold its share within the treatment level, computed
#'   separately for each matching stage. `NULL` unless one of the two
#'   descriptive metrics is included in the `type` argument. When printed, the
#'   counts and percentages of a categorical covariate are compressed into a
#'   single `N (%)` cell, with the percentage always shown to one decimal
#'   place.
#'  * `perc_matched` - A single numeric value indicating the percentage of
#'   observations in the original dataset that were matched.
#'  * `statistic` - A single string defining which statistic will be displayed
#'   in the console.
#'  * `summary_head` - A summary of the matching process. If `max` is included
#'   in the `statistic`, it contains the maximal observed values for each
#'   variable; otherwise, it includes the mean values.
#'  * `n_before` - The number of observations in the dataset before matching.
#'  * `n_after` - The number of observations in the dataset after matching.
#'  * `count_table` - A contingency table showing the distribution of the
#'   treatment variable before and after matching.
#'
#'   The `balqual()` function also prints a well-formatted table with the
#'   defined summary statistics for each variable in the `formula` to the
#'   console.
#'
#' @examples
#' # We try to balance the treatment variable in the cancer dataset based on age
#' # and sex covariates
#' data(cancer)
#'
#' # Then we can estimate the generalized propensity scores
#' gps_cancer <- estimate_gps(formula(status ~ age * sex),
#'   cancer,
#'   method = "multinom",
#'   reference = "control",
#'   verbose_output = TRUE
#' )
#'
#' # ... and drop observations based on the common support region...
#' csr_cancer <- csregion(gps_cancer)
#'
#' # ... to match the samples using `match_gps()`
#' matched_cancer <- match_gps(csr_cancer,
#'   reference = "control",
#'   caliper = 1,
#'   kmeans_cluster = 5,
#'   kmeans_args = list(n.iter = 100),
#'   verbose_output = TRUE
#' )
#'
#' # At the end we can assess the quality of matching using `balqual()`
#' balqual(
#'   matched_data = matched_cancer,
#'   formula = formula(status ~ age * sex),
#'   type = "smd",
#'   statistic = "max",
#'   round = 3,
#'   cutoffs = 0.2
#' )
#'
#' # Adding `desc_full` to `type` additionally reports the descriptive
#' # statistics of every covariate per treatment level, before and after
#' # matching, which allows a direct comparison of the covariate distributions
#' balqual(
#'   matched_data = matched_cancer,
#'   formula = formula(status ~ age * sex),
#'   type = c("smd", "desc_full"),
#'   statistic = "max",
#'   round = 3,
#'   cutoffs = 0.2
#' )
#'
#' # `desc_reduced` gives the same table restricted to the minimum, mean,
#' # median and maximum
#' balqual(
#'   matched_data = matched_cancer,
#'   formula = formula(status ~ age * sex),
#'   type = c("smd", "desc_reduced"),
#'   statistic = "max",
#'   round = 3,
#'   cutoffs = 0.2
#' )
#'
#' @seealso [match_gps()] for matching the generalized propensity scores;
#' [estimate_gps()] for the documentation of the `formula` argument.
#'
#' @references Rubin, D.B. Using Propensity Scores to Help Design Observational
#'   Studies: Application to the Tobacco Litigation. Health Services & Outcomes
#'   Research Methodology 2, 169–188 (2001).
#'   https://doi.org/10.1023/A:1020363010465
#'
#'   Michael J. Lopez, Roee Gutman "Estimation of Causal Effects with Multiple
#'   Treatments: A Review and New Ideas," Statistical Science, Statist. Sci.
#'   32(3), 432-454, (August 2017)
#' @export
balqual <- function(matched_data = NULL,
                    formula = NULL,
                    type = c("smd", "r", "var_ratio"),
                    statistic = c("mean", "max"),
                    cutoffs = NULL,
                    round = 3) {
  ############ ARGUMENT CHECKING AND PROCESSING ################################
  # check data
  .chk_cond(
    is.null(matched_data),
    "The argument `matched_data` is missing with no default!"
  )

  .chk_cond(
    !inherits(matched_data, c("data.frame", "matched")),
    "The argument `matched_data` has to be a data.frame of class
            `matched`!"
  )

  ## get rid of other classes
  old_data <- data.frame(attr(matched_data, "original_data"))
  new_data <- data.frame(matched_data)

  # check and process formula
  data_before <- .process_formula(formula, old_data)
  data_after <- .process_formula(formula, new_data)

  # define all posssible pairwise comaparisons of treatment var
  pairwise_comb <- t(utils::combn(
    unique(as.character(data_before[["treat"]])),
    2
  ))
  pairwise_comb <- data.frame(pairwise_comb)
  pairwise_comb <- pairwise_comb[stats::complete.cases(pairwise_comb), ]
  colnames(pairwise_comb) <- c("group1", "group2")

  # --- check and process `type` ----------------------------------------------
  .chk_cond(
    !.check_vecl(type, check_numeric = FALSE),
    "The `type` argument has to be an atomic character vector."
  )

  # normalise to lower case
  type <- tolower(type)

  allowed_desc <- c("desc_full", "desc_reduced")
  allowed_type <- c("smd", "r", "var_ratio", allowed_desc)

  .chk_cond(
    !all(type %in% allowed_type),
    sprintf(
      "The `type` argument must be one of %s.",
      paste(add_quotes(allowed_type), collapse = ", ")
    )
  )

  # `desc_reduced` is a subset of `desc_full`, so requesting both is ambiguous
  .chk_cond(
    all(allowed_desc %in% type),
    sprintf(
      "The %s and %s values of the `type` argument are mutually exclusive, as
              the latter is a subset of the former.",
      add_quotes("desc_full"), add_quotes("desc_reduced")
    )
  )

  # the descriptive metrics are not pairwise balance metrics, so they are
  # handled separately from the metrics that are aggregated, compared to
  # cutoffs and summarized
  desc_type <- intersect(type, allowed_desc)
  compute_desc <- length(desc_type) > 0L
  balance_type <- setdiff(type, allowed_desc)

  # --- check and process `statistic` -----------------------------------------
  .chk_cond(
    !.check_vecl(statistic, check_numeric = FALSE),
    "The `statistic` argument has to be an atomic character vector."
  )

  statistic <- tolower(statistic)

  allowed_stat <- c("mean", "max")

  .chk_cond(
    !all(statistic %in% allowed_stat),
    sprintf(
      "The `statistic` argument must be one of %s.",
      paste(add_quotes(allowed_stat), collapse = ", ")
    )
  )

  # `statistic` aggregates pairwise metrics, so it does nothing when the only
  # requested metrics are the descriptive ones
  .chk_cond(
    compute_desc && length(balance_type) == 0L &&
      "statistic" %in% names(match.call()),
    sprintf(
      "The `statistic` argument is ignored for the %s type, as descriptive
              statistics are computed per treatment level and not aggregated
              over pairwise comparisons.",
      add_quotes(desc_type)
    ),
    error = FALSE
  )

  # check cutoffs
  if (is.null(cutoffs)) {
    cutoffs <- unlist(lapply(balance_type, function(x) {
      switch(x,
        smd = 0.1,
        r = 0.1,
        var_ratio = 2
      )
    }))
  }

  .chk_cond(
    length(balance_type) > 0L &&
      !.check_vecl(cutoffs, length(balance_type), check_numeric = TRUE),
    "The argument `cutoffs` has to be an numeric vector, with length
            equal to the number of balance metrics in the `type` argument
            (`desc` is not counted)."
  )

  # `create_balqual_output()` expects one cutoff per metric, in the fixed order
  # smd, r, var_ratio. Mapping by name keeps the requested cutoffs aligned with
  # their metric and leaves the defaults in place for the unrequested ones.
  cutoffs_full <- c(smd = 0.1, r = 0.1, var_ratio = 2)

  if (length(balance_type) > 0L) {
    cutoffs_full[balance_type] <- cutoffs
  }

  # check round: single non-negative integer
  .chk_cond(
    length(round) != 1L,
    "The argument `round` must be a single non-negative integer."
  )

  .check_integer(round, "round")

  round <- as.integer(round)

  .chk_cond(
    round < 0L,
    "The argument `round` must be a non-negative integer."
  )

  ############ PERFORMING THE CALCULATIONS #####################################
  # define output list
  quality_list <- .make_list(nrow(pairwise_comb))

  # define output data.frame for smd's (important for the optimization)
  new_colnames <- colnames(data_after[["model_covs"]])

  # Create a data.frame with same number of rows as pairwise_comb
  na_columns <- data.frame(matrix(NA,
    nrow = nrow(pairwise_comb),
    ncol = length(new_colnames)
  ))
  colnames(na_columns) <- new_colnames

  # Combine column-wise
  smd_df <- cbind(pairwise_comb, na_columns)

  # loop over all pairwise combinations to calculate the quality metrics
  for (i in seq_len(nrow(pairwise_comb))) {
    # helper function to recode data to 0 (control) and 1 (treatment)
    .binarize_treat <- function(treatment, comb1, comb2) {
      treat_recoded <- rep(NA_real_, length(treatment))

      # assign 0 to controls
      treat_recoded[treatment == pairwise_comb[i, 1]] <- 0

      # and 1 to treatments
      treat_recoded[treatment == pairwise_comb[i, 2]] <- 1

      treat_recoded
    }

    # recode treatment to binary
    treat_before <- .binarize_treat(
      data_before[["treat"]],
      pairwise_comb[i, 1],
      pairwise_comb[i, 2]
    )

    treat_after <- .binarize_treat(
      data_after[["treat"]],
      pairwise_comb[i, 1],
      pairwise_comb[i, 2]
    )

    # define functions to calculate smd, r and %matched
    # smd
    standardized_bias <- function(treatment, covariate) {
      remove_nas <- which(!is.na(treatment))
      covariate <- covariate[remove_nas]
      treatment <- treatment[remove_nas]

      means <- tapply(
        covariate,
        treatment,
        mean
      )
      sd_treat <- stats::sd(covariate[treatment == 1]) # as suggested by Rubin
      (means[1] - means[2]) / sd_treat
    }

    # r
    wilcox_r <- function(treatment, covariate) {
      data_split <- split(covariate, treatment)

      p_value <- stats::wilcox.test(
        x = data_split[["0"]],
        y = data_split[["1"]],
        paired = FALSE
      )$p.value

      z_stat <- stats::qnorm(p_value / 2)

      z_stat / length(data_split[["0"]])
    }

    # variance ratio
    var_ratio <- function(treatment, covariate) {
      remove_nas <- which(!is.na(treatment))
      covariate <- covariate[remove_nas]
      treatment <- treatment[remove_nas]

      sd_vals <- tapply(
        covariate,
        treatment,
        stats::sd
      )

      max(sd_vals) / min(sd_vals)
    }

    # calculating the statistics
    smd_before <- apply(data_before[["model_covs"]], 2, function(x) {
      standardized_bias(treat_before, x)
    })

    smd_after <- apply(data_after[["model_covs"]], 2, function(x) {
      standardized_bias(treat_after, x)
    })

    r_effsize_before <- apply(data_before[["model_covs"]], 2, function(x) {
      wilcox_r(treat_before, x)
    })

    r_effsize_after <- apply(data_after[["model_covs"]], 2, function(x) {
      wilcox_r(treat_after, x)
    })

    variance_before <- apply(data_before[["model_covs"]], 2, function(x) {
      var_ratio(treat_before, x)
    })

    variance_after <- apply(data_after[["model_covs"]], 2, function(x) {
      var_ratio(treat_after, x)
    })

    variables <- rbind(
      smd_before,
      smd_after,
      r_effsize_before,
      r_effsize_after,
      variance_before,
      variance_after
    )

    variables <- abs(variables)

    # convert to characters for duplicates check
    var_uniques <- data.frame(lapply(t(variables), function(x) {
      as.character(round(x, 8))
    }))

    var_uniques <- duplicated(var_uniques)

    # subset the vars by unique values
    variables <- variables[, !var_uniques, drop = FALSE]
    colnames(variables) <- colnames(data_before[["model_covs"]])[!var_uniques]

    # assign to list
    quality_list[[i]] <- as.data.frame(variables)

    # save the smd's to smd_df for optimization selection
    smd_df[i, 3:ncol(smd_df)] <- abs(smd_after)
  }

  # output data.frames
  quality_dataframe <- data.frame(
    coef_name = rep(c("smd", "r", "var_ratio"),
      each = 2
    ),
    time = rep(c("before", "after"), 3)
  )

  # means
  quality_mean <- create_balqual_output(quality_list,
    quality_dataframe,
    operation = "+",
    round = round,
    which_coefs = balance_type,
    cutoffs = unname(cutoffs_full)
  )

  # maxes
  quality_max <- create_balqual_output(quality_list,
    quality_dataframe,
    operation = "max",
    round = round,
    which_coefs = balance_type,
    cutoffs = unname(cutoffs_full)
  )

  # descriptive statistics of the covariates, per treatment level
  quality_desc <- NULL

  if (compute_desc) {
    # the descriptives describe the covariates as they appear in the `formula`,
    # not the dummy-coded model matrix used for the balance metrics
    desc_long <- rbind(
      .desc_table(
        data_before[["reported_covs"]],
        data_before[["treat"]],
        time = "Before",
        round = round
      ),
      .desc_table(
        data_after[["reported_covs"]],
        data_after[["treat"]],
        time = "After",
        round = round
      )
    )

    # `desc_reduced` keeps the four most commonly used statistics next to the
    # number of observations. The counts of the categorical covariates have
    # nothing to reduce.
    stats_keep <- if (identical(desc_type, "desc_reduced")) {
      c("N", "Min", "Mean", "Median", "Max")
    } else {
      .desc_stat_names
    }

    # the statistics become rows and the matching stages become columns, so
    # that the table reads like the balance tables above it
    quality_desc <- .desc_to_wide(
      desc_long,
      stats_keep = stats_keep,
      var_order = names(data_before[["reported_covs"]]),
      group_order = .desc_groups(data_before[["treat"]])
    )
  }

  # % Matched
  perc_matched <- length(data_after[["treat"]]) /
    length(data_before[["treat"]]) * 100

  perc_matched <- round(perc_matched, 2)

  # calculating total maximas or means
  type_recoded <- lapply(balance_type, function(x) {
    switch(x,
      "smd" = "SMD",
      "r" = "r",
      "var_ratio" = "Var"
    )
  })

  # switch to use different dataframes based on statistics
  summary_head <- lapply(type_recoded, function(x) {
    if ("max" %in% statistic) {
      max(quality_max$After[quality_max$Coef == x])
    } else {
      max(quality_mean$After[quality_mean$Coef == x])
    }
  })

  summary_head <- unlist(summary_head)

  names(summary_head) <- unlist(type_recoded)

  # count table for treatment variable
  times <- c(
    rep("Before", length(data_before[["treat"]])),
    rep("After", length(data_after[["treat"]]))
  )

  count_data <- data.frame(
    Matching = times,
    Treatment = c(
      data_before[["treat"]],
      data_after[["treat"]]
    )
  )

  count_table <- table(
    count_data$Treatment,
    count_data$Matching
  )

  count_table <- cbind(Treatment = rownames(count_table), count_table)
  count_table <- count_table[, c("Treatment", "Before", "After")]
  rownames(count_table) <- NULL

  # defining final object (core list)
  quality_core <- list(
    quality_mean  = quality_mean,
    quality_max   = quality_max,
    quality_desc  = quality_desc,
    perc_matched  = perc_matched,
    statistic     = statistic,
    summary_head  = summary_head,
    n_before      = length(data_before[["treat"]]),
    n_after       = length(data_after[["treat"]]),
    count_table   = count_table,
    type          = type,
    cutoffs       = cutoffs,
    round         = round
  )

  # canonicalise the stored call so equivalent inputs give identical objects
  fc <- match.call()
  if ("type" %in% names(fc)) {
    fc[["type"]] <- type # always lower-case
  }
  if ("statistic" %in% names(fc)) {
    fc[["statistic"]] <- statistic # always lower-case
  }
  if ("round" %in% names(fc)) {
    fc[["round"]] <- round # after you've coerced to integer earlier
  }

  # assemble quality object with attributes and class in one go
  quality_obj <- structure(
    quality_core,
    original_data_before = old_data,
    original_data_after  = new_data,
    formula              = formula,
    function_call        = fc,
    smd_df_combo         = smd_df,
    class                = "quality"
  )

  return(quality_obj)
}

#' @export
print.quality <- function(x, ...) {
  # Helper function to print a table
  print_table <- function(df, colnames_df, ncol) {
    # Define separator based on the number of columns
    separator <- paste(rep("-", ifelse(ncol == 3, 50, 80)), collapse = "")

    # Define column formats based on the number of columns
    col_format <- if (ncol == 3) {
      c("%-25s", "%-10s", "%-10s")
    } else {
      c("%-25s", "%-5s", "%-12s", "%-12s", "%-12s")
    }

    # Print table header
    cat(separator, "\n")
    cat(paste(
      sprintf(col_format[1], colnames_df[1]), " | ",
      sprintf(col_format[2], colnames_df[2]), " | ",
      if (ncol == 3) sprintf(col_format[3], colnames_df[3]),
      if (ncol == 5) paste0(sprintf(col_format[3], colnames_df[3]), " | "),
      if (ncol == 5) paste0(sprintf(col_format[4], colnames_df[4]), " | "),
      if (ncol == 5) sprintf(col_format[5], colnames_df[5]),
      sep = ""
    ), "\n")
    cat(separator, "\n")

    # Print each row of the data frame
    apply(df, 1, function(row) {
      cat(paste(
        sprintf(col_format[1], row[1]), " | ",
        sprintf(col_format[2], row[2]), " | ",
        if (ncol == 3) sprintf(col_format[3], row[3]),
        if (ncol == 5) paste0(sprintf(col_format[3], row[3]), " | "),
        if (ncol == 5) paste0(sprintf(col_format[4], row[4]), " | "),
        if (ncol == 5) sprintf(col_format[5], row[5]),
        sep = ""
      ), "\n")
    })
    cat(separator, "\n")
  }

  # Printer for tables with an arbitrary number of columns. Column widths are
  # derived from the content, so narrow tables stay narrow in the console.
  print_desc_table <- function(df) {
    cells <- lapply(df, function(col) {
      ifelse(is.na(col), "NA", as.character(col))
    })

    widths <- mapply(
      function(col, nm) max(nchar(c(col, nm))),
      cells,
      names(df)
    )

    # 3 characters per separator (" | "), minus the trailing one
    separator <- paste(
      rep("-", sum(widths) + 3L * length(widths) - 3L),
      collapse = ""
    )

    # character columns read better left-aligned, numbers right-aligned
    align <- ifelse(vapply(df, is.numeric, logical(1L)), "", "-")
    col_format <- sprintf("%%%s%ds", align, widths)

    fmt_row <- function(row) {
      paste(mapply(sprintf, col_format, row), collapse = " | ")
    }

    cat(separator, "\n")
    cat(fmt_row(names(df)), "\n")
    cat(separator, "\n")

    for (i in seq_len(nrow(df))) {
      cat(fmt_row(vapply(cells, `[`, character(1L), i)), "\n")
    }

    cat(separator, "\n")
  }

  # Function to print quality statistics tables (mean and max)
  print_quality_table <- function(table, label) {
    cat(label, ":\n")
    colnames_main <- c("Variable", "Coef", "Before", "After", "Quality")
    print_table(
      table,
      colnames_df = colnames_main,
      ncol = 5
    )
    cat("\n")
  }

  # Header for the output
  cat("\nMatching Quality Evaluation\n")
  cat(paste(rep("=", 80), collapse = ""), "\n\n")

  # Print the count table for the treatment variable (3 columns)
  cat("Count table for the treatment variable:\n")
  print_table(
    x$count_table,
    colnames_df = colnames(x$count_table),
    ncol = 3
  )

  # Matching summary statistics
  cat("\n\nMatching summary statistics:\n")
  cat(paste(rep("-", 40), collapse = ""), "\n")
  cat("Total n before matching:\t", format(x$n_before, nsmall = 0), "\n")
  cat("Total n after matching:\t\t", format(x$n_after, nsmall = 0), "\n")
  cat(
    "% of matched observations:\t",
    format(x$perc_matched, nsmall = 2), "%\n"
  )

  # Print summary_head (maximal values for each variable)
  for (i in seq_along(x$summary_head)) {
    sum_text <- "maximal"
    tab_after <- ifelse(
      names(x$summary_head[i]) == "r" && "max" %nin% x$statistic,
      "\t\t",
      "\t"
    )
    cat(
      "Total ", sum_text, " ", names(x$summary_head[i]), "value:",
      tab_after, x$summary_head[i], "\n"
    )
  }

  cat("\n\n")

  # Print mean and/or max quality tables based on the statistic
  if (NROW(x$quality_mean) > 0L) {
    if (length(x$statistic) == 1) {
      if (x$statistic == "mean") {
        print_quality_table(x$quality_mean, "Mean values")
      } else {
        print_quality_table(x$quality_max, "Maximal values")
      }
    } else {
      print_quality_table(x$quality_mean, "Mean values")
      cat("\n")
      print_quality_table(x$quality_max, "Maximal values")
    }
  }

  # Descriptive statistics of the covariates, requested via the `desc_*` types.
  # Continuous and categorical covariates are described by different statistics
  # and are therefore printed as two separate tables.
  if (!is.null(x$quality_desc)) {
    continuous <- x$quality_desc[x$quality_desc$Type == "continuous", ,
      drop = FALSE
    ]
    categorical <- x$quality_desc[x$quality_desc$Type == "categorical", ,
      drop = FALSE
    ]

    if (nrow(continuous) > 0L) {
      cat("Descriptive statistics of the continuous covariates:\n")
      print_desc_table(continuous[, c(
        "Variable", "Group", "Statistic", "Before", "After"
      )])
      cat("\n")
    }

    if (nrow(categorical) > 0L) {
      # the count and the percentage are compressed into a single cell, with
      # the percentage always shown to one decimal to keep the cell narrow
      count_percent <- function(count, percent) {
        ifelse(
          is.na(count),
          "NA",
          sprintf("%s (%.1f%%)", format(count), percent)
        )
      }

      cat("Distribution of the categorical covariates:\n")
      print_desc_table(data.frame(
        Variable = categorical$Variable,
        Group = categorical$Group,
        Level = categorical$Statistic,
        Before = count_percent(
          categorical$Before, categorical$Percent_Before
        ),
        After = count_percent(categorical$After, categorical$Percent_After),
        stringsAsFactors = FALSE
      ))
      cat("\n")
    }
  }

  invisible(x)
}

#' @export
str.quality <- function(object, ...) {
  # attributes from balqual()
  od_before <- attr(object, "original_data_before")
  od_after <- attr(object, "original_data_after")
  smd_df <- attr(object, "smd_df_combo")

  n_before_attr <- object$n_before %||% if (!is.null(od_before)) {
    NROW(od_before)
  } else {
    NA_integer_
  }
  n_after_attr <- object$n_after %||% if (!is.null(od_after)) {
    NROW(od_after)
  } else {
    NA_integer_
  }

  perc <- object$perc_matched %||% if (!is.na(n_before_attr) &&
    n_before_attr > 0) {
    100 * n_after_attr / n_before_attr
  } else {
    NA_real_
  }

  types <- object$type %||% NA_character_
  stats <- object$statistic %||% NA_character_
  cutoffs <- object$cutoffs %||% NA_real_
  q_mean <- object$quality_mean
  q_max <- object$quality_max
  q_desc <- object$quality_desc
  ct <- object$count_table

  # treatment levels from count_table, if available
  tr_levels <- if (!is.null(ct) && "Treatment" %in% names(ct)) {
    unique(as.character(ct$Treatment))
  } else {
    NULL
  }

  ## Header --------------------------------------------------------------------
  cat("quality object: matching diagnostics\n")

  cat(sprintf(
    " Observations (treatment-level panel): before = %s, after = %s",
    if (!is.na(n_before_attr)) n_before_attr else "NA",
    if (!is.na(n_after_attr)) n_after_attr else "NA"
  ))
  cat("\n")

  if (!is.na(perc)) {
    cat(sprintf(" Overall retention: %.2f%%\n", perc))
  } else {
    cat(" Overall retention: NA\n")
  }

  # metrics & summary types
  cat(" Metrics (type): ")
  if (all(is.na(types))) {
    cat("NA\n")
  } else {
    cat(paste(types, collapse = ", "), "\n")
  }

  cat(" Summary statistics: ")
  if (all(is.na(stats))) {
    cat("NA\n")
  } else {
    cat(paste(stats, collapse = ", "), "\n")
  }

  # cutoffs
  cat(" Cutoffs: ")
  if (all(is.na(cutoffs))) {
    cat("NA\n")
  } else {
    cat(paste(cutoffs, collapse = ", "), "\n")
  }

  # treatment levels
  if (!is.null(tr_levels)) {
    cat(" Treatment levels (from count_table): ",
      paste(tr_levels, collapse = ", "),
      "\n",
      sep = ""
    )
  }

  # sizes of key tables
  if (!is.null(q_mean)) {
    cat(sprintf(
      " quality_mean: %d x %d (variables x metrics)\n",
      NROW(q_mean), NCOL(q_mean)
    ))
  } else {
    cat(" quality_mean: <none>\n")
  }

  if (!is.null(q_max)) {
    cat(sprintf(
      " quality_max : %d x %d (variables x metrics)\n",
      NROW(q_max), NCOL(q_max)
    ))
  } else {
    cat(" quality_max : <none>\n")
  }

  if (!is.null(q_desc)) {
    cat(sprintf(
      " quality_desc: %d x %d (variables x groups x stages)\n",
      NROW(q_desc), NCOL(q_desc)
    ))
  } else {
    cat(" quality_desc: <none>\n")
  }

  if (!is.null(ct)) {
    cat(sprintf(
      " count_table : %d x %d (Treatment x Before/After)\n",
      NROW(ct), NCOL(ct)
    ))
  } else {
    cat(" count_table : <none>\n")
  }

  if (!is.null(smd_df)) {
    cat(sprintf(
      " smd_df_combo: %d x %d (pairwise combinations x covariates)\n",
      NROW(smd_df), NCOL(smd_df)
    ))
  } else {
    cat(" smd_df_combo: <none>\n")
  }

  cat("\nUnderlying list structure:\n")
  utils::str(unclass(object), ...)
  invisible(object)
}