# How to visualize interaction effects
# R code from https://philippmasur.de/blog/2018-11-26-visualizing-interaction-effects/

# Loading packages
library(tidyverse)   # includes 'ggplot2' which is probably the best package for data viz
library(ggExtra)     # extends ggplot2 

# Simulating the data
set.seed(100)
m1 <- sample(rep(c(0, 1), 250), 500) # Factor A
m2 <- sample(rep(c(0, 1), 250), 500) # Factor B
y <- 8 + -2*m1 + 2*m2 + 3*(m1*m2) + rnorm(500, 0, 2)
(data <- data.frame(m1 = factor(m1, 
                                label = c("A1", "A2")), 
                    m2 = factor(m2,
                                label = c("B1", "B2")),
                    y) %>%
  as_tibble)

ggplot(data, 
       aes(x = m1,
           fill = m1,  
           y = y)) +
  stat_summary(fun = mean,
               geom = "bar") +
  stat_summary(fun.data = mean_cl_normal, 
               geom="errorbar", 
               width = 0.25) +
  ylim(0, 17) +
  facet_wrap(~m2) +
  labs(x = "Experimental Manipulations",
       y = "Dependent Variable (M)") +
  theme_minimal() +
  scale_fill_brewer(palette = "Blues") +
  theme(legend.position="none")

# Summarizing data
data2 <- data %>% 
  group_by(m1, m2) %>% 
  summarise(y_mean = mean(y),
            y_se = psych::describe(y)$se)

# Creating the plot
data2 %>% 
  ggplot(aes(x = m1, 
             y = y_mean, 
             color = m2)) +
  geom_line(aes(group = m2)) +
  geom_point() +
  geom_errorbar(aes(ymin = y_mean-1.96*y_se, 
                      ymax = y_mean+1.96*y_se),
                width = .1) +
  ylim(0, 17) +
  labs(x = "Factor A",
       color  = "Factor B",
       y = "Dependent Variable") +
  theme_minimal() +
  scale_color_brewer(palette = "Blues")

# Simulating the data
x <- rnorm(500, 2, 3)
m <- rep(c(0,1), 250)
y <- 5 + -4*x + 2*m + 4*(x*m) + rnorm(500, 0, 5)
(data <- data.frame(x, m = factor(m), y) %>%
  as_tibble)

# Setting up the building blocks
basic_plot <- ggplot(data,
       aes(x = x,
           y = y,
           color = m)) +
  theme_bw() +
  labs(x = "Independent variable",
       y = "Dependent variable",
       color = "Moderator")

# Colored scatterplot
basic_plot +
  geom_point()

# Colored scatterplot and regression lines
basic_plot +
  geom_point(alpha = .3, 
             size = .9) +
  geom_smooth(method = "lm")

# ggsave("test2.png",
#        height = 4,
#        width = 5)

# Only regression lines
basic_plot + 
  geom_smooth(method = "lm",
              se = F)

# Simulating the data
n <- 500
x <- runif(n)
m <- rnorm(n, 0.5, 1)
y <- 2*x + -2*m + -6*(x*m) + rnorm(500, 0, 5.88)
(data <- data.frame(x, m, y) %>%
  as_tibble)

data$m_groups <- cut(data$m, breaks = 3) %>% 
  factor(., labels = c("small", "average", "large"))
table(data$m_groups)

ggplot(data,
       aes(x = x,
           y = y,
           color = m_groups)) +
  geom_point(size = .9,
             alpha = .3) +
  geom_smooth(method = "lm") +
  theme_bw() +
  scale_color_brewer(type = "qual", 
                     palette = 3) +
  labs(x = "Independent variable",
       y = "Dependent variable",
       color = "Moderator")

model <- lm(y ~ x + m + x:m, data)
summary(model)

# Function that estimates the conditinal effects for certain values of the moderator
conditional.effects <-  function(model, x, m, quantiles = 10){
  interact = paste0(x,':',m)
  beta.hat = coef(model) 
  covs = vcov(model)
  z0  = quantile(model$model[,m], seq(0 , 1, 1/quantiles))
  dy.dx = beta.hat[x] + beta.hat[interact]*z0
  se.dy.dx = sqrt(covs[x, x] + z0^2*covs[interact, interact] + 2*z0*covs[x, interact])
  upr = dy.dx+1.96*se.dy.dx
  lwr = dy.dx-1.96*se.dy.dx
  data.frame(m = z0, b = dy.dx, lwr, upr)
}

(con_effects <- conditional.effects(model, 
                                   x = "x", 
                                   m = "m", 
                                   quantiles = 500) %>%
  as_tibble)

(mod_plot <- ggplot(con_effects, 
       aes(x = m, 
           y = b, 
           ymin = lwr,
           ymax = upr)) +
  geom_smooth(stat = "identity",
              color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  theme_bw() +
  labs(x = "Moderator",
       y = "Conditional effect of x on y"))

# ggExtra now needs a point layer, so we add invisible points (alpha = 0)
ggMarginal(mod_plot + geom_point(alpha = 0), # function of the package `ggExtra`
           margins = c("x"),
           type = "histogram",
           color = "white",
           fill = "darkgrey",
           size = 3)