Multiverse analysis Steegen et al. (2016)

• Even a simple data analysis is characterized by several (often arbitrary) choices
• The impact of these choices is not always clear, known or easy to predict
• Researchers (with or without awareness) usually report a single set of choices

Some examples

• using a predictor as continous or choosing some thresholds and transforming to categorical
• including or excluding an observation (e.g., outlier)
• including or excluding a covariate (e.g., controlling for age or not)
• using a linear model or an ordinal regression for discrete ordered data
• …

The garden of forking paths

A tree of possibilities…

The garden of forking paths

… Where only some of them produce a certain result.

Only plausible vs all scenarios

Not all scenarios are equally plausible or reasonable. Thus the multiverse is not the entire set but a non-random sample of plausibile choices.

x <- seq(0, 1, 0.0001)
y <- dexp(x, 2)
y <- (y-min(y))/(max(y)-min(y))

plot(x*100, y, type = "l",
     xlab = "All possible scenarios",
     ylab = "Plausibility",
     lwd = 2)

Meta-analysis many choices

Despite useful and very powerful, meta-analysis is characterized by several (arbitrary) choices. For example:

• Should the study \(x\) be excluded for theoretical or statistical (e.g., outliers) reasons?
• Should we use an equal or random-effects model?
• Which value should take the pre-post missing correlation?
• …

An example: Pre-post designs

Let’s imagine that we have \(k\) studies about the efficacy of a certain treatment. They collected a sample of participants measuring a certain variable \(y\) before and after the treatment.

With this design we can summarise the effect computing the difference between the average of the two time points.

An effect size (e.g., Cohen’s \(d\)) can be calculated as:

\[ d = \frac{\overline y_{post} - \overline y_{pre}}{\sigma_p} \]

And \(\sigma_p\) is the pooled pre-post standard deviation.

An example: Pre-post Cohen’s \(d\)

With a pre-post Cohen’s \(d\) we need the pre-post correlation \(\rho\) to calculate the sampling variance:

\[ \sigma^2_{\epsilon_{pp}} = \frac{2(1 - \rho)}{n} + \frac{d^2}{2n} \]

\(\rho\) is usually non reported and need to be chosen from previous literature or a plausible guess.

Pre-post Cohen’s \(d\)

v_dm <- function(r, d, n){
    2*(1-r)/n + d^2/(2*n)
}

vdd <- data.frame(x = seq(-1, 1, 0.001))
vdd$y <- v_dm(vdd$x, 0, 40)

ggplot(vdd, aes(x = x, y = y)) +
    geom_line() +
    xlab(expression(rho)) +
    ylab(expression(sigma^2)) +
    theme_minimal(25)

A simulated example

The data structure: multiple measures of the same outcome within each paper

Simulation approach

Let’s make an example for a paper with \(j = 3\) measures of the outcome:

\[ \begin{bmatrix} y_{11} \\ y_{12} \\ y_{13} \end{bmatrix} = \begin{bmatrix} \mu_{\theta_1} \\ \mu_{\theta_2} \\ \mu_{\theta_3} \end{bmatrix} + \begin{bmatrix} \delta_{\theta_1} \\ \delta_{\theta_2} \\ \delta_{\theta_3} \end{bmatrix} + \begin{bmatrix} \epsilon_{\theta_{11}} \\ \epsilon_{\theta_{12}} \\ \epsilon_{\theta_{13}} \end{bmatrix} \]

\[ \boldsymbol{\delta} \sim \text{MVN}\left( \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \ \begin{matrix} \tau^2_1 & & \\ \rho_{21}\tau_2\tau_1 & \tau^2_1 & \\ \rho_{31}\tau_2\tau_1 & \rho_{32}\tau_2\tau_1 & \tau^2_1 \\ \end{matrix} \right) \]

\[ \boldsymbol{\epsilon} \sim \text{MVN}\left( \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} \ \begin{matrix} \sigma^2_{\epsilon_1} & & \\ \rho_{21}\sigma^2_{\epsilon_2}\sigma^2_{\epsilon_1} & \sigma^2_{\epsilon_2} & \\ \rho_{31}\sigma^2_{\epsilon_3}\sigma^2_{\epsilon_1} & \rho_{32}\sigma^2_{\epsilon_3}\sigma^2_{\epsilon_2} & \sigma^2_{\epsilon_3} \\ \end{matrix} \right) \]

Simulation approach

