Part A

Interpret what ^β0 means in this context.

^β0 is 696.70. This is the vertical intercept of the regression line. It means that a person that is 0 years old earns a $696.70 per week on average. This is often nonsensical, so we don’t often care about the economic meaning of the intercept.

Part B

Interpret what ^β1 means in this context.

^β1 is 9.60 This is the slope of the regression line. It means that for every year older a person is, they can expect their wages to increase by $9.60, on average. This is the marginal effect of Age on AWE (and the causal effect if this model were exogenous).

Part C

The R2=0.023 for this regression. What are the units of the R2, and what does this mean?

R2 has no units, it is the proportion of variation in AWE that is explained by our model, between 0 and 1. This model explains only 2.3% of the variation in AWE, meaning this model is poor, and the line does not fit the data points well.

Part D

The SER,^σu=624.1 for this regression. What are the units of the SER in this context, and what does it mean? Is the SER large in the context of this regression?

SER is measured in the same units as the dependent variable, AWE, so it is measured in dollars. It is the average error or residual for an individual, the difference (in dollars) between OLS’ predicted ^AWE for that person, and their true AWE in the data. This SER is quite big, $624 in average weekly earnings.

Part E

Suppose Maria is 20 years old. What is her predicted ^AWE?

^AWEMaria=696.70+9.60(20)=888.70

Part F

Suppose the data shows her actual AWE is $430. What is her residual? Is this a relatively good or a bad prediction?²

ˆuMaria=YMaria−ˆYMaria=430−888.70=−458.70

This actually a relatively good prediction, as it is much lower than the average prediction error (SER), which was 624.1.

Part G

What does the error term, ui represent in this case? What might individuals have different values of ui?

The error term represents all factors other than age that affects an individual’s average weekly earnings. This could include things like experience, ability, job type, education level, conscienciousness etc.

Part H

Do you think that Age is exogenous? Why or why not? Would we expect ^β1 to be too large or too small?

It’s very unlikely that Age is exogenous. Knowing someone’s age likely gives us information about u: we can guess about their experience or level of education (they are likely higher for older people), and most of these positively affect wages. Thus, we have probably overstimated the effect of age on earnings (i.e. ^β1), and the true β1 is likely smaller.

Part A

What is the OLS estimate of ^β1?

The formula for ^β1=n∑i=1(Xi−ˉX)(Yi−ˉY)n∑i=1(Xi−ˉX)2=cov(X,Y)var(X)=138006900=2

Part B

What is the OLS estimate of ^β0?

The formula for ^β0=ˉY−^β1ˉX=63−30(2)=3

Part C

Suppose the OLS estimate of ^β1 has a standard error of 0.072. Could we probably reject a null hypothesis of H0:β1=0 at the 95% level?

Yes, we could reject the null hypothesis as the estimate of ^β1=2 is more than 2 times its standard error of 0.072. The test-statistic would actually be

t=^β1−β1,0se(^β1)t=2−00.072t≈27.78

This is well beyond the critical value needed to reject H0, and the p-value would be basically 0.

Part D

Calculate the R2 for this model. How much variation in Y is explained by our model?

We know TSS (4 bullet point) and SSE (last bullet point).

R2=1−SSETSS=1−165629000=1−0.057=0.943

This model explains 94.3% of the variation in Yi.

Part E

How large is the average residual?

We need to find the standard error of the regression (SER), but luckily we know the SSE (last bullet point)

SER=√SSEn−2=√165648−2=√36=6

This tells us the average residual is 36 (units of Y).

Part A

Clean up the data a bit by making a new variable to measure home attendance in millions. This will make it easier to interpret your regression later on.

# make home attendance variable in millions
mlb<-mlb %>%
  mutate(home_attend_mil=home_attend/1000000)

Part B

Get the correlation between Runs Scored and Home Attendance.

# summarize and get correlation
mlb %>%
  summarize(Correlation = cor(runs_scored, home_attend_mil))

