6  Correlation

Functions introduced in this chapter

cor, fct_recode

6.1 Introduction

In this chapter, we will learn about the concept of correlation, which is a way of measuring a linear relationship between two numerical variables.

6.1.1 Install new packages

Type the following at the Console:

install.packages("faraway")
install.packages("openintro")

6.1.2 Download the Quarto file

Look at either the top (Posit Cloud) or the upper right corner of the RStudio screen to make sure you are in your intro_stats project.

Then click on the following link to download this chapter as a Quarto file (.qmd).

https://vectorposse.github.io/intro_stats/chapter_downloads/06-correlation.qmd

Once the file is downloaded, move it to your project folder in RStudio and open it there.

6.1.3 Restart R and run all chunks

In RStudio, select “Restart R and Run All Chunks” from the “Run” menu.

6.1.4 Load packages

We load the now-standard tidyverse package. We also include the faraway package to access data about Chicago in the 1970s, the openintro package to access the bdims data set on body dimensions, and the palmerpenguins package to access the familiar penguins data.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(faraway)
library(openintro)

Loading required package: airports
Loading required package: cherryblossom
Loading required package: usdata

Attaching package: 'openintro'

The following objects are masked from 'package:faraway':

    orings, twins

library(palmerpenguins)

6.2 Redlining in Chicago

The data set we will use throughout this chapter is from Chicago in the 1970s studying the practice of “redlining”.

Exercise 1

Do an internet search for “redlining”.

Consult at least two or three sources. Then, in your own words (not copied and pasted from any of the websites you consulted), explain what “redlining” means.

Please write up your answer here.

---

The chredlin data set appears in the faraway package accompanying a book by Julian Faraway (Practical Regression and Anova using R, 2002.) Faraway explains:

“In a study of insurance availability in Chicago, the U.S. Commission on Civil Rights attempted to examine charges by several community organizations that insurance companies were redlining their neighborhoods, i.e. canceling policies or refusing to insure or renew. First the Illinois Department of Insurance provided the number of cancellations, non-renewals, new policies, and renewals of homeowners and residential fire insurance policies by ZIP code for the months of December 1977 through February 1978. The companies that provided this information account for more than 70% of the homeowners insurance policies written in the City of Chicago. The department also supplied the number of FAIR plan policies written and renewed in Chicago by zip code for the months of December 1977 through May 1978. Since most FAIR plan policyholders secure such coverage only after they have been rejected by the voluntary market, rather than as a result of a preference for that type of insurance, the distribution of FAIR plan policies is another measure of insurance availability in the voluntary market.”

In other words, the degree to which residents obtained FAIR policies can be seen as an indirect measure of redlining. This participation in an “involuntary” market is thought to be largely driven by rejection of coverage under more traditional insurance plans.

6.2.1 Exploratory data analysis

Before we learn about correlation, let’s get to know our data a little better.

Type ?chredlin at the Console to read the help file. While it’s not very informative about how the data was collected, it does have crucial information about the way the data is structured.

Here is the data set:

chredlin

Exercise 2

What do each of the rows of this data set represent? You’ll need to refer to the help file. (They are not individual people.)

Please write up your answer here.

Exercise 3

The race variable is numeric. Why? What do these numbers represent? (Again, refer to the help file.)

Please write up your answer here.

---

The glimpse command gives a concise overview of all the variables present.

glimpse(chredlin)

Rows: 47
Columns: 7
$ race     <dbl> 10.0, 22.2, 19.6, 17.3, 24.5, 54.0, 4.9, 7.1, 5.3, 21.5, 43.1…
$ fire     <dbl> 6.2, 9.5, 10.5, 7.7, 8.6, 34.1, 11.0, 6.9, 7.3, 15.1, 29.1, 2…
$ theft    <dbl> 29, 44, 36, 37, 53, 68, 75, 18, 31, 25, 34, 14, 11, 11, 22, 1…
$ age      <dbl> 60.4, 76.5, 73.5, 66.9, 81.4, 52.6, 42.6, 78.5, 90.1, 89.8, 8…
$ involact <dbl> 0.0, 0.1, 1.2, 0.5, 0.7, 0.3, 0.0, 0.0, 0.4, 1.1, 1.9, 0.0, 0…
$ income   <dbl> 11.744, 9.323, 9.948, 10.656, 9.730, 8.231, 21.480, 11.104, 1…
$ side     <fct> n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n…

Exercise 4(a)

Which variable listed above represents participation in the FAIR plan? How is it measured? (Again, refer to the help file.)

Please write up your answer here.

Exercise 4(b)

Why is it important to analyze the number of plans per 100 housing units as opposed to the total number of plans across each ZIP code? (Hint: what happens if some ZIP codes are larger than others?)

