Setting Simulation Parameters

These parameters determine how the simulated sample will be generated. They are based on information that we received from the canteen operator and subject to change. CO2 equivalent measures are based on kgCO2e/kg food.

if (FALSE) {
  # To verify the CO2 values for parameter list
  food_co2 <- read_csv("../data/external/food_co2e.csv", col_types = cols())
  food_co2 %>%  mutate(meat_fish = category %in% c("meat", "fish")) %>%
    group_by(meat_fish) %>%
    summarise(mn_co2 = mean(co2_100g)/100, sd_co2 = sd(co2_100g)/100)
}


parms <- list(
  days = 10,
  treatment_slots_per_day = 3,
  mn_diners_per_day = 3250,
  sd_diners_per_day = 100,
  friday_churn = 0.75,
  monday_churn = 0.9,
  mn_mn_food_weight = 500,
  sd_mn_food_weight = 50,
  mn_sd_food_weight = 100,
  sd_sd_food_weight = 10,
  prob_meat_fish = 0.6,
  sd_meat_fish = 0.1,
  mn_mtfsh_weight = 200,
  sd_mtfsh_weight = 50,
  mn_co2e_veg = 2.2,
  sd_co2e_veg = 2.2,
  mn_co2e_mtfsh = 8.7,
  sd_co2e_mtfsh = 11.0
)

Setting up meal plan

We simulate a meal plan that contains two dishes: Vegetarian and Meat/Fish. Diners can choose whether they add a meat/fish component to their otherwise vegetarian dish. If they do, then the meat_fish indicator variable in the diner dataset is being set to true.

# We need to truncate the mtfsh and vegetable CO2 distribution to enforce
# that mtfsh always has higher CO2 equivalent per 100g than vegetables
# In the field, this will be verified by food menue design.

sim_meal_plan <- function() {
  tibble(
    day = 1:parms$days,
    veg_co2e = rtruncnorm(
      parms$days, 0.2, 3, mean = parms$mn_co2e_veg, sd = parms$sd_co2e_veg
    ),
    mtfsh_co2e = rtruncnorm(
      parms$days, 4, mean = parms$mn_co2e_mtfsh, sd = parms$sd_co2e_mtfsh
    ),
    mtfsh_weight = rtruncnorm(
      parms$days, 100, mean = parms$mn_mtfsh_weight, sd = parms$sd_mtfsh_weight
    )
  )
}

meal_plan <- sim_meal_plan()

Base Sample Descriptives

This is a draw of our base sample (without any treatment)

df <- sim_base_sample(parms, meal_plan)

kable(
  df %>% group_by(day) %>%
    summarise(
      diners = n(),
      meat_fish = mean(meat_fish),
      mn_food_weight = mean(food_weight),
      sd_food_weight = sd(food_weight),
      min_food_weight = min(food_weight),
      max_food_weight = max(food_weight),
      mn_co2e = mean(co2e),
      sd_co2e = sd(co2e),
      min_co2e = min(co2e),
      max_co2e = max(co2e)
    )
) %>% kable_styling()
day diners meat_fish mn_food_weight sd_food_weight min_food_weight max_food_weight mn_co2e sd_co2e min_co2e max_co2e
1 2892 0.5515214 514.7915 97.18396 148.7378 905.1886 1953.215 524.3140 598.07662 3461.395
2 3156 0.6913815 540.8745 94.27001 233.4578 851.7522 1183.907 549.3066 159.71925 1759.195
3 3149 0.7307082 445.3412 88.84015 182.2461 747.7525 3238.920 1780.4283 127.68126 4517.391
4 3235 0.5486862 510.4558 112.72951 108.9215 1035.8452 1515.383 1252.8747 28.76156 2789.702
5 2446 0.5461979 548.7076 100.95546 207.2900 896.9752 2182.359 1345.4919 338.64499 3830.194
6 2935 0.6517888 479.8999 99.93824 148.9004 803.8492 1201.477 453.5745 186.64606 1928.049
7 3216 0.6209577 528.0369 102.82407 188.7875 852.0572 1980.139 1163.1915 177.79470 3190.446
8 3275 0.5013740 462.3673 91.12968 147.0902 771.6921 5144.577 4451.0594 217.88473 10038.007
9 3303 0.6321526 473.0100 103.49692 165.0971 826.3704 2918.806 1497.9840 340.74448 4778.476
10 2415 0.6467909 506.2475 119.53877 146.4754 963.4136 842.937 332.4125 121.89664 1440.293 
ggplot(df) + geom_density(aes(x  = food_weight, color = meat_fish))+
  scale_color_trr266_d() +
  theme_trr() 

