Cluster ANOVA Results and Outcome Distributions

1 Data preparation

The source variables are recoded using the rules in the supplied R script. Binary outcomes are coded 0 = No and 1 = Yes. The severe-weather concern outcome is coded from 1 = Low worry to 3 = High worry. Missing values are excluded separately for each ANOVA and plot.

resolve_input_file <- function(requested_path, alternative_names, file_description) {
  user_profile <- Sys.getenv("USERPROFILE", unset = path.expand("~"))

  input_document <- tryCatch(
    knitr::current_input(dir = TRUE),
    error = function(e) ""
  )

  input_directory <- if (nzchar(input_document)) {
    dirname(normalizePath(input_document, winslash = "/", mustWork = FALSE))
  } else {
    getwd()
  }

  search_directories <- unique(c(
    getwd(),
    input_directory,
    file.path(user_profile, "Downloads"),
    file.path(user_profile, "Desktop"),
    file.path(user_profile, "OneDrive", "Desktop"),
    path.expand("~/Downloads"),
    path.expand("~/Desktop"),
    path.expand("~/OneDrive/Desktop")
  ))

  candidate_paths <- unique(c(
    requested_path,
    unlist(
      lapply(
        search_directories,
        function(directory) file.path(directory, alternative_names)
      ),
      use.names = FALSE
    )
  ))

  matching_paths <- candidate_paths[file.exists(candidate_paths)]

  if (length(matching_paths) == 0) {
    stop(
      paste0(
        "Could not find the ", file_description, ".\n\n",
        "Either place it in the same folder as this R Markdown file, ",
        "or change its path in the params section at the top of the file.\n\n",
        "Paths checked:\n- ",
        paste(candidate_paths, collapse = "\n- ")
      ),
      call. = FALSE
    )
  }

  normalizePath(matching_paths[[1]], winslash = "/", mustWork = TRUE)
}

sav_path <- resolve_input_file(
  requested_path = params$sav_file,
  alternative_names = c("wrp_25+.sav", "wrp_25+(1).sav"),
  file_description = "SPSS data file"
)

cluster_path <- resolve_input_file(
  requested_path = params$cluster_file,
  alternative_names = "clusters.xlsx",
  file_description = "cluster spreadsheet"
)

dataframe <- read_sav(sav_path)
clusters <- read_excel(cluster_path)

required_data_variables <- c(
  "COUNTRYNEW", "WP24198", "WP24199", "WP23345", "WP22252",
  "WP22445", "WP24213", "WP20719", "WP20723", "WP24181",
  "WP24182", "WP24183", "WP24184", "WP24185", "WP24186",
  "WP24187", "WP24188"
)

required_cluster_variables <- c("COUNTRYNEW", "Cluster")

missing_data_variables <- setdiff(required_data_variables, names(dataframe))
missing_cluster_variables <- setdiff(required_cluster_variables, names(clusters))

if (length(missing_data_variables) > 0) {
  stop(
    "The SAV file is missing these required variables: ",
    paste(missing_data_variables, collapse = ", ")
  )
}

if (length(missing_cluster_variables) > 0) {
  stop(
    "The cluster file is missing these required variables: ",
    paste(missing_cluster_variables, collapse = ", ")
  )
}

if (anyDuplicated(clusters$COUNTRYNEW) > 0) {
  stop("The cluster file contains duplicate COUNTRYNEW values.")
}

