• Set Up
    • Question 1
    • Question 2
    • Question 3
    • Question 4
      • Part A
      • Part B
      • Part C
      • Part D
    • Question 5
    • Question 6
    • Question 7
      • Part A
      • Part B
      • Part C
      • Part D
      • Part E
    • Question 8
    • Question 9
    • Question 10
    • Question 11
      • Part A
      • Part B
      • Part C
      • Part D
      • Part E
      • Part F
      • Part G
    • Question 12
      • Part A
      • Part B
      • Part C
      • Part D
    • Question 13

Set Up

To minimize confusion, I suggest creating a new R Project (e.g. regression_practice) and storing any data in that folder on your computer.

Alternatively, I have made a project in R Studio Cloud that you can use (and not worry about trading room computer limitations), with the data already inside (you will still need to assign it to an object).

Question 1

Download and read in (read_csv) the data below.

  • speeding_tickets.csv This data comes from a paper by Makowsky and Strattman (2009) that we will examine later. Even though state law sets a formula for tickets based on how fast a person was driving, police officers in practice often deviate from that formula. This dataset includes information on all traffic stops. An amount for the fine is given only for observations in which the police officer decided to assess a fine. There are a number of variables in this dataset, but the one’s we’ll look at are:
Variable Description
Amount Amount of fine (in dollars) assessed for speeding
Age Age of speeding driver (in years)
MPHover Miles per hour over the speed limit

We want to explore who gets fines, and how much. We’ll come back to the other variables (which are categorical) in this dataset in later lessons.

library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
speed<-read_csv("../data/speeding_tickets.csv")
## Parsed with column specification:
## cols(
##   Black = col_double(),
##   Hispanic = col_double(),
##   Female = col_double(),
##   Amount = col_double(),
##   MPHover = col_double(),
##   Age = col_double(),
##   OutTown = col_double(),
##   OutState = col_double(),
##   StatePol = col_double()
## )

Question 2

How does the age of a driver affect the amount of the fine? Make a scatterplot of the Amount of the fine and the driver’s Age. Add a regression line with an additional layer of geom_smooth(method="lm").

ggplot(data = speed)+
  aes(x = Age,
      y = Amount)+
  geom_point()+
  geom_smooth(method = "lm")
## Warning: Removed 36683 rows containing non-finite values (stat_smooth).
## Warning: Removed 36683 rows containing missing values (geom_point).

Question 3

Find the correlation between Amount and Age.

# get summary of Amount
speed %>%
  select(Amount) %>%
  summary() # see there are 36,683 NA's!
##      Amount     
##  Min.   :  3    
##  1st Qu.: 75    
##  Median :115    
##  Mean   :122    
##  3rd Qu.:155    
##  Max.   :725    
##  NA's   :36683
# so we need to get correlation by dropping the NA's

# method 1: filter by non-NA observations for Amount
speed %>%
  select(Age, Amount) %>%
  filter(!is.na(Amount)) %>%
  cor()
##                Age      Amount
## Age     1.00000000 -0.06546571
## Amount -0.06546571  1.00000000
# method 2: tell cor() to use "pairwise.complete.obs"
speed %>%
  select(Age, Amount) %>%
  cor(., use = "pairwise.complete.obs")
##                Age      Amount
## Age     1.00000000 -0.06546571
## Amount -0.06546571  1.00000000
# note again I am using the tidyverse method of getting correlation

# you can also use summarize():
speed %>% summarize(my_cor = cor(Age, Amount, use = "pairwise.complete.obs"))
ABCDEFGHIJ0123456789
my_cor
<dbl>
-0.06546571
1 row
# or the Base R method:
cor(speed$Age, speed$Amount, use = "pairwise.complete.obs")
## [1] -0.06546571

Question 4

We want to predict the following model:

^Amounti=^β0+^β1Agei

Run a regression, and save it as an object. Now get a summary() of it.

reg1<-lm(Amount~Age, data = speed)
summary(reg1)
## 
## Call:
## lm(formula = Amount ~ Age, data = speed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -123.21  -46.58   -5.92   32.55  600.24 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 131.70665    0.88649  148.57   <2e-16 ***
## Age          -0.28927    0.02478  -11.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56.13 on 31672 degrees of freedom
##   (36683 observations deleted due to missingness)
## Multiple R-squared:  0.004286,   Adjusted R-squared:  0.004254 
## F-statistic: 136.3 on 1 and 31672 DF,  p-value: < 2.2e-16

Part A

Write out the estimated regression equation.

^Amounti=131.71−0.29Agei

# if you want to use equatiomatic
library(equatiomatic)
extract_eq(reg1, use_coefs = TRUE, fix_signs = TRUE, coef_digits = 2)
## $$
## \text{Amount} = 131.71 - 0.29(\text{Age}) + \epsilon
## $$

