Comparison of Recycling Systems in England

Does the collection system have a meaningful impact on recycling rate?

Author

Clare Gibson

Published

July 18, 2026

Introduction

This paper sets out to answer the question: Does the collection system have a meaningful impact on recycling rate?

Local authorities in England can elect to use one of three collection systems for household recycling:

  • Co-Mingled: all non-organic recyclable material is put into a single bin by the householder
  • Two Stream: non-organic recyclable material must be sorted into two bins by the householder (e.g. one bin for glass and metal, and a second bin for paper, card and plastic)
  • Multi-Stream: householders sort their non-organic recyclable material into separate bins for each type of material

Multi-Stream systems place a higher burden on the householder, so this paper asks whether that higher burden is associated with higher recycling rates, and whether any difference is large enough to matter in practice.

The word “meaningful” in the research question is important. A difference in recycling rates can be statistically detectable without being practically meaningful. This analysis therefore reports not only whether the systems differ, but by how much, and with how much confidence.

A note on causation: local authorities choose their own collection systems. Therefore, this is observational data and not a randomised experiment. Throughout this paper the language is deliberately one of association rather than cause.

In this analysis, we set the threshold for statistical significance at p < 0.05.

Setup

Load required packages

Show code 
library(here)
library(tidyverse)
library(knitr)
library(moments)
library(rstatix)
library(broom)
library(showtext)
library(factoextra)
library(cluster)
Show code 
source(here("R/theme_datatranslator.R"))
Show code 
font_add_google("Work Sans", "Work Sans")
font_add_google("Roboto Mono", "Roboto Mono")
showtext_auto()
showtext_opts(dpi = 96)

Import the recycling data

Show code 
recycling <- 
  read_csv(here("data/cln/data_recycling_system.csv")) |> 
  # reorder factor levels to assist with plotting
  mutate(
    collection_type = fct_reorder(collection_type, household_recycle_rate)
  )

head(recycling) |> kable()
lad21cd local_authority population imd_average_score ruc21_settlement_class rurality collection_type household_collected household_recycled household_recycle_rate population_density
E06000001 Hartlepool 92571 35.037 Urban 0.031 Co-Mingled 39552 12893 0.3259759 941.2302
E06000002 Middlesbrough 143734 40.460 Urban 0.002 Co-Mingled 64433 19169 0.2975028 2634.7305
E06000003 Redcar and Cleveland 136616 29.792 Intermediate urban 0.317 Two Stream 57866 22077 0.3815194 538.3132
E06000004 Stockton-on-Tees 197030 25.790 Urban 0.056 Multi-Stream 88050 22592 0.2565815 939.4423
E06000005 Darlington 108222 25.657 Urban 0.131 Multi-Stream 46203 14968 0.3239616 548.0212
E06000006 Halton 128577 32.325 Urban 0.007 Co-Mingled 59958 22137 0.3692084 1423.5478

The dataset includes information about household recycling for 309 local authorities in England in 2021. It contains the unique code and name of each authority, the population of the authority in 2021, the type of collection system in place in 2021, the amount of household waste collected and recycled in 2021 (in tonnes) and the household recycling rate.

Helper code

The chunk below contains code to support consistent plotting.

Show code 
# We'll be creating a number of similar plots, so we can wrap the code into a function
plot_distribution <- function(data, yvar, title, scale = c("none", "percent", "thousands")) {
  scale <- match.arg(scale)

  p <- ggplot(data, aes(x = collection_type, y = {{ yvar }})) +
    geom_violin(aes(fill = collection_type, colour = collection_type),
                linewidth = 1, alpha = 0.5) +
    geom_boxplot(outlier.alpha = 0, coef = 0,
                 colour = dt_navy, width = 0.2) +
    stat_summary(fun = mean, geom = "point", shape = 18, size = 3,
                 colour = dt_navy) +
    coord_flip() +
    scale_colour_dt(guide = "none") +
    scale_fill_dt(guide = "none") +
    labs(x = NULL, y = NULL, title = title) +
    theme_dt() +
    theme(
      panel.grid.major.y = element_blank(),
      panel.grid.minor   = element_blank(),
      plot.background = element_rect(fill = "#ECF0F4")
    )

  if (scale == "percent") {
    p <- p + scale_y_continuous(labels = scales::percent)
  }
  if (scale == "thousands") {
    p <- p + scale_y_continuous(labels = scales::label_number(scale = 1e-3))
  }

  p
}

