# ---
#   output:
#   pdf_document: default
# html_document: default
# ---
  
  
  
  # ANALYSIS
  # Park accessibility and park use behavior
  
  rm(list = ls())

# Packages

# ONCE

# install.packages(c(
#   "readr",
#   "dplyr",
#   "tidyr",
#   "ggplot2",
#   "scales",
#   "clustMixType",
#   "patchwork",
#   "tibble",
#   "rstudioapi",
# ))
# 
  
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(clustMixType)
library(patchwork)

select <- dplyr::select
filter <- dplyr::filter
mutate <- dplyr::mutate
rename <- dplyr::rename
recode <- dplyr::recode

# Data import

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

data <- read_csv("park_accessibility_dataset.csv", show_col_types = FALSE)

# Variable checks

required_vars <- c(
  "id_1",
  "distance_min", "distance_max", "nb_parcs_visites",
  "PARKFREQ", "PARKHAB.ME.", "PARKHAB.PET.", "PARKHAB.GRP.", "PARKHAB.CHILD.",
  "MODEPARK", "TYPE",
  "CIPQAY", "ACCOBJ_CP", "ACCAREA_CP",
  "SEX", "AGE", "ETHNICITE1", "EDUC", "REVENU1", "CHILD1", "AUTO1",
  "PARKIMPT", "PERNUMB", "PEREVA", "PERSATIS", "PERACC", "Grand_parc", "Parc_arrondissement"
)

missing_vars <- setdiff(required_vars, names(data))

if (length(missing_vars) > 0) {
  stop(
    "Variables missing from park_accessibility_dataset.csv: ",
    paste(missing_vars, collapse = ", ")
  )
}

# Recoding

recode_frequency <- function(x) {
  case_when(
    x == "Tous les jours ou presque" ~ 5,
    x == "Quelques fois par semaine" ~ 4,
    x == "Quelques fois par mois" ~ 3,
    x == "Quelques fois durant les trois mois" ~ 2,
    x == "Jamais ou presque jamais" ~ 1,
    TRUE ~ NA_real_
  )
}

data <- data %>%
  mutate(
    PARKFREQ = recode_frequency(PARKFREQ),
    PARKHAB.ME. = recode_frequency(PARKHAB.ME.),
    PARKHAB.PET. = recode_frequency(PARKHAB.PET.),
    PARKHAB.GRP. = recode_frequency(PARKHAB.GRP.),
    PARKHAB.CHILD. = recode_frequency(PARKHAB.CHILD.)
  )

# Park visit context

data <- data %>%
  rowwise() %>%
  mutate(
    PARKHAB.MOYENNE = mean(
      c_across(c(PARKHAB.ME., PARKHAB.PET., PARKHAB.GRP., PARKHAB.CHILD.)),
      na.rm = TRUE
    )
  ) %>%
  ungroup() %>%
  mutate(
    PARKHAB.ECART_ME = PARKHAB.ME. - PARKHAB.MOYENNE,
    PARKHAB.ECART_PET = PARKHAB.PET. - PARKHAB.MOYENNE,
    PARKHAB.ECART_GRP = PARKHAB.GRP. - PARKHAB.MOYENNE,
    PARKHAB.ECART_CHILD = PARKHAB.CHILD. - PARKHAB.MOYENNE
  )

context_data <- data %>%
  select(
    id_1,
    PARKHAB.ECART_ME,
    PARKHAB.ECART_PET,
    PARKHAB.ECART_GRP,
    PARKHAB.ECART_CHILD
  )

set.seed(123)

context_model <- kmeans(
  context_data %>% select(-id_1),
  centers = 6,
  nstart = 25
)

context_data$context_cluster <- context_model$cluster

data <- data %>%
  left_join(context_data %>% select(id_1, context_cluster), by = "id_1") %>%
  mutate(
    CONTEXTE = case_when(
      context_cluster == 1 ~ "Seul",
      context_cluster == 2 ~ "Avec animal et sans enfants",
      context_cluster == 3 ~ "Aucun comportement prédominant",
      context_cluster == 4 ~ "En groupe",
      context_cluster == 5 ~ "Avec enfants",
      context_cluster == 6 ~ "Seul ou en groupe et sans animal",
      TRUE ~ NA_character_
    )
  ) %>%
  select(-context_cluster)