ggplot(df) + geom_density(aes(x  = co2e, color = meat_fish)) + 
  scale_color_trr266_d() +
  theme_trr()

modelsummary(
  feols(co2e ~ food_weight + meat_fish | day + slot, data = df),
  stars = TRUE
)
Model 1
food_weight 1.278***
(0.238)
meat_fishTRUE 2866.155***
(848.273)
Num.Obs. 30022
R2 0.717
R2 Adj. 0.717
R2 Within 0.586
R2 Pseudo
AIC 509282.2
BIC 509398.5
Log.Lik. -254627.093
FE: day X
FE: slot X
Std. errors Clustered (day)
* p < 0.1, ** p < 0.05, *** p < 0.01

Our effect sizes are informed by prior literature.

prior <- read_csv("../data/external/prior_results.csv", col_types = cols())

std_effects <- prior  %>%
  mutate(
    poseffect = sign_dummy * effect,
    std_effect = poseffect/baseline,
    lb_std_effect = poseffect/baseline - 1.96*stderr/baseline,
    ub_std_effect = poseffect/baseline + 1.96*stderr/baseline,
  ) %>%
  arrange(std_effect) %>%
  mutate(
    model = row_number()
  )

if (file.exists("../data/generated/power_runs.rds")) {
  df <- readRDS("../data/generated/power_runs.rds")
  effect_beyer_et_al_h1 <- df %>%
    filter(test == "h1" & mod == "weight") %>%
    summarise(
      area = "Power Results (H1)",
      authors = "Beyer et al. (TBD)",
      std_effect = mean(effect),
      lb_std_effect = std_effect - 1.96*sd(effect),
      ub_std_effect = std_effect + 1.96*sd(effect),
      model = nrow(std_effects) + 1,
    )
} else {
  # Power analysis not run yet - let's use canned values
  effect_beyer_et_al_h1 <- tibble(
      area = "Power Results (H1)",
      authors = "Beyer et al. (TBD)",
      std_effect = -0.0533,
      lb_std_effect = -0.0620,
      ub_std_effect = -0.0446,
      model = 27
    )
}

effect_sim <- bind_rows(std_effects, effect_beyer_et_al_h1) %>%
  arrange(model) %>%
  mutate(model = factor(model))

ggplot(effect_sim, aes(x = model)) +
  geom_pointrange(
    aes(
      y = std_effect, ymin = lb_std_effect, ymax = ub_std_effect,
      color = area
    ),
    fatten = 0.5
  ) + 
  labs(x = "", y = "% Effect to Baseline", color = "Treatment") +
  scale_y_continuous(labels = scales::percent) + 
  scale_x_discrete(labels = effect_sim$authors) +
  scale_color_manual(
    values = c(col_trr266_petrol, col_trr266_yellow, col_trr266_blue)
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )

save_plot("../output/figure_1_eff_sizes", 7, 7)

Design treatment plan

Now we specify a treatment pattern. The idea is to make sure that no treatment can be identified by day and slot fixed effects. This means that each treatment has to be administered at least twice and on varying days and on varying slots. We currently have

Assuming ten days and three slots per day, each treatment can be administered once per slot and on three varying days. We should also make sure that the treatment composition varies across days, meaning that the same treatments are not combined on several days. Finally, we also make sure that each day, both treatment arms are present.

