Moderation and Heterogeneous Effects

Introduction!

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

During this week’s lecture you were introduced to Moderation and Heterogeneous Effects.

In this lab session we will:

• Learn how to perform interaction models with lm()
• Learn how to extract marginal/partial effects with margins::margins() and predictive margins with ggeffects::::ggeffect()
• Learn how to vectorize multiple ifelse() statements with dplyr::case_when()

Packages

set.seed(42) #for consistent results

library(dplyr) # to wrangle our data
library(tidyr) # to wrangle our data - pivot_longer()
library(ggplot2) # to render our graphs
library(readr) # for loading the .csv data
library(janitor) # for data management and tabyl()
library(kableExtra) # to render better formatted tables
library(stargazer) # for formatting your model output

library(margins) #for calculating MARGINAL/PARTIAL EFFECT
library(ggeffects) # easily calculating PREDICTIVE MARGINS

Let’s explore who is in our sample

We can use what we have learned about the janitor::tabyl() function, to check who was in our sample:

moderation_df %>% 
  janitor::tabyl(treatment) %>% 
  knitr::kable(col.names = c("Treatment", "N", "Percent, %")) %>% # create kable table
  kableExtra::kable_styling() # view kable table

moderation_df %>% 
  janitor::tabyl(treatment, victim_verbose) %>% 
  janitor::adorn_totals(c("row", "col")) %>%
  knitr::kable(col.names = c("Treatment", "Not a victim", "Victim", "Total")) %>% # create kable table
  kableExtra::kable_styling() %>% # view kable table
  kableExtra::add_header_above(c("", "Victimization status" = 2, ""))

Let’s explore visually the observed levels of public support for the agreement

ggplot(moderation_df, aes(x = support_thermometer)) +
  geom_density(fill = "#af8dc3", alpha = 0.5) +
  theme_minimal() +
  labs(x = "Support thermometer",
       y = "Density")

What do we see in this graph?

Let’s explore visually the observed levels of public support for the agreement conditional on the treatment status

ggplot(moderation_df, aes(x = support_thermometer, fill = factor(treatment))) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(x = "Support thermometer",
       y = "Density",
       fill = "Treatment") +
  theme(legend.position = "bottom") +
  scale_fill_manual(name = " ", # changes to fill dimension
                     values = c("#a7a8aa", "#cc0055"),
                     labels = c("Control", "Treatment")) 

What do we see through these distributions?

Let’s explore visually the observed levels of public support for the agreement conditional on the treatment and victimhood status

ggplot(moderation_df, aes(x = support_thermometer, fill = factor(treatment))) +
  geom_density(alpha = 0.5) +
  facet_grid(treatment~victim_verbose) +
  theme_bw() +
  labs(x = "Support thermometer",
       y = "Density",
       fill = "Treatment") +
  theme(legend.position = "bottom") +
  scale_fill_manual(name = " ", # changes to fill dimension
                     values = c("#a7a8aa", "#cc0055"),
                     labels = c("Control", "Treatment")) 

What patterns do we see here?

How to estimate heterogenous treatment effects

Heterogeneous treatment effects are usually estimated with regression models that include an interaction between the treatment and the moderator. In our case, the formula would look like this:

\[Y_{i} = β_0 + β_1D_{i} + β_2Victim_{i} + β_3D_i * Victim_{i}+ ϵ_i\]

• \(β_0\): Constant
• \(β_1\): Effect of \(D_i\) on \(Y_i\) if \(Victim_i\) is zero
• \(β_2\): Effect of \(Victim_i\) on \(Y_i\) if \(D_i\) is zero
• \(β_3\): Difference in treatment effects of \(D_i\) depending on \(Victim_i\)

or in lm() terms (same result):

lm(outcome ~ treatment + moderator + (treatment*moderator))
lm(outcome ~ treatment + moderator + treatment:moderator)
lm(outcome ~ treatment * moderator)

Think of the switch logic posed by Prof. Munzert

Let’s model our results

We will move forward by creating two models:

• A naive model, where we will regress support_thermometer on treatment.
• An interaction model, where we will add an interaction between treatment and victim

Naive modeling

naive_model <- lm(support_thermometer ~ treatment, data = moderation_df)

stargazer::stargazer(naive_model, type = "html")

What does this model tell us?

The naive model suggests that the subjects that the support for the peace agreement is about 18.9 percentage points higher on average for the subjects that watched the educational video. We suspect, however, that the video may affect differently victims from non-victims of the conflict.

ggplot(moderation_df, aes(x = support_thermometer, fill = factor(treatment))) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  geom_vline(xintercept = 35.055, linetype = "longdash", color = "#a7a8aa") + #D=0 (just beta0)
  geom_vline(xintercept = 35.055 + 18.898, linetype = "longdash", color = "#cc0055") + #D=1 (beta0+beta1) +
  geom_text(aes(label = "ß0", x = 35.055 + 3, y = 0.04), color = "#a7a8aa") + # we add the 3 to repel from the line
  geom_text(aes(label = "ß0 + ß1", x = 35.055 + 18.898 + 6, y = 0.04 ), color = "#cc0055") + # we add the 6 to repel from the line
  labs(x = "Support thermometer",
       y = "Density",
       fill = "Treatment") +
  theme(legend.position = "bottom") +
  scale_fill_manual(name = " ", # changes to fill dimension
                     values = c("#a7a8aa", "#cc0055"),
                     labels = c("Control", "Treatment")) 