# TABLE 2 - Summary statistics for clustering variables

# Numerical variables

num_stats <- data %>%
  summarise(
    across(
      c(distance_min, distance_max, PARKFREQ, nb_parcs_visites),
      list(
        min = ~ min(.x, na.rm = TRUE),
        median = ~ median(.x, na.rm = TRUE),
        max = ~ max(.x, na.rm = TRUE),
        mean = ~ mean(.x, na.rm = TRUE),
        sd = ~ sd(.x, na.rm = TRUE)
      )
    )
  )

# Categorical variables

format_n_pct <- function(x, category) {
  n <- sum(x == category, na.rm = TRUE)
  pct <- round(n / nrow(data) * 100)
  paste0(n, " (", pct, ")")
}

transport <- c(
  Walking = format_n_pct(data$MODEPARK, "Marche"),
  Cycling = format_n_pct(data$MODEPARK, "Vélo"),
  `Public transport` = format_n_pct(data$MODEPARK, "Transport en commun"),
  Car = format_n_pct(data$MODEPARK, "Automobile")
)

behaviour <- c(
  Alone = format_n_pct(data$CONTEXTE, "Seul"),
  `Alone or in group without pet` =
    format_n_pct(data$CONTEXTE, "Seul ou en groupe et sans animal"),
  `In group` = format_n_pct(data$CONTEXTE, "En groupe"),
  `No predominant behavior` =
    format_n_pct(data$CONTEXTE, "Aucun comportement prédominant"),
  `With children` = format_n_pct(data$CONTEXTE, "Avec enfants"),
  `With pet and without children` =
    format_n_pct(data$CONTEXTE, "Avec animal et sans enfants")
)

park_type <- c(
  Both = format_n_pct(data$TYPE, "Les deux"),
  `Large park` = format_n_pct(data$TYPE, "Grand parc"),
  `Neighborhood park` = format_n_pct(data$TYPE, "Parc de quartier"),
  Neither = format_n_pct(data$TYPE, "Aucun des deux")
)

# Final table

table2 <- data.frame(
  Variable = c(
    "Continuous variables",
    "",
    "Minimum distance (m)",
    "Maximum distance (m)",
    "Park visit frequency",
    "Number of different parks visited",
    "Categorical variables - Number of observations (%)",
    "",
    "Main transport mode",
    "",
    "Visit behaviour",
    "",
    "Type of parks visited"
  ),
  
  C1 = c(
    "", "Min",
    round(num_stats$distance_min_min, 2),
    round(num_stats$distance_max_min, 2),
    round(num_stats$PARKFREQ_min, 2),
    round(num_stats$nb_parcs_visites_min, 2),
    "", "Walking", transport["Walking"],
    "Alone", behaviour["Alone"],
    "Both", park_type["Both"]
  ),
  
  C2 = c(
    "", "Median",
    round(num_stats$distance_min_median, 2),
    round(num_stats$distance_max_median, 2),
    round(num_stats$PARKFREQ_median, 2),
    round(num_stats$nb_parcs_visites_median, 2),
    "", "Cycling", transport["Cycling"],
    "Alone or in group without pet",
    behaviour["Alone or in group without pet"],
    "Large park", park_type["Large park"]
  ),
  
  C3 = c(
    "", "Max",
    round(num_stats$distance_min_max, 2),
    round(num_stats$distance_max_max, 2),
    round(num_stats$PARKFREQ_max, 2),
    round(num_stats$nb_parcs_visites_max, 2),
    "", "Public transport", transport["Public transport"],
    "In group", behaviour["In group"],
    "Neighborhood park", park_type["Neighborhood park"]
  ),
  
  C4 = c(
    "", "Mean",
    round(num_stats$distance_min_mean, 2),
    round(num_stats$distance_max_mean, 2),
    round(num_stats$PARKFREQ_mean, 2),
    round(num_stats$nb_parcs_visites_mean, 2),
    "", "Car", transport["Car"],
    "No predominant behavior",
    behaviour["No predominant behavior"],
    "Neither", park_type["Neither"]
  ),
  
  C5 = c(
    "", "Std dev.",
    round(num_stats$distance_min_sd, 2),
    round(num_stats$distance_max_sd, 2),
    round(num_stats$PARKFREQ_sd, 2),
    round(num_stats$nb_parcs_visites_sd, 2),
    "", "", "",
    "With children", behaviour["With children"],
    "", ""
  ),
  
  C6 = c(
    "", "", "", "", "", "",
    "", "", "",
    "With pet and without children",
    behaviour["With pet and without children"],
    "", ""
  )
)

