---
title: "Box-Whisker Plots"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Box-Whisker Plots}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  dpi = 300,
  fig.show = "hold",
  fig.align = "center",
  fig.width = 7,
  out.width = "80%",
  comment = "#>"
)
```

# 1. Introduction

This vignette documents and illustrates workflows for producing box-and-whisker plots using the `ospsuite.plots` library.

The function for plotting box-whiskers is `plotBoxWhisker`. Basic documentation of the function can be found using: `?plotBoxWhisker`. The output of the function is a `ggplot` object.

## 1.1 Setup

This vignette uses the `{ospsuite.plots}` and `{tidyr}` libraries. We will use the default settings of `{ospsuite.plots}` (see vignette("ospsuite.plots", package = "ospsuite.plots")) but will adjust the legend position and the alignment of the caption.

```{r setup, warning=FALSE, message=FALSE}
options(rmarkdown.html_vignette.check_title = FALSE)

library(ospsuite.plots)
library(tidyr)

# Set Defaults
oldDefaults <- ospsuite.plots::setDefaults()

# Adjust default theme plot for prettier plots
theme_update(legend.position = "top")
theme_update(legend.direction = "horizontal")
theme_update(plot.caption = element_text(hjust = 1))
```

## 1.2 Example Data

This vignette uses a dataset provided by the package:

```{r load-data, results='asis'}
pkRatioData <- exampleDataCovariates |>
  dplyr::filter(SetID == "DataSet1") |>
  dplyr::select(-c("SetID", "gsd", "AgeBin")) |>
  dplyr::mutate(Agegroup = cut(Age, breaks = c(0, 6, 12, 18, 60), include.lowest = TRUE, labels = c("infants", "school children", "adolescents", "adults")))

knitr::kable(head(pkRatioData), digits = 3, caption = "First rows of example data pkRatioData")
```

Metadata is a list that contains dimension and unit information for dataset columns. If available, axis labels are set by this information.

```{r metaData, results='asis'}
metaData <- attr(exampleDataCovariates, "metaData")
knitr::kable(metaData2DataFrame(metaData), digits = 2, caption = "List of meta data")
```

# 2. Examples

## 2.1 Examples for Aggregation of Categorical Data

### 2.1.1 Minimal Example

Age (mapped to y) is aggregated.

```{r minimal-example, fig.alt="Box plot showing age distribution aggregated across all data. The plot displays quartiles with median line, showing the overall age range distribution."}
plotBoxWhisker(data = pkRatioData, mapping = aes(y = Age), metaData = metaData)
```

### 2.1.2 Stratification on X-axis

Age (mapped to y) is aggregated for different countries (mapped to x).

```{r example-x-stratification, fig.alt="Box plots showing age distribution stratified by country. Separate box plots are displayed for each country on the x-axis, allowing comparison of age distributions between countries."}
plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age
  ),
  data = pkRatioData,
  metaData = metaData
)
```

### 2.1.3 Stratification by Color

Age (mapped to y) is aggregated for different countries (mapped to fill).

```{r minimal-example-color, fig.alt="Box plot showing age distribution stratified by country using color coding. Different countries are represented by different fill colors in a single grouped box plot."}
plotBoxWhisker(
  mapping = aes(
    fill = Country,
    y = Age
  ),
  data = pkRatioData,
  metaData = metaData
)
```

### 2.1.4 Stratification by Color and on X-axis

Age (mapped to `y`) is aggregated for different countries (mapped to `x`) and Sex (mapped to `groupby`). `groupby` is an additional aesthetic of `{ospsuite.plots}` that works together with the variable `groupAesthetics`. For the function `plotBoxWhisker()`, `groupAesthetics` is not settable and is fixed to `fill`.

```{r minimal-example-color-x, fig.alt="Box plots showing age distribution stratified by both country (x-axis) and sex (color/grouping). The plot displays separate box plots for each country, with different colors representing sex within each country."}
plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age,
    groupby = Sex
  ),
  data = pkRatioData,
  metaData = metaData
)
```

### 2.1.5 Stratification by Column Combination

Below, `groupby` is mapped to a combination of the columns "Sex" and "Country":

```{r minimal-example-interaction, fig.alt="Box plot showing age distribution grouped by the interaction of country and sex. Each combination of country and sex is represented as a separate group with distinct colors in the legend."}
plotBoxWhisker(
  mapping = aes(
    y = Age,
    groupby = interaction(Country, Sex, sep = "-")
  ),
  data = pkRatioData,
  metaData = metaData
) + theme(legend.title = element_blank())
```

### 2.1.6 Omit Data Points Flagged as Missing Dependent Variable (MDV)

If some of the data should be omitted, we can do this by mapping a logical column to the aesthetic `mdv`. Below, we exclude data from Germany:

```{r minimal-example-mdv, fig.alt="Box plots showing age distribution by country and sex, with data from Germany excluded using the missing dependent variable (MDV) flag. The plot demonstrates how to omit specific data points from the analysis."}
plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age,
    mdv = Country == "Germany",
    groupby = Sex
  ),
  data = pkRatioData,
  metaData = metaData
)
```

## 2.2 Examples for Box-Whisker vs Numeric Data

### 2.2.1 Numeric Data as Factor

In the next example, we added a numeric column as "mean age" of the age group to the dataset. This column is mapped *as factor* to `x`. The values are now displayed as categorical values equidistant:

```{r numeric-example-factor, fig.alt="Box plots showing ratio distribution by mean age treated as categorical factors. Each mean age value is positioned equidistantly on the x-axis, with different age groups represented by different fill colors."}
pkRatioData <- pkRatioData |>
  dplyr::group_by(Agegroup) |>
  dplyr::mutate(meanAge = round(mean(Age), 2))

