Lesson 11 of 17 · Week 05 · Tue Oct 27

Regression in Practice

A question worth a line

Do people in richer countries live longer? Look at one year of gm:

library(ggplot2)
gm07 <- gm[gm$year == 2007, ]
ggplot(gm07, aes(x = gdpPercap, y = lifeExp)) +
  geom_point() +
  scale_x_log10()

On a log scale for GDP, the points fall roughly along a straight line. A linear regression finds the best-fitting line and gives you numbers to describe it: how steep it is, and how closely the points follow it. We'll use log10(gdpPercap) as the predictor, because each step of 1 on that scale is a ten-fold increase in income.

This lesson was written for this course; it has no Carpentries counterpart.

Fitting a model with lm()

fit <- lm(lifeExp ~ log10(gdpPercap), data = gm07)
fit

The formula y ~ x reads “y explained by x”: the response goes on the left of the ~, the predictor on the right. data = says where to find those columns. The result, fit, is a model object: a list holding the estimates, the fitted values, the residuals and more (names(fit)).

Reading the coefficients

coef(fit)
#      (Intercept) log10(gdpPercap)
#            a                b
b <- coef(fit)[2]        # the slope
a <- coef(fit)[1]        # the intercept
a + b * log10(5000)      # predicted life expectancy at $5,000 per person
predict(fit, newdata = data.frame(gdpPercap = 5000))   # the same, with predict()

Slope: for each ten-fold increase in GDP per capita, predicted life expectancy is b years higher. Intercept: the prediction when log10(gdpPercap) is 0, i.e. at $1 per person. That's far outside the data, so on its own it doesn't mean much; it's there to position the line.

Reading summary()

summary(fit)
summary(fit)$r.squared           # R-squared
summary(fit)$coefficients        # estimates, std. errors, t values, p-values

From the output, read three things:

• Estimate: the same numbers as coef().
• Pr(>|t|) for the slope: the p-value. A small value (say below 0.05) says a slope this far from zero would be very unlikely if there were no linear relationship at all. It says nothing about how big or how important the relationship is.
• Multiple R-squared: the proportion of the variation in life expectancy that the line accounts for, from 0 (none) to 1 (all points exactly on the line).

We won't go further into inference (standard errors, confidence intervals, assumptions) in this course. Just read these three numbers carefully.

Plotting the fit

ggplot(gm07, aes(x = gdpPercap, y = lifeExp)) +
  geom_point() +
  scale_x_log10() +
  geom_smooth(method = "lm", se = FALSE)

Because the x-axis is already on a log scale, geom_smooth(method = "lm") fits the same line as lm(lifeExp ~ log10(gdpPercap)).

Residuals: what the line misses

gm07$fitted <- fitted(fit)
gm07$resid  <- residuals(fit)          # observed minus fitted
gm07[order(gm07$resid), c("country", "lifeExp", "fitted", "resid")]
ggplot(gm07, aes(x = fitted, y = resid)) + geom_point() + geom_hline(yintercept = 0)

A residual is how far a point sits above (positive) or below (negative) the line. Large residuals are often the most interesting part of an analysis. Why would a country live much shorter lives than its income predicts? In the residual plot, you hope to see a shapeless cloud around zero. A curve or a funnel means the straight line is missing something.

A line is not a cause. Regression describes an association in these data. Income, health systems, education and history are all tangled together. The model can't tell you that raising GDP will raise life expectancy.

Exercises

Each exercise gives you up to 10 checked attempts. Use Run as many times as you like to experiment first, then click Check answer when you're ready. Your exercises use numbers/subsets unique to you, so answers won't match a classmate's exactly.