analysis.Rmd

---
title: "Behavioral Plasticity Analysis"
date: "`r format(Sys.time(), '%d %B, %Y')`"
author: Fabian Dablander
output:
  html_document:
    toc: true
    theme: united
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE, warning = FALSE, message = FALSE, eval = TRUE,
  fig.align = 'center', fig.width = 10, fig.height = 7, dpi = 300,
  out.width='100%', out.height='100%'
)
```

<!-- avoid border around images -->
<style>
    img {
        border: 0;
    }
</style>

# Data preparation
This document includes code to reproduce all analyses and figures in the paper Nielsen, Dablander, Debnath, Emegor, Ghai, Gwozdz, Hahnel, Hofmann, & Bauer (2024). Perceived plasticity of climate-relevant behaviors and policy support among high- and lower-income individuals.

```{r}
library(rlang)
library(knitr)
library(Hmisc)
library(dplyr)
library(ggplot2)
library(corrplot)
library(tidyverse)
library(doParallel)
library(BayesFactor)
library(RColorBrewer)
source('helpers.R')

df <- readRDS('Data/BPdata.RDS') %>%
  mutate(
    country = factor(country, c('Denmark', 'United States', 'Nigeria', 'India'))
  )
df$climate_concern <- rowMeans(df[,c('cc_worry', 'cc_importance')], na.rm = TRUE)

actions_order <- c(
  # Curtailment
  'Use less energy for heating and cooling your home',
  'Drive fewer kilometers/miles in your car',
  'Take fewer national and international flights',
  'Eat less red meat',
  'Eat less white meat',
  
  # Investment
  'Improve the insulation of your home',
  'Replace older appliances with newer energy efficient models',
  'Install solar panels at home',
  'Purchase an electric vehicle',
  'Move private investments to climate-friendly financial products'
)
```

# Main Figures
## Figure 1: Behavioral plasticity across countries
We look at Behavioral plasticity across countries.

```{r, fig.height = 8}
varnames <- c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))
map_names <- label(df[varnames])

d <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit()

d <- d %>% 
  group_by(country, name, value) %>% 
  summarize(count = n()) %>% 
  left_join(
    d %>% group_by(country, name) %>% summarize(total_count = n()),
    by = c('country', 'name')
  ) %>% 
  mutate(
    prop = count / total_count
  ) %>% 
  mutate(
    value_name = factor(
      ifelse(value == 1, '1: Very unlikely', ifelse(value == 5, '5: Very likely', value)),
      levels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
    ),
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

cols <- c('#ff7f0e', '#A65628', '#029E73', '#56B4E9', '#0072B2')
p_country <- ggplot(d, aes(x = country, y = prop, fill = as.factor(value_name))) +
  geom_bar(stat = 'identity') +
  facet_wrap(~ full_name, ncol = 5) + 
  xlab('') +
  ylab('Proportion') +
  scale_fill_manual(
    values = rev(cols)
    # labels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
  ) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid = element_blank(),
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90, hjust = 1),
    plot.title = element_text(size = 14, hjust = 0.50)
  ) +
  guides(fill = guide_legend(reverse = TRUE))

ggsave('Figures/Figure1.pdf', p_country, width = 9, height = 7)
p_country
```

Calculate proportions for describing the figure.

```{r}
library(pander)

d_comp <- d %>% 
  group_by(country, full_name) %>% 
  summarize(
    unlikely = round(sum(prop[value %in% c(1, 2)]), 2),
    likely = round(sum(prop[value %in% c(4, 5)]), 2)
  ) %>% data.frame

pander(d_comp)
```

We use the `BayesFactor` R package for statistical inference.

## Figure 2: Behavioral plasticity and income
We show the relationship between behavioral plasticity and income across items and countries.

```{r}
df_preds <- do.call('rbind', lapply(varnames, function(varname) get_df_pred_bf(varname, df))) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_bp_income <- ggplot(df_preds, aes(y = ypred, x = combined_income, fill = country)) +
  geom_vline(aes(xintercept = 10), linetype = 'dotted', color = 'gray48') +
  geom_point(aes(y = mean_value, x = combined_income, col = country), size = 1, show.legend = FALSE) +
  geom_line(aes(col = country), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40) + 
  xlab('Income group') +
  ylab('Perceived behavioral plasticity') +
  scale_x_continuous(breaks = c(seq(1, 15, 3), 15), limits = c(1, 15)) +
  scale_y_continuous(breaks = seq(1, 5), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  scale_fill_manual(values = country_colors) +
  facet_wrap(~ full_name) +
  theme_minimal() +
  facet_wrap(~ full_name, ncol = 5) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(size = 8),
    legend.title = element_blank(),
    legend.position = 'top'
  )