Please write up your answer here.

---

We are interested in the association between race and involact. If redlining plays a role in driving people toward FAIR plan policies, we would expect there to be a relationship between the racial composition of a ZIP code and the rate of FAIR plan policies obtained in that ZIP code.

Exercise 5(a)

Since race is a numerical variable, what type of graph or chart is appropriate for visualizing it? (You may need to refer back to the “Numerical data” chapter.)

Please write up your answer here.

Exercise 5(b)

Using ggplot code, create the type of graph you identified above. After creating the initial plot, be sure to go back and set the binwidth and boundary to sensible values. (Refer back to the “Numerical data” chapter for sample code if you’ve forgotten how to make such a graph. If you were unsure about part (a), the instructions about binwidth and boundary should be a pretty big hint.)

# Add code here to create a plot of race

Exercise 5(c)

Describe the shape of the race variable using the three key shape descriptors (modes, symmetry, and outliers).

Please write up your answer here.

Exercise 5(d)

Create the same kind of graph as above, but for involact. (Again, go back and set the binwidth and boundary to sensible values.)

# Add code here to create a plot of involact

Exercise 5(e)

Describe the shape of the involact variable using the three key shape descriptors (modes, symmetry, and outliers).

Please write up your answer here.

Exercise 5(f)

Since both race and involact are numerical variables, what type of graph or chart is appropriate for visualizing the relationship between them?

Please write up your answer here.

Exercise 5(g)

For our research question, is race functioning as a predictor variable or as the response variable? What about involact? Why? Explain why it makes more sense to think of one of them as the predictor and the other as the response.

Please write up your answer here.

Exercise 5(h)

Using ggplot code, create the type of graph you identified above. Be sure to put involact on the y-axis and race on the x-axis. (Again, that’s a hint in case you were confused in part (g).)

# Add code here to create a plot of involact against race

6.3 Correlation

The word correlation describes a linear relationship between two numerical variables. As long as certain conditions are met, we can calculate a statistic called the correlation coefficient, often denoted with a lowercase r.

There are several different ways to compute a statistic that measures correlation. The most common way, and the way we will learn in this chapter, is often attributed to an English mathematician named Karl Pearson. According to his Wikipedia page,

“Pearson was also a proponent of social Darwinism, eugenics and scientific racism.”

Exercise 6

Do an internet search for each of the following terms:

• Social Darwinism
• Eugenics
• Scientific racism

Consult at least two or three sources for each term. Then, in your own words (not copied and pasted from any of the websites you consulted), explain each of these terms.

Please write up your answer here.

---

While Pearson is often credited with its discovery, the so-called “Pearson correlation coefficient” was first developed by a French scientist, Auguste Bravais. Due to the misattribution of discovery, along with the desire to disassociate the useful tool of correlation from its problematic applications to racism and eugenics, we will just refer to it as the correlation coefficient (without a name attached).

The correlation coefficient, r, has some important properties.

• The correlation coefficient is a number between -1 and 1.
• A value close to 0 indicates little or no correlation.
• A value close to 1 indicates strong positive correlation.
• A value close to -1 indicates strong negative correlation.

In between 0 and 1 (or -1), we often use words like weak, moderately weak, moderate, and moderately strong. There are no exact cutoffs for when such words apply. You must learn from experience how to judge scatterplots and r values to make such determinations.

A correlation is positive when low values of one variable are associated with low values of the other value. Similarly, high values of one variable are associated with high values of the other. For example, exercise is positively correlated with burning calories. Low exercise levels will burn a few calories; high exercise levels burn more calories, on average.

A correlation is negative when low values of one variable are associated with high values of the other value, and vice versa. For example, tooth brushing is negatively correlated with cavities. Less tooth brushing may result in more cavities; more tooth brushing is associated with fewer cavities, on average.

6.4 Conditions for correlation

Two variables are considered “associated” any time there is any type of relationship between them (i.e., they are not independent). However, in statistics, we reserve the word “correlation” for situations meeting more stringent conditions:

1. The two variables must be numerical.¹
2. There is a somewhat linear relationship between the variables, as shown in a scatterplot.
3. There are no serious outliers.

For condition (2) above, keep in mind that real data in scatterplots very rarely lines up in a perfect straight line. Instead, you will see a “cloud” of dots. All we want to know is whether that cloud of dots mostly moves from one corner of the scatterplot to the other. Violations of this condition will usually be for one of two reasons:

• The dots are scattered completely randomly with no discernible pattern.
• The dots have a pattern or shape to them, but that shape is curved and not linear.

Exercise 7

