14 Dummy Variables Part 2

There are several places in this file where you need to do things (e.g., write R code, write Latex). Some are labeled “YOUR CODE GOES HERE.” Some places tell you to “UN-COMMENT-OUT” code after you’ve done something else (un-comment-out means remove the hash tags from code in code chunks so that it becomes code instead of comments). Some places say “ReplaceWithCoefficientFromBetaModel”. Others are labeled “PLACEHOLDER”. Make sure to go through the file and find all the places you need to do something.

This chapter goes with LN4.2. In LN4.2 we examined two models of the relationship between education and wages. We examined two different ways we could define sets of dummy variables to represent people with a high school degree, a BA, a Master’s, or a PhD as their highest level of education. We saw how both models could be used to answer questions such as “what is the average wage for someone with a Master’s degree?” or “What is the expected difference in average wage for a person with a Master’s degree compared to a BA?” In this BP chapter you will explore these models empirically using micro-level data (i.e., data on individual people).

Specifically, the data we’ll use comes from the 2019 ACS 1-year Public Use Microdata Sample (PUMS) (the 2021 ACS has not been released yet, and the 2020 data hasn’t either/is marked as experimental due to the Pandemic). This data has the responses of individuals that the Census Bureau uses to calculate the county- and state-level statistics they publish as part of the ACS (e.g., what you’re using for the CP). Formally to analyze this data properly we need to use survey weights to make calculations done using the individual responses representative of the population as a whole. We’re not going to worry about that. However, note that our calculations will be off somewhat from the true average wages for the groups we explore in this chapter.

We’ll obtain our data using the tidycensus package, as described here. If you don’t already have a Census API key, get one and use it, as described here. For our purposes, it’s fine to leave census_api_key("YOUR API KEY GOES HERE") in your code.

library(tidyverse)
library(tidycensus)
library(stargazer)
library(pander)
## This turns off scientific notation when there are fewer then 4 digits. Otherwise pander displays very small p values with way too many decimal places
options(scipen = 4)

14.1 Obtain and prepare data

To work with the PUMS data via tidycensus, you can pull the table of variables to find what variables are included. I’ve set this code chunk to not be evaluated (eval = FALSE) when you build because while doing this is helpful while you’re working with this data, it doesn’t belong in your HTML output

## To work with the PUMS data via tidycensus, you can pull the table of variables to find what variables are included
## I've set this code chunk to not be evaluated (eval = FALSE) when you build. You might want to do this on your own, but it doesn't belong in your HTML output
pums_vars_2019 <- pums_variables %>% 
  filter(year == 2019, survey == "acs1")

Let’s pull the following variables from the 2019 ACS-1 year Public Use Microdata (i.e., data on individual people) for the state of Wisconsin (we’re limiting our sample to Wisconsin so that it doesn’t take too long to load or work with).

Variable	Description
`PUMA`	Public Use Microdata Areas
`WAGP`	WAGe of Person
`AGEP`	AGE of Person
`SCHL`	educational attainment of person (i.e., highest level obtained)
`sex`	SEX of person

The results='hide' code chunk option keeps it from displaying status updates while it downloads (which is sometimes about 100 lines of output, so you definitely don’t want it in your HTML output).

## Download 2019 acs1 for Wisconsin.
wiPUMS <- get_pums(
  variables = c("PUMA","WAGP", "AGEP", "SCHL", "SEX"),
  state = "wi",
  survey = "acs1",
  year = 2019,
  recode = TRUE
)

Let’s rename WAGP as wage so that we remember what it is.

wiPUMS <- rename(wiPUMS, wage = WAGP)

Let’s see what levels of education exist in our data, and how many observations we have for each

wiPUMS %>% count(SCHL,SCHL_label) %>% pander()

