A Multiverse Analysis of Interaction Effects

Arbitrary vs. Non-Arbitrary Decisions

Processing raw data involves many steps, such as cleaning, preparing, and creating analysis variables (among many others). Such decisions eventually produce an analysis-ready dataset, but making (even slightly) different decisions will produce a different dataset. And so emerges the multiverse: all possible datasets that could be generated from different data processing decisions (see also the garden of forking paths).

An important distinction in multiverse analysis is between arbitrary and non-arbitrary decisions. Arbitrary decisions are those that:

1. are not expected to influence results (or their effect is unknown) and,
2. equally defensible/have reasonable justifications

For example, whether to remove outliers at +/- 3 or 2.5 standard deviations could be arbitrary. Neither cutoff is theoretically informed and both are equally defensible.

In contrast, non-arbitrary choices are expected to affect results, either for theoretical or empirical reasons. For example, theory might suggest controlling for a confound is crucial for interpreting results. Empirically, controlling for a confound may aid interpretation because such a confound is known to correlate with the outcome of interest and should be accounted for in our analysis.

Multiverse analysis should be ran over all possible arbitrary decisions. This is important because the multiverse can quickly become unnecessarily large and difficult to interpret. Ideally, the multiverse will comprise all datasets that, by themselves, would make sense to analyze and interpret individually.

Setup

To conduct my multiverse analysis, I will be using a number of packages from the tidyverse. Below are the specific packages I used for my setup:

# Load relevant libraries
library(tidyverse)
library(lmerTest)
library(broom.mixed)
library(tidymodels)
library(ggeffects)
library(interactions)
library(effectsize)
library(showtext)
library(furrr)
# Custom font for ggplot
font_add_google("Lato")
showtext_auto()
# Set the ggplot2 theme
theme_set(
  theme_light() +
    theme(
      text = element_text(family="Lato"),
      panel.grid = element_blank()
    )
)

Simulated Data

The raw data are randomly generated below. I modeled the variables based on an actual project I am working on. The research question is whether or not early life unpredictability¹ (iv_unp) is associated with performance on a cognitive test with two versions, a standard (dv_std) and ecologically relevant (dv_eco) version². We hypothesize that unpredictability negatively impacts performance on the standard version but enhances performance on the ecological version. To that end, I’ve made performance depend on iv_unp with some error.

Exclusion variables are also included in the data.

• ex_sample: The data (not these but those that I will be actually analyzing in a real dataset) come from two samples. The first sample is smaller and considerably more difficult to run our protocol. The second sample is larger and was much easier to control how the protocol was executed. Thus our results could be different depending on whether or not we include the smaller, noisier sample.
• ex_csd: socially desirable responding, which may have undue influence on our predictor. Children who experience adversity, and might want to hide it, may also score high on a measure of socially desirability.
• ex_sped: Some in our sample receive special education. It is possible that some students may have cognitive impairments and shouldn’t be analyzed with rest or the sample.
• ex_att: We included an a few trap questions in our assessment to detect individuals who may not be taking the study seriously.
• ex_imp: Our research assistants noted whether or not participants had trouble reading or seemed impaired somehow.

These variables will be used to construct a multiverse of datasets.

# Set the seed for reproducibility
set.seed(12345)
# detect cores for parrallel processing
cores <- parallel::detectCores()
# Create some data
sim_data <- 
  tibble(
    id          = 1:500,
    intercept   = rnorm(500, mean = 0, sd = 3),
    iv_unp      = rnorm(500, mean = 0, sd = 1),
    ex_sample   = rbinom(500, size = 1, prob = .2),
    ex_csd      = rnorm(500, mean = 0, sd = 1) + iv_unp * (.15) + rnorm(500, mean = 0, 1),
    ex_sped     = rbinom(500, size = 1, prob = .1),
    ex_att      = rbinom(500, size = 1, prob = .1),
    ex_imp      = rbinom(500, size = 1, prob = .05)
  ) %>% 
  mutate(
    intercept   = ifelse(ex_sped == 1, intercept - rnorm(1,.5,.25), intercept),
    intercept   = ifelse(ex_imp == 1, intercept - rnorm(1,.5,.25), intercept),
    intercept   = ifelse(ex_att == 1, intercept - rnorm(1,.25,.25), intercept),
    dv_std      = intercept + iv_unp * -(.2) + rnorm(500, mean = 0, 2),
    dv_eco      = intercept + iv_unp *  (.1) + rnorm(500, mean = 0, 2)
  )