# Helper to set format and structure of stats summaries
summarise_distribution <- function(data, yvar) {
  data |>
    group_by(collection_type) |>
    summarise(
      n      = n(),
      mean   = mean({{ yvar }}, na.rm = TRUE),
      median = median({{ yvar }}, na.rm = TRUE),
      sd     = sd({{ yvar }}, na.rm = TRUE),
      skew   = skewness({{ yvar }})
    )
}

The chunk below contains code to support reporting of results.

Show code 
# Convert a proportion to percentage points
pp <- function(x, digits = 1) round(x * 100, digits)

# Count and share of LAs using a given collection system
n_sys   <- function(x) sum(recycling$collection_type == x)
pct_sys <- function(x) round(mean(recycling$collection_type == x) * 100)

# Pull a single value out of a broom tidy tibble by term name
get_term <- function(tbl, term_name, col = "estimate") {
  tbl[[col]][tbl$term == term_name]
}

Visualise the distributions

Before we run a significance test, let’s view both the number of local authorities employing each type of collection system (see Figure 1), and the distributions of recycling rates by collection type (see Figure 2).

Breakdown of collection systems

Show code 
ggplot(recycling, aes(collection_type, fill = collection_type)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), 
          vjust = -0.5, family = "Work Sans", colour = dt_navy) +
scale_fill_dt(guide = "none") +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
labs(
  x = NULL,
  y = NULL,
  title = "Count of local authorities by recycling collection system"
) +
theme_dt() +
theme(
  axis.text.y  = element_blank(),
  axis.ticks.y = element_blank(),
  plot.background = element_rect(fill = "#ECF0F4"),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank()
)
Figure 1: Count of local authorities by collection system. Group sizes are markedly uneven, which informs the choice of statistical test.

The bar chart in Figure 1 shows that the group sizes are uneven. 157 (51%) authorities use a Co-Mingled system, 112 (36%) use a Two Stream system and 40 (13%) use a Multi-Stream system.

Distribution of recycling rates by system

Show code 
plot_distribution(
  recycling,
  household_recycle_rate,
  title = "Distribution of recycling rates by collection system",
  scale = "percent"
)
Figure 2: Distribution of household recycling rates by collection system. The box shows the median and interquartile range; the diamond shows the mean.
Show code 
summarise_distribution(recycling, household_recycle_rate) |> kable(digits = 3)
Table 1: Summary statistics for household recycling rate by collection system
collection_type n mean median sd skew
Co-Mingled 157 0.407 0.403 0.099 0.215
Two Stream 112 0.432 0.443 0.095 -0.129
Multi-Stream 40 0.454 0.445 0.086 -0.038

As shown in Table 1, we can see that Multi-Stream collection systems achieved the highest mean and median recycling rate in 2021, and Co-Mingled systems achieved the lowest rates. The skew number tells us that the distributions are not normal. Co-Mingled systems show a positive skew, while Two Stream and Multi-Stream both show a negative skew.

Why we analyse rate rather than volume

The raw volumes of household waste recycled by local authority are heavily right-skewed, as shown in the plots below. They are dominated by local authority size, since a large council collects and recycles more of everything. This is why the analysis works on the rate of recycling, which normalises for size. It is also why we will use log(population) when we come to control for size in the regression.

