Skip to content

Instantly share code, notes, and snippets.

@reycn

reycn/ggplot-2-did-trends-interaction

Created January 22, 2024 18:11
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save reycn/33e925eac77f8720dbaec24a97864f21 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save reycn/33e925eac77f8720dbaec24a97864f21 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. reycn created this gist Jan 22, 2024.
    98 changes: 98 additions & 0 deletions ggplot-2-did-trends-interaction
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,98 @@
    # Interaction
    library(ggplot2)
    library(dplyr)
    tips <- read.csv("https://sebastiansauer.github.io/data/tips.csv")
    tips %>%
    group_by(sex, smoker) %>%
    summarise(tip_groups = mean(tip)) -> tips2
    nature_colors <- c("#9B241B", "#014E8D", "#9B7409") # nature brand

    tips2 %>%
    ggplot() +
    aes(x = sex, y = tip_groups, color = smoker) +
    geom_line(aes(group = smoker)) +
    geom_point() +
    geom_hline(
    yintercept = mean(tips$tip),
    linetype = 'dotted',
    col = '#697481'
    ) +
    theme_minimal() +
    theme(
    panel.background = element_rect(fill = "white"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
    ) +
    scale_color_manual(values = nature_colors) +
    labs(
    title = "Moderated Treatment Effects",
    x = "X",
    y = "Y",
    color = "Moderator"
    )

    # Trends
    library(dplyr)
    set.seed(12345)
    T = 100 # no of periods
    N = 40 # no of subjects

    dat = expand.grid(t = 1:T, i = 1:N)

    # Simulate a common AR(1) time trend
    time.trend = as.numeric(arima.sim(n = T, list(
    ar = c(0.4, 0.5), ma = c(0.6, 0.5)
    ))) * 3 + 0.7 * (1:T)

    dat = mutate(
    dat,
    group = ifelse(i > N / 2, "treat", "control"),
    treat = 1L * (group == "treat"),
    exp = 1L * (t > T / 2),
    treat_exp = exp * treat,
    mu.t = time.trend[t],
    eps = rnorm(n()),
    y = mu.t + treat * 40 + treat_exp * 50 + eps
    )
    sample_n(dat, 5)

    show.plot = function(dat,
    label = "",
    show.means = TRUE) {
    library(ggplot2)
    gdat = dat %>%
    group_by(group, t, exp, treat) %>%
    summarize(y = mean(y))

    nature_colors <- c("#9B241B", "#014E8D", "#9B7409") # nature brand
    x_labels <- c(-2, -1, "Treatment", 1, 2)
    gg = ggplot(gdat, aes(y = y, x = t, color = group)) +
    geom_line() +
    stat_smooth(
    linetype = "dashed",
    size = 0.5,
    level = 0.95,
    alpha = 0.1
    ) +
    geom_vline(xintercept = T / 2) +
    theme_bw() +
    annotate("text",
    x = T / 4,
    y = 0.9 * max(gdat$y),
    label = label) +
    scale_color_manual(values = nature_colors) +
    scale_x_continuous(breaks = c(0, 25, 50, 75, 100), labels = x_labels) +
    labs(
    title = "Trends",
    x = "Seasons",
    y = "Y",
    color = "Group"
    )
    # +
    # annotate("label",
    # x = T / 2,
    # y = 150,
    # label = "Treatment")
    gg
    }
    show.plot(dat)