SCHL	SCHL_label	n
01	No schooling completed	1556
02	Nursery school, preschool	749
03	Kindergarten	731
04	Grade 1	576
05	Grade 2	688
06	Grade 3	674
07	Grade 4	713
08	Grade 5	717
09	Grade 6	785
10	Grade 7	738
11	Grade 8	1215
12	Grade 9	959
13	Grade 10	1169
14	Grade 11	1443
15	12th grade - no diploma	864
16	Regular high school diploma	14484
17	GED or alternative credential	1657
18	Some college, but less than 1 year	3783
19	1 or more years of college credit, no degree	6693
20	Associate’s degree	5281
21	Bachelor’s degree	8476
22	Master’s degree	2971
23	Professional degree beyond a bachelor’s degree	685
24	Doctorate degree	478
bb	N/A (less than 3 years old)	1682

We’re interested in examining how wages differ by educational attainment. We’d like to limit our data to people who are old enough to have completed high school (or equivalent) and college, and then be working. So, let’s limit our sample to people age 25 or older. Let’s also only examine people who have a wage.

wiPUMS <- wiPUMS %>% 
            filter(AGEP >= 25 & wage > 0)

We should look at what education levels (SCHL) remain after filtering

wiPUMS %>% count(SCHL,SCHL_label) %>% pander()

SCHL	SCHL_label	n
01	No schooling completed	133
02	Nursery school, preschool	4
03	Kindergarten	3
04	Grade 1	4
05	Grade 2	3
06	Grade 3	7
07	Grade 4	2
08	Grade 5	3
09	Grade 6	30
10	Grade 7	13
11	Grade 8	97
12	Grade 9	88
13	Grade 10	160
14	Grade 11	250
15	12th grade - no diploma	356
16	Regular high school diploma	6583
17	GED or alternative credential	815
18	Some college, but less than 1 year	1952
19	1 or more years of college credit, no degree	3391
20	Associate’s degree	3700
21	Bachelor’s degree	5887
22	Master’s degree	2106
23	Professional degree beyond a bachelor’s degree	458
24	Doctorate degree	334

After filtering, the value “bb” that indicates a person who is under 3 years old and thus can’t have any schooling is now gone. Now all of the levels of SCHL are numeric, which allows us to convert them from a character to a numeric data type. This will then allow us to use inequalities to define dummy variables (e.g., SCHL <16).

wiPUMS$SCHL <- as.numeric(wiPUMS$SCHL)

We’re comparing wages of people who have a high school degree (or equivalent) or higher, so we also need to drop everyone who has less than a high school degree. Because SCHL is now numeric, we can do this using an inequality (rather than listing all of the options separately).

wiPUMS <- wiPUMS %>% filter(SCHL >= 16)

In order to match with LN4.2, we also need to drop people with a professional degree. In practice we would probably want to examine them too, but for the sake of matching with what is in LN4.2 we’re going to drop them.

wiPUMS <- wiPUMS %>% filter(SCHL != 23)

Whenever you do something like filtering data, it’s a very good idea to look at the data and make sure it worked. In previous years of 380, students have wasted many hours trying to get models to work, only to finally go back and look at the data to find that they actually messed up an earlier step that seems easy. So even if it’s something easy you think you know how to do, look at the data. You can display summary measures as we’ll do here, but it’s also a good idea to click on it in the Environment tab and actually scroll through it quickly. Typically you don’t display the results in a paper, but for the purposes of the BP, I want to demonstrate what you might do (e.g., get a count by education levels). Here, we could just re-run this:

wiPUMS %>% count(SCHL,SCHL_label) %>% pander()

SCHL	SCHL_label	n
16	Regular high school diploma	6583
17	GED or alternative credential	815
18	Some college, but less than 1 year	1952
19	1 or more years of college credit, no degree	3391
20	Associate’s degree	3700
21	Bachelor’s degree	5887
22	Master’s degree	2106
24	Doctorate degree	334

Note that this variable is lists the highest level of education attained. That means it generally builds. We can assume someone who reached a 21 has a 16, and a person who reached a 22 has a 21 (and thus also a 16), and a person who reached 24 has a 22 (and thus also a 21 and 16).

14.2 Define dummy variables

In LN4.2 we had an “alpha model” and “beta model” which defined education in different ways.

