ggstatsplot: ggplot2 Based Plots with Statistical Details

# <code>ggstatsplot</code>: <code>ggplot2</code> Based Plots with Statistical Details
## An Introductory Tutorial
### Indrajeet Patil
### 2019-09-11 Follow along: <a href="https://tinyurl.com/y376yz6o" class="uri">https://tinyurl.com/y376yz6o</a>

---

---

**Plan**

---

- Why *ggstatsplot*?

- Current functions

- Benefits and scalability

- Limitations

---

# Why *ggstatsplot*?

---

# Raison d'être

---

.right-column[.font150[
Current count of packages on the Comprehensive R Archive Network (**CRAN**) 
** > 13,000**
]
.footnote[<https://cran.r-project.org/web/packages/>]
]

.right-column[.font150[
Short answer: `ggstatsplot` provides a collection of *information-rich* plots with *statistical details* and is suitable for scholarly publications and quick (exploratory) statistical analysis.
]
]

---

# 1. Simpler data analysis workflow

---

In a typical *exploratory* data analysis workflow, *data visualization* and *statistical modeling* are two different phases: visualization informs modeling,
and modeling in its turn can suggest a different visualization method, and so
on and so forth.

The central idea of **ggstatsplot** is simple: combine these two phases into one
in the form of graphics with statistical details, which makes data exploration
simpler and faster.

---

# 2. Information-rich graphic is worth a thousand words

---

.footnote[[(Matejka & Fitzmaurice, *Autodesk Research*, 2017)](https://www.autodeskresearch.com/publications/samestats)]

---

# 3. Ready-made plot = no customization

The **grammar of graphics** (implemented in `ggplot2`) is an incredibly powerful
framework to prepare graphics and can help you make infinite number of graphics,
each tailored for your specific data visualization problem! But...

---

# 4. Consistent API =  no cognitive fatigue

```r
stats::lm(formula = wt ~ mpg, data = mtcars)
```

```r
stats::cor(x = mtcars$wt, y = mtcars$mpg)
```

```r
stats::cor.test(formula = ~ wt + mpg, data = mtcars)
```

**All** functions in `ggstatsplot`-

1.  consistently rely on a dataframe (e.g., `data, x, y, ...`)
2.  expect the data to be in tidy format (Wickham, [2014](https://vita.had.co.nz/papers/tidy-data.pdf))
3.  accept both `"quoted"` and `unquoted` arguments

---

# 5. Follows best practices for statistical reporting

For all statistical tests reported in the plots, the default template abides
by the [APA](https://my.ilstu.edu/~jhkahn/apastats.html) gold standard for
statistical reporting. For example, here are results from a robust *t*-test:

![](images/stats_reporting_format.png)

---

# Get Started

---

# Installation

Install the **ggstatsplot** package from 
[CRAN](https://cran.r-project.org/web/packages/ggstatsplot/index.html):

```r
utils::install.packages("ggstatsplot")
```
--

Then you can get the development version of the package from
[Github](https://github.com/IndrajeetPatil/ggstatsplot):

```r
library(remotes)
remotes::install_github("IndrajeetPatil/ggstatsplot", dependencies = FALSE)
```

Load the needed packages-

```r
library(ggstatsplot)
library(ggplot2)
```

You are recommended to use the [RStudioIDE](https://www.rstudio.com/products/rstudio/), 
but you don't have to.

---

# ggbetweenstats
For between group/condition comparisons

---

# ggbetweenstats - defaults

---

```r
ggbetweenstats(
  data = movies_long,
* x = mpaa, # > 2 groups
  y = rating,
* type = "p", # default
  messages = FALSE
)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbetweenstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbetweenstats - little code, rich details

---

```r
ggbetweenstats(
  data = movies_long,
  x = mpaa,
  y = rating
)
```

- statistical details

- Bayes Factor

- sample sizes

- distribution summary
]
]

---

# ggbetweenstats - pairwise comparisons

---

```r
ggbetweenstats(
  data = movies_long,
  x = mpaa,
  y = rating,
* type = "np",
* mean.ci = TRUE,
* pairwise.comparisons = TRUE,
* pairwise.display = "ns",
* p.adjust.method = "fdr",
  messages = FALSE
)
```
.font80[
Changing pairwise comparisons displayed

- `"ns"` → only **non-significant**

- `"s"` → only **significant**

- `"all"` → **everything**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbetweenstats_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbetweenstats - (aesthetic) changes & outlier tagging

---

```r
ggbetweenstats(
  data = movies_long,
  x = mpaa,
  y = rating,
* type = "r",
* conf.level = 0.99,
  pairwise.comparisons = TRUE,
  pairwise.annotation = "p", 
* outlier.tagging = TRUE,
* outlier.label = title,
* outlier.coef = 2,
* ggtheme = hrbrthemes::theme_ipsum_tw(),
* palette = "Darjeeling2",
* package = "wesanderson",
  messages = FALSE
)
```

The default palette used
is **colorblind-friendly**.
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbetweenstats_3-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbetweenstats - modification with *ggplot2*

---

```r
ggbetweenstats(
  data = movies_long,
  x = mpaa,
  y = rating,
* type = "bf",
  messages = FALSE
) + 
* scale_y_continuous(sec.axis = dup_axis())
```

.font80[
**Note**: You can modify 
all `ggstatsplot` plots further 
using `ggplot2` functions.
Yaay!
]

]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbetweenstats_4-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbetweenstats - 2 groups

---

```r
ggbetweenstats(
  data = mtcars,
* x = am, # 2 groups
  y = wt,
* type = "p", # default
  messages = FALSE
)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbetweenstats_5-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Summary of tests - *ggbetweenstats*

---

.font100[
Type | No. of groups | Test
----------- | -- | --------------------------
**Parametric** | > 2 | Fisher's or Welch's one-way ANOVA
**Parametric** | 2 | Student's or Welch's *t*-test
 **Non-parametric** | > 2 | Kruskal–Wallis one-way ANOVA
 **Non-parametric** | 2 | Mann–Whitney *U* test
**Robust** | > 2 | Heteroscedastic one-way ANOVA for trimmed means
**Robust** | 2 | Yuen's test for trimmed means
**Bayes Factor** | > 2 | Student's ANOVA
**Bayes Factor** | 2 | Student's *t*-test
]

---

# Effect sizes + CI - *ggbetweenstats*

---

Type | Levels | Test | CI?
------------ | -- | -------------------------- | ---- 
**Parametric** | > 2 | `$\eta_{p}^2$`, `$\eta^2$`, `$\omega_{p}^2$`, `$\omega^2$` | Yes
**Parametric** | 2 | Cohen's *d*, Hedge's *g* (central-and noncentral-*t* distribution based) | Yes
**Non-parametric** | > 2 | `$\epsilon^2$` | Yes
**Non-parametric** | 2 | *r* (computed as `$Z/\sqrt{N_{obs}}$`) | Yes
**Robust** | > 2 | `$\xi$` (Explanatory measure of effect size) | Yes
**Robust** | 2 | `$\xi$` (Explanatory measure of effect size) | Yes
**Bayes Factor** | > 2 | No | No
**Bayes Factor** | 2 | No | No

---

# Pairwise comparion tests - *ggbetweenstats*

---

Type | Equal variance? | Test | *p*-value adjustment? 
----------- | --- | ------------------------- | ---
**Parametric** | No | Games-Howell test | Yes 
**Parametric** | Yes | Student's *t*-test | Yes 
 **Non-parametric** | No | Dwass-Steel-Crichtlow-Fligner test | Yes 
**Robust** | No | Yuen's trimmed means test | Yes 
**Bayes Factor** | No | No | No
**Bayes Factor** | Yes | No | No

---

# ggwithinstats
For within group/condition comparisons

---

# ggwithinstats - repeated measures equivalent

---

```r
ggwithinstats(
  data = WRS2::WineTasting,
  x = Wine,
  y = Taste,
* pairwise.comparisons = TRUE,
  pairwise.annotation = "p",
* ggtheme = hrbrthemes::theme_ipsum_tw(),
  ggstatsplot.layer = FALSE,
  messages = FALSE
)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggwithinstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggwithinstats - little code, rich details

---

```r
ggwithinstats(
  data = WRS2::WineTasting,
* x = Wine, # > 2 groups
  y = Taste,
* pairwise.comparisons = TRUE,
  pairwise.annotation = "p",
* ggtheme = hrbrthemes::theme_ipsum_tw(),
  ggstatsplot.layer = FALSE,
  messages = FALSE
)
```

- statistical details

- pairwise comparisons

- sample sizes

- distribution summary
]
]

---

# ggwithinstats - two groups

---

```r
ggwithinstats(
  data = iris_long,
* x = attribute, # 2 groups
  y = value,
* type = "r",
  messages = FALSE
)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggwithinstats_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Summary of tests - *ggwithinstats*

---

Type | No. of groups | Test
----------- | --- | -------------------------
**Parametric** | > 2 | One-way repeated measures ANOVA
**Parametric** | 2 | Student's *t*-test
**Non-parametric** | > 2 | Friedman rank sum test
**Non-parametric** | 2 | Mann–Whitney *U* test
**Robust** | > 2 | Heteroscedastic one-way repeated measures ANOVA for trimmed means
**Robust** | 2 | Yuen's test on trimmed means for dependent samples
**Bayes Factor** | > 2 | One-way repeated measures ANOVA
**Bayes Factor** | 2 | Student's *t*-test

---

# Effect sizes + CI - *ggwithinstats*

---

Type | No. of groups | Test | CI?
----------- | --- | ------------------------- | --- 
**Parametric** | > 2 | `$\eta_{p}^2$`, `$\omega^2$` | Yes
**Parametric** | 2 | Cohen's *d*, Hedge's *g* (central-and noncentral-*t* distribution based) | Yes 
**Non-parametric** | > 2 | `$W_{Kendall}$` (Kendall's coefficient of concordance) | Yes
**Non-parametric** | 2 | *r* (computed as `$Z/\sqrt{N_{pairs}}$`) | Yes
**Robust** | > 2 | No | Yes
**Robust** | 2 | `$\xi$` (Explanatory measure of effect size) | Yes
**Bayes Factor** | > 2 | No | No
**Bayes Factor** | 2 | No | No

---

# Pairwise comparion tests - *ggwithinstats*

---

Type | Test | *p*-value adjustment? 
----------- | ---------------------------- | ---
**Parametric** | Student's *t*-test | Yes 
 **Non-parametric** | Durbin-Conover test | Yes 
**Robust** | Yuen's trimmed means test | Yes 
**Bayes Factor** | No | No

---

# ggscatterstats
Association between two numeric variables

---

# ggscatterstats - defaults

---

```r
ggscatterstats(
 data = movies_long,
 x = budget,
 y = rating,
 type = "p", # default #<<<
* conf.level = 0.99,
 messages = FALSE
)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggscatterstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggscatterstats - little code, rich details

---

```r
ggscatterstats(
  data = movies_long,
  x = budget,
  y = rating,
  conf.level = 0.99
)
```

- distribution

- Bayes Factor

- statistical details
]
]

---

# ggscatterstats - conditional point tagging

---

```r
ggscatterstats(
  data = movies_long,
  x = budget,
  y = rating,
  type = "r", 
* centrality.para = "mean",
* label.var = title,
* label.expression = budget > 150
*                    & rating > 7.5,
* marginal.type = "density",
  messages = FALSE
)
```

- **histogram**
- **boxplot**
- **density**
- **violin**
- **densigram**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggscatterstats_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Summary of tests - *ggscatterstats*

---

Type | Test | CI?
----------- | ------------------------- | ---
**Parametric** | Pearson's correlation coefficient | Yes
 **Non-parametric** | Spearman's rank correlation coefficient | Yes
**Robust** | Percentage bend correlation coefficient | Yes
**Bayes Factor** | Pearson's correlation coefficient | No

---

# ggscatterstats - changing smoothing functions

---

```r
ggscatterstats(
  data = movies_long,
  x = budget,
  y = rating,
  marginal = FALSE,  
* method = "gam",
* formula = y ~ s(x, k = 3),
* centrality.para = "mean",
  messages = FALSE
)
```

- **mean**

- **median**

]
]

---

# ggcorrmat
Association between multiple numeric variables

---

# ggcorrmat - defaults

---

```r
ggcorrmat(data = dplyr::starwars)
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → not implemented

]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggcorrmat_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggcorrmat - little code, rich details

---

```r
ggcorrmat(data = dplyr::starwars)
```

- statistical details

- sample sizes

- details about test

**Note**: Informative label
about sample sizes in case
`NA`s are present.
]
]

---

# ggcorrmat - changing defaults

---

```r
ggcorrmat(
  data = ggplot2::msleep,
  cor.vars = sleep_cycle:bodywt,
  type = "r",
* matrix.type = "upper",
* p.adjust.method = "holm",
  colors = NULL,
* package = "yarrr",
* palette = "southpark"
)
```

.font80[
In addition to `output = "plot"`, 
this function can also be used to
get a dataframe of results:

- `"r"` → **correlation**

- `"p"` → **p-values**

- `"n"` → **sample sizes**

- `"ci"` → **confidence intervals**

]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggcorrmat_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggcorrmat - dataframe outputs

---

```r
ggcorrmat(
  data = ggplot2::msleep,
  cor.vars = sleep_cycle:bodywt,
* type = "np",
* output = "r",
  p.adjust.method = "fdr",
  messages = FALSE
)
```

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> variable </th>
 <th style="text-align:right;"> sleep_cycle </th>
 <th style="text-align:right;"> awake </th>
 <th style="text-align:right;"> brainwt </th>
 <th style="text-align:right;"> bodywt </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> sleep_cycle </td>
 <td style="text-align:right;"> 1.00 </td>
 <td style="text-align:right;"> 0.49 </td>
 <td style="text-align:right;"> 0.87 </td>
 <td style="text-align:right;"> 0.85 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> awake </td>
 <td style="text-align:right;"> 0.49 </td>
 <td style="text-align:right;"> 1.00 </td>
 <td style="text-align:right;"> 0.59 </td>
 <td style="text-align:right;"> 0.53 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> brainwt </td>
 <td style="text-align:right;"> 0.87 </td>
 <td style="text-align:right;"> 0.59 </td>
 <td style="text-align:right;"> 1.00 </td>
 <td style="text-align:right;"> 0.96 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> bodywt </td>
 <td style="text-align:right;"> 0.85 </td>
 <td style="text-align:right;"> 0.53 </td>
 <td style="text-align:right;"> 0.96 </td>
 <td style="text-align:right;"> 1.00 </td>
 </tr>
</tbody>
</table>

</div>
]
]

```r
ggcorrmat(
  data = ggplot2::msleep,
  cor.vars = sleep_cycle:bodywt,
* type = "np",
* output = "p",
  p.adjust.method = "none",
  messages = FALSE
)
```

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> variable </th>
 <th style="text-align:right;"> sleep_cycle </th>
 <th style="text-align:right;"> awake </th>
 <th style="text-align:right;"> brainwt </th>
 <th style="text-align:right;"> bodywt </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> sleep_cycle </td>
 <td style="text-align:right;"> 0.0000000 </td>
 <td style="text-align:right;"> 0.0045296 </td>
 <td style="text-align:right;"> 0.0e+00 </td>
 <td style="text-align:right;"> 0e+00 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> awake </td>
 <td style="text-align:right;"> 0.0045296 </td>
 <td style="text-align:right;"> 0.0000000 </td>
 <td style="text-align:right;"> 1.4e-06 </td>
 <td style="text-align:right;"> 2e-07 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> brainwt </td>
 <td style="text-align:right;"> 0.0000000 </td>
 <td style="text-align:right;"> 0.0000014 </td>
 <td style="text-align:right;"> 0.0e+00 </td>
 <td style="text-align:right;"> 0e+00 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> bodywt </td>
 <td style="text-align:right;"> 0.0000000 </td>
 <td style="text-align:right;"> 0.0000002 </td>
 <td style="text-align:right;"> 0.0e+00 </td>
 <td style="text-align:right;"> 0e+00 </td>
 </tr>
</tbody>
</table>

</div>
]
]

---

# ggcorrmat - dataframe outputs

---

```r
ggcorrmat(
  data = ggplot2::msleep,
  cor.vars = sleep_cycle:bodywt,
* type = "np",
* output = "n",
  p.adjust.method = "fdr",
  messages = FALSE
)
```

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> variable </th>
 <th style="text-align:right;"> sleep_cycle </th>
 <th style="text-align:right;"> awake </th>
 <th style="text-align:right;"> brainwt </th>
 <th style="text-align:right;"> bodywt </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> sleep_cycle </td>
 <td style="text-align:right;"> 32 </td>
 <td style="text-align:right;"> 32 </td>
 <td style="text-align:right;"> 30 </td>
 <td style="text-align:right;"> 32 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> awake </td>
 <td style="text-align:right;"> 32 </td>
 <td style="text-align:right;"> 83 </td>
 <td style="text-align:right;"> 56 </td>
 <td style="text-align:right;"> 83 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> brainwt </td>
 <td style="text-align:right;"> 30 </td>
 <td style="text-align:right;"> 56 </td>
 <td style="text-align:right;"> 56 </td>
 <td style="text-align:right;"> 56 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> bodywt </td>
 <td style="text-align:right;"> 32 </td>
 <td style="text-align:right;"> 83 </td>
 <td style="text-align:right;"> 56 </td>
 <td style="text-align:right;"> 83 </td>
 </tr>
</tbody>
</table>

</div>
]
]

```r
options(digits = 3)
ggcorrmat(
  data = ggplot2::msleep,
  cor.vars = awake:bodywt,
* type = "np",
* output = "ci",
  p.adjust.method = "fdr",
  messages = FALSE
)
```

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> pair </th>
 <th style="text-align:right;"> r </th>
 <th style="text-align:right;"> lower </th>
 <th style="text-align:right;"> upper </th>
 <th style="text-align:right;"> p </th>
 <th style="text-align:right;"> lower.adj </th>
 <th style="text-align:right;"> upper.adj </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> awake-brainwt </td>
 <td style="text-align:right;"> 0.594 </td>
 <td style="text-align:right;"> 0.392 </td>
 <td style="text-align:right;"> 0.741 </td>
 <td style="text-align:right;"> 0 </td>
 <td style="text-align:right;"> 0.340 </td>
 <td style="text-align:right;"> 0.767 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> awake-bodywt </td>
 <td style="text-align:right;"> 0.535 </td>
 <td style="text-align:right;"> 0.360 </td>
 <td style="text-align:right;"> 0.673 </td>
 <td style="text-align:right;"> 0 </td>
 <td style="text-align:right;"> 0.318 </td>
 <td style="text-align:right;"> 0.698 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> brainwt-bodywt </td>
 <td style="text-align:right;"> 0.957 </td>
 <td style="text-align:right;"> 0.928 </td>
 <td style="text-align:right;"> 0.975 </td>
 <td style="text-align:right;"> 0 </td>
 <td style="text-align:right;"> 0.919 </td>
 <td style="text-align:right;"> 0.978 </td>
 </tr>
</tbody>
</table>

</div>
]
]