color_cond <- c("Neutral", "ColorCoded")
info_cond <- c("NoNumericInfo", "NoContext", "Budget", "Car", "Money")
tments <- expand_grid(
  color_cond = factor(color_cond, levels = color_cond),
  info_cond = factor(info_cond, level = info_cond)
) %>%
  mutate(tment = sprintf("%s_%s", color_cond, info_cond)) 

tplan <- matrix(
  NA_character_, parms$days, parms$treatment_slots_per_day
)

# This could be brute forced or administered via a smart algorithm
# We did it manually for the time being

tment_str <- tments %>% pull(tment)

tplan[1, ] <- tment_str[c(1, 2, 6)]
tplan[2, ] <- tment_str[c(3, 7, 8)]
tplan[3, ] <- tment_str[c(4, 5, 9)]
tplan[4, ] <- tment_str[c(10, 6, 1)]
tplan[5, ] <- tment_str[c(2, 3, 7)]
tplan[6, ] <- tment_str[c(8, 9, 4)]
tplan[7, ] <- tment_str[c(5, 1, 10)]
tplan[8, ] <- tment_str[c(6, 8, 3)]
tplan[9, ] <- tment_str[c(7, 4, 2)]
tplan[10, ] <- tment_str[c(9, 10, 5)]

# Checks - Need each treatment to be present three times and each slot to 
# contain all ten treatments - this is the case

if (FALSE) {
  length(unique(tplan[, 1]))
  length(unique(tplan[, 2]))
  length(unique(tplan[, 3]))
}

table(tplan)
## tplan
##        ColorCoded_Budget           ColorCoded_Car         ColorCoded_Money 
##                        3                        3                        3 
##     ColorCoded_NoContext ColorCoded_NoNumericInfo           Neutral_Budget 
##                        3                        3                        3 
##              Neutral_Car            Neutral_Money        Neutral_NoContext 
##                        3                        3                        3 
##    Neutral_NoNumericInfo 
##                        3
# convert to tibble for later use

tplan <- tibble(
  day = rep(1:parms$days, each = parms$treatment_slots_per_day),
  slot = rep(1:parms$treatment_slots_per_day, parms$days),
  tment = as.vector(t(tplan))
) %>%
  left_join(tments, by = "tment")

Specify Expected Treatment Effects

Here we specify some treatment effect estimates. These can be varied to see how our power depends on these assumptions. We have three outcome measures

# To keep things simple for the time being, we estimate food gram effects in 
# percent consumption to be equal to percentage point effects for the like-
# likelihood to choose the more carbon intensive fish_meat option 
# Also, we set all effect standard deviations to zero which is in line with
# traditional power analyses

nocontext_eff <- -0.02
budget_eff <- -0.08
money_eff <- -0.04
car_eff <- -0.06
color_eff <- -0.05
sd_eff = 0.0

neutral_effects <- c(0, nocontext_eff, budget_eff, car_eff, money_eff)
color_effects <- neutral_effects + color_eff

eff_sizes <- tibble(
  tment = tments$tment,
  mn_weight = c(neutral_effects, color_effects),
  sd_weight = rep(sd_eff, nrow(tments)),
  mn_mtfsh = mn_weight,
  sd_mtfsh = sd_weight
)

Example Analysis

Here we show for one run of the simulation what our main tests would look like in the paper. As we have very strong daily effects, we report the results by day for descriptive purposes.

mp <- sim_meal_plan()
bs <- sim_base_sample(parms, mp)
df <- apply_teffects(bs, mp, tplan, eff_sizes)

kable(
  df %>%
    group_by(day) %>%
    summarise(
      diners = n(),
      meat_fish = mean(meat_fish),
      mn_weight = mean(food_weight),
      sd_weight = sd(food_weight),
      min_weight = min(food_weight),
      max_weight = max(food_weight),
      mn_co2e = mean(co2e),
      sd_co2e = sd(co2e),
      min_co2e = min(co2e),
      max_co2e = max(co2e)
    )
  ) %>% kable_styling()
