The Backdoor Criterion and Basics of Regression in R

Introduction!

Welcome to our fourth tutorial for the Statistics II: Statistical Modeling & Causal Inference (with R) course.

During this week's lecture you reviewed bivariate and multiple linear regressions. You also learned how Directed Acyclic Graphs (DAGs) can be leveraged to gather causal estimates.

In this lab session we will:

• Use the ggdag and dagitty packages to assess your modeling strategy
• Review how to run regression models using R
• Illustrate omitted variable and collider bias

# These are the libraries we will use today. Make sure to install them in your console in case you have not done so previously.
set.seed(42) #for consistent results in randomization
library(wooldridge) # To get our example's dataset 
library(tidyverse) # To use dplyr functions and the pipe operator when needed
library(ggplot2) # To visualize data (this package is also loaded by library(tidyverse))
library(ggdag) # To dagify and plot our DAG objects in R
library(dagitty) # To perform analysis in our DAG objects in R
library(stargazer) # To render nicer regression output
data("wage1") # calls the wage1 dataset from the woorldridge package

Learning about our paths and what adjustments we need

As you have seen, when we dagify a DAG in R a dagitty object is created. These objects tell R that we are dealing with DAGs.

This is very important because in addition to plotting them, we can do analyses on the DAG objects. A package that complements ggdag is the dagitty package.

Today, we will focus on two functions from the dagitty package:

• dagitty::paths(): Returns a list with two components: paths, which gives the actual paths, and open, which shows whether each path is open (d-connected) or closed (d-separated).
• dagitty::adjustmentSets(): Lists the sets of covariates that would allow for unbiased estimation of causal effects, assuming that the causal graph is correct.

We just need to input our DAG object.

dagitty::paths(our_dag)

## $paths
## [1] "d -> y"                "d <- a -> y"           "d <- p -> b <- c -> y"
## 
## $open
## [1]  TRUE  TRUE FALSE

ggdag::ggdag_paths(our_dag) +
  theme_dag()

dagitty::adjustmentSets(our_dag)

## { a }

ggdag::ggdag_adjustment_set(our_dag, shadow = T) +
  theme_dag()

Bivariate regression

In bivariate regression, we are modeling a variable \(y\) as a mathematical function of one variable \(x\). We can generalize this in a mathematical equation as such:

\[y = \beta_{0} + \beta{1}x + ϵ\]

Multiple linear regression

In multiple linear regression, we are modeling a variable \(y\) as a mathematical function of multiple variables \((x, z, m)\). We can generalize this in a mathematical equation as such:

\[y = \beta_{0} + \beta_{1}x + \beta_{2}z + \beta_{3}m + ϵ\]

Visualizing

With regression we can answer EXPLORATORY QUESTIONS. For example:

What is the relationship between education and respondents' salaries?

We can start by exploring the relationship visually with our newly attained ggplot2 skills:

ggplot(wage1, aes(x = educ, y = wage)) +
  geom_point(color = "grey60") +
  geom_smooth(method = "lm", se = F, color = "#CC0055") +
  theme_minimal() +
  labs(x = "Years of education",
       y = "Hourly wage (USD)") 

The lm() function

This question can be formalized mathematically as:

\[Hourly\ wage = \beta_0 + \beta_1Years\ of\ education + ϵ\]

Our interest here would be to build a model that predicts the hourly wage of a respondent (our outcome variable) using the years of education (our explanatory variable). Fortunately for us, R provides us with a very intuitive syntax to model regressions.

The general syntax for running a regression model in R is the following:

your_model_biv <- lm(outcome_variable ~ explanarory_variable, data = your_dataset) #for a bivariate regression
your_model_mult <- lm(outcome_variable ~ explanarory_variable_1 + explanarory_variable_2, data = your_dataset) #for multiple regression

Now let's create our own model and save it into the model_1 object, based on the bivariate regression we specified above in which wage is our outcome variable, educ is our explanatory variable, and our data come from the wage1 object:

model_1 <- lm(wage ~ educ, data = wage1)

summary() and broom::tidy()

We have created an object that contains the coefficients, standard errors and further information from your model. In order to see the estimates, you could use the base R function summary(). This function is very useful when you want to print your results in your console.

Alternatively, you can use the tidy() function from the broom package. The function constructs a data frame that summarizes the model’s statistical findings. You can see what else you can do with broom by running: vignette(“broom”). The broom::tidy() function is useful when you want to store the values for future use (e.g., visualizing them)

Let's try both options in the console up there. You just need to copy this code below the model_1 code.

summary(model_1)
broom::tidy(model_1)

Exercise

\[Hourly\ wage = \beta_0 + \beta_1Years\ of\ education + ϵ\]

How would you interpret the results of our model_1?

• What does the constant mean?
• What does the educ coefficient mean?
• Do these coefficient carry any causal meaning?

Let's explore this visually

What is the relationship between education and respondents' salaries conditional on the sex of the worker?

ggplot(wage1, aes(x = educ, y = wage, color = as.factor(female))) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = F) +
  theme_minimal() +
  labs(x = "Years of education",
       y = "Hourly wage (USD)",
       color = "Female") +
  theme(legend.position = "bottom")

Check what happens when we replace the color = as.factor(female) for color = female

• What insights can we gather from this graph?

Multiple linear regression

This question can be formalized mathematically as:

\[Hourly\ wage = \beta_0 + \beta_1Years\ of\ education + \beta_2Female + ϵ\]

Our interest here would be to build a model that predicts the hourly wage of a respondent (our outcome variable) using the years of education and their sex (our explanatory variables).

Let's remember the syntax for running a regression model in R:

your_model_biv <- lm(outcome_variable ~ explanarory_variable, data = your_dataset) #for a bivariate regression
your_model_mult <- lm(outcome_variable ~ explanarory_variable_1 + explanarory_variable_2, data = your_dataset) #for multiple regression

Now let's create our own model, save it into the model_2 object, and print the results based on the formula regression we specified above in which wage is our outcome variable, educ and female are our explanatory variables, and our data come from the wage1 object:

model_2 <- lm(wage ~ educ + female, data = wage1)
summary(model_2)

## 
## Call:
## lm(formula = wage ~ educ + female, data = wage1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9890 -1.8702 -0.6651  1.0447 15.4998 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.62282    0.67253   0.926    0.355    
## educ         0.50645    0.05039  10.051  < 2e-16 ***
## female      -2.27336    0.27904  -8.147 2.76e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.186 on 523 degrees of freedom
## Multiple R-squared:  0.2588, Adjusted R-squared:  0.256 
## F-statistic: 91.32 on 2 and 523 DF,  p-value: < 2.2e-16

Exercise

\[Hourly\ wage = \beta_0 + \beta_1Years\ of\ education + \beta_2Female + ϵ\]

How would you interpret the results of our model_2?

• What does the constant mean?
• What does the educ coefficient mean?
• What does the female coefficient mean?

Predicting from our models

As we discussed previously, when we do not have our causal inference hats on, the main goal of linear regression is to predict an outcome value on the basis of one or multiple predictor variables.

R has a generic function predict() that helps us arrive at the predicted values on the basis of our explanatory variables.

The syntax of predict() is the following:

predict(name_of_the_model, newdata = data.frame(explanatory1 = value, explanatory2 = value))

Say that based on our model_2, we are interested in the expected average hourly wage of a woman with 15 years of education.

predict(model_2, newdata = data.frame(educ = 15, female = 1))

##        1 
## 5.946237

• What does this result tell us?
• What happens when you change female to 0? What does the result mean?
• Can you think of a way to find the difference in the expected hourly wage between a male with 16 years of education and a female with 17?

predict(model_2, newdata = data.frame(educ = 16, female = 0)) - predict(model_2, newdata = data.frame(educ = 15, female =0))

##         1 
## 0.5064521

True model

true_model_celebrity <- lm(talent ~ beauty, data = celebrity_df)
summary(true_model_celebrity)

## 
## Call:
## lm(formula = talent ~ beauty, data = celebrity_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9225 -0.6588 -0.0083  0.6628  3.5877 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.962209   0.107595  27.531   <2e-16 ***
## beauty      0.006545   0.020755   0.315    0.753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9865 on 998 degrees of freedom
## Multiple R-squared:  9.964e-05,  Adjusted R-squared:  -0.0009023 
## F-statistic: 0.09945 on 1 and 998 DF,  p-value: 0.7526

ggplot(celebrity_df, aes(x=beauty, 
                         y=talent)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  labs(x = "Beauty",
       y = "Talent")

Biased model from previous literature

Let's see:

biased_model_celibrity <- lm(talent ~ beauty + celebrity_status, data = celebrity_df)
summary(biased_model_celibrity)

## 
## Call:
## lm(formula = talent ~ beauty + celebrity_status, data = celebrity_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4244 -0.5394  0.0110  0.5064  2.9429 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.37834    0.13983   38.46   <2e-16 ***
## beauty                        -0.32668    0.02265  -14.43   <2e-16 ***
## celebrity_statusNot Celebrity -1.51375    0.06808  -22.24   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.807 on 997 degrees of freedom
## Multiple R-squared:  0.3316, Adjusted R-squared:  0.3302 
## F-statistic: 247.3 on 2 and 997 DF,  p-value: < 2.2e-16

ggplot(celebrity_df, aes(x=beauty, y=talent, color = celebrity_status)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  labs(x = "Beauty",
       y = "Talent",
       color = "")

As we can see, by controlling for a collider, the previous literature was inducing to a non-existent association between beauty and talent, also known as collider or endogenous bias.

True model

true_model_salary <- lm(salary ~ shoesize + sex, data = salary_df)
summary(true_model_salary)

## 
## Call:
## lm(formula = salary ~ shoesize + sex, data = salary_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2341 -1.3698 -0.0501  1.3595  6.4303 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.73879    1.15886  22.210   <2e-16 ***
## shoesize    -0.02030    0.03044  -0.667    0.505    
## sexMale      5.05924    0.17616  28.720   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.981 on 997 degrees of freedom
## Multiple R-squared:  0.6129, Adjusted R-squared:  0.6121 
## F-statistic: 789.2 on 2 and 997 DF,  p-value: < 2.2e-16

ggplot(salary_df, aes(x=shoesize, y=salary, color = sex)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  labs(x = "Shoe size",
       y = "Salary",
       color = "")

Biased model from previous literature

biased_model_salary <- lm(salary ~ shoesize, data = salary_df)
summary(biased_model_salary)

## 
## Call:
## lm(formula = salary ~ shoesize, data = salary_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8777 -1.9101 -0.0511  1.8496  7.9774 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.68865    1.17280   3.145  0.00171 ** 
## shoesize     0.59429    0.02925  20.319  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.676 on 998 degrees of freedom
## Multiple R-squared:  0.2926, Adjusted R-squared:  0.2919 
## F-statistic: 412.9 on 1 and 998 DF,  p-value: < 2.2e-16

ggplot(salary_df, aes(x=shoesize, y=salary)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  theme_minimal() +
  labs(x = "Shoe size",
       y = "Salary")

As we can see, by failing to control for a confounder, the previous literature was creating a non-existent association between shoe size and salary, incurring in ommited variable bias.