metaData[["meanAge"]] <- metaData[["Age"]]
metaData <- metaData[intersect(names(pkRatioData), names(metaData))]

plotBoxWhisker(
  data = pkRatioData,
  mapping = aes(
    x = as.factor(meanAge),
    y = Ratio,
    fill = Agegroup
  ),
  metaData = metaData
)
```

### 2.2.2 Numeric Data as Distinct Numeric Values

If the column mapped to `x` is numeric and not a factor, the x-position of the boxes corresponds to the numeric value:

```{r numeric-example-distinct, fig.alt="Box plots showing ratio distribution by mean age as continuous numeric values. The x-position of each box plot corresponds to its actual numeric mean age value, with different age groups represented by different fill colors."}
plotBoxWhisker(
  mapping = aes(
    x = meanAge,
    y = Ratio,
    fill = Agegroup
  ),
  data = pkRatioData,
  metaData = metaData
)
```

### 2.2.3 Map Continuous Column with Binning

By mapping `x` with a function, it is possible to aggregate data of a continuous column into bins. Below, we use the function `cut` to aggregate the data into 4 distinct bins (see ?cut).

Attention: for `cut(Age)`, no `metaData` exists. So we have to set the x label manually.

```{r numeric-example-bin4, fig.alt="Box plots showing ratio distribution across four automatically generated age ranges. The continuous age variable is binned into four equal intervals, with each bin represented as a separate box plot on the x-axis."}
plotBoxWhisker(
  mapping = aes(
    x = cut(Age, 4),
    y = Ratio
  ),
  data = pkRatioData,
  metaData = metaData
) + labs(x = "Age range [years]")
```

We can use any function that converts a continuous vector to a factor. Below, we define our own cut function simply as a wrapper around `cut`, with some fixed arguments. This function is then mapped to `x` and `groupby`.

```{r numeric-example-binfunction, fig.alt="Box plots showing ratio distribution across predefined age groups using a custom binning function. Ages are categorized into infants, school children, adolescents, and adults, with each group shown as a separate box plot."}
myCutfun <- function(x) {
  cut(x = x, breaks = c(0, 6, 12, 18, 60), include.lowest = TRUE, labels = c("infants", "school children", "adolescents", "adults"))
}

plotBoxWhisker(
  mapping = aes(
    x = myCutfun(Age),
    y = Ratio,
    groupby = myCutfun(Age)
  ),
  data = pkRatioData,
  metaData = metaData
) +
  theme(legend.title = element_blank()) +
  labs(x = "")
