---
title: "Visualising results from rsimsum"
author: "Alessandro Gasparini"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Visualising results from rsimsum}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options:
  chunk_output_type: console
---

```{r setup, include = FALSE}
options(width = 150)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center", fig.height = 5, fig.width = 7,
  out.width = "90%"
)
```

This vignette requires the following R packages:

```{r packages}
library(rsimsum)
library(ggplot2)
```

# Data

We use data from a simulation study on misspecification of the baseline hazard in survival models. This dataset is included in `rsimsum` and can be loaded with:

```{r data-1}
data("relhaz", package = "rsimsum")
```

Inspecting the structure of the dataset and the first 15 rows of data:

```{r inspect-1}
str(relhaz)
head(relhaz, n = 15)
```

# Summarise results

We use the `simsum` function to summarise results:

```{r summarise-1}
s1 <- simsum(
  data = relhaz, estvarname = "theta", se = "se", true = -0.50,
  methodvar = "model", by = c("n", "baseline"), x = TRUE
)
s1
```

We call `simsum` with `x = TRUE` as that is required for some types of plots (e.g. zip plots, scatter plots, etc.).

```{r summary-summarise-1}
summary(s1)
```

`rsimsum` implements the `autoplot` method for objects of classes: `simsum`, `summary.simsum`, `multisimsum`, `summary.multisimsum`.

See `?ggplot2::autoplot()` for details on the S3 generic function.

# Scatter plots

Scatter plots allow to assess serial trends in estimates and standard errors.

For instance, if we want to compare the point estimates from different methods (across data-generating mechanisms):

```{r plot-est}
autoplot(s1, type = "est")
```

Analogously, if we want to compare standard errors:

```{r plot-se}
autoplot(s1, type = "se")
```

These two plot types allow comparing estimates (and standard errors) obtained from different methods with ease; in ideal settings, the points of the scatterplots should lay on the diagonal (dashed line).
An estimated regression line (of `X` vs `Y`, blue line) is superimposed by default to ease the comparison even more.

In addition to plots comparing estimates and standard errors, Bland-Altman-type plots are supported as well:

```{r plot-est-ba}
autoplot(s1, type = "est_ba")
```

Bland-Altman plots compare the difference between estimates from two competing methods (on the y-axis) with the mean of the estimates from the two methods (x-axis).
In the ideal scenario, there should be no trend and all the points from the scatter plot should lay around the horizontal dashed line.
To ease comparison, a regression line is included here as well.
Bland-Altman plots for standard errors could be obtained by setting the argument `type = "se_ba"`.

## Ridgeline plots