dataframe_combined <- right_join(
  dataframe,
  clusters,
  by = "COUNTRYNEW"
) %>%
  mutate(
    Cluster = factor(Cluster, levels = sort(unique(na.omit(Cluster)))),

    prepared_natgov = case_when(
      WP24198 == 1 ~ 1,
      WP24198 %in% c(2, 3) ~ 0,
      TRUE ~ NA_real_
    ),
    prepared_localgov = case_when(
      WP24199 == 1 ~ 1,
      WP24199 %in% c(2, 3) ~ 0,
      TRUE ~ NA_real_
    ),
    planknown = case_when(
      WP23345 == 1 ~ 1,
      WP23345 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    protectfam = case_when(
      WP22252 == 1 ~ 1,
      WP22252 %in% c(2, 3) ~ 0,
      TRUE ~ NA_real_
    ),
    severeweatherexp = case_when(
      WP22445 %in% c(1, 2, 3) ~ 1,
      WP22445 == 4 ~ 0,
      TRUE ~ NA_real_
    ),
    disasterexp = case_when(
      WP24213 == 1 ~ 1,
      WP24213 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    climateconcern = case_when(
      WP20719 == 1 ~ 1,
      WP20719 %in% c(2, 3, 98) ~ 0,
      TRUE ~ NA_real_
    ),
    weatherconcern = case_when(
      WP20723 == 1 ~ 3,
      WP20723 == 2 ~ 2,
      WP20723 == 3 ~ 1,
      TRUE ~ NA_real_
    ),
    warning_internet = case_when(
      WP24181 == 1 ~ 1,
      WP24181 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_radio = case_when(
      WP24182 == 1 ~ 1,
      WP24182 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_tv = case_when(
      WP24183 == 1 ~ 1,
      WP24183 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_newspapers = case_when(
      WP24184 == 1 ~ 1,
      WP24184 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_sms = case_when(
      WP24185 == 1 ~ 1,
      WP24185 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_whatsapp = case_when(
      WP24186 == 1 ~ 1,
      WP24186 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_billboard = case_when(
      WP24187 == 1 ~ 1,
      WP24187 == 2 ~ 0,
      TRUE ~ NA_real_
    ),
    warning_loudspeaker = case_when(
      WP24188 == 1 ~ 1,
      WP24188 == 2 ~ 0,
      TRUE ~ NA_real_
    )
  )

warning_outcomes <- c(
  warning_internet = "Received warning: Internet",
  warning_radio = "Received warning: Radio",
  warning_tv = "Received warning: TV",
  warning_newspapers = "Received warning: Newspapers",
  warning_sms = "Received warning: SMS",
  warning_whatsapp = "Received warning: Whatsapp",
  warning_billboard = "Received warning: Billboard",
  warning_loudspeaker = "Received warning; Loudspeaker"
)

preparedness_outcomes <- c(
  prepared_natgov = "Prepared for disaster: National govt",
  prepared_localgov = "Prepared for disaster; Local govt",
  planknown = "Plan for disaster known to household members",
  protectfam = "Could protect family in disaster"
)

risk_perception_outcomes <- c(
  climateconcern = paste(
    "Concerned climate change will be a threat to country",
    "in next 20 years"
  ),
  weatherconcern = paste(
    "Worried severe weather events could cause",
    "significant harm"
  )
)

experience_outcomes <- c(
  disasterexp = "Experienced disaster in last 5 years",
  severeweatherexp = "Experiences severe weather event in last 2 years"
)

format_p_value <- function(x) {
  case_when(
    is.na(x) ~ NA_character_,
    x < .001 ~ "< .001",
    TRUE ~ sprintf("%.3f", x)
  )
}

run_anovas <- function(data, outcomes) {
  map_dfr(names(outcomes), function(variable) {
    model_data <- data %>%
      transmute(
        Cluster,
        outcome = .data[[variable]]
      ) %>%
      drop_na()

    if (n_distinct(model_data$Cluster) < 2) {
      return(tibble(
        Outcome = unname(outcomes[[variable]]),
        N = nrow(model_data),
        df_between = NA_real_,
        df_within = NA_real_,
        F = NA_real_,
        p = NA_real_,
        eta_squared = NA_real_
      ))
    }

    model <- aov(outcome ~ Cluster, data = model_data)
    model_table <- broom::tidy(model)

    between_row <- model_table %>% filter(term == "Cluster")
    residual_row <- model_table %>% filter(term == "Residuals")
    total_ss <- sum(model_table$sumsq, na.rm = TRUE)

    tibble(
      Outcome = unname(outcomes[[variable]]),
      N = nrow(model_data),
      df_between = between_row$df,
      df_within = residual_row$df,
      F = between_row$statistic,
      p = between_row$p.value,
      eta_squared = between_row$sumsq / total_ss
    )
  })
}

print_anova_table <- function(results, caption) {
  results %>%
    transmute(
      Outcome,
      N,
      `df between` = as.integer(df_between),
      `df within` = as.integer(df_within),
      `F statistic` = sprintf("%.2f", F),
      `p value` = format_p_value(p),
      `Eta squared` = sprintf("%.3f", eta_squared)
    ) %>%
    knitr::kable(
      caption = caption,
      align = c("l", "r", "r", "r", "r", "r", "r")
    )
}

plot_outcome_distribution <- function(data, variable, outcome_label) {
  plot_data <- data %>%
    filter(
      !is.na(Cluster),
      !is.na(.data[[variable]])
    )

  if (variable == "weatherconcern") {
    plot_data <- plot_data %>%
      mutate(
        Response = factor(
          .data[[variable]],
          levels = c(1, 2, 3),
          labels = c("Low worry", "Moderate worry", "High worry")
        )
      )
  } else {
    plot_data <- plot_data %>%
      mutate(
        Response = factor(
          .data[[variable]],
          levels = c(0, 1),
          labels = c("No", "Yes")
        )
      )
  }

  distribution_data <- plot_data %>%
    count(Cluster, Response, .drop = FALSE) %>%
    group_by(Cluster) %>%
    mutate(Proportion = n / sum(n)) %>%
    ungroup()

  ggplot(
    distribution_data,
    aes(
      x = Cluster,
      y = Proportion,
      fill = Response
    )
  ) +
    geom_col(width = 0.72) +
    scale_y_continuous(
      labels = label_percent(accuracy = 1),
      limits = c(0, 1),
      expand = expansion(mult = c(0, 0.02))
    ) +
    labs(
      title = outcome_label,
      subtitle = paste0(
        "Valid responses: ",
        scales::comma(sum(distribution_data$n))
      ),
      x = "Cluster",
      y = "Percentage of respondents",
      fill = NULL
    ) +
    theme_tufte(base_size = 12) +
    theme(
      legend.position = "top",
      legend.justification = "left",
      plot.title.position = "plot"
    )
}

render_outcome_plots <- function(data, outcomes) {
  purrr::iwalk(outcomes, function(outcome_label, variable) {
    cat("\n\n### ", outcome_label, "\n\n", sep = "")
    print(plot_outcome_distribution(data, variable, outcome_label))
    cat("\n\n")
  })
}

2 Warnings

2.1 ANOVA results

warning_anova_results <- run_anovas(
  dataframe_combined,
  warning_outcomes
)

print_anova_table(
  warning_anova_results,
  "ANOVA results: warnings"
)

ANOVA results: warnings

Outcome	N	df between	df within	F statistic	p value	Eta squared
Received warning: Internet	17867	5	17861	280.01	< .001	0.073
Received warning: Radio	17606	5	17600	123.92	< .001	0.034
Received warning: TV	17819	5	17813	157.29	< .001	0.042
Received warning: Newspapers	17353	5	17347	84.48	< .001	0.024
Received warning: SMS	17454	5	17448	58.83	< .001	0.017
Received warning: Whatsapp	17805	5	17799	311.65	< .001	0.080
Received warning: Billboard	17455	5	17449	60.36	< .001	0.017
Received warning; Loudspeaker	17786	5	17780	196.25	< .001	0.052

2.2 Outcome distributions by cluster

render_outcome_plots(
  dataframe_combined,
  warning_outcomes
)

2.2.1 Received warning: Internet

2.2.2 Received warning: Radio

2.2.3 Received warning: TV

2.2.4 Received warning: Newspapers

2.2.5 Received warning: SMS

2.2.6 Received warning: Whatsapp

2.2.7 Received warning: Billboard

2.2.8 Received warning; Loudspeaker

3 Preparedness

3.1 ANOVA results

preparedness_anova_results <- run_anovas(
  dataframe_combined,
  preparedness_outcomes
)

print_anova_table(
  preparedness_anova_results,
  "ANOVA results: preparedness"
)

ANOVA results: preparedness

Outcome	N	df between	df within	F statistic	p value	Eta squared
Prepared for disaster: National govt	103610	5	103604	2237.98	< .001	0.097
Prepared for disaster; Local govt	103721	5	103715	1731.20	< .001	0.077
Plan for disaster known to household members	112139	5	112133	723.53	< .001	0.031
Could protect family in disaster	111049	5	111043	1041.19	< .001	0.045

3.2 Outcome distributions by cluster

render_outcome_plots(
  dataframe_combined,
  preparedness_outcomes
)

3.2.1 Prepared for disaster: National govt

3.2.2 Prepared for disaster; Local govt

3.2.3 Plan for disaster known to household members

3.2.4 Could protect family in disaster

4 Risk perceptions

4.1 ANOVA results

risk_perception_anova_results <- run_anovas(
  dataframe_combined,
  risk_perception_outcomes
)

print_anova_table(
  risk_perception_anova_results,
  "ANOVA results: risk perceptions"
)

ANOVA results: risk perceptions

Outcome	N	df between	df within	F statistic	p value	Eta squared
Concerned climate change will be a threat to country in next 20 years	115517	5	115511	686.19	< .001	0.029
Worried severe weather events could cause significant harm	114822	5	114816	1154.51	< .001	0.048

4.2 Outcome distributions by cluster

render_outcome_plots(
  dataframe_combined,
  risk_perception_outcomes
)

4.2.1 Concerned climate change will be a threat to country in next 20 years

4.2.2 Worried severe weather events could cause significant harm

5 Experience

5.1 ANOVA results

experience_anova_results <- run_anovas(
  dataframe_combined,
  experience_outcomes
)

print_anova_table(
  experience_anova_results,
  "ANOVA results: experience"
)

ANOVA results: experience

Outcome	N	df between	df within	F statistic	p value	Eta squared
Experienced disaster in last 5 years	115141	5	115135	348.90	< .001	0.015
Experiences severe weather event in last 2 years	114949	5	114943	387.64	< .001	0.017

5.2 Outcome distributions by cluster

render_outcome_plots(
  dataframe_combined,
  experience_outcomes
)

5.2.1 Experienced disaster in last 5 years

5.2.2 Experiences severe weather event in last 2 years