Check the three conditions for the relationship between involact and race. For conditions (2) and (3), you’ll need to check the scatterplot you created above. (You did create a scatterplot for one of the exercises above, right?)

Please write up your answer here.

6.5 Calculating correlation

Since the conditions are met, We calculate the correlation coefficient using the cor command.

cor(chredlin$race, chredlin$involact)

[1] 0.713754

The order of the variables doesn’t matter; correlation is symmetric, so the r value is the same independent of the choice of response and predictor variables.

Since the correlation between involact and race is a positive number and slightly closer to 1 than 0, we might call this a “moderate” positive correlation. You can tell from the scatterplot above that the relationship is not a strong relationship. The words you choose should match the graphs you create and the statistics you calculate.

Exercise 8(a)

Create a scatterplot of income against race. (Put income on the y-axis and race on the x-axis.)

# Add code here to create a scatterplot of income against race

Exercise 8(b)

Check the three conditions for the relationship between income and race. Which condition(s) are seriously violated here?

Exercise 9(a)

Create a scatterplot of theft against fire. (Put theft on the y-axis and fire on the x-axis.)

# Add code here to create a scatterplot of theft against fire

Exercise 9(b)

Check the three conditions for the relationship between theft and fire. Which condition(s) are seriously violated here?

Exercise 9(c)

Even though the conditions are not met, what if you calculated the correlation coefficient anyway? Try it.

# Add code here to calculate the correlation coefficient between theft and fire

Exercise 9(d)

Suppose you hadn’t looked at the scatterplot and you only saw the correlation coefficient you calculated in the previous part. What would your conclusion be about the relationship between theft and fire. Why would that conclusion be misleading?

Please write up your answer here.

The lesson learned here is that you should never try to interpret a correlation coefficient without looking at a plot of the data to assure that the conditions are met and that the result is a sensible thing to interpret.

6.6 Correlation is not causation

When two variables are correlated—indeed, associated in any way, not just in a linear relationship—that means that there is a relationship between them. However, that does not mean that one variable causes the other variable.

For example, we discovered above that there was a moderate correlation between the racial composition of a ZIP code and the new FAIR policies created in those ZIP codes. However, being part of a racial minority does not cause someone to seek out alternative forms of insurance, at least not directly. In this case, the racial composition of certain neighborhoods, through racist policies, affected the availability of certain forms of insurance for residents in those neighborhoods. And that, in turn, caused residents to seek other forms of insurance.

In the Chicago example, there is still likely a causal connection between one variable (race) and the other (involact), but it is indirect. In other cases, there is no causal connection at all. Here are a few of my favorite examples.

Exercise 10

Ice cream sales are positively correlated with drowning deaths. Does eating ice cream cause you to drown? (Perhaps the myth about swimming within one hour of eating is really true!) Do drowning deaths cause ice cream sales to rise? (Perhaps people are so sad about all the drownings that they have to go out for ice cream to cheer themselves up?)

See if you can figure out the real reason why ice cream sales are positively correlated with drowning deaths.

Please write up your answer here.

---

In the Chicago example, the causal effect was indirect. In the example from the exercise above, there is no causation whatsoever between the two variables. Instead, the causal effect was generated by a third factor that caused both ice cream sales to go up, and also happened to cause drowning deaths to go up. (Or, equivalently stated, it caused ice cream sales to be low during certain times of the year and also caused the drowning deaths to be low as well.) Such a factor is called a lurking variable. When a correlation between two variables exists due solely to the intervention of a lurking variable, that correlation is called a spurious correlation. The correlation is real; a scatterplot of ice cream sales and drowning deaths would show a positive relationship. But the reasons for that correlation to exist have nothing to do with any kind of direct causal link between the two.

Here’s another one:

Exercise 11

Most studies involving children create a number of weird correlations. For example, the height of children is very strongly correlated to pretty much everything you can measure about scholastic aptitude. For example, vocabulary count (the number of words children can use fluently in a sentence) is strongly correlated to height. Are tall people just smarter than short people?

The answer is, of course, no. The correlation is spurious. So what’s the lurking variable?

Please write up your answer here.

6.7 Observational studies versus experiments

So when is a statistical finding (like correlation, for example) evidence of a causal relationship? Before we can answer that question, we need a few more definitions.

A lot of data comes from “observational studies” where we simply observe or measure things as they are “in the wild,” so to speak. We don’t interfere in any way. We just write down what we see. Polls are usually observational in that we ask people questions and record their responses. We do not try to manipulate their responses in any way. We just ask the questions and observe the answers. Field studies are often observational. We go out in nature and write stuff down as we observe it.

Another way to gather data is an experiment. In an experiment, we introduce a manipulation or treatment to try to ascertain its effect. For example. if we’re testing a new drug, we will likely give the drug to one group of patients and a placebo to the other.

