Related topics

Section 5.3 Related topics

Subsection 5.3.1 Correlation is not necessarily causation

Throughout this chapter we have been cautious when interpreting regression slope coefficients. We always discussed the “associated” effect of an explanatory variable \(x\) on an outcome variable \(y\text{.}\) For example, our statement from Subsection 5.1.2 that “for every increase of 1 unit in life_exp, there is an associated decrease of on average some number of units of fert_rate.” We include the term “associated” to be extra careful not to suggest we are making a causal statement. So while life_exp is negatively correlated with fert_rate, we can’t necessarily make any statements about life expectancy’s direct causal effect on fertility rates without more information.

Here is another example depicted in Figure 5.3.1: a not-so-great medical doctor goes through medical records and finds that patients who slept with their shoes on tended to wake up more with headaches. So this doctor declares, “Sleeping with shoes on causes headaches!”

Image illustrating the spurious correlation between sleeping with shoes on and waking up with a headache. — Figure 5.3.1. Does sleeping with shoes on cause headaches?

However, there is a good chance that if someone is sleeping with their shoes on, it is potentially because they are intoxicated from alcohol. Furthermore, higher levels of drinking leads to more hangovers, and hence more headaches. The amount of alcohol consumption here is what’s known as a confounding/lurking variable. It “lurks” behind the scenes, confounding the causal relationship (if any) of “sleeping with shoes on” with “waking up with a headache.” We can summarize this in Figure 5.3.2 with a causal graph where Y is a response variable and X is a treatment variable whose causal effect we are interested in.

Causal graph showing X (treatment), Y (response), and Z (confounding variable) with arrows indicating relationships. — Figure 5.3.2. Causal graph.

To study the relationship between Y and X, we could use a regression model where the outcome variable is set to Y and the explanatory variable is set to be X, as you’ve been doing throughout this chapter. However, Figure 5.3.2 also includes a third variable with arrows pointing at both X and Y:

Z is a confounding variable that affects both X and Y, thereby “confounding” their relationship. Here the confounding variable is alcohol.

Alcohol will cause people to be both more likely to sleep with their shoes on as well as be more likely to wake up with a headache. Thus any regression model of the relationship between X and Y should also use Z as an explanatory variable. In other words, our doctor needs to take into account who had been drinking the night before. In the next chapter, we’ll start covering multiple regression models that allow us to incorporate more than one variable in our regression models.

Establishing causation is a tricky problem and frequently takes either carefully designed experiments or methods to control for the effects of confounding variables. Both these approaches attempt, as best they can, either to take all possible confounding variables into account or negate their impact. This allows researchers to focus only on the relationship of interest: the relationship between the outcome variable Y and the treatment variable X.

As you read news stories, be careful not to fall into the trap of thinking that correlation necessarily implies causation. Check out the Spurious Correlations website for some comical examples of variables that are correlated, but are not causally related.

Coming back to our UN member states data, a confounding variable could be the level of economic development of a country. This could affect both life expectancy and fertility rates. A proxy for looking at the level of economic development could be the Human Development Index (HDI). It measures a country’s average achievements in three basic aspects of human development: health (life expectancy), education (mean and expected years of schooling), and standard of living (gross national income per capita). This is stored in the hdi_2022 column of un_member_states_2024. We explore its relationship with life expectancy and fertility rates:

ggplot(data = un_member_states_2024, 
       aes(x = hdi_2022, y = life_expectancy_2022)) +
  geom_point() +
  labs(x = "Human Development Index (HDI)", y = "Life Expectancy")

Scatterplot showing a positive relationship between the Human Development Index (HDI) and life expectancy across UN member states. — Figure 5.3.3. Scatterplot of Human Development Index vs Life Expectancy.

ggplot(data = un_member_states_2024, 
       aes(x = hdi_2022, y = fertility_rate_2022)) +
  geom_point() +
  labs(x = "Human Development Index (HDI)", y = "Fertility Rate")

Scatterplot showing a negative relationship between the Human Development Index (HDI) and fertility rate across UN member states. — Figure 5.3.4. Scatterplot of Human Development Index vs Fertility Rate.

