--- title: "Visualizing Two Numeric Variables" output: html_document: toc: yes code_folding: show code_download: true --- ```{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://www.visualisingdata.com/2013/12/in-praise-of-slopegraphs/ ```{r, echo = FALSE, out.width = "80%"} knitr::include_graphics(IMG("tufteslope.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")) ``` ### 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. #### 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_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/).