bayes factors in model report about report

[mattansb]

• bayesfactor_inclusion
  • text - only results.
  • text_full - Give full explanation about model averaging + posterior probability.
  • table - basically as.data.frame
  • table_full - same as table??
  • docs...
  • tests...
• bayesfactor_models
  • text - method + denominator + {best + worst model}
  • text_full - with BF computation method.
  • table - basically as.data.frame
  • table_full - give also "BF01 compared to best model` column.
  • docs...
  • tests...
• bayesfactor_parameters
  • text - what null (once), BF for each, by name.
  • text_full - same as text??
  • table - basically as.data.frame + row for null and side
  • table_full - same as table??
  • docs...
  • tests...
• bayesfactor_restricted
  • text - ???
  • text_full - ???
  • table - basically as.data.frame
  • table_full - same as table??
  • docs...
  • tests...

Notes to self:

• Should use interpret_bf internally for text, with arg rules.

[mattansb]

@DominiqueMakowski did I miss anything? Any ideas for the "???"?
Any guidelines on how text and text_full / table and table_full should differ?
Is there a timeline on this?
How are BFs dealt with currently in the estimate methods (etc)?

[mattansb]

Just so I know I'm in the right direction, is this what you have in mind? (I'm afraid I haven't been paying close enough attention to report...)

library(bayestestR)
library(report)

mo0 <- lm(Sepal.Length ~ 1, data = iris)
mo1 <- lm(Sepal.Length ~ Species, data = iris)
mo2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
mo3 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)

BFmodels <- bayesfactor_models(mo1, mo2, mo3, denominator = mo0)
inc_bf <- bayesfactor_inclusion(BFmodels, prior_odds = c(1,2,3), match_models = TRUE)

bf_report <- report(inc_bf)

to_table(bf_report) # same for both full = FALSE / TRUE
#> Terms                             | Pr(prior) | Pr(posterior) | Inclusion BF
#> ----------------------------------------------------------------------------
#> Species                           |     0.429 |         0.946 |    3.895e+55
#> Petal.Length                      |     0.286 |         0.946 |    6.891e+26
#> Species:Petal.Length              |     0.429 |         0.054 |        0.038
#>                                   |           |               |             
#> Across matched models only,       |           |               |             
#> with custom prior odds (1, 2, 3). |           |               |

to_text(bf_report, full = FALSE)

We found extreme evidence (BF > 999) in favour of including Species; extreme evidence (BF > 999) in favour of including Petal.Length; strong evidence (BF = 0.04) against including Species:Petal.Length.

to_text(bf_report, full = TRUE)

Bayesian model averaging (BMA) was used to obtain the average evidence for each predictor. Since each model has a prior probability (here we used subjective prior odds of 1, 2, 3), it is possible to sum the prior probability of all models that include a predictor of interest (the prior inclusion probability), and of all models that do not include that predictor (the prior exclusion probability). After the data are observed, we can similarly consider the sums of the posterior models’ probabilities to obtain the posterior inclusion probability and the posterior exclusion probability. The change from prior to posterior inclusion odds is the Inclusion Bayes factor. For each predictor, averaging was done only across models that did not include any interactions with that predictor; additionally, for each interaction predictor, averaging was done only across models that contained the main effect from which the interaction predictor was comprised. This was done to prevent Inclusion Bayes factors from being contaminated with non-relevant evidence (see Mathot, 2017). We found extreme evidence (BF > 999) in favour of including Species; extreme evidence (BF > 999) in favour of including Petal.Length; strong evidence (BF = 0.04) against including Species:Petal.Length.

[DominiqueMakowski]

I was mainly thinking about this case:
report(insight::download_model("stanreg_lm_1"), effsize = "cohen1988", standardize = "refit", test = c("pd", "bf"))

Like implementing the reporting of BFs in the parameters.

Long story short, the table in to_table and to_fulltable (for models) come directly from the parameters table (with one additional column from performance). Then, the text processor picks up whatever indices are present in this table and puts it into text. This is done in the text_parameters.R file. I didn't want to specifically assign this to you it was more like hey we should brainstorm on how to do it best ☺️ so no pressure!

---

I wasn't even thinking about reports for BF inclusion and models yet, but that's awesome!!!

[DominiqueMakowski]

And I believe it will be super useful!

[mattansb]

Wow, then I totally overshot this, hu? 😅

Not sure I follow exactly - what would you like to add to report.stanreg? What is currently missing (to add to the to-do list above)? And will writing report.bayesfactor_parameters help?

I also can't seem to run your example (even after updating all dependencies...):

report::report(insight::download_model("stanreg_lm_1"), effsize = "cohen1988", standardize = "refit", test = c("pd", "bf"))
#> Registered S3 method overwritten by 'xts':
#>   method     from
#>   as.zoo.xts zoo
#> Computation of Bayes factors: sampling priors, please wait...
#> Loading required namespace: logspline
#> Error in UseMethod("format_value"): no applicable method for 'format_value' applied to an object of class "NULL"