We simulated individual participant data, thus:

1. Sampling the true values \(\theta_{ij}\) for each study \(i\) and outcome \(j\) from the multivariate distribution
2. Generating \(n_i\) pre and post data with correlation \(\rho\)
3. Calculating the effect size (imputing the pre-post correlation)
4. Aggregating multiple outcomes within the same paper (imputing the correlation)
5. Fitting the meta-analyis model
6. Calculating the scores
7. Repeating 3-4 for each scenario
8. Using PIMA

Simulation approach

We simulated a relatively simple but plausible multiverse with:

• 4 pre-post correlations
• 4 correlations between multiple measures of the same outcome
• 2 meta-analysis models (fixed and random-effects)

For a total of 32 multiverse scenarios.

Results

p_beta <- ggplot(multi, aes(y = b)) +
    xlim(c(-0.5,0.5)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("$\\mu_{\\,\\, \\theta}$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank())

p_pval <- ggplot(multi, aes(y = -log10(p.value))) +
    xlim(c(-0.5,0.5)) +
    geom_hline(yintercept = -log10(0.05)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("Raw $-log_{10}(p)$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank())

rdata <- cor(zi)
rdata <- data.frame(r = rdata[upper.tri(rdata, diag = FALSE)])

p_cor <- ggplot(rdata, aes(y = r)) +
    xlim(c(-0.5,0.5)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("$\\rho$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank()) +
    ylim(c(-1,1))

cowplot::plot_grid(p_beta, p_cor, p_pval, nrow = 1)

An example from Plessen et al. (2023)

An example from Plessen et al. (2023)

• Over the last four decades, more than 80 meta-analyses have examined the efficacy of psychotherapies for depression
• More than 700 randomised controlled trials (RCTs)
• Not all these studies goes in the same direction

An example from Plessen et al. (2023)

Discrepancies in results could be due to:

Which factors (which data to meta-analyze)

• inclusion/esclusion of a subset of studies (e.g., low quality studies)
• type of control group or control therapy
• …

How factors (how to meta-analyze)

• type of model (e.g., equal vs random)
• model complexity (two-level, three level, robust, etc.)
• correcting for publication bias

Describing the multiverse

Hall et al. (2022) and Liu et al. (2021) proposed several tools to describe the results of a multiverse analysis.

The increase in complexity NEED to be managed using appropriate tools to summarise and visualize the results.

Estimated effect of interest

specif_res |> 
  ggplot(aes(x = b)) +
  geom_histogram(bins = 50,
                 fill = "dodgerblue",
                 col = "black") +
  ylab("Frequency") +
  xlab("Estimated Effect") +
  geom_vline(xintercept = 0, col = "red") +
  geom_vline(xintercept = 0.25, col = "green")

Estimated effect as a function of scenarios¹

above <- specif_res |> 
  arrange(b) |> 
  mutate(id = 1:n()) |> 
  ggplot(aes(x = id, y = b)) +
  geom_segment(aes(xend = id, y = ci.lb, yend = ci.ub)) +
  geom_hline(yintercept = 0, col = "red") +
  geom_hline(yintercept = 0.25, col = "green") +
  theme(axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank())
  
below <- specif_res |> 
  arrange(b) |> 
  mutate(id = 1:n()) |> 
  pivot_longer(c(
    target_group,
    format,
    diagnosis,
    risk_of_bias,
    type,
    method
  )) |> 
  ggplot(aes(x = id, y = value)) +
  geom_tile() +
  theme(axis.title.y = element_blank(),
        axis.text.x = element_blank()) +
  xlab("Specification")

plot_grid(above, below, nrow = 2, align = "hv")

Fixing a specific parameter

We can the effect of a set of scenarios having a certain parameter:

specif_res |> 
  arrange(b) |> 
  mutate(id = 1:n()) |> 
  ggplot(aes(x = b, fill = diagnosis)) +
  geom_density(alpha = 0.5)

Specification Curve Simonsohn, Simmons, and Nelson (2020)

The idea of the specification curve is both descriptive (see previous plots) and inferential.

Inferentially, the approach requires to generate a new dataset under the null hypothesis and compare with the observed value.

Specification Curve in meta-analysis Plessen et al. (2023)

Plessen et al. (2023) describe how to implement the method for meta-analysis.

1. For each dataset/model generate \(k_s\) data (where \(k\) is the number of effect in the specification \(s\)) sampling from \(\mathcal{N}(0,\sigma^2_{\epsilon_i} + \tau^2)\).
2. Fit the model (FE or RE)
3. Get the specification curve
4. Repeat 1-3 for a large number of times (e.g., 1000)
5. Compute the 2.5% and 95.75% quantiles
6. Check if the observed specification is outside the 95% CI

Specification Curve in meta-analysis Plessen et al. (2023)

specif$SD |> 
  filter(spec %in% sample(1:1000, 1)) |> 
  ggplot(aes(x = id, y = b)) +
  geom_line(aes(group = spec))

Specification Curve in meta-analysis Plessen et al. (2023)

specif$SD |> 
  filter(spec %in% sample(1:1000, 1)) |> 
  ggplot(aes(x = id, y = b)) +
  geom_line(aes(group = spec))

Specification Curve in meta-analysis Plessen et al. (2023)

specif$SD |> 
  filter(spec %in% sample(1:1000, 1)) |> 
  ggplot(aes(x = id, y = b)) +
  geom_line(aes(group = spec))

Specification Curve in meta-analysis Plessen et al. (2023)

specif_res <- specif_res |> 
  arrange(b) |> 
  mutate(id = 1:n())

ggplot(data = specif$SD,
       aes(x = id, y = b)) +
  geom_line(aes(group = spec)) +
  geom_line(data = specif_res,
            aes(x = id, y = b),
            col = "firebrick",
            lwd = 2)

PIMA

P-selection Inference in Multiverse Meta-analysis (PIMA) is an inferential approach to multiverse analysis is an inferential framework for:

• obtaining an overall p-value across the multiverse
• obtaining corrected p-values based on resampling methods

PIMA (Girardi et al. 2024)

PIMA on meta-analysis

We implemented PIMA also on meta-analysis but (for the moment) limited to two-level models without moderators.

PIMA Results

multi_flip_overall <- flip::npc(multi_flip, comb.funct = "maxT")
multi_p <- multi_flip_overall@res$`p-value`
fp <- function(p){
    ifelse(p <= 0.001, "p < 0.001", as.character(round(p, 3)))
}

The multiverse is associated with an overall p value of 0.016 ¹. Then we can describe the overall results:

p_beta <- ggplot(multi, aes(y = b)) +
    xlim(c(-0.5,0.5)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("$\\mu_{\\,\\, \\theta}$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank())

p_pval <- ggplot(multi, aes(y = -log10(p.value))) +
    xlim(c(-0.5,0.5)) +
    geom_hline(yintercept = -log10(0.05)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("Raw $-log_{10}(p)$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank())

rdata <- cor(zi)
rdata <- data.frame(r = rdata[upper.tri(rdata, diag = FALSE)])

p_cor <- ggplot(rdata, aes(y = r)) +
    xlim(c(-0.5,0.5)) +
    geom_boxplot(width = 0.5,
                 fill = "dodgerblue",
                 alpha = 0.6) + 
    ylab(latex2exp::TeX("$\\rho$")) +
    theme_minimal(30) +
    theme(axis.text.x = element_blank()) +
    ylim(c(-1,1))

cowplot::plot_grid(p_beta, p_cor, p_pval, nrow = 1)

Impact of multiplicity correction

tp <- function(p){
    -log10(p)
}

lbl <- c("Never $p \\leq 0.05$", 
         "Before correction $p \\leq 0.05$", 
         "After correction $p \\leq 0.05$")


multi$sign <- factor(multi$sign, 
                     labels = latex2exp::TeX(lbl))

ggplot(multi, aes(x = tp(p.value), tp(adjust.maxt), color = sign)) +
    geom_hline(yintercept = tp(0.05), alpha = 0.4) +
    geom_vline(xintercept = tp(0.05), alpha = 0.4) +
    geom_abline(linetype = "dashed", lwd = 0.5) +
    geom_point(size = 5,
               position = position_jitter(width = 0.1),
               alpha = 0.5) +
    xlim(c(0, 3)) +
    ylim(c(0, 3)) +
    xlab(TeX("Raw p value ($-log_{10}$)")) +
    ylab(TeX("maxT p value ($-log_{10}$)")) +
    theme_minimal(base_size = 20) +
    theme(legend.title = element_blank(),
          legend.position = "bottom") +
    scale_color_manual(labels = scales::parse_format(), values = c("#F8766D", "#7CAE00", "#00BFC4"))