---
title: "Visualizing Two Numeric Variables"
output:
  html_document:
    toc: yes
    code_folding: show
    code_download: true
---

<link rel="stylesheet" href="stat4580.css" type="text/css" />
<style type="text/css"> .remark-code { font-size: 85%; } </style>
<!-- title based on Wilke's chapter -->

```{r setup, include = FALSE, message = FALSE}
source(here::here("setup.R"))
knitr::opts_chunk$set(collapse = TRUE, message = FALSE,
                      fig.height = 5, fig.width = 6, fig.align = "center")

set.seed(12345)
library(dplyr)
library(ggplot2)
library(lattice)
library(gridExtra)
library(patchwork)
source(here::here("datasets.R"))
data(father.son, package = "UsingR")
```


## Slope Graphs

The most used graph for visualizing the relationship between two
numeric variables is the _scatter plot_.

But there is one alternative that can be useful and is increasingly
popular: the _slope chart_ or _slope graph_.


### Tufte's Slope Graph

Two articles on slope graphs with examples:

* https://charliepark.org/slopegraphs/
* https://visualisingdata.com/2013/12/in-praise-of-slopegraphs/

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics(IMG("tufteslope.gif"))
```
<!-- http://charliepark.org/images/slopegraphs/slopegraph.gif -->

Tufte showed this example in _The Visual Display of Quantitative
Information_.

Some features of the data that are easy to see:

* order of the countries within each year;

* how each country's values changed;

* how the rates of change compare;

* the country (Britain) that does not fit the general pattern.

The chart uses _no non-data ink_. (Tufte thinks this is important;
others disagree).

The chart in this form is well suited for small data sets or summaries
with modest numbers of categories.

Scalability in this full form is limited, but better if labels and
values are dropped.

The idea can be extended to multiple periods, though two periods or
levels is most common when labeling is used.

Without labeling this becomes a _parallel coordinates plot_.

A multi-period example:
```{r, echo = FALSE, out.width = "60%"}
knitr::include_graphics(IMG("cancer_survival_nash.gif"))
```
<!-- http://charliepark.org/images/slopegraphs/cancer_survival_nash.gif -->


### Barley Mean Yields

A slope graph for the average yields at each experiment station for
the two years 1931 and 1932:

```{r, message = FALSE, class.source = "fold-hide"}
theme_set(theme_minimal() + theme(text = element_text(size = 16)))
library(ggrepel)
library(forcats)
data(barley, package = "lattice")
barley_site_year <- group_by(barley, site, year) |>
    summarize(avg_yield = mean(yield)) |>
    mutate(year = fct_rev(year))
barley_site_year_1932 <- filter(barley_site_year, year == "1932")
ggplot(barley_site_year, aes(x = year, y = avg_yield, group = site)) +
    geom_line() +
    geom_text_repel(aes(label = site),
                    data = barley_site_year_1932,
                    hjust = "left", direction = "y") +
    scale_x_discrete(expand = expansion(mult = c(0.1, .3)),
                     position = "top") +
    labs(x = NULL, y = "Average Yield") +
    theme(axis.line.y.left = element_line(),
          axis.ticks.y.left = element_line())
```

The anomalous result for Morris pops out very clearly.

This graph departs from the classic Tufte style:

* it uses an axis instead of showing the numbers;

* only shows labels on one side.

This is similar to the style used
[here](https://clauswilke.com/dataviz/visualizing-associations.html#associations-paired-data).


#### Creating the Graph

The first step is to compute the averages:

```{r, eval = TRUE}
barley_site_year <- group_by(barley, site, year) |>
    summarize(avg_yield = mean(yield))
head(barley_site_year, 2)
```

The `year` variable is a factor with the levels in the wrong order, so
we need to fix that:

```{r, eval = TRUE}
levels(barley_site_year$year)
barley_site_year <- mutate(barley_site_year, year = fct_rev(year))
levels(barley_site_year$year)
```

Set the default theme to `theme_minimal` with larger text:

```{r, eval = TRUE}
theme_set(
    theme_minimal() +
    theme(text = element_text(size = 16)))
```

The core of a slope graph is produced by

```{r slope-core, eval = FALSE}
p <- ggplot(barley_site_year,
            aes(x = year,
                y = avg_yield,
                group = site)) +
    geom_line()
p
```

```{r slope-core, echo = FALSE}
```

Adding the labels on the 1932 side can be done as

```{r slope-core-label, eval = FALSE}
barley_site_year_1932 <-
    filter(barley_site_year, year == "1932")
p + geom_text(aes(label = site),
              data = barley_site_year_1932,
              hjust = "left")
```

```{r slope-core-label, echo = FALSE}
```

The label positions could use further adjusting.

Using `geom_text_repel` from the `ggrepel` package handles this well:

```{r slope-core-repel, eval = FALSE}
library(ggrepel)
p <- p +
    geom_text_repel(
        aes(label = site),
        data = barley_site_year_1932,
        hjust = "left",
        direction = "y")
