# install.packages(c("stats", "emmeans", "afex", "ggplot2", "tidyverse", "dplyr"))
library(stats)
library(emmeans)
library(afex) # for aov_ev (RM)
library(ggplot2)
library(tidyverse)
library(dplyr)

# the anonymized dataset containing only the variables used in this analysis is data_n1606_anon

# Bicycle Support ---------------------------------------------------------
# select all support types for mode bicycle
bicycle_support <- data_n1606_anon %>% select(c(ResponseID, 
                                                          GeneralSupportTransportBicycle,
                                                          SupportSafetyBicycle,
                                                          SupportSpaceBicycle,
                                                          SupportFinancialBicycle,
                                                          SupportEquityBicycle,
                                                          SupportInfraDevBicycle))
# pivot longer
bicycle_long <- bicycle_support %>%
  select(ResponseID, starts_with("Support"), GeneralSupportTransportBicycle) %>%
  pivot_longer(
    cols = c(GeneralSupportTransportBicycle, starts_with("Support")),
    names_to = "SupportType",
    values_to = "SupportValue"
  ) %>%
  mutate(
    SupportType = recode(SupportType,
                         "GeneralSupportTransportBicycle" = "General",
                         "SupportSafetyBicycle" = "Safety",
                         "SupportSpaceBicycle" = "Space",
                         "SupportFinancialBicycle" = "Financial",
                         "SupportEquityBicycle" = "Equity",
                         "SupportInfraDevBicycle" = "InfraDev"
    ),
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial","InfraDev","Safety","Space"))
  )
# repeated measure ANOVA
aov_bicycle_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "SupportType",
  data   = bicycle_long,
  anova_table = list(correction = "GG", es = "ges")
)
aov_bicycle_rm
# posthoc testing
em <- emmeans(aov_bicycle_rm, ~ SupportType)
posthoc <- pairs(em, adjust = "holm")
posthoc

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(bicycle_long$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# Car Support -----------------------------------------------------------
# select all support types for mode car
car_support <- data_n1606_anon %>% select(c(ResponseID, 
                                                      GeneralSupportTransportCar,
                                                      SupportSafetyCar,
                                                      SupportSpaceCar,
                                                      SupportFinancialCar,
                                                      SupportEquityCar,
                                                      SupportInfraDevCar))
# pivot longer
car_long <- car_support %>%
  select(ResponseID, starts_with("Support"), GeneralSupportTransportCar) %>%
  pivot_longer(
    cols = c(GeneralSupportTransportCar, starts_with("Support")),
    names_to = "SupportType",
    values_to = "SupportValue"
  ) %>%
  mutate(
    SupportType = recode(SupportType,
                         "GeneralSupportTransportCar" = "General",
                         "SupportSafetyCar" = "Safety",
                         "SupportSpaceCar" = "Space",
                         "SupportFinancialCar" = "Financial",
                         "SupportEquityCar" = "Equity",
                         "SupportInfraDevCar" = "InfraDev"
    ),
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial","InfraDev","Safety","Space"))
  )

# run repeated measure ANOVA
aov_car_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "SupportType",
  data   = car_long,
  anova_table = list(correction = "GG", es = "ges")
)
aov_car_rm

