This glossary collects the key concepts, methods and people in this lesson. It can be used as a reference while working through the material or as a review before assessments. Typing in the search box filters the entries.

Lesson 6 · Section 1

Look Before You Model

Plots show the shape of the data, and the shape decides which summaries and models make sense.

Exploratory data analysis

Two kinds of analysis

Exploratory

The analyst asks what the data look like.

Plots and simple summaries show shapes, unusual values and relationships, and they raise new questions.

Confirmatory

The analyst tests a question that was set in advance.

Models give estimates, confidence intervals and p-values for that question.

Tukey (1977) argued that exploration should come first.

The running case

The working group's three questions

• The group asks how lonely people in the survey are.
• The group asks whether loneliness differs by age and gender.
• The group asks how loneliness goes together with social support and mental health.

4,045

respondents in the CSCS 2021 wave

3,083

respondents with complete data on the briefing variables

Anscombe's quartet

Four datasets with the same summary statistics

Statistic	Set 1	Set 2	Set 3	Set 4
Mean of x and of y	9, 7.5	9, 7.5	9, 7.5	9, 7.5
SD of x and of y	3.32, 2.03	3.32, 2.03	3.32, 2.03	3.32, 2.03
Correlation	0.816	0.816	0.816	0.817
Fitted line	y = 3.00 + 0.500 x in every set

Anscombe's quartet

The plots tell four different stories

What to look for

Features that plots reveal

One variable

The analyst checks the shape (symmetry, skew and peaks), unusual and impossible values, heaping on round or extreme values, and missing values.

Two or more variables

The analyst checks the direction of a relationship, its form (straight or curved), its strength, and whether it differs between groups.

Choosing a display

Match the display to the question and the variable types

Recall: Lesson 2, Section 3

A histogram shows one numeric variable, a bar chart the counts of a categorical variable, boxplots a numeric variable across groups, and a scatterplot two numeric variables.

New: many variables

A scatterplot matrix or a correlation heatmap gives an overview, and facets repeat a display for each group.

Check: loneliness across four age groups calls for boxplots side by side, with the points shown.

Poor displays

Four displays that conceal the data

Carry forward

What to take into the next section

• Identical summary statistics can describe very different data, so the analyst plots first.
• The display is chosen from the question and the types of the variables involved.
• Some displays conceal the data, and each has a better alternative.

Introduction and Overview

Lessons 4 and 5 fitted regression models to outcomes of many kinds. Each of those models rested on assumptions about the shape of the data: a straight-line relationship, residuals with an even spread, counts that are not more variable than the model allows, or observations that are independent. This lesson returns to the step that should come before any model is fitted, which is looking carefully at the data. Exploratory data analysis uses plots and simple summaries to find the shape of each variable, unusual or impossible values, and the relationships between variables. It tells the analyst whether the planned model suits the data, and it often suggests questions that the model should answer.

Box 6.1 introduces the running case for the lesson, a one-page visual briefing on social connection that an analyst prepares for a regional working group on loneliness. Each section of the lesson returns to this briefing, so the box sets out the group's three questions and the survey data from which the plots are drawn.

The rest of this section explains what exploratory data analysis is, shows with Anscombe's quartet why summary statistics alone can mislead, matches displays to the questions they answer, and examines four common poor displays.

What Exploratory Data Analysis Is

The statistician John Tukey set out this kind of work in his book Exploratory Data Analysis (Tukey, 1977). He distinguished two activities. Confirmatory analysis tests a question that was set before the data were examined, and it reports estimates, confidence intervals and p-values. Exploratory analysis asks what the data look like, using plots and simple summaries, and it is open to finding things that nobody planned to look for. Tukey compared exploration to detective work: the analyst gathers clues about the data before deciding which formal test or model to apply.

Exploration serves three purposes in an epidemiological analysis. It checks data quality, by revealing values that are impossible, implausible or suspiciously common. It describes distributions, which decides whether a mean or a median is the better summary and whether a variable needs to be transformed or grouped. It shows relationships, including whether they are straight or curved and whether they differ between groups, which decides how variables should enter a model. Lesson 2 introduced numeric summaries, data cleaning and the basic plots for one and two variables (histograms, boxplots, bar charts and scatterplots). This lesson develops those plots into an exploratory workflow and adds displays for many variables.

Table 6.1 sets out five questions that the analyst asks during exploration, what a plot can show in answer to each, and how the answer changes the later analysis.

Table 6.1. Questions asked during exploration, what a plot can show, and what the answer changes in the analysis.

Each question in Table 6.1 is answered by looking at a plot. The next part shows, with four small datasets, why summary statistics cannot take the place of a plot.

Anscombe's Quartet: Why Summary Statistics Can Mislead

In 1973 the statistician Francis Anscombe published four small datasets to argue that graphs are an essential part of statistical analysis and that a computer should produce graphs as well as calculations (Anscombe, 1973). Each dataset has eleven pairs of values, labelled x and y, and the four datasets are supplied with R as a data frame called anscombe. Their summary statistics are almost identical, as Table 6.2 shows.

Table 6.2. Summary statistics for the four datasets in Anscombe's quartet (Anscombe, 1973).

A report that gave only these numbers would describe the four datasets as the same. The plots show that they are very different. Figure 6.1 draws each of the four sets as a scatterplot with its fitted line.

The accordion below takes the four panels of Figure 6.1 in turn and states, for each set, what the plot shows and whether the fitted line is a fair summary.

The same lesson has been repeated with more recent examples. Starting from a dataset drawn in the outline of a dinosaur by Alberto Cairo, Matejka and Fitzmaurice (2017) generated the "Datasaurus Dozen", a set of further datasets that form stars, circles and lines, each with the same means, standard deviations and correlation as the dinosaur to two decimal places. The point is the same as Anscombe's: a summary statistic describes one feature of the data, and only a plot shows whether that feature is the one that matters.

Activity 6.1 reproduces Table 6.2 and Figure 6.1 in R from the anscombe data frame, and its first question asks what an analyst who saw only the summary statistics would conclude.

anscombe                          # four small datasets that come with R
round(colMeans(anscombe), 2)      # the mean of each column
round(sapply(anscombe, sd), 2)    # the standard deviation of each column

round(cor(anscombe$x1, anscombe$y1), 3)   # correlation in set 1
round(cor(anscombe$x2, anscombe$y2), 3)   # set 2
round(cor(anscombe$x3, anscombe$y3), 3)   # set 3
round(cor(anscombe$x4, anscombe$y4), 3)   # set 4
coef(lm(y1 ~ x1, data = anscombe))        # intercept and slope of the line, set 1
coef(lm(y4 ~ x4, data = anscombe))        # the same for set 4

par(mfrow = c(2, 2))    # four plots in one window: two rows, two columns
plot(anscombe$x1, anscombe$y1, main = "Set 1", xlab = "x", ylab = "y")
plot(anscombe$x2, anscombe$y2, main = "Set 2", xlab = "x", ylab = "y")
plot(anscombe$x3, anscombe$y3, main = "Set 3", xlab = "x", ylab = "y")
plot(anscombe$x4, anscombe$y4, main = "Set 4", xlab = "x", ylab = "y")
par(mfrow = c(1, 1))    # back to one plot per window

Anscombe's quartet shows that almost identical summary statistics can describe very different data, so the analyst needs a plot that suits each question. The next part sets out how the number and the type of the variables point to a display.

Matching the Display to the Question