Show code 
plot_distribution(
  recycling,
  household_collected,
  title = "Distribution of household waste collected by collection system",
  scale = "thousands"
)
Figure 3: Distribution of household waste collected by collection system. Values are in thousands of tonnes. The box shows the median and interquartile range; the diamond shows the mean.
Show code 
summarise_distribution(recycling, household_collected) |> kable(digits = 3)
Table 2: Summary statistics for amount of household waste collected (tonnes) by collection system
collection_type n mean median sd skew
Co-Mingled 157 69016.55 54330 45496.03 1.962
Two Stream 112 81700.96 56376 64916.62 2.296
Multi-Stream 40 68661.48 50222 50777.70 2.203
Show code 
plot_distribution(
  recycling,
  household_recycled,
  title = "Distribution of household waste recycled by collection system",
  scale = "thousands"
)
Figure 4: Distribution of household waste recycled by collection system. Values are in thousands of tonnes. The box shows the median and interquartile range; the diamond shows the mean.
Show code 
summarise_distribution(recycling, household_recycled) |> kable(digits = 3)
Table 3: Summary statistics for amount of waste recycled (in tonnes) by collection system
collection_type n mean median sd skew
Co-Mingled 157 27397.06 22379.0 18916.55 2.340
Two Stream 112 35662.92 24325.5 30791.84 2.300
Multi-Stream 40 30769.22 24156.0 22197.09 1.639

Test for significance

Kruskal-Wallis test

The default statistical test would be a one-way ANOVA, which compares group means. But ANOVA assumes each group is roughly normally distributed with similar variance, and our groups are uneven (157 / 112 / 40) and somewhat skewed, as shown above. Instead, we will use the Kruskal-Wallis test: the non-parametric equivalent. It works on the ranks of the values instead of the raw values, so it makes no assumption about the shape of the distribution, and it handles unequal group sizes comfortably. It tests whether at least one group tends to sit higher than the others.

Show code 
# Run the KW test for recycling rate as explained by collection type
kw <- kruskal.test(household_recycle_rate ~ collection_type, data = recycling)
kw

    Kruskal-Wallis rank sum test

data:  household_recycle_rate by collection_type
Kruskal-Wallis chi-squared = 9.8269, df = 2, p-value = 0.007347

The test returns H = 9.83 on 2 degrees of freedom, p = 0.007. At least one collection system therefore differs from the others to a statistically significant degree.

This tells us that a difference exists somewhere, but not where, and not how large. Both of those questions matter more than the first.

Effect size

However, a p-value only tells us whether there is an effect, not how big it is. With 309 authorities, even a trivial difference can come out significant. Our question asks about a meaningful impact, so we need an effect size, which measures how much of the variation in recycling rates is explained by collection system.

Show code 
eff <- recycling |> kruskal_effsize(household_recycle_rate ~ collection_type)
eff |> kable(digits = 3)
Table 4: Kruskal-Wallis effect size for recycling rates as explained by collection system
.y. n effsize method magnitude
household_recycle_rate 309 0.026 eta2[H] small

Eta-squared is 0.026, meaning collection system accounts for 2.6% of the variation in recycling rates between local authorities. By the conventional thresholds (0.01 small, 0.06 medium, 0.14 large) this is a small effect. This means that the collection system is not the biggest contributor to the recycling rate.

Dunn’s test

Finally, we can run Dunn’s test to test the collection systems in a pairwise manner, which allows us to see whether any one of them show a significant difference in recycling rate compared to the others.

Show code 
medians <- recycling |> 
  group_by(collection_type) |> 
  summarise(median_rate = median(household_recycle_rate)) |> 
  mutate(collection_type = as.character(collection_type))

dunn <- recycling |> 
  dunn_test(household_recycle_rate ~ collection_type,
            p.adjust.method = "holm") |> 
  left_join(medians, by = join_by(group1 == collection_type)) |> 
  rename(median1 = median_rate) |> 
  left_join(medians, by = join_by(group2 == collection_type)) |> 
  rename(median2 = median_rate) |> 
  mutate(median_gap_pp = pp(median2 - median1)) |> 
  select(-.y., -median1, -median2)

# Function to pull out data from Dunn test
dunn_gap  <- function(g1, g2) dunn$median_gap_pp[dunn$group1 == g1 & dunn$group2 == g2]
dunn_padj <- function(g1, g2) dunn$p.adj[dunn$group1 == g1 & dunn$group2 == g2]
dunn_p <- function(g1, g2) dunn$p[dunn$group1 == g1 & dunn$group2 == g2]

dunn |> kable(digits = 3)
Table 5: Summary of Dunn’s test results
group1 group2 n1 n2 statistic p p.adj p.adj.signif median_gap_pp
Co-Mingled Two Stream 157 112 2.233 0.026 0.051 ns 4.0
Co-Mingled Multi-Stream 157 40 2.754 0.006 0.018 * 4.2
Two Stream Multi-Stream 112 40 1.149 0.251 0.251 ns 0.2

