R/summary.JointAI.R

In JointAI: Joint Analysis and Imputation of Incomplete Data

#' Summarize the results from an object of class JointAI
#'
#' Obtain and print the \code{summary}, (fixed effects) coefficients
#' (\code{coef}) and credible interval (\code{confint}) for an object of
#' class 'JointAI'.
#'
#' @param digits the minimum number of significant digits to be printed in
#'               values.
#' @param quantiles posterior quantiles
#' @inheritParams sharedParams
#' @param outcome optional; vector identifying for which outcomes the summary
#'                should be given, either by specifying their indices, or their
#'                names (LHS of the respective model formulas as character
#'                string).
#' @param missinfo logical; should information on the number and proportion of
#'                 missing values be included in the summary?
#' @param \dots currently not used
#'
#' @examples
#'
#' \dontrun{
#' mod1 <- lm_imp(y ~ C1 + C2 + M2, data = wideDF, n.iter = 100)
#'
#' summary(mod1, missinfo = TRUE)
#' coef(mod1)
#' confint(mod1)
#' }
#'
#' @seealso The model fitting functions \code{\link{lm_imp}},
#'          \code{\link{glm_imp}}, \code{\link{clm_imp}}, \code{\link{lme_imp}},
#'          \code{\link{glme_imp}}, \code{\link{survreg_imp}} and
#'          \code{\link{coxph_imp}},
#'          and the vignette
#'          \href{https://nerler.github.io/JointAI/articles/SelectingParameters.html}{Parameter Selection}
#'          for examples how to specify the parameter \code{subset}.
#'
#' @export
summary.JointAI <- function(
  object,
  start = NULL,
  end = NULL,
  thin = NULL,
  quantiles = c(0.025, 0.975),
  subset = NULL,
  exclude_chains = NULL,
  outcome = NULL,
  missinfo = FALSE,
  warn = TRUE,
  mess = TRUE,
  ...
) {
  if (is.null(object$MCMC)) {
    errormsg("There is no MCMC sample.")
  }

  cl <- as.list(match.call())[-1]
  autoburnin <- if (is.null(cl$autoburnin)) FALSE else eval(cl$autoburnin)

  MCMC <- prep_MCMC(
    object,
    start = start,
    end = end,
    thin = thin,
    subset = subset,
    exclude_chains = exclude_chains,
    warn = warn,
    mess = mess,
    ...
  )

  # create results matrices
  statnames <- c(
    "Mean",
    "SD",
    paste0(quantiles * 100, "%"),
    "tail-prob.",
    "GR-crit",
    "MCE/SD"
  )

  vars <- if (is.null(outcome)) {
    names(object$coef_list)
  } else {
    names(object$fixed[outcome])
  }
  res_list <- sapply(
    vars,
    function(varname) {
      rdnam <- nlapply(
        names(object$info_list[[varname]]$rd_vcov),
        function(lvl) {
          if (isTRUE(object$info_list[[varname]]$rd_vcov[[lvl]] == "full")) {
            paste0(
              "^",
              c("b", "D", "invD", "RinvD", "KinvD"),
              attr(object$info_list[[varname]]$rd_vcov[[lvl]], "name"),
              "_",
              lvl,
              "\\b"
            )
          }
        }
      )

      MCMCsub <- MCMC[,
        intersect(
          colnames(MCMC),
          c(
            object$coef_list[[varname]]$coef,
            grep(
              paste0("_", object$info_list[[varname]]$varname, "\\b"),
              colnames(MCMC),
              value = TRUE
            ),
            grep(
              paste0("_", object$info_list[[varname]]$varname, "_"),
              colnames(MCMC),
              value = TRUE
            ),
            unlist(lapply(unlist(rdnam), grep, colnames(MCMC), value = TRUE))
          )
        ),
        drop = FALSE
      ]

      if (ncol(MCMCsub) > 0) {
        grcrit <- if (length(object$MCMC) - length(exclude_chains) > 1) {
          GR_crit(
            object = object,
            start = start,
            end = end,
            thin = thin,
            warn = warn,
            mess = FALSE,
            multivariate = FALSE,
            exclude_chains = exclude_chains,
            subset = list(other = colnames(MCMCsub), analysis_main = FALSE),
            autoburnin = autoburnin
          )[[1]][, "Upper C.I."]
        }

        mcerror <- if (length(object$MCMC) - length(exclude_chains) > 1) {
          try(MC_error(
            object,
            subset = list(other = colnames(MCMCsub), analysis_main = FALSE),
            exclude_chains = exclude_chains,
            start = start,
            end = end,
            thin = thin,
            digits = 2,
            warn = FALSE,
            mess = FALSE
          ))
        }

        colnames(MCMCsub)[na.omit(match(
          object$coef_list[[varname]]$coef,
          colnames(MCMCsub)
        ))] <-
          object$coef_list[[varname]]$varnam_print

        stats <- matrix(
          nrow = length(colnames(MCMCsub)),
          ncol = length(statnames),
          dimnames = list(colnames(MCMCsub), statnames)
        )

        stats[, "Mean"] <- apply(MCMCsub, 2, mean)
        stats[, "SD"] <- apply(MCMCsub, 2, sd)
        stats[, paste0(quantiles * 100, "%")] <- t(apply(
          MCMCsub,
          2,
          quantile,
          quantiles
        ))
        stats[, "tail-prob."] <- apply(MCMCsub, 2, computeP)

        if (length(object$MCMC) - length(exclude_chains) > 1) {
          stats[, "GR-crit"] <- grcrit
        }

        if (length(object$MCMC) - length(exclude_chains) > 1) {
          if (!inherits(mcerror, "try-error")) {
            stats[, "MCE/SD"] <- mcerror$data_scale[, "MCSE/SD"]
          }
        }

        regcoef <- stats[
          intersect(rownames(stats), object$coef_list[[varname]]$varnam_print),
          ,
          drop = FALSE
        ]

        sigma <- if (
          object$info_list[[varname]]$family %in%
            c("gaussian", "Gamma", "lognorm") &&
            !is.null(object$info_list[[varname]]$family)
        ) {
          sig <- grep(paste0("sigma_", varname), rownames(stats))

          if (length(sig) > 0) {
            stats[sig, -which(colnames(stats) == "tail-prob."), drop = FALSE]
          }
        }

        intercepts <- if (
          object$info_list[[varname]]$modeltype %in%
            c("clm", "clmm")
        ) {
          get_intercepts(
            stats,
            varname,
            levels(object$Mlist$refs[[varname]]),
            rev = object$info_list[[varname]]$rev
          )
        }

        rd_vcov <- if (!is.null(object$info_list[[varname]]$hc_list)) {
          Dpat <- nlapply(object$Mlist$idvar, function(lvl) {
            if (isTRUE(object$info_list[[varname]]$nranef[lvl] > 0L)) {
              if (
                isTRUE(object$info_list[[varname]]$rd_vcov[[lvl]] == "full")
              ) {
                paste0(
                  "^D[[:digit:]]*_",
                  lvl,
                  "\\[[[:digit:]]+,[[:digit:]]+\\]"
                )
              } else {
                paste0(
                  "^D_",
                  object$info_list[[varname]]$varname,
                  "_",
                  lvl,
                  "\\[[[:digit:]]+,[[:digit:]]+\\]"
                )
              }
            }
          })

          Dpat <- Filter(Negate(is.null), Dpat)

          Ds <- nlapply(names(Dpat), function(lvl) {
            D <- stats[
              grep(Dpat[[lvl]], rownames(stats), value = TRUE),
              ,
              drop = FALSE
            ]

            if (nrow(D) > 0) {
              Ddiag <- sapply(
                strsplit(
                  sub("\\]", "", sub("^[[:print:]]*\\[", "", rownames(D))),
                  ","
                ),
                function(i) length(unique(i)) == 1
              )

              D[Ddiag, "tail-prob."] <- NA
              attr(D, "warnings") <- attr(
                object$info_list[[varname]]$hc_list,
                "warnings"
              )
              attr(D, "rd_vcov") <- object$info_list[[varname]]$rd_vcov[[lvl]]

              for (k in which(names(object$Mlist$rd_vcov[[lvl]]) == "full")) {
                if (varname %in% object$Mlist$rd_vcov[[lvl]][[k]]) {
                  attr(D, "ranef_index") <- attr(
                    object$Mlist$rd_vcov[[lvl]][[k]],
                    "ranef_index"
                  )
                }
              }
            }
            D
          })
        }

        assoc_type <- if (object$info_list[[varname]]$modeltype %in% "JM") {
          object$info_list[[which(
            sapply(object$info_list, "[[", "modeltype") == "JM"
          )]]$assoc_type
        }

        wb_shape <- if (
          object$info_list[[varname]]$modeltype %in% c("survreg")
        ) {
          stats[
            c(paste0("shape_", object$info_list[[varname]]$varname)),
            -which(colnames(stats) == "tail-prob."),
            drop = FALSE
          ]
        }

        events <- if (
          object$info_list[[varname]]$modeltype %in%
            c("survreg", "coxph", "JM")
        ) {
          mat <- object$info_list[[1]]$resp_mat[2]
          col <- object$info_list[[1]]$resp_col[2]
          sum(object$data_list[[mat]][, col])
        }

        other <- setdiff(
          rownames(stats),
          c(
            rownames(regcoef),
            rownames(sigma),
            attr(intercepts, "rownames_orig"),
            if (is.list(rd_vcov)) {
              rownames(do.call(rbind, rd_vcov))
            } else {
              rownames(rd_vcov)
            },
            rownames(wb_shape)
          )
        )

        otherpars <- if (length(other) > 0) {
          stats[other, , drop = FALSE]
        }

        list(
          modeltype = object$info_list[[varname]]$modeltype,
          family = object$info_list[[varname]]$family,
          regcoef = regcoef,
          sigma = sigma,
          intercepts = intercepts,
          rd_vcov = rd_vcov,
          wb_shape = wb_shape,
          assoc_type = assoc_type,
          events = events,
          grcrit = grcrit,
          otherpars = otherpars
        )
      }
    },
    simplify = FALSE
  )

  out <- list()
  out$call <- object$call
  out$start <- ifelse(
    is.null(start),
    start(object$MCMC),
    max(start, start(object$MCMC))
  )
  out$end <- ifelse(is.null(end), end(object$MCMC), min(end, end(object$MCMC)))
  out$thin <- ifelse(is.null(thin), coda::thin(object$MCMC), thin)
  out$nchain <- coda::nchain(object$MCMC) -
    sum(
      exclude_chains %in%
        seq_along(object$MCMC)
    )
  out$res <- res_list
  out$outcome <- outcome
  out$missinfo <- if (missinfo) get_missinfo(object)

  out$analysis_type <- object$analysis_type
  out$size <- object$Mlist$N

  class(out) <- "summary.JointAI"
  return(out)
}


