https://github.com/cran/bayestestR
Raw File
Tip revision: fe07bfa906d7e155439160caee538a3449cd3877 authored by Dominique Makowski on 08 April 2019, 08:42:41 UTC
version 0.1.0
Tip revision: fe07bfa
indicesEstimationComparison.Rmd
---
title: "In-Depth 1: Comparison of Point-Estimates"
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, bayesian, posterior, test]
vignette: >
  \usepackage[utf8]{inputenc}
  %\VignetteIndexEntry{In-Depth 1: Comparison of Point-Estimates}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

This vignette can be referred to by citing the package:

- Makowski, D. \& Lüdecke, D. (2019). *Understand and Describe Bayesian Models and Posterior Distributions using BayestestR*. Available from https://github.com/easystats/bayestestR. DOI: [10.5281/zenodo.2556486](https://zenodo.org/record/2556486).



# Effect Point-Estimates in the Bayesian Framework


## Introduction 

One of the main difference between the Bayesian and the frequentist frameworks is that the former returns a probability distribution of each effect (*i.e.*, parameter of interest of a model, such as a regression slope) instead of a single value. However, there is still a need and demand, for reporting or use in further analysis, for a single value (**point-estimate**) that best characterise the underlying posterior distribution.

There are three main indices used in the literature for effect estimation: the **mean**, the **median** or the **MAP** (Maximum A Posteriori) estimate (roughly corresponding to the mode (the "peak") of the distribution). Unfortunately, there is no consensus about which one to use, as no systematic comparison has ever been done.

In the present work, we will compare these three point-estimates of effect between themselves, as well as with the widely known **beta**, extracted from a comparable frequentist model. With this comparison, we expect to draw bridges and relationships between the two frameworks, helping and easing the transition for the public.

## Methods

### Simulate Regression Data with Noise

We simulated 36000 Bayesian and frequentist linear regression models with systematic variation of "true" effect (1 or 0), sample size (20, 40, 60), levels of noise (0.1, 2.5, 5), and priors (correct - same as true effect or incorrect).

```{r message=FALSE, warning=FALSE}
library(ggplot2)
library(dplyr)
library(tidyr)

df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayes_indices.csv")
```

For the sake of time and computational space, we downloaded the data from github. However, you can find the code to generate it [here](https://easystats.github.io/circus/articles/bayesian_indices.html).




## Results

### Relationship with the Theoretical True Effect

We subtracted the theoretical "true" effect (0 or 1) from the indices, to observe the influence of noise, sample size and priors. 


#### Sensitivity to Noise

```{r, message=FALSE, warning=FALSE}
df %>%
  select(noise, beta, prior_correct, effect, median, mean, map) %>%
  gather(estimate, value, -noise, -effect, -prior_correct) %>%
  mutate(noise = as.factor(noise),
         prior_correct = as.factor(prior_correct),
         value = value-effect) %>%
  ggplot(aes(x = noise, y = value, fill = estimate, color=prior_correct)) +
  geom_boxplot() +
  geom_hline(yintercept = 0) +
  theme_classic() +
  scale_fill_manual(values = c("beta" = "#607D8B", "map" = "#795548", "mean" = "#FF9800", "median" = "#FFEB3B"),
                    name = "Index") +
  scale_color_manual(values = c(`0`="#f44336", `1`="#8BC34A"),
                     name = "Correct Prior") +
  xlab("Point-estimate of the true value 0\n") +
  ylab("\nNoise") +
  coord_cartesian(ylim=c(-1, 1))
```


#### Sensitivity to Sample Size

```{r, message=FALSE, warning=FALSE}
df %>%
  select(sample_size, beta, effect, prior_correct, median, mean, map) %>%
  gather(estimate, value, -sample_size, -effect, -prior_correct) %>%
  mutate(sample_size = as.factor(sample_size),
         prior_correct = as.factor(prior_correct),
         value = value-effect) %>%
  ggplot(aes(x = sample_size, y = value, fill = estimate, color=prior_correct)) +
  geom_boxplot() +
  geom_hline(yintercept = 0) +
  theme_classic() +
  scale_fill_manual(values = c("beta" = "#607D8B", "map" = "#795548", "mean" = "#FF9800", "median" = "#FFEB3B"),
                    name = "Index") +
  scale_color_manual(values = c(`0`="#f44336", `1`="#8BC34A"),
                     name = "Correct Prior") +
  ylab("Point-estimate of the true value 0\n") +
  xlab("\nSample Size") +
  coord_cartesian(ylim=c(-1, 1))
```



#### Statistical Modelling

We fitted a (frequentist) multiple linear regression to statistically test the effects and interactions of noise, sample size and priors on the effect estimates.

```{r, message=FALSE, warning=FALSE}
df %>%
  select(sample_size, beta, effect, prior_correct, median, mean, map, noise) %>%
  gather(estimate, value, -sample_size, -effect, -prior_correct, -noise) %>%
  mutate(noise= scale(noise),
         sample_size = scale(sample_size),
         prior_correct = as.factor(prior_correct)) %>%
  glm(effect ~ estimate/value * noise * sample_size * prior_correct, data=., family="binomial") %>%
  broom::tidy() %>%
  select(term, estimate, p=p.value) %>%
  filter(stringr::str_detect(term, 'value')) %>%
  mutate(term = stringr::str_remove(term, ":value"),
         term = stringr::str_remove(term, "estimate"),
         p = ifelse(p < .001, "< .001***", ifelse(p < .01, "< .01**", ifelse(p < .05, "< .05*", "> .05")))) %>%
  filter(stringr::str_detect(term, ':')) %>%
  knitr::kable(digits=2)
```





### Relationship with the Frequentist Beta

In the next section, we will compare the three Bayesian indices with the frequentist beta.


```{r, message=FALSE, warning=FALSE}
df %>%
  select(sample_size, beta, effect, prior_correct, median, mean, map, noise) %>%
  gather(estimate, value, -beta, -sample_size, -effect, -prior_correct, -noise) %>%
  mutate(effect = as.factor(effect),
         sample_size = as.factor(sample_size),
         estimate = factor(estimate, levels=c("mean", "median", "map"))) %>%
  ggplot(aes(x = beta, y = value, color = effect, shape=sample_size)) +
  geom_point(alpha=0.05) +
  facet_wrap(~estimate, scales = "free") +
  theme_classic() +
  theme(strip.background = element_blank()) +
  scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
  guides(colour = guide_legend(override.aes = list(alpha = 1)),
         shape = guide_legend(override.aes = list(alpha = 1), title="Sample Size"))
```

## Discussion

Conclusions can be found in the [guidelines section](https://easystats.github.io/bayestestR/articles/guidelines.html).
back to top