write.csv(
  table2,
  "Table2_summary_statistics.csv",
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

# TABLE 3 - Descriptive characteristics

# Functions

n_pct <- function(x, value, total = nrow(data)) {
  n <- sum(x == value, na.rm = TRUE)
  paste0(n, " (", round(n / total * 100), ")")
}

stats <- function(x) {
  round(c(
    min(x, na.rm = TRUE),
    median(x, na.rm = TRUE),
    max(x, na.rm = TRUE),
    mean(x, na.rm = TRUE),
    sd(x, na.rm = TRUE)
  ), 2)
}

# Education categories

education <- case_when(
  data$EDUC == "Certificat ou diplôme collégial, cégep ou autre non universitaire" ~ "CEGEP",
  data$EDUC %in% c(
    "Diplôme d'études secondaires ou équivalent",
    "Moins que l'équivalent du secondaire"
  ) ~ "High school",
  data$EDUC == "Certificat ou diplôme d'apprenti enregistré ou autre métier" ~ "Professional program",
  data$EDUC %in% c(
    "Baccalauréat universitaire",
    "Maîtrise (p. ex., M.A., M.Sc., M.Ed., M.B.A.)",
    "Diplôme de doctorat (par exemple, Ph.D.)",
    "Diplôme en médecine, dentisterie, médecine vétérinaire ou optométrie (M.D., D.D.S., D.M.D., D.V.M., O.D.)"
  ) ~ "University",
  TRUE ~ "NA"
)

# Numerical statistics

age <- stats(data$AGE)
cipqay <- stats(data$CIPQAY)
parks <- stats(data$ACCOBJ_CP)
area <- stats(data$ACCAREA_CP)

# Final table

table3 <- data.frame(
  Variable = c(
    "Sociodemographic characteristics - Number of observations (%)",
    "", "Gender",
    "", "Age",
    "", "Number of children",
    "", "Number of cars",
    "", "Household income",
    "", "Education level",
    "", "Ethnic background",
    
    "Importance of parks and green spaces - Number of observations (%)",
    "",
    'Self-declared park importance: "Parks and green spaces are important to me since they connect me to nature"',
    
    "Calculated accessibility*",
    "",
    "Combined Index of Park Quality and Accessibility for Youth (CIPQAY)",
    "Number of accessible parks",
    "Accessible park area (km²)",
    
    "Perceived accessibility - Number of observations (%)",
    "",
    "Easy park access: I have access to parks and green spaces that meet my needs within 5 minutes of walking from my house.",
    "Perceived number of parks accessible: There is a high number of parks and green spaces in my neighborhood.",
    "Satisfied with park access: I am satisfied with my access to parks and green spaces.",
    "",
    "Overall evaluation of park access: How would you rate your overall access to parks and greenspaces from your home?"
  ),
  
  C1 = c(
    "", "Men", n_pct(data$SEX, "Homme"),
    "Min", age[1],
    "0", n_pct(data$CHILD1, "0"),
    "1 car", n_pct(data$AUTO1, "1 automobile"),
    "High", n_pct(data$REVENU1, "Elevé"),
    "CEGEP", n_pct(education, "CEGEP"),
    "Majority group", n_pct(data$ETHNICITE1, "Majorité", sum(!is.na(data$ETHNICITE1))),
    
    "", "Strongly disagree", n_pct(data$PARKIMPT, "Fortement en désaccord"),
    
    "", "Min",
    cipqay[1], parks[1], area[1],
    
    "", "Strongly disagree",
    n_pct(data$PERACC, "Fortement en désaccord"),
    n_pct(data$PERNUMB, "Fortement en désaccord"),
    n_pct(data$PERSATIS, "Fortement en désaccord"),
    "None", n_pct(data$PEREVA, "Nul")
  ),
  
  C2 = c(
    "", "Women", n_pct(data$SEX, "Femme"),
    "Median", age[2],
    "1", n_pct(data$CHILD1, "1"),
    "2+ cars", n_pct(data$AUTO1, "2 automobiles et plus"),
    "Low", n_pct(data$REVENU1, "Faible"),
    "High school", n_pct(education, "High school"),
    "Visible minority", n_pct(data$ETHNICITE1, "Minorité visible", sum(!is.na(data$ETHNICITE1))),
    
    "", "Disagree", n_pct(data$PARKIMPT, "En désaccord"),
    
    "", "Median",
    cipqay[2], parks[2], area[2],
    
    "", "Disagree",
    n_pct(data$PERACC, "En désaccord"),
    n_pct(data$PERNUMB, "En désaccord"),
    n_pct(data$PERSATIS, "En désaccord"),
    "Very low", n_pct(data$PEREVA, "Très faible")
  ),
  
  C3 = c(
    "", "", "",
    "Max", age[3],
    "2", n_pct(data$CHILD1, "2"),
    "None", n_pct(data$AUTO1, "Aucune"),
    "Medium", n_pct(data$REVENU1, "Moyen"),
    "NA", n_pct(education, "NA"),
    "", "",
    
    "", "Neutral", n_pct(data$PARKIMPT, "Neutre"),
    
    "", "Max",
    cipqay[3], parks[3], area[3],
    
    "", "Neutral",
    n_pct(data$PERACC, "Neutre"),
    n_pct(data$PERNUMB, "Neutre"),
    n_pct(data$PERSATIS, "Neutre"),
    "Low", n_pct(data$PEREVA, "Faible")
  ),
  
  C4 = c(
    "", "", "",
    "Mean", age[4],
    "3+", n_pct(data$CHILD1, "3+"),
    "", "",
    "No answer", n_pct(data$REVENU1, "No answer"),
    "Professional program", n_pct(education, "Professional program"),
    "", "",
    
    "", "Agree", n_pct(data$PARKIMPT, "En accord"),
    
    "", "Mean",
    cipqay[4], parks[4], area[4],
    
    "", "Agree",
    n_pct(data$PERACC, "En accord"),
    n_pct(data$PERNUMB, "En accord"),
    n_pct(data$PERSATIS, "En accord"),
    "Moderate", n_pct(data$PEREVA, "Moyen")
  ),
  
  C5 = c(
    "", "", "",
    "Std dev.", age[5],
    "", "",
    "", "",
    "", "",
    "University", n_pct(education, "University"),
    "", "",
    
    "", "Strongly agree", n_pct(data$PARKIMPT, "Fortement en accord"),
    
    "", "Std dev.",
    cipqay[5], parks[5], area[5],
    
    "", "Strongly agree",
    n_pct(data$PERACC, "Fortement en accord"),
    n_pct(data$PERNUMB, "Fortement en accord"),
    n_pct(data$PERSATIS, "Fortement en accord"),
    "High", n_pct(data$PEREVA, "Élevé")
  ),
  
  C6 = c(
    "", "", "",
    "", "",
    "", "",
    "", "",
    "", "",
    "", "",
    "", "",
    
    "", "", "",
    
    "", "",
    "", "", "",
    
    "", "",
    "", "", "",
    "Very high", n_pct(data$PEREVA, "Très élevé")
  )
)

write.csv(
  table3,
  "Table3_descriptive_characteristics.csv",
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

# Clustering

cluster_data <- data %>%
  select(
    id_1,
    distance_min,
    distance_max,
    PARKFREQ,
    nb_parcs_visites,
    MODEPARK,
    CONTEXTE,
    TYPE
  ) %>%
  mutate(
    distance_min = scale(distance_min)[, 1],
    distance_max = scale(distance_max)[, 1],
    PARKFREQ = scale(PARKFREQ)[, 1],
    nb_parcs_visites = scale(nb_parcs_visites)[, 1],
    MODEPARK = as.factor(MODEPARK),
    CONTEXTE = as.factor(CONTEXTE),
    TYPE = as.factor(TYPE)
  )

set.seed(102)

k_values <- 2:10

tot_withinss <- sapply(k_values, function(k) {
  model <- kproto(
    cluster_data %>% select(-id_1),
    k = k
  )
  
  model$tot.withinss
})

plot(
  k_values,
  tot_withinss,
  type = "b",
  pch = 19,
  xlab = "Number of clusters (k)",
  ylab = "Total within-cluster distance",
  main = "Elbow method"
)

set.seed(102)


final_model <- kproto(
  cluster_data %>% select(-id_1),
  k = 6
)

cluster_data$cluster <- final_model$cluster

data_clustered <- data %>%
  left_join(cluster_data %>% select(id_1, cluster), by = "id_1")

write_csv(
  data_clustered,
  "park_accessibility_dataset_clustered.csv"
)

# Cluster summary

effectifs <- cluster_data %>%
  count(cluster, name = "N") %>%
  mutate(Proportion = round(N / sum(N) * 100, 1))

print(effectifs)


# Data preparation

likert_levels <- c(
  "Strongly disagree",
  "Disagree",
  "Neutral",
  "Agree",
  "Strongly agree"
)

evaluation_levels <- c(
  "Null",
  "Very low",
  "Low",
  "Medium",
  "High",
  "Very high"
)

cluster_labels <- c(
  "A" = "A: Frequent large and neighborhood park visitors (n=133)",
  "B" = "B: Frequent single-park visitors (n=78)",
  "C" = "C: Various large park visitors (n=144)",
  "D" = "D: Distant-park car users (n=60)",
  "E" = "E: Infrequent large and neighborhood park visitors (n=158)",
  "F" = "F: Infrequent single-park car users (n=65)"
)

translate_likert <- function(x) {
  case_when(
    x == "Fortement en désaccord" ~ "Strongly disagree",
    x == "En désaccord" ~ "Disagree",
    x == "Neutre" ~ "Neutral",
    x == "En accord" ~ "Agree",
    x == "Fortement en accord" ~ "Strongly agree",
    TRUE ~ NA_character_
  )
}

translate_evaluation <- function(x) {
  case_when(
    x == "Nul" ~ "Null",
    x == "Très faible" ~ "Very low",
    x == "Faible" ~ "Low",
    x == "Moyen" ~ "Medium",
    x == "Élevé" ~ "High",
    x == "Très élevé" ~ "Very high",
    TRUE ~ NA_character_
  )
}

data_clustered <- data_clustered %>%
  mutate(
    PARKIMPT = factor(
      PARKIMPT,
      levels = c(
        "Fortement en désaccord",
        "En désaccord",
        "Neutre",
        "En accord",
        "Fortement en accord"
      )
    ),
    AUTO1 = factor(
      AUTO1,
      levels = c(
        "Aucune",
        "1 automobile",
        "2 automobiles et plus"
      )
    ),
    PERNUMB = factor(
      PERNUMB,
      levels = c(
        "Fortement en désaccord",
        "En désaccord",
        "Neutre",
        "En accord",
        "Fortement en accord"
      )
    ),
    PERSATIS = factor(
      PERSATIS,
      levels = c(
        "Fortement en désaccord",
        "En désaccord",
        "Neutre",
        "En accord",
        "Fortement en accord"
      )
    ),
    PEREVA = factor(
      PEREVA,
      levels = c(
        "Nul",
        "Très faible",
        "Faible",
        "Moyen",
        "Élevé",
        "Très élevé"
      )
    ),
    PERNUMB_en = factor(
      translate_likert(PERNUMB),
      levels = likert_levels
    ),
    PERSATIS_en = factor(
      translate_likert(PERSATIS),
      levels = likert_levels
    ),
    PEREVA_en = factor(
      translate_evaluation(PEREVA),
      levels = evaluation_levels
    ),
    cluster_fac = case_when(
      cluster == 2 ~ "A",
      cluster == 4 ~ "B",
      cluster == 1 ~ "C",
      cluster == 3 ~ "D",
      cluster == 5 ~ "E",
      cluster == 6 ~ "F",
      TRUE ~ NA_character_
    ),
    cluster_fac = factor(
      cluster_fac,
      levels = names(cluster_labels)
    ),
    cluster_fac_name = factor(
      unname(cluster_labels[as.character(cluster_fac)]),
      levels = unname(cluster_labels)
    )
  )

data_clustered %>%
  count(cluster_fac, name = "N") %>%
  mutate(Proportion = round(N / sum(N) * 100, 1))

# Statistical tests

format_p_value <- function(p, digits = 4) {
  if (is.na(p)) {
    return("p-value = NA")
  }
  
  if (p < 0.001) {
    return("p-value < 0.001")
  }
  
  paste0(
    "p-value = ",
    formatC(p, format = "f", digits = digits)
  )
}

auto_p <- chisq.test(table(data_clustered$AUTO1, data_clustered$cluster_fac))$p.value

importance_p <- chisq.test(
  table(data_clustered$PARKIMPT, data_clustered$cluster_fac)
)$p.value

sex_p <- chisq.test(
  table(data_clustered$SEX, data_clustered$cluster_fac)
)$p.value

children_p <- chisq.test(
  table(data_clustered$CHILD1, data_clustered$cluster_fac)
)$p.value

income_p <- chisq.test(
  table(data_clustered$REVENU1, data_clustered$cluster_fac)
)$p.value

education_p <- chisq.test(
  table(data_clustered$EDUC, data_clustered$cluster_fac)
)$p.value

ethnicity_p <- chisq.test(
  table(data_clustered$ETHNICITE1, data_clustered$cluster_fac)
)$p.value

cipqay_model <- aov(CIPQAY ~ cluster_fac, data = data_clustered)
cipqay_p <- summary(cipqay_model)[[1]][["Pr(>F)"]][1]

parks_model <- aov(ACCOBJ_CP ~ cluster_fac, data = data_clustered)
parks_p <- summary(parks_model)[[1]][["Pr(>F)"]][1]

area_model <- aov(ACCAREA_CP ~ cluster_fac, data = data_clustered)
area_p <- summary(area_model)[[1]][["Pr(>F)"]][1]

PERNUMB_p <- chisq.test(
  table(data_clustered$PERNUMB, data_clustered$cluster_fac)
)$p.value

persatis_p <- chisq.test(
  table(data_clustered$PERSATIS, data_clustered$cluster_fac)
)$p.value

pereva_p <- chisq.test(
  table(data_clustered$PEREVA, data_clustered$cluster_fac)
)$p.value

peracc_p <- chisq.test(
  table(data_clustered$PERACC, data_clustered$cluster_fac)
)$p.value
age_model <- aov(AGE ~ cluster_fac, data = data_clustered)

age_p <- summary(age_model)[[1]][["Pr(>F)"]][1]

# Sociodemographic profiles

theme_common <- theme(
  axis.title.x = element_text(size = 18),
  axis.title.y = element_text(size = 18),
  axis.text.x = element_text(size = 16),
  axis.text.y = element_text(size = 16),
  legend.title = element_blank(),
  legend.text = element_text(size = 22),
  legend.key.height = grid::unit(1.2, "cm"),
  legend.key.width = grid::unit(1.2, "cm"),
  plot.title = element_text(size = 22),
  plot.margin = margin(b = 1, unit = "cm")
)


graph_importance <- ggplot(data_clustered, aes(fill = PARKIMPT, x = cluster_fac)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(
    values = c(
      "Fortement en désaccord" = "#2EB62C",
      "En désaccord" = "#57C84D",
      "Neutre" = "#83D475",
      "En accord" = "#ABE098",
      "Fortement en accord" = "#C5E8B7"
    ),
    labels = c(
      "Fortement en désaccord" = "Strongly disagree",
      "En désaccord" = "Disagree",
      "Neutre" = "Neutral",
      "En accord" = "Agree",
      "Fortement en accord" = "Strongly agree"
    ),
    na.value = "gray",
    na.translate = TRUE
  ) +
  labs(
    title = "Self-declared park importance",
    x = paste(
      "Cluster",
      format_p_value(importance_p),
      sep = "\n"
    ),
    y = "Proportion of individuals",
    fill = NULL
  ) +
  theme_common

graph_age <- ggplot(data_clustered, aes(cluster_fac, AGE)) +
  geom_violin() +
  geom_boxplot(width = 0.1) +
  labs(
    title = "Age",
    x = paste(
      "Cluster",
      format_p_value(age_p),
      sep = "\n"
    ),
    y = "Age"
  ) +
  theme_common

# Perceived accessibility

graph_PERNUMB <- ggplot(data_clustered, aes(fill = PERNUMB_en, x = cluster_fac)) +
  geom_bar(position = "fill") +
  scale_fill_brewer(palette = "Blues") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = format_p_value(PERNUMB_p, digits = 2),
    y = "Perceived number of parks accessible",
    fill = NULL
  ) +
  theme_common

graph_persatis <- ggplot(data_clustered, aes(fill = PERSATIS_en, x = cluster_fac)) +
  geom_bar(position = "fill") +
  scale_fill_brewer(palette = "Blues") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = format_p_value(persatis_p),
    y = "Satisfied with park access",
    fill = NULL
  ) +
  theme_common

graph_pereva <- ggplot(data_clustered, aes(fill = PEREVA_en, x = cluster_fac)) +
  geom_bar(position = "fill") +
  scale_fill_brewer(palette = "Reds") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = format_p_value(pereva_p),
    y = "Overall evaluation of park access",
    fill = NULL
  ) +
  theme_common

