Identifying different types of green space deprivation with K means clustering and TidyModels

# reshape data for plotting
foe_green_space_focus_long <- foe_green_space_focus %>%
  pivot_longer(
    cols = population_density:pcnt_bame_pop,
    names_to = "variable",
    values_to = "value"
  )

# get a sense of distributions
ggplot(
  foe_green_space_focus_long,
  aes(value)
) +
  geom_density() +
  facet_wrap(~variable, scales = "free") +
  theme(
    strip.background = element_rect(fill = "white"),
    strip.text = element_text(colour = "black")
  )

#' Get a nice looking table of skew and kurtosis test results
#'
#' @param df a dataframe of variables
#'
#' @return a dataframe of test results
#' @export
#'
#' @examples
test_skew_kurtosis <- function(df) {
  df %>%
    # calculate summary statistics
    summarise(across(dplyr::everything(),
      list(skewness = skewness, kurtosis = kurtosis),
      na.rm = TRUE
    )) %>%
    # rearrange dataframe into a nice format
    pivot_longer(cols = everything()) %>%
    mutate(
      test = if_else(str_detect(name, "skewness"),
        "skewness", "kurtosis"
      ),
      name = str_remove(name, c("_skewness", "_kurtosis"))
    ) %>%
    rename(variable = name) %>%
    pivot_wider(names_from = "test", values_from = "value") %>%
    arrange(desc(skewness), desc(kurtosis))
}

skew_kurtosis_test_res <- test_skew_kurtosis(select(
  foe_green_space_focus,
  -msoa_code
))

kable(skew_kurtosis_test_res)

prop_above_threshold <- function(variable, threshold) {
  obs_above_thres <- length(which(foe_green_space_focus[variable] > threshold))
  total_obs <- length(foe_green_space_focus[[variable]])
  obs_above_thres / total_obs
}

# -----------------------------------------------------------------------------
# garden_area_per_capita: Look at proportions of observations above thresholds
# -----------------------------------------------------------------------------
cut_off_thresholds <- c(1, 10, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000)

prop_omitted <- tibble(threshold = cut_off_thresholds) %>%
  mutate(
    prop_omitted_garden_area_per_capita = map_dbl(
      threshold,
      ~ prop_above_threshold("garden_area_per_capita", .)
    )
  )

kable(prop_omitted)

# -----------------------------------------------------------------------------
# green_space_area_per_capita: Look at proportions of observations above thresholds
# -----------------------------------------------------------------------------

cut_off_thresholds <- c(1, 10, 50, 100, 200, 300, 400, 500, 1000, 5000, 8000, 10000, 100000)

prop_omitted <- tibble(threshold = cut_off_thresholds) %>%
  mutate(
    prop_omitted_green_space_area_per_capita = map_dbl(
      threshold,
      ~ prop_above_threshold("green_space_area_per_capita", .)
    )
  )

kable(prop_omitted)

# retain outliers for inspection and use later on
foe_green_space_outliers <- foe_green_space_focus %>%
  filter(
    garden_area_per_capita > 500 |
      green_space_area_per_capita > 5000
  )

# remove outliers ahead of feature engineering workflow
foe_green_space_focus <- foe_green_space_focus %>%
  filter(
    garden_area_per_capita <= 500,
    green_space_area_per_capita <= 5000
  )


# rerun the skew and kurtosis tests for the variables where I've removed some
# outliers
test_skew_kurtosis(select(
  foe_green_space_focus, garden_area_per_capita, green_space_area_per_capita
)) %>% 
  kable()

set.seed(1289473)

# ------------------------------------------------------------------------------
# key parameters
# ------------------------------------------------------------------------------

#' Run the K means cluster workflow with a single recipe (k is a tuning parameter)
#'
#' @data_df the dataset to be clustered
#' @param num_cluster_levels number of levels created between 1 and 10 in the parameter (i.e. k) tuning grid
#' @param num_cross_validation_folds number of cv folds
#' @recipe the recipe to be used to pre-process the data ahead of loading it into the clustering workflow
#'
#' @return a dataframe of the results of the tuning process
#' @export
#'
#' @examples
run_cluster_wf <- function(data_df, num_cluster_levels = 5,
                           num_cross_validation_folds = 2,
                           recipe) {
  # ---------------------------------------------------------------------------
  # create cross-validation samples
  green_space_cv <- vfold_cv(data_df, v = num_cross_validation_folds)
  # ---------------------------------------------------------------------------

  # ---------------------------------------------------------------------------
  # set up clustering algorithm
  # ---------------------------------------------------------------------------
  kmeans_spec <- k_means(num_clusters = tune()) %>%
    set_engine("stats")

  # ---------------------------------------------------------------------------
  # set up workflow which brings together the clustering algorithm and
  # the data pre-processing (i.e. the recipe) for execution
  # ---------------------------------------------------------------------------
  kmeans_wflow <- workflow(recipe, kmeans_spec)

  # ---------------------------------------------------------------------------
  # create parameter tuning grid - for k (i.e. the number of clusters)
  # ---------------------------------------------------------------------------
  num_clusters_grid <- grid_regular(num_clusters(), levels = num_cluster_levels)

  # ---------------------------------------------------------------------------
  # run tuning
  # ---------------------------------------------------------------------------

  res <- tune_cluster(kmeans_wflow,
    resamples = green_space_cv,
    grid = num_clusters_grid,
    control = control_grid(save_pred = TRUE, extract = identity),
    metrics = cluster_metric_set(
      sse_within_total, sse_total,
      sse_ratio
    )
  )

  return(res)
}


