Data Source

We’ll be using the WHO’s life expectancy dataset, found on Kaggle here.

Libraries

library(tidyverse)
library(skimr)
library(broom)

theme_set(theme_bw())

Data Import and Tidy

Reading in and having a quick look at the dataset:

lifexp_df <- 
  read.csv(file = "./data/who_life_expectancy.csv") %>% 
  janitor::clean_names() 

glimpse(lifexp_df)

## Rows: 2,938
## Columns: 22
## $ country                         <chr> "Afghanistan", "Afghanistan", "Afghani…
## $ year                            <int> 2015, 2014, 2013, 2012, 2011, 2010, 20…
## $ status                          <chr> "Developing", "Developing", "Developin…
## $ life_expectancy                 <dbl> 65.0, 59.9, 59.9, 59.5, 59.2, 58.8, 58…
## $ adult_mortality                 <int> 263, 271, 268, 272, 275, 279, 281, 287…
## $ infant_deaths                   <int> 62, 64, 66, 69, 71, 74, 77, 80, 82, 84…
## $ alcohol                         <dbl> 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.…
## $ percentage_expenditure          <dbl> 71.279624, 73.523582, 73.219243, 78.18…
## $ hepatitis_b                     <int> 65, 62, 64, 67, 68, 66, 63, 64, 63, 64…
## $ measles                         <int> 1154, 492, 430, 2787, 3013, 1989, 2861…
## $ bmi                             <dbl> 19.1, 18.6, 18.1, 17.6, 17.2, 16.7, 16…
## $ under_five_deaths               <int> 83, 86, 89, 93, 97, 102, 106, 110, 113…
## $ polio                           <int> 6, 58, 62, 67, 68, 66, 63, 64, 63, 58,…
## $ total_expenditure               <dbl> 8.16, 8.18, 8.13, 8.52, 7.87, 9.20, 9.…
## $ diphtheria                      <int> 65, 62, 64, 67, 68, 66, 63, 64, 63, 58…
## $ hiv_aids                        <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1…
## $ gdp                             <dbl> 584.25921, 612.69651, 631.74498, 669.9…
## $ population                      <dbl> 33736494, 327582, 31731688, 3696958, 2…
## $ thinness_1_19_years             <dbl> 17.2, 17.5, 17.7, 17.9, 18.2, 18.4, 18…
## $ thinness_5_9_years              <dbl> 17.3, 17.5, 17.7, 18.0, 18.2, 18.4, 18…
## $ income_composition_of_resources <dbl> 0.479, 0.476, 0.470, 0.463, 0.454, 0.4…
## $ schooling                       <dbl> 10.1, 10.0, 9.9, 9.8, 9.5, 9.2, 8.9, 8…

skim(lifexp_df)

Variable type: character

Variable type: numeric

The dataset contains 2938 rows, spanning 22 variables. Every row represents a country and year combination - a total of 193 unique countries for every year from 2000 to 2015. We’ll limit the analysis to one year since having the same country repeated as a new row constitutes a repeated measure, which requires more sophisticated analysis than what we’re doing here. I’m going to (arbitrarily) choose 2012.

We also notice from the glimpse() output that we have quite a few missing values, so we’ll need to drop them from our analysis.

Of the 22 variables available, we’re really only going to focus on the following variables: * life_expectancy: average life expectancy at birth, measured in years * schooling: national average of years of formal education * status: binary variable coded “Developed” or “Developing” (this will allow for some interesting stratification later on)

We will ignore the other variables for this analysis. Thus, our final dataset is as follows:

final_df <-
  lifexp_df %>% 
  filter(year == "2012") %>% 
  drop_na(life_expectancy, schooling, status)

Step 1: Fit and Interpret Model(s)

final_df %>% 
  ggplot(aes(x = schooling, y = life_expectancy)) + 
  geom_point() +
  labs(
    title = "Figure 1: Life Expectancy vs Schooling"
  )

model_1 <-
  lm(data = final_df, life_expectancy ~ schooling)

summary(model_1)