day diners meat_fish mn_weight sd_weight min_weight max_weight mn_co2e sd_co2e min_co2e max_co2e
1 2787 0.5367779 481.1931 116.76963 121.0756 880.4076 1135.0653 943.0112 42.34491 2110.538
2 3328 0.4468149 526.6408 79.00456 211.7414 829.8962 1492.2975 1440.7114 90.71943 3208.592
3 3277 0.3918218 433.8547 91.53410 154.1031 757.1058 1221.1452 1010.7723 189.90937 2783.468
4 3147 0.4696536 447.4115 87.77478 161.2280 765.0852 1465.3187 1269.1851 98.56783 2981.574
5 2493 0.5371039 493.8781 91.06327 141.2710 845.3660 2474.9085 2146.9026 58.52325 4583.175
6 2891 0.4019370 466.8314 84.84336 152.1492 825.0216 998.8698 376.9882 231.88558 1972.395
7 3331 0.6871810 496.7692 118.94360 141.1710 896.8119 1794.1563 625.5598 266.64426 2939.362
8 3428 0.4859977 464.9775 92.97147 184.0874 817.1861 3433.1380 2313.7239 475.02012 6703.977
9 3372 0.7674970 387.6604 98.17851 137.7968 751.2241 755.2638 218.2122 147.48109 1223.321
10 2456 0.5158795 490.7109 90.18368 186.1105 824.0231 2419.5999 974.2702 649.70079 4158.067 
df <- df %>%
  mutate(day = as.factor(day), slot = as.factor(slot))

df_viz <- df %>% 
  mutate(
    day_tment = factor(
      (as.numeric(day) - 1) * 4 + as.numeric(slot),
      levels = 1:39
    )
  )

ylabs <- NULL
for (i in 1:parms$days) {
  tments <- as.character(tplan$info_cond)
  tments[tments == "NoNumericInfo"] <- "No Numeric Info"
  tments[tments == "NoContext"] <- "No Context"
  row_start <- (i - 1)*parms$treatment_slots_per_day + 1
  row_end <- i*parms$treatment_slots_per_day
  ylabs <- c(ylabs, tments[row_start:row_end])
  if (i < parms$days) ylabs <- c(ylabs, "")
}

fig <- ggplot(df_viz, aes(x = food_weight, color = color_cond, y = day_tment)) + 
  geom_density_ridges(fill = NA) + 
  labs(y = "", x = "Food purchased [grams]") +
  scale_y_discrete(labels = ylabs, drop = FALSE) + 
  scale_x_continuous(breaks = c(250, 500, 750))

style_ridges_plot(fig)

save_plot("../output/figure_2_food_weight", 8, 6)


df_viz2 <- df %>% 
  group_by(day, tment, info_cond, color_cond) %>%
  summarise(
    meat_fish = mean(meat_fish),
    .groups = "drop"
  ) %>%
  mutate(
    label = case_when(
      info_cond == "NoNumericInfo" ~  fontawesome("fa-question"),
      info_cond == "NoContext" ~ fontawesome("fa-cloud"),
      info_cond == "Budget" ~ fontawesome("fa-pie-chart"),
      info_cond == "Car" ~ fontawesome("fa-car"),
      info_cond == "Money" ~ fontawesome("fa-eur")
    ),
    nudge_left = case_when(
      day == 1 & info_cond == "NoContext" ~ TRUE,
      day == 2 & info_cond == "NoContext" ~ TRUE,
      day == 6 & meat_fish < 0.4 & info_cond == "Car" ~ TRUE,
      day == 10 & info_cond == "Money" & color_cond == "ColorCoded" ~ TRUE,
      TRUE ~ FALSE
    ),
    nudge_right = case_when(
      day == 1 & info_cond == "NoNumericInfo" & color_cond == "ColorCoded" ~ TRUE,
      day == 2 & info_cond == "Budget" & color_cond == "Neutral" ~ TRUE,
      day == 6 & meat_fish < 0.4 & info_cond == "Budget" ~ TRUE,
      day == 10 & info_cond == "Money" & color_cond == "Neutral"  ~ TRUE,
      TRUE ~ FALSE
    )
  )