p
```

```{r slope-core-repel, echo = FALSE}
```

Adjust the `x` scale:

```{r slope-core-x, eval = FALSE}
p <- p +
    scale_x_discrete(
        expand = expansion(mult = c(0.1, .3)),
        position = "top")
p
```

```{r slope-core-x, echo = FALSE}
```

Final theme adjustments:

```{r scale-core-theme, eval = FALSE}
p + labs(x = NULL,
         y = "Average Yield") +
    theme(axis.line.y.left = element_line(),
          axis.ticks.y.left = element_line())
```

```{r scale-core-theme, echo = FALSE}
```


### Father-Son Heights

The `father.son` data set has `r nrow(father.son)` observations, which
is too large for the labeled slope graph, but the basic representation
is useful:

```{r, fig.height = 4.5, class.source = "fold-hide"}
library(tidyr)
fs <-
    mutate(father.son,
           id = seq_len(nrow(father.son))) |>
    pivot_longer(1 : 2,
                 names_to = "which",
                 values_to = "height")
ggplot(fs, aes(x = which, y = height)) +
    geom_line(aes(group = id), alpha = 0.1) +
    scale_x_discrete(
        expand = expansion(mult = c(.1, .1)),
        labels = c("Father", "Son"),
        position = "top") +
    labs(x = NULL, y = "Height (Inches)")
```

This very clearly shows the famous _regression to the mean_ effect:

* taller parents tend to be taller than their children;
* shorter parents tend to be shorter than their children.

Conversely,

* taller children tend to be taller than their parents;
* shorter children tend to be shorter than their parents.

<!--
fs <- mutate(fs, TorSF = abs(height - mean(height)) > 1.5 * sd(height) &
                             which == "fheight")
ggplot(fs) +
geom_line(aes(x = which, y = height, group = id, color = TorSF), alpha = 0.1) +
scale_x_discrete(expand = c(.1, 0)) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black"))
-->


#### Creating the Graph

To make creating the graph easier we can convert the data frame into
a longer form with variables

* `height`, the height measurement
* `which`, `fheight` or `sheight`
* `id`, identifying the pair:

Add the `id` variable:
    
```{r, eval = TRUE}
fs <- mutate(father.son,
             id = seq_len(nrow(father.son)))
head(fs)
```

Pivot to the longer format:

```{r, eval = TRUE}
fs <- pivot_longer(fs, 1 : 2,
                   names_to = "which",
                   values_to = "height")
head(fs)
```

The basic plot is quite simple:

```{r fs-basic, eval = FALSE}
ggplot(fs, aes(x = which,
               y = height,
               group = id)) +
    geom_line()
```

```{r fs-basic, echo = FALSE}
```

With an `alpha` adjustment to reduce over-plotting:

```{r fs-alpha, eval = FALSE}
ggplot(fs, aes(x = which,
               y = height,
               group = id)) +
    geom_line(alpha = 0.1)
```

```{r fs-alpha, echo = FALSE}
```

With scale and label adjustments:

```{r fs-theme, eval = FALSE}
ggplot(fs, aes(x = which,
                    y = height,
                    group = id)) +
    geom_line(alpha = 0.1) +
    scale_x_discrete(
        expand = expansion(mult = c(.1, .1)),
        labels = c("Father", "Son"),
        position = "top") +
    labs(x = NULL, y = "Height (Inches)")
```

```{r fs-theme, echo = FALSE}
```


## Scatter Plots

```{r, include = FALSE}
theme_set(theme_minimal() +
          theme(text = element_text(size = 16)) +
          theme(panel.background = element_rect(color = "grey", fill = NA)))
```
The most common way to show the relationship between two variables is a
_scatter plot_.

A scatter plot of two variables maps the values of one variable
to the vertical axis and the other to the horizontal axis of a
cartesian coordinate system and places a mark for each observation
at the resulting point.

Some conventions:

* Plot `A` versus/against `B` means `A` is mapped to the vertical, or
  $y$, axis, and `B` to the horizontal, or $x$ axis.

* If we can think of variation in `A` as being partly explained by `B`
  then we usually plot `A` against `B`.

* If we can think of `B` as helping to predict `A`, then we usually
  plot `A` against `B`.


### Barley Yields

For a scatter plot of mean yield in 1932 against mean yield in 1931
for the different sites it is useful to have a data frame containing
variables for each year.

This requires converting the data frame to a wider format.

```{r}
wide_barley_site_year <- pivot_wider(barley_site_year,
                                     names_from = "year",
                                     names_prefix = "avg_yield_",
                                     values_from = "avg_yield")
head(wide_barley_site_year)
```

The basic scatter plot of `y = avg_yield_1932` against `x =
avg_yield_year1931`:

```{r, class.source = "fold-hide"}
p <- ggplot(wide_barley_site_year,
            aes(x = avg_yield_1931,
                y = avg_yield_1932)) +
    geom_point()
p
```

Adding labels using `geom_text_repel` identifies the Morris site:

```{r, class.source = "fold-hide"}
p <- p + geom_text_repel(aes(label = site))
p
```

To recognize the reversal for Morris we can add the 45 degree line:

```{r, class.source = "fold-hide"}
p + geom_abline(aes(intercept = 0, slope = 1),
                linetype = 2)
