Who Earns What in England? A Visual Investigation of the 2021 Household Census Snapshot

library(readxl); library(dplyr); library(tidyr); library(ggplot2)
library(forcats); library(scales); library(patchwork)
theme_set(theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"), legend.position = "top",
        panel.grid.minor = element_blank(), plot.title.position = "plot"))
sex_cols <- c(Female = "#C2185B", Male = "#1565C0")
ed_lv <- c("Less than Secondary School", "Secondary School", "Some HE",
           "Bachelors Degree", "Masters or higher")
raw <- read_xlsx("census.xlsx", na = c("", "NA"))  # missing income is stored as text "NA"
census <- raw |>
  rename(Highest_Ed = `Highest Ed`) |>
  mutate(Highest_Ed = factor(trimws(Highest_Ed), levels = ed_lv),
         Mar_Stat = fct_recode(Mar_Stat, Married = "married"),
         Sex = factor(if_else(Female == 1, "Female", "Male")), # see Section 2
         Eth = fct_infreq(Eth),
         Age_band = cut(Age, c(-1, 14, 24, 34, 44, 54, 64, 74, 120),
           labels = c("0-14","15-24","25-34","35-44","45-54","55-64","65-74","75+")))
earners <- filter(census, INC > 0)  # adults reporting a positive income

1. Introduction

The data are a modified snapshot of the 2021 household census for England: 27,410 people living in 10,565 households, described by age, sex, marital status, ethnicity, education, income and one housing question (H8: whether the household has private, lockable use of all its rooms). Rather than describing each variable in turn, this report follows one question — what shapes annual income? — and, along the way, checks whether the data can be trusted. Two of the most interesting findings turn out to be about the data itself, not the population.

2. Getting to know the data: a coding problem hiding in plain sight

The codebook says Female = 0 means yes, female. Before relying on this, a basic demographic sanity check is worth doing: women in England live longer, so they should dominate the oldest ages and the widowed.

census |> count(Age_band, Sex) |>
  mutate(n = if_else(Sex == "Male", -n, n)) |>
  ggplot(aes(n, Age_band, fill = Sex)) + geom_col(width = 0.85) +
  scale_fill_manual(values = sex_cols) +
  scale_x_continuous(labels = \(x) comma(abs(x))) +
  labs(title = "Older age bands lean heavily towards the group coded Female = 1",
       x = "Number of people", y = "Age band", fill = NULL)

Figure 1: Figure 1. Age-sex pyramid, with sex taken from Female = 1 (see text).

census |> filter(!is.na(Mar_Stat)) |> count(Mar_Stat, Sex) |>
  group_by(Mar_Stat) |> mutate(p = n / sum(n)) |>
  ggplot(aes(p, fct_reorder2(Mar_Stat, Sex, p), fill = Sex)) +
  geom_col(width = 0.7) + geom_vline(xintercept = 0.5, linetype = 2) +
  scale_fill_manual(values = sex_cols) + scale_x_continuous(labels = percent) +
  labs(title = "Four in five widowed people are coded Female = 1",
       x = "Share of people in each status", y = NULL, fill = NULL)

Figure 2: Figure 2. Sex composition of each marital status (adults only).

Commentary.

I therefore treat Female = 1 as female throughout. The point matters: under the documented coding, every later conclusion about the gender pay gap would be reversed.

3. Missing values are structural, not random

census |> group_by(Age_band) |>
  summarise(across(c(INC, Mar_Stat, Highest_Ed), \(x) mean(is.na(x)))) |>
  pivot_longer(-Age_band, names_to = "Variable", values_to = "Missing") |>
  ggplot(aes(Age_band, Variable, fill = Missing)) +
  geom_tile(colour = "white", linewidth = 1) +
  geom_text(aes(label = percent(Missing, 1)), size = 3.4) +
  scale_fill_gradient(low = "#F1F8E9", high = "#E65100", labels = percent) +
  labs(title = "Missingness is almost entirely a children's phenomenon",
       x = "Age band", y = NULL, fill = "Missing")

Figure 3: Figure 3. Share of missing values by age band and variable.

Commentary. Income and marital status are blank for essentially all children under 15 and complete for adults; education is missing only for the under-3s. This is by design (the census does not ask these questions of children), so deleting or imputing these rows would be wrong. Income analysis below is restricted to people with a reported income above zero (n = 18,658); 2,579 adults report exactly £0 and are discussed separately where relevant.

4. Education is the strongest visual predictor of income

med_ed <- earners |> filter(!is.na(Highest_Ed)) |>
  group_by(Highest_Ed) |> summarise(m = median(INC))
earners |> filter(!is.na(Highest_Ed)) |>
  ggplot(aes(INC, Highest_Ed)) +
  geom_violin(fill = "#90CAF9", colour = NA, alpha = 0.6) +
  geom_boxplot(width = 0.18, outlier.shape = NA, fill = "white") +
  geom_label(data = med_ed, aes(m, label = label_currency(prefix = "£")(m)),
             nudge_y = 0.38, size = 3.2, label.size = 0) +
  scale_x_log10(labels = label_currency(prefix = "£", scale_cut = cut_short_scale())) +
  labs(title = "Each step up the education ladder roughly raises the median income",
       x = "Annual income (log scale)", y = NULL)

Figure 4: Figure 4. Income by highest qualification (log scale); labels show medians.