Choosing a display begins with two questions about the request: how many variables it involves, and what type each variable is. The variable types are those from Lesson 2. A numeric variable records an amount and is continuous (any value in a range, such as age) or discrete (whole numbers, such as the number of people in a household). A categorical variable places each person into a group and is binary, ordinal (ordered groups, such as self-rated mental health from poor to excellent) or nominal (unordered groups, such as province). Scale scores such as the loneliness score, which takes the whole-number values 3 to 9, are numeric for plotting purposes, although their small number of distinct values affects how they look, as later sections show.

Box 6.2 recalls the guide to the basic plots from Lesson 2 and poses a retrieval question about comparing loneliness across age groups.

Figure 6.2 extends that guide into a decision tree, which starts from the number of variables in the question and then branches on their types.

A histogram divides the range of a numeric variable into intervals (bins) and draws a bar for the number of people in each one, so the bars touch and their order follows the number line. A bar chart draws one bar for each category of a categorical variable, with gaps between the bars, and the categories can be placed in any sensible order, such as the order of an ordinal scale or from most to least common. Students often confuse the two because both use bars. The test is the variable on the horizontal axis: a histogram has a numeric axis, and a bar chart has categories.

For two variables, the display follows the pair of types. A scatterplot places one point per person for two numeric variables. Boxplots side by side compare a numeric variable across the groups of a categorical variable, and adding the individual points shows how many people are in each group and how their values are spread. For two categorical variables, a bar chart of percentages (for example, the percentage of each gender in each mental health category) compares the groups fairly when the groups differ in size. For many variables, a scatterplot matrix or a correlation heatmap gives an overview of every pair, and facets repeat one display for each group in a grid of small panels, an idea that Tufte (1983) called small multiples.

Interactive 6.1 applies the decision tree in Figure 6.2 to six requests from the working group. Each request is answered by choosing a display, and the feedback under the options explains why each choice suits the request or fails it.

Choosing the right type of display is the first step. A display of the right type can still mislead when it is drawn poorly, and the next part examines four common examples.

A Gallery of Poor Displays

Research on graphical perception has found that people judge positions along a common scale more accurately than lengths, angles or areas (Cleveland & McGill, 1984). A good display therefore encodes the comparison that matters as a position or a length on a shared axis, and it shows the data themselves wherever possible. The four displays below break these principles in ways that are common in reports and presentations. Each tab shows the poor version beside a better one, drawn from the briefing file.

The four tabs hold Figure 6.3 (a truncated axis), Figure 6.4 (means without spread), Figure 6.5 (a pie chart with many slices) and Figure 6.6 (overplotting), and the text under each figure states what the poor version conceals.

Truncated axis Means without spread Pie with many slices Overplotting

What it conceals. A bar encodes a value by its length, so a bar chart whose axis starts above zero misrepresents every comparison. On the left, the bar for men looks about one tenth of the height of the bar for non-binary respondents, although the means are 5.36 and 5.93, a ratio of about nine to ten. The right panel starts the axis at zero, so the lengths are proportional to the values. A related distortion is the dual axis, in which two series share a plot with two different vertical scales; the apparent crossing points and relative heights depend entirely on how the two scales were chosen, and two separate panels are usually clearer.

Note on axis ranges. The zero rule applies to bars. A dot or line display may use an axis that covers the possible range of the scale (3 to 9 for this loneliness score), provided the range is labelled.

What it conceals. The bars show four means between 5.39 and 6.01 and nothing else. They hide the spread within each group (scores from 3 to 9 in every group), the shape of each distribution, and the number of people in each group, which ranges from 366 people aged 65 and older to 1,079 people aged 16 to 29. Weissgerber and colleagues (2015) reviewed 703 research articles in leading physiology journals, found that bar graphs were the most common way of presenting continuous data (85.6% of the articles included at least one), and recommended displays that show the distribution of the data, particularly dot plots of the individual values when samples are small. The right panel shows that the differences between age groups are small compared with the variation within each group.

What it conceals. A pie chart asks the reader to compare angles and areas, which people judge less accurately than positions along a common scale. With fourteen slices, several of them thin, it is hard to tell whether Manitoba or New Brunswick has more respondents, or to read the smaller territories at all. The sorted bar chart places every category on a common scale, so the order and the size of the differences are clear. A pie chart can be acceptable for two or three categories that make up a whole; for more categories, a bar chart is clearer.

What it conceals. The loneliness score takes only seven values, so 3,083 points fall on seven horizontal lines and lie on top of one another. The left panel cannot show where most people are, and a row with five people looks the same as a row with five hundred. The right panel adds a small random vertical shift to each point (jitter) and draws the points with transparency, so dense areas look darker. The fitted line makes the downward trend easy to see: people with more social support tend to report less loneliness.

Several of these displays compare groups, such as genders or age groups. Box 6.3 sets out when groups are better shown by colour within one panel and when they are better shown in separate panels.

This section has explained why the analyst looks at the data before fitting a model, how the question and the variable types point to a display, and what four common poor displays conceal. The knowledge check and the reflection that follow review these ideas, and Section 2 draws the displays for the briefing file in base R.

Lesson 6 · Section 2

Quick Looks in Base R

Base R draws a useful diagnostic plot with one line of code.

Base R plotting

Four familiar functions and one new one

Function	What it draws	Variables
hist(x)	Histogram	One numeric
barplot(table(x))	Bar chart of counts	One categorical
boxplot(y ~ group)	Boxplots side by side	Numeric across groups
plot(x, y)	Scatterplot	Two numeric
pairs(data)	Scatterplot matrix	Several numeric

Recall: the first four were used in the Lesson 2, Section 3 R activity, and pairs() is new.

Preparing the data

One briefing file for every plot

• The code loads the CSCS file and keeps the 2021 wave, which has 4,045 people.
• Short names replace the long survey names, for example loneliness and support.
• Bracket recoding sets "Presented but no response" to missing and builds four age groups.
• na.omit() keeps the 3,083 people with complete data on all nine briefing variables.

One variable

A histogram and a bar chart

Finding a problem

A spike at the largest value

Responding to the problem

Check, decide and report

Check

The codebook and the source questions show that 20 means 20 or more.

Decide

A count analysis sets 20 to missing, and a grouped variable absorbs the problem.

Report

The methods section states the problem, the decision and the number of people affected.

Two variables

Boxplots and a scatterplot

Layouts and saving

Many plots at once, and plots saved to files

Code	What it does
pairs(brief[, c(...)])	It draws a scatterplot for every pair of the chosen variables.
par(mfrow = c(1, 2))	It places the next plots in one row of two panels.
png("file.png", width, height, res)	It sends the next plots to an image file.
dev.off()	It closes the file, which saves it.

Carry forward

What to take into the next section

• Base R draws a quick diagnostic plot for each kind of question.
• Quick plots reveal problems such as spikes at a top value and heaping on round numbers.
• Jitter and transparency reveal patterns when many points share the same values.

Introduction and Overview

Section 1 explained why the analyst looks before modelling and how the question and the variable types point to a display. This section draws those displays for the survey data with base R, the plotting functions that are part of R itself. Base R plots need no extra packages and usually one line of code, which makes them the analyst's quick diagnostic tool: they are drawn, read and set aside, and only the useful ones are later rebuilt as polished figures (Section 3). The section also prepares the briefing file that every later plot uses, and it uses a quick plot to find a data problem that the cleaning lesson left behind.

The section begins with the labelling arguments that every base R plot accepts. It then prepares the briefing file, reads plots of one variable and of two variables, and ends with layouts of several plots and plots saved to files.

Base R as the Analyst's Quick Diagnostic Tool

