Nothing
#' Convert (refit) a Bayesian model to frequentist
#'
#' Refit Bayesian model as frequentist. Can be useful for comparisons.
#'
#' @param model A Bayesian model.
#' @param data Data used by the model. If `NULL`, will try to extract it
#' from the model.
#' @param REML For mixed effects, should models be estimated using
#' restricted maximum likelihood (REML) (`TRUE`, default) or maximum
#' likelihood (`FALSE`)?
#' @examplesIf require("rstanarm")
#' \donttest{
#' # Rstanarm ----------------------
#' # Simple regressions
#' model <- rstanarm::stan_glm(Sepal.Length ~ Species,
#' data = iris, chains = 2, refresh = 0
#' )
#' bayesian_as_frequentist(model)
#'
#' model <- rstanarm::stan_glm(vs ~ mpg,
#' family = "binomial",
#' data = mtcars, chains = 2, refresh = 0
#' )
#' bayesian_as_frequentist(model)
#'
#' # Mixed models
#' model <- rstanarm::stan_glmer(
#' Sepal.Length ~ Petal.Length + (1 | Species),
#' data = iris, chains = 2, refresh = 0
#' )
#' bayesian_as_frequentist(model)
#'
#' model <- rstanarm::stan_glmer(vs ~ mpg + (1 | cyl),
#' family = "binomial",
#' data = mtcars, chains = 2, refresh = 0
#' )
#' bayesian_as_frequentist(model)
#' }
#'
#' @export
convert_bayesian_as_frequentist <- function(model, data = NULL, REML = TRUE) {
if (is.null(data)) {
data <- insight::get_data(model)
}
info <- insight::model_info(model, verbose = FALSE)
model_formula <- insight::find_formula(model)
model_family <- insight::get_family(model)
# fix exception: The 0 + Intercept syntax in brms can be used to facilitate
# prior specification for the intercept, but but it leads to issues where it
# wrongly can be believed that Intercept is a variable and not a special term.
f_string <- insight::safe_deparse(model_formula$conditional)
if (grepl("0 + Intercept", f_string, fixed = TRUE)) {
model_formula$conditional <- stats::as.formula(gsub("0 + Intercept", "1", f_string, fixed = TRUE))
}
if (inherits(model_family, "brmsfamily")) {
insight::check_if_installed("glmmTMB")
# exception: ordbetareg()
if ("custom" %in% model_family$family && all(model_family$name == "ord_beta_reg")) {
model_family <- glmmTMB::ordbeta()
} else {
# not all families return proper objects from "get", so we capture
# some families via switch here...
model_family <- .safe(switch(model_family$family,
beta = glmmTMB::beta_family(link = model_family$link),
beta_binomial = glmmTMB::betabinomial(link = model_family$link),
negbinomial = glmmTMB::nbinom1(link = model_family$link),
lognormal = glmmTMB::lognormal(link = model_family$link),
student = glmmTMB::t_family(link = model_family$link),
get(model_family$family)(link = model_family$link)
))
}
}
# if family could not be identified, stop here
if (is.null(model_family)) {
insight::format_error("Model could not be automatically converted to frequentist model.")
}
# first attempt
freq <- tryCatch(.convert_bayesian_as_frequentist(
info = info, formula = model_formula, data = data, family = model_family, REML = REML
), error = function(e) e)
if (inherits(freq, "error")) {
# try again to extract family, using generic approach
model_family <- get(model_family$family)(link = model_family$link)
freq <- .convert_bayesian_as_frequentist(
info = info, formula = model_formula, data = data, family = model_family, REML = REML
)
}
if (inherits(freq, "error")) {
insight::format_error("Model could not be automatically converted to frequentist model.")
}
freq
}
# internal
.convert_bayesian_as_frequentist <- function(info, formula, data, family, REML = TRUE) {
# TODO: Check for
# nonlinear formulas,
# correlation structures,
# weights,
# offset,
# subset,
# knots,
# meta-analysis
if (info$is_dispersion || info$is_orderedbeta || info$is_beta || info$is_betabinomial || info$is_zero_inflated || info$is_zeroinf || info$is_hurdle || info$is_negbin) { # nolint
insight::check_if_installed("glmmTMB")
cond_formula <- .rebuild_cond_formula(formula)
dispformula <- formula$dispersion
if (is.null(dispformula)) dispformula <- formula$sigma
if (is.null(dispformula)) dispformula <- ~1
ziformula <- formula$zero_inflated
if (is.null(ziformula)) ziformula <- formula$zi
if (is.null(ziformula)) ziformula <- ~0
freq <- tryCatch(
glmmTMB::glmmTMB(
formula = cond_formula,
ziformula = ziformula,
dispformula = dispformula,
family = family,
data = data,
REML = REML
),
error = function(e) e
)
} else if (info$is_gam) {
insight::check_if_installed("gamm4")
freq <- tryCatch(
gamm4::gamm4(
formula = formula$conditional,
random = formula$random,
family = family,
data = data
),
error = function(e) e
)
} else if (info$is_mixed) {
insight::check_if_installed("lme4")
insight::check_if_installed("glmmTMB")
cond_formula <- .rebuild_cond_formula(formula)
if (info$is_linear) {
freq <- tryCatch(
lme4::lmer(
formula = cond_formula,
data = data
),
error = function(e) e
)
} else {
## TODO: check if beta/Gamma are correctly captured
freq <- tryCatch(
lme4::glmer(
formula = cond_formula,
family = family,
data = data
),
error = function(e) e
)
if (inherits(freq, "error")) {
freq <- tryCatch(
glmmTMB::glmmTMB(
formula = cond_formula,
family = family,
data = data
),
error = function(e) e
)
}
}
} else if (info$is_linear) {
freq <- stats::lm(formula$conditional, data = data)
} else {
freq <- stats::glm(formula$conditional, data = data, family = family)
}
freq
}
.rebuild_cond_formula <- function(formula) {
if (is.null(formula$random)) {
return(formula$conditional)
}
if (is.list(formula$random)) {
random_formula <- paste(
lapply(
formula$random, function(x) {
paste0("(", as.character(x)[-1], ")")
}
),
collapse = " + "
)
} else {
random_formula <- paste0("(", as.character(formula$random)[-1], ")")
}
fixed_formula <- paste(as.character(formula$conditional)[c(2, 1, 3)], collapse = " ")
stats::as.formula(paste(fixed_formula, random_formula, sep = " + "))
}
#' @rdname convert_bayesian_as_frequentist
#' @export
bayesian_as_frequentist <- convert_bayesian_as_frequentist
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.