#' @rdname summary.JointAI
#' @param x an object of class \code{summary.JointAI} or \code{JointAI}
#' @export
print.summary.JointAI <- function(
  x,
  digits = max(3, .Options$digits - 4),
  ...
) {
  if (!inherits(x, "summary.JointAI")) {
    errormsg("Use only with objects.", sQuote("summary.JointAI"))
  }

  cat("\n")

  if (sum(!sapply(x$res, is.null)) > 1 | !is.null(x$outcome)) {
    cat("Bayesian joint model fitted with JointAI", "\n")
  } else {
    cat(
      "Bayesian",
      print_type(x$res[[1]]$modeltype, x$res[[1]]$family),
      "fitted with JointAI\n"
    )
  }

  cat(
    "\nCall:\n",
    paste(deparse(x$call), sep = "\n", collapse = "\n"),
    "\n",
    sep = ""
  )

  for (k in seq_along(x$res)) {
    if (!is.null(x$res[[k]])) {
      cat("\n\n")
      if (sum(!sapply(x$res, is.null)) > 1 | !is.null(x$outcome)) {
        cat(paste0(
          "# ",
          paste0(c(rep("-", 69)), collapse = ""),
          " #\n",
          "  ",
          "Bayesian ",
          print_type(x$res[[k]]$modeltype, x$res[[k]]$family),
          " for ",
          dQuote(names(x$res)[k]),
          "\n",
          "# ",
          paste0(c(rep("-", 35)), collapse = " "),
          " #\n\n"
        ))
      }

      if (!is.null(x$res[[k]]$events)) {
        cat("Number of events:", x$res[[k]]$events, "\n\n")
      }

      if (!is.null(x$res[[k]]$regcoef)) {
        cat("Posterior summary:\n")
        print(x$res[[k]]$regcoef, digits = digits, na.print = "")
      }

      if (!is.null(x$res[[k]]$intercepts)) {
        # remove the attributes to avoid printing them
        attr(x$res[[k]]$intercepts, "rownames_orig") <- NULL

        cat("\nPosterior summary of the intercepts:\n")
        print(x$res[[k]]$intercepts, digits = digits, na.print = "")
      }

      if (!is.null(x$res[[k]]$rd_vcov)) {
        cat("\n\nPosterior summary of random effects covariance matrix:\n")
        for (lvl in names(x$res[[k]]$rd_vcov)) {
          warnings <- attr(x$res[[k]]$rd_vcov[[lvl]], "warnings")
          rd_vcov <- attr(x$res[[k]]$rd_vcov[[lvl]], "rd_vcov")
          ranef_index <- attr(x$res[[k]]$rd_vcov[[lvl]], "ranef_index")
          attr(x$res[[k]]$rd_vcov[[lvl]], "warnings") <- NULL
          attr(x$res[[k]]$rd_vcov[[lvl]], "rd_vcov") <- NULL
          attr(x$res[[k]]$rd_vcov[[lvl]], "ranef_index") <- NULL

          if (length(names(x$res[[k]]$rd_vcov)) > 1) {
            cat(paste0("\n* For level ", dQuote(lvl), ":", "\n"))
          }

          w <- if (isTRUE(rd_vcov == "full")) {
            which(names(x$res)[k] == names(ranef_index))
          }

          print_d <- isTRUE(rd_vcov == "full") &&
            (w == 1 ||
              if (w > 1) {
                all(lvapply(x$res[names(ranef_index)[seq_len(w - 1)]], is.null))
              })

          if (print_d || isTRUE(rd_vcov != "full")) {
            if (!is.null(ranef_index)) {
              ranef_index <- lapply(ranef_index, function(nr) {
                eval(parse(text = nr))
              })

              cat(
                paste0("\r", tab(), "Indices:"),
                paste0(
                  names(ranef_index),
                  ": ",
                  paste0(ranef_index),
                  collapse = "; "
                ),
                "\n"
              )
            }
            print(x$res[[k]]$rd_vcov[[lvl]], digits = digits, na.print = "")

            if (!is.null(unlist(warnings))) {
              warnmsg(warnings)
            }
          } else {
            w_ref <- min(which(
              !lvapply(x$res[names(ranef_index)[seq_len(w - 1)]], is.null)
            ))
            cat(paste0(
              tab(),
              "(see model for ",
              dQuote(names(ranef_index)[w_ref]),
              ")"
            ))
          }
        }
        cat("\n")
      }

      if (!is.null(x$res[[k]]$sigma)) {
        cat("\nPosterior summary of residual std. deviation:\n")
        print(x$res[[k]]$sigma, digits = digits, na.print = "")
      }

      if (!is.null(x$res[[k]]$wb_shape)) {
        cat("\nPosterior summary of the shape of the Weibull distribution:\n")
        print(x$res[[k]]$wb_shape, digits = digits, na.print = "")
      }

      if (!is.null(x$res[[k]]$assoc_type)) {
        cat("\nAssociation types:\n")
        cat(
          paste0(
            names(x$res[[k]]$assoc_type),
            ": ",
            sapply(x$res[[k]]$assoc_type, function(i) {
              switch(
                i,
                "underl.value" = "underlying value",
                "obs.value" = "observed value"
              )
            }),
            collapse = "\n"
          ),
          "\n"
        )
      }

      if (!is.null(x$res[[k]]$otherpars)) {
        cat("\nPosterior summary of other parameters:\n")
        print(x$res[[k]]$otherpars, digits = digits, na.print = "")
      }
    }
  }

  cat("\n\n")
  if (sum(!sapply(x$res, is.null)) > 1) {
    cat("#", paste0(c(rep("-", 59)), collapse = ""), "#\n\n")
  }

  cat("MCMC settings:\n")
  cat("Iterations = ", x$start, ":", x$end, "\n", sep = "")
  cat("Sample size per chain =", floor((x$end - x$start + 1) / x$thin), "\n")
  cat("Thinning interval =", x$thin, "\n")
  cat("Number of chains =", x$nchain, "\n")
  cat("\n")
  cat("Number of observations:", x$size["lvlone"], "\n")
  if (length(x$size) > 1) {
    i <- which(!names(x$size) %in% "lvlone")
    cat(
      "Number of groups:\n",
      paste0("- ", names(x$size)[i], ": ", x$size[i], "\n")
    )
  }

  if (!is.null(x$missinfo)) {
    cat("\n\n")
    cat("Number and proportion of complete cases:\n")
    print(x$missinfo$complete_cases, digits = digits)
    cat("\nNumber and proportion of missing values:\n")
    for (k in seq_along(x$missinfo$miss_list)) {
      print(x$missinfo$miss_list[[k]], digits = digits)
      cat("\n")
    }
  }

  invisible(x)
}