And here are the first 10 rows of the simulated data:

Check Sample Sizes

One issue that could arise in multiverse analysis is that the sample sizes of some individual datasets may be very small, at least compared to the full dataset. If this occurs, it may not make sense to analyze the datasets with small sample sizes. To check this issue, I can plot the distribution of sample sizes as a function of the decisions in our decision grid:

# create a sample size data.frame from our data list
mulitverse_sample_size <- 
  map_df(multiverse_data, function(x){
    # Grab the decisions flatten them out
    decisions <- x$decisions %>% 
      gather(key,value,-decision,-variable) %>% 
      unite(grid_vars, variable, key) %>% 
      spread(grid_vars, value)
    
    # Add the decisions to the model table
    bind_cols(n = x$n, decisions)
  })
# histogram of sample sizes in our multiverse
mulitverse_sample_size %>% 
  ggplot(aes(x = n)) +
  geom_histogram(color = "black", bins = 15) +
  labs(title = "Histogram of Multiverse Sample Sizes", x = "\nSample Size", y = "Number of Datasets\n")

In my case, the minimum sample size was 304 and the maximum was 500. We can also look at which data processing decisions affected sample size the most:

# boxplots showing how sample sizes look across different decisions
mulitverse_sample_size %>% 
  select(-ends_with("exp")) %>% 
  gather(variable, exclusion, -decision, -n) %>%
  mutate(variable = str_replace(variable, "_(name)$|_(variable_group)$", ".\\1\\2")) %>% 
  separate(variable, c("var","group"), sep = "\\.") %>% 
  spread(group, exclusion) %>% 
  arrange(decision, variable_group) %>% 
  ggplot() +
  geom_boxplot(aes(x = name, y = n), notch = T, width = .5)  +
  scale_x_discrete("Exclusion Criteria") +
  scale_y_continuous("Sample Size (max = 500)\n") +
  ggtitle("Sample Sizes in the Multiverse") +
  facet_wrap(~variable_group, scales = "free", ncol = 2)

Multiverse Effect Sizes

With the multiverse of results, I can extract the coefficients (i.e. regression slopes) for plotting and interpretation. Here, I just loop through my multiverse of results to get this information (i.e., fixed main effects and interaction term).

# Get the model estimates out of the giant list
multiverse_estimates <- 
  multiverse_models %>% 
  map_df(function(x){
    # Grab the decisions flatten them out
    decisions <- 
      x$decisions %>% 
      gather(key,value,-decision,-variable) %>% 
      unite(grid_vars, variable, key) %>% 
      spread(grid_vars, value)
    
    # Add the decisions to the model table
    bind_cols(x$mod_tidy, n = x$n, decisions)
  })
# Only look at our main effects and interaction terms
multiverse_plot_data <- 
  multiverse_estimates %>% 
  filter(str_detect(mod_term,"^iv_unp_z$|^iv_unp_resid_z$|task")) %>% 
  mutate(
    term = str_remove(mod_term,"_resid"),
    term = factor(term, 
                  levels = c("iv_unp_z","task_type","task_type:iv_unp_z"), 
                  labels = c("Main Effect Unp","Main Effect Task","Interaction")),
    p_value = case_when(mod_p.value <= .05 ~ "p <= .05",
                        mod_p.value >  .05 ~ "p > .05"),
    p_value = factor(p_value, levels = c("p > .05","p <= .05"))
  ) %>% 
  group_by(term) %>% 
  mutate(multiverse_set = as.numeric(fct_reorder(decision, mod_estimate)))