panel_accessibility <- (
  graph_persatis +
    plot_spacer() +
    graph_pereva +
    guide_area()
) +
  plot_layout(
    ncol = 4,
    widths = c(1, 0.10, 1, 1.2),
    guides = "collect"
  ) &
  theme(
    legend.title = element_blank(),
    legend.text = element_text(size = 22),
    legend.key.height = grid::unit(1.2, "cm"),
    legend.key.width = grid::unit(1.2, "cm")
  )

# Standardized cluster profiles

vars_std <- c(
  "distance_min",
  "distance_max",
  "PARKFREQ",
  "nb_parcs_visites"
)

# Create respondent-level binary indicators
# and standardize all variables at the respondent level

data_std <- data_clustered %>%
  mutate(
    car = as.numeric(MODEPARK == "Automobile"),
    contexte = as.numeric(
      CONTEXTE == "Aucun comportement prédominant"
    ),
    grandparc = as.numeric(
      Grand_parc == 1 | Grand_parc == "1"
    ),
    parcarr = as.numeric(
      Parc_arrondissement == 1 | Parc_arrondissement == "1"
    )
  ) %>%
  mutate(
    across(
      c(
        all_of(vars_std),
        car,
        contexte,
        grandparc,
        parcarr
      ),
      ~ as.numeric(scale(.x))
    )
  )