The Kruskal-Wallis test is an omnibus test: it reports that some difference exists without identifying which systems differ. Dunn’s test makes the three pairwise comparisons, using the same pooled ranks. Because running three tests inflates the chance of a false positive, the p-values are adjusted using the Holm correction, and it is the adjusted values (p.adj in the table above) that should be interpreted.

Only one pair differs significantly after correction. Co-Mingled authorities are associated with lower recycling rates than Multi-Stream authorities (p.adj = 0.018), a median gap of 4.2pp.

The other two comparisons are more nuanced. The Co-Mingled versus Two Stream comparison narrowly misses our stated threshold for significance (p.adj = 0.051), and yet the median gap is 4pp, which is comparable in size to the significant result above. The evidence here is suggestive but not conclusive.

The Two Stream versus Multi-Stream comparison is different. It is not significant (p.adj = 0.251) and the median gap is a negligible 0.2pp. Here both the test and the estimate point the same way, towards there being little or no real difference between the two.

Taken together, the pattern is consistent with the burden of separation mattering mainly at the first step, going from mixing everything into one bin to separation into two bins, with little further gain from separating into more streams. The data support this interpretation but do not establish it. With only 40 Multi-Stream authorities, this analysis has limited power to detect modest differences.

Controlling for authority size

Significance testing has now been taken about as far as it can usefully go. Every result has been interpreted as either significant or not significant, with one landing almost exactly on the threshold. We now have two further questions to answer.

Firstly, how large is each difference, and how certain are we? We can answer this question with regression by reporting effect estimates with confidence intervals, which say more than a p-value alone.

Secondly, is the collection system really doing the work? If Multi-Stream authorities happen to be systematically smaller, then what looks like a system effect may be a size effect in disguise. This is an example of a confounding variable. Regression addresses this by comparing systems among authorities of similar size, holding size constant.

Population is used on a log scale. As the volume distributions in Figure 3 and Figure 4 showed, authority size is heavily right-skewed, with a long tail of very large councils. Logging the variable means the model reads proportional differences in size, such as a doubling of population, rather than treating each additional resident as equivalent whether it is the two-thousandth or the millionth.

Co-Mingled is the reference level, so each coefficient represents the difference in recycling rate relative to a Co-Mingled authority. Because the recycling rate is a proportion, a coefficient of 0.04 corresponds to 4pp.

Show code 
mod_1 <- lm(household_recycle_rate ~ collection_type, data = recycling)
mod_2 <- lm(household_recycle_rate ~ collection_type + log(population), data = recycling)

tidy_1   <- tidy(mod_1, conf.int = TRUE)
tidy_2   <- tidy(mod_2, conf.int = TRUE)
glance_1 <- glance(mod_1)
glance_2 <- glance(mod_2)
Show code 
tidy_2 |> kable(digits = 4)
Table 6: Summary of regression for data controlled for population
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.6733 0.1054 6.3858 0.0000 0.4658 0.8808
collection_typeTwo Stream 0.0272 0.0118 2.3100 0.0216 0.0040 0.0504
collection_typeMulti-Stream 0.0451 0.0168 2.6745 0.0079 0.0119 0.0782
log(population) -0.0223 0.0088 -2.5315 0.0119 -0.0397 -0.0050
Show code 
glance_2 |> kable(digits = 4)
Table 7: Further summary of regression for data controlled for population
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.0498 0.0404 0.0951 5.327 0.0014 3 290.6516 -571.3031 -552.6364 2.7572 305 309
Show code 
tidy_1 |> kable(digits = 4)
Table 8: Summary of regression for data uncontrolled for population
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 0.4071 0.0077 53.1791 0.0000 0.3920 0.4221
collection_typeTwo Stream 0.0252 0.0119 2.1213 0.0347 0.0018 0.0485
collection_typeMulti-Stream 0.0465 0.0170 2.7361 0.0066 0.0131 0.0799
Show code 
glance_1 |> kable(digits = 4)
Table 9: Further summary of regression for data uncontrolled for population
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.0298 0.0235 0.0959 4.7032 0.0097 2 287.4389 -566.8779 -551.9445 2.8152 306 309