Exercise 12

Here’s another internet rabbit hole for you. First, look up the definition of placebo. You do not need to write up your own version of that definition here; just familiarize yourself with the term if you’re not already familiar with it. Next, find some websites about the placebo effect and read those.

Given what you have learned about the placebo effect, why is it important to have a placebo group in a drug trial? Why not just give one set of patients the drug and compare them to another group that takes no pill at all?

Please write up your answer here.

---

The goal of the experiment is to learn whether the treatment (in this example, the drug) is effective when compared to the control (in this example, the placebo).

Note that the word “effective” implies a causal claim. We want to know if the drug causes patients to get better.

Unlike an observational study, in which the relationship between variables can be caused by a lurking variable, in an experiment, we purposefully manipulate one of the variables and try to control all others. For example, we manipulate the drug variable (we purposefully give some people the drug and others the placebo). But we control the amount of the drug given and the schedule on which patients are required to take the pills.

There are lots of things we cannot control. For example, it would be very difficult to control the diet of every person in the experiment. Could diet play a role in whether a patient gets better? Sure, so how do we know diet is not a lurking variable? In the context of an experiment, lurking variables are often called confounders or confounding variables.

One way to mitigate the effect of confounders that we cannot directly control is to randomize the patients into the treatment and control groups. With random selection, there will likely be people who have relatively healthy diets in both the control and treatment groups. If the drugs work, in theory they should still work better for the treatment group than for those taking the placebo. And likewise, patients with less healthy diets will generally be mixed up in both groups, and the drug should also work better for them.

The mantra of experimental design is, “Control as much as you can. Randomize to take care of the rest.”

There are lots of aspects of experimental design that we will not go into here (for example, blinding and blocking). But we will continue to mention the differences between observational studies and experiments in future chapters as we exercise caution in making causal claims.

6.8 Prediction versus explanation

Even when claims are not causal, we can use associations (and correlations more specifically) for purposes of prediction.

Exercise 13

If I tell you that ice cream sales are high right now, can you make a reasonable prediction about the relative number of drowning deaths this month (high or low)? Why or why not?

Please write up your answer here.

---

So even when there is no direct causal link between two variables, if they are positively correlated, then large values of one variable are associated with large values of the other variable. So if I tell you one value is large, it is reasonable to predict that the other value will be large as well.

We use the language “predictor” variable and “response” variable to reinforce this idea.

In a properly designed and controlled experiment, we can use different language. In this case, we can explain the outcome using the treatment variable. If we’ve controlled for everything else, the only possible explanation for a difference between the treatment and control groups must be the treatment variable. If the patients get better on the drug (more so than those on the placebo) and we’ve controlled for every other possible confounding variable, the only possible explanation is that the drug works. The drug “explains” the difference in the response variable.

Be careful, as sometimes statisticians use the term “explanatory variable” to mean any kind of variable that predicts or explains. In this course, we will try to use the term “predictor variable” exclusively.

6.9 Visualizing lurking variables

When we create a scatterplot, we can visualize associations between the two numerical variables. Is there a way to see lurking variables in the scatterplot as well?

One simple case is when the lurking variable is a categorical variable. We saw several examples of that in Chapters 3 and 4 in the penguins data. The association (or lack thereof) between variables was often misleading when we failed to take into account the fact that there were three different species of penguin.

Here are a few more interesting examples. The bdims data (hosted in the openintro package) consists of many body measurements taken from 507 physically active individuals. Type ?bdims at the Console to read the help file.

Here is the data:

bdims

glimpse(bdims)