Part B

What is ^β0? What does it mean in the context of our question?

^β0=131.71. It means when Age is 0, the fine Amount is expected to be $131.71.

Part C

What is ^β1? What does it mean in the context of our question?

^β1=−0.29. It means for every additional year old someone is, their expected Amount of the fine decreases by $0.29.

Part D

What is the marginal effect of age on amount?

The marginal effect is ^β1, described in the previous part. Just checking to make sure you know it’s the marginal effect!

Question 5

Redo question 4 with the broom package. Try out tidy() and glance(). This is just to keep you versatile.

library(broom)
# tidy() to get coefficients
reg1_tidy<-tidy(reg1)

reg1_tidy
ABCDEFGHIJ0123456789
term
<chr>
estimate
<dbl>
std.error
<dbl>
statistic
<dbl>
p.value
<dbl>
(Intercept)131.70665080.88649485148.570130.000000e+00
Age-0.28927120.02477541-11.675741.967241e-31
2 rows
# glance() at original reg for measures of fit
glance(reg1)
ABCDEFGHIJ0123456789
r.squared
<dbl>
adj.r.squared
<dbl>
sigma
<dbl>
statistic
<dbl>
p.value
<dbl>
df
<int>
logLik
<dbl>
AIC
<dbl>
0.0042857590.00425432156.1254136.32281.967241e-312-172512.3345030.6
1 row | 1-8 of 11 columns

Question 6

How big would the difference in expected fine be for two drivers, one 18 years old and one 40 years old?

18-year-old driver:

^Amounti=^β0+^β1Agei=131.71−0.29(18)=126.49

40-year-old driver:

^Amounti=^β0+^β1Agei=131.71−0.29(40)=120.11

The difference is 126.49−120.11=6.38 only!

# we can use R as a calculator and use the regression coefficients

# extract beta0 and save it
beta_0<-reg1_tidy %>%
  filter(term=="(Intercept)") %>%
  select(estimate)

# extract beta1 and save it
beta_1<-reg1_tidy %>%
  filter(term=="Age") %>%
  select(estimate)

# predicted amount for 18 year old
amount_18yr<-beta_0+beta_1*18

# predicted amount for 40 year old
amount_40yr<-beta_0+beta_1*40

# difference
amount_18yr-amount_40yr
ABCDEFGHIJ0123456789
estimate
<dbl>
6.363966
1 row

Question 7

Now run the regression again, controlling for speed (MPHover).

reg2<-lm(Amount~Age+MPHover, data=speed)
summary(reg2)
## 
## Call:
## lm(formula = Amount ~ Age + MPHover, data = speed)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -308.763  -19.783    7.682   25.757  226.295 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.49386    0.95428   3.661 0.000251 ***
## Age          0.02496    0.01760   1.418 0.156059    
## MPHover      6.89175    0.03869 178.130  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.67 on 31671 degrees of freedom
##   (36683 observations deleted due to missingness)
## Multiple R-squared:  0.5026, Adjusted R-squared:  0.5026 
## F-statistic: 1.6e+04 on 2 and 31671 DF,  p-value: < 2.2e-16

Part A

Write the new regression equation.

^Amounti=3.49+0.02Agei+6.89MPHoveri

# if you want to use equatiomatic
extract_eq(reg2, use_coefs = TRUE, fix_signs = TRUE, coef_digits = 2)
## $$
## \text{Amount} = 3.49 + 0.02(\text{Age}) + 6.89(\text{MPHover}) + \epsilon
## $$

Part B

What is the marginal effect of Age on Amount? What happened to it?

^β1=0.02. For every additional year old someone is, holding their speed constant, their expected fine increases by $0.02. The magnitude of the effect shrank and it switched from negative to positive!

Part C

What is the marginal effect of MPHover on Amount?

^β2=6.89. For every additional MPH someone drives above the speed limit, holding their Age constant, their expected fine increases by $6.89

Part D

What is ^β0, and what does it mean?

^β0=3.49. This is the expected fine for a driver of Age =0 and MPHover =0.

Part E

What is the adjusted ˉR2? What does it mean?

It is 0.5026. We can explain 50% of the variation in Amount from variation in our model.

Question 8

Now suppose both the 18 year old and the 40 year old each went 10 MPH over the speed limit. How big would the difference in expected fine be for the two drivers?

18-year-old driver:

^Amounti=^β0+^β1Agei+^β2MPHoveri=3.49+0.02(18)+6.89(10)=72.75

40-year-old driver:

^Amounti=^β0+^β1Agei+^β2MPHoveri=3.49+0.02(40)+6.89(10)=73.19