Effect Size Curve

multiverse_effect_curve <- 
  multiverse_plot_data %>% 
  ggplot(aes(x = as.numeric(multiverse_set), y = mod_estimate, color = p_value)) +
  geom_point(size = .5) +
  geom_errorbar(
    aes(x = as.numeric(multiverse_set), ymin = mod_estimate - mod_std.error, ymax = mod_estimate + mod_std.error, color = p_value),
    size = .25,
    width = 0,
    alpha = .25, 
    inherit.aes = F, 
    show.legend = T
  ) +
  scale_y_continuous("Estimate") +
  scale_color_manual(values = c("gray","#f0506e")) +
  scale_fill_manual(values = c("gray","#f0506e")) +
  facet_wrap(~term) +
  theme(
    axis.text.x = element_blank(), 
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank()
  )
multiverse_effect_curve

Specification Grid

multiverse_spec_grid <- 
  multiverse_plot_data %>% 
  select(-matches("exp$")) %>% 
  gather(variable, spec, matches("name|group$"), -starts_with("mod"), -decision, -multiverse_set) %>%
  mutate(variable = str_replace(variable, "_(name$)|_(variable_group$)",".\\1\\2")) %>% 
  separate(variable, c("variable","group"), sep = "\\.") %>% 
  spread(group,spec) %>% 
  ggplot(aes(x = multiverse_set, y = name, color = p_value)) +
  geom_point(size = .75, shape = 124, show.legend = F) +
  scale_y_discrete("Specification\n") +
  scale_x_continuous("\nMultiverse Dataset") +
  scale_color_manual(values = c("gray","#f0506e")) +
  facet_grid(variable_group~term, scales = "free_y") +
  theme(strip.background.x = element_blank(), strip.text.x = element_blank(), strip.placement = "outside")
multiverse_spec_grid

Creating the Full Multiverse Plot

cowplot::plot_grid(
  multiverse_effect_curve,
  multiverse_spec_grid,
  ncol = 1, 
  align = 'v', 
  axis = "lr",
  rel_heights = c(.25,.75)
)

The Histogram of p-values

multiverse_plot_data %>% 
  ggplot(aes(x = mod_p.value)) +
  geom_histogram(color = "black") +
  geom_vline(aes(xintercept = .05), color = "red", linetype = "dotted") +
  scale_x_continuous("\nMultiverse p-value") +
  scale_y_continuous("Count") +
  facet_wrap(~term, scales = "free_y")

Simple Slopes Effect Curves

These curves are, in principle, the same as the effect size curves above. The difference is that the effect size of interest is a simple slope instead of a regression coefficient from the mixed model. There are four simple effects of interest:

• The effect of task type when unpredictability is low
• The effect of task type when unpredictability is high
• The effect of unpredictability on the standard version
• The effect of unpredictability on the ecological version

The code below grabs this information from our multiverse analysis so I can visualize these effects. Note that I could also make a specification grid for each of these plots but I opted not to for brevity.

# Get the simple slopes tests from our giant list
multiverse_ss <- map_df(multiverse_models, function(x){
  decisions <- x$decisions %>% 
    gather(key,value,-decision,-variable) %>% 
    unite(grid_vars, variable, key) %>% 
    spread(grid_vars, value)
  
  simple_slopes <- 
    bind_rows(
      x$mod_sim_slopes$unp$slopes %>% 
        rename_all(~c("value", "slope", "se","lower","upper","t.value","p.value")) %>% 
        mutate(group = "unp"),
      x$mod_sim_slopes$task$slopes %>% 
        rename_all(~c("value", "slope", "se","lower","upper","t.value","p.value")) %>% 
        mutate(group = "task"),
    )
  
  bind_cols(simple_slopes, decisions)
}) %>% 
  mutate(simple_slope = case_when(group=="unp" & value ==  1 ~ "High Unp",
                                  group=="unp" & value == -1 ~ "Low Unp",
                                  group=="task" & value ==  1 ~ "Ecological",
                                  group=="task" & value == -1 ~ "Standard"),
         simple_slope = factor(simple_slope, levels = c("Low Unp","High Unp", "Standard", "Ecological")),
         p_value = case_when(p.value <= .05 ~ "p <= .05",
                             p.value >  .05 ~ "p > .05"),
         p_value = factor(p_value, levels = c("p > .05","p <= .05")))

