The ecological context

The Irish hare (Lepus timidus hibernicus) is a subspecies endemic to Ireland — found nowhere else on Earth. It has been present since before the last Ice Age, making it one of the island’s oldest residents and a genuine conservation priority.

The European (brown) hare (Lepus europaeus), by contrast, is an invasive species in Northern Ireland, introduced sometime in the 19th century, probably, like so many ecological misfortunes, by people who thought it was a good idea at the time. L. europaeus poses a threat to L. t. hibernicus through competition and hybridisation, and its management has been actively debated.

The question is not only can we cull brown hares, but should we and, critically, will people support it? That last consideration is where social science enters the room, takes off its coat, and starts asking uncomfortable questions.

Learning objectives

By the end of this tutorial you will be able to:

1. Load, clean, and recode a real-world survey dataset using tidyverse
2. Compute descriptive statistics
3. Apply chi-square tests and Fisher’s exact test with appropriate multiple-comparison corrections
4. Compute pairwise phi coefficients and visualise response linkages as networks
5. Fit and interpret a binary logistic regression model
6. Perform quantitative thematic analysis of multi-select responses
7. Identify respondent classes and interpret them in meaningful ways

The study

This tutorial uses data from:

Caravaggi A, Montgomery WI, Reid N (2017). Management and control of invasive brown hares (Lepus europaeus): contrasting attitudes of selected environmental stakeholders and the wider rural community. Biology and Environment: Proceedings of the Royal Irish Academy, 117B, 53–63. https://doi.org/10.3318/bioe.2017.08

The study surveyed two groups:

• CAI: Members of Countryside Alliance Ireland — an organisation representing rural and field sports interests
• Public: The wider rural community (non-members)

Packages

# Install any missing packages by uncommenting the line below
# install.packages(c("tidyverse","vcd","ca","igraph","ggraph",
#                    "corrplot","broom","scales","patchwork","janitor","ggrepel", "polCA"))

library(vcd)         # Visualising categorical data (mosaic plots, association stats)
library(ca)          # Correspondence analysis
library(igraph)      # Network/graph objects
library(ggraph)      # ggplot2-style network visualisation
library(corrplot)    # Correlation/association matrix heatmaps
library(broom)       # Converts messy model output into tidy data frames
library(scales)      # Axis label formatting helpers
library(patchwork)   # Combine multiple ggplots into one figure
library(janitor)     # Data cleaning utilities (clean_names, tabyl, etc.)
library(ggrepel)     # Non-overlapping text labels on plots
library(poLCA)       # Latent Class Analysis
library(tidyverse)   # The Swiss Army knife: dplyr, ggplot2, tidyr, readr, etc.

set.seed(42)  # Reproducibility: ensures any randomness gives the same result
              # 42 chosen for reasons that will be obvious to anyone who has
              # read the right books.

Loading the data

Download https://zenodo.org/records/891251 and place it in your working directory. The dataset is published under a CC BY 4.0 licence, meaning you can use it freely as long as you cite it.

# Load the raw data
# Adjust the path if your file is in a subfolder, e.g. "data/Hare_data.csv"
hare_raw <- read.csv("data.csv",
                     stringsAsFactors = FALSE,
                     na.strings        = c("", "NA"))

# clean_names() from {janitor} converts column names to consistent snake_case:
# "Own/rent land?" becomes "own_rent_land", "L. europaeus" becomes "l_europaeus", etc.
# This saves enormous amounts of grief when typing column names.
hare_raw <- clean_names(hare_raw)

dim(hare_raw)

## [1] 341  44

Awareness by geographic zone

One of the central ecological questions in this study is whether proximity to established L. europaeus populations affects awareness. Respondents from the CORE zone live alongside the invasive species; those in WIDER NI may never have encountered one. Do their perceptions differ?

aware_labels <- c(
  seen_hares                           = "Seen Irish hares",
  seen_l_europaeus                     = "Seen European hares",
  aware_of_2_sp_in_ni                  = "Aware of 2 species in NI",
  is_threat_of_l_europaeus_important   = "Considers European hare threat important"
)