ggplot() + 
  geom_text(
    data = df_viz2 %>% filter(!nudge_left & !nudge_right),
    aes(x = day, y = meat_fish, color = color_cond, label = label), 
    family='fontawesome-webfont', size = 6
  ) +
  geom_text(
    data = df_viz2 %>% filter(nudge_left),
    aes(x = day, y = meat_fish, color = color_cond, label = label), 
    family='fontawesome-webfont', size = 6,
    position = position_nudge(x = -0.15)
  ) +
  geom_text(
    data = df_viz2 %>% filter(nudge_right),
    aes(x = day, y = meat_fish, color = color_cond, label = label), 
    family='fontawesome-webfont', size = 6,
    position = position_nudge(x = +0.15)
  ) +
  scale_color_manual(
    values = c("ColorCoded" = "limegreen", "Neutral" = "black"),
    guide = FALSE
  ) +
  scale_y_continuous(labels = scales::label_percent(accuracy = 1L)) + 
  labs(x = "Day", y = "Share of meat or fish dishes") +
  theme_minimal()

save_plot("../output/figure_3_food_choice", 6, 6)

fig <- ggplot(df_viz, aes(x = co2e/1000, color = color_cond, y = day_tment)) + 
  geom_density_ridges(fill = NA) + 
  labs(y = "", x = "Carbon footprint [kgCO\u2082e]") +
  scale_y_discrete(labels = ylabs, drop = FALSE) + 
  scale_x_continuous(breaks = c(0, 2.5, 5.0, 7.5))

style_ridges_plot(fig)

save_plot("../output/figure_4_food_co2e", 8, 6)


weight_mod <- feols(
  log(food_weight) ~ info_cond*color_cond| day + slot, 
  df, cluster = c("day", "slot")
)
mtfsh_mod <- femlm(
  meat_fish ~ info_cond*color_cond| day + slot, 
  family = "logit", df, cluster = c("day", "slot")
)
mtfsh_mod_ols <- feols(
  meat_fish ~ info_cond*color_cond| day + slot, 
  df, cluster = c("day", "slot")
)
co2e_mod <- feols(
  log(co2e) ~ info_cond*color_cond| day + slot, 
  df, cluster = c("day", "slot")
)


modelsummary(
  list(
    weight = weight_mod, meat_fish_logit = mtfsh_mod, co2e = co2e_mod
  ),
  stars = TRUE
)
weight meat_fish_logit co2e
info_condNoContext -0.029*** -0.063 -0.044
(0.003) (0.043) (0.032)
info_condBudget -0.085*** -0.355*** -0.199***
(0.004) (0.050) (0.034)
info_condCar -0.066*** -0.339*** -0.167***
(0.006) (0.040) (0.048)
info_condMoney -0.042*** -0.294*** -0.133**
(0.003) (0.051) (0.057)
color_condColorCoded -0.042*** -0.190*** -0.116***
(0.003) (0.029) (0.016)
info_condNoContext × color_condColorCoded -0.010 -0.074 0.011
(0.007) (0.076) (0.035)
info_condBudget × color_condColorCoded -0.014 -0.045 -0.029
(0.005) (0.090) (0.052)
info_condCar × color_condColorCoded -0.012 -0.060*** -0.006
(0.007) (0.014) (0.015)
info_condMoney × color_condColorCoded -0.010** 0.045 0.020
(0.002) (0.061) (0.028)
Num.Obs. 30510 30510 30510
R2 0.158 0.205
R2 Adj. 0.157 0.205
R2 Within 0.021 0.006
R2 Pseudo 0.045
AIC -6346.7 40352.4 86263.5
BIC -6171.8 40527.3 86438.4
Log.Lik. 3194.329 -20155.219 -43110.770
FE: day X X X
FE: slot X X X
Std. errors Two-way (day & slot) Two-way (day & slot) Two-way (day & slot)
* p < 0.1, ** p < 0.05, *** p < 0.01 
test_hypothesis <- function(mod, test) {
  rv <- car::linearHypothesis(mod, test)
  tibble(
    effect = unname(attr(rv, "value")[1,1]),
    pvalue = rv[2, 4],
    starred = sprintf(
      "%.3f%s", effect, 
      case_when(
        pvalue < 0.01 ~ "***",
        pvalue < 0.05 ~ "**",
        pvalue < 0.1 ~ "*",
        TRUE ~ ""
      )
    )
  )
}