ggsave('Figures/Figure2.pdf', p_bp_income, width = 10, height = 6)
p_bp_income
```

We calculate inclusion Bayes factors in favor of an effect.

```{r}
df_inclusion_bfs <- get_bf_inclusion_fig2_df(actions_order)
pander(df_inclusion_bfs %>% select(-rscaleFixed, -rscaleCont))
```

## Figure 3: Behavioral plasticity and domain-matched policy support
We run the relevant three-way interaction models.

```{r}
df_meat_narm <- df %>%
  filter(!is.na(support_7), !is.na(curtailment_2)) %>% 
  mutate(
    curtailment_2 = curtailment_2 - mean(curtailment_2),
    combined_income = combined_income - mean(combined_income)
  )

df_flights_narm <- df %>%
  filter(!is.na(support_8), !is.na(curtailment_1)) %>% 
  mutate(
    curtailment_1 = curtailment_1 - mean(curtailment_1),
    combined_income = combined_income - mean(combined_income)
  )
  
df_finance_narm <- df %>%
  filter(!is.na(support_6), !is.na(investment_4)) %>% 
  mutate(
    investment_4 = investment_4 - mean(investment_4),
    combined_income = combined_income - mean(combined_income)
  )

rscaleFixed <- sqrt(1/2)
rscaleCont <- sqrt(1/2) / 4

# Lose 422 rows
fit_meat <- generalTestBF(
  support_7 ~ curtailment_2 * combined_income * country,
  data = df_meat_narm, whichModels = 'withmain', 
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 443 rows
fit_flights <- generalTestBF(
  support_8 ~ curtailment_1 * combined_income * country,
  data = df_flights_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 530 rows
fit_finance <- generalTestBF(
  support_6 ~ investment_4 * combined_income * country,
  data = df_finance_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Get the model that predicts the data best
fit_meat_best <- fit_meat[which.max(fit_meat@bayesFactor[, 1])]
fit_flights_best <- fit_flights[which.max(fit_flights@bayesFactor[, 1])]
fit_finance_best <- fit_finance[which.max(fit_finance@bayesFactor[, 1])]

## Used below
countries <- c('Denmark', 'United States', 'Nigeria', 'India')
newdat <- expand.grid(
  country = countries,
  combined_income = seq(15) - mean(seq(15)),
  bp_value = seq(5) - mean(seq(5))
)

map_back <- function(x) {
  as.numeric(factor(x, labels = seq(5)))
}
```

We look at model predictions using the models with the highest marginal likelihood.

```{r}
matches <- list(
  c('support_8', 'curtailment_1'),
  c('support_7', 'curtailment_2'),
  c('support_6', 'investment_4')
)

df_preds <- bind_rows(
  get_df_pred_policy(fit_meat_best, 'curtailment_2', 'support_7', df_meat_narm),
  get_df_pred_policy(fit_flights_best, 'curtailment_1', 'support_8', df_flights_narm),
  get_df_pred_policy(fit_finance_best, 'investment_4', 'support_6', df_finance_narm)
) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(map_names)
    ),
    
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
  )
```

```{r, fig.height = 8}
matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_matching <- ggplot(df_preds, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_y_continuous(breaks = seq(1, 7, 1)) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/Figure3.pdf', p_matching, width = 9, height = 8)
p_matching
```

```{r}
library(VGAM)

model_size_prior <- dbetabinom.ab(seq(0, 7), size = 7, shape1 = 1, shape2 = 1)
x <- model_size_prior[1]

# We specify the priors here according to the sequence from the output of the BayesFactor package
# This is if we wanted to use the betabinomial prior, but for the paper we report a uniform model prior
# There's barely any difference
model_prior <- c(
  # Size 0: 1 --> grand mean only
  # Size 1: 3 --> one main effect
  # Size 2: 3 --> two main effects
  # Size 3: 4 --> all main effects or two main effects and one interaction
  # Size 4: 3 --> all main effects and one interactions
  # Size 5: 3 --> all main effects and two two-way interactions
  # Size 6: 1 --> all main effects and all two-way interactions
  # Size 7: 1 --> all main effects and all interactions
  x,
  rep(x / 3, 3),
  rep(x / 3, 3),
  rep(x / 4, 4),
  rep(x / 3, 3),
  rep(x / 3, 3),
  x,
  x
)

set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4')
))

colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)
```

Not so different with uniform prior, which we report in the main paper.

```{r}
set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
))
colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)
```

# Extended Data Figures
## Figure ED1: Already did the action or cannot do it
We visualize the proportion of people who already did a particular action or cannot do it.

```{r, fig.height = 8}
varnames_peb <- paste0(varnames, '_peb')
varnames_cannot <- paste0(varnames, '_cannot')
map_names_cannot <- map_names
names(map_names_cannot) <- varnames_cannot

df_cannot <- df %>%
  select(country, all_of(varnames_cannot)) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == 1), # Only 0 / 1 exists
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_cannot[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    prop = ifelse(grepl('curtailment', name), NA, prop),
    prop_se = ifelse(grepl('curtailment', name), NA, prop_se)
  )

varnames_peb <- paste0(varnames, '_peb')
map_names_peb <- map_names
names(map_names_peb) <- varnames_peb

df_peb <- df %>%
  select(country, all_of(varnames_peb)) %>%
  pivot_longer(-country) %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df %>%
      select(country, all_of(varnames_peb)) %>%
      pivot_longer(-country) %>% 
      mutate(value = as.numeric(as.character(value))) %>% 
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_both <- bind_rows(
  df_cannot %>% mutate(type = 'Cannot do'),
  df_peb %>% mutate(type = 'Already did / Never do it')
) %>% 
  mutate(
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)],
        'Do not eat white meat',
        'Do not eat red meat',
        'Never fly',
        'Do not own a car',
        'Do not heat or cool their home'
      )
    )
  )

pwidth <- 0.96
p_both <- ggplot(df_both, aes(x = prop, y = full_name, fill = country, col = country)) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
    geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ type) +
  ylab('') +
  xlab('Proportion') +
  scale_fill_manual(values = country_colors) +
  scale_color_manual(values = country_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 10),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureED1.pdf', p_both, width = 9, height = 8)
p_both
```

## Figure ED2: Already did and behavioral plasticity
Visualize relationship between exclusion rate and mean plasticity rating for behaviors across countries.

```{r, fig.height = 8}
df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    se_bp = sd_bp / sqrt(n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_mean_exclusion <- df_mean %>% 
  left_join(
    df_both %>%
      filter(type == 'Already did / Never do it') %>% 
      mutate(
        name = str_extract(name, "^[^_]+_[^_]+")
      ) %>% 
      select(country, name, prop, prop_se)
    , by = c('country', 'name')
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = rev(wrap_label_vectorized(actions_order))
    )
  )
  
cols <- brewer.pal(10, 'Set3')
p_mean_exclusion <- ggplot(df_mean_exclusion, aes(x = prop, y = mean_bp, col = full_name)) +
  geom_point(size = 2) +
  geom_errorbar(
    aes(ymin = mean_bp - 1.96 * se_bp, ymax = mean_bp + 1.96 * se_bp),
    width = 0.01
  ) +
  # Horizontal error bars
  geom_errorbarh(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    height = 0.1
  ) +
  facet_wrap(~ country) +
  labs(
    x = 'Proportion: Already do / never do it',
    y = 'Mean perceived behavioral plasticity'
  ) +
  scale_color_manual(name = '', values = cols) +
  scale_x_continuous(breaks = seq(0, 0.3, 0.05)) +
  scale_y_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.10, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureED2.pdf', p_mean_exclusion, width = 7, height = 7)
p_mean_exclusion
```

## Figure ED3: Already did the action across income groups
```{r}
df_prep <- df %>%
  select(country, combined_income, varnames_peb) %>%
  pivot_longer(cols = -c(country, combined_income)) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    )
  )