Part C

Plot a scatterplot of Runs Scored (y) on Home Attendance (x). Add a regression line.

# create scatterplot with regression line 
scatter<-ggplot(data = mlb)+
  aes(x = home_attend_mil,
      y = runs_scored)+
  geom_point()+
  geom_smooth(method="lm")
scatter

Part D

Run a regression of Runs Scored on Home Attendance. What are β0 and ^β1? Interpret them in the context of our question.

# run regression, save as reg
reg<-lm(runs_scored ~ home_attend_mil, data = mlb)

# get summary of reg
summary(reg)

## 
## Call:
## lm(formula = runs_scored ~ home_attend_mil, data = mlb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -295.566  -52.754    1.414   63.769  271.377 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      572.618      8.081   70.86   <2e-16 ***
## home_attend_mil   68.798      4.183   16.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 91.47 on 836 degrees of freedom
## Multiple R-squared:  0.2445, Adjusted R-squared:  0.2436 
## F-statistic: 270.5 on 1 and 836 DF,  p-value: < 2.2e-16

# Here I'm going to save beta 0 hat and beta 1 hat
# as objects to call up in the text of the markdown document
# We'll need broom and tidy() first
library(broom)
reg_tidy<-tidy(reg)

reg_tidy

# first beta 0 hat 

beta_0_hat<-reg_tidy %>%
  filter(term=="(Intercept)") %>% # look at intercept row
  pull(estimate) %>% # extract beta 0 hat
  round(., 2) # round to 2 decimal places

beta_0_hat # check 

## [1] 572.62

# now beta 1 hat 

beta_1_hat<-reg_tidy %>%
  filter(term=="home_attend_mil") %>% # look at X-variable row
  pull(estimate) %>% # extract beta 1 hat
  round(., 2) # round to 2 decimal places

beta_1_hat

## [1] 68.8

^β0= 572.62 ^β1: 68.8

Teams that have a home attendance of 0 over their season score 572.618 runs. For every additional 1 million fans attending at home games over the season, a team score s68.798 more runs.

Part E

Write out the estimated regression equation.

^Runs scoredi=572.618−68.798 Home attendance (mil)

# if you are using markdown, try out the equatiomatic package
#install.packages("equatiomatic")
library(equatiomatic)
extract_eq(reg, # the regression
           use_coefs = TRUE, # use the estimated numbers
           coef_digits = 3, # how many digits to show
           fix_signs = TRUE) # fix negatives

## $$
## \text{runs\_scored} = 572.618 + 68.798(\text{home\_attend\_mil}) + \epsilon
## $$

Part F

Make a regression table of the output.

# load huxtable
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

## The following object is masked from 'package:pander':
## 
##     wrap

huxreg(reg, # this is sufficient, the rest is decoration
       coefs = c("Constant" = "(Intercept)",
                 "Home Attendance (Millions)" = "home_attend_mil"),
       statistics = c("N" = "nobs",
                      "R-Squared" = "r.squared",
                      "SER" = "sigma"),
       number_format = 3)

Part G

Now let’s start running some diagnostics of the regression. Make a histogram of the residuals. Do they look roughly normal?

# here we need broom's augment() command to add residuals to the data

# load broom
library(broom)
# augment the regression, save as reg_aug
reg_aug<-reg %>%
  augment()

# now we use this as the data in our histogram plot in ggplot, (x is .resid)
ggplot(data = reg_aug)+
  aes(x = .resid)+
  geom_histogram(color="white")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

They do look roughly normal.

Part H

Make a residual plot.

# this is another plot from reg_aug, where x is home attendance and y is .fitted
ggplot(data = reg_aug)+
  aes(x = home_attend_mil,
      y = .resid)+
  geom_point()+
  geom_hline(yintercept=0, color="red")

Part I