Interpreting the regression models

The regression analysis showed two results.

Firstly, the system effect was robust to authority size. The collection type coefficients barely changed when population was factored in: Two Stream moved from 2.5pp to 2.7pp, and Multi-Stream from 4.6pp to 4.5pp. If size had been confounding the comparison, these would have reduced materially. The association between collection system and recycling rate is not an artefact of larger or smaller councils favouring particular systems.

Secondly, larger authorities recycle less. The coefficient on log(population) is -0.022 (p = 0.012), which corresponds to a fall of roughly 1.5pp in the recycling rate for each doubling of population. This is a modest but real association, and it is independent of the collection system in use.

The confidence intervals are more informative than the p-values. Relative to Co-Mingled, Two Stream is associated with an increase of 2.7pp (95% CI: 0.4pp to 5pp) and Multi-Stream with an increase of 4.5pp (95% CI: 1.2pp to 7.8pp). Both intervals exclude zero, so both effects are likely real, but both are wide. The Multi-Stream interval spans everything from a barely noticeable 1pp to a substantial 8pp. A fair interpretation of this is that a real effect exists and its likely size is modest, but this dataset cannot pin it down precisely.

There is one apparent discrepancy that we should explain. In the regression, Two Stream differs significantly from Co-Mingled (p = 0.022), whereas Dunn’s test found the same comparison borderline (p.adj = 0.051). These are not contradictory. Dunn’s raw p-value for that pair was 0.026, and the Holm correction for three comparisons roughly doubled it. The regression applies no such correction. The two results agree on the substance: there is a difference of around 3 to 4pp, and the evidence for it is real but not overwhelming.

Finally, the model explains very little. The R-squared is 0.05. This means that 95% of the variation in recycling rates between English local authorities lies outside collection system and population entirely.

Hierarchical clustering

If collection system alone cannot fully explain the difference in recycling rates, there must be other factors at play. We can use a clustering technique to define strata of similar authorities, within which we can test the effect of the collection system on recycling rates.

We can use hierarchical agglomerative clustering with Ward’s D2 linkage on a set of standardised features, with k-means as a stability cross-check. This model is suitable because our n of 309 is small. The method is deterministic, so requires no seed to reproduce the results. Ward’s minimises within-cluster variance, so it tends to give compact, reasonably balanced groups, which is what this stratification needs. Every peer group will be populated enough to compare systems inside it.

The features we have pre-selected for this cluster analysis are:

  • population: the 2021 population estimate for the local authority
  • population_density: population per km^2
  • rurality: the proportion of the population living in rural areas
  • imd_average_score: the 2019 index of deprivation

Building peer groups

The first step is to get a clean, comparably-scaled feature matrix and measure the correlation between them. Since population_density and rurality both measure the urban-rural axis, they are likely to be strongly negatively correlated; if that is around -0.85 or stronger then we are double-weighting one dimension, and we should elect to drop one of the features.

Show code 
# context features only; rate and system deliberately excluded
cluster_data <- recycling |>
  select(lad21cd, local_authority,
         population, population_density, imd_average_score, rurality)

cluster_data |> head() |> kable(digits = 2)
Table 10: Pre-selected features for hierarchical clustering (first 5 rows only)
lad21cd local_authority population population_density imd_average_score rurality
E06000001 Hartlepool 92571 941.23 35.04 0.03
E06000002 Middlesbrough 143734 2634.73 40.46 0.00
E06000003 Redcar and Cleveland 136616 538.31 29.79 0.32
E06000004 Stockton-on-Tees 197030 939.44 25.79 0.06
E06000005 Darlington 108222 548.02 25.66 0.13
E06000006 Halton 128577 1423.55 32.33 0.01

Check for completeness

Show code 
# distance methods cannot take NAs, so check completeness first
cluster_data |>
  summarise(across(everything(), ~ sum(is.na(.x))))
# A tibble: 1 × 6
  lad21cd local_authority population population_density imd_average_score
    <int>           <int>      <int>              <int>             <int>