# Cluster means

profil_cluster <- data_std %>%
  group_by(cluster_fac_name) %>%
  summarise(
    across(
      c(
        all_of(vars_std),
        car,
        contexte,
        grandparc,
        parcarr
      ),
      ~ mean(.x, na.rm = TRUE)
    ),
    n = n(),
    .groups = "drop"
  ) %>%
  mutate(
    prop = n / sum(n),
    label_prop = percent(prop, accuracy = 0.1)
  ) %>%
  rename(
    prop_car = car,
    prop_contexte = contexte,
    prop_grandparc = grandparc,
    prop_parcarr = parcarr
  )

# Cluster profiles


profil_long <- profil_cluster %>%
  dplyr::select(-n, -prop, -label_prop) %>%
  pivot_longer(
    cols = -cluster_fac_name,
    names_to = "variable",
    values_to = "value"
  ) %>%
  mutate(
    variable = factor(
      variable,
      levels = c(
        "PARKFREQ",
        "nb_parcs_visites",
        "distance_min",
        "distance_max",
        "prop_car",
        "prop_contexte",
        "prop_grandparc",
        "prop_parcarr"
      )
    ),
    value_plot = if_else(
      variable == "prop_contexte",
      -value,
      value
    )
  )

plot_standard <- ggplot(
  profil_long,
  aes(x = variable, y = value_plot, fill = variable)
) +
  geom_col(width = 0.8) +
  facet_wrap(
    ~ cluster_fac_name,
    nrow = 1,
    labeller = label_wrap_gen(width = 18)
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed",
    colour = "grey40"
  ) +
  scale_y_continuous(limits = c(-3, 3)) +
  scale_fill_manual(
    values = c(
      "distance_min" = "#3A50F0",
      "distance_max" = "#8494F5",
      "PARKFREQ" = "#EDC055",
      "nb_parcs_visites" = "#997317",
      "prop_car" = "#82E0C6",
      "prop_contexte" = "#08805E",
      "prop_grandparc" = "#E67780",
      "prop_parcarr" = "#7A010C"
    ),
    labels = c(
      "distance_min" = "Minimum distance",
      "distance_max" = "Maximum distance",
      "PARKFREQ" = "Park visit frequency",
      "nb_parcs_visites" = "Number of different parks visited",
      "prop_car" = "Car use*",
      "prop_contexte" = "Predominant behavior**",
      "prop_grandparc" = "Large park***",
      "prop_parcarr" = "Neighborhood park***"
    )
  ) +
  labs(
    x = NULL,
    y = "Standardized mean",
    fill = NULL
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "bottom",
    panel.spacing = grid::unit(1, "cm"),
    axis.text.y = element_text(size = 26),
    axis.title.y = element_text(size = 26),
    legend.text = element_text(size = 28),
    legend.key.height = grid::unit(1.2, "cm"),
    legend.key.width = grid::unit(1.2, "cm"),
    strip.text.x = element_text(size = 24, face = "bold"),
    text = element_text(family = "Arial")
  )