#' Run the K means cluster workflow for multiple recipes (k is a tuning parameter)
#'
#' #' @data_df the dataset to be clustered
#' @param num_cluster_levels number of levels created between 1 and 10 in the parameter (i.e. k) tuning grid
#' @param num_cross_validation_folds number of cv folds
#' @param recipes a named list of recipes
#'
#' @return a ggplot graph showing the performance (i.e. sse_ratio) of clustering with each recipe
#' @export
#'
#' @examples
run_wf_with_recipes <- function(data_df,
                                num_cluster_levels = 5,
                                num_cross_validation_folds = 2,
                                recipes) {
  # ---------------------------------------------------------------------------
  # create a dataframe with each recipe
  # to store the objects create at each stage of the clustering process
  # ---------------------------------------------------------------------------

  clustering_df <- tibble(
    recipe_id = names(recipes), # for meaningful names in plotting
    recipe = recipes
  )

  # ---------------------------------------------------------------------------
  # Run the workflow with each recipe and get the performance metrics for
  # each
  # ---------------------------------------------------------------------------

  metric_df <- clustering_df %>%
    mutate(
      # run a clustering workflow for each recipe
      workflow = map(recipe, ~ run_cluster_wf(
        recipe = .,
        data_df = data_df,
        num_cluster_levels = num_cluster_levels,
        num_cross_validation_folds = num_cross_validation_folds
      )),

      # get the performance metrics for each workflow
      metrics = map(workflow, ~ collect_metrics(.))
    )

  # ---------------------------------------------------------------------------
  # transform the metric data for use in ggplot2
  # including selecting the metric of interest
  # ---------------------------------------------------------------------------

  metric_df_focus <- metric_df %>%
    # expand the nested metrics tables
    select(recipe_id, metrics) %>%
    unnest(metrics) %>%
    # focus on metric of interest
    filter(.metric == "sse_ratio") %>%
    # so that different coloured lines can be plotted for each recipe
    mutate(recipe_id = factor(recipe_id))

  # ---------------------------------------------------------------------------
  # plot the performance metric for each recipe
  # ---------------------------------------------------------------------------

  p <- ggplot(
    data = metric_df_focus,
    mapping = aes(
      x = num_clusters, y = mean,
      colour = recipe_id, group = recipe_id
    )
  ) +
    geom_line() +
    scale_color_viridis_d() +
    theme_minimal() +
    ylab("Performance (mean WSS/TSS ratio)") +
    xlab("Number of clusters") +
    scale_x_continuous(breaks = 1:10)

  # make plot interactive
  plotly::ggplotly(p)
}

#' Having decided on a recipe (set of pre-processing steps) and k, run the
#' clustering workflow
#'
#' @param k the number of clusters to identified by the k means clustering
#' algorithm
#' @param recipe the set of data pre-processing steps (a TidyModels object)
#' @param data_df the main dataset used by the clustering algorithm
#' @param outlier_df any outlier observations that you do not want to influence 
#' the clustering algorithm. The clusters for these observations will be 
#' 'predicted' once the cluster algorithm has been run.
#' 
#'
#' @return a dataframe including observation from both `data_df` and 
#' `outlier_df` including a cluster assignment for each observation

get_clustering_results <- function(k, recipe, data_df, outlier_df) {
  # ---------------------------------------------------------------------------
  # set up clustering algorithm
  # ---------------------------------------------------------------------------
  kmeans_spec <- k_means(num_clusters = k) %>%
    set_engine("stats")

  # ---------------------------------------------------------------------------
  # set up workflow which brings together the clustering algorithm and
  # the data pre-processing (i.e. the recipe) for execution
  # ---------------------------------------------------------------------------
  kmeans_wflow <- workflow(recipe, kmeans_spec)

  # ---------------------------------------------------------------------------
  # Main data used for determining how observations should be classified
  # ---------------------------------------------------------------------------
  kmeans_fit <- kmeans_wflow %>%
    fit(data = select(data_df, -msoa_code))

  cluster_assignments <- extract_cluster_assignment(kmeans_fit)

  main <- data_df %>%
    bind_cols(cluster_assignments)

  # ---------------------------------------------------------------------------
  # Outliers data not involved in determining how observations should be
  # classified
  # ---------------------------------------------------------------------------
  outliers <- kmeans_fit %>%
    augment(new_data = outlier_df) %>%
    rename(.cluster = .pred_cluster)

  # ---------------------------------------------------------------------------
  # return all the data together including cluster classifications
  # ---------------------------------------------------------------------------
  bind_rows(main, outliers)
}