Interaction model

Following any of the different formats renders the same results.

interaction_model <- lm(support_thermometer ~ treatment + victim + (treatment*victim), data = moderation_df)
interaction_model_2 <- lm(support_thermometer ~ treatment + victim + treatment:victim, data = moderation_df)
interaction_model_3 <- lm(support_thermometer ~ treatment*victim, data = moderation_df)

stargazer::stargazer(interaction_model, interaction_model_2, interaction_model_3, type = "html")

What does this model tell us?

• \(β_0\): Constant = The average support for non-victims who were not exposed to the video was 27.92
• \(β_1\): Effect of \(D_i\) on \(Y_i\) if \(Victim_i\) is zero = The educational video results in an increase of about 8 percentage points of the peace agreement for the non-victims
• \(β_2\): Effect of \(Victim_i\) on \(Y_i\) if \(D_i\) is zero = On average, victims’ support for the peace agreement is 14.3 percentage points higher than that of the non-victims in the control group
• \(β_3\): Difference in treatment effects of \(D_i\) depending on \(Victim_i\) = The educational video results in an additional increase for victims of about 21.8 percentage points, compared to the effect for non-victims

moderation_df %>%
  dplyr::group_by(treatment,victim_verbose) %>%
  dplyr::mutate(avg_support = mean(support_thermometer)) %>%
ggplot(., aes(x = support_thermometer, fill = factor(treatment))) +
  geom_density(alpha = 0.5) +
  geom_vline(aes(xintercept = avg_support), linetype = "longdash") +
  facet_grid(treatment~victim_verbose) +
  theme_bw() +
  labs(x = "Support thermometer",
       y = "Density",
       fill = "Treatment") +
  theme(legend.position = "bottom") +
  scale_fill_manual(name = " ", # changes to fill dimension
                     values = c("#a7a8aa", "#cc0055"),
                     labels = c("Control", "Treatment"))

Let’s plot the marginal/partial effect of the treatment for victims and non-victims

pe_margins <- margins::margins(interaction_model, variables = "treatment", at = list(victim = c(0,1)))

pe_plotting <- summary(pe_margins) %>% #NOTE we use the summary output, instead of the margins object
  dplyr::select(victim, AME, lower, upper) %>% # you will need to adapt this based your moderator
  dplyr::mutate(victim_labels = ifelse(victim == 1, "Victim", "No Victim")) # you may not need this to create labels for the assignment

ggplot(pe_plotting, aes(x = victim_labels,
                        y = AME)) +
  geom_point(size = 1.5) +
  geom_segment(aes(x = victim_labels, xend = victim_labels, y = lower, yend = upper), size = 0.5) + # render whiskers from confidence intervals
  theme_minimal() +
  scale_y_continuous(limits = c(0,45)) + # you may need to change the limits for your plots based on the specific effects of your application
  labs(x = "Respondent status\n",
       y = "\nPartial effect of educational video", 
       caption = "Note: You can utilize it to describe what the plot illustrates during your assignment.") + 
  coord_flip() # to flip the plot

Let’s plot these predictive margins

In this exercise we will meet a very important function from dplyr, dplyr::case_when(). This function allows us to vectorize multiple ifelse() statements. The syntax is the following:

dplyr::case_when(condition ~ what to do if met)

Let’s see it at play.

extracted_me <-  ggeffects::ggeffect(interaction_model, terms = c("victim","treatment")) %>%
  dplyr::mutate(group_labels = dplyr::case_when(x == 1 & group == 1 ~ "Victim (1) - Treatment (1)",
                                                x == 1 & group == 0 ~ "Victim (1) - Treatment (0)",
                                                x == 0 & group == 1 ~ "Victim (0) - Treatment (1)",
                                                x == 0 & group == 0 ~ "Victim (0) - Treatment (0)"
  )) # adding labels to (x,group) combinations for the plot

extracted_me

ggplot(extracted_me, aes(x = group_labels, 
                         y= predicted, 
                         fill = x,
                         alpha = group
                         )) +
  geom_bar(stat = "identity") + 
  geom_point(size = 1.5) +
  geom_segment(aes(x = group_labels, 
                   xend = group_labels, 
                   y = conf.high, 
                   yend = conf.low), size = 0.5) + # render whiskers from confidence intervals
  theme_minimal() +
  labs(x = "\nRespondent status",
       y = "Predictive margins",
       fill = "Treatment",
       caption = "Note: You can utilize it to describe what the plot illustrates during your assignment.") +
  theme(legend.position = "note") +
  scale_alpha_manual(values=c(0.6, 1))