Average Weekly Wage SDID and SC

Functions

The functions to be used for the SDID and SC Estimates were created.

with.overlay <- function(est, s) { attr(est,'overlay') = s; est }

estimators <- function(sdid,sc,s) {
  estimator.list = list(with.overlay(sdid, s), sc)
  names(estimator.list)=c('SDID', 'SC')
  estimator.list
}

plot.estimators <- function(ests, alpha.multiplier) {
  p = synthdid_plot(ests, se.method='none',
                    alpha.multiplier=alpha.multiplier, facet=rep(1,length(ests)),
                    trajectory.linetype = 1, effect.curvature=-.4,
                    trajectory.alpha=.5, effect.alpha=.5, diagram.alpha=1)
  suppressMessages(p + scale_alpha(range=c(0,1), guide='none'))
}

Data Preparation

The data to be used was prepared and all MSA in Louisiana were excluded except the New Orleans MSA

# Excluding all LA MSAs except New Orleans
# then Balancing data
df <- read.csv("Average Weekly Wage.csv") %>%
    mutate(Date = case_when(Qtr == 1 ~ paste0(Year, "-01-01"),
        Qtr == 2 ~ paste0(Year, "-04-01"),
        Qtr == 3 ~ paste0(Year, "-07-01"),
        Qtr == 4 ~ paste0(Year, "-10-01"),
        TRUE ~ "NA"),
        Date = as.Date(Date)) %>%
    filter(State != "LA" | grepl("New Orleans", ignore.case = TRUE, AreaName))

dt <- seq.Date(from = as.Date("1998-01-01"),
    to = as.Date("2015-01-01"), by = "quarter")

df.fix <- data.frame(AreaName = rep(unique(df$AreaName), length(dt))) %>%
    arrange(AreaName) %>%
    mutate(Date = rep(dt, length(unique(df$AreaName))))

# Y Axis Label
dt1 <- df %>%
    distinct(Date, time) %>%
    filter(Date >= as.Date("2000-01-01") & Date <= as.Date("2015-01-01")) %>%
    mutate(month = month(Date),
        time = ifelse(month <= 9, paste0(year(Date), "-0", month), time)) %>%
    filter(month == 1 & year(Date) %% 5 == 0)

Total Employment

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    rename(Total.Covered = Total.Covered.Total..all.industries) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Total.Covered)) & Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Total.Covered

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Total Employment", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl

Goods Producing

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Goods.producing)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Goods.producing

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl1 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Goods Producing", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl1

Natural Resources and Mining

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Natural.resources.and.mining)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Natural.resources.and.mining

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl2 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Natural Resources and Mining", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl2

Construction

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Construction)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Construction

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl3 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Construction", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl3

Manufactruing

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Manufacturing)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Manufacturing

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl4 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Manufacturing", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl4

Service Providing

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Service.providing)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Service.providing

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl5 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Service Providing", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl5

Trade, Transportation and Utilities

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Trade..transportation..and.utilities)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Trade..transportation..and.utilities

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl6 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Trade, Transportation and Utilities", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl6

Financial Activities

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Financial.activities)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Financial.activities

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl7 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"),
        linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date),
        max(dt1$Date), by = "5 year"), labels = dt1$time) +
    labs(title = "",
        y = "Financial Activities", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl7

Professional and Business Services

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Professional.and.business.services)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Professional.and.business.services

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl8 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"), linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date), max(dt1$Date), by = "5 year"),
        labels = dt1$time) +
    labs(title = "",
        y = "Professional and Business Services", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl8

Education and Health Services

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Education.and.health.services)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Education.and.health.services

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl9 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"),
        linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date),
        max(dt1$Date), by = "5 year"), labels = dt1$time) +
    labs(title = "",
        y = "Education and Health Services", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl9

Leisure and Hospitality

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Private.Leisure.and.hospitality)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Private.Leisure.and.hospitality

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl10 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"),
        linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date),
        max(dt1$Date), by = "5 year"), labels = dt1$time) +
    labs(title = "",
        y = "Leisure and Hospitality", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl10

Local Government

bal_df <- left_join(df.fix, df, by = c("AreaName", "Date")) %>%
    mutate(Date = as.Date(Date)) %>%
    group_by(AreaName) %>%
    filter(!any(is.na(Local.Government.Total..all.industries)) &
        Date >= as.Date("2000-01-01")) %>%
    ungroup() %>%
    mutate(treat = ifelse(State == "LA" & Date >= as.Date("2005-07-01"), 1, 0))

bal_df$AreaName <- factor(bal_df$AreaName)
bal_df <- as.data.frame(bal_df)
bal_df$lg <- bal_df$Local.Government.Total..all.industries

setup <- panel.matrices(bal_df,
    unit = "AreaName", time = "Date", outcome = "lg", treatment = "treat")

sdid.df <- synthdid_estimate(setup$Y, setup$N0, setup$T0)
sc.df <- sc_estimate(setup$Y, setup$N0, setup$T0)

p4 <- plot.estimators(estimators(sdid = sdid.df, sc = sc.df, s = 1),
    alpha.multiplier = c(1, 1, 1))

plot.theme <- theme(legend.position = "bottom",
    legend.direction = 'horizontal',
    legend.key = element_blank(),
    legend.background = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12))

pl11 <- p4 + plot.theme +
    geom_vline(xintercept = as.Date("2005-07-01"),
        linetype = "dashed", color = "black") +
    scale_y_continuous(labels = function(x) format(x, scientific = FALSE)) +
    scale_x_continuous(breaks = seq(min(dt1$Date),
        max(dt1$Date), by = "5 year"), labels = dt1$time) +
    labs(title = "",
        y = "Local Government", x = "Date")

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

pl11

All Plots

# Get the legends from each plot
legends <- lapply(list(pl, pl1, pl2, pl3, pl4, pl5,
    pl6, pl7, pl8, pl9, pl10, pl11),
    function(p) {
  # Extract the legend from the plot
  g_legend <- function(p) {
    tmp <- ggplotGrob(p)
    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
    legend <- tmp$grobs[[leg]]
    return(legend)
  }
  g_legend(p)
})

# Select the legend you want to display at the bottom
lgd <- legends[[3]]

# Plots without Legend
main_plot <- grid.arrange(pl + theme(legend.position = "none"),
    pl1 + theme(legend.position = "none"),
    pl2 + theme(legend.position = "none"),
    pl3 + theme(legend.position = "none"),
    pl4 + theme(legend.position = "none"),
    pl5 + theme(legend.position = "none"),
    pl6 + theme(legend.position = "none"),
    pl7 + theme(legend.position = "none"),
    pl8 + theme(legend.position = "none"),
    pl9+ theme(legend.position = "none"),
    pl10 + theme(legend.position = "none"),
    pl11 + theme(legend.position = "none"), nrow = 4,
    left = textGrob("
    Average Weekly Wage
    ", rot = 90,
        gp = gpar(fontface = "bold", fontsize = 15)))

All Plots with Legend

grid.arrange(main_plot, lgd, nrow = 2, heights = c(13, 1))