#' @export
formula.JointAI <- function(x, ...) {
  if (!(inherits(x, "JointAI") | inherits(x, "JointAI_errored"))) {
    errormsg("Use only %s with objects.", sQuote("JointAI"))
  }

  x$formula

  #   if (inherits(x$call, "call")) {
  #   if (is.null(x$call$formula)) {
  #     as.formula(x$call$fixed)
  #   } else {
  #     fmla <- eval(x$call$formula)
  #     if (inherits(eval(x$call$formula), "list")) {
  #       x$call$formula
  #     } else {
  #       as.formula(x$call$formula)
  #     }
  #   }
  # } else if (inherits(x$call, "list")) {
  #   fmla_list <- lapply(x$call, function(k) {
  #     if (!is.null(k$formula)) {
  #       k$formula
  #     } else {
  #       k$fixed
  #     }
  #   })
  #   fmla_list <- lapply(fmla_list[lvapply(fmla_list, inherits, "call")], as.formula)
  #   if (length(fmla_list) == 1L) {
  #     fmla_list[[1]]
  #   } else {
  #     fmla_list
  #   }
  # }
}


#' @rdname summary.JointAI
#' @export
coef.JointAI <- function(
  object,
  start = NULL,
  end = NULL,
  thin = NULL,
  subset = NULL,
  exclude_chains = NULL,
  warn = TRUE,
  mess = TRUE,
  ...
) {
  if (!inherits(object, "JointAI")) {
    errormsg("Use only with %s objects.", sQuote("JointAI"))
  }

  if (is.null(object$MCMC)) {
    errormsg("There is no MCMC sample.")
  }

  MCMC <- prep_MCMC(
    object,
    start,
    end,
    thin,
    subset,
    exclude_chains = exclude_chains,
    mess = mess,
    warn = warn
  )

  params <- parameters(object)

  nlapply(names(object$fixed), function(k) {
    x <- params[params$outcome == internal_clean_survname(k), , drop = FALSE]
    rev <- object$info_list[[k]]$rev

    cols <- unlist(
      lapply(x$coef, function(var) {
        grep(
          paste0(
            glob2rx(var),
            "|^",
            gsub("\\[", "\\\\[", var),
            "\\[[[:digit:]]+\\]$"
          ),
          colnames(MCMC),
          value = TRUE
        )
      })
    )

    cfs <- colMeans(MCMC[, cols, drop = FALSE])

    # replace the regression coefficient parameters with the corresponding
    # variable names
    reg_coefs <- intersect(names(cfs), x$coef[!is.na(x$varname)])
    names(cfs)[match(reg_coefs, names(cfs))] <- x$varname[match(
      reg_coefs,
      x$coef
    )]

    if (object$info_list[[k]]$modeltype %in% c("clm", "clmm")) {
      interc <- grep(paste0("gamma_", k, "\\["), names(cfs))

      lvl <- levels(object$Mlist$refs[[k]])
      if (isTRUE(rev)) {
        names(cfs)[interc] <- paste(k, "\u2264", lvl[-length(lvl)])
      } else {
        names(cfs)[interc] <- paste(k, ">", lvl[-length(lvl)])
      }
    }
    cfs
  })
}