Rows: 507
Columns: 25
$ bia_di <dbl> 42.9, 43.7, 40.1, 44.3, 42.5, 43.3, 43.5, 44.4, 43.5, 42.0, 40.…
$ bii_di <dbl> 26.0, 28.5, 28.2, 29.9, 29.9, 27.0, 30.0, 29.8, 26.5, 28.0, 29.…
$ bit_di <dbl> 31.5, 33.5, 33.3, 34.0, 34.0, 31.5, 34.0, 33.2, 32.1, 34.0, 33.…
$ che_de <dbl> 17.7, 16.9, 20.9, 18.4, 21.5, 19.6, 21.9, 21.8, 15.5, 22.5, 20.…
$ che_di <dbl> 28.0, 30.8, 31.7, 28.2, 29.4, 31.3, 31.7, 28.8, 27.5, 28.0, 30.…
$ elb_di <dbl> 13.1, 14.0, 13.9, 13.9, 15.2, 14.0, 16.1, 15.1, 14.1, 15.6, 13.…
$ wri_di <dbl> 10.4, 11.8, 10.9, 11.2, 11.6, 11.5, 12.5, 11.9, 11.2, 12.0, 10.…
$ kne_di <dbl> 18.8, 20.6, 19.7, 20.9, 20.7, 18.8, 20.8, 21.0, 18.9, 21.1, 19.…
$ ank_di <dbl> 14.1, 15.1, 14.1, 15.0, 14.9, 13.9, 15.6, 14.6, 13.2, 15.0, 14.…
$ sho_gi <dbl> 106.2, 110.5, 115.1, 104.5, 107.5, 119.8, 123.5, 120.4, 111.0, …
$ che_gi <dbl> 89.5, 97.0, 97.5, 97.0, 97.5, 99.9, 106.9, 102.5, 91.0, 93.5, 9…
$ wai_gi <dbl> 71.5, 79.0, 83.2, 77.8, 80.0, 82.5, 82.0, 76.8, 68.5, 77.5, 81.…
$ nav_gi <dbl> 74.5, 86.5, 82.9, 78.8, 82.5, 80.1, 84.0, 80.5, 69.0, 81.5, 81.…
$ hip_gi <dbl> 93.5, 94.8, 95.0, 94.0, 98.5, 95.3, 101.0, 98.0, 89.5, 99.8, 98…
$ thi_gi <dbl> 51.5, 51.5, 57.3, 53.0, 55.4, 57.5, 60.9, 56.0, 50.0, 59.8, 60.…
$ bic_gi <dbl> 32.5, 34.4, 33.4, 31.0, 32.0, 33.0, 42.4, 34.1, 33.0, 36.5, 34.…
$ for_gi <dbl> 26.0, 28.0, 28.8, 26.2, 28.4, 28.0, 32.3, 28.0, 26.0, 29.2, 27.…
$ kne_gi <dbl> 34.5, 36.5, 37.0, 37.0, 37.7, 36.6, 40.1, 39.2, 35.5, 38.3, 38.…
$ cal_gi <dbl> 36.5, 37.5, 37.3, 34.8, 38.6, 36.1, 40.3, 36.7, 35.0, 38.6, 40.…
$ ank_gi <dbl> 23.5, 24.5, 21.9, 23.0, 24.4, 23.5, 23.6, 22.5, 22.0, 22.2, 23.…
$ wri_gi <dbl> 16.5, 17.0, 16.9, 16.6, 18.0, 16.9, 18.8, 18.0, 16.5, 16.9, 16.…
$ age    <int> 21, 23, 28, 23, 22, 21, 26, 27, 23, 21, 23, 22, 20, 26, 23, 22,…
$ wgt    <dbl> 65.6, 71.8, 80.7, 72.6, 78.8, 74.8, 86.4, 78.4, 62.0, 81.6, 76.…
$ hgt    <dbl> 174.0, 175.3, 193.5, 186.5, 187.2, 181.5, 184.0, 184.5, 175.0, …
$ sex    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …

Most physical body measurements are known to be correlated; this makes sense because when one part of the body is larger, we expect lots of other body parts to be larger as well (and similarly for smaller individuals).

For example, it’s no surprise that shoulder girth (sho_gi) and chest girth (che_gi) are strongly correlated:

ggplot(bdims, aes(y = sho_gi, x = che_gi)) +
  geom_point()

Is there a possible lurking variable here, though? You may wonder about sex. (In this data set, the sex variable is presumed to be biological sex assigned at birth.)

Before we go any further, go back to the help file and the glimpse output above and note that sex is coded as an integer (a whole number). We’ll use the mutate and as_factor commands—illustrated in Chapters 3 and 5—to make a new factor variable.

bdims <- bdims |>
  mutate(sex_fct = as_factor(sex))
glimpse(bdims)

