Lecture

Video on YouTube

Script

LinearRegression.R

Linear regression

Testing a continuous response variable against a continuous predictor variable is called linear regression.

To present linear model fits graphically, start with a scatterplot of the data…

sp_gg <- 
  ggplot(mtcars, aes(x=hp, y=mpg)) + theme_bw(14) +
    geom_point()

… and add a trendline with geom_smooth:

sp_gg + geom_smooth(method='lm', se=FALSE)

geom_smooth is powerful. It can fit a line for almost any statistical model one can make in \({\bf\textsf{R}}\). By default, the trendline is a curvy loess function, so we need to specify method='lm'; this actually sends geom_smooth off to calculate a linear regression using lm and then it plots these results, similar to how stat_functionchugs a function defined by fun= and parameterized with args=. Indeed, geom_smooth will pass any optional arguments to the function specified in method= via a similar args=list() construction.

By default, geom_smooth adds a 95% confidence interval around the trendline. Here’s a totally default call:

sp_gg + geom_smooth( )

I don’t know why the default is to add the shaded 95% confidence interval, but you should pretty much always turn it off with se = F. ¯\_(ツ)_/¯.

We use lm() to test for a linear relationship between the (log) response variable, lmpg, and the predictor variable, hp:

hp_lm <- lm(lmpg ~ hp, mtcars)

There are several options view the results of the linear model.

coef() returns just the \(\beta\) coefficients:

coef(hp_lm)

##  (Intercept)           hp 
##  3.460466874 -0.003428734

summary() returns the full model results, including both the results of the linear model, goodness-of-fit statistics (\(R^2\)), and an \(F\) test. One can focus on the linear model results, stored in the coefficients matrix:

summary(hp_lm)$coefficients

##                 Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  3.460466874 0.0785838302 44.035355 8.047220e-29
## hp          -0.003428734 0.0004866904 -7.045001 7.853353e-08

If one is interested in the \(F\) test, one can call anova() to get the full ANOVA table:

anova(hp_lm)

## Analysis of Variance Table
## 
## Response: lmpg
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## hp         1 1.7132 1.71320  49.632 7.853e-08 ***
## Residuals 30 1.0355 0.03452                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note that although the \(P\) values are the same, these tables perform different calculations. The \(t\) statistic is the ratio of the \(\beta\) coefficient to the error associated with the estimate, or

\(t = -6.74 = \frac{-0.0682282}{0.010119}\)

It maintains the sign of the \(\beta\) estimate because the \(t\) distribution is two-sided.

In the ANOVA, however, the \(F\) statistic is calculated as the dividend of the mean sum of squares of the term over the residual error:

\(F = \frac{Mean~Sum~Squares_{Group}}{Mean~Sum~Squares_{Residuals}}\)

\(45.46 = \frac{678.37}{14.92}\)

Because these values are derived from the sums of squared error, the sign of negative effects are lost, thus the \(F\) distribution is one-sided. Therefore one would need to see the data to determine the direction of a significant effect.

To circle back to the discussion of weak alternative hypotheses \(H_a\) vs. strong \(H_A\), \(t\) statistics are better at informing the latter because the two-sided \(t\) distribution can indicate whether the slope is positive or negative. \(F\) statistics can only inform the weaker \(H_a\), reporting whether there is any evidence against the null hypothesis. But the predicted direction of the effect should still be stated a priori; don’t just state a lame \(H_a\) and rely upon the \(t\) statistic to tell you which direction it goes.

Fitting multiple regression models

Multiple regression is testing a single response variable against two or more predictor variables. The additional variables can be continuous or categorical.
We’ll start with two continuous variables.

Our trusty response variable is, of course, fuel economy, measured as miles per gallon and stored as mpg. But there are several variables in the dataset that could potentially explain variability in fuel economy: More powerful cars might use more fuel (larger engines with more cylinders, gear ratios suited more for power than economy); larger, heavier cars might take more fuel to drag around.

We can use multiple regression to compare the relative effect of several potentially explanatory variables:

• hp: horsepower, a US Customary measure of power produced by the engine
• wt: how heavy is the car that the engine has to pull around?
• drat: How many times the axle turns for every turn of the driveshaft

Recall the statistical formula for a line:

\(y_i = \beta_0 + \beta_1 x_i + \varepsilon\)

where \(\beta_0\) is the Y axis intercept, \(\beta_1\) is the slope, and \(\varepsilon\) is the residual error.

The formula for multiple regression modifies it only slightly:

\(y = \beta_0 + \beta_1 x_1 + ... \beta_n x_n + \varepsilon\)

where \(\beta_0\) is the Y axis intercept, \(\beta_1\) is the slope of the first term, \(x_1\), with additional variables added – each with their own \(\beta\) coefficient, up to \(x_n\).

The formula for our example would look something like

\(mpg = \beta_0 + \beta_1 hp + \beta_2 wt + \beta_3 drat + \varepsilon\)

In R’s y ~ x formula notation, we add additional predictor variables to the RHS separated by \(+\):

mr_lm <- lm(mpg ~ hp + wt + drat, mtcars)

Assessing model results

summary(mr_lm)

## 
## Call:
## lm(formula = mpg ~ hp + wt + drat, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3598 -1.8374 -0.5099  0.9681  5.7078 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29.394934   6.156303   4.775 5.13e-05 ***
## hp          -0.032230   0.008925  -3.611 0.001178 ** 
## wt          -3.227954   0.796398  -4.053 0.000364 ***
## drat         1.615049   1.226983   1.316 0.198755    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.561 on 28 degrees of freedom
## Multiple R-squared:  0.8369, Adjusted R-squared:  0.8194 
## F-statistic: 47.88 on 3 and 28 DF,  p-value: 3.768e-11

summary(mr_lm)$coefficients %>%
  pander::pander() 

coef(mr_lm)

## (Intercept)          hp          wt        drat 
##  29.3949345  -0.0322304  -3.2279541   1.6150485

sc1_lm <- 
  mtcars %>%
    select(mpg, hp, wt, drat) %>%
      mutate_all(~scale(.)) %>%
        lm(mpg ~ hp + wt + drat, .)
coef(sc1_lm) %>% 
  pander::pander() 

sc2_lm <- lm(scale(mpg) ~ scale(hp) + scale(wt) + scale(drat), mtcars)
coef(sc2_lm) %>% 
  pander::pander() 

confint(mr_lm) %>%
  pander::pander() 

Visualization

Visualizing three continuous predictor variables is difficult (impossible?) in a scatterplot. To see the pairwise relationships in the raw data, one can use a scatterplot matrix:

pacman::p_load(GGally)
ggpairs(mtcars, columns = c("mpg","hp", "wt", "drat"))

Once the multiple regression has been fit, the coefficients and the confidence intervals can be combined into a table for ggplot.

It is a bit tricky, though, because coef() prints in wide format and confint() prints in long format. Here’s a simple function called coefR() that combines them into a tidy object:

coefR <- function(x) {
              require(tidyverse)
              confint(x) %>%
              as.data.frame() %>%
              rownames_to_column("name") %>% 
              full_join(enframe(coef(x))) %>%
              setNames(c("term", "lwr", "upr", "est")) %>%
              select(term, lwr, est, upr)
            }

Now call the function on the lm object:

terms <- coefR(mr_lm)
terms %>% pander::pander() 

We can easily plot the data in this table to visualize the regression results. Unless one has a specific way to interpret it, the intercept term, \(\beta_0\), isn’t much use in a multiple regression model, so we’ll simply remove the top row using slice:²

terms %>%
  slice(-1) %>%
ggplot(aes(x = term)) + theme_bw(14) +
  geom_hline(yintercept = 0, lty = 2, color="grey70") + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.1, size = 1, col="blue") + 
  geom_point(aes(y = est), shape = 21, 
             col = "blue", fill = "lightblue", 
             size = 4, stroke = 1.25) +
  coord_flip() 

One can see the effect of scale mis-matches; the interval around hp essentially disappears. The size of the point is such that it looks like it overlaps zero.

Let’s re-plot with the model fit to scaled data:

coefR(sc1_lm) %>% 
  slice(-1) %>%