un_member_states_2024 |> 

  get_correlation(life_expectancy_2022 ~ hdi_2022, na.rm = TRUE)

# A tibble: 1 × 1
    cor
  <dbl>
1 0.889

un_member_states_2024 |> 

  get_correlation(fertility_rate_2022 ~ hdi_2022, na.rm = TRUE)

# A tibble: 1 × 1
     cor
   <dbl>
1 -0.849

Looking at both of the scatterplots above, we see a strong positive linear relationship between the Human Development Index (HDI) and life expectancy, as well as a strong negative correlation between the Human Development Index (HDI) and fertility rates. The correlation coefficients between the Human Development Index (HDI) and life expectancy, as well as the Human Development Index (HDI) and fertility rates, are also given. These findings as well as some additional understanding of socioeconomic factors suggest that the Human Development Index (HDI) is a confounding variable in the relationship between life expectancy and fertility rates.

Subsection 5.3.2 Best-fitting line

Regression lines are also known as “best-fitting” lines. But what do we mean by “best”? We unpack the criteria that is used in regression to determine “best.” Recall Figure 5.1.5, where for Bosnia and Herzegovina we marked the observed value \(y\) with a circle, the fitted value \(\widehat{y}\) with a square, and the residual \(y - \widehat{y}\) with an arrow. We re-display Figure 5.1.5 alongside three more arbitrarily chosen countries (Chad, India, and the Solomon Islands), each with their own residuals marked:

Four scatterplots progressively adding residuals for Bosnia and Herzegovina, Chad, India, and the Solomon Islands. — Figure 5.3.5. Example of observed value, fitted value, and residual for four countries.

The four plots refer to:

The country of Bosnia and Herzegovina had a life expectancy \(x = 77.98\) and fertility rate \(y = 1.3\text{.}\) The residual in this case is \(1.3 - 1.894 = -0.594\text{,}\) which we mark with a blue arrow in the top-left plot.
The addition of Chad, which had a life expectancy \(x = 59.15\) and fertility rate \(y = 6\text{.}\) The residual in this case is \(6 - 4.479 = 1.521\text{,}\) which we mark with a new blue arrow in the top-right plot.
The addition of India, which had a life expectancy \(x = 67.22\) and fertility rate \(y = 2\text{.}\) The residual in this case is \(2 - 3.371 = -1.371\text{,}\) which we mark with a new blue arrow in the bottom-left plot.
The addition of the Solomon Islands, which had a life expectancy \(x = 76.7\) and fertility rate \(y = 3.8\text{.}\) The residual in this case is \(3.8 - 2.069 = 1.731\text{,}\) which we mark with a new blue arrow in the bottom-right plot.

Now say we repeated this process of computing residuals for all 181 countries with complete information, then we squared all the residuals, and then we summed them. We call this quantity the sum of squared residuals; it is a measure of the lack of fit of a model. Larger values of the sum of squared residuals indicate a bigger lack of fit. This corresponds to a worse-fitting model.

If the regression line fits all the points perfectly, then the sum of squared residuals is 0. This is because if the regression line fits all the points perfectly, then the fitted value \(\widehat{y}\) equals the observed value \(y\) in all cases, and hence the residual \(y - \widehat{y} = 0\) in all cases, and the sum of even a large number of 0’s is still 0.

Furthermore, of all possible lines we can draw through the cloud of points, the regression line minimizes this value. In other words, the regression and its corresponding fitted values \(\widehat{y}\) minimizes the sum of the squared residuals:

\begin{equation*} \sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 \end{equation*}

We use our data-wrangling tools to compute the sum of squared residuals exactly:

# Fit regression model and regression points
demographics_model <- lm(fert_rate ~ life_exp, data = UN_data_ch5)
regression_points <- get_regression_points(demographics_model)

# Compute sum of squared residuals
regression_points |>
  mutate(squared_residuals = residual^2) |>
  summarize(sum_of_squared_residuals = sum(squared_residuals))