## 
## Call:
## lm(formula = life_expectancy ~ schooling, data = final_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8459  -2.5066   0.3488   3.2111  12.4766 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  41.8211     1.7213   24.30   <2e-16 ***
## schooling     2.2932     0.1319   17.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.003 on 171 degrees of freedom
## Multiple R-squared:  0.6386, Adjusted R-squared:  0.6365 
## F-statistic: 302.1 on 1 and 171 DF,  p-value: < 2.2e-16

final_df %>% 
  ggplot(aes(x = schooling, y = life_expectancy)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Figure 2: Fitted Model"
  )

model_1_df = 
  broom::augment(model_1)
model_1_df

## # A tibble: 173 × 8
##    life_expectancy schooling .fitted .resid    .hat .sigma   .cooksd .std.resid
##              <dbl>     <dbl>   <dbl>  <dbl>   <dbl>  <dbl>     <dbl>      <dbl>
##  1            59.5       9.8    64.3 -4.79  0.0117    5.00 0.00551      -0.964 
##  2            76.9      14.2    74.4  2.52  0.00729   5.01 0.000936      0.505 
##  3            75.1      14.4    74.8  0.257 0.00773   5.02 0.0000104     0.0516
##  4            56        10.3    65.4 -9.44  0.00987   4.96 0.0179       -1.90  
##  5            75.9      13.8    73.5  2.43  0.00658   5.01 0.000789      0.488 
##  6            75.9      17.2    81.3 -5.36  0.0197    5.00 0.0118       -1.08  
##  7            74.4      12.7    70.9  3.46  0.00578   5.01 0.00139       0.693 
##  8            82.3      20.1    87.9 -5.61  0.0436    5.00 0.0300       -1.15  
##  9            88        15.7    77.8 10.2   0.0119    4.96 0.0253        2.05  
## 10            71.9      11.8    68.9  3.02  0.00637   5.01 0.00118       0.605 
## # … with 163 more rows

model_1_df %>% 
  ggplot(aes(x = schooling, y = life_expectancy)) + 
  geom_point() +
  geom_segment(aes(xend = schooling, yend = .fitted), color = "red") +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Figure 3: Model Residuals"
  )

final_df %>% 
  mutate(index = row_number()) %>% 
  ggplot(aes(x = index, y = life_expectancy, group = status)) + 
  geom_point(aes(color = status)) + 
  labs(
    title = "Figure 4: Life Expectancy Difference by Status "
  )

model_2 = 
  lm(data = final_df, life_expectancy ~ schooling + status)
summary(model_2)

## 
## Call:
## lm(formula = life_expectancy ~ schooling + status, data = final_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7545  -3.1520   0.5462   3.5773  12.4017 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       45.4937     2.7312  16.657   <2e-16 ***
## schooling          2.1420     0.1578  13.578   <2e-16 ***
## statusDeveloping  -2.1010     1.2176  -1.725   0.0863 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.975 on 170 degrees of freedom
## Multiple R-squared:  0.6448, Adjusted R-squared:  0.6406 
## F-statistic: 154.3 on 2 and 170 DF,  p-value: < 2.2e-16

model_3 = lm(data = final_df, life_expectancy ~ status)
summary(model_3)

## 
## Call:
## lm(formula = life_expectancy ~ status, data = final_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.408  -5.408   1.307   5.592  14.892 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        80.393      1.330  60.456  < 2e-16 ***
## statusDeveloping  -11.285      1.458  -7.742 8.17e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.161 on 171 degrees of freedom
## Multiple R-squared:  0.2596, Adjusted R-squared:  0.2552 
## F-statistic: 59.94 on 1 and 171 DF,  p-value: 8.169e-13

final_df %>% 
  ggplot(aes(x = status, y = life_expectancy)) + 
  geom_point() + 
  stat_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Figure 5: Status vs Life Expectancy"
  )

final_df %>% 
  ggplot(aes(x = schooling, y = life_expectancy, group = status)) + 
  geom_point(aes(color = status))

Step 2: Figure out how “good” your model is