The difference is 72.75−73.19=−0.44 only!

# we can use R as a calculator and use the regression coefficients

# first we need to tidy reg2
reg2_tidy<-tidy(reg2)

# extract beta0 and save it
multi_beta_0<-reg2_tidy %>%
  filter(term=="(Intercept)") %>%
  select(estimate)

# extract beta1 and save it
multi_beta_1<-reg2_tidy %>%
  filter(term=="Age") %>%
  select(estimate)

# extract beta2 and save it
multi_beta_2<-reg2_tidy %>%
  filter(term=="MPHover") %>%
  select(estimate)

# predicted amount for 18 year old going 10 MPH over
amount_18yr_10mph<-multi_beta_0+multi_beta_1*18+multi_beta_2*10

# predicted amount for 40 year old going 10 MPH over
amount_40yr_10mph<-multi_beta_0+multi_beta_1*40+multi_beta_2*10

# difference
amount_18yr_10mph-amount_40yr_10mph # close, we have some rounding error!
ABCDEFGHIJ0123456789
estimate
<dbl>
-0.5492232
1 row

Question 9

How about the difference in expected fine between two 18 year olds, one who went 10 MPH over, and one who went 30 MPH over?

18-year-old driver, 10 MPH over: we saw was $72.75.

18-year-old driver, 30 MPH over:

^Amounti=^β0+^β1Agei+^β2MPHoveri=3.49+0.02(18)+6.89(30)=210.55

The difference is now 210.55−72.75=−137.80!

# we can use R as a calculator and use the regression coefficients

# predicted amount for 40 year old going 30 MPH over
amount_18yr_30mph<-multi_beta_0+multi_beta_1*18+multi_beta_2*30

# difference
amount_18yr_30mph-amount_18yr_10mph
ABCDEFGHIJ0123456789
estimate
<dbl>
137.835
1 row

So clearly the speed plays a much bigger role than age does!

Question 10

Use the huxtable package’s huxreg() command to make a regression table of your two regressions: the one from question 4, and the one from question 7.

library(huxtable)
## 
## Attaching package: 'huxtable'
## The following object is masked from 'package:dplyr':
## 
##     add_rownames
## The following object is masked from 'package:purrr':
## 
##     every
## The following object is masked from 'package:ggplot2':
## 
##     theme_grey
huxreg(reg1,reg2) # the default
(1) (2)
(Intercept) 131.707 *** 3.494 ***
(0.886)    (0.954)   
Age -0.289 *** 0.025    
(0.025)    (0.018)   
MPHover          6.892 ***
         (0.039)   
N 31674         31674        
R2 0.004     0.503    
logLik -172512.295     -161520.081    
AIC 345030.591     323048.163    
*** p < 0.001; ** p < 0.01; * p < 0.05. 
# some customization
huxreg("Simple OLS"=reg1,
       "Mulitvariate OLS"=reg2,
       coefs = c("Constant" = "(Intercept)",
                 "Age" = "Age",
                 "MPH Over" = "MPHover"),
       statistics = c("N" = "nobs",
                      "R-Squared" = "r.squared",
                      "SER" = "sigma"),
       number_format = 2)
Simple OLS Mulitvariate OLS
Constant 131.71 *** 3.49 ***
(0.89)    (0.95)   
Age -0.29 *** 0.02    
(0.02)    (0.02)   
MPH Over         6.89 ***
        (0.04)   
N 31674        31674       
R-Squared 0.00     0.50    
SER 56.13     39.67    
*** p < 0.001; ** p < 0.01; * p < 0.05.

Look at that change in R2!

Question 11

Are our two independent variables multicollinear? Do younger people tend to drive faster?

Part A

Get the correlation between Age and MPHover.

# same alternative methods as question 3

# fortunately no NA's here!

speed %>%
  select(Age, MPHover) %>%
  cor()
##                Age    MPHover
## Age      1.0000000 -0.1035001
## MPHover -0.1035001  1.0000000

Part B

Make a scatterplot of MPHover on Age.

ggplot(data = speed)+
  aes(x = Age,
      y = MPHover)+
  geom_point()+
  geom_smooth(method = "lm")

Part C

Run an auxiliary regression of MPHover on Age.

reg_aux<-lm(MPHover~Age, data = speed)
summary(reg_aux)
## 
## Call:
## lm(formula = MPHover ~ Age, data = speed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.683  -4.020  -0.488   2.512  59.551 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.541094   0.054399  304.07   <2e-16 ***
## Age         -0.039007   0.001434  -27.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.056 on 68355 degrees of freedom
## Multiple R-squared:  0.01071,    Adjusted R-squared:  0.0107 
## F-statistic: 740.2 on 1 and 68355 DF,  p-value: < 2.2e-16