Another way of visually comparing estimates from different methods is given by ridgeline plots.
According to the documentation of the `ggridges` package ([source](https://cran.r-project.org/package=ggridges/vignettes/introduction.html)),

> Ridgeline plots are partially overlapping line plots that create the impression of a mountain range.
> They can be quite useful for visualizing changes in distributions over time or space.

In the settings of simulation studies, we aim to visualise changes in distribution over data-generating mechanisms.

For instance, say we want to compare estimates across data-generating mechanisms and methods:

```{r ridge-est-1}
autoplot(s1, type = "est_ridge")
```

This allows us to see how the estimates from the `Exp` method are different from the other methods in two of the four data-generating mechanisms: `50, Weibull` and `250, Weibull`.

To obtain a similar plot for standard errors, call the `autoplot` method with `type = "se_ridge` instead.

# Lolly plots

Lolly plots are used to present estimates for a given summary statistic with confidence intervals based on Monte Carlo standard errors (if calling the `autoplot` method on `summary` objects).
They allow to easily compare methods.

Say we are interested in bias:

```{r lolly-bias-1}
autoplot(summary(s1), type = "lolly", stats = "bias")
```

It is straightforward to identify the exponential model as yielding biased results when the true baseline hazard is Weibull, irrespectively of the sample size.

On a relative scale, we can plot _relative bias_ as well:

```{r lolly-bias-1.1}
autoplot(summary(s1), type = "lolly", stats = "rbias")
```

If confidence intervals based on Monte Carlo errors are not required, it is sufficient to call the `autoplot` method on the `simsum` object:

```{r lolly-bias-2}
autoplot(s1, type = "lolly", stats = "bias")
```

Analogously, for coverage:

```{r lolly-coverage-1}
autoplot(summary(s1), type = "lolly", stats = "cover")
```

# Forest plots

Forest plots could be an alternative to lolly plots, with similar interpretation:

```{r forest-bias-1}
autoplot(s1, type = "forest", stats = "bias")
```

```{r forest-bias-2}
autoplot(summary(s1), type = "forest", stats = "bias")
```

# Zipper plots

Zipper plots (or zip plots), introduced in Morris _et al_. (2019), help understand coverage by visualising the confidence intervals directly. For each data-generating mechanism and method, the confidence intervals are centile-ranked according to their significance against the null hypothesis \(H_0: \theta\) = \(\theta_{\text{true}}\), assessed via a Wald-type test. This ranking is used for the vertical axis and is plotted against the intervals themselves.

When a method has 95% coverage, the colour of the intervals switches at 95 on the vertical axis. Finally, the horizontal lines represent confidence intervals for the estimated coverage based on Monte Carlo standard errors.

```{r zipper-1}
autoplot(s1, type = "zip")
```

The zipper plot for the exponential model with n = 50 and a true Weibull baseline hazard shows that coverage is approximately 95%; however, there are more intervals to the right of \(\theta\) = -0.50 than to the left: this indicates that the model standard errors must be overestimating the empirical standard error, because coverage is appropriate despite bias.

It is also possible to _zoom_ on the top x% of the zip plot to increase readability, e.g. on the top 30%:

```{r zipper-zoom}
autoplot(s1, type = "zip", zoom = 0.3)
```

# Heat plots

Heat plots are a new visualisation that we suggest and include here for the first time. With heat plots, we produce a heat-map-like plot where the filling of each tile represents a given summary statistic, say bias:

```{r heat-bias-1}
autoplot(s1, type = "heat", stats = "bias")
```

This visualisation type automatically puts all data-generating mechanisms on the y-axis; by default, the methods are included on the `x` axis. Therefore, this plot is most useful when a simulation study includes different methods to be compared and many data-generating mechanisms. Using a heat plot, it is immediate to identify visually which method performs better and under which data-generating mechanisms.

By default, heat plots use the default `ggplot` scale for the filling aesthetic.
It is recommended to use a different colour palette with better characteristics, e.g. the viridis colour palette from `matplotlib`; see the next section for details on how to do this, and [here](https://cran.r-project.org/package=viridis/vignettes/intro-to-viridis.html) for details on the viridis colour palette.

# Contour plots and hexbin plots

Individual point estimates and standard errors could also be plotted using contour plots or hexbin plots.

Contour plots represent a 3-dimensional surface by plotting constant z slices (called contours) on a 2-dimensional format.
That is, given a value for z, lines are drawn for connecting the (x, y) coordinates where the value of z is (relatively) homogenous.

Hexbin plots are useful to represent the relationship of 2 numerical variables when you have a lot of data points: instead of overlapping, the plotting window is split in several hexbins, and the number of points per hexbin is counted.
The colour filling denotes then the number of points.

Both plots provide an alternative to scatter plots when there is a large number of data points that overlap.
Contour plots and hexbin plots can be easily obtained using the `autoplot` method once again, and using the argument `type = "est_density"`, `type = "se_density"`, `type = "est_hex"`, or `type = "se_hex"`.
For instance, focussing on point estimates:

```{r plot-est-density}
autoplot(s1, type = "est_density")
```

```{r plot-est-hex}
autoplot(s1, type = "est_hex")
```

Of course, analogous plots could be obtained for standard errors.

# Custom plotting

All plots produced by `rsimsum` are meant to be quick explorations of results from Monte Carlo simulation studies: they are not meant to be final manuscript-like-quality plots (although they can be useful as a starting point).

Generally, the output of all types of `autoplot` calls are `ggplot` objects; hence, it should be generally straightforward to customise all plots.
For instance, say we want to add a custom theme:

```{r add-theme}
autoplot(summary(s1), type = "lolly", stats = "bias") +
  ggplot2::theme_bw()
```

Compare to the default plot:

```{r add-theme-default}
autoplot(summary(s1), type = "lolly", stats = "bias")
```

We also mentioned before that the colour palette of heat plots should be customised. Say we want to use the default viridis colour palette:

```{r heat-bias-viridis}
autoplot(s1, type = "heat", stats = "bias") +
  ggplot2::scale_fill_viridis_c()
```

Analogously, say we want to customise the colour palette of ridgeline plots, again using the viridis colour palette:

```{r colour-palette}
autoplot(s1, type = "est_ridge") +
  scale_fill_viridis_d() +
  scale_colour_viridis_d()
```