multiverse_ss %>% 
  filter(group == "unp") %>% 
  group_by(group, value) %>% 
  mutate(multiverse_set = as.numeric(fct_reorder(decision, slope))) %>% 
  ggplot(aes(x = multiverse_set, y = slope, color = p_value)) +
  geom_point() +
  scale_x_continuous("\nMultiverse") +
  scale_y_continuous("Estimate\n") +
  scale_color_manual(values = c("gray","#f0506e")) +
  scale_fill_manual(values = c("gray","#f0506e")) +
  facet_wrap(~simple_slope)

multiverse_ss %>% 
  filter(group == "task") %>% 
  group_by(group, value) %>% 
  mutate(multiverse_set = as.numeric(fct_reorder(decision, slope))) %>% 
  ggplot(aes(x = multiverse_set, y = slope, color = p_value)) +
  geom_point() +
  scale_x_continuous("\nMultiverse") +
  scale_y_continuous("Estimate\n") +
  scale_color_manual(values = c("gray","#f0506e")) +
  scale_fill_manual(values = c("gray","#f0506e")) +
  facet_wrap(~simple_slope)

Interaction plots

Here, I simply plot interactions based on high versus low unpredictability and task-type. More specifically, the predicted level of performance when unpredictability is high versus low and on the standard versus ecological version.

Each line represents an interaction from a single dataset in the multiverse. I’ve broken the plots out into panels so the overall pattern is more visible. I’ve also made one plot with unpredictability on the x-axis and separate panels for task type. I made another plot with these flipped.

# Get the predicted points out of our giant list
multiverse_effects <- map_df(multiverse_models, function(x){
  decisions <- x$decisions %>% 
    gather(key,value,-decision,-variable) %>% 
    unite(grid_vars, variable, key) %>% 
    spread(grid_vars, value)
  
  bind_cols(x$mod_effects, n = x$n, decisions)
}) %>% 
  left_join(
    multiverse_plot_data %>% filter(str_detect(mod_term,"\\:")) %>% select(term, p_value, decision),
    by = "decision"
  )

# Plot the interactions
multiverse_effects %>% 
  mutate(group = factor(group, levels = c(-1,1), labels = c("Standard Version","Ecological Version"))) %>% 
  ggplot(aes(x = x, y = predicted, group = decision, color = p_value)) +
  geom_line(alpha = .25,size = .25, show.legend = F) +
  scale_x_continuous("\nEnvironmental Unpredictability") +
  scale_y_continuous("Performance\n") +
  scale_color_manual(values = c("gray","#f0506e")) +
  facet_wrap(~group) 

# Plot the interactions
multiverse_effects %>% 
  mutate(
    group_jitter = jitter(as.numeric(group), amount = .1),
    group        = factor(group, levels = c(-1,1), labels = c("Standard Version","Ecological Version")),
    x            = factor(x, levels = c(-1,1), labels = c("Low Unpredictability","High Unpredictability"))
  ) %>% 
  ggplot(aes(x = group_jitter, y = predicted, group = decision, color = p_value)) +
  geom_line(alpha = .25,size = .25, show.legend = F) +
  geom_point(alpha = .25,size = .25, show.legend = F) +
  scale_x_continuous("\nTask Type", breaks = c(1,2), labels = c("Standard","Ecological")) +
  scale_y_continuous("Performance\n") +
  scale_color_manual(values = c("gray","#f0506e")) +
  facet_wrap(~x)