Are natural green spaces less accessible to people from minority ethnic backgrounds?

# prepare data for plotting
plot_df <- pan_clean %>% 
  
  count(ethnicity) %>% 
  
  # to order groups within the plot
  mutate(ethnicity = fct_reorder(ethnicity, n))

# will be useful for ordering facets in a plot later on
ethnicity_levels <- rev(levels(plot_df$ethnicity))

# plot frequency distribution
ggplot(plot_df, aes(ethnicity, n)) +
  
  geom_col(fill = colours[1], width = 0.8) +
  geom_text(aes(label = n), nudge_y = 2000) +
  
  coord_flip() +
  
  labs(y = "Number of survey respodents",
       x = NULL)

# format strings for table output
levels(plot_df$ethnicity) <- rev(str_replace(ethnicity_levels, "\n", " "))

# look at proportion of survey sample in each ethnicity group
plot_df %>% 
  mutate(prop = round(n / sum(n), 3)) %>% 
  kable()

# -----------------------------------------------------------------------------
# compare weighted and unweighted counts
# -----------------------------------------------------------------------------

# unweighted and weighted group counts
ethnicity <- pan_clean %>% 
  count(ethnicity) %>% 
  mutate(weighted = FALSE)

ethnicity_wt <- pan_clean %>% 
  count(ethnicity, wt = weight_percent) %>% 
  mutate(weighted = TRUE)

# combine unweighted and weighted group counts for plotting
plot_df <- bind_rows(ethnicity, ethnicity_wt) 
  
# produce the plot
ggplot(plot_df, aes(ethnicity, n, 
                    fill = weighted)) +
  
  geom_col(position = "dodge", width = 0.6) +
  
  scale_fill_manual(values = colours) +
  
  coord_flip() +
  
  labs(x = NULL,
       y= "Number of survey respondents",
       fill = "With weigthing?") +
  
  facet_wrap(~factor(ethnicity, ethnicity_levels), scales = "free") +
  theme(axis.text.y = element_blank()) +
  style_facets()

# any_visits
any_visits <- pan_clean %>% 
  count(any_visits_14) %>% 
  mutate(weighted = FALSE,
         prop = round(n / sum(n), 3))

# plot frequency distribution
ggplot(any_visits, aes(fct_reorder(any_visits_14, n),
                    prop)) +
  
  geom_col(fill = colours[1], width = 0.6) +
  geom_text(aes(label = prop), nudge_y = 0.03) +
  
  coord_flip() +
  
  labs(y = "Proportion of survey respodents",
       x = NULL,
       subtitle = "Have you visited green space in the last 14 days?")

# show table including proportions
select(any_visits, -weighted) %>% 
  kable()

# ----------------------------------------------------------------------------
# compare weighted and unweighted counts
# ----------------------------------------------------------------------------
any_visits_wt <- pan_clean %>% 
  count(any_visits_14, wt = weight_percent) %>% 
  mutate(weighted = TRUE)

plot_df <- bind_rows(any_visits, any_visits_wt)

ggplot(plot_df, aes(any_visits_14, n, fill = weighted)) +
  geom_col(position = "dodge", width = 0.6) +
  
  scale_fill_manual(values = colours) +
  
  labs(x = NULL,
       y= "Number of survey respondents",
       fill = "With weigthing?") +
  
  coord_flip() +
  
  facet_wrap(~any_visits_14, scales = "free") +
  style_facets() +
  theme(axis.text.y = element_blank())

# check why data has been read in as character
unique(pan_clean$no_of_visits) # "Don't know" etc. values

# remove text based responses
pan_clean_no_of_visits <- pan_clean %>%
  filter(!(no_of_visits %in% c("Don’t Know", "Prefer not to say"))) %>% 
  mutate(no_of_visits = as.numeric(no_of_visits))

p <- pan_clean_no_of_visits %>% 
  # calculate frequencies
  count(no_of_visits) %>% 
  
  # produce plot
  ggplot(aes(no_of_visits, n)) +
    
    geom_col(fill = colours[1]) +
    
    labs(x = "Number of visits to green space in the last 14 days",
         subtitle = "Number of survey respondents",
         y = NULL) +
  
    ggforce::facet_zoom(x = no_of_visits < 20)

p

s <- skim(select(pan_clean_no_of_visits, 
                  no_of_visits))

s

# ---------------------------------------------------------------------------------
# Prepare the data
# ---------------------------------------------------------------------------------

diff_in_props_df <- pan_clean %>% 
  
  # remove ambiguous answers
  filter(ethnicity != "Don't know", ethnicity != "Prefer not to say") %>% 
  filter(any_visits_14 != "Don't know", any_visits_14 != "Prefer not to say") %>% 
  
  mutate(ethnicity = if_else(ethnicity == "White", ethnicity, "Minority Ethnic Background")) 

# ---------------------------------------------------------------------------------
# simulate a representative sample based on the survey weightings
# ---------------------------------------------------------------------------------

set.seed(1683763) # so we get the same sample each time the code block runs

# what would group proprtions look like if we applied the weighting factor
group_props_weighted <- diff_in_props_df %>% 
  count(ethnicity, any_visits_14,
        wt = weight_percent) %>% 
  mutate(prop = n / sum(n), 3)


num_resp <- nrow(diff_in_props_df)

sim_diff_in_props_df <- group_props_weighted %>% 
  
  # simulate respondents based on weighted group proportions
  slice_sample(weight_by = prop, n = num_resp, replace = TRUE) %>% 
  select(-c(n, prop))
  
  
sim_diff_in_props_counts <- sim_diff_in_props_df %>% 
  # calculated proportions based on simulation of a sample using weightings
  count(ethnicity, any_visits_14) %>% 
  mutate(prop_weighted = round(n / sum(n), 3)) %>% 
  rename(n_weighted = n)