---

# gghistostats
Distribution of a numeric variable

---

# gghistostats - defaults

---

```r
gghistostats(
  data = movies_long,    
  x = budget,
* test.value = 50,
  messages = FALSE
) 
```

- `"p"` → **parametric**

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/gghistostats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# gghistostats - further customization

---

```r
gghistostats(
  data = movies_long,    
  x = budget,
  effsize.type = "d",
  test.value = 50,
* bar.measure = "mix",
  centrality.para = "median", 
* test.value.line = TRUE,
* normal.curve = TRUE,
  ggtheme = hrbrthemes::theme_ipsum_tw(), 
  ggstatsplot.layer = FALSE, 
  messages = FALSE
) 
```

- **count**
- **proportion**
- **both** (of the above)
- **density**

]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/gghistostats_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# gghistostats - little code, rich details

---

```r
gghistostats(
  data = movies_long,    
  x = budget,
  effsize.type = "d",
  test.value = 50,
  test.value.size = TRUE, 
  bar.measure = "mix", 
  centrality.para = "median", 
  test.value.line = TRUE,
  normal.curve = TRUE
)
```

- statistical details

- Bayes Factor

- frequency

- distribution summary
]
]

---

# ggdotplotstats
Distribution of a numeric variable with labels

---

# ggdotplotstats - defaults

