Part A

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

# I just want to show you how to extract these and put it into text with markdown
# This relies on the code in Question 6 Part B
library(broom)
reg_tidy<-reg %>%
  tidy()

beta_0_hat<-reg_tidy %>% # save this as beta_0_hat
  slice(1) %>% # look only at first row (intercept)
  pull(estimate) %>% # extract the value of the estimate
  round(.,2) # round to 2 decimal places

beta_0_hat # look at it

## [1] 0.56

^β0= 0.56, meaning the rate at which Members of Congress agree with the President when the Net Trump Vote was 0 (i.e. both candidates got 50%) is 0.56 (out of 1).¹

Part B

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

beta_1_hat<-reg_tidy %>% # save this as beta_0_hat
  slice(2) %>% # look only at second row (net_trump_score)
  pull(estimate) %>% # extract the value of the estimate
  round(.,2) # round to 2 decimal places

beta_1_hat # look at it

## [1] 0.01

^β1= 0.01, meaning for every additional Net Trump Vote % in a Member of Congress’ district, that Member will agree 0.01 more with the President (out of 1).

Part C

What is R2? What does it mean?

# this uses the glance command in question 6 part C. 

r_sq<-glance(reg)$r.squared %>%
  round(.,2) # round to 2 decimals

r_sq # look at it

## [1] 0.61

R2= 0.61, meaning we have explained 61%² of the variation in agree_pct.

Part D

What is the SER? What does it mean?

# this uses the glance command in question 6 part C. 

ser<-glance(reg)$sigma %>%
  round(.,2) # round to 2 decimals

ser # look at it

## [1] 0.25

SER= 0.25, meaning the average size of the residuals or errors (distance from an actual data point to the predicted regression line) is 0.25 (out of 1).

Part A

Install and then load broom.

#install.packages("broom")
library(broom)

Part B

Run the function tidy() on your regression object (saved in question 5). Save this as an object and then look at it.

reg_tidy<-reg %>%
  tidy()
reg_tidy

Part C

Run the glance() function on your original regression object. What does it show you?

reg %>%
  glance()

It shows you summary statistics of the regression, particularly its goodness of fit.

• r.squared is the R2 value
• sigma is the SER(^σu)

That’s all we need for now.

Part D

Now run the augment() function on your original regression object, and save this as an object. Look at it. What does it show you?

reg_aug<-reg %>%
  augment()
reg_aug

It takes the original data points for Y and X and adds new variables with values for those observations: - .fitted is ^Yi|X, the predicted value for a specific net_trump_vote - .resid is ^ui, the residual or error for a specific net_trump_vote ^ui=Yi−^Yi.resid=agree_pct−.fitted

That’s all we need right now.

Part A

Take the augmented regression object from Question 6-D and use it as the source of your data to create a histogram, where x is .resid. Does it look roughly normal?

ggplot(data = reg_aug)+
  aes(x = .resid)+
  geom_histogram(color="white")+
  theme_classic()

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

It does look rougly normal.

Part B

Take the augmented regression object and make a residual plot, which is a scatterplot where x is the normal x variable, and y is the .resid. Feel free to add a horizontal line at 0 with geom_hline(yintercept=0).

ggplot(data = reg_aug)+
  aes(x = net_trump_vote, y=.resid)+
  geom_point()+
  geom_hline(yintercept = 0, color = "red")+
  theme_classic()

Part A

Looking at the scatterplot and residual plots in Questions 3 and 7B, do you think the errors are heteroskedastic or homoskedastic?

It does seem like it would be, looking at the residual plot, since there are clear downward-moving clusters of points that are at times above the line, then below the line, then above and below the line.

Part B

Install and load the lmtest package and run bptest on your regression object. Was the data heteroskedastic or homoskedastic?

#install.packages("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 = 6.0092, df = 1, p-value = 0.01423

Since the p-value was less than 0.05, we can reject the null hypothesis H0 (that the errors are homoskedastic), and conclude that the errors are heteroskedastic. We will need to fix them if we want accurate standard errors.

Part C