1       0               0          0                  0                 0
# ℹ 1 more variable: rurality <int>

We have no missing data, so there is no need to drop any rows.

Check for skewness

Show code 
# distribution shape of each feature
cluster_data |>
  summarise(across(c(population, population_density, imd_average_score, rurality),
                   skewness))
# A tibble: 1 × 4
  population population_density imd_average_score rurality
       <dbl>              <dbl>             <dbl>    <dbl>
1       2.74               2.68             0.608    0.870

Scale the features

As expected, population and population_density show strong positive skew and should therefore be log-transformed. imd_average_score and rurality both show a milder positive skew, imd_average_score is mild enough not to require transformation and rurality is bounded at zero, and so cannot be log-transformed.

Once the log transformations are complete, we can then scale all four features so that they are directly comparable. the scale() function in R converts each raw value into its z-score, which is the distance from the mean measured in standard deviations.

Show code 
cluster_scaled <- cluster_data |>
  mutate(
    population         = log10(population),
    population_density = log10(population_density)
    # imd_average_score and rurality stay on their original scale
  ) |>
  select(lad21cd, population, population_density, imd_average_score, rurality) |>
  column_to_rownames("lad21cd") |>
  scale()

cluster_scaled |> head() |> kable(digits = 2)
Table 11: Scaled features for hierarchical clustering (first 5 rows only)
population population_density imd_average_score rurality
E06000001 -0.82 0.18 1.92 -0.79
E06000002 -0.10 0.92 2.60 -0.91
E06000003 -0.19 -0.22 1.26 0.41
E06000004 0.41 0.18 0.76 -0.69
E06000005 -0.56 -0.21 0.74 -0.37
E06000006 -0.28 0.48 1.58 -0.89

Correlation check

Show code 
cor(cluster_scaled) |> round(2)
                   population population_density imd_average_score rurality
population               1.00               0.41              0.35    -0.41
population_density       0.41               1.00              0.42    -0.89
imd_average_score        0.35               0.42              1.00    -0.41
rurality                -0.41              -0.89             -0.41     1.00

As expected, population_density and rurality show a strong negative correlation. Since it is stronger than our pre-defined threshold of -0.85, we should consider that these features are near-duplicates and omit one of them. We will keep rurality as the more direct measure, and drop population_density from the feature set.

Final feature set

Show code 
cluster_final <- cluster_scaled[, c("population", "imd_average_score", "rurality")]

head(cluster_final) |> kable()
Table 12: Final feature set (first 5 rows only)
population imd_average_score rurality
E06000001 -0.8172692 1.9186685 -0.7917248
E06000002 -0.1032319 2.5993062 -0.9136364
E06000003 -0.1856571 1.2603715 0.4105759
E06000004 0.4086029 0.7580828 -0.6866286
E06000005 -0.5637651 0.7413900 -0.3713399
E06000006 -0.2840770 1.5782869 -0.8926171

How many peer groups?

We now need to decide how many clusters k to divide the data into. We will use a dendrogram and a silhouette plot to aid this decision but we must also keep in mind that with 309 authorities and three unevenly sized collection systems, fine-grained clusters risk leaving peer groups with zero or one authority of some system. Stratification only controls confounding if every stratum still holds a spread of systems to compare. So k is bounded from above by hese factors, and not just by dendrogram or silhouette results.

Dendrogram

Show code 
# distance matrix on the final feature set
cluster_dist <- dist(cluster_final, method = "euclidean")

# Ward's hierarchical clustering
cluster_hc <- hclust(cluster_dist, method = "ward.D2")

# dendrogram
fviz_dend(cluster_hc, show_labels = FALSE)

In the dendrogram, at approximately height 5, k = 10. At approximately height 10, k = 5.

Silhouette

Show code 
# silhouette width across candidate k
fviz_nbclust(cluster_final, FUNcluster = hcut, method = "silhouette")

k = 2 is by far the most optimized number of clusters according to the silhouette, but that is a usual result. Silhouette rewards the cleanest split, but two groups of roughly 150 authorities each is barely a peer group. k = 3, 7, 9 are reasonably well supported by the silhouette.