construct_lin_hypothesis <- function(df) {
  no_info <- info_cond[1]
  info_weights <- pull(
    df %>% filter(info_cond != no_info) %>% group_by(info_cond) %>% 
      summarise(share = n()/nrow(df))
  ) 
  
  pres_weight <- sum(df$color_cond == color_cond[1])/nrow(df)
  
  h1_test <- paste0(
    paste0(
      info_weights/sum(info_weights), " * info_cond", info_cond[-1], collapse = " + "
    ), " + ", 
    paste0(
      info_weights/sum(info_weights) * pres_weight, " * ", 
      "info_cond", info_cond[-1], ":color_cond", color_cond[2], collapse = " + "
    ),
    " = 0"
  )
  
  h2_test <- paste0(
    "- info_cond", info_cond[2], " ", -pres_weight, 
    " * info_cond", info_cond[2], ":color_cond", color_cond[2], " + ", 
    paste0(
      info_weights[2:4]/sum(info_weights[2:4]), " * info_cond", info_cond[3:5], 
      collapse = " + "
    ), 
    " + ",
    paste0(
      info_weights[2:4]/sum(info_weights[2:4]) * pres_weight,  
      " * info_cond", info_cond[3:5], ":color_cond", color_cond[2], collapse = " + "
    ),
    " = 0"
  )
  
  h3a_test <- paste0(
    "info_cond", info_cond[3], " + ", pres_weight, 
    " * info_cond", info_cond[3], ":color_cond", color_cond[2], " - ", 
    "info_cond", info_cond[4], " - ", pres_weight,  
    " * info_cond", info_cond[4], ":color_cond", color_cond[2], 
    " = 0"
  )
  
  h3b_test <- paste0(
    "info_cond", info_cond[4], " + ", pres_weight, 
    " * info_cond", info_cond[4], ":color_cond", color_cond[2], " - ", 
    "info_cond", info_cond[5], " - ", pres_weight,  
    " * info_cond", info_cond[5], ":color_cond", color_cond[2], 
    " = 0"
  )
  
  h4_test <- paste0(
    "color_cond", color_cond[2], " + ", 
    paste0(
      paste0(pull(
        df %>% filter(info_cond != no_info) %>% group_by(info_cond) %>% 
          summarise(share = n()/nrow(df))
      ) , " * "), 
      "info_cond", info_cond[-1], ":color_cond", color_cond[2], collapse = " + "
    ),
    " = 0"
  )
  return(c(
    h1 = h1_test,
    h2 = h2_test,
    h3a = h3a_test,
    h3b = h3b_test,
    h4 = h4_test
  ))
} 

lh <- construct_lin_hypothesis(df)

hresults <- matrix(NA_character_, 5, 3)

col <- 0
for (mod in list(weight_mod, mtfsh_mod, co2e_mod)) {
  col <- col + 1
  hresults[1, col] <- test_hypothesis(mod, lh["h1"])$starred
  hresults[2, col] <- test_hypothesis(mod, lh["h2"])$starred
  hresults[3, col] <- test_hypothesis(mod, lh["h3a"])$starred
  hresults[4, col] <- test_hypothesis(mod, lh["h3b"])$starred
  hresults[5, col] <- test_hypothesis(mod, lh["h4"])$starred
}
rownames(hresults) <- c("H1", "H2", "H3a", "H3b", "H4")
colnames(hresults) <- c("weight_mod", "mtfsh_mod", "co2e_mod")