```

## 2.3 Aggregation Function

By default, the data is aggregated by percentiles defined by the default option `ospsuite.plots.percentiles`. The percentiles can be customized for a specific plot by the input variable `percentiles` or for all plots generated by `plotBoxWhisker()` by changing the default options (`setOspsuite.plots.option(optionKey = OptionKeys$percentiles, value = c(0.05, 0.25, 0.5, 0.75, 0.95))`).

It is also possible to use a customized function via the input variable `statFun`. If `statFun` is not `NULL`, it will overwrite the percentiles.

**Important**: If you override the defaults this way, please make sure to specify this in the plot annotations, as you are essentially redefining a box plot, and the reader might misinterpret it.

### 2.3.1 Example for Customized Percentiles

In the example below, we do not show the whiskers by setting whisker percentiles to box percentiles.

```{r customized-percentiles, fig.alt="Box plots showing age distribution by country with customized percentiles. The whiskers are removed by setting both whisker percentiles equal to the box percentiles, showing only the interquartile range."}
# B No whiskers
plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age
  ),
  data = pkRatioData,
  metaData = metaData,
  percentiles = c(0.25, 0.25, 0.5, 0.75, 0.75)
)
```

### 2.3.2 Example for Customized Aggregation Function

To customize the aggregation, provide a function that has as input the vector to aggregate and as output a named list with entries `ymin`, `lower`, `middle`, `upper`, and `ymax`. In the example below, the function `myStatFun` is provided, which uses mean and standard deviation to aggregate:

```{r aggregation-functions, fig.alt="Box plots showing age distribution by country using a custom aggregation function based on mean and standard deviation. The middle line represents the mean, box edges show mean ± SD, and whiskers extend to mean ± 1.96 SD."}
myStatFun <- function(y) {
  r <- list(
    ymin = mean(y) - 1.96 * stats::sd(y),
    lower = mean(y) - stats::sd(y),
    middle = mean(y),
    upper = mean(y) + stats::sd(y),
    ymax = mean(y) + 1.96 * stats::sd(y)
  )
  return(r)
}

plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age
  ),
  data = pkRatioData,
  metaData = metaData,
  statFun = myStatFun
) + labs(caption = "Mean for the middle line, mean +/- standard deviation for the box edges, and mean +/- 1.96 standard deviation for the whiskers.")
```

## 2.3.3 Tables Corresponding to Plot

As there already exist many possibilities to aggregate data in R, this package does not provide an extra function to return the aggregated data in a tabular form.

However, the `plotBoxWhisker()` adds the aggregation function to the plot object. So the aggregation can be easily done, e.g., with the `{data.table}` package using the same function as used for the plot.

```{r aggregationTable, fig.alt="Box plots showing age distribution by country and sex. This plot demonstrates how the aggregation function can be extracted from the plot object for creating corresponding data tables."}
# Generate plot object
plotObject <- plotBoxWhisker(
  mapping = aes(
    x = Country,
    y = Age,
    fill = Sex
  ),
  data = pkRatioData,
  metaData = metaData
)

plot(plotObject)

# Convert data to data.table and use statFun saved in plotObject for aggregation
# Make sure to add all relevant aesthetics to "by" columns
dt <- plotObject$data |>
  data.table::setDT() |>
  (\(x) x[, as.list(plotObject$statFun(Age)), by = c("Country", "Sex")])()

knitr::kable(dt)
```

## 2.4 Show Outliers

### 2.4.1 Outliers with Default Settings

Outliers are displayed if the variable `outliers` is TRUE. Default outliers are flagged when outside the range from the "25th" percentile - 1.5 x IQR to the "75th" percentile + 1.5 x IQR.

```{r example-outliers, fig.alt="Box plots showing ratio distribution by sex and country with outliers displayed. Outliers are shown as individual points outside the whiskers, which indicate the 90% range (5th-95th percentile). Default outlier detection uses 1.5 × IQR beyond the quartiles."}
plotBoxWhisker(
  mapping = aes(
    x = Sex,
    y = Ratio,
    fill = Country
  ),
  data = pkRatioData,
  metaData = metaData,
  outliers = TRUE
) + labs(tag = "A", caption = "Default settings: Whiskers indicate 90% range (5th - 95th percentile) and outliers indicate all measurements outside 25th percentiles - 1.5 x IQR to 75th percentiles + 1.5 x IQR.")
```

### 2.4.2 Outliers with Customized Settings

In the following example, the aggregation is customized to set the whiskers to 10% and 90% via the input variable `percentiles`, and the outlier range is customized to show all points outside the whiskers. For that, we define a function `myStatFunOutlier`, which has as input the vector of the values and as output the vector of values outside the outlier range.

```{r example-outliers-custom, fig.alt="Box plots showing ratio distribution by sex and country with custom outlier detection. The whiskers indicate the 80% range (10th-90th percentile), and outliers show all measurements outside the whisker range."}
myStatFunOutlier <- function(y) {
  q <- stats::quantile(y, probs = c(0.1, 0.9), names = FALSE, na.rm = TRUE)
  yOutsideRange <- subset(y, y < q[1] | y > q[2])

  if (length(yOutsideRange) < 1) {
    return(as.double(NA))
  } else {
    return(yOutsideRange)
  }
}