hare |>
  dplyr::select(zone, all_of(aware_cols)) |>
  pivot_longer(-zone) |>
  group_by(zone, name) |>
  summarise(pct = mean(value, na.rm = TRUE) * 100, .groups = "drop") |>
  mutate(name = factor(name, levels = names(aware_labels), labels = aware_labels)) |>
  ggplot(aes(x = zone, y = pct, fill = name)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 50, linetype = "dashed", colour = "grey50", linewidth = 0.5) +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = percent_format(scale = 1), limits = c(0, 105)) +
  labs(title    = "Hare awareness by geographic zone",
       subtitle = "Dashed line = 50%",
       x = "Zone", y = "% respondents answering Yes", fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 9))

Conservation concern: a mosaic plot

A mosaic plot is the multivariate extension of a stacked bar chart. The width of each column represents the proportion of respondents in that group; the height of each cell within the column represents the proportion choosing each response. Colour shading (Pearson residuals) shows which cells deviate from what you would expect if group membership and conservation concern were independent.

# Compute residuals
tbl <- table(Group   = hare$group,
             Concern = hare$conservation_concern) |>
  (\(t) t[rowSums(t) > 0, colSums(t) > 0])()

res_df <- chisq.test(tbl)$residuals |>
  as.data.frame() |>
  rename(group = Group, concern = Concern, residual = Freq)

# Compute proportions for tile sizing
prop_df <- as.data.frame(prop.table(tbl, margin = 1)) |>
  rename(group = Group, concern = Concern, prop = Freq)

plot_df <- left_join(res_df, prop_df, by = c("group", "concern"))