Variable	Alpha Model Definition	Beta Model Definition
`BA`	=1 if have B.A. degree, =0 if don’t	=1 if highest degree is B.A. degree, =0 otherwise
`Masters`	=1 if have Master’s degree, =0 if don’t	=1 if highest degree is Master’s, =0 otherwise
`PhD`	=1 if have PhD, =0 if don’t	=1 if highest degree is PhD, =0 otherwise

14.2.1 Alpha Model

In LN4.2, we defined the “Alpha Model” as

$wage = \alpha_0 + \alpha_1 BA + \alpha_2 Masters + \alpha_3 PhD + u$

We need to create the dummy variables used for the “Alpha Model”. We’ll prefix the variables with “a” for “Alpha”. Later you’ll define “Beta Model” variables and prefix them with “b”.

In our data, some people have some college or an Associates Degree. While we might be interested in differences between people with some college and a only a high school degree, for the purposes of what we’re doing here (learning about dummy variables), let’s classify anyone without a BA as having high school as their highest degree (even if they have some college or an Associates Degree).

# Create dummy variables for "alpha model" (starting with a)
wiPUMS <- wiPUMS %>%
            mutate(aBA = ifelse(SCHL >= 21, 1, 0),
                   aMasters = ifelse(SCHL >= 22, 1, 0),
                   aPhD = ifelse(SCHL == 24, 1, 0)
            )

Now let’s estimate the regression shown as “Alpha Model” in LN4

alphaModel <- lm(wage ~ aBA + aMasters + aPhD, data = wiPUMS)
pander(summary(alphaModel))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	42649	373.5	114.2	0
aBA	20853	727.4	28.67	7.28e-178
aMasters	11197	1216	9.208	3.568e-20
aPhD	19184	2821	6.801	1.059e-11

Fitting linear model: wage ~ aBA + aMasters + aPhD
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
24768	47889	0.0637	0.06359

14.2.2 Beta Model

In LN4.2, we defined the “Beta Model” (version 1…we won’t use the “version 2” at all) as:

$wage = \beta_0 + \beta_1 BA + \beta_2 Masters + \beta_3 PhD + u$

We need to create the dummy variables for the “Beta Model” from LN4.2. Use the same categories as you did for the Alpha Model (so we’re considering anyone with some college or an associates degree as having high school as their highest education). Prefix these variables with a “b” for “Beta”.

YOUR CODE GOES HERE: add variables bBA, bMasters, and bPhD to wiPUMS

# Create dummy variables for "Beta Model" (prefix with b)
wiPUMS <- wiPUMS %>%
            mutate(bBA = ifelse(SCHL == 21, 1, 0),
                   bMasters = ifelse(SCHL == 22, 1, 0),
                   bPhD = ifelse(SCHL == 24, 1, 0)
            )

Now let’s estimate the regression shown as “Beta Model” in LN4.2

UN-COMMENT-OUT the following code after you estimate the Beta Model

betaModel <- lm(wage ~ bBA + bMasters + bPhD, data = wiPUMS)
pander(summary(betaModel))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	42649	373.5	114.2	0
bBA	20853	727.4	28.67	7.28e-178
bMasters	32050	1108	28.92	7.426e-181
bPhD	51234	2647	19.36	7.366e-83

Fitting linear model: wage ~ bBA + bMasters + bPhD
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
24768	47889	0.0637	0.06359

14.3 Compare the regressions side-by-side

In LN4.2 we talk about how different models can lead to equivalent results theoretically. Here we want to examine that theoretical equivalence numerically. Specifically, we want to write out various conditional expectations that should be equivalent theoretically and show that the estimated values are indeed equivalent.

Let’s start by displaying the two models on the same table. You have to be very careful doing this (so you don’t accidentally mis-label variables), but we can re-name the variable labels for each model so that they are the same (e.g., rename aBA and bBA to both be BA) and thus display on the same row of the stargazer table. We’ll make a copy of the model results that we use specifically for this purpose. Then we’ll rename the coefficients to have the same names. Then we’ll display them using stargazer.

UN-COMMENT-OUT the following code after you estimate the Beta Model