kable(hresults) %>% kable_styling()
weight_mod mtfsh_mod co2e_mod
H1 -0.062*** -0.282** -0.137***
H2 -0.036*** -0.240** -0.132***
H3a -0.021** -0.009 -0.044
H3b -0.025** -0.098** -0.046*
H4 -0.052*** -0.217** -0.117***

Calculating Power

Now we estimate the power for our hypotheses tests by

test_hypotheses <- function(df, level) {
  weight_mod <- feols(
    log(food_weight) ~ info_cond*color_cond| day + slot, 
    cluster = c("day", "slot"), df
  )
  mtfsh_mod <- femlm(
    meat_fish ~ info_cond*color_cond| day + slot, family = "logit", 
    cluster = c("day", "slot"), df
  )
  co2e_mod <- feols(
    log(co2e) ~ info_cond*color_cond| day + slot, 
    cluster = c("day", "slot"), df
  )

  bind_cols(
    expand_grid(
      mod = c("weight", "mtfsh", "co2e"),
      test = c("h1", "h2", "h3a", "h3b", "h4")
    ),
    bind_rows(lapply(
    list(weight_mod, mtfsh_mod, co2e_mod),
    function(mod) {
      bind_rows(lapply(
        construct_lin_hypothesis(df),
        function(test) test_hypothesis(mod, test)
      ))
    }
    ))
  )
}

run_sim <- function(p, tp, te) {
  mp <- sim_meal_plan()
  bs <- sim_base_sample(p, mp)
  df <- apply_teffects(bs, mp, tp, te)
  test_hypotheses(df)
}

if (file.exists("../data/generated/power_runs.rds")) {
  df <- readRDS("../data/generated/power_runs.rds")
} else {
  df <- do.call(
    rbind, 
    replicate(1000, run_sim(parms, tplan, eff_sizes), simplify = FALSE)
  )
  saveRDS(df, "../data/generated/power_runs.rds")
}

sim_results <-  df %>% group_by(mod, test) %>%
  filter(!is.na(pvalue) & pvalue != 1) %>%
  summarise(
    n = n(),
    mn_effect = mean(effect),
    sd_effect = sd(effect),
    pct_sig_neg = mean(pvalue < 0.05 & effect < 0),
    pct_sig_pos = mean(pvalue < 0.05 & effect > 0)
  ) %>%
  pivot_wider(
    id_cols = test, names_from = mod,
    values_from = c(n, mn_effect, sd_effect, pct_sig_pos, pct_sig_neg),
    names_glue = "{mod}_{.value}"
  ) %>%
  select(test, starts_with("weight"), starts_with("mtfsh"), starts_with("co2e"))
  
kable(sim_results) %>% kable_styling()
test weight_n weight_mn_effect weight_sd_effect weight_pct_sig_pos weight_pct_sig_neg mtfsh_n mtfsh_mn_effect mtfsh_sd_effect mtfsh_pct_sig_pos mtfsh_pct_sig_neg co2e_n co2e_mn_effect co2e_sd_effect co2e_pct_sig_pos co2e_pct_sig_neg
h1 991 -0.0532858 0.0044594 0 1.0000000 987 -0.2148750 0.0440937 0.000000 0.9118541 990 -0.1096010 0.0244533 0.000000 0.9030303
h2 989 -0.0431221 0.0045963 0 1.0000000 988 -0.1677387 0.0443400 0.000000 0.8248988 989 -0.0876007 0.0214355 0.000000 0.8685541
h3a 978 -0.0221133 0.0059461 0 0.7934560 976 -0.0874515 0.0564016 0.000000 0.3032787 980 -0.0467585 0.0279185 0.000000 0.2887755
h3b 995 -0.0216296 0.0051758 0 0.7959799 997 -0.0826989 0.0513707 0.001003 0.2808425 998 -0.0428989 0.0274234 0.001002 0.2294589
h4 1000 -0.0534723 0.0026821 0 1.0000000 1000 -0.2104588 0.0257848 0.000000 1.0000000 1000 -0.1089842 0.0149074 0.000000 0.9990000