Now let’s make some heteroskedasticity-robust standard errors. Install and load the estimatr package and use the lm_robust() command (instead of the lm() command) to run a new regression (and save it). Make sure you add se_type="stata" inside the command to calculate robust SEs. Look at it. What changes?

#install.packages("estimatr")
library(estimatr)

reg_rse<-lm_robust(agree_pct ~ net_trump_vote, data = congress, se_type="stata")
reg_rse

##                   Estimate   Std. Error  t value Pr(>|t|)    CI Lower
## (Intercept)    0.561645133 0.0060117565 93.42447        0 0.549854072
## net_trump_vote 0.009921844 0.0001493938 66.41403        0 0.009628833
##                  CI Upper   DF
## (Intercept)    0.57343619 1733
## net_trump_vote 0.01021486 1733

The standard errors on both estimates (^β0 and ^β1) change, but the estimates themselves do not change!

Note t values and Pr(>|t|) will change, and this command adds confidence intervals (CI) and degrees of freedom for t, DF, but we have not covered any of these yet!

Part D

Now let’s see this in a nice regression table. Use huxreg() again, but add both your original regression and your regression saved in part C. Notice any changes?

huxreg("Original" = reg, "Robust SEs" = reg_rse,
       number_format = 6) # need to see more decimal places!

Note that the standard errors are too small to register anything within 3 decimal places (the default). If we were to make a true percent variable out of 100% instead of out of 1 (as I have you do in Question 13), we would see more of a difference. Look carefully, the standard error on net_trump_vote is 0.000192 (from Question 5), and the robust standard error is 0.000149 from Part C.

Part A

Just looking at the scatterplot in Question 3, do you see any outliers?

It looks like there might be one - for the Member whose district’s Net vote was somewhere around −60 and s/he agrees with the President about 0.6.

Part B

Install and load the car package. Run the outlierTest command on your regression object. Does it detect any outliers?

#install.packages("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
## 1708 -2.886125          0.0039485           NA

Yes, surprisingly it’s not the point that looks like it stands out on the graph (middle-left), it’s a different point.

Part C

Look in your original data to match this outlier with an observation. Hint: use the slice() command, as the outlier test gave you an observation (row) number!

congress %>%
  slice(1708)

Part A

First, the variable congress (session of Congress) seems a bit off. Get a count() of congress.

congress %>%
  count(congress)

Part B

Let’s get rid of the 0 values for congress (someone made a mistake coding this, probably). Also, while we’re at it, let’s take agree_pct and mutate a variable that is a proper percentage (i.e. *100).

congress_tidy<-congress %>%
  filter(congress!=0) %>%
  mutate(agree_pct = agree_pct * 100)

Part C

The variable party is also quite a mess. count() by party to see. Then let’s mutate a variable to make a better measure of political party - just "Republican", "Democrat", and "Independent". Try doing this with the case_when() command (as your mutate formula).³ When you’re done count() by your new party variable to make sure it worked.

congress_tidy <- congress_tidy %>%
  mutate(pol_party = case_when(
    party %in% c("Democrat", "D") ~ "Democrat",
    party %in% c("Republican", "R") ~ "Republican",
    party %in% c("Independent", "I") ~ "Independent"
  ))

congress_tidy %>%
  count(pol_party)

Part D

Now plot a scatterplot (same as Question 3) and set color to your party variable. Notice R uses its own default colors, which don’t match to the actual colors these political parties use! Make a vector where you define the party colors as follows: party_colors <- c("Democrat" = "blue", "Republican" = "red", "Independent" = "gray"). Then, run your plot again, adding the following line to customize the colors +scale_colour_manual("Parties", values = party_colors).⁴

party_colors <- c("Democrat" = "blue", "Republican" = "red", "Independent" = "gray")

p<-ggplot(data = congress_tidy)+
  aes(x = net_trump_vote, y=agree_pct)+
  geom_point(aes(color = pol_party))+
  geom_smooth(method="lm", color = "black")+
  scale_colour_manual("Parties", values = party_colors)+
  theme_light()
p

Part E

Now facet your scatterplot by chamber.

p+facet_wrap(~chamber)

# just for kicks
p+facet_grid(chamber~pol_party)