df_selection <- df_prep %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  group_by(country, income_bracket, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df_prep %>%
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, income_bracket, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    ),
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        'Do not heat or cool their home',
        'Do not own a car',
        'Never fly',
        'Do not eat red meat',
        'Do not eat white meat',
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)]
      )
    )
  )

p_selection <- ggplot(
    df_selection,
    aes(y = prop, x = country, fill = income_bracket)
  ) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
  geom_point(
    aes(x = country, y = prop), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = country, y = prop, col = income_bracket), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se, col = income_bracket),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('Proportion') +
  xlab('') +
  scale_fill_manual(values = matching_colors) +
  scale_color_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureEx3.pdf', p_selection, width = 9, height = 7)
p_selection
```

Let's run a proportion test to assess the evidence for differences.

```{r}
df_bfs <- c()
for (country in unique(df_prep$country)) {
  for (name in unique(df_prep$name)) {
    bf <- test_prop(df, country, name)
    df_bfs <- rbind(df_bfs, cbind(bf, country, name))
  }
}

df_bfs <- data.frame(df_bfs) %>% 
  mutate(
    name = gsub('_peb$', '', name),
    log_bf = as.numeric(bf),
    full_name = map_names[name]
  )

pander(df_bfs %>% select(country, full_name, log_bf))
```

# Supplementary Figures
## Figure S1: Mean behavioral plasticity across countries
Here we visualize mean behavioral plasticity across all countries.

```{r}
df_mean_sd_beh <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  mutate(curtailment = as.numeric(grepl('curtailment', name))) %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n)
  ) %>% 
  # Compute the mean of the mean /sd behavioral plasticity of behaviors
  group_by(country) %>% 
  summarize(
    mean_bp = mean(mean_bp),
    sd_bp = mean(sd_bp)
  )