# calculate actual proportions in the raw data for comaparision
acutal_diff_in_props_counts <- diff_in_props_df %>% 
  count(ethnicity, any_visits_14) %>% 
  mutate(prop = round(n / sum(n), 3))

acutal_diff_in_props_counts %>% 
  left_join(sim_diff_in_props_counts) %>% 
  kable()

sim_diff_in_props_df %>% 
  count(ethnicity, any_visits_14) %>% 
  group_by(ethnicity) %>% 
  mutate(prop_of_ethnicity_group = n / sum(n)) %>% 
  kable()

set.seed(16833) # so we get the same results from the permuation test each time

# calculate the observed difference in the proportions of people with and without 
# Minority Ethnic backgrounds who have visited green space in the last 14 days
d_hat <- sim_diff_in_props_df %>% 
  specify(any_visits_14 ~ ethnicity, success = "Any visits") %>%
  calculate(stat = "diff in props", order = c("Minority Ethnic Background", 
                                              "White"))
# generate the null distribution
null_dist <- sim_diff_in_props_df %>%
  specify(any_visits_14 ~ ethnicity, success = "Any visits") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "diff in props", order = c("Minority Ethnic Background", 
                                              "White"))

# annotation
ann_text<-data.frame(
  x = 4, y = 100,
  label="geeks for geeks"
)

label_text <- "Red line shows\ndifference in proportions\nin the actual survey sample"

# visualise the null distibution and the observed differnce
visualise(null_dist) +
  shade_p_value(obs_stat = d_hat, direction = "two-sided") +
  
  geom_label(x = -0.0185, y = 100, 
            label = label_text, 
            colour = "red") +
  
  labs(x = "Difference in proportions between groups",
       y = "Number of simulations")

null_dist %>%
  get_p_value(obs_stat = d_hat, direction = "two-sided") %>% 
  kable()

bootstrap_dist <- sim_diff_in_props_df %>%
  specify(any_visits_14 ~ ethnicity, success = "Any visits") %>%
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "diff in props", order = c("Minority Ethnic Background", 
                                              "White"))
# get confidence intervals
percentile_ci <- get_ci(bootstrap_dist, level = 0.99)

# output and visualise confidence intervals
percentile_ci %>% 
  mutate(observed_diff_in_proportions = d_hat$stat) %>% 
  relocate(observed_diff_in_proportions) %>% 
  kable()

visualize(bootstrap_dist) +
  shade_confidence_interval(endpoints = percentile_ci) +
  
  geom_label(x = -0.0185, y = 150, 
            label = "99% CI\nin blue", 
            colour = "grey20") +
  
  labs(x = "Difference in proportions between groups",
       y = "Number of simulations")

raw_means_comp <- pan_clean_no_of_visits %>% 
  
  # remove outliers
  # assume than visiting green space more than twice per day
  # i.e. 28 times in the last 14 days is an outlier
  filter(no_of_visits <= 28) %>% 
  
  # remove ambiguous survey responses
  filter(!(ethnicity %in% c("Don’t know", "Prefer not to say")))

skim_df <- raw_means_comp %>%
  select(ethnicity, no_of_visits) %>% 
  
  # for ouput table formatting
  mutate(ethnicity = str_replace(ethnicity, '\n', ' ')) %>%
  
  # for group-wise summary statistics
  group_by(ethnicity) %>% 
  skim()

# reorder columns for readability
raw_means <- skim_df %>% 
  select(skim_variable:numeric.sd) %>% 
  rename(mean = numeric.mean,
         sd = numeric.sd) %>% 
  relocate(mean, .after = ethnicity)

# output to notebook
kable(raw_means)

# create plot
raw_means_comp %>% 
  ggplot(aes(no_of_visits, ethnicity)) +
    geom_violin() +
  
    labs(y = NULL, 
         x = "Number of Visits to Green Space in the last 14 days",
         subtitle = "Probability densities")

# create model
model <- lm(no_of_visits ~ ethnicity, 
            data = raw_means_comp)

#calculate marginal means
marginal_means <- estimate_means(model, ci = 0.99)

# format marginal means for notebook output
marginal_means %>% 
  format_ethnicity_str() %>% 
  kable()

# format ethnicity strings for plotting
raw_means <- raw_means %>% 
  mutate(ethnicity = if_else(ethnicity == 'Any other ethnic group or background',
                               'Any other ethnic group\nor background',
                               ethnicity))


# create plot 
ggplot() +
  
  # raw means
  geom_point(data = raw_means,
             mapping = aes(mean, ethnicity),
             colour = "red",
             position = position_nudge(y = -0.1)) +
  
  # marginal means
  geom_pointrange(data = marginal_means,
                 mapping = aes(Mean, ethnicity,
                               xmin = CI_low, xmax = CI_high),
                 position = position_nudge(y = 0.1)) +
  
  # tidy up text
  labs(subtitle = "Red points are raw means.
Black point-lines show marginal means 
    plus confidence intervals",
       y = NULL,
      x = "Mean number of visits in the last 14 days")

contrasts <- estimate_contrasts(model, contrast = "ethnicity")

contrasts_full_res <- contrasts %>% 
  # format strings for table output
  mutate(across(where(is.character), str_replace, '\n', ' ')) %>% 
  
  mutate(across(where(is.numeric), 
                ~round(., 3))
         )

contrasts_full_res %>% 
  kable()

contrasts_full_res %>% 
  filter(p < 0.01) %>% 
  kable()

plot(filter(contrasts, p < 0.01), estimate_means(model)) +
  coord_flip() +
  labs(x = NULL,
       y = "Mean number of visits in the last 14 days")