---

```r
ggdotplotstats(
  data = movies_long,
  x = budget,
  y = genre,
  effsize.type = "d", 
* test.value = 52,
  centrality.para = "median", 
* test.value.line = TRUE,
  test.value.color = "red", 
  ggtheme = ggthemes::theme_par(),
  messages = FALSE
) 
```

- `"np"` → **non-parametric**

- `"r"` → **robust**

- `"bf"` → **bayes factor**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggdotplotstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggdotplotstats - little code, rich details

---

```r
ggdotplotstats(
  data = movies_long,
  x = budget,
  y = genre,
  effsize.type = "d",
  test.value = 52,     
  centrality.para = "median", 
  test.value.line = TRUE, 
  test.value.color = "red"
) 
```

- statistical details

- Bayes Factor

- distribution summary
]
]

---

# Summary of tests - *gghistostats*/*ggdotplotstats*

---

Type | Test
------------------ | -------------------------
**Parametric** | One-sample Student's *t*-test
**Non-parametric** | One-sample Wilcoxon test 
**Robust** | One-sample percentile bootstrap
**Bayes Factor** | One-sample Student's *t*-test

Type | Effect size | CI?
----------- | ------------------------- | --- 
**Parametric** | Cohen's *d*, Hedge's *g* (central-and noncentral-*t* distribution based) | Yes
**Non-parametric** | *r* (computed as `$Z/\sqrt{N_{obs}}$`) | Yes
**Robust** | `$M_{robust}$` (a robust location measure) | Yes
**Bayes Factor** | No | No