Rows: 507
Columns: 26
$ bia_di  <dbl> 42.9, 43.7, 40.1, 44.3, 42.5, 43.3, 43.5, 44.4, 43.5, 42.0, 40…
$ bii_di  <dbl> 26.0, 28.5, 28.2, 29.9, 29.9, 27.0, 30.0, 29.8, 26.5, 28.0, 29…
$ bit_di  <dbl> 31.5, 33.5, 33.3, 34.0, 34.0, 31.5, 34.0, 33.2, 32.1, 34.0, 33…
$ che_de  <dbl> 17.7, 16.9, 20.9, 18.4, 21.5, 19.6, 21.9, 21.8, 15.5, 22.5, 20…
$ che_di  <dbl> 28.0, 30.8, 31.7, 28.2, 29.4, 31.3, 31.7, 28.8, 27.5, 28.0, 30…
$ elb_di  <dbl> 13.1, 14.0, 13.9, 13.9, 15.2, 14.0, 16.1, 15.1, 14.1, 15.6, 13…
$ wri_di  <dbl> 10.4, 11.8, 10.9, 11.2, 11.6, 11.5, 12.5, 11.9, 11.2, 12.0, 10…
$ kne_di  <dbl> 18.8, 20.6, 19.7, 20.9, 20.7, 18.8, 20.8, 21.0, 18.9, 21.1, 19…
$ ank_di  <dbl> 14.1, 15.1, 14.1, 15.0, 14.9, 13.9, 15.6, 14.6, 13.2, 15.0, 14…
$ sho_gi  <dbl> 106.2, 110.5, 115.1, 104.5, 107.5, 119.8, 123.5, 120.4, 111.0,…
$ che_gi  <dbl> 89.5, 97.0, 97.5, 97.0, 97.5, 99.9, 106.9, 102.5, 91.0, 93.5, …
$ wai_gi  <dbl> 71.5, 79.0, 83.2, 77.8, 80.0, 82.5, 82.0, 76.8, 68.5, 77.5, 81…
$ nav_gi  <dbl> 74.5, 86.5, 82.9, 78.8, 82.5, 80.1, 84.0, 80.5, 69.0, 81.5, 81…
$ hip_gi  <dbl> 93.5, 94.8, 95.0, 94.0, 98.5, 95.3, 101.0, 98.0, 89.5, 99.8, 9…
$ thi_gi  <dbl> 51.5, 51.5, 57.3, 53.0, 55.4, 57.5, 60.9, 56.0, 50.0, 59.8, 60…
$ bic_gi  <dbl> 32.5, 34.4, 33.4, 31.0, 32.0, 33.0, 42.4, 34.1, 33.0, 36.5, 34…
$ for_gi  <dbl> 26.0, 28.0, 28.8, 26.2, 28.4, 28.0, 32.3, 28.0, 26.0, 29.2, 27…
$ kne_gi  <dbl> 34.5, 36.5, 37.0, 37.0, 37.7, 36.6, 40.1, 39.2, 35.5, 38.3, 38…
$ cal_gi  <dbl> 36.5, 37.5, 37.3, 34.8, 38.6, 36.1, 40.3, 36.7, 35.0, 38.6, 40…
$ ank_gi  <dbl> 23.5, 24.5, 21.9, 23.0, 24.4, 23.5, 23.6, 22.5, 22.0, 22.2, 23…
$ wri_gi  <dbl> 16.5, 17.0, 16.9, 16.6, 18.0, 16.9, 18.8, 18.0, 16.5, 16.9, 16…
$ age     <int> 21, 23, 28, 23, 22, 21, 26, 27, 23, 21, 23, 22, 20, 26, 23, 22…
$ wgt     <dbl> 65.6, 71.8, 80.7, 72.6, 78.8, 74.8, 86.4, 78.4, 62.0, 81.6, 76…
$ hgt     <dbl> 174.0, 175.3, 193.5, 186.5, 187.2, 181.5, 184.0, 184.5, 175.0,…
$ sex     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ sex_fct <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

If you look at the glimpse output above, you see that we do have a new variable called sex_fct and it is properly coded as a factor variable. However, the labels 0 and 1 (for females and males, respectively) are not very helpful. Can we change them? Yes, the forcats package has a fct_recode function that does just that. Here is what it looks like:

bdims <- bdims |>
  mutate(sex_fct = fct_recode(sex_fct, "female" = "0", "male" = "1"))
glimpse(bdims)

