Multiple Regression

As we broaden our scope beyond the relationship of only two variables (the dependent variable and a single regressor), some new issues can be seen, such as multicollinearity and omitted variable bias (OVB). Let us investigate this using R.

Let’s use a dataset that contains data on test performance, school characteristics and student demographic backgrounds for school districts in California, and define variables for the student-to-teacher ratio STR and test scores score.

# load the data set and define variables
if(!require(AER)) install.packages("AER")

## Loading required package: AER

## Loading required package: car

## Loading required package: carData

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

if(!require(MASS)) install.packages("MASS")

## Loading required package: MASS

if(!require(mvtnorm)) install.packages("mvtnorm")

## Loading required package: mvtnorm

library(AER)
library(MASS)
data(CASchools)   


CASchools$size <- CASchools$students/CASchools$teachers       
CASchools$score <- (CASchools$read + CASchools$math)/2

Performing a multiple regression in R is straightforward. One can simply add additional variables to the right hand side of the formula argument of the function lm() by using their names and the + operator. Let us estimate a model with and without adding a control for a the percent of english learner in a district english.

model1 <- lm(score ~ size, data = CASchools)
model2 <- lm(score ~ size + english, data = CASchools)

summary(model1)

## 
## Call:
## lm(formula = score ~ size, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.727 -14.251   0.483  12.822  48.540 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 698.9329     9.4675  73.825  < 2e-16 ***
## size         -2.2798     0.4798  -4.751 2.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared:  0.05124,    Adjusted R-squared:  0.04897 
## F-statistic: 22.58 on 1 and 418 DF,  p-value: 2.783e-06

summary(model2)

## 
## Call:
## lm(formula = score ~ size + english, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
## size         -1.10130    0.38028  -2.896  0.00398 ** 
## english      -0.64978    0.03934 -16.516  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

Examples of Perfect Multicollinearity

Let’s try adding a variable to our regression FracEL that is perfectly collinear with english.

# define the fraction of English learners        
CASchools$FracEL <- CASchools$english / 100

# estimate the model

return <- lm(score ~ size + english + FracEL, data = CASchools)

# obtain a summary of the model
    
summary(return)

## 
## Call:
## lm(formula = score ~ size + english + FracEL, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients: (1 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
## size         -1.10130    0.38028  -2.896  0.00398 ** 
## english      -0.64978    0.03934 -16.516  < 2e-16 ***
## FracEL             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

Graded Question 1

When happened in the output of our regression due to perfectly correlated regressors?

Answer: It gets assigned a NA value.

The Distribution of the OLS Estimators in Multiple Regression

# load packages
library(MASS)
library(mvtnorm)

# set sample size
n <- 50

# initialize vector of coefficients
coefs <- cbind("hat_beta_1" = numeric(10000), "hat_beta_2" = numeric(10000))

# set seed for reproducibility
set.seed(1)

# loop sampling and estimation
for (i in 1:10000) {
  
  X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
  u <- rnorm(n, sd = 5)
  Y <- 5 + 2.5 * X[, 1] + 3 * X[, 2] + u
  coefs[i,] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
  
}

# compute density estimate
kde <- kde2d(coefs[, 1], coefs[, 2])

# plot density estimate
persp(kde, 
      theta = 310, 
      phi = 30, 
      xlab = "beta_1", 
      ylab = "beta_2", 
      zlab = "Est. Density",
      main = "2D Kernel Density Estimate")

Omitted Variable Bias

Graded Question 2

Use the cor function to see if income is positively or negatively associated with class size and with test scores.

cor(CASchools$income, CASchools$size)

## [1] -0.2321937

cor(CASchools$income, CASchools$score)

## [1] 0.7124308

Now let’s regress test scores onto class size, with and without controlling for income.

#Estimate model with only size
SizeOnly_model <- lm(score ~ size, data = CASchools)

