---
title: "Effect sizes for Bayesian models"
output: 
  github_document:
    toc: true
    fig_width: 10.08
    fig_height: 6
  rmarkdown::html_vignette:
    toc: true
    fig_width: 10.08
    fig_height: 6
tags: [r, effect size, rules of thumb, bayesian, effect size]
vignette: >
  \usepackage[utf8]{inputenc}
  %\VignetteIndexEntry{Effect sizes for Bayesian models}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
bibliography: bibliography.bib
---


```{r message=FALSE, warning=FALSE, include=FALSE}
library(knitr)
options(knitr.kable.NA = '')
options(digits = 2)
knitr::opts_chunk$set(comment = ">")

if (!requireNamespace("parameters", quietly = TRUE) ||
    !requireNamespace("rstanarm", quietly = TRUE) ||
    !requireNamespace("bayestestR", quietly = TRUE) ||
    !requireNamespace("ppcor", quietly = TRUE)) {
  knitr::opts_chunk$set(eval = FALSE)
}

library(effectsize)
```


**Obtaining indices of effect size for Bayesian model is currently an issue**, as no test statistics are present to help us compute such indices. The two alternatives (both available in **effectsize**) is 1) to compute [**standardized parameters**](https://easystats.github.io/effectsize/articles/standardize_parameters.html) or 2) obtain effect sizes by **test-statistic approximation**. This second method is described here.


## Bayesian Regression Coefficients to Correlation via Test-Statistic Approximation

### Overview
```{r, warning=FALSE, message=FALSE}
library("ppcor")
df <- iris[, 1:4]  # Remove the Species factor
ppcor::pcor(df)$estimate[2:4, 1]  # Select the rows of interest
```

The goal is to retrieve coefficients similar to the above **(partial) correlations** for this multiple regression model: `Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width`. This can easily be achieved for frequentist models by converting the *t* statistic into a correlation:


```{r, warning=FALSE, message=FALSE}
library(effectsize)
library(parameters)

model <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = df) 

parameters <- model_parameters(model)[2:4,]
convert_t_to_r(parameters$t, parameters$df_residual)
```

Note that these are *not* equivalent to *standardized parameters* from multiple regressions (from which coefficients can be higher than 1).


```{r, warning=FALSE, message=FALSE}
standardize_parameters(model)$Std_Coefficient[2:4]
```

Let's start by fitting the Bayesian regression:

```{r, warning=FALSE, message=FALSE, eval = FALSE}
library("rstanarm")
model <- stan_glm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = df) 
```
```{r, warning=FALSE, message=FALSE, echo = FALSE}
library("rstanarm")
model <- stan_glm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = df, refresh = 0, chains = 2) 
```


The **effectsize** package provides the `posteriors_to_r()` function, which attempts to estimate the posterior of the (partial) correlation coefficient by approximating the *t* statistic (deviding the coefficient by the SD of the posterior) and using frequentist degrees of freedom. It is an hybrid method that needs to be validated. However, the results appear as very close:

```{r, warning=FALSE, message=FALSE, eval=FALSE}
library(bayestestR)

r <- convert_posteriors_to_r(model)
bayestestR::describe_posterior(r)$Median[2:4]
```

### Applicability

Does it work in all cases?

#### Logistic Models

Let's start with the frequentist model:

```{r, warning=FALSE, message=FALSE, eval=FALSE}
model <- glm(vs ~ cyl + disp + drat, data = mtcars, family = "binomial")

parameters <- model_parameters(model)
parameters$r <- convert_z_to_r(parameters$z, n = insight::n_obs(model))
parameters
```

However, as logistic models return log-odds, these can be directly converted to *r*:

```{r, warning=FALSE, message=FALSE, eval=FALSE}
parameters$r_from_odds <- convert_odds_to_r(parameters$Coefficient, log = TRUE)
parameters
```

And now the Bayesian:

```{r, warning=FALSE, message=FALSE, eval = FALSE}
model <- stan_glm(vs ~ cyl + disp + drat,, data = mtcars, family = "binomial") 
parameters <- model_parameters(model)
r <- convert_posteriors_to_r(model)
parameters$r <- bayestestR::describe_posterior(r)$Median
parameters$r_from_odds <- convert_odds_to_r(parameters$Median, log = TRUE)
parameters
```
```{r, warning=FALSE, message=FALSE, echo = FALSE, eval=FALSE}
model <- stan_glm(vs ~ cyl + disp + drat,, data = mtcars, family = "binomial", refresh = 0) 
parameters <- model_parameters(model)
r <- convert_posteriors_to_r(model)
parameters$r <- bayestestR::describe_posterior(r)$Median
parameters$r_from_odds <- convert_odds_to_r(parameters$Median, log = TRUE)
parameters
```


## References