## ----set-options, echo = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  fig.width = 7,
  fig.height = 3.5,
  message = FALSE,
  warning = FALSE,
  package.startup.message = FALSE
)
options(width = 800)
arrow_color <- "#FF00cc"
my_predictions <- NULL

pkgs <- c("ggplot2", "marginaleffects", "see", "parameters")

options(modelbased_join_dots = FALSE)
options(modelbased_select = "minimal")

if (!all(insight::check_if_installed(pkgs, quietly = TRUE))) {
  knitr::opts_chunk$set(eval = FALSE)
}

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(modelbased)
library(parameters)
library(ggplot2)

set.seed(123)
n <- 200
d <- data.frame(
  outcome = rnorm(n),
  grp = as.factor(sample(c("treatment", "control"), n, TRUE)),
  episode = as.factor(sample(1:3, n, TRUE)),
  sex = as.factor(sample(c("female", "male"), n, TRUE, prob = c(0.4, 0.6)))
)
model1 <- lm(outcome ~ grp + episode, data = d)
model_parameters(model1)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
my_predictions <- estimate_means(model1, "episode")
my_predictions

plot(my_predictions)

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
p <- plot(my_predictions, pointrange = list(position = position_dodge(width = 0.4)))
line_data <- as.data.frame(my_predictions)[1:2, ]
p + annotate(
  "segment",
  x = as.numeric(line_data$episode[1]) + 0.06,
  xend = as.numeric(line_data$episode[2]) - 0.06,
  y = line_data$Mean[1],
  yend = line_data$Mean[2],
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
  ggtitle("Within \"episode\", do levels 1 and 2 differ?")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# argument `comparison` defaults to "pairwise"
estimate_contrasts(model1, "episode")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(model1, contrast = "episode=c(1,2)")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
model2 <- lm(outcome ~ grp * episode, data = d)
model_parameters(model2)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
my_predictions <- estimate_means(model2, by = c("episode", "grp"))
my_predictions

plot(my_predictions)

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
p <- plot(my_predictions, pointrange = list(position = position_dodge(width = 0.4)))
line_data <- as.data.frame(my_predictions)[3:4, 1:3]
line_data$group_col <- "control"
p + annotate(
  "segment",
  x = as.numeric(line_data$episode[1]) - 0.06,
  xend = as.numeric(line_data$episode[2]) + 0.06,
  y = line_data$Mean[1],
  yend = line_data$Mean[2],
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
  ggtitle("Within level 2 of \"episode\", do treatment and control group differ?")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# we want "episode = 2-2" and "grp = control-treatment"
estimate_contrasts(model2, contrast = c("episode", "grp"))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_means(model2, by = c("episode", "grp"))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# compute specific contrast directly
estimate_contrasts(model2, contrast = c("episode", "grp"), comparison = "b2 = b5")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# return pairwise comparisons for specific values, in
# this case for episode = 2 in both groups
estimate_contrasts(model2, contrast = c("episode=2", "grp"))

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
p <- plot(my_predictions)
line_data <- as.data.frame(my_predictions)[c(2, 5), 1:3]
line_data$group_col <- "treatment"
p + annotate(
  "segment",
  x = as.numeric(line_data$episode[1]) + 0.06,
  xend = as.numeric(line_data$episode[2]) - 0.06,
  y = line_data$Mean[1],
  yend = line_data$Mean[2],
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
  ggtitle("Do different episode levels differ between groups?")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(model2, contrast = c("episode", "grp"), comparison = "b4 = b3")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_means(model2, by = c("episode=c(1,3)", "grp"))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(
  model2,
  contrast = c("episode = c(1, 3)", "grp"),
  comparison = "b3 = b2"
)

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
my_predictions <- estimate_means(model2, c("grp", "episode"))
p <- plot(my_predictions, pointrange = list(position = position_dodge(width = 0.4)))
line_data <- as.data.frame(my_predictions)[, 1:3, ]
line_data$group_col <- "1"
p + annotate(
  "segment",
  x = as.numeric(line_data$grp[1]) - 0.05,
  xend = as.numeric(line_data$grp[1]) - 0.05,
  y = line_data$Mean[1],
  yend = line_data$Mean[2],
  color = "orange",
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40, type = "closed")
) + annotate(
  "segment",
  x = as.numeric(line_data$grp[4]) - 0.05,
  xend = as.numeric(line_data$grp[4]) - 0.05,
  y = line_data$Mean[4],
  yend = line_data$Mean[5],
  color = "orange",
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40, type = "closed")
) + annotate(
  "segment",
  x = as.numeric(line_data$grp[1]) - 0.05,
  xend = as.numeric(line_data$grp[4]) - 0.05,
  y = (line_data$Mean[1] + line_data$Mean[2]) / 2,
  yend = (line_data$Mean[4] + line_data$Mean[5]) / 2,
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40, type = "closed")
) +
  ggtitle("Differnce-in-differences")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_means(model2, c("episode", "grp"))

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(
  model2,
  c("episode", "grp"),
  comparison = "(b1 - b2) = (b4 - b5)"
)