Test the regression for heteroskedasticity. Are the errors homoskedastic or heteroskedastic? Generate robust standard errors. Make a regression output table, with one column with regular standard errors and another with robust standard errors.

# this requires the lmtest package for the bptest() command

# install.packages("lmtest")
# load lmtest
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

bptest(reg)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg
## BP = 0.22515, df = 1, p-value = 0.6351

The null hypothesis H0 is that the errors are homoskedastic. The p-value for this test is very large, so we cannot reject the null hypothesis.

This is good, it means the errors are homoskedastic, and our OLS estimators’ standard errors are accurate and do not need to be corrected for heteroskedasticity.

Part J

Test the data for outliers. If there are any, identify which team(s) and season(s) are outliers.

# this requires the car package for the outlierTest() command

# install.packages("car")
# load car

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

outlierTest(reg)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##      rstudent unadjusted p-value Bonferroni p
## 816 -3.255201          0.0011787      0.98779

This test detected one outlier, which is observation (row) number 816. Let’s look it up:

mlb %>%
  slice(816)

The Toronto Blue Jays’ 1981 season is an outlier. Just for kicks, let’s point it out on the scatterplot.

outlier<-mlb %>%
  slice(816)

library(ggrepel)
scatter+ # our scatterplot saved from part c
  geom_point(data = outlier,
             aes(x = home_attend_mil,
                 y = runs_scored),
             color = "red")+
  geom_text_repel(data = outlier,
             aes(x = home_attend_mil,
                 y = runs_scored),
             label = "1981 Blue Jays",
             color = "red")

Part K

What is the marginal effect of home attendance on runs scored? Is this statistically significant? Why or why not?

This is an interpretation question, no need to calculate anything. The marginal effect is ^β1: for every 1 additional million fans attending home games over a team’s season, the team scores 69 more runs.

Looking back at the regression output in part c, the t-score for the hypothesis test H0:β1=0, H1:β1≠0 is 16.45, yielding a very very small p-value. We have sufficient evidence to reject H0 in favor of our alternative hypothesis, that there is a relationship between home attendance and runs scored over a season.

Part L

Now we’ll try out the infer package to understand the p-value and confidence interval for our observed slope in our regression model. Save the (value of) our sample ^β1 from your regression in Part D as an object. Then, install and load the infer package, and then calculate the slope³ under the null hypothesis that there is no connection between attendance and runs.⁴ for 1000 additional simulated samples⁵, and save this as an object (it’s a tibble). Then, use this to get_p_value()⁶. Compare to the p-value given by lm() and tidy() above.

# save slope from our regression as sample_slope

sample_slope<-reg_tidy %>%
  filter(term=="home_attend_mil") %>%
  pull(estimate)

# double check it worked
sample_slope

## [1] 68.79777

# note I could also just use the beta_1_hat I optionally made in Part D

# make 1000 simulations of sample slopes under null hypothesis that slope = 0

#install.packages(infer)
library(infer) # load infer
slope_simulations<-mlb %>% # save as a tibble
  specify(runs_scored ~ home_attend_mil) %>% # our lm model
  hypothesize(null = "independence") %>% # null hypothesis, slope = 0 (X and Y independent)
  generate(reps = 1000, type = "permute") %>% # make 1000 permutations
  calculate(stat = "slope") # calculate the sample slope of each permutation

# make sure it worked
slope_simulations %>%
  head(., n = 10) # there are a LOT! I only print the first 10 for space

slope_simulations %>%
  get_p_value(obs_stat = sample_slope,
              direction = "both")

The p-value is basically 0. According to the regression output in part D, it is smaller than 0.00000000000000002.

Part M

Make a histogram of the simulated slopes, and plot our sample slope on that histogram, shading the p-value.⁷

# now make a histogram of this
slope_simulations %>%
  visualize(obs_stat = sample_slope)+
  shade_p_value(obs_stat = sample_slope, # set the obs_stat equal to our saved slope from above
                direction = "both") # two-sided test, shade both sides