# posthoc testing 
em <- emmeans(aov_car_rm, ~ SupportType)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_car_rm$anova_table # ges = 0.0066387

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(car_long$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# Walking Support ---------------------------------------------------------
# select all support types for mode walking
walking_support <- data_n1606_anon %>% select(c(ResponseID, 
                                                          GeneralSupportTransportWalking,
                                                          SupportSafetyWalking,
                                                          SupportSpaceWalking,
                                                          SupportFinancialWalking,
                                                          SupportEquityWalking,
                                                          SupportInfraDevWalking))
# pivot longer
walking_long <- walking_support %>%
  select(ResponseID, starts_with("Support"), GeneralSupportTransportWalking) %>%
  pivot_longer(
    cols = c(GeneralSupportTransportWalking, starts_with("Support")),
    names_to = "SupportType",
    values_to = "SupportValue"
  ) %>%
  mutate(
    SupportType = recode(SupportType,
                         "GeneralSupportTransportWalking" = "General",
                         "SupportSafetyWalking" = "Safety",
                         "SupportSpaceWalking" = "Space",
                         "SupportFinancialWalking" = "Financial",
                         "SupportEquityWalking" = "Equity",
                         "SupportInfraDevWalking" = "InfraDev"
    ),
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial","InfraDev","Safety","Space"))
  )
# repeated measure ANOVA
aov_walking_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "SupportType",
  data   = walking_long,
  anova_table = list(correction = "GG", es = "ges")
)
aov_walking_rm
# posthoc testing
em <- emmeans(aov_walking_rm, ~ SupportType)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_walking_rm$anova_table # ges = 0.096967

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(walking_long$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# Bus Support ---------------------------------------------------------
# select all support types for mode bus
bus_support <- data_n1606_anon %>% select(c(ResponseID, 
                                                      GeneralSupportTransportBus,
                                                      SupportSafetyBus,
                                                      SupportSpaceBus,
                                                      SupportFinancialBus,
                                                      SupportEquityBus,
                                                      SupportInfraDevBus))
# pivot longer
bus_long <- bus_support %>%
  select(ResponseID, starts_with("Support"), GeneralSupportTransportBus) %>%
  pivot_longer(
    cols = c(GeneralSupportTransportBus, starts_with("Support")),
    names_to = "SupportType",
    values_to = "SupportValue"
  ) %>%
  mutate(
    SupportType = recode(SupportType,
                         "GeneralSupportTransportBus" = "General",
                         "SupportSafetyBus" = "Safety",
                         "SupportSpaceBus" = "Space",
                         "SupportFinancialBus" = "Financial",
                         "SupportEquityBus" = "Equity",
                         "SupportInfraDevBus" = "InfraDev"
    ),
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial","InfraDev","Safety","Space"))
  )
# repeated measure ANOVA
aov_bus_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "SupportType",
  data   = bus_long,
  anova_table = list(correction = "GG", es = "ges")
)
aov_bus_rm
# posthoc testing
em <- emmeans(aov_bus_rm, ~ SupportType)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_bus_rm$anova_table # ges = 0.036195

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(bus_long$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# Rail Support ---------------------------------------------------------
# select all support types for mode rail
rail_support <- data_n1606_anon %>% select(c(ResponseID, 
                                                       GeneralSupportTransportRail,
                                                       SupportSafetyRail,
                                                       SupportSpaceRail,
                                                       SupportFinancialRail,
                                                       SupportEquityRail,
                                                       SupportInfraDevRail))
# pivot longer
rail_long <- rail_support %>%
  select(ResponseID, starts_with("Support"), GeneralSupportTransportRail) %>%
  pivot_longer(
    cols = c(GeneralSupportTransportRail, starts_with("Support")),
    names_to = "SupportType",
    values_to = "SupportValue"
  ) %>%
  mutate(
    SupportType = recode(SupportType,
                         "GeneralSupportTransportRail" = "General",
                         "SupportSafetyRail" = "Safety",
                         "SupportSpaceRail" = "Space",
                         "SupportFinancialRail" = "Financial",
                         "SupportEquityRail" = "Equity",
                         "SupportInfraDevRail" = "InfraDev"
    ),
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial","InfraDev","Safety","Space"))
  )
# repeated measure ANOVA
aov_rail_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "SupportType",
  data   = rail_long,
  anova_table = list(correction = "GG", es = "ges")
)
aov_rail_rm
# posthoc testing
em <- emmeans(aov_rail_rm, ~ SupportType)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_rail_rm$anova_table # ges = 0.016515

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(rail_long$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)


# VISUALIZE TYPES OF SUPPORT ACROSS MODES ------------------------------------------------------


# make df with just the support variables for each mode
modal_support <- data_n1606_anon %>% select(c(
  ResponseID,
  GeneralSupportTransportBus:GeneralSupportTransportCar,
  SupportSafetyBus:SupportInfraDevCar
))

# pivot longer
modal_long <- modal_support %>%
  pivot_longer(
    cols = -ResponseID,
    names_to = c("SupportType", "Mode"),
    names_pattern = "(.*)(Bus|Walking|Rail|Bicycle|Car)$",
    values_to = "SupportValue"
  )

# pivot long & make changes for cleaner plotting
modal_long2 <- modal_support %>%
  pivot_longer(
    cols = -ResponseID,
    names_to = c("SupportType", "Mode"),
    names_pattern = "(.*)(Car|Bus|Rail|Walking|Bicycle)$",
    values_to = "SupportValue"
  ) %>%
  mutate(
    # clean labels for plotting
    SupportType = recode(SupportType,
                         "GeneralSupportTransport" = "General",
                         "SupportSafety"          = "Safety",
                         "SupportSpace"           = "Space",
                         "SupportFinancial"       = "Financial",
                         "SupportEquity"          = "Equity",
                         "SupportInfraDev"        = "Infrastructure\nDevelopment",
                         .default = SupportType
    ),
    Mode = factor(Mode,
                  ordered = TRUE,
                  levels = c("Car", "Bicycle", "Walking", "Bus", "Rail")),
    # set a readable order for the x-axis
    SupportType = factor(SupportType,
                         levels = c("General","Equity","Financial", "Infrastructure\nDevelopment", "Safety","Space"))
  )