Part D

Interpret the coefficient on Age.

In our notation, this is ^δ1=−0.03, which says for every 1 year older someone is, they drive 0.03 MPH less over the speed limit.

Part E

Look at your regression table in question 10. What happened to the standard error on Age? Why (consider the formula for variation in ^β1)

The formula for var[^β1]=11−R21×SER2n×var[Age]

The standard error of the regression (SER) went from 56.13 to 39.67!1 Adding the second variable created a much better fit, which lowered variance on the betas (and hence, standard error, the square root of variance).

Variance might have slightly increased to more than it otherwise would be, due to multicollinearity between Age and MPHover, which we investigate next.

Part F

Calculate the Variance Inflation Factor (VIF) using the car package’s vif() command.2

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
vif(reg2)
##      Age  MPHover 
## 1.010149 1.010149

The VIF is 1.01 for ^β1 and ^β2. Variance increases only by 1% due to (weak) multicollinearity between Age and MPHover.

Part G

Calculate the VIF manually, using what you learned in this question.

VIF=11−R21, where R21 comes from the auxiliary regression of MPHover on Age (Part C). The R2 from that regression was 0.011, so:

VIF=11−R21=11−0.011=10.989≈1.011

# we can do this in R

# first we need to extract and save the R-squared from the aux reg
# use broom's glance()

aux_r_sq<-glance(reg_aux)$r.squared %>%
  round(.,3) # round to 3 decimals

vif<-1/(1-aux_r_sq)
vif
## [1] 1.011122

Again, we have some slight rounding error.

Question 12

Let’s now think about the omitted variable bias. Suppose the “true” model is the one we ran from Question 7.

Part A

Do you suppose that MPHover fits the two criteria for omitted variable bias?

  1. MPHover probably affects Amount
  2. cor(MPHover,Age)≠0

Possibly, though it is only weakly correlated with Age (−0.10). Speed clearly has a big role to play in affecting and predicting fines, but probably does not add very much bias to the marginal effect of Age on Amount.

Part B

Look at the regression we ran in Question 4. Consider this the “omitted” regression, where we left out MPHover. Does our estimate of the marginal effect of Age on Amount overstate or understate the true marginal effect?

Since there is negative correlation between MPHover and Age, ^β1 likely understates the effect of Age on Amount.

Part C

Use the “true” model (Question 7), the “omitted” regression (Question 4), and our “auxiliary” regression (Question 11) to identify each of the following parameters that describe our biased estimate of the marginal effect of Age on Amount:3 α1=β1+β2δ1

Using our notation from class, and the regressions from questions 4, 7, and 11, we have three regressions:

  • True Model: ^Amounti=3.49+0.02Agei+6.89MPHoveri
  • Omitted Reg: ^Amounti=131.71−0.29Agei
  • Auxiliary Reg: ^MPHoveri=16.54−0.04Agei

^α1=β1+β2δ1−0.29=0.02+6.89(−0.04)

Where:

  • β1=0.02: the true marginal effect of Age on Amount (holding MPHover constant)
  • β2=6.89: the true marginal effect of MPHover on Amount (holding Age constant)
  • δ1=6.89: the (auxiliary) effect of MPHover on Age

We have some slight rounding error, but this is the relationship.

Part D

From your answer in part C, how large is the omitted variable bias from leaving out MPHover?

^α1=β1+β2δ1⏟Bias−0.29=0.02+6.89(−0.04)⏟biasBias=−0.2756

Question 13

Make a coefficient plot of your coefficients from the regression in Question 7.

# first we have to make sure to tidy our reg with confidence intervals!
reg2_tidy<-tidy(reg2, conf.int=TRUE)

ggplot(data = reg2_tidy)+
  aes(x = estimate,
      y = term,
      color = term)+
    geom_point()+ # point for estimate
    # now make "error bars" using conf. int.
    geom_segment(aes(x = conf.low,
                     xend = conf.high,
                     y=term,
                     yend=term))+
  # add line at 0
  geom_vline(xintercept=0, size=1, color="black", linetype="dashed")+
  #scale_x_continuous(breaks=seq(-1.5,0.5,0.5),
  #                   limits=c(-1.5,0.5))+
  labs(x = expression(paste("Marginal effect, ", hat(beta[j]))),
       y = "Variable")+
  guides(color=F)+ # hide legend
  theme_classic(base_family = "Fira Sans Condensed",
           base_size=20)

  1. Note: you need to be sure to show sigma in your huxreg() or look at each of your regression results with broom’s glance() or Base R’s summary() output.↩︎

  2. Run it on your multivariate regression object from Question 7.↩︎

  3. See the notation from class 3.2.↩︎