---

# ggpiestats
For composition of categorical variables

---

# ggpiestats - defaults

---

```r
# let's use subset of data
ggpiestats(
  data = dplyr::filter(.data = movies_long, 
  genre %in% c("Drama", "Comedy", "Animated")), 
  x = genre,
  y =  mpaa,
* paired = FALSE, # default
* conf.level = 0.99,
  package = "ggsci",
  palette = "default_ucscgb",
  messages = FALSE
)
```

Test by design

- `paired = FALSE` → Pearson's `$\chi^2$`

- `paired = TRUE` → McNemar

]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggpiestats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggpiestats - little code, rich details

---

```r
# let's use subset of data
ggpiestats(
  data = dplyr::filter(movies_long, 
  genre %in% c("Drama", "Comedy", "Animated")), 
  x = genre,
  y =  mpaa,
* conf.level = 0.99
)
```

- statistical details

- Bayes Factor

- sample sizes

- proportion test results
]
]

---

# ggpiestats - proportion test

---

```r
ggpiestats(
  data = as.data.frame(Titanic), 
* x = Survived,
* counts = Freq,
* slice.label = "both",
  messages = FALSE
)
```

Test by analysis

- `condition != NULL` → contingency table

- `y = = NULL` → goodness of fit

]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggpiestats_2-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbarstats
For composition of categorical variables