# compute robust summary
coeftest(SizeOnly_model, 
         vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 698.93295   10.36436 67.4362 < 2.2e-16 ***
## size         -2.27981    0.51949 -4.3886 1.447e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Estimate model with income and size
SizeAndIncome_model <- lm(score ~ size + income, data = CASchools)

# compute robust summary
coeftest(SizeAndIncome_model, 
         vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 638.72916    7.30124  87.482  < 2e-16 ***
## size         -0.64874    0.35334  -1.836  0.06707 .  
## income        1.83911    0.11473  16.029  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Graded Question 3

Is the bias on the coefficient for size from omitting income positive or negative? What about income relationship to class size and test scores might explain this?

Answer: The bias on the coefficient for ‘size’ is negative. For higher income relationships, they usually have a negative correlation between income and size. and a positive correlation relationship between income and score.

Nonlinear Regression

Let’s now consider the relationship between test scores and income.

# fit a simple linear model
linear_model<- lm(score ~ income, data = CASchools)
summary(linear_model)

## 
## Call:
## lm(formula = score ~ income, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.574  -8.803   0.603   9.032  32.530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 625.3836     1.5324  408.11   <2e-16 ***
## income        1.8785     0.0905   20.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.39 on 418 degrees of freedom
## Multiple R-squared:  0.5076, Adjusted R-squared:  0.5064 
## F-statistic: 430.8 on 1 and 418 DF,  p-value: < 2.2e-16

Here, income and test scores are positively related: school districts with above average income tend to achieve above average test scores. Does a linear regression function model the data adequately? Let us plot the data and add a linear regression line.

# plot the observations
plot(CASchools$income, CASchools$score,
     col = "steelblue",
     pch = 20,
     xlab = "District Income (thousands of dollars)", 
     ylab = "Test Score",
     cex.main = 0.9,
     main = "Test Score vs. District Income and a Linear OLS Regression Function")

# add the regression line to the plot
abline(linear_model, 
       col = "red", 
       lwd = 2)
legend("bottomright", legend="linear fit",lwd=2,col="red")

Fortunately, OLS does not only handle linear functions of the regressors. We can, for example, model test scores as a function of income and the square of income. The corresponding regression model is \[ TestScore_i = \beta_0 + \beta_1 \times income_i + \beta_2 \times income_i^2 + u_i \] Hence, the quadratic model is a special case of a multivariate regression model. When fitting the model with lm() we have to use the ^ operator in conjunction with the function I() to add the quadratic term as an additional regressor to the argument formula. This is because the regression formula we pass to formula is converted to an object of the class formula. For objects of this class, the operators +, -, * and ^ have a nonarithmetic interpretation. I() ensures that they are used as arithmetical operators.

Fit the quadratic Model

quadratic_model <- lm(score ~ income + I(income^2), data = CASchools)

# obtain the model summary
coeftest(quadratic_model, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) 607.3017435   2.9017544 209.2878 < 2.2e-16 ***
## income        3.8509939   0.2680942  14.3643 < 2.2e-16 ***
## I(income^2)  -0.0423084   0.0047803  -8.8505 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot Models

# draw a scatterplot of the observations for income and test score
plot(CASchools$income, CASchools$score,
     col  = "steelblue",
     pch = 20,
     xlab = "District Income (thousands of dollars)",
     ylab = "Test Score",
     main = "Estimated Linear and Quadratic Regression Functions")

# add a linear function to the plot
abline(linear_model, col = "green", lwd = 2)

# add quatratic function to the plot
order_id <- order(CASchools$income)

lines(x = CASchools$income[order_id], 
      y = fitted(quadratic_model)[order_id],
      col = "red", 
      lwd = 2) 
legend("bottomright",legend=c("Linear Line","Quadratic Line"),
       lwd=2,col=c("green","red"))

Logarithms

Another way to specify a nonlinear regression function is to use the natural logarithm of Y and/or X. Logarithms convert changes in variables into percentage changes. This is convenient as many relationships are naturally expressed in terms of percentages. There are three different cases in which logarithms might be used.

Transform X to its logarithm, but not Y.
Analogously we can transform Y to its logarithm but leave X at level.
Both Y and X are transformed to their logarithms.

The interpretation of the regression coefficients is different in each case.

Linear-Log Case:

\[Score_i = \beta_0 + \beta_1 \log(income_i) + u_i \]

Similar to polynomial regressions, we do not have to create a new variable before using lm(). We can simply adjust the formula argument of lm() to tell R that the log-transformation of a variable should be used.

LinearLog_model <- lm(score ~ log(income), data = CASchools)

# compute robust summary
coeftest(LinearLog_model, 
         vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 557.8323     3.8399 145.271 < 2.2e-16 ***
## log(income)  36.4197     1.3969  26.071 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Logarithm Plot

# draw a scatterplot
plot(score ~ income, 
     col = "steelblue",
     pch = 20,
     data = CASchools,
     ylab="Score",
     xlab="Income",
     main = "Linear-Log Regression Line")

# add the linear-log regression line
order_id  <- order(CASchools$income)

lines(CASchools$income[order_id],
      fitted(LinearLog_model)[order_id], 
      col = "red", 
      lwd = 2)
legend("bottomright",legend = "Linear-log line",lwd = 2,col ="red")

Graded Questions 4-6

Based on the slope coefficient, a 1% increase in income is associated with how much of an increase in test scores?

Answer: A one percent increase would result in an increase in test scores that are equal to approximately .36 points

Rerun the regression using log(score) and income instead. How do we interpret the the slope coefficient now?

return <- lm(log(score) ~ income, data = CASchools)
coeftest(return, vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 6.43936234 0.00289382 2225.210 < 2.2e-16 ***
## income      0.00284407 0.00017509   16.244 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Answer: A increase in income would result in a 100xB1 percet increase in test scores.

Repeat Question 5 but now with log(score) and log(income).

return <- lm(log(score) ~ log(income), data = CASchools)
coeftest(return, vcov = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 6.3363494  0.0059246 1069.501 < 2.2e-16 ***
## log(income) 0.0554190  0.0021446   25.841 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Answer: A 1% increase in district income would result in a B1% increase in test scores.

Interactions between Independent Variables

There are research questions where it is interesting to learn how the effect on Y of a change in an independent variable depends on the value of another independent variable. For example, we may ask if districts with many English learners benefit differently from a decrease in class sizes to those with few English learning students.

To assess this using a multiple regression model, we include an interaction term. We consider three cases:

Interactions between two binary variables.
Interactions between a binary and a continuous variable.
Interactions between two continuous variables.

Example with Interactions Between a Binary and Continuous Variable

Let’s make a binary variable HiEL to indicate if percent of English learners is at least 10.

# append HiEL to CASchools
CASchools$HiEL <- as.numeric(CASchools$english >= 10)

We proceed by estimating the model: \[TestScores_i = \beta_0 + \beta_1 HiEL + \beta_2 Size + \beta_3 HiEL \times Size + u_i\] There are several ways to add the interaction term to the formula argument when using lm()

# estimate the model with an interaction term
inter_model <- lm(score ~ HiEL * size, data = CASchools)
#bi_model <- lm(score ~ HiSTR +  HiEL + HiSTR:HiEL, data = CASchools)


# print a robust summary of the coefficients
coeftest(inter_model, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 682.24584   11.86781 57.4871   <2e-16 ***
## HiEL          5.63914   19.51456  0.2890   0.7727    
## size         -0.96846    0.58910 -1.6440   0.1009    
## HiEL:size    -1.27661    0.96692 -1.3203   0.1875    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The estimated regression model is ˆ \[TestScores_i = 682.2 + 5.6 × HiEL − 1.0 × Score - 1.3 × (HiEL × Score)\]

plot(CASchools$size, CASchools$score, xlab = "Class Size", ylab = "Test Score")
abline(a = inter_model$coefficients[1], b = inter_model$coefficients[3], col = "red")
abline(a = inter_model$coefficients[1]+inter_model$coefficients[2], b = inter_model$coefficients[3]+inter_model$coefficients[4], col = "blue")

Graded Question 7

Is the interaction between the student-to-teacher ratio and percent English learners significant at the 0.05 level?

Answer:No, because the p-value is greater than .05, we fail to reject the null hypothesis.

Graded Question 8

Interpret the coefficient on size.

Answer:The base coefficient shows that for districts with fewer than approximately 10% of English learners, an increase of 1 student per teacher means an approximate decrease by 1 point in test scores.

R Day 3 - Multiple Regression

Dev Amin