reycn · January 22, 2024 18:11
diff --git a/ggplot-2-did-trends-interaction b/ggplot-2-did-trends-interaction
 # 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)
	# 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)