---

# ggbarstats - defaults

---

```r
ggbarstats(
  data = movies_long, 
  x = genre,
  y =  mpaa,
* paired = FALSE, # default
  package = "ggsci",
  palette = "default_igv",
* caption = substitute(
*   paste(italic("Source"), ": www.imdb.com")
* ),
  messages = FALSE
)
```

.font70[
**Note**: Even if you display 
Bayes Factor message in a caption, 
you can still use the `caption`
argument.

Label information

- **percentage** (default) 
- **counts**
- **both** (of the above)
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggbarstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggbarstats - little code, rich details

---

```r
ggbarstats(
  data = movies_long, 
  x = genre,
  y =  mpaa
)
```

- statistical details

- Bayes Factor

- sample sizes

- proportion test results
]
]

---

# Test summary - *ggpiestats*/*ggbarstats*

---

**Tests**

Data | Design | Test
----------- | ------------  | -------------------------
Unpaired | `$n \times p$` contingency table | Pearson's `$\chi^2$` test
Paired | `$n \times p$` contingency table | McNemar's test
Frequency | `$n \times 1$` contingency table | Goodness of fit

**Effect sizes + CI**

Test | Effect size | CI?
----------------------- | ----------------------------- | -----
Pearson's `$\chi^2$` test | Cramer's `$V$` | Yes
McNemar's `$\chi^2$` test | Cohen's `$g$` | Yes
Goodness of fit `$\chi^2$` test | Cramer's `$V$` | Yes