Base R draws a quick diagnostic plot with one line of code. Every plotting function accepts the same labelling arguments: main for the title, xlab and ylab for the axis labels. Every plot in the briefing should carry these labels from the start, because a plot without labels is easily misread, even by the person who drew it.

Box 6.4 recalls the four standard plotting functions from the Lesson 2 activity and names the briefing variables to which this section applies them.

This section adds three tools that Lessons 1 and 2 did not use. The function pairs(data) draws a scatterplot matrix, one small scatterplot for every pair of several numeric variables, such as loneliness, support, depression and age. The setting par(mfrow = ...) places several plots in one window, and png() sends plots to a file; both are described at the end of this section.

Before any of these functions can be applied, the survey data must be prepared as a single analysis file, which the next part describes.

Preparing the Briefing File

All the plots in this lesson use one analysis file, called brief. The code that builds it follows the pattern of the course's data cleaning and regression walkthroughs. It loads the CSCS file from the course's GitHub repository, keeps the 2021 wave (which gives one row per person), and gives short names to six variables. It then uses bracket recoding, in which a condition inside square brackets selects the rows to change, to set the answer "Presented but no response" to missing (NA) for gender and self-rated mental health, and to build four age groups from age in years. The function factor() with a levels argument fixes the order in which categories appear in tables and plots. Finally, na.omit() keeps the people who have no missing value on any of the nine briefing variables.

Table 6.3 lists the nine briefing variables, the survey measure behind each one, its possible values and the way it is treated for plotting.

Table 6.3. The nine variables in the briefing file.

Two cautions apply to every result drawn from this file. Box 6.5 describes how the CSCS recruited its respondents and what its sample can support, and Box 6.6 describes the people whom the complete-case file leaves out.

Activity 6.2 carries out these steps in R. It loads the survey, builds the briefing file, draws a histogram of the loneliness score and a bar chart of self-rated mental health, and summarizes the household size variable that a later part of this section examines.

github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData"
load(url(github_url))
data <- data[data$SURVEY_collection_year == 2021, ]   # the 2021 wave: one row per person
nrow(data)                                             # how many people?

# Short names for the briefing variables
data$loneliness <- data$LONELY_ucla_loneliness_scale_score   # 3 to 9, higher = lonelier
data$support    <- data$PSYCH_zimet_multidimensional_social_support_scale_score   # 1 to 7
data$depression <- data$WELLNESS_phq_score                   # 0 to 6, higher = more symptoms
data$anxiety    <- data$WELLNESS_gad_score                   # 0 to 6, higher = more symptoms
data$life_sat   <- data$WELLNESS_life_satisfaction_num       # 1 to 10, higher = more satisfied
data$age        <- data$DEMO_age                             # years

# Gender: the answer "Presented but no response" becomes missing (NA)
data$gender <- as.character(data$DEMO_gender)
data$gender[data$gender == "Presented but no response"] <- NA
data$gender <- factor(data$gender, levels = c("Woman", "Man", "Non-binary"))

# Self-rated mental health, in order from Poor to Excellent
data$mental_health <- as.character(data$WELLNESS_self_rated_mental_health)
data$mental_health[data$mental_health == "Presented but no response"] <- NA
data$mental_health <- factor(data$mental_health,
  levels = c("Poor", "Fair", "Good", "Very good", "Excellent"))

# Four age groups
data$age_group[data$age < 30] <- "16 to 29"
data$age_group[data$age >= 30 & data$age < 45] <- "30 to 44"
data$age_group[data$age >= 45 & data$age < 65] <- "45 to 64"
data$age_group[data$age >= 65] <- "65 and older"
data$age_group <- factor(data$age_group,
  levels = c("16 to 29", "30 to 44", "45 to 64", "65 and older"))

# The briefing file: the briefing variables for people with no missing values
brief <- na.omit(data[, c("loneliness", "support", "depression", "anxiety",
                          "life_sat", "age", "age_group", "gender", "mental_health")])
nrow(brief)                                            # people in the briefing file

summary(brief$loneliness)              # a numeric summary first
hist(brief$loneliness,
     breaks = seq(2.5, 9.5, by = 1),       # one bar for each score from 3 to 9
     main = "Loneliness (UCLA 3-item scale)",
     xlab = "Loneliness score (3 to 9)", ylab = "Number of people")

table(brief$mental_health)             # the counts that the bars will show
barplot(table(brief$mental_health),
        main = "Self-rated mental health",
        xlab = "Rating", ylab = "Number of people")

summary(data$GEO_housing_household_size)   # people each respondent lives with
hist(data$GEO_housing_household_size,
     breaks = seq(-0.5, 20.5, by = 1),             # one bar for each value from 0 to 20
     main = "Household size",
     xlab = "Number of people the respondent lives with", ylab = "Number of people")
sum(data$GEO_housing_household_size == 20, na.rm = TRUE)   # how many at the top value?

With the briefing file in place, the next part turns to plots of one variable.

Looking at One Variable

The first look at any variable is a plot of its distribution. Figure 6.7 shows the histogram of the loneliness score and the bar chart of self-rated mental health drawn in Activity 6.2, and the two paragraphs that follow read each plot in turn.

Reading the histogram. The argument breaks = seq(2.5, 9.5, by = 1) sets the edges of the bins at 2.5, 3.5, 4.5 and so on, so each bar holds exactly one whole score. Without it, R chooses its own bins, which can merge two scores into one bar or leave empty gaps, and the picture of a variable with only seven values can change a good deal. The loneliness histogram rises from 433 people at the lowest score of 3 to a peak of 755 at a score of 6, and then falls to 179 at a score of 8. The bar at 9, with 270 people (8.8%), is taller than the bar at 8. A rise at the top of a scale is called a ceiling pile-up: a group of people chose the most extreme answer to all three items, and the scale cannot record anything lonelier. The mean (5.57) is a little above the median (5), which reflects the longer upper tail.

Reading the bar chart. The bars follow the order set by factor(), from Poor to Excellent. Very good is the most common rating (1,024 people, 33.2%), followed by Good (928). Together, Poor (230) and Fair (467) make up 697 people, or 22.6% of the briefing file. If the levels had not been set, R would have placed the categories in alphabetical order (Excellent, Fair, Good, Poor, Very good), which hides the shape of an ordinal scale.

Box 6.7 gathers the questions behind these two readings into a checklist that applies to any plot of one variable.

The checklist asks about values outside the possible range and about answers that pile up on particular values. The next part applies it to a variable whose numeric summary looked unremarkable.

Finding a Problem That Cleaning Left Behind

The household size variable, GEO_housing_household_size, records the number of people each respondent lives with, not counting the respondent, so 0 means living alone (503 people). The course's cleaning walkthrough used it to build a living-alone indicator, and its summary in Activity 6.2 looks unremarkable: a median of 3, a mean of 4.46 and a maximum of 20. The histogram tells a different story. Panel A of Figure 6.8 shows that histogram, and panel B shows a second variable from the 2021 wave, the hours worked per week.

In panel A, the bars fall smoothly from one other person toward larger households, and then a bar of 116 people appears at exactly 20, the largest value in the data. The real distribution of household sizes does not rise again at its maximum, so the spike signals that something about the value 20 differs from the other values. Tracing the variable back to the survey explains it. Household size is the sum of eight questions that ask how many partners, children, grandchildren, parents, in-laws, siblings, roommates and other people the respondent lives with, and totals above 20 were recorded as 20. For 105 of the 116 people at 20, the eight answers added up to more than 20, and one total reached 180. Some of these answers are implausible, such as five grandchildren, five in-laws and seven siblings in one household, and they may reflect careless or mistaken responses. The summary statistics did not reveal this problem, and the histogram made it obvious.