#' @export
coef.summary.JointAI <- function(
  object,
  start = NULL,
  end = NULL,
  thin = NULL,
  subset = NULL,
  exclude_chains = NULL,
  warn = TRUE,
  mess = TRUE,
  ...
) {
  if (!inherits(object, "summary.JointAI")) {
    errormsg("Use only with %s objects.", sQuote("summary.JointAI"))
  }

  Filter(Negate(is.null), lapply(object$res, "[[", "regcoef"))
}


#' @rdname summary.JointAI
#' @param parm same as \code{subset} (for consistency with \code{confint}
#'             method for other types of objects)
#' @param level confidence level (default is 0.95)
#' @export
confint.JointAI <- function(
  object,
  parm = NULL,
  level = 0.95,
  quantiles = NULL,
  start = NULL,
  end = NULL,
  thin = NULL,
  subset = NULL,
  exclude_chains = NULL,
  warn = TRUE,
  mess = TRUE,
  ...
) {
  if (!inherits(object, "JointAI")) {
    errormsg("Use only with %s objects.", sQuote("JointAI"))
  }

  if (is.null(object$MCMC)) {
    errormsg("There is no MCMC sample.")
  }

  if (is.null(subset) & !is.null(parm)) {
    subset <- parm
  }

  if (!is.null(subset) & !is.null(parm)) {
    errormsg(
      "At least one of %s and should be NULL.",
      dQuote("parm"),
      dQuote("subset")
    )
  }

  if (is.null(quantiles) & !is.null(level)) {
    quantiles <- c((1 - level) / 2, 1 - (1 - level) / 2)
  }

  MCMC <- prep_MCMC(
    object,
    start,
    end,
    thin,
    subset,
    exclude_chains = exclude_chains,
    mess = mess,
    warn = warn
  )

  params <- parameters(object)

  nlapply(names(object$fixed), function(k) {
    x <- params[params$outcome == internal_clean_survname(k), , drop = FALSE]
    rev <- object$info_list[[k]]$rev

    cols <- unlist(
      lapply(x$coef, function(var) {
        grep(
          paste0(
            glob2rx(var),
            "|^",
            gsub("\\[", "\\\\[", var),
            "\\[[[:digit:]]+\\]$"
          ),
          colnames(MCMC),
          value = TRUE
        )
      })
    )

    MCMC_sub <- MCMC[, cols, drop = FALSE]

    quants <- t(apply(MCMC_sub, 2, quantile, quantiles))

    # replace the regression coefficient parameters with the corresponding
    # variable names
    reg_coefs <- intersect(rownames(quants), x$coef[!is.na(x$varname)])
    rownames(quants)[match(reg_coefs, rownames(quants))] <-
      x$varname[match(reg_coefs, x$coef)]

    if (object$info_list[[k]]$modeltype %in% c("clm", "clmm")) {
      lvl <- levels(object$Mlist$refs[[k]])
      interc <- grep(paste0("gamma_", k, "\\["), rownames(quants))

      if (isTRUE(rev)) {
        rownames(quants)[interc] <- paste(k, "\u2264", lvl[-length(lvl)])
      } else {
        rownames(quants)[interc] <- paste(k, ">", lvl[-length(lvl)])
      }
    }
    quants
  })
}