---

# ggcoefstats
Displaying results from regression analyses

---

# ggcoefstats - defaults

---

```r
# model
mod <- stats::aov(
 formula = rating ~ mpaa * genre,
 data = movies_long
)

# plot
ggcoefstats(x = mod)
```

.font80[
In addition to `output = "plot"`, 
this function can also be used to
get a dataframe of results:

- `"tidy"` → **estimates**

- `"glance"` → **model summary**

- `"augment"` → **predictions**
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/ggcoefstats_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# ggcoefstats - little code, rich details

---

```r
# model
mod <- stats::aov(
 formula = rating ~ mpaa * genre,
 data = movies_long
)

# plot
ggcoefstats(x = mod)
```

- estimate + 95% CI

- model summary

- statistical details
]
]

---

# ggcoefstats - dataframe outputs

---

```r
library(lme4)

# model
mod1 <- 
 lme4::lmer(formula = Reaction ~ Days + (Days | Subject),
 data = sleepstudy)

# dataframe
ggcoefstats(x = mod1, 
*           output = "glance")
```

<table>
 <thead>
 <tr>
 <th style="text-align:right;"> sigma </th>
 <th style="text-align:right;"> logLik </th>
 <th style="text-align:right;"> AIC </th>
 <th style="text-align:right;"> BIC </th>
 <th style="text-align:right;"> REMLcrit </th>
 <th style="text-align:right;"> df.residual </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:right;"> 25.59181 </td>
 <td style="text-align:right;"> -871.8141 </td>
 <td style="text-align:right;"> 1755.628 </td>
 <td style="text-align:right;"> 1774.786 </td>
 <td style="text-align:right;"> 1743.628 </td>
 <td style="text-align:right;"> 174 </td>
 </tr>
</tbody>
</table>

</div>
]
]

```r
library(ordinal)

# model
mod2 <- clm(formula = rating ~ temp * contact, 
 data = wine)

# dataframe
ggcoefstats(x = mod2, 
*           output = "augment") %>%
  head(5)
```

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> rating </th>
 <th style="text-align:left;"> temp </th>
 <th style="text-align:left;"> contact </th>
 <th style="text-align:right;"> .fitted </th>
 <th style="text-align:right;"> .se.fit </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> 2 </td>
 <td style="text-align:left;"> cold </td>
 <td style="text-align:left;"> no </td>
 <td style="text-align:right;"> 0.5622964 </td>
 <td style="text-align:right;"> 0.0884961 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 3 </td>
 <td style="text-align:left;"> cold </td>
 <td style="text-align:left;"> no </td>
 <td style="text-align:right;"> 0.2086491 </td>
 <td style="text-align:right;"> 0.0788338 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 3 </td>
 <td style="text-align:left;"> cold </td>
 <td style="text-align:left;"> yes </td>
 <td style="text-align:right;"> 0.4346731 </td>
 <td style="text-align:right;"> 0.0837035 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 4 </td>
 <td style="text-align:left;"> cold </td>
 <td style="text-align:left;"> yes </td>
 <td style="text-align:right;"> 0.0893885 </td>
 <td style="text-align:right;"> 0.0436085 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 4 </td>
 <td style="text-align:left;"> warm </td>
 <td style="text-align:left;"> no </td>
 <td style="text-align:right;"> 0.1902823 </td>
 <td style="text-align:right;"> 0.0711286 </td>
 </tr>
</tbody>
</table>

</div>
]
]

---

# *ggcoefstats*: Supported models

---

.font90[
.pull-left[
.pull-left[
- `aareg`
- `anova`
- `aov`
- `aovlist`
- `Arima` 
- `bglmerMod`
- `biglm`
- `blmerMod`
- `brmsfit`
- `btergm`
- `cch`
- `clm`
- `clmm`
- `confusionMatrix`
- `coxph`
- `data.table`
- `drc`
- `epi.2by2`
- `ergm`
]
.pull-right[
- `felm`
- `fitdistr`
- `glmerMod`
- `glmmPQL`
- `glmmTMB`
- `gls`
- `gam`
- `Gam`
- `gamlss`
- `garch`
- `glm` 
- `glmmadmb` 
- `glmmTMB`
- `glmrob`
- `glmRob`
- `gmm`
- `ivreg`
]
]

.pull-right[
.pull-left[
- `lm`
- `lm.beta`
- `lme`
- `lmerMod`
- `lmodel2`
- `lmrob`
- `lmRob`
- `mclogit` 
- `mcmc`
- `MCMCglmm`
- `mediate` 
- `mjoint`
- `mle2`
- `mlm`
- `mmclogit`
- `multinom`
- `negbin`
- `nlmerMod`
- `nlrq`
- `nls`
]
.pull-right[
- `orcutt`
- `plm`
- `polr`
- `ridgelm`
- `rjags`
- `rlm`
- `rlmerMod`
- `rq`
- `speedglm`
- `speedlm`
- `stanreg`
- `survreg`
- `svyglm`
- `svyolr`
- `svyglm`
- `tobit`
- `wblm`
]
]
]