The analyst's response has three steps. The first is to check the codebook and the source questions, which is how the meaning of 20 was established. The second is to decide how the variable will be used. Used as a count, the value 20 cannot be trusted and can be set to missing with a bracket recode (Activity 6.2, question 3). Used as a grouped variable, such as living alone versus living with others, or 0, 1, 2 to 4 and 5 or more, the problem is absorbed by the top category and matters much less. The third step is to report the decision in the methods section, with the number of people affected, and, where the decision could change a result, to repeat the analysis both ways as a sensitivity analysis.

Panel B shows two further problems that the same plot can reveal. Answers to the question on hours worked per week heap on round numbers: 340 of the 2,978 people who answered reported exactly 40 hours, compared with one person at 39 and one at 41, and smaller peaks appear at 20, 30 and 35. Heaping of this kind, which Lesson 2 called digit preference, usually reflects rounding and standard working weeks, and it means that small differences in reported hours carry little information. The long right tail reaches 140 hours, and 13 people reported more than 112 hours a week, which would mean working more than 16 hours a day, every day. These values are implausible and would be checked and, in most analyses, set to missing.

Household size and working hours show that a histogram can reveal problems that summary statistics hide. Activity 6.3 draws boxplots and scatterplots of two variables and a scatterplot matrix, and it saves a two-panel figure to a file.

boxplot(loneliness ~ age_group, data = brief,
        main = "Loneliness by age group",
        xlab = "Age group", ylab = "Loneliness score (3 to 9)")
tapply(brief$loneliness, brief$age_group, median)   # the thick line in each box

plot(brief$support, brief$loneliness,              # every person as one point
     xlab = "Social support (1 to 7)", ylab = "Loneliness score (3 to 9)")
set.seed(2021)                                       # the same jitter every time
plot(jitter(brief$support), jitter(brief$loneliness),   # nudge each point a little
     pch = 16, col = rgb(0, 0, 0, 0.15),             # small grey see-through dots
     xlab = "Social support (1 to 7)", ylab = "Loneliness score (3 to 9)")

pairs(brief[, c("loneliness", "support", "depression", "age")],
      pch = 16, col = rgb(0, 0, 0, 0.1))   # every pair of variables in one grid

png("briefing_base_r.png", width = 1600, height = 800, res = 200)   # open a file
par(mfrow = c(1, 2))                       # one row, two plots
hist(brief$loneliness, breaks = seq(2.5, 9.5, by = 1),
     main = "A. Loneliness", xlab = "Loneliness score", ylab = "Number of people")
boxplot(loneliness ~ age_group, data = brief,
        main = "B. Loneliness by age group", xlab = "Age group", ylab = "Loneliness score")
par(mfrow = c(1, 1))
dev.off()                                  # close the file, which saves it

The next part turns from single variables to pairs of variables.

Looking at Two Variables

A plot of two variables shows how the variables go together. Figure 6.9 shows the boxplots of loneliness by age group and the two scatterplots of loneliness against social support from Activity 6.3, and the paragraphs below read each panel.

Boxplots. A boxplot summarizes a numeric variable with five numbers. The thick line is the median. The box runs from the lower to the upper quartile, so it holds the middle half of the people, and its length is the interquartile range (IQR). The whiskers extend to the most extreme values that lie within 1.5 times the IQR of the box, and any value beyond the whiskers is drawn as a separate point. Table 6.4 gives the numbers behind panel A.

Table 6.4. The numbers behind the boxplots in panel A of Figure 6.9: the number of people, the whiskers, the quartiles and the median in each age group.

The two older groups have a median of 6 and the two younger groups a median of 5, and the box for people aged 45 to 64 is the widest, so loneliness varies most in that group. The box for people aged 16 to 29 is narrow, because most of them scored 5 or 6. Its whiskers therefore stop at 4 and 7, and the scores 3, 8 and 9 in this group are drawn as separate points (212 people in all). These points are flagged by the 1.5 × IQR rule, and they are ordinary answers on the scale. With a variable that takes only a few whole values, the points beyond the whiskers say more about the narrowness of the box than about unusual people.

Scatterplots and overplotting. Panel B draws one solid point for each of the 3,083 people. Because the loneliness score takes only seven values, the points fall on seven horizontal lines and cover one another, a problem called overplotting, and the plot cannot show where most people are. Panel C applies two remedies. The function jitter() adds a small random amount to each value so that points sharing a value spread out a little, and col = rgb(0, 0, 0, 0.15) draws each point in black at 15% opacity, so that areas where many points overlap look darker. The function set.seed() fixes the random numbers, so the jittered plot is the same each time the code runs. Jitter changes only the picture: the data used in any calculation are unchanged. Panel C shows that most people have support scores between 4 and 6 and loneliness scores of 5 or 6, and that people with more support tend to report less loneliness.

The scatterplot matrix. The function pairs() draws a small scatterplot for every pair of the chosen variables, with the variable names on the diagonal. Each panel is a quick look at one pair, and the panel above the diagonal shows the same pair as the panel below it with the axes swapped. With 3,083 people and scores that take few values, the panels are crowded, so the matrix is best used to spot strong patterns or unexpected shapes, after which a single scatterplot or a correlation heatmap (Section 4) gives a clearer view.

Figure 6.10 shows the scatterplot matrix that Activity 6.3 draws for loneliness, support, depression and age.

Boxplots, scatterplots and the scatterplot matrix complete the base R displays that this section draws. The last part of this section shows how to arrange several plots in one window and how to keep them in a file.

Several Plots in One Window, and Saving Plots

The setting par(mfrow = c(rows, columns)) divides the plotting window into a grid, and each new plot fills the next cell, row by row. par(mfrow = c(1, 2)) places two plots side by side, and par(mfrow = c(2, 2)), used in the Anscombe activity, places four in a square. The setting stays in force until it is changed, so the code resets it with par(mfrow = c(1, 1)) afterwards.

Plots in the RStudio Plots pane disappear when the session ends. To keep a plot, the analyst sends it to a file. The function png("name.png", width = 1600, height = 800, res = 200) opens an image file 1,600 pixels wide and 800 pixels high, at a resolution of 200 pixels per inch, which prints at 8 by 4 inches. Every plot drawn afterwards goes into the file and does not appear on screen. The function dev.off() closes the file, which saves it in the working directory. The functions pdf() and jpeg() work in the same way for other formats, and the Export button in the Plots pane saves a plot by hand. Saving with code is better practice, because the code records exactly how the figure was made and remakes it in one step if the data change.

The card below links to a narrated R walkthrough of the base R plots in this section and the ggplot2 plots in Section 3.

This section has prepared the briefing file and used base R to read one variable, to find a data problem that cleaning left behind, to show two variables together and to save a figure. The knowledge check and the reflection that follow review these steps, and Section 3 rebuilds the most useful plots as polished figures with ggplot2.

Lesson 6 · Section 3

The Grammar of Graphics with ggplot2

Every ggplot2 plot is built from the same parts, added in layers.

The layered grammar

Six parts of every plot

Data

A data frame with one row per person.

Aesthetics

Mappings from variables to x, y, colour or fill.

Geoms

The marks: bars, boxes, points and lines.

Scales

Rules that turn values into positions and colours.

Facets

Panels that repeat the plot for each group.

Theme

Fonts, grid lines and background.

The template

Data, a mapping, a geom and labels

Line of code	Part of the grammar
ggplot(brief, aes(x = loneliness)) +	Data and the aesthetic mapping
geom_histogram(binwidth = 1, fill = "#0B7B6B", colour = "white") +	The geom, with fixed (set) colours
labs(title = "...", x = "...", y = "...")	Labels

