Extracting and Visualising the results of a multiverse analysis

Model #2: Effect of Fertility and Relationship status on Religiosity

Steegen et al. analyse the data using 6 models. All the models except the first use this dataset. For the visualisations, we will focus on the effect of Fertility and Relationship status on Religiosity.

The authors perform an ANOVA to study the effect of Fertility, Relationship and their interaction term, on the composite Religiosity score (RelComp). We perform this analysis and extract the results from the multiverse into a tidy data frame.

inside(M, {
  fit_RelComp <- lm( RelComp ~ Fertility * RelationshipStatus, data = df )
  
  fit_FiscConsComp <- lm( FiscConsComp ~ Fertility * RelationshipStatus, data = df)
  
  fit_SocConsComp <- lm( SocConsComp ~ Fertility * RelationshipStatus, data = df)

  fit_Donate <- glm( Donate ~ Fertility * Relationship, data = df, family = binomial(link = "logit") )
  
  fit_Vote <- glm( Vote ~ Fertility * Relationship, data = df, family = binomial(link = "logit") )
})

In most analysis, you would then look at the results using the summary function. However, when performing a multiverse analysis, the summary would be executed in each universe of the multiverse, and would not be accessible outside each universe. However, there are ways to easily extract the results from each universe, as long as the results within each universe are as tidy data.

Computing the results in each multiverse

The broom package allows extracting the summary of most linear models using the tidy() function. The tidy method summarises information about the components of a model including the estimates, std. error, t-statistic, p-value for each coefficient and returns a tibble (a tidy data frame). Additional arguments can be passed to obtain the upper and lower 95% confidence limits, or to change the confidence level.

inside(M, {
  summary_RelComp <- fit_RelComp %>% 
    broom::tidy( conf.int = TRUE )
  
  summary_FiscConsComp <- fit_FiscConsComp %>% 
    broom::tidy( conf.int = TRUE )
  
  summary_SocConsComp <- fit_SocConsComp %>% 
    broom::tidy( conf.int = TRUE )
  
  summary_Donate <- fit_Donate %>% 
    broom::tidy( conf.int = TRUE )
  
  summary_Vote <- fit_Vote %>% 
    broom::tidy( conf.int = TRUE )
})

Finally, we will execute all the generated analysis combinations in the multiverse, since so far we have only defined them.

execute_multiverse(M)

Extracting the results

All the computed summaries (summary_RelComp, summary_FiscConsComp, summary_SocConsComp, summary_Donate and summary_Vote) are stored within .results column (a separate environment) in the multiverse table. We can extract these summaries as a tibble using and store it in a separate column, and then unnest this new column:

expand(M) %>%
  mutate( summary = map(.results, "summary_RelComp") ) %>%
  unnest( summary )

## # A tibble: 840 × 17
##    .universe cycle_length certainty   menstrual_calcul… fertile relationship_st…
##        <int> <chr>        <chr>       <chr>             <chr>   <chr>           
##  1         1 cl_option1   cer_option1 mc_option1        fer_op… rs_option1      
##  2         1 cl_option1   cer_option1 mc_option1        fer_op… rs_option1      
##  3         1 cl_option1   cer_option1 mc_option1        fer_op… rs_option1      
##  4         1 cl_option1   cer_option1 mc_option1        fer_op… rs_option1      
##  5         2 cl_option1   cer_option1 mc_option1        fer_op… rs_option2      
##  6         2 cl_option1   cer_option1 mc_option1        fer_op… rs_option2      
##  7         2 cl_option1   cer_option1 mc_option1        fer_op… rs_option2      
##  8         2 cl_option1   cer_option1 mc_option1        fer_op… rs_option2      
##  9         3 cl_option1   cer_option1 mc_option1        fer_op… rs_option3      
## 10         3 cl_option1   cer_option1 mc_option1        fer_op… rs_option3      
## # … with 830 more rows, and 11 more variables: .parameter_assignment <list>,
## #   .code <list>, .results <list>, .errors <list>, term <chr>, estimate <dbl>,
## #   std.error <dbl>, statistic <dbl>, p.value <dbl>, conf.low <dbl>,
## #   conf.high <dbl>

We can then use various R packages to visualise the results. Below, we show some standard visualisations.

Animating the coefficients of the model from each multiverse

One way to visualise over the multiverse is to animate over the results from each universe in the multiverse. This approach, inspired by the concept of Hypothetical Outcome plots was described by Dragicevic et al. in Increasing the Transparency of Research Papers with Explorable Multiverse Analyses. This approach allows us to quickly see the robustness of a result — if a particular result is consistent across all analysis paths or idiosyncratic to a specific analysis path.