---

# *ggcoefstats*: If not implemented, use a dataframe

---

```r
# dataframe with results
df <- tibble::tribble(
 ~term, ~estimate, ~std.error, ~statistic, ~p.value,
 "(Intercept)", 3.77, 0.165, 22.9, 1.49e-20,
 "x", -1.36, 0.258, -5.26, 1.13e-5
)

# plot
# `statistic` argument decides label format
ggcoefstats(
  x = df,
* statistic = "z",
  exclude.intercept = FALSE
)
```

At the minimum, two
columns needed - term and estimate.
]

]

---

# *ggcoefstats*: You can also do meta-analysis!

---

```r
# made up data
meta_df <- tibble::tribble(
 ~term, ~estimate, ~std.error,
 "study_1", 0.111, 0.065,
 "study_2", -0.003, 0.258,
 "study_3", 0.001, 0.120,
 "study_4", 0.032, 0.022,
 "study_5", -0.765, 0.650,
 "study_6", -0.032, 0.058
)

# plot
ggcoefstats(
  x = meta_df,
* meta.analytic.effect = TRUE,
* bf.message = TRUE,
  xlab = "estimate"
)
```

.font80[
- Frequentist random-effects meta-analysis from `metafor`
- Bayesian random-effects meta-analysis from `metaBMA`
]

]

---

# *grouped_* variants of all functions
Running the same function for 
all levels of a single grouping variable

---

# *grouped_* functions

You can repeat the same analysis for all levels of a single grouping variable.

---

```r
# only one additional argument
grouped_ggpiestats(
  data = mtcars, 
  x = cyl,
* grouping.var = am,
* results.subtitle = FALSE,
  messages = FALSE
)
```