Bars

geom_bar() counts; geom_col() draws values

Groups

A boxplot with every person shown

A smoother

Loneliness across age, with a loess curve

Facets

Small multiples with facet_grid()

Colour

A colour-blind-safe palette

Saving

ggsave() writes a figure at a chosen size

Argument	What it controls
"loneliness_age_gender.png"	The file name; the extension sets the file type.
p	The plot to save; without it, the last plot shown is saved.
width = 7, height = 4.5	The size in inches.
dpi = 300	The resolution in dots per inch, the usual standard for print.

Carry forward

What to take into the next section

• A ggplot2 plot combines data, aesthetic mappings and geoms, with scales, facets and a theme.
• Layers can be stacked, such as a boxplot with jittered points or a scatterplot with a smoother.
• Groups are shown with colour or facets, using a palette that people with colour-vision deficiency can read.

Introduction and Overview

Base R plots are quick to draw, and Section 2 used them to check the data. A figure for the working group's briefing, or for a manuscript, needs more care: consistent labels, readable colours, panels that share a scale, and a file of the right size and resolution. The ggplot2 package (Wickham, 2016) is the most widely used tool in R for such figures. It belongs to the tidyverse, a family of packages that Lesson 1 loaded when setting up R. The walkthroughs and this lesson keep data preparation in base R, exactly as in Section 2, and ggplot2 is the tidyverse package that this lesson uses for figures.

The section first describes the grammar on which ggplot2 is built. It then draws displays of one variable, of a numeric variable across groups and of two numeric variables, compares groups with facets and with colour, and ends with saving a figure to a file.

The Layered Grammar of Graphics

Leland Wilkinson described a grammar of graphics in which every statistical graphic is assembled from a small set of independent parts, in the way that a sentence is assembled from parts of speech (Wilkinson, 2005). Hadley Wickham adapted the idea into a layered grammar, in which a plot is built by adding layers one at a time, and implemented it as ggplot2 (Wickham, 2010). The practical benefit is that one template covers almost every display: once the parts are familiar, a histogram, a boxplot and a faceted scatterplot differ only in a line or two of code.

Figure 6.11 lists the six parts of a ggplot2 plot, each with an example taken from the briefing.

The six flip cards below repeat the parts in Figure 6.11 and serve as a reference for the code in the activities of this section.

Click each card for a short description of the part and the code that controls it.

Equation 6.1 shows how the parts are joined into a single block of code.

Every ggplot2 plot follows the same template. The function ggplot() receives the data and the aesthetic mappings, and each further part is added with a plus sign at the end of the line before it:

The accordion below describes three errors that beginners meet first when they use the template in Equation 6.1, with the remedy for each.

With the template in hand, the following parts apply it to the displays drawn in base R in Section 2, beginning with one variable.

One Variable in ggplot2

The histogram in Activity 6.4 maps loneliness to x and adds geom_histogram(binwidth = 1), which gives one bar per whole score, as breaks did in base R. The fill colour and the white bar outlines are set outside aes(), so they apply to every bar. Figure 6.12 carries the same information as the base R histogram in Section 2 (Figure 6.7), with labels that can be reused in the briefing.

Two geoms draw bars, and they are easily confused. geom_bar() counts the rows in each category, so it needs only an x aesthetic, such as mental_health. geom_col() draws bars whose heights are supplied in the data, so it needs both x and y. To draw mean loneliness by age group, the means are first computed with the base R function aggregate(), which returns a small data frame with one row per age group, and geom_col() then draws those values. Figure 6.13 draws one chart with each geom.

The right-hand panel is the "means without spread" display from Section 1, and it is shown here to illustrate geom_col(). For the briefing, loneliness by age group is better shown with boxplots and the individual points.

Activity 6.4 draws the plots in Figure 6.12 and Figure 6.13 and a boxplot with jittered points.

# install.packages("ggplot2")   # run once, if not yet installed
library(ggplot2)
ggplot(brief, aes(x = loneliness)) +          # data and the aesthetic mapping
  geom_histogram(binwidth = 1,                # the geom: one bar per score
                 fill = "#0B7B6B", colour = "white") +
  labs(title = "Loneliness scores",           # labels
       x = "Loneliness score (3 to 9)", y = "Number of people")

ggplot(brief, aes(x = mental_health)) +
  geom_bar() +                                # geom_bar() counts the rows in each category
  labs(x = "Self-rated mental health", y = "Number of people")

means <- aggregate(loneliness ~ age_group, data = brief, FUN = mean)
means                                         # one row per age group
ggplot(means, aes(x = age_group, y = loneliness)) +
  geom_col() +                                # geom_col() draws the values it is given
  labs(x = "Age group", y = "Mean loneliness score")

set.seed(2021)                                # the same jitter every time
ggplot(brief, aes(x = age_group, y = loneliness)) +
  geom_boxplot(outlier.shape = NA) +          # the box summarizes each group
  geom_jitter(width = 0.2, height = 0.2, alpha = 0.1) +   # every person as a faint dot
  labs(x = "Age group", y = "Loneliness score (3 to 9)")

The last plot in Activity 6.4 combines a boxplot with the individual points, a display that the next part examines more closely.

A Numeric Variable Across Groups

A boxplot with jittered points uses two geoms on the same axes. geom_boxplot(outlier.shape = NA) draws the median, the quartiles and the whiskers for each age group and suppresses the separate outlier points, because geom_jitter(width = 0.2, height = 0.2, alpha = 0.1) already draws every person. The width and height arguments set how far each point may be nudged sideways and up or down, and alpha sets the opacity, from 0 (invisible) to 1 (solid). Figure 6.14 shows the result for the four age groups.

The combined display shows three things that bars of means would hide. The medians differ by one point (5 in the two younger groups and 6 in the two older groups). Every group includes people across the whole range of the scale, from 3 to 9. And the groups differ in size, from 1,079 people aged 16 to 29 to 366 aged 65 and older, which the density of the points makes visible.

Boxplots compare a numeric variable across the groups of a categorical one. When both variables are numeric, the display is a scatterplot, which the next part extends with a smoother.

Two Numeric Variables and a Smoother

A scatterplot maps two numeric variables to x and y. Because the loneliness score takes only seven values, the points are jittered vertically (height = 0.2) and drawn with transparency. The layer geom_smooth() adds a curve that summarizes the average of y across the values of x, with a grey 95% confidence band. With method = "loess", the curve is a loess smoother (locally estimated scatterplot smoothing), which fits many small weighted regressions along the axis and joins them, so it can bend. With method = "lm", the curve is a straight least-squares line. When it draws the curve, ggplot2 prints the message `geom_smooth()` using formula = 'y ~ x', which only reports the formula it used and is not an error. Figure 6.15 adds a loess smoother to a scatterplot of loneliness against age.

The curve tells a story that a straight line would miss. Predicted loneliness is highest among the youngest respondents (about 6.4 at age 16), falls to about 5.3 around age 30, rises to about 6.0 between ages 50 and 60, and declines after 70 (about 5.4 at age 80). The band is narrow where most people are (ages 20 to 40) and widens at both ends, where the curve rests on few people: 48 respondents are under 20 and 29 are 80 or older. The ends of a smoother deserve the least trust, and a single respondent aged 100 sits at the far right. A smoother is a description of the data. It suggests that age might enter a later model as a curve or in groups, and it gives no reason, by itself, for the pattern.

The smoother summarizes one relationship across the whole briefing file. The next part asks whether a pattern differs between groups by drawing one panel for each group.