Power by Effect Size

How does the power of a single treatment vary relative to its (differential) effect size?

test_diff <- function(p, tp, te) {
  mp <- sim_meal_plan()
  bs <- sim_base_sample(p, mp)
  df <- apply_teffects(bs, mp, tp, te)
 
  weight_mod <- feols(
    log(food_weight) ~ info_cond*color_cond| day + slot, 
    cluster = c("day", "slot"), df
  )
  mtfsh_mod <- femlm(
    meat_fish ~ info_cond*color_cond| day + slot, family = "logit", 
    cluster = c("day", "slot"), df
  )
  co2e_mod <- feols(
    log(co2e) ~ info_cond*color_cond| day + slot, 
    cluster = c("day", "slot"), df
  )

  
  pres_weight <- sum(df$color_cond == "ColorCoded")/nrow(df)
  
  diff_test <- paste0(
    "info_condNoContext + ", pres_weight, " * info_condNoContext:color_condColorCoded = 0"
  )
  
  bind_cols(
    mod = c("weight", "mtfsh", "co2e"),
    bind_rows(lapply(
      list(weight_mod, mtfsh_mod, co2e_mod),
      function(mod) {
        test_hypothesis(mod, diff_test)
      }
    ))
  )   
}

if (file.exists("../data/generated/power_differential_effects.rds")) {
  df <- readRDS("../data/generated/power_differential_effects.rds")
} else {
  prot <- expand_grid(
    run = 1:1000,
    diff_effect = seq(-0.0, -0.10, length.out = 21)
  )
  
  run_protocol <- function(n) {
    df <- prot[n, ]
    te <- eff_sizes
    te[te$tment == "Neutral_NoContext", "mn_weight"] <- df$diff_effect
    te[te$tment == "ColorCoded_NoContext", "mn_weight"] <- df$diff_effect + color_eff 
    te[te$tment == "Neutral_NoContext", "mn_mtfsh"] <- df$diff_effect
    te[te$tment == "ColorCoded_NoContext", "mn_mtfsh"] <- df$diff_effect + color_eff
    message(sprintf(
      "Running run %d with parameter %.3f", df$run, df$diff_effect
    ))
    bind_cols(df, test_diff(parms, tplan, te))
  }
  
  df <- bind_rows(
    lapply(1:nrow(prot), run_protocol)
  )
  
  saveRDS(df, "../data/generated/power_differential_effects.rds")
}

df %>%
  group_by(mod, diff_effect) %>%
  summarise(
    alpha = mean(effect < 0 & pvalue < 0.05),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = diff_effect, y = alpha, color = mod, shape = mod)) + 
  geom_point() +
  geom_line() +
  labs(y = "Power", x = "Treatment effect differential [percentage points]")+
  scale_y_continuous(breaks = seq(0, 1, by = 0.2)) +
    geom_hline(yintercept = 0.8, color = col_trr266_red, lty = 2) +
    scale_shape_manual(
      values = c("co2e" = 16, "mtfsh" = 15, "weight" = 17),
      name = "Dependent variables",
      labels = c("CO\u2082 emission", "Choice for meat or fish", "Food weight")
      ) +
    scale_color_manual(
      values = c("co2e" = col_trr266_petrol, "mtfsh" = col_trr266_yellow, "weight" = col_trr266_blue),
      name = "Dependent variables",
      labels = c("CO\u2082 emission", "Choice for meat or fish", "Food weight")
      )+
    theme_minimal() +
    theme(
      legend.position = "bottom"
    )

save_plot("../output/figure_5_power", 6, 6)