pander(df_mean_sd_beh)
```

```{r}
df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n),
    se_sd_bp = sd_bp / sqrt(2 * n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

pwidth <- 0.91
p_mean <- ggplot(df_mean, aes(x = mean_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_bp - 1.96 * se_bp, xmax = mean_bp + 1.96 * se_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS1.pdf', p_mean, width = 9, height = 7)
p_mean
```

## Figure S1: Standard deviation of behavioral plasticity across countries
```{r}
pwidth <- 0.91
p_sd <- ggplot(df_mean, aes(x = sd_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = sd_bp - 1.96 * se_sd_bp, xmax = sd_bp + 1.96 * se_sd_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Standard deviation of perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(0, 2, 0.50), limits = c(0, 2)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS2.pdf', p_sd, width = 9, height = 7)
p_sd
```

## Figure S3: Mean policy support across countries
Here we visualize mean policy support across all countries.

```{r}
support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_mean_ps <- df %>%
  select(country, starts_with('support_')) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_ps = mean(value),
    sd_ps = sd(value),
    n = n(),
    se_ps = sd_ps / sqrt(n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(support_labels[name]),
      levels = wrap_label_vectorized(support_labels[name])
    )
  )

pwidth <- 0.91
country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_mean_ps <- ggplot(df_mean_ps, aes(x = mean_ps, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_ps - 1.96 * se_ps, xmax = mean_ps + 1.96 * se_ps),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean policy support') +
  scale_x_continuous(breaks = seq(1, 7, 1), limits = c(1, 7)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 6),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS3.pdf', p_mean_ps, width = 9, height = 8)
p_mean_ps
```

## Figure S4: Differences in high and low income groups
Here we visualize the differences in mean behavioral plasticity between the top 10% income earners and the rest.

```{r}
set.seed(1)
df_diffs <- do.call(
  'rbind', lapply(varnames, function(varname) get_top10_difference(varname, df))
) %>% 
  mutate(
    country = factor(country, levels = rev(countries)),
    full_name = factor(
      map_names[name],
      levels = rev(actions_order)
    )
  )
```

```{r, fig.height = 10}
cols <- brewer.pal(10, 'Set3')
p_estimates <- ggplot(df_diffs, aes(x = mean_diff, y = country, color = full_name)) +
  geom_vline(aes(xintercept = 0), linetype = 'dotted', color = 'gray48') +
  geom_errorbar(
    aes(xmin = mean_diff_lo, xmax = mean_diff_hi),
    position = position_dodge(width = 0.70), width = 0.80, size = 1.5
  ) +
  geom_point(position = position_dodge(width = 0.70), size = 3) + 
  scale_x_continuous(breaks = round(seq(-0.80, 1.2, 0.2), 3), limits = c(-0.90, 1.3)) +
  theme_minimal() +
  scale_color_manual(name = '', values = cols) +
  labs(
    x = 'Mean difference in perceived behavioral plasticity',
    y = ''
  ) +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.06, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureS4.pdf', p_estimates, width = 8, height = 10)
p_estimates
```

## Figures S5 & S6: Sensitivity analysis I
We assess how sensitive the results are to changes in the prior width of the income effect and the interaction of income and country.

```{r}
scale_lo <- 0.05
scale_hi <- 0.30
scales_cont <- seq(scale_lo, scale_hi, 0.0125)

df_inclusion_bfs_sens <- do.call('rbind', lapply(scales_cont, function(scale) {
  get_bf_inclusion_fig2_df(actions_order, rscaleCont = scale)
}))
end <- Sys.time()

df_inclusion_bfs_sens <- df_inclusion_bfs_sens %>% 
    mutate(
    full_name = factor(
      wrap_label_vectorized(full_name),
      levels = wrap_label_vectorized(actions_order)
    )
  )
```

```{r, fig.height = 6}
library(scales)

# Custom function to apply scientific notation for numbers >= 10
custom_label <- function(x) {
  ifelse(x > 100, scientific_format()(x), format(x, scientific = FALSE))
}

df_inclusion_bfs_sens2 <- df_inclusion_bfs_sens %>% 
  mutate(income = ifelse(income <= 1, 1 / income, income))

p_sens_fig2_income <- ggplot(df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income))) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS5.pdf',
  p_sens_fig2_income,
  width = 10, height = 5
)
p_sens_fig2_income
```

```{r, fig.height = 6}
p_sens_fig2_income_country <- ggplot(
  df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income_x_country))
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS6.pdf',
  p_sens_fig2_income_country,
  width = 10, height = 5
)
p_sens_fig2_income_country
```

## Figure S7: Sensitivity analysis II
We run sensitivity analyses for the results reported in Figure 3.

```{r}
if (!file.exists('Data/df_sensitivity_matched.csv')) {
  
  rscaleFixed <- sqrt(1/2)
  scales_cont <- seq(0.05, 0.30, 0.0125)
  
  df_inclusion_bfs <- do.call('rbind', lapply(scales_cont, function(scale) {
    
    set.seed(1)
    fit_meat <- generalTestBF(
      support_7 ~ curtailment_2 * combined_income * country,
      data = df_meat_narm, whichModels = 'withmain', 
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_flights <- generalTestBF(
      support_8 ~ curtailment_1 * combined_income * country,
      data = df_flights_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_finance <- generalTestBF(
      support_6 ~ investment_4 * combined_income * country,
      data = df_finance_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    df_inclusion_bfs <- data.frame(rbind(
      get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
    ))
    
    colnames(df_inclusion_bfs) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
    df_inclusion_bfs$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]
    df_inclusion_bfs$rscaleCont <- scale
    df_inclusion_bfs
    
  }))
  
  df_inclusion_bfs_sens2 <- df_inclusion_bfs %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(full_name),
        levels = wrap_label_vectorized(actions_order)
      )
    )
  
  df_inclusion_bfs_sens2_long <- df_inclusion_bfs_sens2 %>%
    pivot_longer(
      cols = c(bp, bp_country, bp_income, bp_income_country),
      names_to = 'type',
      values_to = 'value'
    ) %>% 
    mutate(
      type = case_when(
        type == 'bp' ~ 'BP',
        type == 'bp_country' ~ 'BP x Country',
        type == 'bp_income' ~ 'BP x Income',
        type == 'bp_income_country' ~ 'BP x Income x Country'
      )
    )
} else {
  df_inclusion_bfs_sens2_long <- read.csv('Data/df_sensitivity_matched.csv') %>% 
    mutate(
      full_name = factor(
        full_name, levels = c(
          'Take fewer national and\ninternational flights', 'Eat less red meat',
          'Move private investments to\nclimate-friendly financial products'
        )
      )
    )
}
```

```{r, fig.height = 9}
p_sens_fig3 <- ggplot(
  df_inclusion_bfs_sens2_long, aes(x = rscaleCont, y = value)
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = label_scientific()) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  # facet_nested(~ type + full_name, scales = 'free_y') +
  facet_wrap(~ type + full_name, ncol = 3, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS7.pdf',
  p_sens_fig3,
  width = 8, height = 8
)
p_sens_fig3
```


## Figures S8 & S9: More loosely domain-matched analyses
We estimate the relationship between more loosely domain-matched policies and behavioral plasticity items.

```{r}
get_best <- function(fit_all) fit_all[which.max(fit_all@bayesFactor[, 1])]