^{Created on 2019-08-15 by the reprex package (v0.3.0)}

[DominiqueMakowski]

I also can't seem to run your example (even after updating all dependencies...):

have you updated the latest github's performance?

report::report(insight::download_model("stanreg_lm_1"), test = c("pd", "bf"))
#> Registered S3 method overwritten by 'xts':
#>   method     from
#>   as.zoo.xts zoo
#> Computation of Bayes factors: sampling priors, please wait...
#> Loading required namespace: logspline
#> We fitted a Bayesian linear model to predict mpg with wt. The model's explanatory power is substantial (R2's median = 0.75, adj. R2 = 0.71). The model's intercept is at 37.25. Within this model:
#> 
#>   - The effect of wt has a probability of 100% of being negative and can be considered as large (median = -5.34, 89% CI [-6.23, -4.46], std. median = -0.87).

^{Created on 2019-08-15 by the reprex package (v0.3.0)}

Currently, the text doesn't mention the BFs that are present in the table. I'd like to neatly integrate it with the current textual output... cf. the discussion on the best reporting sentence in the Bayesian indices paper

[mattansb]

Ah I see. I think for BF reporting could be put the in the parentheses, as I mostly see done in the psych papers I read, and reported as:

rep_bf <- function(bf, ...){
  ori_bf <- bf
  dir <- ifelse(log(bf) < 0,"BF01", "BF10")
  
  bf[bf < 1] <- 1 / bf[bf < 1]
  
  paste0(
    parameters::format_bf(bf, name = dir),
    ", considered ",
    report::interpret_bf(ori_bf, include_value = FALSE, ...),
    " an effect."
  )
}

rep_bf(exp(-4:4)) # is vectorized
#> [1] "BF01 = 54.60, considered very strong evidence against an effect."     
#> [2] "BF01 = 20.09, considered strong evidence against an effect."          
#> [3] "BF01 = 7.39, considered moderate evidence against an effect."         
#> [4] "BF01 = 2.72, considered anecdotal evidence against an effect."        
#> [5] "BF10 = 1.00, considered anecdotal evidence in favour of an effect."   
#> [6] "BF10 = 2.72, considered anecdotal evidence in favour of an effect."   
#> [7] "BF10 = 7.39, considered moderate evidence in favour of an effect."    
#> [8] "BF10 = 20.09, considered strong evidence in favour of an effect."     
#> [9] "BF10 = 54.60, considered very strong evidence in favour of an effect."

^{Created on 2019-08-15 by the reprex package (v0.3.0)}

If you agree, feel free to add this in the appropriate place - I don't seem to understand the coding structure of the .report_regressions function/s...

[DominiqueMakowski]

the structure is a bit convoluted indeed (but still simpler than it was ^^). Basically, there are primitives in text_parameters that are called in text_parameters_regression, which is used in report_regression (which is just an internal that can be applied to all regressions to avoid having the same code everywhere). Let me try to provide a start

[DominiqueMakowski]

• pd and ROPE and BF

There is moderate evidence (BF = 0.29) in favour of an absence of effect of X, which has a probability of 90.14% of being negative (Median = -0.03, 89% CI [-0.05, 0.01]), and can be considered as very small (Std. Median = -0.09) and not-significant (6.42% in ROPE).

• pd and ROPE

The effect of X has a probability of 90.14% of being negative (Median = -0.03, 89% CI [-0.05, 0.01]) and can be considered as very small (Std. Median = -0.09) and not-significant (6.42% in ROPE).

• pd and BF

There is moderate evidence (BF = 0.29) in favour of an absence of effect of X, which has a probability of 90.14% of being negative (Median = -0.03, 89% CI [-0.05, 0.01]) and can be considered as very small (Std. Median = -0.09).

• ROPE and BF

There is moderate evidence (BF = 0.29) in favour of an absence of effect of X (Median = -0.03, 89% CI [-0.05, 0.01]) which can be considered as very small (Std. Median = -0.09) and not-significant (6.42% in ROPE).

• Only pd

The effect of X has a probability of 90.14% of being negative (Median = -0.03, 89% CI [-0.05, 0.01]) and can be considered as very small (Std. Median = -0.09).

• Only ROPE

The effect of X (Median = -0.03, 89% CI [-0.05, 0.01]) can be considered as very small (Std. Median = -0.09) and not-significant (6.42% in ROPE).

• Only BF

There is moderate evidence (BF = 0.29) in favour of an absence of effect of X (Median = -0.03, 89% CI [-0.05, 0.01]) which can be considered as very small (Std. Median = -0.09).

[mattansb]

