This case study is adapted from the multiverse package case study of the same analysis.

Initiate a multiverse

library(dplyr)
library(tidyr)
library(purrr)
library(tidybayes)
library(ggplot2)
library(mverse)
data('hurricane')
# read and process data
hurricane <- hurricane %>%
  rename(
    year = Year,
    name = Name,
    dam = NDAM,
    death = alldeaths,
    female = Gender_MF,
    masfem = MasFem,
    category = Category,
    pressure = Minpressure_Updated_2014,
    wind = HighestWindSpeed
  ) %>%
  mutate(
    post = if_else(year > 1979, 1, 0),
    zcat = as.numeric(scale(category)),
    zmin = -scale(pressure),
    zwin = as.numeric(scale(wind)),
    z3 = as.numeric((zmin + zcat + zwin) / 3)
  )
mv <- create_multiverse(hurricane)

Define filter branches

death_outliers <- filter_branch(
  TRUE,
  name != 'Katrina',
  !(name %in% c('Katrina', 'Audrey'))
)
damage_outliers <- filter_branch(
  TRUE,
  name != 'Sandy',
  !(name %in% c('Sandy', 'Andrew')),
  !(name %in% c('Sandy', 'Andrew', 'Donna'))
)
mv <- mv %>%
  add_filter_branch(death_outliers, damage_outliers)
summary(mv)
## # A tibble: 12 x 3
##    universe death_outliers_branch           damage_outliers_branch              
##    <fct>    <fct>                           <fct>                               
##  1 1        "TRUE"                          "TRUE"                              
##  2 2        "TRUE"                          "name != \"Sandy\""                 
##  3 3        "TRUE"                          "!(name %in% c(\"Sandy\", \"Andrew\…
##  4 4        "TRUE"                          "!(name %in% c(\"Sandy\", \"Andrew\…
##  5 5        "name != \"Katrina\""           "TRUE"                              
##  6 6        "name != \"Katrina\""           "name != \"Sandy\""                 
##  7 7        "name != \"Katrina\""           "!(name %in% c(\"Sandy\", \"Andrew\…
##  8 8        "name != \"Katrina\""           "!(name %in% c(\"Sandy\", \"Andrew\…
##  9 9        "!(name %in% c(\"Katrina\", \"… "TRUE"                              
## 10 10       "!(name %in% c(\"Katrina\", \"… "name != \"Sandy\""                 
## 11 11       "!(name %in% c(\"Katrina\", \"… "!(name %in% c(\"Sandy\", \"Andrew\…
## 12 12       "!(name %in% c(\"Katrina\", \"… "!(name %in% c(\"Sandy\", \"Andrew\…

Define mutate branches

femininity <- mutate_branch(
  female,
  masfem
)
damage <- mutate_branch(
  dam,
  log(dam)
)
mv <- mv %>%
  add_mutate_branch(femininity, damage)
summary(mv)
## # A tibble: 48 x 5
##    universe death_outliers_b… damage_outliers_br… femininity_bran… damage_branch
##    <fct>    <fct>             <fct>               <fct>            <fct>        
##  1 1        TRUE              "TRUE"              female           dam          
##  2 2        TRUE              "TRUE"              female           log(dam)     
##  3 3        TRUE              "TRUE"              masfem           dam          
##  4 4        TRUE              "TRUE"              masfem           log(dam)     
##  5 5        TRUE              "name != \"Sandy\"" female           dam          
##  6 6        TRUE              "name != \"Sandy\"" female           log(dam)     
##  7 7        TRUE              "name != \"Sandy\"" masfem           dam          
##  8 8        TRUE              "name != \"Sandy\"" masfem           log(dam)     
##  9 9        TRUE              "!(name %in% c(\"S… female           dam          
## 10 10       TRUE              "!(name %in% c(\"S… female           log(dam)     
## # … with 38 more rows

Define model branches

mverse methods aren’t implemented yet

# wrapper for model branching to be implemented
multiverse::inside(
  mv, {
    fit <- glm(
      branch(
        model,
        'linear' ~ log(death + 1),
        'poisson' ~ death
      ) ~ branch(
        femininity_damage_interaction,
        'yes' ~ femininity * damage,
        'no' ~ femininity + damage
      ) + branch(
        other_predictors,
        'pressure' ~ femininity * zmin,
        'wind' ~ femininity * zwin,
        'category' ~ femininity * zcat,
        'all' ~ femininity * z3,
        'all_no_interaction' ~ z3,
        'none' ~ NULL
      ) + branch(
        covariates,
        '1' ~ NULL,
        '2' ~ year:damage,
        '3' ~ post:damage
      ), family = branch(
        model,
        'linear' ~ gaussian,
        'poisson' ~ poisson
      ), data = data
    )
    # attach fitted values and standard errors to the original data
    pred <- predict(fit, se.fit = TRUE, type = "response")
    disagg_fit <- data %>%
      mutate(
        fitted = pred$fit,
        se.fit = pred$se.fit,
        df = df.residual(fit),
        sigma = sigma(fit),
        se.residual = sqrt(sum(residuals(fit)^2))/ df
      )
    pred2expectation <- function(mu, sigma) {
      branch(model, "linear" ~ exp(mu + sigma^2/2) - 1, "poisson" ~ mu)
    }
    expectation <- disagg_fit %>%
      mutate(expected_deaths = pred2expectation(fitted, sigma)) %>%
      group_by(female) %>%
      summarise(mean_deaths = mean(expected_deaths), .groups = "drop_last") %>%
      compare_levels(mean_deaths, by = female)
  }
)

Execute and Inspect Results

# execute the multiverse
mv <- mv %>%
  execute_multiverse()

mverse methods aren’t implemented yet (will wait until multiverse::extract_variables is available in the released version of multiverse)

# a wrapper function to be implemented using multiverse::extract_variables
fem_coefs <- multiverse::extract_variables(mv) %>%
  mutate(res = map(.results, 'expectation')) %>%
  unnest(res) %>%
  rename(universe = .universe)
fem_coefs %>%
  arrange(mean_deaths) %>%
  mutate(index = 1:nrow(.)) %>%
  ggplot(aes(y = mean_deaths, x = index, color = model)) +
  geom_point(alpha = 0.1) +
  xlab('Universe') + ylab('Mean Difference in Expected Deaths') +
  theme_minimal()
Estimated coefficient of femininity of a hurricane's name.

Estimated coefficient of femininity of a hurricane’s name.