p <- expand(M) %>%
  mutate( summary_RelComp = map(.results, "summary_RelComp") ) %>%
  unnest( cols = c(summary_RelComp) ) %>%
  mutate( term = recode( term, 
                 "RelationshipStatusSingle" = "Single",
                 "Fertilitylow:RelationshipStatusSingle" = "Single:Fertility_low"
  ) ) %>%
  filter( term != "(Intercept)" ) %>%
  ggplot() + 
  geom_vline( xintercept = 0,  colour = '#979797' ) +
  geom_point( aes(x = estimate, y = term)) +
  geom_errorbarh( aes(xmin = conf.low, xmax = conf.high, y = term), height = 0) +
  theme_minimal() +
  transition_manual( .universe )

animate(p, nframes = 210, fps = 2)

Animation of linear regression coefficients

Histogram of p-values of the Fertility × Relationship status interaction

Another approach would be to summarise the results. Steegen et al. which depicts a histogram of p values of the Fertility × Relationship status interaction on religiosity for the multiverse of 210 data sets in Study 2 (summary_RelComp), on fiscal and social political attitudes for the multiverse of 210 data sets in Study 2 (summary_FiscConsComp and summary_SocConsComp), and on voting and donation preferences for the multiverse of 210 data sets in Study 2 (summary_Donate and summary_Vote). The dashed line indicates p = 0.05.

Below we show how to re-create the visualisation.

expand(M) %>%
  mutate( index = seq(1:nrow(.)) ) %>%
  mutate( 
    summary_RelComp = map(.results, "summary_RelComp" ),
    summary_FiscConsComp = map(.results, "summary_FiscConsComp" ),
    summary_SocConsComp = map(.results, "summary_SocConsComp" ),
    summary_Donate = map(.results, "summary_Donate" ),
    summary_Vote = map(.results, "summary_Vote" )
  ) %>%
  select( summary_RelComp:summary_Vote ) %>%
  gather( "analysis", "result" ) %>%
  unnest(result) %>%
  filter( term == "Fertilitylow:RelationshipStatusSingle" | term == "Fertilitylow:Relationship") %>%
  ggplot() +
  geom_histogram(aes(x = p.value), bins = 100, fill = "#ffffff", color = "#333333") +
  geom_vline( xintercept = 0.05, color = "red", linetype = "dashed") +
  facet_wrap(~ analysis, scales = "free", nrow = 3)

Specification curve

Simonsohn et al. propose visualising results from the multiverse as a specification curve, which consists of two panels. The top panel shows the effect size (here the value of a coefficient) for each universe (or specification). The bottom panel shows which specification of the parameters results in that result.

data.spec_curve <- expand(M) %>%
  mutate( summary_RelComp = map(.results, "summary_RelComp") ) %>%
  unnest( cols = c(summary_RelComp) ) %>%
  filter( term == "Fertilitylow:RelationshipStatusSingle" ) %>%
  select( .universe, !! names(parameters(M)), estimate, p.value ) %>%
  arrange( estimate ) %>%
  mutate( .universe = 1:nrow(.))

p1 <- data.spec_curve %>%
  gather( "parameter_name", "parameter_option", !! names(parameters(M)) ) %>%
  select( .universe, parameter_name, parameter_option) %>%
  mutate(
    parameter_name = factor(str_replace(parameter_name, "_", "\n"))
  ) %>%
  ggplot() +
  geom_point( aes(x = .universe, y = parameter_option, color = parameter_name), size = 0.5 ) +
  labs( x = "universe #", y = "option included in the analysis specification") + 
  facet_grid(parameter_name ~ ., space="free_y", scales="free_y", switch="y") +
  theme(strip.placement = "outside",
        strip.background = element_rect(fill=NA,colour=NA),
        panel.spacing.x=unit(0.15,"cm"), 
        strip.text.y = element_text(angle = 180, face="bold", size=10), 
        panel.spacing = unit(0.25, "lines")
      )

p2 <- data.spec_curve %>%
  ggplot() +
  geom_point( aes(.universe, estimate, color = (p.value < 0.05)), size = 0.25) +
  labs(x = "", y = "coefficient of\ninteraction term:\nfertility x relationship")
  
cowplot::plot_grid(p2, p1, axis = "bltr",  align = "v", ncol = 1, rel_heights = c(1, 3))

Specification curve visualisation