Nothing
# Part of the rstanarm package for estimating model parameters
# Copyright (C) 2015, 2016, 2017 Trustees of Columbia University
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#' Print method for stanreg objects
#'
#' The \code{print} method for stanreg objects displays a compact summary of the
#' fitted model. See the \strong{Details} section below for descriptions of the
#' different components of the printed output. For additional summary statistics
#' and diagnostics use the \code{\link[=summary.stanreg]{summary}} method.
#'
#' @export
#' @method print stanreg
#' @templateVar stanregArg x
#' @template args-stanreg-object
#' @param detail Logical, defaulting to \code{TRUE}. If \code{FALSE} a more
#' minimal summary is printed consisting only of the parameter estimates.
#' @param digits Number of digits to use for formatting numbers.
#' @param ... Ignored.
#' @return Returns \code{x}, invisibly.
#' @details
#' \subsection{Point estimates}{
#' Regardless of the estimation algorithm, point estimates are medians computed
#' from simulations. For models fit using MCMC (\code{"sampling"}) the posterior
#' sample is used. For optimization (\code{"optimizing"}), the simulations are
#' generated from the asymptotic Gaussian sampling distribution of the
#' parameters. For the \code{"meanfield"} and \code{"fullrank"} variational
#' approximations, draws from the variational approximation to the posterior are
#' used. In all cases, the point estimates reported are the same as the values
#' returned by \code{\link[=coef.stanreg]{coef}}.
#' }
#' \subsection{Uncertainty estimates (MAD_SD)}{
#' The standard deviations reported (labeled \code{MAD_SD} in the print output)
#' are computed from the same set of draws described above and are proportional
#' to the median absolute deviation (\code{\link[stats]{mad}}) from the median.
#' Compared to the raw posterior standard deviation, the MAD_SD will be
#' more robust for long-tailed distributions. These are the same as the values
#' returned by \code{\link[=se.stanreg]{se}}.
#' }
#' \subsection{Additional output}{
#' \itemize{
#' \item For GLMs with group-specific terms (see \code{\link{stan_glmer}}) the printed
#' output also shows point estimates of the standard deviations of the group
#' effects (and correlations if there are both intercept and slopes that vary by
#' group).
#'
#' \item For analysis of variance models (see \code{\link{stan_aov}}) models, an
#' ANOVA-like table is also displayed.
#'
#' \item For joint longitudinal and time-to-event (see \code{\link{stan_jm}}) models
#' the estimates are presented separately for each of the distinct submodels.
#' }
#' }
#'
#' @seealso \code{\link{summary.stanreg}}, \code{\link{stanreg-methods}}
#'
print.stanreg <- function(x, digits = 1, detail = TRUE, ...) {
if (detail) {
cat(x$stan_function)
cat("\n family: ", family_plus_link(x))
cat("\n formula: ", formula_string(formula(x)))
cat("\n observations:", nobs(x))
if (isTRUE(x$stan_function %in%
c("stan_glm", "stan_glm.nb", "stan_lm", "stan_aov"))) {
cat("\n predictors: ", length(coef(x)))
}
if (!is.null(x$call$subset)) {
cat("\n subset: ", deparse(x$call$subset))
}
cat("\n------\n")
}
mer <- is.mer(x)
gamm <- isTRUE(x$stan_function == "stan_gamm4")
ord <- is_polr(x) && !("(Intercept)" %in% rownames(x$stan_summary))
aux_nms <- .aux_name(x)
if (!used.optimizing(x)) {
if (isTRUE(x$stan_function %in% c("stan_lm", "stan_aov"))) {
aux_nms <- c("R2", "log-fit_ratio", aux_nms)
}
mat <- as.matrix(x$stanfit) # don't used as.matrix.stanreg method b/c want access to mean_PPD
nms <- setdiff(rownames(x$stan_summary), c("log-posterior", aux_nms))
if (gamm) {
smooth_sd_nms <- grep("^smooth_sd\\[", nms, value = TRUE)
nms <- setdiff(nms, smooth_sd_nms)
smooth_sd_mat <- mat[, smooth_sd_nms, drop = FALSE]
smooth_sd_estimates <- .median_and_madsd(smooth_sd_mat)
}
if (mer) {
nms <- setdiff(nms, grep("^b\\[", nms, value = TRUE))
}
if (ord) {
cut_nms <- grep("|", nms, fixed = TRUE, value = TRUE)
nms <- setdiff(nms, cut_nms)
cut_mat <- mat[, cut_nms, drop = FALSE]
cut_estimates <- .median_and_madsd(cut_mat)
}
ppd_nms <- grep("^mean_PPD", nms, value = TRUE)
nms <- setdiff(nms, ppd_nms)
coef_mat <- mat[, nms, drop = FALSE]
estimates <- .median_and_madsd(coef_mat)
if (mer) {
estimates <- estimates[!grepl("^Sigma\\[", rownames(estimates)),, drop=FALSE]
}
.printfr(estimates, digits, ...)
if (length(aux_nms)) {
aux_estimates <- .median_and_madsd(mat[, aux_nms, drop=FALSE])
cat("\nAuxiliary parameter(s):\n")
.printfr(aux_estimates, digits, ...)
}
if (ord) {
cat("\nCutpoints:\n")
.printfr(cut_estimates, digits, ...)
}
if (gamm) {
cat("\nSmoothing terms:\n")
.printfr(smooth_sd_estimates, digits, ...)
}
if (mer) {
cat("\nError terms:\n")
print(VarCorr(x), digits = digits + 1, ...)
cat("Num. levels:",
paste(names(ngrps(x)), unname(ngrps(x)), collapse = ", "), "\n")
}
if (is(x, "aov")) {
print_anova_table(x, digits, ...)
}
} else {
# used optimization
nms <- names(x$coefficients)
ppd_nms <- grep("^mean_PPD", rownames(x$stan_summary), value = TRUE)
estimates <- x$stan_summary[nms, 1:2, drop=FALSE]
.printfr(estimates, digits, ...)
if (length(aux_nms)) {
cat("\nAuxiliary parameter(s):\n")
.printfr(x$stan_summary[aux_nms, 1:2, drop=FALSE], digits, ...)
}
}
if (detail) {
cat("\n------\n")
cat("* For help interpreting the printed output see ?print.stanreg\n")
cat("* For info on the priors used see ?prior_summary.stanreg\n")
}
invisible(x)
}
#' @rdname print.stanreg
#' @export
#' @method print stanmvreg
print.stanmvreg <- function(x, digits = 3, ...) {
M <- x$n_markers
mvmer <- is.mvmer(x)
surv <- is.surv(x)
jm <- is.jm(x)
stubs <- paste0("(", get_stub(x), 1:M, "):")
cat(x$stan_function)
if (mvmer) {
for (m in 1:M) {
cat("\n formula", stubs[m], formula_string(formula(x, m = m)))
cat("\n family ", stubs[m], family_plus_link(x, m = m))
}
}
if (surv) {
cat("\n formula (Event):", formula_string(formula(x, m = "Event")))
cat("\n baseline hazard:", x$basehaz$type_name)
}
if (jm) {
sel <- grep("^which", rownames(x$assoc), invert = TRUE, value = TRUE)
assoc <- lapply(1:M, function(m) {
vals <- sel[which(x$assoc[sel,m] == TRUE)]
paste0(vals, " (Long", m, ")")
})
cat("\n assoc: ", paste(unlist(assoc), collapse = ", "))
}
cat("\n------\n")
mat <- as.matrix(x$stanfit)
nms <- collect_nms(rownames(x$stan_summary), M,
stub = get_stub(x), value = TRUE)
# Estimates table for longitudinal submodel(s)
if (mvmer) {
link <- sapply(1:M, function(m) x$family[[m]]$link)
for (m in 1:M) {
terms_m <- terms(x)[[m]]
sel <- attr(terms_m, "response")
yvar <- rownames(attr(terms_m, "factors"))[sel]
if (is.jm(x)) {
cat(paste0("\nLongitudinal submodel", if (M > 1) paste0(" ", m),
": ", yvar,"\n"))
} else {
cat(paste0("\nSubmodel for y", m, ": ", yvar,"\n"))
}
coef_mat <- mat[, c(nms$y[[m]], nms$y_extra[[m]]), drop = FALSE]
# Calculate median and MAD
estimates <- .median_and_madsd(coef_mat)
# Add column with eform
if (link[m] %in% c("log", "logit"))
estimates <- cbind(estimates,
"exp(Median)" = c(exp(estimates[nms$y[[m]], "Median"]),
rep(NA, length(nms$y_extra[[m]]))))
# Print estimates
rownames(estimates) <-
gsub(paste0("^", get_stub(x), m, "\\|"), "", rownames(estimates))
.printfr(estimates, digits, ...)
}
}
# Estimates table for event submodel
if (surv) {
cat("\nEvent submodel:\n")
coef_mat <- mat[, c(nms$e, nms$a, nms$e_extra), drop = FALSE]
# Calculate median and MAD
estimates <- .median_and_madsd(coef_mat)
# Add column with eform
estimates <- cbind(estimates,
"exp(Median)" = c(exp(estimates[c(nms$e, nms$a), "Median"]),
rep(NA, length(nms$e_extra))))
rownames(estimates) <- gsub("^Event\\|", "", rownames(estimates))
rownames(estimates) <- gsub("^Assoc\\|", "", rownames(estimates))
.printfr(estimates, digits, ...)
}
# Estimates table for group-level random effects
if (mvmer) {
cat("\nGroup-level error terms:\n")
print(VarCorr(x), digits = digits + 1, ...)
cat("Num. levels:", paste(names(ngrps(x)), unname(ngrps(x)),
collapse = ", "), "\n")
# Sample average of the PPD
ppd_mat <- mat[, nms$ppd, drop = FALSE]
ppd_estimates <- .median_and_madsd(ppd_mat)
cat("\nSample avg. posterior predictive distribution \nof",
if (is.jm(x)) "longitudinal outcomes:\n" else "y:\n")
.printfr(ppd_estimates, digits, ...)
}
cat("\n------\n")
cat("For info on the priors used see help('prior_summary.stanreg').")
invisible(x)
}
#' Summary method for stanreg objects
#'
#' Summaries of parameter estimates and MCMC convergence diagnostics
#' (Monte Carlo error, effective sample size, Rhat).
#'
#' @export
#' @method summary stanreg
#'
#' @templateVar stanregArg object
#' @template args-stanreg-object
#' @template args-regex-pars
#'
#' @param ... Currently ignored.
#' @param pars An optional character vector specifying a subset of parameters to
#' display. Parameters can be specified by name or several shortcuts can be
#' used. Using \code{pars="beta"} will restrict the displayed parameters to
#' only the regression coefficients (without the intercept). \code{"alpha"}
#' can also be used as a shortcut for \code{"(Intercept)"}. If the model has
#' varying intercepts and/or slopes they can be selected using \code{pars =
#' "varying"}.
#'
#' In addition, for \code{stanmvreg} objects there are some additional shortcuts
#' available. Using \code{pars = "long"} will display the
#' parameter estimates for the longitudinal submodels only (excluding group-specific
#' pparameters, but including auxiliary parameters).
#' Using \code{pars = "event"} will display the
#' parameter estimates for the event submodel only, including any association
#' parameters.
#' Using \code{pars = "assoc"} will display only the
#' association parameters.
#' Using \code{pars = "fixef"} will display all fixed effects, but not
#' the random effects or the auxiliary parameters.
#' \code{pars} and \code{regex_pars} are set to \code{NULL} then all
#' fixed effect regression coefficients are selected, as well as any
#' auxiliary parameters and the log posterior.
#'
#' If \code{pars} is \code{NULL} all parameters are selected for a \code{stanreg}
#' object, while for a \code{stanmvreg} object all
#' fixed effect regression coefficients are selected as well as any
#' auxiliary parameters and the log posterior. See
#' \strong{Examples}.
#' @param probs For models fit using MCMC or one of the variational algorithms,
#' an optional numeric vector of probabilities passed to
#' \code{\link[stats]{quantile}}.
#' @param digits Number of digits to use for formatting numbers when printing.
#' When calling \code{summary}, the value of digits is stored as the
#' \code{"print.digits"} attribute of the returned object.
#'
#' @return The \code{summary} method returns an object of class
#' \code{"summary.stanreg"} (or \code{"summary.stanmvreg"}, inheriting
#' \code{"summary.stanreg"}), which is a matrix of
#' summary statistics and
#' diagnostics, with attributes storing information for use by the
#' \code{print} method. The \code{print} method for \code{summary.stanreg} or
#' \code{summary.stanmvreg} objects is called for its side effect and just returns
#' its input. The \code{as.data.frame} method for \code{summary.stanreg}
#' objects converts the matrix to a data.frame, preserving row and column
#' names but dropping the \code{print}-related attributes.
#'
#' @details
#' \subsection{mean_PPD diagnostic}{
#' Summary statistics are also reported for \code{mean_PPD}, the sample
#' average posterior predictive distribution of the outcome. This is useful as a
#' quick diagnostic. A useful heuristic is to check if \code{mean_PPD} is
#' plausible when compared to \code{mean(y)}. If it is plausible then this does
#' \emph{not} mean that the model is good in general (only that it can reproduce
#' the sample mean), however if \code{mean_PPD} is implausible then it is a sign
#' that something is wrong (severe model misspecification, problems with the
#' data, computational issues, etc.).
#' }
#'
#' @seealso \code{\link{prior_summary}} to extract or print a summary of the
#' priors used for a particular model.
#'
#' @examples
#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") {
#' if (!exists("example_model")) example(example_model)
#' summary(example_model, probs = c(0.1, 0.9))
#'
#' # These produce the same output for this example,
#' # but the second method can be used for any model
#' summary(example_model, pars = c("(Intercept)", "size",
#' paste0("period", 2:4)))
#' summary(example_model, pars = c("alpha", "beta"))
#'
#' # Only show parameters varying by group
#' summary(example_model, pars = "varying")
#' as.data.frame(summary(example_model, pars = "varying"))
#' }
#' @importMethodsFrom rstan summary
summary.stanreg <- function(object,
pars = NULL,
regex_pars = NULL,
probs = c(0.1, 0.5, 0.9),
...,
digits = 1) {
mer <- is.mer(object)
pars <- collect_pars(object, pars, regex_pars)
if (!used.optimizing(object)) {
args <- list(object = object$stanfit, probs = probs)
out <- do.call("summary", args)$summary
if (is.null(pars) && used.variational(object)) {
out <- out[!rownames(out) %in% "log-posterior", , drop = FALSE]
}
if (!is.null(pars)) {
pars <- allow_special_parnames(object, pars)
out <- out[rownames(out) %in% pars, , drop = FALSE]
}
out <- out[!grepl(":_NEW_", rownames(out), fixed = TRUE), , drop = FALSE]
stats <- colnames(out)
if ("n_eff" %in% stats) {
out[, "n_eff"] <- round(out[, "n_eff"])
}
if ("se_mean" %in% stats) {# So people don't confuse se_mean and sd
colnames(out)[stats %in% "se_mean"] <- "mcse"
}
} else { # used optimization
if (!is.null(probs)) {
stanmat <- object$asymptotic_sampling_dist
object$stan_summary <- cbind(Median = apply(stanmat, 2L, median),
MAD_SD = apply(stanmat, 2L, mad),
t(apply(stanmat, 2L, quantile, probs)))
}
object$stan_summary <- cbind(object$stan_summary, object$diagnostics)
if (is.null(pars)) {
famname <- family(object)$family
mark <- names(object$coefficients)
if (is.gaussian(famname))
mark <- c(mark, "sigma")
if (is.nb(famname))
mark <- c(mark, "reciprocal_dispersion")
} else {
mark <- NA
if ("alpha" %in% pars)
mark <- c(mark, "(Intercept)")
if ("beta" %in% pars)
mark <- c(mark, setdiff(names(object$coefficients), "(Intercept)"))
mark <- c(mark, setdiff(pars, c("alpha", "beta")))
mark <- mark[!is.na(mark)]
}
out <- object$stan_summary[mark, , drop=FALSE]
}
structure(
out,
call = object$call,
algorithm = object$algorithm,
stan_function = object$stan_function,
family = family_plus_link(object),
formula = formula(object),
posterior_sample_size = posterior_sample_size(object),
nobs = nobs(object),
npreds = if (isTRUE(object$stan_function %in% c("stan_glm", "stan_glm.nb", "stan_lm")))
length(coef(object)) else NULL,
ngrps = if (mer) ngrps(object) else NULL,
print.digits = digits,
priors = object$prior.info,
no_ppd_diagnostic = no_mean_PPD(object),
class = "summary.stanreg"
)
}
#' @rdname summary.stanreg
#' @export
#' @method print summary.stanreg
#'
#' @param x An object of class \code{"summary.stanreg"}.
print.summary.stanreg <-
function(x, digits = max(1, attr(x, "print.digits")),
...) {
atts <- attributes(x)
cat("\nModel Info:")
cat("\n function: ", atts$stan_function)
cat("\n family: ", atts$family)
cat("\n formula: ", formula_string(atts$formula))
cat("\n algorithm: ", atts$algorithm)
if (!is.null(atts$posterior_sample_size) && atts$algorithm == "sampling") {
cat("\n sample: ", atts$posterior_sample_size,
"(posterior sample size)")
}
cat("\n priors: ", "see help('prior_summary')")
cat("\n observations:", atts$nobs)
if (!is.null(atts$npreds)) {
cat("\n predictors: ", atts$npreds)
}
if (!is.null(atts$call$subset)) {
cat("\n subset: ", deparse(atts$call$subset))
}
if (!is.null(atts$ngrps)) {
cat("\n groups: ", paste0(names(atts$ngrps), " (",
unname(atts$ngrps), ")",
collapse = ", "))
}
cat("\n\nEstimates:\n")
if (used.optimizing(atts) || used.variational(atts)) {
hat <- "khat"
str_diag <- "Monte Carlo diagnostics"
str1 <- "and khat is the Pareto k diagnostic for importance sampling"
str2 <- " (perfomance is usually good when khat < 0.7).\n"
} else {
hat <- "Rhat"
str_diag <- "MCMC diagnostics"
str1 <- "and Rhat is the potential scale reduction factor on split chains"
str2 <- " (at convergence Rhat=1).\n"
}
sel <- which(colnames(x) %in% c("mcse", "n_eff", hat))
has_mc_diagnostic <- length(sel) > 0
if (has_mc_diagnostic) {
xtemp <- x[, -sel, drop = FALSE]
colnames(xtemp) <- paste(" ", colnames(xtemp))
} else {
xtemp <- x
}
ppd_nms <- grep("^mean_PPD", rownames(x), value = TRUE)
has_ppd_diagnostic <- !atts$no_ppd_diagnostic && length(ppd_nms) > 0
if (has_ppd_diagnostic) {
ppd_estimates <- xtemp[rownames(xtemp) %in% ppd_nms, , drop=FALSE]
} else {
ppd_estimates <- NULL
}
xtemp <- xtemp[!rownames(xtemp) %in% c(ppd_nms, "log-posterior"), , drop=FALSE]
# print table of parameter stats
.printfr(xtemp, digits)
if (has_ppd_diagnostic) {
cat("\nFit Diagnostics:\n")
.printfr(ppd_estimates, digits)
cat("\nThe mean_ppd is the sample average posterior predictive ",
"distribution of the outcome variable ",
"(for details see help('summary.stanreg')).\n",
sep = '')
}
if (has_mc_diagnostic) {
cat("\n", str_diag, "\n", sep = '')
mcse_hat <- format(round(x[, c("mcse", hat), drop = FALSE], digits),
nsmall = digits)
n_eff <- format(x[, "n_eff", drop = FALSE], drop0trailing = TRUE)
print(cbind(mcse_hat, n_eff), quote = FALSE)
cat("\nFor each parameter, mcse is Monte Carlo standard error, ",
"n_eff is a crude measure of effective sample size, ",
str1,
str2, sep = '')
}
invisible(x)
}
#' @rdname summary.stanreg
#' @method as.data.frame summary.stanreg
#' @export
as.data.frame.summary.stanreg <- function(x, ...) {
as.data.frame(unclass(x), ...)
}
#' @rdname summary.stanreg
#' @export
#' @method summary stanmvreg
summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL,
probs = NULL, ..., digits = 3) {
pars <- collect_pars(object, pars, regex_pars)
M <- object$n_markers
mvmer <- is.mvmer(object)
surv <- is.surv(object)
jm <- is.jm(object)
if (mvmer) {
# Outcome variable for each longitudinal submodel
y_vars <- sapply(1:M, function(m, object) {
terms_m <- terms(object)[[m]]
sel <- attr(terms_m, "response")
ret <- rownames(attr(terms_m, "factors"))[sel]
}, object = object)
# Family and link for each longitudinal submodel
fam <- lapply(1:M, function(m) family_plus_link(object, m = m))
}
if (jm) {
# Association structure
sel <- grep("^which", rownames(object$assoc), invert = TRUE, value = TRUE)
assoc <- list_nms(lapply(1:M, function(m)
sel[which(object$assoc[sel,m] == TRUE)]), M)
}
# Construct summary table
args <- list(object = object$stanfit)
if (!is.null(probs))
args$probs <- probs
out <- do.call("summary", args)$summary
nms <- collect_nms(rownames(object$stan_summary), M,
stub = get_stub(object), value = TRUE)
if (!is.null(pars)) {
pars2 <- NA
if ("alpha" %in% pars) pars2 <- c(pars2, nms$alpha)
if ("beta" %in% pars) pars2 <- c(pars2, nms$beta)
if ("long" %in% pars) pars2 <- c(pars2, unlist(nms$y), unlist(nms$y_extra))
if ("event" %in% pars) pars2 <- c(pars2, nms$e, nms$a, nms$e_extra)
if ("assoc" %in% pars) pars2 <- c(pars2, nms$a)
if ("fixef" %in% pars) pars2 <- c(pars2, unlist(nms$y), nms$e, nms$a)
if ("b" %in% pars) pars2 <- c(pars2, nms$b)
pars2 <- c(pars2, setdiff(pars,
c("alpha", "beta", "varying", "b",
"long", "event", "assoc", "fixef")))
pars <- pars2[!is.na(pars2)]
} else {
pars <- rownames(object$stan_summary)
pars <- setdiff(pars, b_names(pars, value = TRUE))
if (used.variational(object))
pars <- setdiff(pars, "log-posterior")
}
out <- out[rownames(out) %in% pars, , drop = FALSE]
out <- out[!grepl(":_NEW_", rownames(out), fixed = TRUE), , drop = FALSE]
stats <- colnames(out)
if ("n_eff" %in% stats)
out[, "n_eff"] <- round(out[, "n_eff"])
if ("se_mean" %in% stats) # So people don't confuse se_mean and sd
colnames(out)[stats %in% "se_mean"] <- "mcse"
# Reorder rows of output table
nms_tmp <- rownames(out)
nms_tmp_y <- lapply(1:M, function(m)
grep(paste0("^", get_stub(object), m, "\\|"), nms_tmp, value = TRUE))
nms_tmp_e <- grep("^Event\\|", nms_tmp, value = TRUE)
nms_tmp_a <- grep("^Assoc\\|", nms_tmp, value = TRUE)
nms_tmp_b <- b_names(nms_tmp, value = TRUE)
nms_tmp_Sigma <- grep("^Sigma", nms_tmp, value = TRUE)
nms_tmp_lp <- grep("^log-posterior$", nms_tmp, value = TRUE)
out <- out[c(unlist(nms_tmp_y), nms_tmp_e, nms_tmp_a, nms_tmp_b,
nms_tmp_Sigma, nms_tmp_lp), , drop = FALSE]
# Output object
if (mvmer)
out <- structure(
out, y_vars = y_vars, family = fam, n_markers = object$n_markers,
n_yobs = object$n_yobs, n_grps = object$n_grps)
if (surv)
out <- structure(
out, n_subjects = object$n_subjects, n_events = object$n_events,
basehaz = object$basehaz)
if (jm)
out <- structure(
out, id_var = object$id_var, time_var = object$time_var, assoc = assoc)
structure(
out, formula = object$formula, algorithm = object$algorithm,
stan_function = object$stan_function,
posterior_sample_size = posterior_sample_size(object),
runtime = object$runtime, print.digits = digits,
class = c("summary.stanmvreg", "summary.stanreg"))
}
#' @rdname summary.stanreg
#' @export
#' @method print summary.stanmvreg
print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")),
...) {
atts <- attributes(x)
mvmer <- atts$stan_function %in% c("stan_mvmer", "stan_jm")
jm <- atts$stan_function == "stan_jm"
tab <- if (jm) " " else ""
cat("\nModel Info:\n")
cat("\n function: ", tab, atts$stan_function)
if (mvmer) {
M <- atts$n_markers
stubs <- paste0("(", if (jm) "Long" else "y", 1:M, "):")
for (m in 1:M) {
cat("\n formula", stubs[m], formula_string(atts$formula[[m]]))
cat("\n family ", stubs[m], atts$family[[m]])
}
}
if (jm) {
cat("\n formula (Event):", formula_string(atts$formula[["Event"]]))
cat("\n baseline hazard:", atts$basehaz$type_name)
assoc_fmt <- unlist(lapply(1:M, function(m)
paste0(atts$assoc[[m]], " (Long", m, ")")))
cat("\n assoc: ", paste(assoc_fmt, collapse = ", "))
}
cat("\n algorithm: ", tab, atts$algorithm)
cat("\n priors: ", tab, "see help('prior_summary')")
if (!is.null(atts$posterior_sample_size) && atts$algorithm == "sampling")
cat("\n sample: ", tab, atts$posterior_sample_size, "(posterior sample size)")
if (mvmer) {
obs_vals <- paste0(atts$n_yobs, " (", if (jm) "Long" else "y", 1:M, ")")
cat("\n num obs: ", tab, paste(obs_vals, collapse = ", "))
}
if (jm) {
cat("\n num subjects: ", atts$n_subjects)
cat(paste0("\n num events: ", atts$n_events, " (",
round(100 * atts$n_events/atts$n_subjects, 1), "%)"))
}
if (!is.null(atts$n_grps))
cat("\n groups: ", tab,
paste0(names(atts$n_grps), " (", unname(atts$n_grps), ")", collapse = ", "))
if (atts$algorithm == "sampling") {
maxtime <- max(atts$runtime[, "total"])
if (maxtime == 0) maxtime <- "<0.1"
cat("\n runtime: ", tab, maxtime, "mins")
}
cat("\n\nEstimates:\n")
sel <- which(colnames(x) %in% c("mcse", "n_eff", "Rhat"))
if (!length(sel)) {
.printfr(x, digits)
} else {
xtemp <- x[, -sel, drop = FALSE]
colnames(xtemp) <- paste(" ", colnames(xtemp))
.printfr(xtemp, digits)
cat("\nDiagnostics:\n")
mcse_rhat <- format(round(x[, c("mcse", "Rhat"), drop = FALSE], digits),
nsmall = digits)
n_eff <- format(x[, "n_eff", drop = FALSE], drop0trailing = TRUE)
print(cbind(mcse_rhat, n_eff), quote = FALSE)
cat("\nFor each parameter, mcse is Monte Carlo standard error, ",
"n_eff is a crude measure of effective sample size, ",
"and Rhat is the potential scale reduction factor on split chains",
" (at convergence Rhat=1).\n", sep = '')
}
invisible(x)
}
# internal ----------------------------------------------------------------
.printfr <- function(x, digits, ...) {
print(format(round(x, digits), nsmall = digits), quote = FALSE, ...)
}
.median_and_madsd <- function(x) {
cbind(Median = apply(x, 2, median), MAD_SD = apply(x, 2, mad))
}
# equivalent to isFALSE(object$compute_mean_PPD)
no_mean_PPD <- function(object) {
x <- object$compute_mean_PPD
is.logical(x) && length(x) == 1L && !is.na(x) && !x
}
# Allow "alpha", "beta", "varying" as shortcuts
#
# @param object stanreg object
# @param pars result of calling collect_pars(object, pars, regex_pars)
allow_special_parnames <- function(object, pars) {
pars[pars == "varying"] <- "b"
pars2 <- NA
if ("alpha" %in% pars)
pars2 <- c(pars2, "(Intercept)")
if ("beta" %in% pars) {
beta_nms <- if (is.mer(object))
names(fixef(object)) else names(object$coefficients)
pars2 <- c(pars2, setdiff(beta_nms, "(Intercept)"))
}
if ("b" %in% pars) {
if (is.mer(object)) {
pars2 <- c(pars2, b_names(rownames(object$stan_summary), value = TRUE))
pars[pars == "b"] <- NA
} else {
warning("No group-specific parameters. 'varying' ignored.",
call. = FALSE)
}
}
pars2 <- c(pars2, setdiff(pars, c("alpha", "beta", "varying")))
pars2[!is.na(pars2)]
}
# Family name with link in parenthesis
# @param x stanreg object
# @param ... Optionally include m to specify which submodel for stanmvreg models
family_plus_link <- function(x, ...) {
fam <- family(x, ...)
if (is.character(fam)) {
stopifnot(identical(fam, x$method))
fam <- paste0("ordered [", fam, "]")
} else if (inherits(x, "betareg")) {
fam <- paste0("beta [",
x$family$link,
", link.phi=",
x$family_phi$link,
"]")
} else {
fam <- paste0(fam$family, " [", fam$link, "]")
}
return(fam)
}
# @param formula formula object
formula_string <- function(formula, break_and_indent = TRUE) {
coll <- if (break_and_indent) "--MARK--" else " "
char <- gsub("\\s+", " ", paste(deparse(formula), collapse = coll))
if (!break_and_indent)
return(char)
gsub("--MARK--", "\n\t ", char, fixed = TRUE)
}
# get name of aux parameter based on family
.aux_name <- function(object) {
aux <- character()
if (!is_polr(object)) {
aux <- .rename_aux(family(object))
if (is.na(aux)) {
aux <- character()
}
}
return(aux)
}
# print anova table for stan_aov models
# @param x stanreg object created by stan_aov()
print_anova_table <- function(x, digits, ...) {
labels <- attributes(x$terms)$term.labels
patterns <- gsub(":", ".*:", labels)
dnms <- dimnames(extract(x$stanfit, pars = "beta",
permuted = FALSE))$parameters
groups <- sapply(patterns, simplify = FALSE, FUN = grep, x = dnms)
names(groups) <- gsub(".*", "", names(groups), fixed = TRUE)
groups <- groups[sapply(groups, length) > 0]
effects_dim <- dim(x$effects)
effects <- x$effects^2
effects <- sapply(groups, FUN = function(i) {
apply(effects[, , i, drop = FALSE], 1:2, mean)
})
dim(effects) <- c(effects_dim[-3], ncol(effects))
dim(effects) <- c(nrow(effects) * ncol(effects), dim(effects)[3])
colnames(effects) <- paste("Mean Sq", names(groups))
anova_table <- .median_and_madsd(effects)
cat("\nANOVA-like table:\n")
.printfr(anova_table, digits, ...)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.