Comparing Across Groups with Facets

Facets draw the same plot in one panel per group. facet_wrap(~ age_group) takes one categorical variable and wraps its panels into a grid, which suits a single grouping variable with several levels. All panels share the same axes by default, so heights and positions can be compared across panels. Figure 6.16 applies facet_wrap() to the histogram of loneliness, with one panel for each age group.

The shared vertical axis shows both the shape and the size of each group. The two younger groups, with more than 1,000 people each, have tall bars at 5 and 6. The two older groups have shorter, flatter histograms with more people at both ends of the scale, and they hold most of the ceiling pile-up seen in Section 2: in the 45 to 64 group, 111 people (19.2%) scored 9, as many as scored 6 (112), and in the 65 and older group the most common score is 3, the lowest possible (85 people, 23.2%). The older groups are therefore more divided, with more people at the least lonely and the most lonely ends of the scale. If the groups' shapes matter more than their sizes, facet_wrap(~ age_group, scales = "free_y") gives each panel its own vertical axis, at the cost of making heights incomparable across panels.

facet_grid(gender ~ age_group) takes two variables, the first for the rows and the second for the columns, and draws one panel for every combination. Figure 6.17 shows social support against loneliness, with a straight line in each panel.

The lines slope downward in all twelve panels, so higher support goes with lower loneliness in every combination of gender and age group. The lines look steeper in the two older age groups: within the age groups as a whole, the slopes are about −0.40 and −0.54 points of loneliness per point of support in the two younger groups and about −0.76 and −0.84 in the two older groups. The table of counts printed in Activity 6.5 is the necessary companion to the grid. The non-binary panels hold between 6 and 26 people, so their lines and bands are very uncertain, and any difference they appear to show would need a much larger sample to assess. An exploratory figure should show small groups, and the text should say how small they are.

Table 6.5 gives these counts, the number of people behind each panel of Figure 6.17.

Table 6.5. Number of people in each panel of Figure 6.17, by gender and age group.

Facets give each group its own panel. Groups can also share one panel and be told apart by colour, and the next part turns to colour, palettes and themes.

Colour, Palettes and Themes

Mapping a variable to colour (for points and lines) or fill (for bars and boxes) shows groups within one panel. Colour and facets suit different situations. Colour lets the reader compare groups directly on the same axes, and it works well for two or three groups. Facets give each group its own panel, and they work better when there are many groups, when points overlap heavily, or when one group is much smaller than the others and would be hidden. The two can also be combined, with colour for one grouping variable and facets for another.

The colours themselves need care. About 8% of men and 0.5% of women of northern European ancestry have some form of red-green colour-vision deficiency, so a palette that relies on telling red from green fails a sizeable part of any audience. The Okabe-Ito palette (Okabe & Ito, 2008) is a set of eight colours chosen to remain distinguishable for people with the common forms of colour-vision deficiency. Its first three colours after black, orange (#E69F00), sky blue (#56B4E9) and bluish green (#009E73), are used below with scale_fill_manual(). The viridis palettes, available in ggplot2 as scale_fill_viridis_d() for categories, are another safe choice, and they also remain readable when printed in greyscale. Figure 6.18 maps fill to gender in boxplots of loneliness by age group, with the three Okabe-Ito colours and a minimal theme.

The theme controls the appearance of everything except the data. The default theme, theme_grey(), uses a grey panel with white grid lines. theme_minimal() removes the grey panel, theme_bw() adds a black border, and theme_classic() keeps only the axis lines. The argument base_size sets the size of all text in points, and 11 or 12 suits most journal figures. Labels are set with labs(), which also names the legend through the aesthetic it describes, as in fill = "Gender".

Once a figure has its colours, theme and labels, it is saved to a file at the size at which it will be printed or shown, which the next part describes.

Saving a Figure with ggsave()

The function ggsave() writes a ggplot2 plot to a file. Its first argument is the file name, and the extension sets the format: .png for most uses, .pdf or .svg for figures that must stay sharp at any size, and .tiff when a journal asks for it. The second argument is the plot object, here p, and without it ggsave() saves the last plot shown. The arguments width and height are in inches by default, and dpi sets the resolution, with 300 dots per inch the usual requirement for print. A figure of 7 by 4.5 inches at 300 dpi is 2,100 by 1,350 pixels. Text size is relative to the physical size of the figure, so a plot saved at a large width with a small base_size will have text that is too small when the figure is shrunk to fit a page.

Activity 6.5 draws the plots in Figure 6.15, Figure 6.16, Figure 6.17 and Figure 6.18, prints the counts in Table 6.5, and saves the last plot with ggsave().

set.seed(2021)
ggplot(brief, aes(x = age, y = loneliness)) +
  geom_jitter(height = 0.2, alpha = 0.1) +
  geom_smooth(method = "loess") +             # a smooth curve through the average score
  labs(x = "Age (years)", y = "Loneliness score (3 to 9)")

ggplot(brief, aes(x = loneliness)) +
  geom_histogram(binwidth = 1, fill = "#0B7B6B", colour = "white") +
  facet_wrap(~ age_group) +                   # one panel for each age group
  labs(x = "Loneliness score (3 to 9)", y = "Number of people")

set.seed(2021)
ggplot(brief, aes(x = support, y = loneliness)) +
  geom_jitter(height = 0.2, alpha = 0.2) +
  geom_smooth(method = "lm") +                # a straight line in each panel
  facet_grid(gender ~ age_group) +            # rows: gender; columns: age group
  labs(x = "Social support (1 to 7)", y = "Loneliness score (3 to 9)")

table(brief$gender, brief$age_group)         # people in each panel

p <- ggplot(brief, aes(x = age_group, y = loneliness, fill = gender)) +
  geom_boxplot() +
  scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73")) +  # Okabe-Ito colours
  labs(x = "Age group", y = "Loneliness score (3 to 9)", fill = "Gender") +
  theme_minimal(base_size = 12)               # a plain theme with larger text
p                                             # show the plot
ggsave("loneliness_age_gender.png", p, width = 7, height = 4.5, dpi = 300)

The card below links again to the narrated walkthrough introduced in Section 2.

This section has described the grammar on which ggplot2 is built and used it to draw the briefing plots for one variable, for a numeric variable across groups and for two numeric variables, to compare groups with facets and colour, and to save a figure. The knowledge check and the reflection that follow review these steps, and Section 4 adds displays for many variables and combines the briefing plots into one figure.

Lesson 6 · Section 4

Many Variables and Combined Figures

An overview of many variables, and one figure that brings the briefing together.

Correlation

Pearson's r for every pair of variables

Loneliness with	r	Direction
Depressive symptoms (PHQ-2)	0.47	Positive
Anxiety symptoms (GAD-2)	0.43	Positive
Social support	−0.40	Negative
Life satisfaction	−0.37	Negative
Age	0.07	Close to zero

The strongest correlation in the matrix is between depression and anxiety (r = 0.59).

Heatmap

A correlation matrix drawn as colours

Interpretation

What a correlation heatmap cannot show

Cause

A correlation fits several causal stories, including reverse causation and common causes.

Curves

Pearson's r measures straight-line association, so a curved pattern can give r near zero.

Shape

One number per pair hides outliers, clusters and the number of people involved.

Marginal histograms

A relationship and two distributions in one display

Combining panels

The briefing as one figure

Publication figures

Conventions for a manuscript figure

• Each axis has a label in words, with units or the range of the scale.
• Plain labels replace R variable names such as life_sat.
• Panels carry letters, and legends say what colours and marks mean.
• Colours remain readable with colour-vision deficiency and in greyscale.
• Text is legible at the printed size, and the file is saved at about 300 dpi.

Captions

What a figure caption states

Element	Example from the briefing figure
What the figure shows	Loneliness and social connection
Data source and year	Canadian Social Connection Survey, 2021
People included	3,083 respondents with complete data
Measures and ranges	UCLA 3-item scale (3 to 9); support (1 to 7)
What the marks mean	Boxes, whiskers, jitter, line and 95% band

Lesson summary

From a first look to a finished figure

Step	Tool
Look before modelling	Plots of every variable and relationship
Choose the display	The question and the variable types
Check the data quickly	Base R: hist(), boxplot(), plot(), pairs()
Build the figure	ggplot2: data, aesthetics, geoms, facets, themes
Combine and caption	patchwork, ggsave() and a caption that stands alone

Introduction and Overview

The working group's third question asks how loneliness goes together with social support and mental health, which involves several variables at once. One scatterplot per pair would answer it slowly, because six variables have fifteen pairs. This section shows how to summarize many pairs in one display with a correlation matrix and a heatmap, and it explains the limits of that summary. It then adds marginal histograms to a scatterplot, combines the briefing plots into one multi-panel figure, and sets out the conventions of a publication figure and its caption. The finished figure is the kind of exploratory figure that a research paper based on survey data often includes.

The section begins with the correlation matrix, the numeric summary from which the heatmap is drawn.

Looking at Many Variables at Once

The Pearson correlation coefficient, written r, measures how closely two numeric variables follow a straight line. It ranges from −1 (all points on a line sloping downward) through 0 (no straight-line relationship) to +1 (all points on a line sloping upward). Its sign gives the direction and its size gives the strength. Cohen (1988) suggested that values of about 0.1, 0.3 and 0.5 be called small, medium and large in the behavioural sciences, and these labels are a rough guide only, since what counts as a strong association depends on the field and the question.

The function cor(), given a data frame of numeric columns, returns a correlation matrix with one row and one column per variable. Each cell holds the correlation between its row variable and its column variable. The diagonal holds the correlation of each variable with itself, which is always 1, and the matrix is symmetric, so the values above the diagonal repeat those below it. Six variables give 36 cells and 15 distinct correlations.

Table 6.6 gives the correlation matrix for the six numeric briefing variables.

Table 6.6. Pearson correlations among the six numeric variables in the briefing file (n = 3,083).

In the briefing file, loneliness has moderate positive correlations with depressive symptoms (0.47) and anxiety symptoms (0.43), and moderate negative correlations with social support (−0.40) and life satisfaction (−0.37). Its correlation with age is close to zero (0.07). The strongest correlation in the matrix is between depressive and anxiety symptoms (0.59), two measures that often move together.

Two practical points apply to cor(). First, it returns NA for any pair with a missing value unless told how to handle them, for example with use = "complete.obs". The briefing file has no missing values, so this does not arise here. Second, Pearson's r assumes numeric variables and is sensitive to outliers and to skew. For ordinal scores or skewed variables, Spearman's rank correlation, computed with cor(..., method = "spearman"), uses the ranks of the values and is less affected by extreme values. For loneliness and social support, Spearman's correlation is −0.36, close to the Pearson value of −0.40.

Activity 6.6 computes Table 6.6 with cor(), reshapes it into one row per pair of variables, and draws a heatmap of the correlations.

vars <- c("loneliness", "support", "depression", "anxiety", "life_sat", "age")
cor_mat <- round(cor(brief[, vars]), 2)     # Pearson correlations, rounded
cor_mat

cor_long <- as.data.frame(as.table(cor_mat))   # one row for each pair
names(cor_long) <- c("var1", "var2", "r")
head(cor_long)

ggplot(cor_long, aes(x = var1, y = var2, fill = r)) +
  geom_tile(colour = "white") +                     # one coloured square per pair
  geom_text(aes(label = r), size = 3.5) +           # print the value on each square
  scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       limits = c(-1, 1)) +         # blue negative, red positive
  labs(x = NULL, y = NULL, fill = "r") +
  theme_minimal()

