R/print.model_set.R

In modelbpp: Model BIC Posterior Probability

#' @title Print a `model_set`-Class Object
#'
#' @description Print the content of
#' a `model_set`-class object.
#'
#' @details
#' It is the print method of the
#' output of [model_set()].
#'
#' @return
#' `x` is returned invisibly.
#' Called for its side effect.
#'
#' @param x A `model_set`-class object.
#'
#' @param bic_digits The number of
#' decimal places to be displayed
#' for BIC.
#' Default is 3.
#'
#' @param bpp_digits The number of
#' decimal places to be displayed
#' for BIC posterior probability
#' and prior probabilities.
#' Default is 3.
#'
#' @param sort_models Whether the models
#' will be sorted by BIC posterior
#' probability.
#' Default is `TRUE`.
#'
#' @param max_models The maximum number
#' of models to be printed. Default is
#' 20.
#'
#' @param bpp_target The desired
#' BIC probability. Used to compute
#' and print
#' the minimum prior probability
#' of the target model required to
#' achieve `bpp_target`. Default
#' is `NULL`.
#'
#' @param target_name The name of the
#' target model as appeared in the
#' model list. Default is `"original"`.
#' Used if `bpp_target` is not `NULL`.
#'
#' @param more_fit_measures Character
#' vector. To be passed to
#' [lavaan::fitMeasures()]. Default
#' is `c("cfi", "rmsea")`. Set it to
#' `NULL` to disable printing additional
#' fit measures.
#'
#' @param fit_measures_digits The number of
#' decimal places to be displayed
#' for additional fit measures, if
#' requested. Default is 3.
#'
#' @param short_names If `TRUE`,
#' then simple short names will be
#' printed
#' along with full model names.
#' Default is `FALSE`. Short names
#' can be used when interpreting
#' the graph from `model_graph()`
#' if short names are used in the graph.
#'
#' @param cumulative_bpp If `TRUE` and
#' the models are sorted by BPPs,
#' cumulative BPPs will be printed.
#' Default is `FALSE`.
#'
#' @param ...  Optional arguments.
#' Ignored.
#'
#' @author Shu Fai Cheung <https://orcid.org/0000-0002-9871-9448>
#'
#' @seealso A `model_set`-class object
#' is generated by [model_set()].
#'
#' @examples
#'
#' library(lavaan)
#'
#' dat <- dat_path_model
#'
#' mod <-
#' "
#' x3 ~ a*x1 + b*x2
#' x4 ~ a*x1
#' ab := a*b
#' "
#'
#' fit <- sem(mod, dat_path_model, fixed.x = TRUE)
#'
#' out <- model_set(fit)
#' out
#'
#' @export