As mentioned earlier, neither plot is the decision on its own, because the binding constraint that the selected k has to leave every peer group with a workable spread of all three systems, and Multi-Stream is the rare category, which caps how fine we can go.

Based on the constraint and the interpretation of the two charts, we will shortlist k = 5 and k = 3, and run a crosstab to make the final decision.

Crosstab

Show code 
# Assign local authorities into either 3 or 5 clusters using cutree
recycling <- recycling |>
  mutate(
    peer_k3 = factor(cutree(cluster_hc, k = 3)[lad21cd]),
    peer_k5 = factor(cutree(cluster_hc, k = 5)[lad21cd])
  )
Show code 
# does every peer group retain all three systems?
recycling |> count(peer_k5, collection_type) |>
  pivot_wider(names_from = collection_type, values_from = n, values_fill = 0) |> 
  mutate(n = rowSums(across(where(is.numeric))), .before = `Co-Mingled`) |> 
  kable()
Table 13: Peer group statistics for k = 5
peer_k5 n Co-Mingled Two Stream Multi-Stream
1 58 29 22 7
2 65 33 19 13
3 60 31 24 5
4 42 20 17 5
5 84 44 30 10
Show code 
recycling |> count(peer_k3, collection_type) |>
  pivot_wider(names_from = collection_type, values_from = n, values_fill = 0) |> 
  mutate(n = rowSums(across(where(is.numeric))), .before = `Co-Mingled`) |> 
  kable()
Table 14: Peer group statistics for k = 3
peer_k3 n Co-Mingled Two Stream Multi-Stream
1 58 29 22 7
2 149 77 49 23
3 102 51 41 10

The crosstab shows that both k = 3 and k = 5 produce clusters that contain all three collection systems.

We will select k = 5 for further analysis, for the following reasons:

  • It produces finer stratification on the features we selected.
  • It produces more balanced totals.

Peer signature

Let’s review summary statistics for k = 5.

Show code 
peer_signature <- cluster_final |> 
  as.tibble(rownames = "lad21cd") |> 
  left_join(
    select(
      recycling,
      lad21cd,
      peer_k5
    ),
    by = "lad21cd"
  ) |> 
  group_by(peer_k5) |> 
  summarise(
    across(c(population, imd_average_score, rurality), mean),
    n = n(),
    .groups = "drop"
  )

peer_signature |> kable(digits = 2)
Table 15: Peer signature for k = 5
peer_k5 population imd_average_score rurality n
1 -0.47 0.59 -0.74 58
2 -0.86 -0.31 1.42 65
3 1.12 1.14 -0.81 60
4 0.96 -0.23 -0.31 42
5 -0.30 -0.87 0.14 84

In order to label the clusters, we look at the mean z-score of each feature.

Show code 
peer_labels <- tribble(
  ~peer_k5, ~peer_label,
  "1", "Urban, deprived",
  "2", "Rural, small",
  "3", "Urban, large, high deprivation",
  "4", "Urban, large, low deprivation",
  "5", "Semi-Rural, low deprivation"
)

# Add the labels to recycling data
recycling <- recycling |>
  mutate(peer_k5 = as.character(peer_k5)) |>
  left_join(peer_labels, by = "peer_k5")

recycling |> 
  group_by(peer_label) |> 
  summarise(median(population), mean(imd_average_score), mean(rurality), n()) |> 
  kable(digits = 2)
Table 16: Labelled clusters
peer_label median(population) mean(imd_average_score) mean(rurality) n()
Rural, small 99435.0 17.30 0.56 65
Semi-Rural, low deprivation 130034.0 12.82 0.25 84
Urban, deprived 112352.0 24.47 0.04 58
Urban, large, high deprivation 288597.0 28.82 0.03 60
Urban, large, low deprivation 275111.5 17.94 0.15 42

Comparing systems within peer groups

We can now build a regression model to investigate the effect of collection system on recycling rate within the five peer groups.

The model

Peer group goes in as a fixed factor, not a random effect. A mixed model wants many more than five groups to estimate a variance reliably, so with k = 5 we choose fixed blocking: four dummy coefficients that we treat as nuisance adjustment.