[1] 81.3

Any other straight line drawn in the figure would yield a sum of squared residuals greater than this value. This is a mathematically guaranteed fact that you can prove using calculus and linear algebra. That’s why alternative names for the linear regression line are the best-fitting line and the least-squares line.

Why do we square the residuals (i.e., the arrow lengths)? So that both positive and negative deviations of the same amount are treated equally. (That being said, while taking the absolute value of the residuals would also treat both positive and negative deviations of the same amount equally, squaring the residuals is used for reasons related to calculus: taking derivatives and minimizing functions. To learn more we suggest you consult a textbook on mathematical statistics.)

Exercises Exercises

1.

Note in Figure 5.3.6 there are 3 points marked with dots and:

The “best” fitting solid regression line in blue
An arbitrarily chosen dotted red line
Another arbitrarily chosen dashed green line

Scatterplot with three points and three lines: the best-fit regression line (blue solid), an arbitrary red dotted line, and a green dashed line. — Figure 5.3.6. Regression line and two others.

Compute the sum of squared residuals by hand for each line. Show that the regression line in blue has the smallest value of these three lines.

Answer.

The regression line in blue:

\begin{equation*} \sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-1.5)^2+(1.0-2.0)^2+(3.0-2.5)^2 = 1.5 \end{equation*}

The arbitrarily chosen dotted red line (horizontal at \(y = 2.5\)):

\begin{equation*} \sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-2.5)^2+(1.0-2.5)^2+(3.0-2.5)^2 = 2.75 \end{equation*}

The arbitrarily chosen dashed green line (\(y = 2 - x\)):

\begin{equation*} \sum_{i=1}^{n}(y_i - \widehat{y}_i)^2 = (2.0-2.0)^2+(1.0-1.5)^2+(3.0-1.0)^2 = 4.25 \end{equation*}

As calculated, \(1.5 \lt 2.75 \lt 4.25\text{.}\) Therefore, we show that the regression line in blue has the smallest value of the residual sum of squares.

Subsection 5.3.3 get_regression_x() functions

Recall in this chapter we introduced a wrapper function from the moderndive package:

get_regression_points() that returns point-by-point information from a regression model in Subsection 5.1.3.

What is going on behind the scenes with the get_regression_points() function? We mentioned in Subsection 5.1.2 that this was an example of a wrapper function. Such functions take other pre-existing functions and “wrap” them into single functions that hide the user from their inner workings. This way all the user needs to worry about is what the inputs look like and what the outputs look like. In this subsection, we’ll “get under the hood” of these functions and see how the “engine” of these wrapper functions works.

The get_regression_points() function is a wrapper function, returning information about the individual points involved in a regression model like the fitted values, observed values, and the residuals. get_regression_points() uses the augment() function in the broom package to produce the data shown in Table 5.3.7. Additionally, it uses clean_names() from the janitor package [23] to clean up the variable names.

library(broom)
library(janitor)

demographics_model |>
  augment() |>
  mutate_if(is.numeric, round, digits = 3) |>
  clean_names() |>
  select(-c("std_resid", "hat", "sigma", "cooksd", "std_resid"))

Table 5.3.7. Regression points using augment() from broom package

fert_rate	life_exp	fitted	resid
4.3	53.6	5.23	-0.934
1.4	79.5	1.69	-0.289
2.7	78.0	1.89	0.813
5.0	62.1	4.07	0.928
1.6	77.8	1.92	-0.318
1.9	78.3	1.85	0.052
1.6	76.1	2.15	-0.548
1.6	83.1	1.19	0.408
1.5	82.3	1.30	0.195
1.6	74.2	2.42	-0.819

In this case, it outputs only the variables of interest to students learning regression: the outcome variable \(y\) (fert_rate), all explanatory/predictor variables (life_exp), all resulting fitted values \(\hat{y}\) used by applying the equation of the regression line to life_exp, and the residual \(y - \hat{y}\text{.}\)

If you are even more curious about how these and other wrapper functions work, take a look at the source code for these functions on GitHub.

Prev Top Next