tidy_model <- function(
model, ci.lvl, tf, type, bpe, robust, facets, show.zeroinf, p.val = NULL,
standardize = FALSE, bootstrap = FALSE, iterations = 1000, seed = NULL,
p_adjust = NULL, keep = NULL, drop = NULL, std.response = TRUE, ...) {
if (!is.logical(standardize) && standardize == "") standardize <- NULL
if (is.logical(standardize) && standardize == FALSE) standardize <- NULL
if (is.stan(model)) {
out <- tidy_stan_model(model, ci.lvl, tf, type, bpe, show.zeroinf, facets, ...)
} else {
if (!is.null(standardize)) {
if (isTRUE(standardize)) standardize <- "std"
model <- datawizard::standardize(
model,
two_sd = isTRUE(standardize == "std2"),
include_response = isTRUE(std.response)
)
}
if (!is.null(seed)) {
set.seed(seed)
}
minfo <- insight::model_info(model)
if (!is.null(minfo) && ((show.zeroinf && minfo$is_zero_inflated) || minfo$is_dispersion)) {
component <- "all"
} else {
component <- "conditional"
}
if (is.null(p.val)) {
if (inherits(model, c("glm", "polr")) && !inherits(model, "svyglm")) {
p.val <- "profile"
} else {
p.val <- "wald"
}
}
ci_method <- switch(
p.val,
"r" = ,
"residual" = "residual",
"wald" = "wald",
"kr" = ,
"kenward" = "kenward",
"s" = ,
"satterthwaite" = "satterthwaite",
"n" = ,
"normal" = "normal",
"profile" = "profile",
p.val
)
if (isTRUE(bootstrap)) {
ci_method <- "eti"
}
if (!insight::is_empty_object(insight::compact_list(robust))) {
if (!is.null(robust$vcov.type)) {
robust$vcov.args[["type"]] <- robust$vcov.type
}
model_params <- parameters::model_parameters(model, ci = ci.lvl, component = component, bootstrap = bootstrap, iterations = iterations, vcov = robust$vcov.fun, vcov_args = robust$vcov.args, ci_method = ci_method, p_adjust = p_adjust, effects = "fixed", keep = keep, drop = drop, verbose = FALSE)
} else {
model_params <- parameters::model_parameters(model, ci = ci.lvl, component = component, bootstrap = bootstrap, iterations = iterations, ci_method = ci_method, p_adjust = p_adjust, effects = "fixed", keep = keep, drop = drop, verbose = FALSE)
}
out <- insight::standardize_names(model_params, style = "broom")
# warning for p-values?
tryCatch({
if (insight::model_info(model)$is_mixed && ci_method == "kenward" && insight::find_algorithm(model)$algorithm != "REML") {
warning("Model was not fitted by REML. Re-fitting model using REML, but p-values, df, etc. still might be unreliable.", call. = FALSE)
}
},
error = function(e) { NULL }
)
column <- which(colnames(out) == "response")
if (length(column)) colnames(out)[column] <- ifelse(isTRUE(facets), "facet", "response.level")
column <- which(colnames(out) == "component")
if (length(column)) colnames(out)[column] <- "wrap.facet"
if (!is.null(out$effect) && "random" %in% out$effect) {
out$group <- NULL
}
column <- which(colnames(out) == "group")
if (length(column)) colnames(out)[column] <- "wrap.facet"
# remove duplicated column names
dupl_cols <- duplicated(colnames(out))
if (any(dupl_cols)) {
out <- out[!dupl_cols]
}
if ("component" %in% colnames(out)) {
out$component[out$component == "zero_inflated"] <- "Zero-Inflated Model"
out$component[out$component == "zi"] <- "Zero-Inflated Model"
out$component[out$component == "conditional"] <- "Conditional Model"
out$component[out$component == "count"] <- "Conditional Model"
}
attr(out, "pretty_names") <- attributes(model_params)$pretty_names
}
out
}
#' @importFrom rlang .data
tidy_stan_model <- function(model, ci.lvl, tf, type, bpe, show.zeroinf, facets, ...) {
# set defaults
p.inner <- .5
p.outer <- ci.lvl
# get model information
modfam <- insight::model_info(model)
if (insight::is_multivariate(model))
modfam <- modfam[[1]]
# additional arguments for 'effects()'-function?
add.args <- lapply(match.call(expand.dots = F)$`...`, function(x) x)
# check whether we have "prob.inner" and "prob.outer" argument
# and if so, use these for CI and Bayesian point estimate
if ("prob.inner" %in% names(add.args)) p.inner <- eval(add.args[["prob.inner"]])
if ("prob.outer" %in% names(add.args)) p.outer <- eval(add.args[["prob.outer"]])
# get two CI-intervals
if (type == "re")
ty <- "random"
else
ty <- "fixed"
d1 <- bayestestR::ci(model, ci = p.outer, effects = ty, component = "all")
d2 <- bayestestR::ci(model, ci = p.inner, effects = ty, component = "all")
if (!is.null(tf)) {
funtrans <- match.fun(tf)
d1$CI_low <- funtrans(d1$CI_low)
d1$CI_high <- funtrans(d1$CI_high)
d2$CI_low <- funtrans(d2$CI_low)
d2$CI_high <- funtrans(d2$CI_high)
}
# bind columns, so we have inner and outer hdi interval
dat <- d2 %>%
dplyr::select(.data$CI_low, .data$CI_high) %>%
sjmisc::var_rename(CI_low = "conf.low50", CI_high = "conf.high50") %>%
sjmisc::add_columns(d1) %>%
sjmisc::var_rename(CI_low = "conf.low", CI_high = "conf.high", Parameter = "term") %>%
dplyr::select(-.data$CI, -.data$Effects, -.data$Component)
# for brmsfit models, we need to remove some columns here to
# match data rows later
mod.dat <- as.data.frame(model, optional = FALSE)
if (inherits(model, "brmsfit")) {
re.sd <- string_starts_with("sd_", x = colnames(mod.dat))
re.cor <- string_starts_with("cor_", x = colnames(mod.dat))
lp <- string_starts_with("lp__", x = colnames(mod.dat))
resp.cor <- string_starts_with("rescor__", x = colnames(mod.dat))
priors <- string_starts_with("prior_", x = colnames(mod.dat))
xme <- string_starts_with(pattern = "Xme_me", x = colnames(mod.dat))
xme.sd <- string_starts_with(pattern = "sdme_me", x = colnames(mod.dat))
brmsfit.removers <- unique(c(re.sd, re.cor, lp, resp.cor, priors, xme, xme.sd))
if (!sjmisc::is_empty(brmsfit.removers))
mod.dat <- dplyr::select(mod.dat, !! -brmsfit.removers)
# also clean prepared data frame
resp.cor <- string_starts_with("rescor__", x = dat$term)
if (!sjmisc::is_empty(resp.cor))
dat <- dplyr::slice(dat, !! -resp.cor)
}
# do transformation on posterior samples first,
# then summarize (see https://discourse.mc-stan.org/t/monotonic-effects-in-non-gaussian-models/6353/5)
# need to transform point estimate as well
if (!is.null(tf)) {
funtrans <- match.fun(tf)
all.cols <- sjmisc::seq_col(mod.dat)
simp.pars <- string_starts_with("simo_mo", colnames(mod.dat))
if (!sjmisc::is_empty(simp.pars)) all.cols <- all.cols[-simp.pars]
for (i in all.cols) mod.dat[[i]] <- funtrans(mod.dat[[i]])
}
# add bayesian point estimate
est <- purrr::map_dbl(mod.dat, ~ sjmisc::typical_value(.x, fun = bpe))
dat <- data_frame(
term = names(est),
estimate = est,
p.value = 0,
std.error = purrr::map_dbl(mod.dat, stats::mad)
) %>%
dplyr::inner_join(
dat,
by = "term"
)
# sort columns, for tab_model()
sorted_columns <- intersect(
c("term", "estimate", "std.error", "conf.low", "conf.high", "conf.low50", "conf.high50", "p.value"),
colnames(dat)
)
dat <- dat[, sorted_columns]
# remove some of the information not needed for plotting
if ("sigma" %in% dat$term) dat <- dplyr::filter(dat, .data$term != "sigma")
if ("lp__" %in% dat$term) dat <- dplyr::filter(dat, .data$term != "lp__")
if ("shape" %in% dat$term) dat <- dplyr::filter(dat, .data$term != "shape")
# remove sd_c and cor_ row
re <- string_starts_with("sd_", x = dat$term)
if (!sjmisc::is_empty(re)) dat <- dplyr::slice(dat, !! -re)
re <- string_starts_with("cor_", x = dat$term)
if (!sjmisc::is_empty(re)) dat <- dplyr::slice(dat, !! -re)
# check if we need to keep or remove random effects
re <- string_starts_with("b[", x = dat$term)
re.s <- string_starts_with("Sigma[", x = dat$term)
re.i <- intersect(
string_starts_with("r_", x = dat$term),
string_ends_with(".", x = dat$term)
)
# and all random effect error terms
if (!sjmisc::is_empty(re.s)) dat <- dplyr::slice(dat, !! -re.s)
if (type == "est") {
# remove all random effect intercepts
if (!sjmisc::is_empty(re)) dat <- dplyr::slice(dat, !! -re)
# remove random effects from brmsfit-models
if (!sjmisc::is_empty(re.i)) dat <- dplyr::slice(dat, !! -re.i)
} else if (type == "re") {
# remove all random effect intercepts
if (!sjmisc::is_empty(re)) dat <- dplyr::slice(dat, !! re)
# remove random effects from brmsfit-models
if (!sjmisc::is_empty(re.i)) dat <- dplyr::slice(dat, !! re.i)
}
# for plot-type random effects, make sure that the random effects
# are plotted as facet grid, grouped by groups
if (type == "re") {
dat$facet <- "re"
# find random intercepts
ri <- grep("b\\[\\(Intercept\\) (.*)\\]", dat$term)
if (!sjmisc::is_empty(ri)) {
dat$facet[ri] <- "(Intercept)"
dat$term[ri] <- gsub("b\\[\\(Intercept\\) (.*)\\]", "\\1", dat$term[ri])
}
# find random intercepts
ri1 <- grep("r_(.*)\\.(.*)\\.", dat$term)
ri2 <- which(gsub("r_(.*)\\.(.*)\\.", "\\2", dat$term) == "Intercept")
if (!sjmisc::is_empty(ri1)) {
ri <- intersect(ri1, ri2)
dat$facet[ri] <- "(Intercept)"
dat$term[ri] <- gsub("r_(.*)\\.(.*)\\.", "\\1", dat$term[ri])
}
# fix multiple random intercepts
if (inherits(model, "brmsfit")) {
pattern <- "(.*)\\.(.*)"
} else {
pattern <- "(.*)\\:(.*)"
}
interc <- which(dat$facet == "(Intercept)")
if (!sjmisc::is_empty(interc)) {
interc.grps <- gsub(pattern, "\\1", dat$term[interc])
resp.lvl <- gsub(pattern, "\\2", dat$term[interc])
if (!sjmisc::is_empty(interc.grps) && dplyr::n_distinct(interc.grps) > 1) {
dat$facet[interc] <- sprintf("(Intercept: %s)", interc.grps)
dat$term[interc] <- resp.lvl
}
}
# find random slopes
rs1 <- grep("b\\[(.*) (.*)\\]", dat$term)
rs2 <- which(gsub("b\\[(.*) (.*)\\]", "\\1", dat$term) != "(Intercept)")
if (!sjmisc::is_empty(rs1)) {
rs <- intersect(rs1, rs2)
rs.string <- gsub("b\\[(.*) (.*)\\]", "\\1", dat$term[rs])
dat$facet[rs] <- rs.string
dat$term[rs] <- gsub("b\\[(.*) (.*)\\]", "\\2", dat$term[rs])
}
# find random slopes
rs1 <- grep("r_(.*)\\.(.*)\\.", dat$term)
rs2 <- which(gsub("r_(.*)\\.(.*)\\.", "\\2", dat$term) != "Intercept")
if (!sjmisc::is_empty(rs1)) {
rs <- intersect(rs1, rs2)
rs.string <- gsub("r_(.*)\\.(.*)\\.", "\\2", dat$term[rs])
dat$facet[rs] <- rs.string
dat$term[rs] <- gsub("r_(.*)\\.(.*)\\.", "\\1", dat$term[rs])
}
}
# categorical model?
if (inherits(model, "brmsfit") && modfam$is_categorical) {
# terms of categorical models are prefixed with "mu"
if (length(string_starts_with("b_mu", x = dat$term)) == nrow(dat)) {
dat$term <- substr(dat$term, 5, max(nchar(dat$term)))
# create "response-level" variable
dat <- sjmisc::add_variables(dat, response.level = "", .before = 1)
dat$response.level <- gsub("(.*)\\_(.*)", "\\1", dat$term)
dat$term <- gsub("(.*)\\_(.*)", "\\2", dat$term)
}
}
# multivariate-response model?
if (inherits(model, "brmsfit") && insight::is_multivariate(model)) {
# get response variables
responses <- stats::formula(model)$responses
# also clean prepared data frame
resp.sigma1 <- string_starts_with("sigma_", x = dat$term)
resp.sigma2 <- string_starts_with("b_sigma_", x = dat$term)
resp.sigma <- c(resp.sigma1, resp.sigma2)
if (!sjmisc::is_empty(resp.sigma))
dat <- dplyr::slice(dat, !! -resp.sigma)
# create "response-level" variable
dat <- sjmisc::add_variables(dat, response.level = "", .before = 1)
# copy name of response into new character variable
# and remove response name from term name
for (i in responses) {
m <- grep(pattern = sprintf("_%s_", i), x = dat$term)
dat$response.level[intersect(which(dat$response.level == ""), m)] <- i
dat$term <- gsub(sprintf("b_%s_", i), "", dat$term, fixed = TRUE)
dat$term <- gsub(sprintf("s_%s_", i), "", dat$term, fixed = TRUE)
}
# check whether each category should be printed in facets, or
# in a single graph (with dodged geoms)
if (!missing(facets) && isTRUE(facets))
colnames(dat)[1] <- "facet"
else
colnames(dat)[1] <- "response.level"
}
# do we have a zero-inflation model?
if (modfam$is_zero_inflated || sjmisc::str_contains(dat$term, "b_zi_", ignore.case = T)) {
dat$wrap.facet <- "Conditional Model"
# zero-inflated part
zi <- string_starts_with("b_zi_", x = dat$term)
# check if zero-inflated part should be shown or removed
if (show.zeroinf) {
dat$wrap.facet[zi] <- "Zero-Inflated Model"
dat$term[zi] <- sub(pattern = "b_zi_", replacement = "b_", x = dat$term[zi], fixed = T)
} else {
if (!sjmisc::is_empty(zi)) dat <- dplyr::slice(dat, !! -zi)
}
}
# check model for monotonic effects
simplex.terms <- string_starts_with(pattern = "simo_mo", x = dat$term)
if (!sjmisc::is_empty(simplex.terms)) {
if (!obj_has_name(dat, "wrap.facet")) {
dat$wrap.facet <- ""
dat$wrap.facet[simplex.terms] <- "Simplex Parameters"
} else {
dat$wrap.facet[simplex.terms] <- sprintf(
"%s (Simplex Parameters)",
dat$wrap.facet[simplex.terms]
)
}
}
# remove facet column if not necessary
if (!show.zeroinf && obj_has_name(dat, "wrap.facet"))
dat <- dplyr::select(dat, -.data$wrap.facet)
dat
}
.get_confint <- function(ci.lvl = .95) {
if (!is.null(ci.lvl) && !is.na(ci.lvl))
(1 + ci.lvl) / 2
else
.975
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.