I would change "in favour of an absence of effect of X" simply to "against an effect of X", and have BF01 / BF10 accordingly to also reflect that more explicitly (or just allow values smaller than 1).

[DominiqueMakowski]

agreed

[DominiqueMakowski]

The tricky part is to implement these different cases in a not-horrible way

[mattansb]

Can you point me to the code lines where this is implemented? I'll give it a go (:

…

-- Mattan S. Ben-Shachar, PhD student Department of Psychology & Zlotowski Center for Neuroscience Ben-Gurion University of the Negev The Developmental ERP Lab

On Fri, Aug 16, 2019, 07:27 Dominique Makowski ***@***.***> wrote: The tricky part is to implement these different cases in a not-horrible way — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#38?email_source=notifications&email_token=AINRP6AN7L62Y3YXL4DNC7LQEYUENA5CNFSM4IL2WBB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD4NUAOA#issuecomment-521879608>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AINRP6EQYDH24RLME66MNE3QEYUENANCNFSM4IL2WBBQ> .

[DominiqueMakowski]

Currently, the assemblage is made here:

report/R/text_parameters_regression.R

Lines 46 to 53 in d8212bc

	.text_parameters <- function(parameters, ci, effsize, type = "d", prefix = " - ", bayesian_diagnostic = FALSE) {
	text <- .text_parameters_combine(
	direction = .text_parameters_direction(parameters),
	size = .text_parameters_size(parameters, effsize = effsize, type = type),
	significance = .text_parameters_significance(parameters),
	indices = .text_parameters_indices(parameters, ci = ci),
	bayesian_diagnostic = .text_parameters_bayesian_diagnostic(parameters, bayesian_diagnostic = bayesian_diagnostic)
	)

Which calls this:

report/R/text_parameters.R

Lines 345 to 352 in d8212bc

	.text_parameters_combine <- function(direction = "", size = "", significance = "", indices = "", bayesian_diagnostic = "") {
	text <- direction
	text <- ifelse(significance != "" & size != "", paste0(text, " and can be considered as ", size, " and ", significance),
	ifelse(significance != "" & size == "", paste0(text, " and can be considered as ", significance),
	ifelse(significance == "" & size != "", paste0(text, " and can be considered as ", size), text)
	)
	)

Which get its content, for instance for significance, from here:

report/R/text_parameters.R

Lines 99 to 125 in d8212bc

	.text_parameters_significance <- function(parameters, rope_ci = 1) {
	text <- ""
	if ("p" %in% names(parameters)) {
	text <- interpret_p(parameters$p)
	}
	if ("ROPE_Percentage" %in% names(parameters) & !is.null(rope_ci)) {
	text <- interpret_rope(parameters$ROPE_Percentage, ci = rope_ci)
	}
	text
	}
	#' @keywords internal
	.text_parameters_indices <- function(parameters, ci = 0.89) {
	text <- ""
	if ("Coefficient" %in% names(parameters)) {
	text <- paste0(
	text,
	"beta = ",
	parameters::format_value(parameters$Coefficient)
	)
	}

[mattansb]

Okay... I think I figured it out(?)
I found the best (easiest?) place to put it was inside the parentheses, next to the ROPE - both moved to the end of the perentheses. So it looks like this:

library(report)
mod <- insight::download_model("stanreg_lm_1")
#> Registered S3 method overwritten by 'xts':
#>   method     from
#>   as.zoo.xts zoo
report(mod, test = c("pd", "bf","ROPE"))
#> Computation of Bayes factors: sampling priors, please wait...
#> Loading required namespace: logspline

We fitted a Bayesian linear model to predict mpg with wt. The model's explanatory power is substantial (R2's median = 0.75, adj. R2 = 0.71). The model's intercept is at 37.25. Within this model:
- The effect of wt has a probability of 100% of being negative and can be considered as large and significant (median = -5.34, 89% CI [-6.23, -4.46], 0% in ROPE, std. median = -0.87, BF10 > 999, considered extreme evidence in favour of the effect).

^{Created on 2019-08-16 by the reprex package (v0.3.0)}

What do you think?

[mattansb]

(currently BF is not in .text_parameters_significance... seemed too clumsy to write something like:

The effect of wt has a probability of 100% of being negative and can be considered as large, significant and strongly not supported by the data (median = -5.34, 89% CI [-6.23, -4.46], 0% in ROPE, std. median = -0.87, BF01 > 999, considered extreme against in favour of the effect).

[DominiqueMakowski]

I am not a fan of the interpretation within the parenthesis :)
overall I think it's a good path to try as much as possible to dissociate indices and their interpretation, rather than cramming all interpretation on one side, and all values within brackets, which can lead to confusion about which conclusion is based on what number

let me see if what can I do here

[DominiqueMakowski]

This is not in because it's currently only sexist, but we'll need to re-add that to allow more flexibility