.font70[
Available `grouped_` variants
- `grouped_ggdotplotstats`
- `grouped_ggbarstats`    
- `grouped_ggscatterstats`
- `grouped_gghistostats`  
- `grouped_ggpiestats`    
- `grouped_ggbetweenstats`
- `grouped_ggwithinstats`
- `grouped_ggcorrmat`  
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/grouped_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Utility beyond *ggstatsplot*
What if I don't like the default plots 
but still want to display statistical results?

---

# Using as helper functions

---

`ggstatsplot` can also be used just to get the statistical details.

```r
# using `ggstatsplot` for stats
results <- 
 ggstatsplot::ggpiestats( 
 data = Titanic_full,
 x = Survived,
 y = Sex,
* return = "subtitle",
 messages = FALSE
)

# using `ggiraphExtra` for plot
*ggiraphExtra::ggSpine(
  data = Titanic_full,
  aes(x = Sex, fill = Survived),
  addlabel = TRUE,
  interactive = FALSE
*) + labs(subtitle = results)
```

.font70[
 
**All** included analyses have their
corresponding [helper functions](https://indrajeetpatil.github.io/ggstatsplot/reference/index.html#section-helper-functions-for-preparing-statistics-subtitles) for preparing 
subtitles with statistical details.
]
]

.right-plot[
<img src="ggstatsplot_presentation_files/figure-html/subtitle_1-1.png" width="100%" style="display: block; margin: auto;" />
]

---

# Glossary
Statistical reporting in *ggstatsplot*

---

# Best practices in reporting statistical details

- As discussed before, the details included in statistical analyses follow the
APA gold standard.

- The default tests follow the best practices. For example, `ggbetweenstats`
function by default runs Welch's *t*-test and Welch's ANOVA - and not Student's *t*-test and Fisher's
ANOVA - based on recent work (Delacre et al.,
[2017](https://www.rips-irsp.com/article/10.5334/irsp.82/),
[2018](https://psyarxiv.com/wnezg)).

- No *p*-value error 
(Lilienfeld et al.,
[2015](https://www.frontiersin.org/articles/10.3389/fpsyg.2015.01100/full))

---

# Avoiding errors

Since the plot and the statistical analysis are yoked together, the chances of
making an error in reporting the results are minimized. You never have to write
the results manually or copy-paste them from someplace else.

![](images\statcheck.png)

.footnote[[(Nuijten et al., *Behavior Research Methods*, 2016)](https://link.springer.com/article/10.3758/s13428-015-0664-2)]

---

# Making sense of null results

Combination of frequentist and Bayesian statistics for each analysis to properly
interpret the null results.

![](images/null_results.png)

.footnote[[(Aczel et al., *AMPPS*, 2018)](https://journals.sagepub.com/doi/pdf/10.1177/2515245918773742)]

---

# Toggling between type of statistics

```r
# anova
ggbetweenstats(
  data = mtcars,
  x = cyl,
  y = wt,
* type = "p"
)

# correlation analysis
ggscatterstats(
  data = mtcars,
  x = wt,
  y = mpg,
* type = "p"
)

# t-test
gghistostats(
  data = mtcars,
  x = wt,
  test.value = 2,
* type = "p"
)
```
]

```r
# anova
ggbetweenstats(
  data = mtcars,
  x = cyl,
  y = wt,
* type = "np"
)

# correlation analysis
ggscatterstats(
  data = mtcars,
  x = wt,
  y = mpg,
* type = "np"
)

# t-test
gghistostats(
  data = mtcars,
  x = wt,
  test.value = 2,
* type = "np"
)
```
]

---

# Glossary
Summary of statistical tests included
(as of 0.1.1.9000)

---

# Types of statistical analyses supported

Functions | Description | **Parametric** | **Non-parametric** | **Robust** | **Bayes Factor**
------- | ------------------ | ---- | ----- | ----| ----- 
`ggbetweenstats` | Between group/condition comparisons | Yes | Yes | Yes | Yes
`ggwithinstats` | Within group/condition comparisons | Yes | Yes | Yes | Yes
`gghistostats`, `ggdotplotstats` | Distribution of a numeric variable | Yes | Yes | Yes | Yes
`ggcorrmat` | Correlation matrix | Yes | Yes | Yes | No
`ggscatterstats` | Correlation between two variables | Yes | Yes | Yes | Yes
`ggpiestats`, `ggbarstats` | Association between categorical variables | Yes | `NA` | `NA` | Yes
`ggcoefstats` | Regression model coefficients | Yes | Yes| Yes | Yes

---

# Effect sizes + CI available?

Test | **Parametric** | **Non-parametric** | **Robust** | **Bayes Factor**
---------------------------- | ------ | ------ | ------ | ------
one-sample *t*-test | Yes | Yes | Yes | No
two-sample *t*-test (between) | Yes | Yes | Yes | No
two-sample *t*-test (within) | Yes | Yes | Yes | No
One-way ANOVA (between) | Yes | Yes | Yes | No
One-way ANOVA (within) | Yes | Yes | No | No
correlations | Yes | Yes | Yes | No
contingency table | Yes | `NA` | `NA` | No
goodness of fit | Yes | `NA` | `NA` | No
regression | Yes | Yes | Yes | No

---

# Why use *ggstatsplot*? 
Summary of benefits

---

# Benefits of using *ggstatsplot*

---

- Truly makes your figures worth a thousand words.

- No need to copy-paste results to the text editor (MS-Word, e.g.).

- Disembodied figures stand on their own and are easy to evaluate for the reader.

- More breathing room for theoretical discussion and other text.

- No need to worry about updating figures and statistical details separately
  if something about the data changes.

- Minimal amount of code needed for all functions (typically only `data`, `x`,
  and `y`). This minimizes chances of error.

---

# Example - before and after *ggstatsplot*

---

**Before ggstatsplot**

![](images/before_ggstats.PNG)

![](images/no_joy.jpg)

]

**After ggstatsplot**

![](images/after_ggstats.PNG)
]

---

# Quality assurance

---

`ggstatsplot` includes over **1400** tests to make
sure that all the functions are working properly.

If you use `ggstatsplot` to report results in your publication and find an error
that traces back to code in this package, I will personally write an apology
letter to you and to the editor of the journal in question.

![](images/codecov.png)

---

# Informative plot designs

---

.pull-left[
.font150[
The default plots in *ggstatsplot* are
**opinionated**, yes, but an attempt has been made to make sure that they
follow best principles outlined in the data visualization research.
]
]

---

# Documentation

---

Exhaustive documentation at the dedicated website- 
<https://indrajeetpatil.github.io/ggstatsplot/>

![](images/website.PNG)

---

# Limitations

---

# Limitations

---

- Limited kinds of plots available.

- Limited number of statistical tests (and effect
sizes) available.

- Faceting (or small multiples) not implemented.

- Default plots can be too complicated for effectively communicating results in
time-constrained presentation settings (e.g., conference talks).

- Bulky API (in terms of number of function arguments to keep in mind). 
(Saving grace: Defaults are sufficient most of the time.)

---

# Overcoming these limitations

---

![](images/needs_you.jpg)

]

- Read and correct any inconsistencies in the
[documentation](https://indrajeetpatil.github.io/ggstatsplot/)

- Raise issues about bugs or wanted features

- Review code

- Add new functionality (in the form of new plotting functions or helpers for
preparing subtitles)
]
]

---

# Not the only game in town
other options that pack informative graphics with statistical details

---

# GUI softwares

---

**jamovi**

![](images/jamovi.PNG)
]

**JASP**

![](images/jasp.PNG)
]

---

# Acknowledgments

Contributor to *ggstatsplot*: [Chuck Powell](https://github.com/ibecav)

Advisors: [Mina Cikara](http://www.intergroupneurosciencelaboratory.com/) and [Fiery Cushman](http://cushmanlab.fas.harvard.edu/index.php)

Slides created via the R package [xaringan](https://github.com/yihui/xaringan).

The CSS template comes from [Garrick Aden-Buie](https://github.com/gadenbuie/gentle-ggplot2).

---

# Thanks!

.font150.text-white[
[@patilindrajeets](https://twitter.com/patilindrajeets) 
[github.com/IndrajeetPatil](https://github.com/IndrajeetPatil) 
patilindrajeet.science@gmail.com
]