The heatmap at the end of Activity 6.6 displays the same matrix in colour, and the next part explains how it is built and how it is read.

Correlation Heatmaps

A matrix of 36 numbers is slow to read. A correlation heatmap draws each cell as a coloured square, so that the pattern can be seen at a glance. Drawing it in ggplot2 takes two steps. The matrix is first turned into a long data frame with one row per pair: as.table() treats the matrix as a table, and as.data.frame() turns that table into three columns, here renamed var1, var2 and r. Then geom_tile() draws one square per row, with var1 mapped to x, var2 to y and r to fill, and geom_text() prints the value on each square. Figure 6.19 shows the heatmap for the six briefing variables.

The colour scale is a diverging scale, set with scale_fill_gradient2(): blue for negative values, white at zero and red for positive values, with limits = c(-1, 1) so that the colours keep the same meaning in any heatmap. A diverging scale suits correlations because zero is a meaningful midpoint. Blue and red remain distinguishable for people with red-green colour-vision deficiency, which a red-to-green scale would not.

The heatmap shows a clear structure. Loneliness, depression and anxiety form a block of positive correlations, and each of the three is negatively correlated with social support and life satisfaction, which are positively correlated with each other (0.35). Age is close to white in every square. For the briefing, this panel gives the working group an overview of the third question in one picture.

The overview has limits, which the following subsection sets out.

What a Heatmap Shows and What It Cannot

A correlation heatmap shows the strength and direction of straight-line association between each pair of variables. Three common readings go beyond what it shows.

The accordion below explains each of these readings and why the heatmap cannot support it.

A heatmap therefore points to the pairs of variables that deserve a closer look. The next part shows one way to take that closer look: a scatterplot of one pair with the distribution of each variable along its margins.

Scatterplots with Marginal Histograms

A scatterplot shows how two variables go together, and a histogram shows how one variable is distributed. A scatterplot with marginal histograms shows both at once: the scatterplot sits in the middle, with a histogram of the horizontal variable along the top and a histogram of the vertical variable along the right side. The function ggMarginal() from the ggExtra package adds the margins to a finished ggplot2 scatterplot. Its argument type can be "histogram", "density" or "boxplot". Figure 6.20 shows loneliness against social support with a histogram of each variable along its margin.

The marginal histogram of social support shows a distribution skewed to the left: most people have scores between 4 and 6, and a long, thin tail runs toward low support. It also shows spikes at whole-number scores. In the briefing file, 170 people have a score of exactly 5, 136 exactly 4 and 117 exactly 7. The support score is the mean of twelve items, so a whole-number score is what a person who gives the same answer to every item receives, and a score of 7 requires the highest answer on all twelve. Such answers can reflect genuine views or a respondent moving quickly through the questions (sometimes called straight-lining), and the plot is a reason to look more closely before analyzing the score. The marginal histogram of loneliness shows the seven separate values of that score.

The object that ggMarginal() returns is a finished drawing, and it is not a ggplot2 object, so further layers cannot be added to it and patchwork cannot place it in a grid without extra steps. It can be saved with ggsave(). The combined briefing figure (Figure 6.21) therefore uses the scatterplot without its margins.

The marginal histograms reveal spikes in the support score at whole-number values. The next part places the scatterplot beside three other briefing plots in a single figure made with patchwork.

Combining Panels with patchwork