plotBoxWhisker(
  data = pkRatioData,
  metaData = metaData,
  mapping = aes(
    x = Sex,
    y = Ratio,
    fill = Country
  ),
  outliers = TRUE,
  percentiles = c(0.1, 0.25, 0.5, 0.75, 0.9),
  statFunOutlier = myStatFunOutlier
) + labs(caption = "Whiskers indicate 80% range (10th - 90th percentile) and outliers indicate all measurements outside whiskers.")
```

### 2.4.3 Retrieve Outliers as Table

The `plotBoxWhisker()` adds the outlier function to the plot object. So outliers can be easily retrieved, e.g., with the `{data.table}` package.

```{r outlierTable}
# Generate plotObject
plotObject <- plotBoxWhisker(
  data = pkRatioData,
  metaData = metaData,
  mapping = aes(
    x = Sex,
    y = Ratio,
    fill = Country
  ),
  outliers = TRUE,
  percentiles = c(0.1, 0.25, 0.5, 0.75, 0.9),
  statFunOutlier = myStatFunOutlier
)

# Convert data to data.table and use statFun saved in plotObject for aggregation
# Make sure to add all relevant aesthetics to "by" columns
dt <- plotObject$data |>
  data.table::setDT() |>
  (\(x) x[, .(outliers = plotObject$statFunOutlier(Ratio)), by = c("Country", "Sex")])()

knitr::kable(dt)
```

# 3. Plot Configuration of Boxplots

Below you see an example of the same plot:

- A default layout
- B customized layout
- The variable `geomBoxplotAttributes` now has an entry for color 'orange'.
- The variable `geomPointAttributes` now has an entry for size = 4, color 'orange', and shape = 'diamond'.
- Fill is added to the mapping, so the outlier points are also colored according to 'Sex'.
- With `scale_fill_manual`, the colors for 'Sex' are defined.

```{r boxplot-update-configuration, fig.alt="Two side-by-side box plots (A and B) showing age distribution by country and sex. Plot A shows default styling, while Plot B demonstrates custom styling with orange colors for boxes and diamond-shaped outlier points, plus custom fill colors for sex groups."}
# A default layout
plotBoxWhisker(
  data = pkRatioData,
  metaData = metaData,
  mapping = aes(
    x = Country,
    y = Age,
    fill = Sex
  ),
  outliers = TRUE,
  percentiles = c(0.1, 0.25, 0.5, 0.75, 0.9),
  statFunOutlier = myStatFunOutlier
) + labs(tag = "A")

# B customized layout
plotBoxWhisker(
  data = pkRatioData,
  metaData = metaData,
  mapping = aes(
    x = Country,
    y = Age,
    fill = Sex
  ),
  outliers = TRUE,
  percentiles = c(0.1, 0.25, 0.5, 0.75, 0.9),
  statFunOutlier = myStatFunOutlier,
  geomBoxplotAttributes = list(color = "orange", position = position_dodge(width = 1)),
  geomPointAttributes = list(position = position_dodge(width = 1), size = 4, color = "orange", shape = "diamond")
) + scale_fill_manual(values = c(Female = "pink", Male = "dodgerblue")) + labs(tag = "B")
```

```{r cleanup, echo=FALSE}
resetDefaults(oldDefaults)
```