#' @rdname summary.JointAI
#' @export
print.JointAI <- function(x, digits = max(4, getOption("digits") - 4), ...) {
  if (!inherits(x, "JointAI")) {
    errormsg("Use only with %s objects.", sQuote("JointAI"))
  }

  MCMC <- if (!is.null(x$MCMC)) {
    prep_MCMC(
      x,
      start = NULL,
      end = NULL,
      thin = NULL,
      subset = NULL,
      exclude_chains = NULL,
      mess = TRUE,
      warn = TRUE,
      ...
    )
  }

  cat("\nCall:\n")
  print(x$call)

  if (!is.null(MCMC)) {
    coefs <- coef(x)

    for (k in seq_along(coefs)) {
      varname <- names(coefs)[k]

      intercepts <- if (
        x$info_list[[varname]]$modeltype %in%
          c("clm", "clmm")
      ) {
        stats <- t(t(colMeans(MCMC[, grep(
          paste0("gamma_", varname),
          colnames(MCMC)
        )])))
        get_intercepts(
          stats,
          varname,
          lvls = levels(x$Mlist$refs[[varname]]),
          rev = x$info_list[[varname]]$rev
        )
      }

      cat(
        "\n",
        "Bayesian",
        print_type(
          x$info_list[[varname]]$modeltype,
          x$info_list[[varname]]$family
        ),
        "for",
        dQuote(varname),
        "\n"
      )
      if (
        x$info_list[[names(coefs)[k]]]$modeltype %in%
          c("glmm", "clmm", "mlogitmm")
      ) {
        if (
          !is.null(intercepts) |
            length(coefs[[k]][x$coef_list[[varname]]$varname]) > 0
        ) {
          cat("\nFixed effects:\n")
          print(
            c(
              setNames(c(intercepts), rownames(intercepts)),
              coefs[[k]][x$coef_list[[varname]]$varname]
            ),
            digits = digits
          )
        }

        cat("\n\nRandom effects covariance matrix:\n")
        print(get_Dmat(object = x, varname = varname), digits = digits)
      } else {
        if (
          length(coefs[[k]][x$coef_list[[varname]]$varname] > 0) |
            !is.null(intercepts)
        ) {
          cat("\n\nCoefficients:\n")
          print(
            c(
              setNames(c(intercepts), rownames(intercepts)),
              coefs[[k]][x$coef_list[[varname]]$varname]
            ),
            digits = digits
          )
        }
      }

      if (paste0("sigma_", varname) %in% colnames(MCMC)) {
        cat("\n\nResidual standard deviation:\n")
        print(
          colMeans(MCMC[, paste0("sigma_", names(coefs)[k]), drop = FALSE]),
          digits = digits
        )
      }
    }
  } else {
    cat("\n(The object does not contain an MCMC sample.)")
  }

  invisible(x)
}