print.model_set <- function(x,
                            bic_digits = 3,
                            bpp_digits = 3,
                            sort_models = TRUE,
                            max_models = 20,
                            bpp_target = NULL,
                            target_name = "original",
                            more_fit_measures = c("cfi", "rmsea"),
                            fit_measures_digits = 3,
                            short_names = FALSE,
                            cumulative_bpp = FALSE,
                            ...) {
    fit_n <- length(x$models)
    fit_names <- names(x$models)
    models_fitted <- !is.null(x$fit)
    model_set_call <- x$model_set_call
    if (!models_fitted) {
        fit_many_call <- NULL
        k_converged <- NA
        k_post_check <- NA
        all_converged <- NA
        all_post_checked <- NA
      } else {
        fit_many_call <- x$call
        k_converged <- sum(sapply(x$converged, isTRUE))
        k_post_check <- sum(sapply(x$post_check, isTRUE))
        all_converged <- all(x$converged)
        all_post_checked <- all(x$post_check)
      }
    if (!models_fitted) {
        change_tmp <- rep(NA, fit_n)
        model_df_tmp <- rep(NA, fit_n)
        prior_tmp <- rep(NA, fit_n)
        bic_tmp <- rep(NA, fit_n)
        postprob_tmp <- rep(NA, fit_n)
      } else {
        change_tmp <- x$change
        if (!is.null(x$model_df)) {
            model_df_tmp <- x$model_df
          } else {
            model_df_tmp <- sapply(x$fit, lavaan::fitMeasures,
                                  fit.measures = "df")
          }
        model_df_tmp <- unname(model_df_tmp)
        prior_tmp <- x$prior
        bic_tmp <- x$bic
        postprob_tmp <- x$bpp
      }
    out_table <- data.frame(modification = fit_names,
                            model_df = model_df_tmp,
                            df_diff = change_tmp)
    if (models_fitted && !all_converged) {
        out_table$Converged <- ifelse(x$converged,
                                      "Yes",
                                      "No")
      }
    if (models_fitted && !all_post_checked) {
        out_table$Check <- ifelse(x$post_check,
                                      "Passed",
                                      "Failed")
      }
    out_table$Prior <- prior_tmp
    out_table$BIC <- bic_tmp
    out_table$BPP <- postprob_tmp
    if (!is.null(more_fit_measures)) {
        fit_fm <- sapply(x$fit,
                         function(xx) {
                            out <- tryCatch(lavaan::fitMeasures(xx,
                                            fit.measures = more_fit_measures,
                                            output = "vector"),
                                            error = function(e) e)
                            if (inherits(out, "error")) {
                                out <- rep(NA, length(more_fit_measures))
                              }
                            out
                          })
        fit_fm <- t(fit_fm)
        out_table <- cbind(out_table, fit_fm)
      }
    if (sort_models && models_fitted && all_converged) {
        i <- order(out_table$BPP,
                   decreasing = TRUE)
        out_table <- out_table[i, ]
        if (cumulative_bpp) {
            tmp <- round(cumsum(out_table$BPP), bpp_digits)
            tmp <- formatC(tmp,
                          digits = bpp_digits,
                          format = "f")
            out_table["Cumulative"] <- tmp
          }
      }
    out_table_print <- out_table
    out_table_print$Prior <- round(out_table_print$Prior,
                                 digits = bpp_digits)
    out_table_print$Prior <- formatC(out_table_print$Prior,
                                   digits = bpp_digits,
                                   format = "f")
    out_table_print$BIC <- round(out_table_print$BIC,
                                 digits = bic_digits)
    out_table_print$BIC <- formatC(out_table_print$BIC,
                                   digits = bic_digits,
                                   format = "f")
    out_table_print$BPP <- round(out_table_print$BPP,
                                 digits = bpp_digits)
    out_table_print$BPP <- formatC(out_table_print$BPP,
                                   digits = bpp_digits,
                                   format = "f")
    if (!is.null(more_fit_measures)) {
        fm_names <- colnames(fit_fm)
        for (xx in fm_names) {
            out_table_print[xx] <- formatC(out_table_print[, xx, drop = TRUE],
                                           digits = fit_measures_digits,
                                           format = "f")
          }
      }
    if (short_names) {
        if (is.character(x$short_names)) {
            tmp1 <- x$short_names
            tmp2 <- out_table_print$modification
            if (isTRUE(setequal(names(tmp1),
                                tmp2))) {
                xx <- unname(tmp1[tmp2])
                out_table_print <- cbind(Short = xx,
                                         out_table_print)
              }
          }
      }
    if (isTRUE(fit_n > max_models)) {
        gt_max_models <- TRUE
        x_tmp <- out_table_print[seq_len(max_models), ]
      } else {
        gt_max_models <- FALSE
        x_tmp <- out_table_print
      }
    cat("\n")
    cat("Call:\n")
    print(model_set_call)
    cat("\n")
    if (models_fitted) {
        tmp1 <- c("Number of model(s) fitted",
                  "Number of model(s) converged",
                  "Number of model(s) passed post.check")
        tmp2 <- c(fit_n,
                  k_converged,
                  k_post_check)
        tmp3 <- auto_tab(tmp1,
                         tmp2,
                         between = ": ")
        cat(tmp3, sep = "\n")
        if (!is.null(bpp_target)) {
            bpp_min <- min_prior(x$bic,
                                 bpp_target = bpp_target,
                                 target_name = target_name)
            tmp <- data.frame(x = c(
                      formatC(bpp_target, digits = bpp_digits, format = "f"),
                      formatC(bpp_min, digits = bpp_digits, format = "f"),
                      formatC(x$bpp[target_name], digits = bpp_digits, format = "f")))
            colnames(tmp) <- paste0("Target Model: ",
                                    target_name)
            rownames(tmp) <- c("Desired minimum BPP:",
                               "Required minimum prior probability:",
                               "Current BPP:")
            print(tmp)
          }
      } else {
        cat("Models are not fitted.")
        cat("\n")
      }
    cat("\n")
    tmp1 <- ifelse(sort_models && all(!is.na(out_table$BPP)),
                   " (sorted by BPP)",
                   "")
    if (!models_fitted) tmp1 <- ""
    if (gt_max_models) {
        cat("The first ",
            max_models,
            " model(s)",
            tmp1,
            ":\n", sep = "")
      } else {
        cat("The models",
            tmp1,
            ":\n", sep = "")
      }
    rownames(x_tmp) <- x_tmp$modification
    x_tmp$modification <- NULL
    print(x_tmp)

    if (models_fitted && !all_converged) {
        x_tmp2 <- out_table_print[out_table_print$Converged != "Yes", ]
        rownames(x_tmp2) <- x_tmp2$modification
        x_tmp2$modification <- NULL
        cat("\nModel(s) not converged:\n")
        print(x_tmp2)
      }

    if (models_fitted && !all_post_checked) {
        x_tmp2 <- out_table_print[out_table_print$Check != "Passed", ]
        rownames(x_tmp2) <- x_tmp2$modification
        x_tmp2$modification <- NULL
        cat("\nModel(s) failed lavaan's post.check:\n")
        print(x_tmp2)
      }

    if (!is.null(x$equivalent_clusters)) {
        tmp <- as.vector(sapply(x$equivalent_clusters, paste, collapse = ", "))
        tmp <- data.frame(Cluster = tmp)
        tmp2 <- names(x$models_equivalent)
        cat("\nModels that are equivalent:\n")
        print(tmp, right = FALSE)
        cat("\nEquivalent model(s) excluded from the analysis:\n")
        catwrap(paste(tmp2, collapse = ", "))
      }

    cat("\nNote:\n")
    cat("- BIC: Bayesian Information Criterion.\n")
    cat("- BPP: BIC posterior probability.\n")
    cat("- model_df: Model degrees of freedom.\n")
    cat("- df_diff: Difference in df compared to the original/target model.\n")
    if (sort_models) {
        if (("Cumulative" %in% colnames(x_tmp))) {
            cat("- Cumulative: Cumulative BIC posterior probability.\n")
          } else {
            cat("- To show cumulative BPPs, call print() with 'cumulative_bpp = TRUE'.\n")
          }
      }
    if (gt_max_models) {
        tmp <- paste(fit_n,
                   "models were fitted but",
                   max_models,
                   "were printed. Call print() and",
                   "set 'max_models' to a larger number",
                   "to print more models, or set it to",
                   "NA to print all models.")
        catwrap(tmp, initial = "- ", exdent = 2)
      }
    if (models_fitted &&
        (k_converged != fit_n) &&
        any(is.na(postprob_tmp))) {
        tmp <- "BPP and/or prior not computed because one or more models not converged."
        catwrap(tmp, initial = "- ", exdent = 2)
      }
    if (models_fitted &&
        (k_post_check != fit_n)) {
        tmp <- "Interpret with caution. One or more models failed lavaan's post.check."
        catwrap(tmp, initial = "- ", exdent = 2)
      }
    if (x$drop_equivalent_models && x$fixed_x) {
        tmp <- "At least one model has fixed.x = TRUE. The models are not checked for equivalence."
        catwrap(tmp, initial = "- ", exdent = 2)
      }
    tmp <- paste0("Since Version 0.1.3.5, the default values of ",
                  "exclude_feedback and exclude_xy_cov changed to TRUE. ",
                  "Set them to FALSE to reproduce results from previous versions.")
    catwrap(tmp, initial = "- ", exdent = 2)

    invisible(x)
  }