Smith et al (2017) describe regression diagnostics as an “art”, referring to the absence of a structured process for model evaluation. After looking through a digital bushel of papers, I would have to agree. If you spend any time looking at scientific work employing linear regression models, you’ll be hard-pressed to find any consistency in the measures scientists use to evaluate their models (if they use such measures at all).

However, it’s clear that simply fitting a linear model is not enough. While there are many different constructs and measures to evaluate model performance, I would like to, at minimum, understand the magnitude of a model’s error, its predictive capacity, the impact of outliers and influential observations, and whether it meets the assumptions for linear modeling in the first place. We’ll look at this last topic in the next section.

First, let’s look at perhaps the most common model diagnostic, $R^2$, or, more formally, the coefficient of determination. For a single-variable model, $R^2$, is just the square of our friend, the Pearson correlation coefficient. It turns out that this is actually a much more useful quantity:

$$R^2=\frac{{\Sigma }_{i=1}^n(y_i-{\hat{y}}_i{)}^2}{{\Sigma }_{i=1}^n(y_i-{\overline{y}}_i{)}^2}$$

$R^2$ is powerful because it’s completely intuitive - it equals the percentage of variance in the outcome explained by the predictor(s). This is part of the summary() output:

summary(model_1)

## 
## Call:
## lm(formula = life_expectancy ~ schooling, data = final_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8459  -2.5066   0.3488   3.2111  12.4766 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  41.8211     1.7213   24.30   <2e-16 ***
## schooling     2.2932     0.1319   17.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.003 on 171 degrees of freedom
## Multiple R-squared:  0.6386, Adjusted R-squared:  0.6365 
## F-statistic: 302.1 on 1 and 171 DF,  p-value: < 2.2e-16

For the single-variable model, $R^2=0.6386$, so 63.86% of the variance in a nation’s life expectancy can be explained by its linear relationship to education. For a single variable model, this is a pretty high $R^2$. Note that when we fit model_2, the $R^2$ only went up to 64.48%, all while taking away a degree of freedom from the model - further evidence that status is not a worthwhile predictor in this context.

To quantify the model’s error, which, again, comes down to residuals, we can look at Root Mean Square Error (RMSE), another highly common regression diagnostic. It is defined as:

$$RMSE=\sqrt{\frac{{\Sigma }_{i=1}^n({\hat{y}}_i-y_i{)}^2}{n}}$$

and as you might intuit from the formula, RMSE is a measure of the standard deviation of residuals. RMSE can be computed using the metrics package, or just using a quick manual calculation:

sqrt(mean(model_1_df$.resid^2))

## [1] 4.97418

Our calculated RMSE of 4.974 indicates that actual life expectancy deviates from life expectancy predicted by the model by about 5 years, on average.

Another way to assess model performance is to figure out whether it was influenced by a small set of influential observations. Perhaps our model started out as a perfectly nice model, chugging along, predicting stuff with few mistakes. But then it came across some highly influential characters - data points that don’t hang with the pack, non-conformers, and they applied heavy influence and ultimately swayed our model off its course.

In truth, the study of outliers and influential observation is a whole thing and worthy of its own project. For now, let’s do two things. First, let’s look at Figure 2 and acknowledge that there’s not much visual evidence of outliers. Second, let’s do a quick check using Cook’s distance, .cooksd in our augment() dataframe.

Cook’s distance is a metric based on deleted residuals and is calculated for each data point in the set. It is a measure of the difference that would occur in our predicted values if we were to re-run the regression model without that point. Let’s look at the Cook’s distances in model_1:

model_1_df %>% 
  arrange(desc(.cooksd)) %>% 
  inner_join(lifexp_df, on = life_expectancy) %>% 
  head()