## Warning: `visualize()` should no longer be used to plot a p-value.
## Arguments `obs_stat`, `obs_stat_color`, `pvalue_fill`, and `direction` are
## deprecated. Use `shade_p_value()` instead.

# infer no longer shades the p-value, so I'll do it manually

slope_simulations %>%
  ggplot(data = .)+
  aes(x = stat)+ # slope from simulations
  geom_histogram(color = "white", # the histogram
                 fill = "indianred")+
  # add "shading" for p-value's two sides
  # right side
  geom_rect(xmin=68.798,
            xmax=100,
            ymin=0,
            ymax=450,
            fill = "pink",
            alpha=0.4)+
  # left side
  geom_rect(xmin=-100,
            xmax=-68.798,
            ymin=0,
            ymax=450,
            fill = "pink",
            alpha=0.4)+
  # add vertical line for our sample slope
  geom_vline(xintercept = sample_slope,
             color = "red",
             size = 2)+
  # add label
  geom_label(x = sample_slope,
             y = 200,
             color = "red",
             label = expression(paste("Our ", hat(beta[1]))))+
  scale_x_continuous(breaks=c(-75,-50,-25,0,25,50,75),
                     limits=c(-75,75),
                     expand=c(0,0))+
  scale_y_continuous(limits=c(0,450),
                       expand=c(0,0))+
  labs(x = expression(paste("Distribution of ", hat(beta[1]), " under ", H[0], " that ", beta[1]==0)),
       y = "Samples")+
    theme_classic(base_family = "Fira Sans Condensed",
           base_size=20)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'

Note this is the sampling distribution of ^β1 under the null hypothesis (the true β1=0). Values on the horizontal axis are values of ^β1, not the number of standard deviations away from the null hypothesis.

What the test-statistic does is standardize this distribution much like we would standardize a distribution to the standard normal distribution via calculating Z-scores.

Let me visualize what would happen (and this is what R does to calculate the t value and associationed p-value.)

slope_simulations %>%
  summarize(mean = mean(stat), # get the mean slope of the simulated distribution of slopes
            se = (sd(stat))) # get the standard error of the slope 

# standardize slopes to t-statistics
tstatistics<-slope_simulations %>%
  mutate(tscores = ((stat - mean(stat))/sd(stat)))

our_t<-((sample_slope-mean(slope_simulations$stat))/sd(slope_simulations$stat))

# what is our t-statistic?
our_t

## [1] 14.14881

# now plot t-statistics %>%
ggplot(data = tstatistics)+
  aes(x = tscores)+ # slope from simulations
  geom_histogram(color = "white", # the histogram
                 fill = "indianred")+
  # add vertical line for our sample slope's t-statistic
  geom_vline(xintercept = our_t,
             color = "red",
             size = 2)+
  #add label
  geom_label(x = our_t,
             y = 50,
             color = "red",
             label = "Our t")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

This is off slightly (see the t-statistic R calculated in part C), but it gives you an idea. Our test-statistic is somewhere around 14-16, this means it is 14-16 standard deviations away from 0, so again, a very very low p-value!

Part N

Get the 95% confidence interval for your slope estimate,⁸ and then make a histogram of the simulated slopes (like part L), but instead, add +shade_confidence_interval().⁹ Compare this to what you get with tidy() above.

# save as a tibble called "ci_values"
ci_values<-slope_simulations %>%
  get_confidence_interval(level = 0.95,
                          type = "se",
                          point_estimate = sample_slope)

# see what we made
ci_values

# we'll use this for the endpoints in the shade_ci() command 

slope_simulations %>%
  visualize(obs_stat = sample_slope)+
  shade_confidence_interval(endpoints = ci_values)

## Warning: `visualize()` should no longer be used to plot a p-value.
## Arguments `obs_stat`, `obs_stat_color`, `pvalue_fill`, and `direction` are
## deprecated. Use `shade_p_value()` instead.