alphaModelStargazer <- alphaModel
betaModelStargazer <- betaModel
names(alphaModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD")
names(betaModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD")
stargazer(alphaModelStargazer, betaModelStargazer,
          type = "html",
          report=('vc*p'),
          notes = "<em>&#42;p&lt;0.1;&#42;&#42;p&lt;0.05;&#42;&#42;&#42;p&lt;0.01</em>",
          notes.append = FALSE,
          model.numbers = FALSE,
          column.labels = c("Alpha","Beta"))


	Dependent variable:

	wage
	Alpha	Beta

BA	20,853.350^***	20,853.350^***
	p = 0.000	p = 0.000

Masters	11,196.800^***	32,050.150^***
	p = 0.000	p = 0.000

PhD	19,183.600^***	51,233.750^***
	p = 0.000	p = 0.000

Constant	42,649.000^***	42,649.000^***
	p = 0.000	p = 0.000


Observations	24,768	24,768
R²	0.064	0.064
Adjusted R²	0.064	0.064
Residual Std. Error (df = 24764)	47,889.000	47,889.000
F Statistic (df = 3; 24764)	561.593^***	561.593^***

Note:	p<0.1;p<0.05;**p<0.01

14.4 Compare the predictions of each model

In LN4.2 we claimed that the Alpha Model and Beta Model were theoretically equivalent. Here we are demonstrating that they are also empirically equivalent. The first thing you notice is that the $R^2$ and adjusted $R^2$ are identical. We haven’t talked about the other statistics (Residual Std. Error and F Statistic), but they’re identical too. That’s what we would expect if the models are indeed equivalent. Recall the definition of $R^2$ . You’d expect equivalent models to explain the same fraction of the variation in wages.

Now let’s check if real-world situations that are supposed to be equivalent in theory are indeed equivalent numerically.

First, let’s store our coefficients in variables so they’re easier to work with:

## Alpha Model
a0 <- coef(alphaModel)["(Intercept)"]
a1 <- coef(alphaModel)["aBA"]
a2 <- coef(alphaModel)["aMasters"]
a3 <- coef(alphaModel)["aPhD"]

## Beta Model
## FILL THESE IN WITH COEFFICIENTS FROM YOUR BETA MODEL
b0 <- coef(betaModel)["(Intercept)"]
b1 <- coef(betaModel)["bBA"]
b2 <- coef(betaModel)["bMasters"]
b3 <- coef(betaModel)["bPhD"]

## I'll use this to let you know where you need to fill in answers
PLACEHOLDER <- "PLACEHOLDER"

We’ll check a variety of situations explored by these models. For each situation, we’ll write out the following 6 steps (if the steps aren’t clear, also look at the examples below):

First, we’ll write out the situation in English (e.g., “Expected/average wage for a person whose highest degree is high school”)
Next, we’ll write out the conditional expectation corresponding to that situation (e.g., $E(wage|HS)$ ) on the left hand side of the equation. To te right of the | in the E(...|...), write it out as a short description of the situation in the world (e.g., HS) instead of writing out all the variable values (e.g., BA=0, Masters=0, PhD=0).
Then, on the right hand side of the equation, we’ll write out the conditional expectation in terms of model parameters using parameters from the Alpha Model (e.g., $\alpha_0$ )…
…and in terms of model parameters using parameters from the Beta Model (e.g., $\beta_0$ ).
Last, for the third and fourth rows on the right hand side of the equation, we’ll calculate the value estimated using the Alpha Model (e.g., 42649.00)…
…and the value estimated using the Beta Model (e.g., 42649.00).

Note that we’re using format(...,nsmall=2) instead of round(...,2) to display the coefficients to two decimal places.

I’ve filled in a few parts for you below to show you what you’re supposed to do. Anyplace you see “PLACEHOLDER” you need to fill in whatever should go there. (Note that I filled in the coefficient b0 for you, but until you do the parts above, it will display as ReplaceWithCoefficientFromBetaModel.)

Use LaTex for parameters (e.g., $\alpha_0$ , $\beta_0$ ). For values, each value should be displayed to 2 decimal places using format(...,nsmall=2). Note that you can add together multiple coefficients inside the format function (e.g. format(a0 + a1,nsmall=2)). This structure is already set up for you below (i.e., you just need to replace everywhere that says PLACEHOLDER).

When you start filling it in, make sure to build frequently so you don’t mess it up and then are unable to figure out what isn’t working. For the first few, I’d build each time you make a change (remember to set it in _bookdown.yml to only build this chapter while you’re working on it, and then just build the entire book when you’re finished before committing/pushing to GitHub).

For each scenrio, your parameters in rows one and two on the right hand side should match the slides, and the values in rows three and four should be identical. If the values aren’t identical, you did something wrong.

Here are the scenarios:

Expected/average wage for a person whose highest degree is high school (HS) $\begin{aligned} E(wage|HS) &= \alpha_0 \\&= \beta_0 \\&= 42649.00 \\&= 42649.00 \end{aligned}$

Expected/average wage for a person whose highest degree is BA $\begin{aligned} E(wage|BA) &= \alpha_0 + \alpha_1 \\&= \beta_0 + \beta_1 \\&= 63502.35 \\&= 63502.35 \end{aligned}$

Expected/average wage for a person whose highest degree is Master’s $\begin{aligned} E(wage|Masters) &= \alpha_0 + \alpha_1 + \alpha_2 \\&= \beta_0 + \beta_2 \\&= 74699.15 \\&= 74699.15 \end{aligned}$

Expected/average wage for a person whose highest degree is PhD $\begin{aligned} E(wage|PhD) &= \alpha_0 + \alpha_1 + \alpha_2 + \alpha_3 \\&= \beta_0 + \beta_3 \\&= 93882.75 \\&= 93882.75 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a BA compared to someone whose highest degree is high school (HS) $\begin{aligned} E(wage|BA) - E(wage|HS) &= (\alpha_0 + \alpha_1) - \alpha_0 \\&= (\beta_0 + \beta_1) - \beta_0 \\&= 20853.35 \\&= 20853.35 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a Master’s compared to someone whose highest degree is high school (HS) $\begin{aligned} E(wage|Masters) - E(wage|HS) &= \alpha_0 + \alpha_1 + \alpha_2 - (\alpha_0) = \alpha_1 + \alpha_2 \\&= \beta_0 + \beta_2 - (\beta_0) = \beta_2 \\&= 32050.15 \\&= 32050.15 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a PhD compared to someone whose highest degree is high school (HS) $\begin{aligned} E(wage|PhD) - E(wage|HS) &= (\alpha_0 + \alpha_1 + \alpha_2 + \alpha_3) - \alpha_0 \\&= (\beta_0 + \beta_3) - \beta_0 \\&= 51233.75 \\&= 51233.75 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a Master’s compared to someone whose highest degree is a BA $\begin{aligned} E(wage|Masters) - E(wage|BA) &= (\alpha_0 + \alpha_1 + \alpha_2) - (\alpha_0 + \alpha_1) \\&= (\beta_0 + \beta_2) - (\beta_0 + \beta_1) \\&= 11196.80 \\&= 11196.80 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a PhD compared to someone whose highest degree is a BA $\begin{aligned} E(wage|PhD) - E(wage|BA) &= (\alpha_0 + \alpha_1 + \alpha_2 + \alpha_3) - (\alpha_0 + \alpha_1) \\&= (\beta_0 + \beta_3) - (\beta_0 + \beta_1) \\&= 30380.40 \\&= 30380.40 \end{aligned}$

How much higher do you expect the average wage to be for someone whose highest degree is a PhD compared to someone whose highest degree is a Master’s $\begin{aligned} E(wage|PhD) - E(wage|Masters) &= (\alpha_0 + \alpha_1 + \alpha_2 + \alpha_3) - (\alpha_0 + \alpha_1 + \alpha_2) \\&= (\beta_0 + \beta_3) - (\beta_0 + \beta_2) \\&= 19183.60 \\&= 19183.60 \end{aligned}$

14.5 Other things you should think about (and possibly do)

Doing the things below (and doing them correctly) is required to get a check plus. If you don’t do them, you can still get a check.

Remember that a check plus is a 100 and a check is a 93, so the difference is small. Only do the following if you actually have time (e.g., the priority is doing well on the group projects). But if you have time, do them. It will help you better understand.

14.5.1 Group averages

Linear regression when all explanatory variables are dummy variables produces the same estimates of average values for different groups as if you simply took the average of the $y$ value for each group. Try confirming this using group_by() and summarize() to get average wage for each education level (make sure not to save wiPUMS as grouped it still works as expected below), and then compare them with what we found above in the regressions (i.e., make sure they are the same).

wiPUMS %>%
  group_by(SCHL) %>%
    summarize(averagewage = mean(wage)) %>%
      filter(SCHL %in% c(18, 21, 22, 24))

## # A tibble: 4 × 2
##    SCHL averagewage
##   <dbl>       <dbl>
## 1    18      42821.
## 2    21      63502.
## 3    22      74699.
## 4    24      93883.

They are the same values as the estimates from the regression.

Why then do we use regression? One reason is because we easily get standard errors, confidence intervals, and p-values for hypothesis tests of the group averages and differences between groups. Another reason is because we can also control for other variables (such as age…see below).

14.5.2 Causal estimates?

Make sure you can explain why we should not interpret our results as causal. For example, above we estimated that the average wage is 20853.35 higher for someone whose highest degree is a BA compared to someone whose highest degree is high school (HS). Why can we not claim that getting a BA causes this increase in wage? In other words, why is $\hat{\alpha}_1$ (or $\hat{\beta}_1$ ) not an unbiased estimate of the true effect of getting a BA on wages?

There are many other factors that have an effect on wage and level of education, including age. This model likely fails the ZCM assumption.

14.5.3 What about age?

How do the results change if you control for age (i.e., add age as a variable)? You can use either the alpha or beta model and display the results side-by-side (i.e., the original model and the new model including age) using stargazer.

ageModel <- lm(wage ~ bBA + bMasters + bPhD + AGEP, data = wiPUMS)

ageModelStargazer <- ageModel
betaModelStargazer <- betaModel
names(ageModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD", "Age")
names(betaModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD")
stargazer(ageModelStargazer, betaModelStargazer,
          type = "html",
          report=('vc*p'),
          notes = "<em>&#42;p&lt;0.1;&#42;&#42;p&lt;0.05;&#42;&#42;&#42;p&lt;0.01</em>",
          notes.append = FALSE,
          model.numbers = FALSE,
          column.labels = c("Age","Beta"))


	Dependent variable:

	wage
	Age	Beta

BA	20,892.150^***	20,853.350^***
	p = 0.000	p = 0.000

Masters	32,062.110^***	32,050.150^***
	p = 0.000	p = 0.000

PhD	51,229.020^***	51,233.750^***
	p = 0.000	p = 0.000

Age	10.522
	p = 0.653

Constant	42,134.440^***	42,649.000^***
	p = 0.000	p = 0.000


Observations	24,768	24,768
R²	0.064	0.064
Adjusted R²	0.064	0.064
Residual Std. Error	47,889.770 (df = 24763)	47,889.000 (df = 24764)
F Statistic	421.232^*** (df = 4; 24763)	561.593^*** (df = 3; 24764)

Note:	p<0.1;p<0.05;**p<0.01

For the model with age, some of the effects of the other variables decrease, like with the effect of PhD, while BA and Master’s go up slightly. This is somewhat interesting because adding additional explanatory variables to a model generally causes a decrease in the magnitude of effect that the other variables have. In this case, there was very little change.

14.5.4 What about sex?

How might you add sex to the model?

In the data, male is coded as sex=1 and female is coded as sex=2 (see the PUMS data dictionary). Create a female dummy variable (i.e., female=1 indicates female, female=0 indicates male because that’s the only other option in the data).

wiPUMS <- wiPUMS %>%
            mutate(female = ifelse(SEX == 2, 1, 0))

Estimate two models (explained below) and display the results side-by-side using stargazer, along with the original model for comparison (either alpha or beta).

For the first model, simply add female as a variable.

female1 <- lm(wage ~ bBA + bMasters + bPhD + female, data = wiPUMS)
pander(female1)

Fitting linear model: wage ~ bBA + bMasters + bPhD + female
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	52574	456.6	115.2	0
bBA	22498	710.5	31.66	9.491e-216
bMasters	34998	1084	32.3	3.083e-224
bPhD	51720	2580	20.05	1.127e-88
female	-21466	595.8	-36.03	4.273e-277

For the second model, allow how education affects wages to be different for males versus females.

wiPUMS <- wiPUMS %>%
  mutate(bBAf = bBA * female, bMastersf = bMasters * female, bPhDf = bPhD * female)

female2 <- lm(wage ~ bBA + bMasters + bPhD + female +  bBAf + bMastersf + bPhDf, data = wiPUMS)
pander(female2)

Fitting linear model: wage ~ bBA + bMasters + bPhD + female + bBAf + bMastersf + bPhDf
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	50550	495.4	102	0
bBA	28382	1022	27.77	4.02e-167
bMasters	44678	1679	26.61	7.555e-154
bPhD	59331	3586	16.55	3.684e-61
female	-17089	728.6	-23.45	2.446e-120
bBAf	-11539	1419	-8.131	4.45e-16
bMastersf	-17143	2196	-7.807	6.093e-15
bPhDf	-15895	5151	-3.086	0.002033

female1ModelStargazer <- female1
female2ModelStargazer <- female2
betaModelStargazer <- betaModel

names(female1ModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD", "Female")
names(female2ModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD", "Female", "BA * Female", "Masters * Female", "PhD * Female")
names(betaModelStargazer$coefficients) <- c("(Intercept)", "BA", "Masters", "PhD")
stargazer(female1ModelStargazer, female2ModelStargazer, betaModelStargazer,
          type = "html",
          report=('vc*p'),
          notes = "<em>&#42;p&lt;0.1;&#42;&#42;p&lt;0.05;&#42;&#42;&#42;p&lt;0.01</em>",
          notes.append = FALSE,
          model.numbers = FALSE,
          column.labels = c("Female 1", "Female 2", "Beta"))


	Dependent variable:

	wage
	Female 1	Female 2	Beta

BA	22,497.730^***	28,381.850^***	20,853.350^***
	p = 0.000	p = 0.000	p = 0.000

Masters	34,998.190^***	44,678.180^***	32,050.150^***
	p = 0.000	p = 0.000	p = 0.000

PhD	51,719.950^***	59,330.580^***	51,233.750^***
	p = 0.000	p = 0.000	p = 0.000

Female	-21,466.070^***	-17,088.530^***
	p = 0.000	p = 0.000

BA * Female		-11,539.230^***
		p = 0.000

Masters * Female		-17,143.450^***
		p = 0.000

PhD * Female		-15,895.490^***
		p = 0.003

Constant	52,574.500^***	50,550.400^***	42,649.000^***
	p = 0.000	p = 0.000	p = 0.000


Observations	24,768	24,768	24,768
R²	0.110	0.114	0.064
Adjusted R²	0.110	0.114	0.064
Residual Std. Error	46,682.100 (df = 24763)	46,578.050 (df = 24760)	47,889.000 (df = 24764)
F Statistic	767.763^*** (df = 4; 24763)	456.936^*** (df = 7; 24760)	561.593^*** (df = 3; 24764)

Note:	p<0.1;p<0.05;**p<0.01

What does each model allow you to say about differences in wages by sex? In your answer, include specific examples using your coefficients.

The Female 1 Model does not allow for differences in the effect being a female versus being a male has on wage across difference education levels. The second female model does, however. For example, in Model 1, the difference between a female with a BA and a male with a BA is -21466.07, while the difference between a female with a Masters and a male with a Masters is also -21466.07. In Model 2, though, the difference between a female with a BA and a male with a BA is -11539.23, while the difference between a female with a Masters and a male with a Masters is -17143.45. It appears the largest disparity between male’s and female’s wages is at the Master’s level, and the smallest dipsarity is at the Bachelor’s level.