Simple Linear Regression with GDP Data

Often, when conducting quantitative research, you need to run regression to examine a theory or to explore the relationship between variables. In this section, we use the GDP per capita and life expectancy data from Gapminder again to study how to estimate a simple linear model in R. We will also do some test in the theory of economic convergence and implement what data visualization and transformation we have learned in the previous two sections.

Economic Convergence

Before studying linear regression in R, let’s first try to examine a phenomenon in economic growth - convergence growth. The idea of convergence growth is that poor countries should be able to grow faster than rich countries. According to the Solow growth model, this is because the capital per capita or capital per efficiency unit of labor in poor countries grows faster in poor countries (far from steady state). Their low levels of capital suggests that there should be high returns to investment. Further, since they can copy frontier technologies, they should be able to improve their technology levels faster.

Prerequisites

Again, we will use GDP per capita data and life expectancy from Gapminder and the Penn World Table. First, let’s load required data and packages into current R session.

library(tidyverse)
library(gapminder)
gapminder<-gapminder

Because we want to to focus on long term economic growth instead of business cycles, we selected the very first year and the very last year available in our data: 1952 and 2007.

gapminder5207<-filter(gapminder,year==1952|year==2007)

Create variable g to represents the average growth rate of GDP per capita between 1952 and 2007.

gapminder5207<-gapminder5207 %>%
  group_by(country) %>%
  mutate(g = log(lead(gdpPercap)/gdpPercap)/(2007-1952)) %>%
  filter(year == 1952)

Rename gdpPercap as y to make our coding easier.

gapminder5207<-rename(gapminder5207, y=gdpPercap)

Convergence Fails?

Let’s first check whether poor countries, in general, grew faster in the period 1952-2007 time periods. To test this, draw a graph of initial GDP on the x-axis and subsequent growth on the y-axis. Use natural log to rescale the x-axis when draw the graph.

ggplot(gapminder5207)+
  geom_text(mapping = aes(x=log(y), y=g, label=country), size = 3)

Apparently, there is very little connection between initial GDP per capita and subsequent growth. If there was convergence, this graph should slope downwards: rich countries should grow slower. Thus, convergence fails! We actually see some rich countries grow really fast.

Does that me our the prediction from convergence theory and Solow model fails? The answer is NO! An obvious issue is omitted variable bias. Countries that were poor in 1960 had other problems that might have made them less likely to grow. To test if this is a relevant concern, plot 1952 GDP per capita on the x-axis, and the life expectancy on 1952 on the y-axis as we did in the Data Visualization section. Remember to use natural log to rescale the x-axis.

ggplot(gapminder5207)+
  geom_text(mapping = aes(x=log(y), y=lifeExp, label=country), size = 3)+
  geom_smooth(mapping = aes(x=log(y), y=lifeExp), method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

As you see when you do the plot, there is an extremely strong relationship between income levels and the people’s life expectancy in 1952. Since there are fewer people who are able to work healthily, the low life expectancy is a likely candidate for subsequent low growth.

A natural way to check how much this could matter is to plot the life expectancy against subsequent growth. Do this by plotting life expectancy on the x-axis and g on the y-axis. Add both a linear plot and a smooth plot using geom_smooth.

ggplot(gapminder5207)+
  geom_text(mapping = aes(x=lifeExp, y=g, label=country), size = 3)+
  geom_smooth(mapping = aes(x=lifeExp, y=g), method = "lm", se = FALSE)+
  geom_smooth(mapping = aes(x=lifeExp, y=g), linetype = "dashed", se=FALSE)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

If you did this correctly, you should find a weak positive relationship between life expectancy and subsequent growth, but it does not jump out of the data. Indeed, if you allow for a smooth fit, you will find a non-monotonic relationship, with growth turning down as life expectancy years become sufficiently high.

What is going on here? Is longer life expectancy not good for growth? Or is it first good, and then bad? Not so fast. Remember that when we use life expectancy, we now have a new omitted variable: the initial GDP per capita. Since the countries with low life expectancy also have low incomes, their growth will be higher thanks to convergence dynamics. Even if there was a positive relationship between life expectancy and growth, this might be hid by the convergence dynamics.

So how do we get out of this bind? The key insight is that we should look whether countries have high life expectancy conditional on their income level (like the second graph). So in the next part of this handout, we fit a linear model of life expectancy on GDP per capita.

Linear Model

A linear model has the general form y = a_1 + a_2 * x_1 + a_3 * x_2 + ... + a_n * x_(n - 1), where y is the dependent variable and there are n independent variables (one constant and n-1 other regressors). We want regress life expectancy on income per capita, so the simple model here is equivalent to a general linear model where n is 2 and x_1 is x. Our model to be estimated is \[\text{Life Expectancy}_i=a_0+a_1 Income_i\]

R has a tool specifically designed for fitting linear models called lm(). lm() has a special way to specify the model family: formulas. Formulas look like y ~ x, which lm() will translate to a function like y = a_1 + a_2 * x, thus regressing y on a constant (notice that R does this automatically) and the regressor x.

We can fit the model and look at the output:

model <- lm(lifeExp ~ log(y), data=gapminder5207)
summary(model)
## 
## Call:
## lm(formula = lifeExp ~ log(y), data = gapminder5207)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.9571  -5.7319   0.7517   6.5770  13.7361 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.8457     5.0668  -3.522 0.000578 ***
## log(y)        8.8298     0.6626  13.326  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.146 on 140 degrees of freedom
## Multiple R-squared:  0.5592, Adjusted R-squared:  0.556 
## F-statistic: 177.6 on 1 and 140 DF,  p-value: < 2.2e-16

We can also create a table for our estimates using a package called stargazer.

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
stargazer(model, type = "html", dep.var.labels = "Life Expectancy", 
          covariate.labels = c("GDP per capita"), out = "model2.htm")
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Life Expectancy</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">GDP per capita</td><td>8.830<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.663)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-17.846<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(5.067)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>142</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.559</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.556</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>8.146 (df = 140)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>177.587<sup>***</sup> (df = 1; 140)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Predicted Values

Add fitted values (called predictions in R) to gapminder5207. You should use modelr::add_predictions() which takes a data frame and a model as arguments. It adds the predictions from the model to a new column in the data frame.

library(modelr)
gapminder5207 <- add_predictions(gapminder5207,model)

We can further plot the predicted line from the regression of life expectancy on log of income per capita, overlaying it on a scatter plot of the data.

ggplot(gapminder5207)+
  geom_point(mapping = aes(x=log(y), y=lifeExp))+
  geom_smooth(mapping = aes(x=log(y), y=pred))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Residuals

The flip-side of predictions are residuals. The predictions tells you the pattern that the model has captured, and the residuals tell you what the model has missed. The residuals are just the distances between the observed and predicted values in the fitted plot that you drew. You can add residuals to the data with add_residuals(), which works much like add_predictions().

gapminder5207 <- add_residuals(gapminder5207,model)

There are a few different ways to understand what the residuals tell us about the model. One way is to simply draw a frequency polygon to help us understand the spread of the residuals:

ggplot(gapminder5207) + 
  geom_freqpoly(mapping=aes(resid))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

This helps you calibrate the quality of the model: how far away are the predictions from the observed values? Note that the average of the residual will always be 0. You’ll often want to recreate plots using the predicted values and residuals.