## # A tibble: 6 × 28
##   life_exp…¹ schoo…² .fitted .resid    .hat .sigma .cooksd .std.…³ country  year
##        <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>   <dbl>   <dbl> <chr>   <int>
## 1       63.6     5      53.3  10.3  0.0473    4.95  0.111     2.11 Eritrea  2012
## 2       63       5.1    53.5   9.48 0.0462    4.96  0.0912    1.94 Niger    2012
## 3       49.7     9.1    62.7 -13.0  0.0149    4.92  0.0518   -2.62 Sierra…  2012
## 4       77       9.9    64.5  12.5  0.0113    4.92  0.0360    2.51 Bangla…  2012
## 5       52.2    11      67.0 -14.8  0.00785   4.89  0.0351   -2.98 Lesotho  2012
## 6       52.7     9.7    64.1 -11.4  0.0121    4.94  0.0321   -2.29 Congo    2001
## # … with 18 more variables: status <chr>, adult_mortality <int>,
## #   infant_deaths <int>, alcohol <dbl>, percentage_expenditure <dbl>,
## #   hepatitis_b <int>, measles <int>, bmi <dbl>, under_five_deaths <int>,
## #   polio <int>, total_expenditure <dbl>, diphtheria <int>, hiv_aids <dbl>,
## #   gdp <dbl>, population <dbl>, thinness_1_19_years <dbl>,
## #   thinness_5_9_years <dbl>, income_composition_of_resources <dbl>, and
## #   abbreviated variable names ¹​life_expectancy, ²​schooling, ³​.std.resid

We can see that our top five most influential observations are Eritrea, Niger, Sierra Leone, Bangladesh, and Lesotho. To complement using a quick visualization:

model_1_df %>% 
  mutate(id = row_number()) %>% 
  ggplot(aes(x = id, y = .cooksd)) +
  geom_point() + 
  theme_bw()

We see three points (Eritrea, Niger, Sierra Leone) that stand out from the pack. These observations are not just outliers, but influential outliers, meaning they swayed our model towards themselves, compromising its predictive capability.

So what do we do with them? You can find lots of different “rules of thumb” online, dictating which Cook’s distance might prompt removal of a data point from your set, but what makes most sense to me in this situation is to accept the data as they are. We don’t have any super high leverage points here, and even if we did, we would need justification for removing them. As you can probably intuit, few natural processes are simple enough to be accurately explained by a linear model. Hence the discovery/invention of more sophisticated probabilistic models and machine learning.

Step 3: Validate Regression Assumptions

Now that we’ve done all this work, we need to figure out whether any of it was worthwhile. This is a bit of a nuisance with linear regression - you can’t really check the assumptions before you fit the model because you need the model to know if the assumptions were satisfied.

You’ll often see regression assumptions summarized as follows: 1. Linearity: If the data aren’t linear, don’t fit a linear model. We looked at the scatterplot in step 1 and found it was pretty linear. 2. Independence: Observations shouldn’t be clustered or have any influence on other observations. It’s hard to say this is the case in our situation, since it’s pretty clear that developing countries and developed countries cluster together and through political and economic means influence their neighbors’ policies and practices. 3. Normality of residuals: If you look at the residuals in Figure 3 and superimpose little sideways bell curves along the regression line, you should see that the residuals are normally distributed. Basically, the majority of the residuals are located close to the line and a few are found further away. This is best explained visually:

It looks like this is generally the case for our model, using the trusty eyeball method. 4. Homoscedasticity ¹: Homoscedasticity means that the variance of the residual terms is somewhat constant. If your data are heteroscedastic, linear regression modeling is probably not a good choice. Again, this is better visualized:

When we first looked at our data in scatterplot form above, we noted some definite heteroscedasticity. Let’s look at it again:

final_df %>% 
  ggplot(aes(x = schooling, y = life_expectancy)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Life Expectancy vs Schooling"
  )

Below 10 years of schooling, our residuals are much wider than on the other side of 10. There might be some intuitive explanation for this observation. For instance, it’s possible that there’s a distinction between countries where a majority of the population finishes what is considered high school in the US and in countries where that isn’t the case, years of schooling doesn’t matter as much as other factors - like access to clean water and nutritious foods.

All in all, checking regression assumptions convinced me that linear regression is not really the best choice here. But on the bright side, we learned some good fundamental principles about how statistical modeling works.