Rows: 507
Columns: 26
$ bia_di  <dbl> 42.9, 43.7, 40.1, 44.3, 42.5, 43.3, 43.5, 44.4, 43.5, 42.0, 40…
$ bii_di  <dbl> 26.0, 28.5, 28.2, 29.9, 29.9, 27.0, 30.0, 29.8, 26.5, 28.0, 29…
$ bit_di  <dbl> 31.5, 33.5, 33.3, 34.0, 34.0, 31.5, 34.0, 33.2, 32.1, 34.0, 33…
$ che_de  <dbl> 17.7, 16.9, 20.9, 18.4, 21.5, 19.6, 21.9, 21.8, 15.5, 22.5, 20…
$ che_di  <dbl> 28.0, 30.8, 31.7, 28.2, 29.4, 31.3, 31.7, 28.8, 27.5, 28.0, 30…
$ elb_di  <dbl> 13.1, 14.0, 13.9, 13.9, 15.2, 14.0, 16.1, 15.1, 14.1, 15.6, 13…
$ wri_di  <dbl> 10.4, 11.8, 10.9, 11.2, 11.6, 11.5, 12.5, 11.9, 11.2, 12.0, 10…
$ kne_di  <dbl> 18.8, 20.6, 19.7, 20.9, 20.7, 18.8, 20.8, 21.0, 18.9, 21.1, 19…
$ ank_di  <dbl> 14.1, 15.1, 14.1, 15.0, 14.9, 13.9, 15.6, 14.6, 13.2, 15.0, 14…
$ sho_gi  <dbl> 106.2, 110.5, 115.1, 104.5, 107.5, 119.8, 123.5, 120.4, 111.0,…
$ che_gi  <dbl> 89.5, 97.0, 97.5, 97.0, 97.5, 99.9, 106.9, 102.5, 91.0, 93.5, …
$ wai_gi  <dbl> 71.5, 79.0, 83.2, 77.8, 80.0, 82.5, 82.0, 76.8, 68.5, 77.5, 81…
$ nav_gi  <dbl> 74.5, 86.5, 82.9, 78.8, 82.5, 80.1, 84.0, 80.5, 69.0, 81.5, 81…
$ hip_gi  <dbl> 93.5, 94.8, 95.0, 94.0, 98.5, 95.3, 101.0, 98.0, 89.5, 99.8, 9…
$ thi_gi  <dbl> 51.5, 51.5, 57.3, 53.0, 55.4, 57.5, 60.9, 56.0, 50.0, 59.8, 60…
$ bic_gi  <dbl> 32.5, 34.4, 33.4, 31.0, 32.0, 33.0, 42.4, 34.1, 33.0, 36.5, 34…
$ for_gi  <dbl> 26.0, 28.0, 28.8, 26.2, 28.4, 28.0, 32.3, 28.0, 26.0, 29.2, 27…
$ kne_gi  <dbl> 34.5, 36.5, 37.0, 37.0, 37.7, 36.6, 40.1, 39.2, 35.5, 38.3, 38…
$ cal_gi  <dbl> 36.5, 37.5, 37.3, 34.8, 38.6, 36.1, 40.3, 36.7, 35.0, 38.6, 40…
$ ank_gi  <dbl> 23.5, 24.5, 21.9, 23.0, 24.4, 23.5, 23.6, 22.5, 22.0, 22.2, 23…
$ wri_gi  <dbl> 16.5, 17.0, 16.9, 16.6, 18.0, 16.9, 18.8, 18.0, 16.5, 16.9, 16…
$ age     <int> 21, 23, 28, 23, 22, 21, 26, 27, 23, 21, 23, 22, 20, 26, 23, 22…
$ wgt     <dbl> 65.6, 71.8, 80.7, 72.6, 78.8, 74.8, 86.4, 78.4, 62.0, 81.6, 76…
$ hgt     <dbl> 174.0, 175.3, 193.5, 186.5, 187.2, 181.5, 184.0, 184.5, 175.0,…
$ sex     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ sex_fct <fct> male, male, male, male, male, male, male, male, male, male, ma…

This will be a lot more helpful!

Now, back to the scatterplots.

One way we learned (in Chapters 3 and 4) to incorporate a third variable into the analysis is through the use of color as an additional aesthetic element. We’ll use our new sex_fct variable. Also, don’t forget to use the Viridis color palette and the black-and-white theme.

ggplot(bdims, aes(y = sho_gi, x = che_gi, color = sex_fct)) +
  geom_point() +
  scale_color_viridis_d() +
  theme_bw()

In this example, there is a strong correlation between shoulder girth and chest girth, but females and males lie in completely different parts of the graph. Having said that, if you focus on the females separately, you can still see a strong positive correlation, and if you focus on males separately, there is also a strong positive correlation there. So the inclusion of sex didn’t really change much about the nature of the correlation in this example. Even still, the correlation coefficients do change a little depending on whether we look at the whole data set versus females/males separately:

cor(bdims$sho_gi, bdims$che_gi)

[1] 0.9271923

bdims |>
  group_by(sex_fct) |>
  summarise(corr = cor(sho_gi, che_gi))

Exercise 14

Why would the correlation coefficient be stronger for the whole data set and slightly less strong for the sexes separately? (Hint: think about sample size.)

Please write up your answer here.

---

In the previous example, sex was a lurking variable, but it did not radically alter the nature of the association. What about the examples in these next two sets of exercises?

Exercise 15(a)

Create a scatterplot of thigh girth against weight (put thi_gi on the y-axis and wgt on the x-axis).

# Add code here to create a scatterplot of thigh girth against weight.

Exercise 15(b)

Change the scatterplot above to include sex_fct as a color aesthetic. (Use the Viridis color palette and theme_bw.)

# Add code here to add color for sex_fct.

Exercise 15(c)