A multi-panel figure lets a reader see several related displays together, with consistent styling and one caption. The patchwork package combines ggplot2 plots with arithmetic-like operators. A plus sign (+) or a vertical bar (|) places plots side by side, a forward slash (/) stacks them, and parentheses group them. The expression (p_hist + p_box) / (p_scatter + p_heat) therefore gives two plots in the top row and two in the bottom row. The function plot_annotation() adds a title, a subtitle or a caption to the whole figure, and its argument tag_levels = "A" labels the panels A, B, C and D in reading order, which the caption can then refer to. The function plot_layout() controls the number of columns, the relative widths and heights of panels, and whether legends are collected in one place. Figure 6.21 shows the draft briefing figure that this expression produces from the histogram, the boxplots, the scatterplot and the heatmap.

The figure is saved with ggsave("briefing_figure.png", figure, width = 11, height = 8.5, dpi = 300). A four-panel figure needs a larger size than a single plot, so that each panel keeps legible text, and 11 by 8.5 inches suits a full page turned sideways or a slide. When the figure is placed in a narrower space, such as one column of a journal, it is better to save it at that width and adjust the text size with base_size than to shrink a large image.

Activity 6.7 draws the scatterplot with marginal histograms in Figure 6.20, and then assembles and saves the four-panel briefing figure in Figure 6.21.

# install.packages(c("ggExtra", "patchwork"))   # run once, if not yet installed
library(ggExtra)     # ggMarginal() adds histograms to the margins of a plot
library(patchwork)   # + and / combine ggplot2 plots into one figure
set.seed(2021)
p_scatter <- ggplot(brief, aes(x = support, y = loneliness)) +
  geom_jitter(height = 0.2, width = 0, alpha = 0.15) +
  geom_smooth(method = "lm", colour = "#CC0033") +
  labs(x = "Social support (1 to 7)", y = "Loneliness score (3 to 9)")
ggMarginal(p_scatter, type = "histogram")    # a histogram on each axis

p_hist <- ggplot(brief, aes(x = loneliness)) +
  geom_histogram(binwidth = 1, fill = "#0B7B6B", colour = "white") +
  labs(x = "Loneliness score (3 to 9)", y = "Number of people")
p_box <- ggplot(brief, aes(x = age_group, y = loneliness)) +
  geom_boxplot(fill = "#E6F3F0") +
  labs(x = "Age group", y = "Loneliness score (3 to 9)")
p_heat <- ggplot(cor_long, aes(x = var1, y = var2, fill = r)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = r), size = 3) +
  scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       limits = c(-1, 1)) +
  labs(x = NULL, y = NULL, fill = "r")

figure <- (p_hist + p_box) / (p_scatter + p_heat) +   # + side by side, / stacked
  plot_annotation(tag_levels = "A",                     # label the panels A to D
    title = "Loneliness and social connection, CSCS 2021 (n = 3,083)")
figure
ggsave("briefing_figure.png", figure, width = 11, height = 8.5, dpi = 300)

The draft figure still shows R variable names in panel D. The next part sets out the conventions that a publication figure follows.

Conventions of a Publication Figure

An exploratory figure in a manuscript has a different audience from the analyst's quick looks. Readers see it once, often without the text nearby, and they should be able to understand it without access to the code. Tufte (1983) argued that a graphic should spend its ink on the data and remove decoration that carries no information, which he called chartjunk. Table 6.7 lists the conventions most journals expect and checks the draft briefing figure against them.

Table 6.7. Conventions of a publication figure, checked against the draft briefing figure (Figure 6.21).

A figure that meets the conventions in Table 6.7 still needs a caption, which the next part describes.

Writing a Figure Caption

A caption makes a figure stand on its own. Most journals place the caption below the figure, beginning with the figure number and a short title, followed by sentences that tell the reader what is shown. A good caption states what the figure displays and for whom: the data source and year, the number of people included and the reason others were excluded, the measures and their ranges, and the meaning of every panel and mark, such as what the boxes, whiskers, lines and bands represent and whether points were jittered. The caption describes the figure; the interpretation, such as which differences matter, belongs in the results text that refers to it.

Worked Example 6.1 applies these rules to the briefing. It lists the changes made to the draft in Figure 6.21, gives the full caption of the finished figure, and shows the three sentences of text that accompany the figure in the briefing.

The caption and the three sentences of text complete the briefing. The last part of this section places a figure of this kind within the results of a research paper.

From Analysis to Manuscript

An exploratory figure is one part of the results of a research paper. STROBE, the reporting guideline for cohort, case-control and cross-sectional studies (von Elm et al., 2007), lists what the results section should report. Item 13 asks for the number of people at each stage of the study and the reasons they were not included, which the caption of the briefing figure gives (3,083 of 4,045 respondents, with 962 excluded for missing data). Item 14 asks for the characteristics of the participants and the number with missing data for each variable, usually in a Table 1, and item 15 asks for summaries of the outcome, such as the distribution in panel A. Item 16 asks for the main results: unadjusted and adjusted estimates with their confidence intervals, and a statement of which confounders were adjusted for and why. These estimates are usually set out in a results table with one column for the crude estimates and one for the adjusted estimates, as Lesson 4, Section 1 showed for a logistic model.

The figure and the tables then support the rest of the paper. The abstract reports the main estimate with its confidence interval in one or two sentences. The discussion interprets the results, compares them with earlier studies and states the limitations; for the CSCS, these include the cross-sectional design, which limits causal claims, and the online recruitment, which limits generalization to the Canadian population. For optional reading, HSCI 207 Lesson 12 (Writing Up Research) treats results tables, the discussion and limitations in more detail.

This section has summarized many variables with a correlation matrix and a heatmap, set out the limits of that summary, added marginal histograms to a scatterplot, combined the briefing plots into one figure, and described the conventions of a publication figure and its caption. The knowledge check and the reflection that follow review these ideas, and the final page of the lesson gathers the key takeaways, a closing reflection and the final assessment.

Bringing It All Together

This lesson followed one request from start to finish: a one-page visual briefing on social connection for a regional loneliness working group, drawn from the 2021 wave of the Canadian Social Connection Survey. Section 1 set out why the analyst looks before modelling. Anscombe's quartet showed four datasets with the same means, standard deviations, correlation and fitted line but completely different patterns, and a decision guide matched the display to the number and types of the variables in a question. A gallery of poor displays showed what a truncated axis, a bar chart of means, a crowded pie chart and an overplotted scatterplot conceal.

Section 2 prepared the briefing file in base R (3,083 people with complete data on nine variables) and used quick base R plots to check it. A histogram revealed a spike of 116 people at the top value of household size, where larger totals had been recorded as 20, and heaping at 40 hours in reported hours of work, problems that the numeric summaries did not show. Section 3 rebuilt the plots with ggplot2, using the layered grammar of data, aesthetics, geoms, scales, facets and themes. A loess smoother showed that loneliness dips around age 30 and rises to a second high point in the fifties, facets compared the relationship between social support and loneliness across gender and age groups, and the Okabe-Ito palette kept the colours readable. Section 4 summarized many variables with a correlation heatmap, explained why a heatmap shows association and cannot show cause, curves or shape, added marginal histograms to a scatterplot, and combined four panels into one figure with a caption.

The same habits apply to every analysis in the rest of the course: each variable is plotted before it is summarized, each relationship is plotted before it is modelled, and each figure carries labels and a caption that let it stand on its own. The next lesson, Measurement and Psychometrics, takes up a question that this lesson left open: whether the scale scores plotted here, such as the loneliness and social support scores, measure what they are meant to measure, and how reliably they do so.

The final assessment covers all four sections. All 15 questions must be answered correctly (100%), and the final reflection completed, to finish the lesson.

Final Knowledge Assessment