# plot boxplots
ggplot(modal_long2, aes(x = SupportType, y = SupportValue)) +
  geom_boxplot() +
  facet_wrap(~ Mode, scales = "free_y") +
  stat_summary(fun.y=mean, geom="point", shape=20, size=2, color="red", fill="red") +
  theme_bw() +
  labs(
    title = "Support Types across Modes",
    x = "Support Type",
    y = "Support Value"
  ) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))



# GENERAL -----------------------------------------------------------------
# select all modes for general support
general_data <- modal_long %>% filter(SupportType == "GeneralSupportTransport")
# run RM ANOVA
aov_general_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = general_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_general_rm
# posthoc testing
em <- emmeans(aov_general_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_general_rm$anova_table # ges = 0.15254

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(general_data$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)


# SPACE -------------------------------------------------------------------
# select all modes for space support
space_data <- modal_long %>% filter(SupportType == "SupportSpace")
# RM ANOVA
aov_space_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = space_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_space_rm
# posthoc testing
em <- emmeans(aov_space_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_space_rm$anova_table # ges = 0.10135

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(space_data$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)


# SAFETY ------------------------------------------------------------------
# select all modes for safety support
safety_data <- modal_long %>% filter(SupportType == "SupportSafety")
# RM ANOVA
aov_safety_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = safety_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_safety_rm
# posthoc testing
em <- emmeans(aov_safety_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_safety_rm$anova_table # ges = 0.0909

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(safety_data$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)


# FINANCIAL ---------------------------------------------------------------
# select all modes for financial support
financial_data <- modal_long %>% filter(SupportType == "SupportFinancial")
# RM ANOVA
aov_financial_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = financial_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_financial_rm
# posthoc testing
em <- emmeans(aov_financial_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_financial_rm$anova_table # ges = 0.12264

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(financial_data$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)


# EQUITY ------------------------------------------------------------------
# select all modes for equity support
equity_data <- modal_long %>% filter(SupportType == "SupportEquity")
# RM ANOVA
aov_equity_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = equity_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_equity_rm
# posthoc testing
em <- emmeans(aov_equity_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_equity_rm$anova_table # ges = 0.19687

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(equity_data$SupportValue) 
posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# INFRASTRUCTURE DEVELOPMENT ----------------------------------------------
# select all modes for infra support
infra_data <- modal_long %>% filter(SupportType == "SupportInfraDev")
# RM ANOVA
aov_infra_rm <- aov_ez(
  id     = "ResponseID",
  dv     = "SupportValue",
  within = "Mode",
  data   = infra_data,
  anova_table = list(correction = "GG", es = "ges")
)
aov_infra_rm
# posthoc testing
em <- emmeans(aov_infra_rm, ~ Mode)
posthoc <- pairs(em, adjust = "holm")
posthoc

# check effect size (provided in anova_table)
aov_infra_rm$anova_table # ges = 0.10922

# cohen's d
posthoc_summary <- summary(posthoc)

pooled_sd <- sd(infra_data$SupportValue)  

posthoc_summary %>%
  mutate(cohens_d = estimate / pooled_sd)

# VISUALIZE MODAL SUPPORT ACROSS TYPES OF SUPPORT ---------------------------------------------------------------

# order modes in ascending order (generally) for visualization
modal_long$Mode <-  factor(modal_long$Mode,
                           ordered = TRUE,
                           levels = c("Car", "Bicycle", "Walking", "Bus", "Rail"))
# clean up support type labels for plotting
modal_long <- modal_long %>% mutate(SupportType = recode(SupportType,
                                                         "GeneralSupportTransport" = "General",
                                                         "SupportSafety"          = "Safety",
                                                         "SupportSpace"           = "Space",
                                                         "SupportFinancial"       = "Financial",
                                                         "SupportEquity"          = "Equity",
                                                         "SupportInfraDev"        = "Infrastructure\nDevelopment",
                                                         .default = SupportType
),
SupportType = factor(SupportType, ordered = T,
                     levels = c("General","Equity","Financial", "Infrastructure\nDevelopment", "Safety","Space")))


# visualize
ggplot(modal_long, aes(x = Mode, y = SupportValue)) +
  geom_boxplot() +
  facet_wrap(~ SupportType, scales = "free_y") +
  stat_summary(fun.y=mean, geom="point", shape=20, size=2, color="red", fill="red") +
  theme_bw() +
  labs(title = "Support by Mode across Support Types",
       y = "Support Value")