#' @export
print.modelstring <- function(x, ...) {
  if (!inherits(x, "modelstring")) {
    errormsg("Use only with %s objects.", sQuote("modelstring"))
  }

  cat(x)
}


#' Obtain a summary of the missing values involved in an object of class JointAI
#'
#' This function returns a \code{data.frame} or a \code{list} of
#' \code{data.frame}s per grouping level. Each of the \code{data.frames}
#' has columns \code{variable}, \code{#NA} (number of missing values) and
#' \code{%NA} (proportion of missing values in percent).
#'
#' @param object object inheriting from class JointAI
#'
#' @export
#'
#' @examples
#' mod <-  lm_imp(y ~ C1 + B2 + C2, data = wideDF, n.iter = 100)
#' get_missinfo(mod)
#'
#'
get_missinfo <- function(object) {
  if (!(inherits(object, "JointAI") | inherits(object, "JointAI_errored"))) {
    errormsg("Use only with 'JointAI' objects.")
  }

  allvars <- all_vars(
    object$fixed,
    object$random,
    object$Mlist$auxvars,
    object$Mlist$timevar
  )

  groups <- object$Mlist$groups

  # data_lvls <- cvapply(object$data[, allvars], check_varlevel, groups = groups)
  data_lvls <- get_datlvls(object$data[, allvars, drop = FALSE], groups)

  complcases <- lapply(names(groups), function(k) {
    cc <- complete.cases(object$data[
      match(unique(groups[[k]]), groups[[k]]),
      names(data_lvls[data_lvls == k])
    ])

    as.data.frame(
      Filter(
        Negate(is.null),
        list(
          level = if (length(object$Mlist$groups) > 1) k,
          "#" = sum(cc),
          "%" = mean(cc) * 100
        )
      ),
      check.names = FALSE,
      row.names = k
    )
  })

  # dat_lvls <- sapply(object$data[allvars], check_varlevel,
  #                    groups = object$Mlist$groups)
  dat_lvls <- get_datlvls(object$data[allvars], object$Mlist$groups)

  miss_list <- sapply(
    unique(dat_lvls),
    function(lvl) {
      subdat <- object$data[
        match(unique(object$Mlist$groups[[lvl]]), object$Mlist$groups[[lvl]]),
        names(dat_lvls)[dat_lvls == lvl],
        drop = FALSE
      ]
      missinfo <- as.data.frame(
        Filter(
          Negate(is.null),
          list(
            level = if (length(unique(dat_lvls)) > 1) lvl,
            "# NA" = colSums(is.na(subdat)),
            "% NA" = colMeans(is.na(subdat)) * 100
          )
        ),
        check.names = FALSE
      )
      missinfo[order(missinfo$`# NA`), ]
    },
    simplify = FALSE
  )

  list("complete_cases" = do.call(rbind, complcases), miss_list = miss_list)
}