```

A 45 degree line also helps when viewing the full data:

```{r, class.source = "fold-hide"}
bw <- pivot_wider(barley,
                  names_from = "year",
                  names_prefix = "yield_",
                  values_from = "yield")
ggplot(bw, aes(x = yield_1931,
               y = yield_1932)) +
    geom_point() +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) +
    geom_point(data = filter(bw, site == "Morris"),
               color = "red") +
    ggtitle("Barley Yields",
            "Values for Morris are shown in red.") +
    labs(x = "1931", y = "1932")
```

A recent [blog
post](https://eagereyes.org/blog/2020/in-praise-of-the-diagonal-reference-line)
discusses the value of reference lines as plot annotations.

If the primary goal is to show the change from one year to the next
then a _mean-difference plot_ is a good choice:

```{r, class.source = "fold-hide"}
ggplot(wide_barley_site_year,
       aes(x = (avg_yield_1932 + avg_yield_1931) / 2,
           y = avg_yield_1932 - avg_yield_1931)) +
    geom_point() +
    geom_text_repel(aes(label = site)) +
    geom_abline(aes(intercept = 0, slope = 0), linetype = 2)
```

For the full data:

```{r, class.source = "fold-hide"}
ggplot(bw,
       aes(x = (yield_1932 + yield_1931) / 2,
           y = yield_1932 - yield_1931)) +
    geom_point() +
    geom_point(data = filter(bw, site == "Morris"),
               color = "red") +
    geom_abline(aes(intercept = 0,
                    slope = 0),
                linetype = 2)
```

The comparison of changes is now an aligned axis comparison.

Mean-difference plots are also known as

* Tukey mean-difference plots;
* MA-plots;
* [Bland-Altman plots](https://en.wikipedia.org/wiki/Bland%E2%80%93Altman_plot).

Plotting the difference against the `x` variable is also often useful:

```{r, class.source = "fold-hide"}
ggplot(bw,
       aes(x = yield_1931,
           y = yield_1932 - yield_1931)) +
    geom_point() +
    geom_point(
        data = filter(bw, site == "Morris"),
        color = "red") +
    geom_abline(aes(intercept = 0,
                    slope = 0),
                linetype = 2) +
    ggtitle(
        "Barley Yield Differences",
        "Values for Morris are shown in red.") +
    labs(x = "Yield in 1931",
         y = "Difference in Yield for 1932")
```


### Father and Son Heights

The basic scatter plot:

```{r, class.source = "fold-hide"}
p0 <- ggplot(father.son,
             aes(x = fheight,
                 y = sheight))
p1 <- p0 + geom_point()
p1
```

Adding a line with slope one helps identify the regression to the mean
phenomenon:

```{r, class.source = "fold-hide"}
p2 <- p1 +
    geom_abline(
        aes(intercept =
                mean(sheight) - mean(fheight),
            slope = 1),
        color = "red",
        linewidth = 1.5)
p2
```

Adding a regression line helps further:

```{r, class.source = "fold-hide"}
p2 + geom_smooth(method = "lm")
```

But for showing the regression effect it is hard to beat the scatter
plot of `sheight - fheight` against `fheight`:

```{r, class.source = "fold-hide"}
ggplot(father.son) +
    geom_point(aes(x = fheight,
                   y = sheight - fheight)) +
    geom_hline(aes(yintercept = 0),
               linetype = 2)
```


### Old Faithful Eruptions

A scatter plot of the waiting times until the next eruption against
the duration of the current eruption for the `faithful` data set shows
the two clusters corresponding to the short and long eruptions:

```{r, class.source = "fold-hide"}
ggplot(faithful) +
    geom_point(aes(x = eruptions,
                   y = waiting))
```

For the `geyser` data set from the `MASS` package a plot of the two
variables shows a different pattern:

```{r, class.source = "fold-hide"}
data(geyser, package = "MASS")
ggplot(geyser) +
    geom_point(aes(x = duration,
                   y = waiting))
```

The reason for the difference is that in the `geyser` data set the
waiting time reflects the time since the _previous_ eruption, not the
time until the _next_ one.

For this ordering it is more natural to plot `duration` against `waiting`:

```{r, class.source = "fold-hide"}
ggplot(geyser) +
    geom_point(aes(x = waiting,
                   y = duration))
```

How well does the waiting time predict whether the duration will be
longer or shorter?

The question the park service is more interested in: How well
does duration predict waiting time until the next eruption?

We can adjust these data to pair durations with waiting times until
the next eruption using the `lag` function from `dplyr`.

This produces the same basic pattern as for the `faithful` data set:

```{r, class.source = "fold-hide"}
ggplot(geyser) +
    geom_point(aes(x = lag(duration),
                   y = waiting),
               na.rm = TRUE)
```


## Reading

Chapter [_Visualizing associations among two or more quantitative
  variables_](https://clauswilke.com/dataviz/visualizing-associations.html)
  in [_Fundamentals of Data
  Visualization_](https://clauswilke.com/dataviz/).