# a quick test things are working as expected
clustering_res <- get_clustering_results(
  k = 6,
  recipe = pca_recipe_2,
  data_df = foe_green_space_focus,
  outlier_df = foe_green_space_outliers
)

# confirm the returned dataframe has the expected number of rows
assertthat::assert_that(nrow(clustering_res) == nrow(foe_green_space_focus) +
  nrow(foe_green_space_outliers))

#' View some visualisations of FoE green space clustering results
#'
#' @param df a dataframe of observations including a `.cluster` variable
#' @param k an integer the number of clusters identified 
#' (for inclusion in plot titles)
#' @param recipe_name a string naming the recipe used to preprocess the clustering
#' input data (for inclusion in plot titles)
#'
#' @return a list of three ggplotly objects (interactive ggplot charts created
#' using plotly::ggplotly(p)) showing clustering results.

view_cluster_dist <- function(df, k, recipe_name) {
  
  #-----------------------------------------------------------------------------
  # set up common features used across multiple plots
  #-----------------------------------------------------------------------------
  
  common_plot_elements <- function(facet = TRUE) {
    
    background_grey <- "grey90"
    
    list(
      
      # colour palette
      scale_colour_viridis_d(),
      
      # facet styling
      if(facet) 
        facet_wrap(~variable, scales = "free", ncol = 2),
      theme(
        strip.background = element_rect(fill = background_grey),
        strip.text = element_text(colour = "black"),
        
        # add a neutral background (rather than stark white)
        plot.background = element_rect(fill = background_grey),
        panel.background = element_rect(fill = background_grey),
        legend.background = element_blank()
      ),
      
      # plot labeling
      ggtitle(glue::glue("{k} clusters. Preprocessing: {recipe_name}"))
    )
    
  }

  # ---------------------------------------------------------------------------
  # process data for plotting
  # ---------------------------------------------------------------------------
  
  # put data in long format - required for facetted plots
  plot_df <- df %>%
    select(-c(msoa_code)) %>%
    pivot_longer(
      cols = population_density:pcnt_bame_pop,
      names_to = "variable"
    )

  # split data for use in different plots
  
  standard_scale_vars <-  c(
      "income_index",
      "pcnt_pop_without_go_space_access"
    )
  
  # plots with a standard x axis scale
  standard_plot <- plot_df %>%
    filter(variable %in% standard_scale_vars)

  # plots with a log x axis scale
  log_plot <- plot_df %>%
    filter(!(variable %in% standard_scale_vars) & variable != "green_space_area_per_capita")
  
  gsapc_plot <- plot_df %>% # more customisation is need to get good view of the
                            # distribution of this variable
    filter(variable == "green_space_area_per_capita")

  #-----------------------------------------------------------------------------
  # produce plots
  #-----------------------------------------------------------------------------
  
  # on standard x axis scale
  p1 <- standard_plot %>%
    ggplot(aes(value, colour = .cluster, group = .cluster)) +
    geom_density() +
    common_plot_elements() +
    labs(x = NULL)


  # on log x axis scale
  p2 <- log_plot %>%
    # + 1e-10 an offset to avoid log(0)
    ggplot(aes(value + 1e-10, colour = .cluster)) +
    geom_density() +
    scale_x_log10() +
    common_plot_elements() +
    labs(x = NULL)
  
  # more customisation is need to get good view of the
  # distribution of this variable
  p3 <- gsapc_plot %>% 
    # + 1e-10 an offset to avoid log(0)
    ggplot(aes(value + 1e-10, colour = .cluster)) +
    geom_density() +
    scale_x_log10(limits = c(1, 1000)) +
    common_plot_elements(facet = FALSE) +
    labs(x = "log green_space_area_per_capita")

  # return plots
  list(p1, 
       p2,
       p3)
}

# set number of clusters as used in two functions below
k <- 6

res <- get_clustering_results(
  k = k,
  recipe = pca_recipe_2,
  data_df = foe_green_space_focus,
  outlier_df = foe_green_space_outliers
) 

# output for inspection
kable(head(res))

# view results charts
res %>%
  view_cluster_dist(k = k, 
                    # get recipe_name as string
                    recipe = deparse(substitute(pca_recipe_2)))

# output for use in other notebooks
readr::write_csv(res, "../data_out/clusters_cv_5_k_6_pca_recipe_2.csv")