matches_loose <- list(
  
  # Increasing energy efficiency requirements for buildings & insulating your home
  c('support_5', 'investment_3'),
  
  # Increase price of electricity during peak hours & replace older with more energy efficient appliances
  c('support_3', 'investment_2'),
  
  # increasing subsidies for renewable energy projects & install solar panels at home
  c('support_4', 'investment_5'),
  
  # Ban sale of diesel & purchasing an EV
  c('support_10', 'investment_1'),
  
  # Requires banks to disclose & move investments
  c('support_6', 'investment_4'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_2'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_3'),
  
  # Expand public transport & drive fewer kilometres with the car
  c('support_2', 'curtailment_4'),
  
  # Strengthen requirements for energy efficiency & use less energy for heating and cooling
  c('support_5', 'curtailment_5')
)

#' Prepare data.frame for loosely domain-matched policy analysis
prepare_df <- function(df, bp_name, policy_name) {
  bp_name <- ensym(bp_name)
  policy_name <- ensym(policy_name)
  
  df %>% 
    filter(!is.na(!!bp_name), !is.na(!!policy_name)) %>% 
    mutate(
      !!bp_name := !!bp_name - mean(!!bp_name, na.rm = TRUE),
      combined_income = combined_income - mean(combined_income, na.rm = TRUE)
    )
}

if (!file.exists('Data/df_preds_loose.csv')) {
  
  fit_cars <- get_bf_fit_loose(df, 'support_10', 'investment_1')
  fit_electricity <- get_bf_fit_loose(df, 'support_3', 'investment_2')
  fit_banks <- get_bf_fit_loose(df, 'support_6', 'investment_4')
  fit_replace <- get_bf_fit_loose(df, 'support_5', 'investment_3')
  fit_renewables <- get_bf_fit_loose(df, 'support_4', 'investment_5')
  
  fit_food_red <- get_bf_fit_loose(df, 'support_12', 'curtailment_2')
  fit_food_white <- get_bf_fit_loose(df, 'support_12', 'curtailment_3')
  fit_public <- get_bf_fit_loose(df, 'support_2', 'curtailment_4')
  fit_efficiency <- get_bf_fit_loose(df, 'support_5', 'curtailment_5')
  
  df_preds_loose <- bind_rows(
    get_df_pred_policy(
      get_best(fit_cars), 'investment_1', 'support_10',
      prepare_df(df, 'investment_1', 'support_10')
    ) %>% mutate(type = 'cars'),
    
    get_df_pred_policy(
      get_best(fit_electricity), 'investment_2', 'support_3',
      prepare_df(df, 'investment_2', 'support_3')
    ) %>% mutate(type = 'electricity'),
    
    get_df_pred_policy(
      get_best(fit_banks), 'investment_4', 'support_6',
      prepare_df(df, 'investment_4', 'support_6')
    ) %>% mutate(type = 'banks'),
    
    get_df_pred_policy(
      get_best(fit_replace), 'investment_3', 'support_5',
      prepare_df(df, 'investment_3', 'support_5')
    ) %>% mutate(type = 'replace'),
    
    get_df_pred_policy(
      get_best(fit_renewables), 'investment_5', 'support_4',
      prepare_df(df, 'investment_5', 'support_4')
    ) %>% mutate(type = 'renewables'),
    
    get_df_pred_policy(
      get_best(fit_food_red), 'curtailment_2', 'support_12',
      prepare_df(df, 'curtailment_2', 'support_12')
    ) %>% mutate(type = 'red_meat'),
    
    get_df_pred_policy(
      get_best(fit_food_white), 'curtailment_3', 'support_12',
      prepare_df(df, 'curtailment_3', 'support_12')
    ) %>% mutate(type = 'white_meat'),
    
    get_df_pred_policy(
      get_best(fit_public), 'curtailment_4', 'support_2',
      prepare_df(df, 'curtailment_4', 'support_2')
    ) %>% mutate(type = 'public'),
    
    get_df_pred_policy(
      get_best(fit_efficiency), 'curtailment_5', 'support_5',
      prepare_df(df, 'curtailment_5', 'support_5')
    ) %>% mutate(type = 'public')
  ) %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
    )
} else {
  df_preds_loose <- read.csv('Data/df_preds_loose.csv') %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
      income_bracket = factor(income_bracket, levels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)'))
    )
}