Commentary. Median income rises at every level, from about £10.5k for those without secondary education to £50k with a Masters or higher, a roughly five-fold difference. A log scale is essential here: income is heavily right-skewed (maximum £720,000), and on a linear axis every box collapses to the left edge. The wide lower tails of the violins show that a qualification raises the typical income but does not guarantee it.

5. Income across the life course

earners |> filter(!is.na(Highest_Ed), Age >= 18) |>
  group_by(Age_band, Highest_Ed) |> summarise(m = median(INC), n = n(), .groups = "drop") |>
  filter(n >= 30) |>
  ggplot(aes(Age_band, m, colour = Highest_Ed, group = Highest_Ed)) +
  geom_line(linewidth = 1.1) + geom_point(size = 2) +
  scale_colour_viridis_d(end = 0.9, direction = -1) +
  scale_y_continuous(labels = label_currency(prefix = "£", scale_cut = cut_short_scale())) +
  guides(colour = guide_legend(nrow = 2)) +
  labs(title = "Graduates' earnings climb steeply to a peak in mid-career",
       x = "Age band", y = "Median income", colour = NULL)

Figure 5: Figure 5. Median income by age band and education.

Commentary. All groups follow the familiar hump-shaped life-cycle profile, rising to a peak between 35 and 64 and falling with retirement. The gap between groups is not constant, though: lines fan out between 25 and 54 as higher qualifications compound with experience, then narrow after 65 when pension income replaces earnings. Education matters most during working life.

6. The gender gap persists at every level of education

gap <- earners |> filter(!is.na(Highest_Ed), between(Age, 25, 64)) |>
  group_by(Highest_Ed, Sex) |> summarise(m = median(INC), .groups = "drop")
gap_lab <- gap |> pivot_wider(names_from = Sex, values_from = m) |>
  mutate(pct = percent(1 - Female / Male, 1))
ggplot(gap, aes(m, Highest_Ed)) +
  geom_line(aes(group = Highest_Ed), colour = "grey70", linewidth = 2) +
  geom_point(aes(colour = Sex), size = 4) +
  geom_text(data = gap_lab, aes(x = Male, label = pct), nudge_x = 5500, size = 3.4) +
  scale_colour_manual(values = sex_cols) +
  scale_x_continuous(labels = label_currency(prefix = "£", scale_cut = cut_short_scale()),
                     expand = expansion(mult = c(0.05, 0.12))) +
  labs(title = "At every qualification level, women earn less than men",
       subtitle = "Numbers show the percentage by which women's median falls short",
       x = "Median annual income", y = NULL, colour = NULL)

Figure 6: Figure 6. Median income of working-age earners (25-64) by sex and education.

Commentary. Comparing like with like (same qualification, working age) removes the excuse that the gap simply reflects education: women’s median income is lower at every level. Notably, Figure 2’s coding check is what makes this reading possible; under the codebook’s literal coding the chart would wrongly show women out-earning men by a similar margin.

7. A housing question that only Asian households answered “no”

p1 <- census |> group_by(Eth) |> summarise(p = mean(H8 == 1)) |>
  ggplot(aes(p, Eth)) + geom_col(fill = "#6A1B9A", width = 0.6) +
  geom_text(aes(label = percent(p, 1)), hjust = -0.15, size = 3.4) +
  scale_x_continuous(labels = percent, limits = c(0, 1.12)) +
  labs(title = "H8 = 'no' by ethnicity", x = "Share of people", y = NULL)
p2 <- earners |> ggplot(aes(INC, Eth)) +
  geom_boxplot(fill = "#CE93D8", outlier.alpha = 0.08) +
  scale_x_log10(labels = label_currency(prefix = "£", scale_cut = cut_short_scale())) +
  labs(title = "Income by ethnicity", x = "Income (log)", y = NULL)
p1 + p2

Figure 7: Figure 7. Share of people answering H8 = ‘no’ (shared rooms), and income, by ethnicity.

Commentary. This is the most surprising result in the dataset. Every Asian respondent answered “no” to H8 and no one else did: the housing variable is a perfect copy of the Asian ethnicity category. In real census data, shared accommodation is spread across all groups, so this is almost certainly an artefact of how the snapshot was modified. The practical consequence is that H8 and ethnicity cannot be separated: a naive chart of “income by H8” (earner medians of £15k vs £23k) would look like a housing effect but is really an ethnicity comparison, and vice versa. Any model including both would be perfectly collinear. Ethnic income differences themselves are visible in the right panel, with White respondents having the highest median income and Asian and Other the lowest, but these should be read alongside the H8 caveat.

8. Conclusions

• Check the codebook against the data. Demographic patterns (Figures 1-2) show the Female coding is inverted relative to its documentation; using it as written would reverse every finding about sex.
• Missingness is structural. Blanks belong to children, not careless recording (Figure 3), so they should be filtered, not imputed.
• Education is the clearest driver of income, with a five-fold difference in medians (Figure 4) that widens through mid-career (Figure 5).
• A gender gap remains within every education level (Figure 6).
• H8 is perfectly confounded with Asian ethnicity (Figure 7), an artefact of the modified snapshot that makes the two variables inseparable.

The wider lesson is that a visual investigation is not only a way to present results; plotting simple demographic relationships first exposed two data problems that would otherwise have silently corrupted the conclusions.