ggplot(plot_df, aes(x = concern, y = group, fill = residual)) +
  geom_tile(colour = "white", linewidth = 3) +
  geom_text(aes(label = round(residual, 2)),
            size = 4.5, fontface = "bold",
            colour = "white") +
  scale_fill_gradient2(
    low      = "#d73027",
    mid      = "white",
    high     = "#2166ac",
    midpoint = 0,
    name     = "Pearson\nresidual"
  ) +
  labs(
    title    = "Conservation concern by respondent group",
    subtitle = "Cell colour = Pearson residual.Blue = more than expected. Red = fewer",
    x        = "Conservation concern",
    y        = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(panel.grid = element_blank())

The shading uses Pearson residuals — the difference between observed and expected cell counts, standardised by the square root of the expected count. A residual greater than ±2 is conventionally considered noteworthy. Which cells, if any, show strong deviations? What might explain them?

Management support: bar chart and heatmap

We want to see not just overall support levels, but how much the two groups differ. A bar chart gives the absolute picture; a heatmap lets the eye quickly spot patterns across many variables simultaneously.

mgmt_labels <- c(
  support_management_of_l_europaeus       = "General management\nsupport",
  habitat_management_for_l_t_hibernicus   = "Habitat management\n(non-lethal)",
  lethal_culling                          = "Lethal culling (any)",
  support_management_via_petition         = "Sign a petition",
  support_cull_with_direct_involvement    = "Direct involvement\nin cull",
  allow_cull_on_land                      = "Allow cull on land",
  support_gov_t_cull_decision             = "Support\nGovernment cull"
)

plot_df <- hare |>
  dplyr::select(group, all_of(mgmt_cols)) |>
  pivot_longer(-group) |>
  group_by(group, name) |>
  summarise(pct = mean(value, na.rm = TRUE) * 100, .groups = "drop") |>
  mutate(name = factor(name, levels = names(mgmt_labels), labels = mgmt_labels))

p_bar <- ggplot(plot_df, aes(x = pct, y = fct_rev(name), fill = group)) +
  geom_col(position = position_dodge(0.7), width = 0.6) +
  geom_text(aes(label = paste0(round(pct), "%")),
            position = position_dodge(0.7), hjust = -0.1, size = 3) +
  scale_fill_manual(values = c("CAI" = "#2166ac", "Public" = "#d73027")) +
  scale_x_continuous(limits = c(0, 115), labels = percent_format(scale = 1)) +
  labs(title    = "Support for European hare management strategies",
       subtitle = "By respondent group",
       x = "% respondents", y = NULL, fill = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top")

print(p_bar)

ggplot(plot_df, aes(x = name, y = group, fill = pct)) +
  geom_tile(colour = "white", linewidth = 0.8) +
  geom_text(aes(label = paste0(round(pct), "%")),
            size = 3.5, colour = "grey20") +
  scale_fill_gradient2(
    low      = "#d73027",
    mid      = "white",
    high     = "#2166ac",
    midpoint = 50,
    limits   = c(0, 100),
    name     = "% Yes"
  ) +
  labs(title = "Management support patterns",
       subtitle = "Heatmap view — each cell shows % respondents answering Yes",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 40, hjust = 1, size = 9))

The bar chart is better for reading exact values and comparing groups side-by-side. The heatmap is better for spotting patterns across many variables at once. Neither is objectively superior — it depends on what question you are trying to answer. In a real report, you would probably use one or the other, not both.

Theory

The chi-square test asks a simple question: is there a statistically significant association between two categorical variables? It does so by comparing the counts of actual observations to the counts you would expect if the variables were completely independent.

When any expected cell count falls below 5, the chi approximation becomes unreliable. Fisher’s exact test calculates the exact probability rather than an approximation so is always valid.

A chi-square test does not tell Whether the association is large or small, however. A statistically significant result only means that the effect is probably not zero and not that it is meaningful. To explore meaning we need a measure of effect size.

Cramér’s V: measuring effect size

Cramér’s V scales the \(\chi^2\) statistic to a 0–1 range, correcting for sample size and table dimensions. Typically we interpret V as ~0.1 = small, ~0.3 = medium, ~0.5 = large.

Applying to our data

We test each management variable against group membership (CAI vs Public). In plain terms we’re asking whether each group answers differently on this issue. We’re asking this of 10 variables and if you test things many times then the likelihood of false positives increases considerably. To mitigate this we also apply a Bonferroni correction that divides our usual cutoff for ‘significant’ by the number of tests. We have 10 tests so our P threshold is now 0.005.

# A function that runs the appropriate test for a single variable
compare_groups <- function(var, data = hare) {
  tbl    <- table(data$group, data[[var]])
  # Suppress warnings about small expected values (we check manually)
  exp_ct <- suppressWarnings(chisq.test(tbl)$expected)

  if (any(exp_ct < 5, na.rm = TRUE)) {
    test   <- fisher.test(tbl)
    method <- "Fisher"
    stat   <- NA_real_
    p      <- test$p.value
    OR     <- as.numeric(test$estimate)
  } else {
    test   <- suppressWarnings(chisq.test(tbl))
    method <- "Chi-sq"
    stat   <- as.numeric(test$statistic)
    p      <- test$p.value
    # Cramér's V for chi-square
    n_obs  <- sum(tbl)
    OR     <- sqrt(stat / (n_obs * (min(dim(tbl)) - 1)))  # repurpose column for V
  }

  tibble(
    variable  = var,
    method    = method,
    statistic = round(stat, 3),
    p_value   = round(p, 4),
    cramers_V_or_OR = round(OR, 3),
    effect_note = if_else(method == "Chi-sq", "Cramer's V", "Odds ratio"),
    sig = case_when(p < 0.001 ~ "***",
                    p < 0.01  ~ "**",
                    p < 0.05  ~ "*",
                    TRUE       ~ "")
  )
}

comparison_table <- map_dfr(mgmt_cols, compare_groups) |>
  arrange(p_value) |>
  mutate(
    p_bonferroni = p.adjust(p_value, method = "bonferroni"),
    sig_bonf     = case_when(p_bonferroni < 0.001 ~ "***",
                              p_bonferroni < 0.01  ~ "**",
                              p_bonferroni < 0.05  ~ "*",
                              TRUE                  ~ "ns")
  )

comparison_table |>
  dplyr::select(variable, method, statistic, p_value, cramers_V_or_OR,
         effect_note, sig, sig_bonf) |>
  knitr::kable(
    caption = "Table 2. Statistical comparisons of management support between CAI members and Public respondents. Sig. = before Bonferroni correction; Sig. (Bonf.) = after correction.",
    col.names = c("Variable","Test","Statistic","p-value",
                  "Effect size","Measure","Sig.","Sig. (Bonf.)")
  )

Does species awareness predict management support?

An ecologically interesting sub-question: do respondents who know there are two hare species in Northern Ireland show different levels of support for management?

hare |>
  filter(!is.na(aware_of_2_sp_in_ni),
         !is.na(support_management_of_l_europaeus)) |>
  count(aware_of_2_sp_in_ni, support_management_of_l_europaeus) |>
  mutate(
    aware_of_2_sp_in_ni              = if_else(aware_of_2_sp_in_ni, "Aware", "Not aware"),
    support_management_of_l_europaeus = if_else(support_management_of_l_europaeus,
                                                 "Supports", "Opposes")
  ) |>
  knitr::kable(caption = "Table 3. Species awareness vs management support crosstab.")

fisher.test(table(hare$aware_of_2_sp_in_ni,
                  hare$support_management_of_l_europaeus))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(hare$aware_of_2_sp_in_ni, hare$support_management_of_l_europaeus)
## p-value = 0.5904
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.4839208 1.4835166
## sample estimates:
## odds ratio 
##  0.8457105

If people who are more aware of the species distinction are less likely to support management, what might this suggest? Could ecological knowledge cut both ways - increasing concern for the Irish hare but also increasing reluctance to sanction lethal control of the other? This is a genuinely interesting question that has real implications for how conservation communication should be framed.

Theory: phi coefficients for binary data

So far we have been looking at one variable at a time, or comparing one variable against group membership. But survey responses are not independent of each other. Someone who says they would sign a petition may also be likely to allow a cull on their land. Mapping these internal associations reveals the structure of opinion.

For two binary (Y/N) variables, the appropriate correlation coefficient is phi (\(\phi\)), which is conceptually identical to Pearson’s \(r\) applied to 0/1-coded data.

Computing the phi matrix

all_binary_cols <- c(aware_cols, mgmt_cols, "member_of_conservation_org", activity_cols)

binary_mat <- hare |>
  dplyr::select(all_of(all_binary_cols)) |>
  mutate(across(everything(), as.integer))

phi_mat <- cor(binary_mat, use = "pairwise.complete.obs")

## Warning in cor(binary_mat, use = "pairwise.complete.obs"): the standard
## deviation is zero

phi_mat[is.na(phi_mat)] <- 0  # replace NA correlations with 0

# Readable labels
nice_names <- c(
  seen_hares                           = "Seen Irish hare",
  seen_l_europaeus                     = "Seen Eur. hare",
  aware_of_2_sp_in_ni                  = "Aware 2 spp",
  is_threat_of_l_europaeus_important   = "Threat important",
  support_management_of_l_europaeus    = "Support mgmt",
  habitat_management_for_l_t_hibernicus= "Habitat mgmt",
  lethal_culling                       = "Lethal cull",
  support_management_via_petition      = "Petition",
  support_cull_with_direct_involvement = "Direct involve",
  allow_cull_on_land                   = "Allow cull",
  support_gov_t_cull_decision          = "Govt cull",
  member_of_conservation_org           = "Consv org",
  hunt                                 = "Hunt",
  game_birds                           = "Game birds",
  rabbits                              = "Rabbits",
  wildfowl                             = "Wildfowl",
  fishing_angling                      = "Fishing"
)
rownames(phi_mat) <- colnames(phi_mat) <- nice_names[colnames(phi_mat)]

Heatmap of phi coefficients

# Make the plot device tall and wide before calling corrplot
par(oma = c(0, 0, 2, 0))  # outer margin for title

corrplot(
  phi_mat,
  method       = "color",
  type         = "lower",          # full matrix reads better at this size
  order        = "hclust",        # group correlated variables together
  hclust.method = "ward.D2",      # matches your clustering analysis
  tl.cex       = 0.85,            # larger labels
  tl.col       = "black",
  tl.srt       = 45,              # angled top labels, less cramped
  addCoef.col  = "black",
  number.cex   = 0.65,
  number.digits = 1,              # 1 decimal place keeps cells readable
  col          = colorRampPalette(c("#d73027", "white", "#2166ac"))(200),
  mar          = c(0, 0, 1, 0),
  diag         = FALSE            # hide the diagonal 1s
)
title("Phi coefficients: binary response associations",
      outer = TRUE, cex.main = 1.1)

Hierarchical clustering: what clusters with what?

Hierarchical clustering groups variables by similarity. Here we treat \(1 - |\phi|\) as a distance (so highly associated variables, whether positively or negatively, are close together), and use Ward’s D2 linkage to merge groups.

dist_mat   <- as.dist(1 - abs(phi_mat))
hclust_res <- hclust(dist_mat, method = "ward.D2")

plot(hclust_res,
     main  = "Hierarchical clustering of questionnaire responses",
     sub   = "Distance = 1 - |φ|  (Ward's D2 linkage)",
     xlab  = "Response variable",
     ylab  = "Distance",
     cex   = 0.85)
rect.hclust(hclust_res, k = 4,
            border = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3"))

Variables that merge at a low height are closely associated. Variables that only merge near the top of the tree are weakly associated. The coloured rectangles suggest a four-cluster partition but selecting \(k\) is a judgement call, not an algorithm output. Try k = 3 and k = 5 and see how the interpretation changes.

Association network

We can go one step further and visualise the associations as a network where variables are nodes, and we draw an edge between two nodes if their \(|\phi| > 0.3\) (a conventional threshold for a “meaningful” binary association). Edge thickness and colour encode the strength and direction of the association.

phi_long <- as_tibble(phi_mat, rownames = "from") |>
  pivot_longer(-from, names_to = "to", values_to = "phi") |>
  filter(from < to,        # lower triangle only — avoid duplicate edges
         abs(phi) > 0.3,
         !is.na(phi))

if (nrow(phi_long) > 0) {
  g <- graph_from_data_frame(
    d        = phi_long |> mutate(weight = abs(phi)),
    directed = FALSE,
    vertices = tibble(name = colnames(phi_mat))
  )
  E(g)$sign <- ifelse(phi_long$phi > 0, "positive", "negative")

  set.seed(42)
  ggraph(g, layout = "fr") +
    geom_edge_link(aes(width = weight, colour = sign), alpha = 0.8) +
    geom_node_point(size = 7, colour = "#444444", fill = "#dddddd",
                    shape = 21, stroke = 1) +
    geom_node_text(aes(label = name), size = 2.8, repel = TRUE,
                   max.overlaps = 20) +
    scale_edge_width(range = c(0.5, 3.5), name = "|φ|") +
    scale_edge_colour_manual(
      values = c(positive = "#2166ac", negative = "#d73027"),
      name   = "Direction"
    ) +
    labs(
      title    = "Association network: questionnaire responses",
      subtitle = "Edges where |φ| > 0.3 · blue = positive · red = negative"
    ) +
    theme_graph(base_family = "sans") +
    theme(legend.position = "right")
} else {
  cat("No pairs with |phi| > 0.3 in this sample. Try lowering the threshold.\n")
}

In a well-designed attitude survey, you often see two major clusters: one for knowledge/awareness items and one for action/support items. Do these emerge here? What does it mean if awareness and support cluster together — or separately?

Fitting the model

hare_mod <- hare |>
  mutate(
    group                = relevel(group, ref = "Public"),
    conservation_concern = relevel(conservation_concern, ref = "Not concerned")
  ) |>
  filter(!is.na(support_gov_t_cull_decision))

glm_fit <- glm(
  support_gov_t_cull_decision ~
    group + zone + own_rent_land + conservation_concern +
    aware_of_2_sp_in_ni + consider_hares_pests +
    is_threat_of_l_europaeus_important + hunt,
  data   = hare_mod,
  family = binomial(link = "logit")
)

summary(glm_fit)

## 
## Call:
## glm(formula = support_gov_t_cull_decision ~ group + zone + own_rent_land + 
##     conservation_concern + aware_of_2_sp_in_ni + consider_hares_pests + 
##     is_threat_of_l_europaeus_important + hunt, family = binomial(link = "logit"), 
##     data = hare_mod)
## 
## Coefficients:
##                                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                            -0.25476    0.51240  -0.497   0.6191  
## groupCAI                                0.06835    0.46442   0.147   0.8830  
## zonePERIPHERY                           0.24281    0.37854   0.641   0.5212  
## zoneWIDER NI                            0.09899    0.38234   0.259   0.7957  
## own_rent_landTRUE                      -0.28734    0.35088  -0.819   0.4128  
## conservation_concernConcerned           0.25992    0.39922   0.651   0.5150  
## aware_of_2_sp_in_niTRUE                -0.76739    0.38291  -2.004   0.0451 *
## consider_hares_pestsTRUE               -0.20057    0.41016  -0.489   0.6248  
## is_threat_of_l_europaeus_importantTRUE  0.82124    0.32678   2.513   0.0120 *
## huntTRUE                                1.20848    0.48102   2.512   0.0120 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 262.50  on 199  degrees of freedom
## Residual deviance: 241.55  on 190  degrees of freedom
##   (121 observations deleted due to missingness)
## AIC: 261.55
## 
## Number of Fisher Scoring iterations: 4

Odds ratios and confidence intervals

We calculate odds ratios in logistic regression to express how much the odds of an outcome change when a predictor increases or differs between groups. They make results easier to interpret by turning model coefficients into a clear “how many times more or less likely” comparison.

glm_tidy <- tidy(glm_fit, conf.int = TRUE, exponentiate = TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(
    sig = case_when(p.value < 0.001 ~ "***",
                    p.value < 0.01  ~ "**",
                    p.value < 0.05  ~ "*",
                    TRUE             ~ "")
  )

glm_tidy |>
  dplyr::select(term, estimate, conf.low, conf.high, p.value, sig) |>
  knitr::kable(
    digits  = 3,
    caption = "Table 4. Logistic regression results. Odds ratios (ORs) with 95% confidence intervals and p-values. Reference group: Public respondents; Not concerned about conservation.",
    col.names = c("Predictor", "OR", "95% CI lower", "95% CI upper", "p-value", "Sig.")
  )

ggplot(glm_tidy, aes(x = estimate, y = reorder(term, estimate))) +
  geom_point(size = 3.5, colour = "#2166ac") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 height = 0.3, colour = "#2166ac", linewidth = 0.8) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "#e74c3c",
             linewidth = 0.8) +
  scale_x_log10() +
  labs(
    title    = "Predictors of government cull support",
    subtitle = "Odds ratios with 95% CI (log scale)",
    x = "Odds ratio", y = NULL
  ) +
  theme_minimal(base_size = 12)

## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `height` was translated to `width`.

Model fit

cat("AIC:", round(AIC(glm_fit), 1), "\n")

## AIC: 261.6

cat("Null deviance:", round(glm_fit$null.deviance, 1),
    "on", glm_fit$df.null, "df\n")

## Null deviance: 262.5 on 199 df

cat("Residual deviance:", round(glm_fit$deviance, 1),
    "on", glm_fit$df.residual, "df\n")

## Residual deviance: 241.6 on 190 df

cat("McFadden pseudo-R:",
    round(1 - glm_fit$deviance / glm_fit$null.deviance, 3), "\n")

## McFadden pseudo-R: 0.08

Unlike the usual R^2 in linear regression, McFadden’s pseudo-R^2 does not tell you how much of the outcome is “explained” by the model. Instead, it tells you how much better your model fits the data compared to a baseline model with only an intercept. Values around 0.2–0.4 are often considered quite good in social science, but it should be interpreted cautiously rather than taken as a strict measure of model quality.

Which predictors have ORs substantially different from 1.0? Are the confidence intervals narrow (precise estimates) or wide (high uncertainty)? How might you improve the model? What would an interaction between group and aware_of_2_sp_in_ni mean conceptually?

Theory

Some questions in the survey were multi-select, i.e., respondents could tick all that applied. These cannot be treated as simple binary variables as the combination of choices carries meaning.

In this section we’ll use a qualitative thematic approach. We identify the “themes” (choice profiles) present in the data and examine which groups favour which themes. This allows us to be fully quantitative and use pattern analysis on binary data.

Culling method profiles

cull_mat <- hare |>
  dplyr::select(group, all_of(cull_methods)) |>
  filter(!is.na(shooting) | !is.na(netting) | !is.na(trapping)) |>
  mutate(across(all_of(cull_methods), ~replace_na(as.integer(.), 0L))) |>
  mutate(
    n_methods = shooting + netting + trapping,
    profile   = case_when(
      n_methods == 0              ~ "None selected",
      shooting & !netting & !trapping ~ "Shooting only",
      !shooting & netting & !trapping ~ "Netting only",
      !shooting & !netting & trapping ~ "Trapping only",
      shooting & netting & !trapping  ~ "Shooting + Netting",
      shooting & !netting & trapping  ~ "Shooting + Trapping",
      !shooting & netting & trapping  ~ "Netting + Trapping",
      TRUE                            ~ "All three"
    )
  )

count(cull_mat, group, profile) |> arrange(group, desc(n))

##     group             profile  n
## 1     CAI       Shooting only 34
## 2     CAI           All three 29
## 3     CAI        Netting only 21
## 4     CAI  Shooting + Netting  7
## 5     CAI  Netting + Trapping  5
## 6     CAI       Trapping only  4
## 7     CAI Shooting + Trapping  2
## 8  Public       Shooting only 33
## 9  Public        Netting only 27
## 10 Public           All three 13
## 11 Public  Netting + Trapping 10
## 12 Public       Trapping only  9
## 13 Public  Shooting + Netting  2

cull_mat |>
  count(group, profile) |>
  group_by(group) |>
  mutate(pct = n / sum(n) * 100) |>
  ggplot(aes(x = pct, y = reorder(profile, pct), fill = group)) +
  geom_col(position = position_dodge(0.7), width = 0.65) +
  geom_text(aes(label = paste0(round(pct), "%")),
            position = position_dodge(0.7), hjust = -0.15, size = 3) +
  scale_fill_manual(values = c("CAI" = "#2166ac", "Public" = "#d73027")) +
  scale_x_continuous(limits = c(0, 75), labels = percent_format(scale = 1)) +
  labs(title    = "Preferred culling method combinations",
       subtitle = "Multi-select question: tick all that apply",
       x = "% respondents", y = NULL, fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

Linking themes to attitudes

We can use chi-square to quantitatively test whether culling method preference differs between groups.

profile_tbl <- table(cull_mat$group, cull_mat$profile)
chisq.test(profile_tbl)

## Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
## incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  profile_tbl
## X-squared = 14.926, df = 6, p-value = 0.02084

Activity co-occurrence: the Jaccard index

For the rural activities question, we want to know which activities tend to co-occur, i.e., which activities are commonly practised together. We use the Jaccard similarity index, which measures the overlap between two sets. A Jaccard of 1 means the two activities always occur together; 0 means they never do.

act_mat <- hare |>
  dplyr::select(all_of(activity_cols)) |>
  mutate(across(everything(), ~replace_na(as.integer(.), 0L))) |>
  as.matrix()

co_occur <- t(act_mat) %*% act_mat   # n(i AND j) for all pairs
diag(co_occur) <- 0

n_each      <- colSums(act_mat)
jaccard_mat <- matrix(NA_real_, nrow = ncol(act_mat), ncol = ncol(act_mat),
                      dimnames = list(activity_cols, activity_cols))

for (i in seq_len(ncol(act_mat))) {
  for (j in seq_len(ncol(act_mat))) {
    union_ij          <- n_each[i] + n_each[j] - co_occur[i, j]
    jaccard_mat[i, j] <- if (union_ij > 0) co_occur[i, j] / union_ij else 0
  }
}

round(jaccard_mat, 2)

##                 hunt game_birds rabbits wildfowl fishing_angling
## hunt            0.00       0.66    0.54     0.54            0.44
## game_birds      0.66       0.00    0.56     0.71            0.53
## rabbits         0.54       0.56    0.00     0.54            0.38
## wildfowl        0.54       0.71    0.54     0.00            0.53
## fishing_angling 0.44       0.53    0.38     0.53            0.00

g_act <- graph_from_adjacency_matrix(
  jaccard_mat, mode = "undirected", weighted = TRUE, diag = FALSE
)

set.seed(7)
ggraph(g_act, layout = "circle") +
  geom_edge_link(aes(width = weight, alpha = weight), colour = "#27ae60") +
  geom_node_point(size = 14, colour = "#1a7837", fill = "#a8ddb5",
                  shape = 21, stroke = 1.5) +
  geom_node_text(aes(label = name), size = 3.2, colour = "black") +
  scale_edge_width(range = c(0.3, 5), name = "Jaccard") +
  scale_edge_alpha(range = c(0.2, 1), guide = "none") +
  labs(title    = "Co-occurrence of rural activities",
       subtitle = "Edge width proportional to Jaccard similarity") +
  theme_graph(base_family = "sans")

Respondent typology

Finally, we create a simple classification by combining two measures: awareness (how much a respondent knows) and support (how strongly they support active management). This is similar to stakeholder mapping, where people are grouped based on their knowledge and attitudes to help guide conservation decisions.

hare <- hare |>
  mutate(
    awareness_index = as.integer(seen_hares) +
                      as.integer(aware_of_2_sp_in_ni) +
                      as.integer(is_threat_of_l_europaeus_important),

    support_index   = as.integer(support_management_of_l_europaeus) +
                      as.integer(lethal_culling) +
                      as.integer(support_management_via_petition) +
                      as.integer(support_cull_with_direct_involvement) +
                      as.integer(allow_cull_on_land) +
                      as.integer(support_gov_t_cull_decision),

    activity_breadth = as.integer(hunt) + as.integer(game_birds) +
                       as.integer(rabbits) + as.integer(wildfowl) +
                       as.integer(fishing_angling),

    typology = case_when(
      support_index >= 4 & awareness_index >= 2 ~ "Informed advocate",
      support_index >= 4 & awareness_index <  2 ~ "Pragmatic advocate",
      support_index <  2 & awareness_index >= 2 ~ "Aware but opposed",
      support_index <  1 & awareness_index <  1 ~ "Disengaged",
      TRUE                                      ~ "Ambivalent"
    ),
    typology = factor(typology,
                      levels = c("Informed advocate", "Pragmatic advocate",
                                 "Ambivalent", "Aware but opposed", "Disengaged"))
  )

count(hare, group, typology) |> arrange(group, typology)

##    group           typology   n
## 1    CAI  Informed advocate  27
## 2    CAI Pragmatic advocate   6
## 3    CAI         Ambivalent  83
## 4    CAI  Aware but opposed  23
## 5 Public  Informed advocate  32
## 6 Public Pragmatic advocate   8
## 7 Public         Ambivalent 130
## 8 Public  Aware but opposed  31
## 9 Public         Disengaged   1

It’s important to note that the thresholds used to define these categories are based on subjective judgement rather than fixed rules. This means that different reasonable cut-offs could change how individuals are grouped, so the typology should be interpreted as a useful simplification rather than a precise classification. Where possible, these thresholds should be checked against empirical evidence or validated with expert judgement to ensure they are meaningful and robust.

hare |>
  filter(!is.na(typology)) |>
  count(group, typology) |>
  group_by(group) |>
  mutate(pct = n / sum(n) * 100) |>
  ggplot(aes(x = pct, y = fct_rev(typology), fill = group)) +
  geom_col(position = position_dodge(0.7), width = 0.65) +
  geom_text(aes(label = paste0(round(pct), "%")),
            position = position_dodge(0.7), hjust = -0.1, size = 3.2) +
  scale_fill_manual(values = c("CAI" = "#2166ac", "Public" = "#d73027")) +
  scale_x_continuous(limits = c(0, 75), labels = percent_format(scale = 1)) +
  labs(
    title    = "Respondent typology by group",
    subtitle = "Based on awareness index (0–4) and support index (0–6)",
    x = "% respondents", y = NULL, fill = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

How would you use these typologies in conservation practice? The “informed advocate” already knows and supports so do you need to spend resources on them? The “pragmatic advocate” supports but doesn’t really know why; are they a potential liability if the evidence changes? The “aware but opposed” understand the issue and still say no. What would it take to change their mind, and should you try?