# Pairwise age comparisons

significance_code <- function(p) {
  case_when(
    is.na(p) ~ "",
    p < 0.001 ~ "***",
    p < 0.01 ~ "**",
    p < 0.05 ~ "*",
    p < 0.1 ~ ".",
    TRUE ~ ""
  )
}

tukey_pairs <- TukeyHSD(
  age_model,
  "cluster_fac"
)$cluster_fac %>%
  as.data.frame() %>%
  tibble::rownames_to_column("comparison") %>%
  separate(
    comparison,
    into = c("g1", "g2"),
    sep = "-"
  ) %>%
  mutate(
    significance = significance_code(`p adj`)
  )

clusters <- sort(unique(c(tukey_pairs$g1, tukey_pairs$g2)))

tukey_matrix <- matrix(
  "",
  nrow = length(clusters),
  ncol = length(clusters),
  dimnames = list(clusters, clusters)
)

for (i in seq_len(nrow(tukey_pairs))) {
  tukey_matrix[tukey_pairs$g1[i], tukey_pairs$g2[i]] <-
    tukey_pairs$significance[i]
}

diag(tukey_matrix) <- "-"

tukey_age <- as.data.frame(tukey_matrix) %>%
  tibble::rownames_to_column("group")

# Save figures

ggsave(
  filename = "cluster_profiles_standardized.jpg",
  plot = plot_standard,
  width = 22,
  height = 8,
  units = "in",
  dpi = 300
)

ggsave(
  filename = "park_importance_by_cluster.jpg",
  plot = graph_importance,
  width = 9,
  height = 7,
  units = "in",
  dpi = 300
)


ggsave(
  filename = "perceived_accessibility_by_cluster.jpg",
  plot = panel_accessibility,
  width = 14,
  height = 7,
  units = "in",
  dpi = 300
)

ggsave(
  filename = "age_by_cluster.jpg",
  plot = graph_age,
  width = 9,
  height = 7,
  units = "in",
  dpi = 300
)

write.csv(
  tukey_age,
  file = "tukey_age.csv",
  row.names = FALSE
)