ggplot(aes(x = term)) + theme_bw(14) +
  geom_hline(yintercept = 0, lty = 2, color="grey70") + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), 
                width = 0.1, size = 1, col="blue") + 
  geom_point(aes(y = est), shape = 21, 
             col = "blue", fill = "lightblue", 
             size = 4, stroke = 1.25) +
  coord_flip() +
  labs(y = "Regression coefficient", 
       x = "Model term")

This gives us a much different picture of the relative effect of weight and horsepower. When plotting CIs, it is more important to consider scaling data prior to fitting the model, because the visual impact of the plot will be relative values, not just rank.

Continuous + categorical variables

While multiple continuous explantory variables are often tested together to determine their relative effects, additional categorical variables might be included to determine how the effects of predictor variables are conditioned by other variables, i.e., how does the response variable \(y\) change with \(x_1\) if \(x_2\) is held at a constant value? (One can also test how \(y\) changes as both \(x_1\) and \(x_2\) change by adding an interaction term; see below.)

For example, we know that fuel economy is affected by both transmission type and vehicle weight:

We might wonder, does the weight trend hold within each transmission type? Which variable is more important to fuel economy–a lighter car, or one with a manual transmission?

We can visualize the first question by adding the categorical variable as a conditional aesthetic in the scatterplot:

ggplot(mtcars, 
           aes(x = wt, 
               y = mpg)) + theme_bw(16) +
    geom_smooth(aes(color = transmission), 
                method="lm", se = F) +
    geom_point(aes(fill=transmission), 
               pch = 21, size =2, stroke = 1.1) +
    labs(y = "Fuel economy (mpg)", 
         x = "Vehicle weight (x 1000 lbs)")

The plot indicates that for both transmission types, fuel economy declines as cars get heavier. The decline seems more drastic, however, for cars with manual transmissions.

Fit a linear model:

cat_lm <- lm(mpg ~ wt + transmission, mtcars)
summary(cat_lm)$coefficients %>%
  pander::pander() 

These results indicate that in a model that includes vehicle weight, transmission type does not have a significant effect on fuel economy.

Another way to determine whether a term is worth including in statistical models is to compare models with and without the term using the anova() function. This method compares the Residual Sums of Squares (RSS) between the models; the lower the RSS, the better the model fits.

If the additional term is meaningful in explaining variation in the response variable, the models will be statistically-signficantly different. If the additional term does not add meaningful value, the models will not be different:

wt_lm <- lm(mpg ~ wt, mtcars) 
anova(wt_lm, cat_lm) 

## Analysis of Variance Table
## 
## Model 1: mpg ~ wt
## Model 2: mpg ~ wt + transmission
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1     30 278.32                          
## 2     29 278.32  1 0.0022403 2e-04 0.9879

These results confirm the above: adding the transmission term does not explain signficantly more variation in fuel economy than vehicle weight alone.

NOTE

Comparing linear models with anova() should be used sparingly.

Do not compare several models at once, nor models with substantially different numbers of terms. Adding terms changes the relationships between the models and their relative probabilities of including a particularly highly-explantory variable by chance, simply because there are many variables. Proper model selection methods use an information-theoretic criterion to rank models, and those methods typically have a means for taking the number of terms into account. We don’t cover it here, but it is covered elsewhere.

So why do cars with manual transmissions appear to have better fuel economy?

After considering all of these graphs, it appears that the reason is not that having the driver change gears saves energy, but that vehicles with manual transmissions simply tend to be lighter. This illustrates an issue with testing hypotheses from observational data–the mtcars dataset does not provide a complete factorial design in which all vehicle weights are represented by both manual and automatic transmissions. Really only then could we make a conclusive statement about which is more important about cars in general (beyond this specific dataset).

lm(y ~ x1 * x2, data)

lm(y ~ x1 + x2 + x1:x2, data)

int_lm <- lm(mpg ~ wt + transmission + wt:transmission, mtcars)
anova(wt_lm, cat_lm, int_lm)

## Analysis of Variance Table
## 
## Model 1: mpg ~ wt
## Model 2: mpg ~ wt + transmission
## Model 3: mpg ~ wt + transmission + wt:transmission
##   Res.Df    RSS Df Sum of Sq       F   Pr(>F)   
## 1     30 278.32                                 
## 2     29 278.32  1     0.002  0.0003 0.985556   
## 3     28 188.01  1    90.312 13.4502 0.001017 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1