Show code 
# peer group as a fixed blocking factor
mod_3 <- lm(household_recycle_rate ~ collection_type + peer_label, data = recycling)

tidy_3   <- tidy(mod_3, conf.int = TRUE)
glance_3 <- glance(mod_3)

Testing the system effect

Show code 
mod_peer <- lm(household_recycle_rate ~ peer_label, data = recycling)
anova(mod_peer, mod_3)   # what does adding collection_type buy over peer group alone?
Analysis of Variance Table

Model 1: household_recycle_rate ~ peer_label
Model 2: household_recycle_rate ~ collection_type + peer_label
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    304 2.1677                                
2    302 2.0863  2   0.08139 5.8908 0.003093 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusion

Does the collection system have a meaningful impact on recycling rates?

The answer has three parts, and they pull in different directions.

Yes, there is a real difference, and it is not trivial in size. Authorities that require householders to separate their recycling are associated with higher recycling rates than those that collect everything co-mingled. Adjusting for authority size, the gap is 2.7pp for Two Stream and 4.5pp for Multi-Stream, relative to Co-Mingled. A 4pp shift in the national recycling rate would be a substantial policy outcome, so this is not a difference that can be dismissed as academic.

But the effect appears to be a step, not a ladder. The evidence does not support the intuition that asking householders to sort into ever more streams yields ever higher rates. Two Stream and Multi-Stream authorities are statistically indistinguishable from one another, with a median gap of 0.2pp. What appears to matter is whether households separate their recycling at all, not how finely they separate it. If this interpretation is correct, it has a clear practical implication: the additional burden that Multi-Stream places on householders, relative to Two Stream, gives no detectable improvement. This is the most interesting finding in the analysis, and it is also the one held most tentatively, resting as it does on only 40 Multi-Stream authorities.

And most of what drives recycling rates lies elsewhere. Collection system alone explains 3% of the variation between authorities. Adding population raises this only to 5%. The other 95% belongs to factors this analysis cannot see. A local authority hoping to transform its recycling performance by changing its bin system alone would, on this evidence, be disappointed. The system is one lever among many, and not the largest.

There are also sound reasons why an authority might choose a system irrespective of its effect on rates. Without access to a Materials Recovery Facility, where co-mingled waste can be mechanically sorted, a council may have little practical alternative to kerbside sorting. The choice of system is constrained by infrastructure, cost and geography, not made freely on the basis of expected recycling rate.

Data sources and limitations

  1. Observational, not experimental. Local authorities chose their own collection systems. Any of the associations reported here could reflect the kind of authority that adopts a given system rather than the system itself. Nothing in this analysis establishes cause.

  2. Residual confounding is likely. The only control available in this dataset is population. The factors that the recycling literature would flag most strongly are absent: deprivation, rural versus urban character, and housing stock. Housing type in particular is a known driver, since flats recycle markedly less well than houses whatever bins are provided, and it is plausibly correlated with the choice of system. The regression should therefore be read as adjusting for authority size only.

  3. A single year. All figures relate to 2021. Rates fluctuate, and 2021 was likely affected by altered patterns of household waste generation during the pandemic. A multi-year panel would be more robust and would allow authorities that changed system to be observed before and after, which is a considerably stronger design.

  4. Unequal and small groups. Only 40 of the 309 authorities operate Multi-Stream systems. Comparisons involving this group have limited statistical power, and the finding that Two Stream and Multi-Stream cannot be distinguished may reflect that limitation rather than a true absence of difference.

  5. System classification. Collection system data was obtained by Freedom of Information request and reflects the system in place at the time of response. Authorities operating hybrid or transitional arrangements are assigned to a single category, which introduces some measurement error.

Further work

  • Test Two Stream against Multi-Stream directly by releveling the factor, rather than inferring it from two comparisons against a shared reference.
  • Join contextual variables (deprivation, rural or urban classification, housing type) and re-fit the model. If the system effect survives these controls, the finding is considerably stronger.
  • Use those contextual variables to construct peer groups of comparable authorities through cluster analysis, deliberately excluding both the recycling rate and the collection system, and then compare systems within peer groups. This controls for confounding by stratification.
  • Extend to multiple years to exploit authorities that changed system.