df_preds_curtailment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'curtailment'))

df_preds_investment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'investment'))
```

```{r, fig.height = 10}
matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_investment_loose <- ggplot(
  df_preds_investment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

p_curtailment_loose <- ggplot(
  df_preds_curtailment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Behavioral plasticity') +
  ylab('Policy support') +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS8.pdf', p_curtailment_loose, width = 9, height = 10)
p_curtailment_loose
```

```{r, fig.height = 10}
ggsave('Figures/FigureS9.pdf', p_investment_loose, width = 9, height = 10)
p_investment_loose
```

## Figure S10: Mean behavioral plasticy and mean policy support
We run the model for mean support across all policies and mean of all behavioral plasticity items.

```{r, fig.height = 4}
df_avg <- df %>% 
  rowwise() %>% 
  mutate(
    support_avg = mean(c_across(starts_with('support')), na.rm = TRUE),
    bp_value = mean(c_across(c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))), na.rm = TRUE)
  ) %>% data.frame

fit_all <- generalTestBF(
  support_avg ~ bp_value * combined_income * country,
  data = df_avg, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

fit_all_best <- fit_all[which.max(fit_all@bayesFactor[, 1])]

df_preds_avg <- get_df_pred_policy(
  fit_all_best, 'bp_value', 'support_avg', df_avg, type = 'averaged',
  join = FALSE, use_empirical = FALSE
) %>% 
  mutate(
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_emp_avg <- df_avg %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    ),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_both <- df_emp_avg %>% 
  left_join(df_preds_avg, by = c('country', 'bp_value', 'income_bracket'))

p_mean_matched <- ggplot(df_both, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = support_avg, col = income_bracket),
    size = 1, alpha = 0.30, width = 0.10, show.legend = FALSE
  ) +
  geom_line(data = df_preds_avg, aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(data = df_preds_avg, aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  scale_y_continuous(breaks = seq(1, 7, 1), limits = c(0.9, 7.2)) +
  xlab('Mean perceived behavioral plasticity') +
  ylab('Mean policy support') +
  facet_wrap(~ country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS10.pdf', p_mean_matched, width = 9, height = 4)
p_mean_matched
```

## Supplementary Figures 11 - 14: Correlation plots
We calculate Kendall's correlation coefficients.

```{r}
support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_cors <- df %>% select(all_of(varnames), starts_with('support_'), country)
labels <- c(wrap_label_vectorized(map_names), wrap_label_vectorized(support_labels))

res_india <- get_cors_all(df_cors, labels, 'India')
res_denmark <- get_cors_all(df_cors, labels, 'Denmark')
res_nigeria <- get_cors_all(df_cors, labels, 'Nigeria')
res_usa <- get_cors_all(df_cors, labels, 'United States')
```

### Denmark
```{r, fig.width = 10, fig.height = 10}
pdf('Figures/FigureS11.pdf', width = 10, height = 10)
get_corrplot_all(res_denmark)
dev.off()
get_corrplot_all(res_denmark)
```

### United States
```{r, fig.width = 10, fig.height = 10}
pdf('Figures/FigureS12.pdf', width = 10, height = 10)
get_corrplot_all(res_usa)
dev.off()
get_corrplot_all(res_usa)
```

### Nigeria
```{r, fig.width = 10, fig.height = 10}
pdf('Figures/FigureS13.pdf', width = 10, height = 10)
get_corrplot_all(res_nigeria)
dev.off()
get_corrplot_all(res_nigeria)
```

### India
```{r, fig.width = 10, fig.height = 10}
pdf('Figures/FigureS14.pdf', width = 10, height = 10)
get_corrplot_all(res_india)
dev.off()
get_corrplot_all(res_india)
```