From 46f822fa3239b7c9d0329fd890b72be99d68d052 Mon Sep 17 00:00:00 2001 From: mutlusun Date: Fri, 22 Jan 2021 09:20:02 +0100 Subject: [PATCH] add initial support for brms models - models with multiple dependent variables are not supported yet --- DESCRIPTION | 6 +- NAMESPACE | 12 +++ R/report.brmsfit.R | 109 +++++++++++++++++++++++++++ R/report_random.R | 6 ++ man/report.brmsfit.Rd | 67 ++++++++++++++++ man/report_random.Rd | 6 ++ tests/testthat/test-report.brmsfit.R | 27 +++++++ 7 files changed, 231 insertions(+), 2 deletions(-) create mode 100644 R/report.brmsfit.R create mode 100644 man/report.brmsfit.Rd create mode 100644 tests/testthat/test-report.brmsfit.R diff --git a/DESCRIPTION b/DESCRIPTION index 3359d2c1..cfe5f26a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,6 +37,7 @@ Imports: effectsize (>= 0.4.0), performance (>= 0.5.0) Suggests: + brms, BayesFactor, dplyr, httr, @@ -69,6 +70,9 @@ Collate: 'report.htest.R' 'report.aov.R' 'report.bayesfactor_models.R' + 'report.lme4.R' + 'report.stanreg.R' + 'report.brmsfit.R' 'report.character.R' 'report.compare_performance.R' 'report.data.frame.R' @@ -78,10 +82,8 @@ Collate: 'report.glmmTMB.R' 'report.lavaan.R' 'report.lme.R' - 'report.lme4.R' 'report.numeric.R' 'report.sessionInfo.R' - 'report.stanreg.R' 'report.survreg.R' 'report.test_performance.R' 'report.zeroinfl.R' diff --git a/NAMESPACE b/NAMESPACE index 58eeda5c..0922d74f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -31,6 +31,7 @@ S3method(report,aov) S3method(report,aovlist) S3method(report,bayesfactor_inclusion) S3method(report,bayesfactor_models) +S3method(report,brmsfit) S3method(report,character) S3method(report,compare_performance) S3method(report,data.frame) @@ -55,6 +56,7 @@ S3method(report_effectsize,MixMod) S3method(report_effectsize,anova) S3method(report_effectsize,aov) S3method(report_effectsize,aovlist) +S3method(report_effectsize,brmsfit) S3method(report_effectsize,default) S3method(report_effectsize,glm) S3method(report_effectsize,glmmTMB) @@ -69,6 +71,7 @@ S3method(report_info,MixMod) S3method(report_info,anova) S3method(report_info,aov) S3method(report_info,aovlist) +S3method(report_info,brmsfit) S3method(report_info,default) S3method(report_info,glm) S3method(report_info,glmmTMB) @@ -80,6 +83,7 @@ S3method(report_info,stanreg) S3method(report_info,survreg) S3method(report_info,zeroinfl) S3method(report_intercept,MixMod) +S3method(report_intercept,brmsfit) S3method(report_intercept,default) S3method(report_intercept,glm) S3method(report_intercept,glmmTMB) @@ -93,6 +97,7 @@ S3method(report_model,MixMod) S3method(report_model,anova) S3method(report_model,aov) S3method(report_model,aovlist) +S3method(report_model,brmsfit) S3method(report_model,default) S3method(report_model,glm) S3method(report_model,glmmTMB) @@ -107,6 +112,7 @@ S3method(report_parameters,MixMod) S3method(report_parameters,anova) S3method(report_parameters,aov) S3method(report_parameters,aovlist) +S3method(report_parameters,brmsfit) S3method(report_parameters,character) S3method(report_parameters,compare_performance) S3method(report_parameters,data.frame) @@ -127,6 +133,7 @@ S3method(report_parameters,survreg) S3method(report_parameters,test_performance) S3method(report_parameters,zeroinfl) S3method(report_performance,MixMod) +S3method(report_performance,brmsfit) S3method(report_performance,default) S3method(report_performance,glm) S3method(report_performance,glmmTMB) @@ -137,9 +144,11 @@ S3method(report_performance,merMod) S3method(report_performance,stanreg) S3method(report_performance,survreg) S3method(report_performance,zeroinfl) +S3method(report_priors,brmsfit) S3method(report_priors,default) S3method(report_priors,stanreg) S3method(report_random,MixMod) +S3method(report_random,brmsfit) S3method(report_random,default) S3method(report_random,glmmTMB) S3method(report_random,lme) @@ -149,6 +158,7 @@ S3method(report_statistics,MixMod) S3method(report_statistics,anova) S3method(report_statistics,aov) S3method(report_statistics,aovlist) +S3method(report_statistics,brmsfit) S3method(report_statistics,character) S3method(report_statistics,compare_performance) S3method(report_statistics,data.frame) @@ -172,6 +182,7 @@ S3method(report_table,aov) S3method(report_table,aovlist) S3method(report_table,bayesfactor_inclusion) S3method(report_table,bayesfactor_models) +S3method(report_table,brmsfit) S3method(report_table,character) S3method(report_table,compare_performance) S3method(report_table,data.frame) @@ -198,6 +209,7 @@ S3method(report_text,aov) S3method(report_text,aovlist) S3method(report_text,bayesfactor_inclusion) S3method(report_text,bayesfactor_models) +S3method(report_text,brmsfit) S3method(report_text,character) S3method(report_text,compare_performance) S3method(report_text,data.frame) diff --git a/R/report.brmsfit.R b/R/report.brmsfit.R new file mode 100644 index 00000000..b47bf2e8 --- /dev/null +++ b/R/report.brmsfit.R @@ -0,0 +1,109 @@ +#' Reporting Bayesian Models from brms +#' +#' Create reports for Bayesian models. The description of the parameters +#' follows the Sequential Effect eXistence and sIgnificance Testing framework +#' (see \link[bayestestR:sexit]{SEXIT documentation}). +#' +#' @inheritParams report.lm +#' @inherit report return seealso +#' +#' @examples +#' library(report) +#' +#' # Bayesian models +#' if(require("brms")){ +#' model <- brm(mpg ~ qsec + wt, data = mtcars, refresh=0, iter=300) +#' r <- report(model) +#' r +#' summary(r) +#' as.data.frame(r) +#' summary(as.data.frame(r)) +#' } +#' @include report.lm.R report.stanreg.R report.lme4.R +#' @export +report.brmsfit <- function(x, ...) { + table <- report_table(x, include_effectsize=FALSE, ...) + text <- report_text(x, table = table, ...) + + as.report(text, table = table, ...) +} + +#' @export +report_effectsize.brmsfit <- report_effectsize.lm + +#' @export +report_table.brmsfit <- report_table.lm + +#' @export +report_performance.brmsfit <- report_performance.lm + +#' @export +report_statistics.brmsfit <- report_statistics.lm + +#' @export +report_random.brmsfit <- report_random.merMod + +#' @export +report_model.brmsfit <- report_model.lm + +#' @export +report_text.brmsfit <- report_text.lm + + +# ==================== Specific to Bayes =================================== + + + +# report_priors ----------------------------------------------------------- + +#' @export +report_priors.brmsfit <- function(x, ...) { + params <- bayestestR::describe_prior(x) + params <- params[params$Parameter != "(Intercept)", ] + + # Return empty if no priors info + if (!"Prior_Distribution" %in% names(params) | + nrow(params) == 0 | + all(is.na(params$Prior_Scale))) { + return("") + } + + values <- ifelse(params$Prior_Distribution == "normal", + paste0("mean = ", insight::format_value(params$Prior_Location), ", SD = ", insight::format_value(params$Prior_Scale)), + paste0("location = ", insight::format_value(params$Prior_Location), ", scale = ", insight::format_value(params$Prior_Scale)) + ) + + values <- paste0(params$Prior_Distribution, " (", values, ")") + + if (length(unique(values)) == 1 & nrow(params) > 1) { + text <- paste0("all set as ", values[1]) + } else { + text <- paste0("set as ", format_text(values)) + } + + text <- paste0("Priors over parameters were ", text, " distributions") + as.report_priors(text) +} + + +# report_parameters ------------------------------------------------------- + + +#' @export +report_parameters.brmsfit <- report_parameters.stanreg + + + +# report_intercept -------------------------------------------------------- + + +#' @export +report_intercept.brmsfit <- report_intercept.stanreg + + + +# report_info ------------------------------------------------------------- + + +#' @export +report_info.brmsfit <- report_info.stanreg diff --git a/R/report_random.R b/R/report_random.R index fd236613..427a522b 100644 --- a/R/report_random.R +++ b/R/report_random.R @@ -27,6 +27,12 @@ #' r #' summary(r) #' } +#' if(require("brms")){ +#' model <- brm(mpg ~ disp + (1 | cyl), data = mtcars, refresh=0, iter=1000) +#' r <- report_random(model) +#' r +#' summary(r) +#' } #' @export report_random <- function(x, ...) { UseMethod("report_random") diff --git a/man/report.brmsfit.Rd b/man/report.brmsfit.Rd new file mode 100644 index 00000000..c21842d8 --- /dev/null +++ b/man/report.brmsfit.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/report.brmsfit.R +\name{report.brmsfit} +\alias{report.brmsfit} +\title{Reporting Bayesian Models from brms} +\usage{ +\method{report}{brmsfit}(x, ...) +} +\arguments{ +\item{x}{Object of class \code{lm} or \code{glm}.} + +\item{...}{Arguments passed to or from other methods.} +} +\value{ +A list-object of class \code{report}, which contains further list-objects +with a short and long description of the model summary, as well as a short +and long table of parameters and fit indices. +} +\description{ +Create reports for Bayesian models. The description of the parameters +follows the Sequential Effect eXistence and sIgnificance Testing framework +(see \link[bayestestR:sexit]{SEXIT documentation}). +} +\examples{ +library(report) + +# Bayesian models +if(require("brms")){ + model <- brm(mpg ~ qsec + wt, data = mtcars, refresh=0, iter=300) + r <- report(model) + r + summary(r) + as.data.frame(r) + summary(as.data.frame(r)) +} +} +\seealso{ +Specific components of reports (especially for stats models): +\itemize{ + \item \code{\link{report_table}} + \item \code{\link{report_parameters}} + \item \code{\link{report_statistics}} + \item \code{\link{report_effectsize}} + \item \code{\link{report_model}} + \item \code{\link{report_priors}} + \item \code{\link{report_random}} + \item \code{\link{report_performance}} + \item \code{\link{report_info}} + \item \code{\link{report_text}} +} +Other types of reports: +\itemize{ + \item \code{\link{report_system}} + \item \code{\link{report_packages}} + \item \code{\link{report_participants}} + \item \code{\link{report_sample}} + \item \code{\link{report_date}} +} +Methods: +\itemize{ + \item \code{\link{as.report}} +} +Template file for supporting new models: +\itemize{ + \item \code{\link{report.default}} +} +} diff --git a/man/report_random.Rd b/man/report_random.Rd index 2e045ab3..89783c02 100644 --- a/man/report_random.Rd +++ b/man/report_random.Rd @@ -35,4 +35,10 @@ if(require("rstanarm")){ r summary(r) } +if(require("brms")){ + model <- brm(mpg ~ disp + (1 | cyl), data = mtcars, refresh=0, iter=1000) + r <- report_random(model) + r + summary(r) +} } diff --git a/tests/testthat/test-report.brmsfit.R b/tests/testthat/test-report.brmsfit.R new file mode 100644 index 00000000..62a35e4a --- /dev/null +++ b/tests/testthat/test-report.brmsfit.R @@ -0,0 +1,27 @@ +if (require("testthat") && require("brms")) { + model <- brm(mpg ~ qsec + wt, data = mtcars, refresh=0, iter=300) + + testthat::test_that("model", { + r <- report(model) + testthat::expect_is(summary(r), "character") + testthat::expect_is(as.data.frame(r), "data.frame") + + expect_equal( + as.data.frame(r)$Parameter, + c( + "(Intercept)", "qsec", "wt", NA, "ELPD", "LOOIC", "WAIC", "R2", + "R2 (adj.)", "Sigma" + ) + ) + expect_equal( + as.data.frame(r)$Median, + c(19.906865, 0.930295, -5.119548, rep(NA, 7)), + tolerance = 1e-1 + ) + expect_equal( + as.data.frame(r)$pd, + c(rep(1, 3), rep(NA, 7)), + tolerance = 1e-1 + ) + }) +}