Calculate the correlation coefficients for thigh girth and weight, once for the whole data set, and again and for the data split by sex_fct (as above).

# Add code here to calculate the correlation coefficient
# between thigh girth and weight.

# Add code here to calculate the correlation coefficient
# between thigh girth and weight split by sex.

Exercise 15(d)

Explain how sex is a lurking variable here. In other words, how did ignoring/considering sex alter the way we perceived the correlation between thigh girth and weight? What changed about the nature of the correlation within each sex category?

Please write up your answer here.

Exercise 16(a)

The help file for the bia_di variable describes it as the “respondent’s biacromial diameter in centimeters.” What is “biacromial diameter”?

Please write up your answer here.

Exercise 16(b)

Create a scatterplot of biacromial diameter against weight (put bia_di on the y-axis and wgt on the x-axis).

# Add code here to create a scatterplot of biacromial diameter against weight.

Exercise 16(c)

Change the scatterplot above to include sex_fct as a color aesthetic. (Use the Viridis color palette and theme_bw.)

# Add code here to add color for sex_fct.

Exercise 16(d)

Calculate the correlation coefficients for biacromial diameter and weight, once for the whole data set, and again and for the data split by sex_fct (as above).

# Add code here to calculate the correlation coefficient
# between biacromial diameter and weight

# Add code here to calculate the correlation coefficient
# between biacromial diameter and weight split by sex

Exercise 16(e)

Explain how sex is a lurking variable here. In other words, how did ignoring/considering sex alter the way we perceived the correlation between biacromial diameter and weight? What changed about the nature of the correlation within each sex category?

Please write up your answer here.

---

The take-home message here is that lurking variables can change the strength of the correlation between two variables, making it appear stronger or weaker. In more extreme cases, it’s even possible to change the direction of the correlation altogether! There isn’t an example of this phenomenon in the bdims data, but we do find one in the penguins data.

Here is a scatterplot of bill depth against bill length.

ggplot(penguins, aes(y = bill_depth_mm, x = bill_length_mm)) +
  geom_point()

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

There is not much correlation between bill depth and bill length, but if anything, it looks like there might be a slightly negative association. (In the following code chunk, the cor command uses a different method for dealing with missing data.)

cor(penguins$bill_depth_mm, penguins$bill_length_mm,
    use = "complete.obs")

[1] -0.2350529

Now split by species:

ggplot(penguins, aes(y = bill_depth_mm, x = bill_length_mm,
                     color = species)) +
  geom_point() +
  scale_color_viridis_d() +
  theme_bw()

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

penguins |>
  group_by(species) |>
  summarise(corr = cor(bill_depth_mm, bill_length_mm,
                       use = "complete.obs"))

There was a very weak negative correlation in the full data set, but, behold, bill depth and bill length are positive correlated within each species!

The phenomenon of an association between two variables “reversing” direction when considering a third variable is often called “Simpson’s Paradox”.² We’ll revisit Simpson’s Paradox in a future chapter.

6.10 Conclusion

If we have two numerical variables that have a linear association between them (also assuming there are no serious outliers), we can compute the correlation coefficient that measures the strength and direction of that linear association.

Keep in mind that in an observational study, the correlation coefficient is a measure of association but it does not signify that one variable causes the other. It’s possible that one variable causes the other, but it’s also possible that a third “lurking” variable is responsible for all or part of the association. Either way, the fact that a relationship exists means it is possible to use values of one variable to make reasonable predictions about the values of the other variable.

In a properly designed experiment, the manipulation of one variable while controlling for others (and randomizing to take care of other confounders) ensures that there is a causal link between the treatment variable and the response of interest. In this case, the treatment can “explain” the response, not just predict it.

Finally, when there is a lurking variable (confounding), we can visualize it by adding it to our scatterplots as a color aesthetic. There, we often find that associations change; they can get stronger, weaker, or even change direction altogether in the presence of that third variable.

6.10.1 Preparing and submitting your assignment

1. From the “Run” menu, select “Restart R and Run All Chunks”.
2. Deal with any code errors that crop up. Repeat steps 1–2 until there are no more code errors.
3. Spell check your document by clicking the icon with “ABC” and a check mark.
4. Hit the “Render” button one last time to generate the final draft of the HTML file. (If there are errors here, you may need to go back and fix broken inline code or other markdown issues.)
5. Proofread the HTML file carefully. If there are errors, go back and fix them, then repeat steps 1–5 again.

If you have completed this chapter as part of a statistics course, follow the directions you receive from your professor to submit your assignment.

---

1. There are other ways of measuring association for variables that are not numerical, but these aren’t covered in this course.

2. Just like for every other named concept in statistics, Simpson’s Paradox wasn’t first observed by Simpson.