Unemployment Rates by Occupation: Raw vs. Shrunken

A hierarchical empirical Bayes adjustment of LFS occupational unemployment rates, 2024–2025

Author

Richard Martin

Published

July 21, 2026

Code
library(tidyverse)
library(here)
library(janitor)
library(readxl)
library(vroom)
library(lme4)
library(conflicted)
conflicts_prefer(dplyr::filter)
#functions-----------------------
fit_at_rm <- function(rm) {
  # Sample membership is fixed: lfs was already filtered on the labour force
  # count, which does not depend on rm. Only the weights change.
  d <- lfs |> mutate(n = lf / sum(lf) * rm)

  m <- glmer(urate ~ 1 + (1 | noc_2/noc_5),
             weights = n,
             family = binomial, data = d,
             control = glmerControl(optimizer = "bobyqa",
                                    optCtrl = list(maxfun = 1e5)))

  d |>
    mutate(rm        = rm,
           urate_raw = urate,
           urate_eb  = fitted(m))
}
#constants------------------
max_year  <- year(today()) - 1   # last complete year
years_use <- (max_year - 1):max_year
cut_n <- 125 #HOO size
# BC respondent-months IN THE LABOUR FORCE, 2024-25.
#
# COUNTED, not reconstructed. BC labour-force records in the LFS public use
# microdata files:
#
#   Jan 2024   8,249   (2015 sample design)
#   Jan 2025   8,754   (2015 sample design)
#   Dec 2025   8,058   (2025 sample design)
#   mean      ~8,350   x 24 months = ~200,000 respondent-months
#
# The three months bracket the April-September 2025 design change and show no
# systematic shift, so the spread is month-to-month noise rather than design.
# That ~9% spread is the practical precision here; nothing finer is meaningful.
resp_months <- 200000
#read the data--------------------
noc21_descriptions <- readxl::read_excel(
  here("data", "2021_noc_descriptions.xlsx"),
  col_types = "text") |>
  mutate(noc_5 = str_pad(noc_5, width = 5, side = "left", pad = "0"))

status_by_noc <- vroom(
  list.files("data", pattern = "stat", full.names = TRUE),
  col_types = vroom::cols(
    SYEAR    = vroom::col_double(),
    NOC_5    = vroom::col_character(),  # character so leading zeros not stripped
    LF_STAT  = vroom::col_character(),
    `_COUNT_` = vroom::col_double())) |>
  clean_names() |>
  filter(noc_5 != "missi") |>
  mutate(noc_5 = str_pad(noc_5, width = 5, pad = "0")) |>
  full_join(noc21_descriptions) |>
  filter(
    syear %in% years_use,
    lf_stat %in% c("Employed", "Unemployed")) |>
  summarize(count = sum(count), .by = c(lf_stat, noc_5, class_title))

allocate_to <- status_by_noc |>
  filter(noc_5 %in% c(41220, 41221)) |>
  mutate(prop = count / sum(count), .by = lf_stat)

allocate_from <- status_by_noc |>
  filter(noc_5 == 41229) |>
  select(lf_stat, to_be_allocated = count)

allocated <- full_join(allocate_to, allocate_from) |>
  mutate(count = count + prop * to_be_allocated) |>
  select(-prop, -to_be_allocated)

noc_5 <- status_by_noc |>
  filter(!noc_5 %in% c(41220, 41221, 41229)) |>
  bind_rows(allocated) |>
  arrange(noc_5) |>
  pivot_wider(names_from = "lf_stat", values_from = "count") |>
  mutate(noc_2 = str_sub(noc_5, 1, 2), .after = "noc_5") |>
  mutate(urate = Unemployed / (Employed + Unemployed))

lfs <- noc_5 |>
  filter(Employed + Unemployed > 0) |>
  mutate(
    lf    = Employed + Unemployed,
    # Each occupation's share of the survey's respondents is taken to equal its
    # share of employment. Not rounded: the model does not need integers.
    n     = lf / sum(lf) * resp_months,
    urate = Unemployed / lf,
    y     = urate * n
  )

rm_grid <- c(150000, 175000, resp_months, 225000, 250000)

# --- rounding analysis -------------------------------------------------------
# RTRA rounds every published count to the nearest 250. A cell therefore reports
# zero unemployed whenever the true count is anywhere in 0-124. For an occupation
# with labour force lf, that bounds the true rate at 124/lf: a large lf makes a
# recorded zero an informative statement (the rate is provably tiny), not an
# uninformative one.
round_base <- 250
zero_cells <- lfs |>
  filter(Unemployed == 0) |>
  mutate(rate_ceiling = (round_base/2 - 1) / lf)

zero_bands <- zero_cells |>
  mutate(band = cut(rate_ceiling,
    breaks = c(0, 0.005, 0.02, Inf),
    labels = c("Provably below 0.5%", "Provably below 2%",
               "Genuinely uncertain"))) |>
  count(band, name = "Occupations")

n_zero        <- nrow(zero_cells)
n_provably_lo <- sum(zero_cells$rate_ceiling < 0.005)
biggest_zero  <- zero_cells |> slice_max(lf, n = 1)


sweep <- map_dfr(rm_grid, fit_at_rm)

ranked <- sweep |>
  mutate(rk = rank(urate_eb), .by = rm) |>
  summarize(
    times_in = sum(rk <= cut_n),
    best_rk  = min(rk),
    worst_rk = max(rk),
    swing    = max(rk) - min(rk),
    .by = c(noc_5, class_title, urate_raw)
  )

# The two occupations used to illustrate the mechanism below. Picked from the
# data rather than named, so they survive a change to the grid or the counts.
#   - a zero occupation that moves a lot: little own evidence, so borrowing
#     decides where it sits, and that depends on the respondent count
#   - a well-measured occupation with a real rate: its own estimate barely
#     moves, but its rank does, because the zeros climb past it
illus_zero <- ranked |>
  filter(urate_raw == 0) |>
  slice_max(swing, n = 1)

illus_solid <- ranked |>
  filter(urate_raw > 0) |>
  left_join(
    sweep |> filter(rm == resp_months) |> select(noc_5, n),
    by = "noc_5"
  ) |>
  filter(swing >= quantile(swing, 0.9)) |>
  slice_max(n, n = 1)

illus_pair <- c(illus_zero$class_title, illus_solid$class_title)

illus_tbl <- sweep |>
  mutate(rk = rank(urate_eb), .by = rm) |>
  filter(class_title %in% illus_pair)

# n varies with rm, so report it at the adopted count rather than arbitrarily.
illus_range <- function(occ) {
  d <- illus_tbl |> filter(class_title == occ)
  tibble(
    lo = min(d$urate_eb),
    hi = max(d$urate_eb),
    n  = d$n[d$rm == resp_months][1]
  )
}
illus_zero_range  <- illus_range(illus_zero$class_title)
illus_solid_range <- illus_range(illus_solid$class_title)

cutoff_stability <- ranked |>
  filter(times_in > 0) |>
  summarize(
    `Ever selected`          = dplyr::n(),
    `Selected at every value` = sum(times_in == length(rm_grid)),
    `Movers`                 = sum(times_in > 0 & times_in < length(rm_grid))
  )

fit <- glmer(
  urate ~ 1 + (1 | noc_2 / noc_5),
  weights = n,
  family = binomial,
  data = lfs,
  control = glmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 1e5)
  )
)

vc <- as.data.frame(VarCorr(fit))

sd_noc2 <- vc |> filter(grp == "noc_2") |> pull(sdcor)
sd_noc5 <- vc |> filter(grp == "noc_5:noc_2") |> pull(sdcor)
v_noc2  <- vc |> filter(grp == "noc_2") |> pull(vcov)
v_noc5  <- vc |> filter(grp == "noc_5:noc_2") |> pull(vcov)
centre  <- fixef(fit)[["(Intercept)"]]

grand_logodds <- fixef(fit)[["(Intercept)"]]
grand_rate    <- plogis(grand_logodds)

re_noc2 <- ranef(fit)$noc_2 |>
  rownames_to_column("noc_2") |>
  rename(offset_noc2 = `(Intercept)`)

size_noc2 <- lfs |>
  summarize(
    n_noc2        = sum(n),
    y_noc2        = sum(y),
    n_occupations = dplyr::n(),
    .by = noc_2
  ) |>
  mutate(urate_pooled_noc2 = y_noc2 / n_noc2)

middle_layer <- size_noc2 |>
  left_join(re_noc2, by = "noc_2") |>
  mutate(
    logodds_noc2  = grand_logodds + offset_noc2,
    urate_eb_noc2 = plogis(logodds_noc2)
  ) |>
  arrange(desc(n_noc2))

re_noc5 <- ranef(fit)$`noc_5:noc_2` |>
  rownames_to_column("key") |>
  rename(offset_noc5 = `(Intercept)`) |>
  separate(key, into = c("noc_5", "noc_2"), sep = ":")

eb <- lfs |>
  mutate(
    urate_raw = urate,
    urate_eb  = fitted(fit)
  ) |>
  left_join(re_noc5, by = c("noc_5", "noc_2")) |>
  left_join(
    middle_layer |> select(noc_2, urate_eb_noc2, n_noc2),
    by = "noc_2"
  )

# --- the worked example: NOC major group 32 ---
focus_group <- "32"

# A single position for every zero, set below the smallest non-zero rate in the
# group so the zeros form a clean floor rather than overlapping real estimates.
group_floor <- eb |>
  filter(noc_2 == focus_group, urate_raw > 0) |>
  pull(urate_raw) |>
  min() / 2

# Only label breaks that fall inside the plotted range.
occ_breaks <- c(0.002, 0.005, 0.01, 0.02, 0.05, 0.10, 0.20)
occ_breaks <- occ_breaks[occ_breaks > group_floor]

one_group <- eb |>
  filter(noc_2 == focus_group) |>
  mutate(
    # A zero-unemployment occupation has no finite own-data log-odds. Place ALL
    # zeros at a single common floor, below the smallest non-zero rate in the
    # group. Using 0.5/n instead would scatter them by sample size, which is
    # exactly the information the raw rate does not contain: every one of these
    # occupations reported the same raw number, zero. The differences belong on
    # the fitted side, not the raw side.
    zero_raw = Unemployed == 0,
    lo_raw   = qlogis(if_else(zero_raw, group_floor, urate_raw)),
    lo_eb    = qlogis(urate_eb),
    label    = str_trunc(class_title, 38)
  )

group_target <- middle_layer |>
  filter(noc_2 == focus_group) |>
  pull(logodds_noc2)

zeros    <- one_group |> filter(Unemployed == 0)
tiny_zero <- zeros |> slice_min(n, n = 1)
big_zero  <- zeros |> slice_max(n, n = 1)

The problem in one sentence

The current High Opportunity Occupations (HOO) methodology puts 10% weight on the unemployment rate, but a quarter of occupations report a raw rate of exactly 0% — artifacts of small samples and of rounding, not genuinely unemployment-free occupations — so that weight lands largely on noise.

The problem in one paragraph

We want a true-ish current unemployment rate for every one of roughly 500 detailed occupations (five-digit NOC codes), because the current unemployment rate is one of seven indicators feeding a measure of occupational opportunity. The Labour Force Survey is our only source of (timely) labour force data, and it was never designed to carry that much detail. Some occupations are represented by thousands of survey responses; others by very few. The rate we compute directly from the survey — the raw rate — is trustworthy for the first group and close to noise for the second. This document describes a method for producing a rate for every occupation that is honest about how much evidence stands behind each one.

Why the raw rate misleads for small occupations

Code
thirty_zero <-  dbinom(
     x = 0,      # observed successes
     size = 30,  # sample size
     prob = 1/30 # true success probability
)

thirty_one <-  dbinom(
     x = 1,      # observed successes
     size = 30,  # sample size
     prob = 1/30 # true success probability
)

sixty_zero <-  dbinom(
     x = 0,      # observed successes
     size = 60,  # sample size
     prob = 1/30 # true success probability
)

one_twenty_zero <-  dbinom(
     x = 0,      # observed successes
     size = 120,  # sample size
     prob = 1/30 # true success probability
 )

Suppose an occupation’s true unemployment rate is 3.33%, or one worker in thirty. If your sample size happens to be 30, you are almost equally likely to observe zero unemployment (36.17%) as the true rate of 3.33% (37.41%). If you double the sample size, the probability of observing zero unemployment (if the truth is still 3.33%) drops to 13.08%. Doubling the sample size again, the probability of observing zero unemployment becomes negligible at 1.71%. If we use two years of monthly LFS data, 223 occupations do not reach the above “negligible” threshold of 120 observations.

The problem with ranking occupations by their raw unemployment rate is that small sample sizes mechanically produce many extreme values, including numerous zero rates. As a result, the occupations at the extremes of the ranking are often the smallest — not because their true unemployment rates are exceptional, but because sampling variation is greatest when samples are small. A ranking based on raw unemployment rates therefore becomes, in part, a ranking of sample sizes rather than of the underlying measure correlated with opportunity: the true unemployment rate.

Sampling is not the only thing that manufactures zeros, though, and for the low end it is not even the main one. Statistics Canada rounds its published counts, which drives many large occupations to a recorded zero as well — a second mechanism, with the opposite implication, that the section on the zero problem takes up below.

Shrinkage you can see: one major group

The clearest view of the model is a single major group’s occupations, each shown twice: where its own data puts it, and where the model puts it. Every occupation lies between its own rate and its group’s rate, and the only question is how far along.

Code
ggplot(one_group, aes(x = reorder(label, lo_eb))) +
  geom_hline(yintercept = group_target, colour = "grey40",
             linetype = "dashed", linewidth = 0.8) +
  geom_segment(
    aes(xend = reorder(label, lo_eb), y = lo_raw, yend = lo_eb),
    arrow = arrow(length = unit(0.1, "cm"), type = "closed"),
    colour = "grey55"
  ) +
  geom_point(aes(y = lo_raw, size = n, shape = zero_raw), alpha = 0.4) +
  geom_point(aes(y = lo_eb), colour = "firebrick", size = 1.8) +
  scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 1),
                     labels = c(`FALSE` = "some unemployed", `TRUE` = "none unemployed"),
                     name = NULL) +
  scale_y_continuous(breaks = qlogis(occ_breaks),
                     labels = scales::percent(occ_breaks, accuracy = 0.1)) +
  coord_flip() +
  labs(
    x = NULL,
    y = "Unemployment rate",
    size = "Respondent-months",
    title = paste0("Shrinkage within major group ", focus_group),
    subtitle = "Grey = own data only; red = model estimate.\nDashed line = the group's fitted rate, which is what they are pulled toward."
  ) +
  theme_minimal()

Major group 32 (technical health occupations). Grey = the occupation’s own data; red = the model’s estimate. The dashed line is the group’s fitted rate. The ten occupations with no unemployed respondents all sit at the same floor — they reported the same raw number — but fan out to different fitted rates according to how much each zero was worth.

Group 32 is a good illustration because ten of its eighteen occupations recorded no unemployed respondents at all. Taken at face value, each of those has an unemployment rate of exactly zero, and all ten would tie for the lowest rate in the province. They appear as hollow points, all at the same position on the floor of the plot — because they all reported the same raw number.

The model does not treat them alike, and the reason is the whole point of the exercise. Denturists has only about 12 respondent-months behind its zero — barely any evidence — and is pulled almost all the way to the group’s rate. Massage therapists has about 451, and its zero survives partly intact: it stays visibly below the group, because observing no unemployment across that many respondent-months is real evidence of low unemployment.

Both had a raw rate of zero. They end up in different places, and the distance each one travelled is the model’s statement about how much its zero was worth. That is exactly what an opportunity measure needs: the zeros that mean something stay low, and the zeros that mean nothing do not get to claim first place on the strength of being small.

NoteWhy the axis is spaced oddly

The tick labels are unemployment rates, but they are not evenly spaced — the gap from 1% to 2% takes as much room as the gap from 5% to 10%. That is deliberate.

The model works in log-odds, and shrinkage is a straight pull toward the dashed line on that scale. If we plotted rates on an ordinary axis, the pull would look crooked: some arrows would appear to overshoot the line and a few to move away from it, purely because converting log-odds to rates is a curve rather than a straight line. Nothing would be wrong with the model — the picture would just be lying about the geometry.

The zero problem

For this work the failure has a specific and awkward shape, because we are looking for occupations with low unemployment — signals of opportunity — rather than high.

The lowest possible raw rate is zero, and zero is common. Not because those occupations have solved unemployment, but for two distinct reasons that a raw ranking cannot tell apart. The first is sampling: if you observe thirty people in an occupation and none happens to be unemployed, the arithmetic returns 0%. Over the 2024:2025 window 130 occupations record no unemployed respondents at all. Taken at face value, every one of those zeros ties for first place. A ranking of raw rates would hand the top of the opportunity list to occupations we often know very little about. That is precisely backwards, and it is the failure this method exists to prevent.

The second reason is rounding, and it turns out to matter more than the first — enough to warrant its own treatment below.

Most of the zeros are not what they seem

Statistics Canada rounds every count it releases to the nearest 250. An occupation is therefore recorded with zero unemployed whenever its true count of unemployed people is anywhere from 0 to 124 — the rounding carries it down to zero. This is a different mechanism from the sampling story above, and it has the opposite implication.

The reason it matters is that the same rounding puts a ceiling on the true rate. If an occupation has a labour force of 177,250 and reports zero unemployed, its true unemployment count is at most 124, so its true rate is at most 0.07%. That is not an uninformative zero. It is a near-certain statement that the occupation’s unemployment rate is very low — exactly the signal an opportunity measure is trying to find, disguised as a tie for first place with occupations we truly know nothing about.

Sorting the 130 zero occupations by how tightly the rounding pins their rate:

Code
zero_bands |>
  knitr::kable(
    caption = paste0("The ", n_zero, " occupations recorded with zero unemployed, ",
                     "classified by the rate ceiling that rounding to ", round_base,
                     " implies (", round_base/2 - 1, "/lf)")
  )
The 130 occupations recorded with zero unemployed, classified by the rate ceiling that rounding to 250 implies (124/lf)
band Occupations
Provably below 0.5% 45
Provably below 2% 57
Genuinely uncertain 28

The reading is stark. Only the last group is genuinely uncertain — small enough that a recorded zero really could sit anywhere from zero up to a few percent. The other two groups, 102 occupations in all, are provably low-rate: rounding disguised occupations we can prove have low unemployment as occupations we supposedly know nothing about. 45 of them are provably below half a percent.

This sharpens the case against ranking raw rates, and it also bounds what the shrinkage below can achieve. Shrinkage corrects for sampling noise, and it does that well. It does not know about the rounding, and the consequence is subtle. A large occupation with a recorded zero has a large n, so the model shrinks it only a little — but a little is not nothing, and the pull is toward the group rate, which is positive. So these cells do not stay at zero; they lift slightly off the floor. The trouble is that they lift toward the group, not toward what the rounding actually implies. Rounding tells us the true rate lies between zero and 124/lf — a sliver near zero for a large occupation. Shrinkage aims instead at the group rate, several times higher, and stops wherever n happens to arrest it. For the very largest occupations that pull is so slight it lands inside the permitted sliver by luck; for merely large ones it can overshoot, producing a fitted rate the rounding proves impossible.

So shrinkage partially mitigates these zeros as a side effect — it moves them off the artificial floor — but it targets the wrong quantity and cannot be relied on to land them correctly.

It is worth being clear that this is not a problem more data would solve, because rounding and sampling are independent sources of error. Imagine the extreme case: the survey covers every one of the 177,250 people for the full window, so there is no sampling error at all and the true unemployment count is known exactly — and then that exact count is rounded to the nearest 250 before release. A recorded zero still means only that the true count was somewhere in 0–124; you still cannot say where. Shrinkage has nothing to do in this world — every occupation is measured perfectly — yet the rounding ambiguity is entirely untouched. That is the proof that the two problems are orthogonal: the largest possible sample eliminates the one shrinkage is built for and leaves the other exactly as it was. Worse, a large sample makes shrinkage trust the number more, so it preserves the rounded zero rather than correcting it. The one consolation is that shrinkage’s pull and the needed correction happen to point the same way — both upward, off the floor — so the model under-corrects these cells rather than sending them the wrong direction.

The right repair is to use the bound directly: treat the zero as an interval (true count in 0–124) rather than a point, which is a change to the likelihood this document does not make. The point here is narrower and prior: whatever the downstream index does with these occupations, it should not treat a provably-low occupation and a genuinely-uncertain one as the same number.

What the model actually estimates

The model has three unknowns:

  • The overall (grand) log-odds of unemployment — the intercept. One number describing the centre of the whole system.
  • How much major groups vary around that centre — one variance.
  • How much occupations vary around their own major group — a second variance.

Those two variances set the geometry: how far apart the major groups sit, and how much occupations scatter inside them. The model estimates both from the data rather than being told.

What they do not do is decide how much any given occupation moves. That is governed by its sample size. The variances say where the pull points; the sample size says how hard it pulls.

The trade-off for any single occupation is roughly:

Pull toward the group is large when the occupation’s own sample is small, or when the occupations in its group are all similar to each other. Pull is small when its own sample is large, or when its group’s occupations are genuinely spread out.

Everything below is a working out of that sentence.

Teacher reallocation

NOC 41229 is a residual “other” category for teachers not used in the LMO. Its counts are distributed across 41220 and 41221 in proportion to their sizes, separately for the employed and the unemployed.

Recovering sample size

Here is the difficulty at the centre of the exercise. Shrinkage needs to know how many survey respondents stand behind each occupation, because that is what determines how noisy the raw rate is. But the tabulations available to us give weighted population counts — estimates of how many British Columbians are in each occupation, not how many people were interviewed. However, the LFS public use microdata files are one record per respondent. Counting the British Columbia labour-force records gives the monthly total directly:

Month BC labour-force records Sample design
January 2024 8,249 2015
January 2025 8,754 2015
December 2025 8,058 2025

About 8,350 per month, or roughly 200,000 respondent-months across the 24 months of 2024–25.

The three months were chosen to bracket the sample redesign that was phased in between April and September 2025. They show no systematic shift across it: the December 2025 figure, fully under the new design, sits below both Januaries under the old one. Whatever changed nationally did not change BC’s respondent count. The 9% spread across the three is month-to-month variation, and it is the practical precision of this figure — nothing finer is meaningful.

The distribution: still an assumption

The total is measured. How it distributes across occupations is not, and cannot be — that is precisely what RTRA declines to return.

We assume that an occupation’s share of the survey’s respondents equals its share of employment. An occupation holding 1% of BC’s jobs is taken to be 1% of the people the survey talked to.

\[n_{\text{occ}} = \frac{\text{count}_{\text{occ}}}{\sum \text{count}} \times \text{respondent-months}\]

Where that assumption comes from, and when it fails

The LFS does not sample occupations. It samples dwellings, stratified geographically — economic regions, census metropolitan areas, employment insurance regions. Occupation is observed after the fact, from whoever happens to live in the selected dwellings. Nothing in the design targets occupations at all.

That is what makes proportionality a reasonable default: with no occupational targeting, an occupation’s respondents are roughly its share of the population that got sampled. But it fails whenever an occupation’s geography differs from the population’s:

  • Occupations concentrated in oversampled areas have more respondents than their employment share implies. British Columbia carries provincially-funded oversamples targeting specific sub-populations.
  • Occupations concentrated in Vancouver have fewer: a large population sampled at a low rate to meet CMA quality targets.

For those occupations we are misjudging how much evidence exists — under-crediting the first group and over-crediting the second. The correction requires respondent-level weights, which RTRA will not return.

This is now the only inferred step in the chain, and it is worth being clear that the sensitivity analysis below does not test it. Varying the total tests whether the overall scale matters. Proportionality is about how that total is divided, and nothing here can bound it.

NoteSources
  • Respondent counts: Statistics Canada, Labour Force Survey public use microdata files, January 2024, January 2025, and December 2025.
  • Sample design, redesign phase-in dates, and population coverage: Statistics Canada, Guide to the Labour Force Survey, 2025. Catalogue no. 71-543-G, released 4 April 2025, Sections 4 and 5.

n and y are in respondent-months: one person observed in one month. An occupation with n = 600 might be 50 people observed across twelve months, or 600 people observed once. The model treats these identically.

WarningWhat remains uncertain

Measuring the total closed most of the questions that used to sit here. Two things about the level survive, and both are small:

  1. Three months is not twenty-four. The counts vary by about 9% across the months sampled, so 200,000 is a mean with a margin around it. The sensitivity analysis below spans considerably more than that.

  2. The design change could still matter, though the evidence says not. The sample redesign phased in between April and September 2025 falls inside the window. A reconstruction from the Guide predicted BC’s respondent count would rise with the larger national sample. The December 2025 microdata says it did not — that month, fully under the new design, has fewer BC labour-force records than either January. One month is not proof, but it is evidence against, and it is more than the reconstruction had.

Both affect the level: how much total evidence the model thinks it has. Getting the level wrong scales every n by a common factor, so the ordering of shrinkage survives; what changes is the amount. Too large a respondent count makes the data look more informative than it is, and the model shrinks too little.

The assumption that does not appear on this list, because the sensitivity analysis cannot touch it, is proportionality — how the measured total distributes across occupations. It is discussed above and it is now the weakest link.

WarningA second problem, with a known direction

Respondent-months are not independent observations. The LFS uses a rotating panel in which dwellings remain in the sample for six consecutive months, producing a five-sixths month-to-month overlap. One person unemployed for six months contributes six correlated respondent-months, not six independent pieces of evidence.

The starkest case: to limit response burden, answers from people aged 70 and older are carried forward from the initial interview to the subsequent five months. Those are not merely correlated observations — they are the same response, repeated. Our n counts each of them six times. (Most are outside the labour force, so the effect on unemployment rates is limited, but it illustrates the point exactly.)

This makes n overstate the true information content, which makes the model shrink less than it should. Unlike the uncertainty in the total — which is now measured, and small — this one has a known sign: the panel structure always inflates the apparent evidence, never deflates it.

The standard correction is an effective sample size, n_eff = (Σw)²/Σw², which also handles the fact that weights vary across occupations. Computing it requires the sum of squared weights per cell, which RTRA will not return. It remains out of reach.

So on this count the estimates are conservative: closer to the raw rates than a fully correct treatment would put them. This is the largest known problem with the analysis, and unlike the respondent count it cannot be fixed by counting something.

How much does the respondent count actually matter?

Below we test the sensitivity of the amount of shrinkage with respect to our assumption of 200,000 respondent-months. The grid runs from 150,000 to 250,000 — roughly ±25% around 200,000, which is far wider than the 9% month-to-month variation in the counts themselves.

Code
sweep_breaks <- c(0.005, 0.01, 0.02, 0.05, 0.10, 0.20)

sweep |>
  ggplot(aes(rm, qlogis(urate_eb), group = noc_5)) +
  geom_line(alpha = 0.15) +
  geom_vline(xintercept = resp_months, colour = "firebrick", linetype = "dashed") +
  scale_y_continuous(breaks = qlogis(sweep_breaks),
                     labels = scales::percent(sweep_breaks, accuracy = 0.1)) +
  labs(
    x = "BC respondent-months assumed",
    y = "Shrunken unemployment rate",
    title = "Sensitivity of shrunken rates to the assumed respondent count",
    subtitle = paste0("Red line = adopted value (",
                      scales::comma(resp_months), ", measured).\n",
                      "Grid spans +/-25% around the measured figure.")
  ) +
  theme_minimal()

Each line is one occupation. Flat lines are insensitive to the assumed respondent count; steep ones are not.
Code
# The occupations of interest are those with the LOWEST unemployment rates:
# this analysis feeds an opportunity measure, not a distress measure. cut_n
# is the index's own selection size, set in the constants chunk above.
cutoff_stability |>
  knitr::kable(
    caption = paste0("Stability of the ", cut_n,
                     " lowest-unemployment occupations across the respondent-count range")
  )
Stability of the 125 lowest-unemployment occupations across the respondent-count range
Ever selected Selected at every value Movers
139 111 28
Code
ranked |>
  filter(times_in > 0, times_in < length(rm_grid)) |>
  arrange(desc(swing)) |>
  slice_head(n = 10) |>
  mutate(`Raw rate` = scales::percent(urate_raw, accuracy = 0.1)) |>
  select(Occupation = class_title, `Raw rate`,
         `Best rank` = best_rk, `Worst rank` = worst_rk, `Swing` = swing) |>
  knitr::kable(
    digits = 0,
    caption = "Occupations that move in and out of the selection as the assumed respondent count varies"
  )
Occupations that move in and out of the selection as the assumed respondent count varies
Occupation Raw rate Best rank Worst rank Swing
Pest controllers and fumigators 0.0% 123 181 58
Postal services representatives 0.0% 112 169 57
Theatre, fashion, exhibit and other creative designers 0.0% 92 141 49
Medical laboratory technologists 1.4% 100 148 48
Technical occupations in geomatics and meteorology 0.0% 103 149 46
Elementary school and kindergarten teachers 1.3% 94 138 44
Assemblers and inspectors of other wood products 0.0% 96 140 44
Retail and wholesale trade managers 1.3% 92 135 43
Nurse aides, orderlies and patient service associates 1.4% 120 162 42
Assemblers and inspectors, electrical appliance, apparatus and equipment manufacturing 0.0% 93 134 41

That second table looks alarming, and taken alone it would be. The largest movers swing by dozens of places — comfortably inside the selection at one respondent count, clearly outside it at another. Those are not boundary cases.

But the table is misleading about why, and the reason matters more than the finding. Look at the raw-rate column: most of the big movers report zero.

It is the zeros that move

Code
ranked |>
  mutate(`Raw rate` = if_else(urate_raw == 0, "Zero", "Non-zero")) |>
  summarize(
    `Occupations`  = dplyr::n(),
    `Median swing` = median(swing),
    `Largest swing`= max(swing),
    .by = `Raw rate`
  ) |>
  knitr::kable(
    digits = 1,
    caption = "Rank movement across the respondent-count range, split by whether the occupation recorded any unemployed respondents"
  )
Rank movement across the respondent-count range, split by whether the occupation recorded any unemployed respondents
Raw rate Occupations Median swing Largest swing
Zero 130 19.5 85
Non-zero 383 7.0 71

Occupations with a raw rate of zero swing far more than the rest. That is not a malfunction — it is the design working. A zero carries only weak information about where its occupation sits: observing no one unemployed is genuine evidence of a low rate, but on a thin sample it is not much evidence, so the estimate leans heavily on the group. How heavily is exactly what the respondent count controls — raise it and the same zero counts for more, lower it and the model leans harder on the group. Lower the count, then, and the evidence looks thinner, the zeros get pulled further toward their group rates, and they rise.

The consequence is the part worth understanding. Around 25% of all occupations are zeros. When they rise, everything they pass falls in rank — without moving at all. Two occupations make the mechanism visible:

Code
illus_tbl |>
  mutate(
    `Respondent-mo.` = round(n),
    `Raw rate`       = scales::percent(urate_raw, accuracy = 0.01),
    `Fitted rate`    = scales::percent(urate_eb, accuracy = 0.01),
    Rank             = rk
  ) |>
  select(Occupation = class_title, `Respondent-months assumed` = rm,
         `Respondent-mo.`, `Raw rate`, `Fitted rate`, Rank) |>
  arrange(Occupation, `Respondent-months assumed`) |>
  knitr::kable(caption = "Two occupations across the respondent-count range: one moving, one being moved past.")
Two occupations across the respondent-count range: one moving, one being moved past.
Occupation Respondent-months assumed Respondent-mo. Raw rate Fitted rate Rank
Oil and gas well drilling and related workers and services operators 150000 25 0.00% 3.34% 326
Oil and gas well drilling and related workers and services operators 175000 29 0.00% 2.89% 297
Oil and gas well drilling and related workers and services operators 200000 34 0.00% 2.59% 275
Oil and gas well drilling and related workers and services operators 225000 38 0.00% 2.38% 260
Oil and gas well drilling and related workers and services operators 250000 42 0.00% 2.18% 241
Retail and wholesale trade managers 150000 2829 1.27% 1.29% 92
Retail and wholesale trade managers 175000 3300 1.27% 1.29% 109
Retail and wholesale trade managers 200000 3772 1.27% 1.28% 117
Retail and wholesale trade managers 225000 4243 1.27% 1.28% 126
Retail and wholesale trade managers 250000 4715 1.27% 1.28% 135

Oil and gas well drilling and related workers and services operators reported no unemployed respondents at all, on a thin sample. Its fitted rate runs from 2.18% to 3.34% across the range — it is genuinely moving, because on so thin a sample its own zero counts for little and the estimate leans mostly on the group, and the assumed respondent count decides how hard.

Retail and wholesale trade managers has a raw rate of 1.27% on about 3772 respondent-months. Its fitted rate moves from 1.28% to 1.29% — for practical purposes immovable. Its rank swings 43 places anyway, because the zeros are climbing past it.

Fitting the model

Reading the formula:

  • urate — the outcome is the occupation’s unemployment rate: a proportion between 0 and 1.
  • weights = n — how much that proportion counts for, in respondent-months. This is how the model knows an occupation with n = 30 carries less weight than one with n = 3000. The sample size is not a side note; it is what makes shrinkage work.
  • ~ 1 — the only fixed effect is the intercept: the overall rate.
  • (1 | noc_2 / noc_5) — the hierarchy. This expands to (1 | noc_2) + (1 | noc_2:noc_5): a random intercept for each major group, and a random intercept for each occupation within its major group. The / is what makes this nested rather than crossed, and it is the whole point.
  • family = binomial — the outcome is a proportion bounded by 0 and 1, so the model works in log-odds and converts back at the end.

The two variance components

Code
vc |>
  select(grp, vcov, sdcor) |>
  mutate(
    level = case_when(
      grp == "noc_2"        ~ "Between major groups",
      grp == "noc_5:noc_2"  ~ "Between occupations, within major group",
      TRUE                  ~ grp
    )
  ) |>
  select(level, variance = vcov, sd = sdcor) |>
  knitr::kable(digits = 4, caption = "Estimated variance components (log-odds scale)")
Estimated variance components (log-odds scale)
level variance sd
Between occupations, within major group 1.0001 1.0001
Between major groups 0.3352 0.5790

These two numbers describe how spread out the hierarchy is at each level. They are on the log-odds scale, which is hard to read directly, so the table below translates them into rates.

Code
tibble(
  Level = c("A typical major group, relative to the overall rate",
            "A typical occupation, relative to its major group"),
  `One SD below` = plogis(centre + c(-sd_noc2, -sd_noc5)),
  `Centre`       = plogis(centre + c(0, 0)),
  `One SD above` = plogis(centre + c(sd_noc2, sd_noc5))
) |>
  mutate(across(where(is.numeric), ~ scales::percent(.x, accuracy = 0.1))) |>
  knitr::kable(
    caption = "The variance components translated into unemployment rates, for readability"
  )
The variance components translated into unemployment rates, for readability
Level One SD below Centre One SD above
A typical major group, relative to the overall rate 1.4% 2.5% 4.3%
A typical occupation, relative to its major group 0.9% 2.5% 6.4%

Both rows are now centred on the same overall rate and show one standard deviation either side, so they can be compared directly.

The first row says major groups are not interchangeable: a group one standard deviation from the centre sits at a materially different unemployment rate from one on the other side. There is real signal in knowing which group an occupation belongs to — which is what makes borrowing from the group worth doing at all.

The second row is wider. The within-group variance (1) exceeds the between-group variance (0.335): occupations inside a major group differ from each other more than the groups differ from one another. Knowing an occupation’s major group tells you something real, but it is a long way from telling you everything.

NoteWhat these numbers do and don’t determine

They set the targets and the scale: where an occupation with no data of its own will land, and how far apart the levels sit.

They do not determine how much any given occupation gets shrunk. That is governed by sample size — an occupation with 3,000 respondent-months keeps its own rate almost regardless of these variances, and one with 30 gets pulled toward its group almost regardless. The variance components say where the pull points; n says how hard it pulls.

It is also tempting to read the larger within-group variance as meaning the major groups are barely distinguishable and the model is leaning on the overall rate. That is not what it means. Compare the two rows above: the first shows real separation between groups, and the comparison between the rows says something narrower than “groups don’t matter” — occupations inside a group are heterogeneous, so a group’s rate is a useful but imprecise predictor of any one occupation in it. The borrowing is real, and it is noisy.

Seeing the two levels at work

The model borrows strength at two levels: occupations lean on their major group, and major groups lean on the overall rate. The variance components above are the evidence that both are happening. This section shows the borrowing where it is directly visible — at the occupation level, where each estimate lies between the occupation’s own data and its major group’s rate.

The overall rate

Code
tibble(
  quantity = c("Grand mean (log-odds)", "Grand mean (as a rate)"),
  value    = c(round(grand_logodds, 4), scales::percent(grand_rate, accuracy = 0.01))
) |>
  knitr::kable(caption = "The centre of the whole system")
The centre of the whole system
quantity value
Grand mean (log-odds) -3.6798
Grand mean (as a rate) 2.46%
Note

This is not the same as the provincial unemployment rate, and shouldn’t be expected to match it. It is the centre of the distribution of occupations, with each occupation contributing on its own terms rather than in proportion to how many workers it contains. A small occupation and a huge one both inform it.

The major-group rates

Each major group’s fitted rate is the overall log-odds plus that group’s estimated offset, converted back to a rate. These are the targets that detailed occupations are pulled toward, and they are what the rest of this section uses.

Code
middle_layer |>
  select(
    `Major group`    = noc_2,
    `Occupations`    = n_occupations,
    `Respondent-mo.` = n_noc2,
    `Fitted rate`    = urate_eb_noc2
  ) |>
  mutate(`Fitted rate` = scales::percent(`Fitted rate`, accuracy = 0.01)) |>
  knitr::kable(
    caption = "Fitted unemployment rate for each major group — the target that its occupations are pulled toward."
  )
Fitted unemployment rate for each major group — the target that its occupations are pulled toward.
Major group Occupations Respondent-mo. Fitted rate
65 16 16500.16467 3.51%
72 53 15638.37962 2.67%
64 18 13805.02799 3.38%
41 23 13474.31112 1.69%
21 36 11917.47173 1.96%
11 8 9910.52805 1.81%
31 19 8229.49830 0.75%
14 18 8016.11044 3.47%
73 18 7681.96289 3.63%
63 10 7334.77879 2.36%
13 8 7022.58755 2.70%
62 12 6724.11900 2.21%
12 17 6724.11900 1.78%
75 9 6223.92689 4.16%
60 5 6097.67812 1.85%
42 7 5972.11549 2.10%
22 27 5205.01702 3.43%
33 5 4328.82314 2.13%
32 18 3665.33099 1.05%
70 5 3663.27259 1.86%
94 40 3206.30695 3.10%
10 9 3173.37249 2.16%
43 6 2776.10056 1.75%
52 9 2425.48578 5.14%
40 9 2118.78362 1.69%
54 1 1917.74619 3.08%
51 11 1917.74619 2.20%
20 3 1836.78230 2.74%
85 9 1695.43858 7.83%
74 9 1656.32891 2.54%
53 12 1534.19695 3.67%
92 13 1149.27544 1.81%
95 9 1080.66198 5.62%
90 2 867.27412 2.28%
44 3 852.86530 3.03%
83 5 693.68207 4.73%
82 5 674.47030 3.22%
80 4 578.41146 1.86%
30 1 408.93622 2.10%
00 6 377.37403 1.56%
84 6 357.47612 5.97%
50 3 244.95005 2.94%
93 4 157.81096 2.93%
45 1 89.19750 3.19%
55 1 74.10254 3.65%

The occupation layer

Now the same picture one level down. Each occupation’s fitted rate is the overall log-odds, plus its major group’s offset, plus its own offset within that group.

TipThe decomposition, made explicit

For any occupation, the fitted log-odds is a sum of three pieces:

Code
bind_rows(
  eb |> filter(Unemployed == 0) |> arrange(n) |> head(4),
  eb |> filter(Unemployed == 0) |> arrange(desc(n)) |> head(4)
) |>
  mutate(
    `Overall`        = grand_logodds,
    `+ Group offset` = qlogis(urate_eb_noc2) - grand_logodds,
    `+ Own offset`   = offset_noc5,
    `= Fitted`       = qlogis(urate_eb)
  ) |>
  select(
    Occupation = class_title,
    `Resp-mo.` = n,
    Overall, `+ Group offset`, `+ Own offset`, `= Fitted`,
    `Raw rate` = urate_raw,
    `Fitted rate` = urate_eb
  ) |>
  mutate(
    `Resp-mo.` = round(`Resp-mo.`),
    across(c(`Raw rate`, `Fitted rate`),
           ~ scales::percent(.x, accuracy = 0.1))
  ) |>
  knitr::kable(
    digits = 3,
    caption = "Occupations with a raw rate of zero: the four smallest and four largest by sample. Log-odds columns sum left to right."
  )
Occupations with a raw rate of zero: the four smallest and four largest by sample. Log-odds columns sum left to right.
Occupation Resp-mo. Overall + Group offset + Own offset = Fitted Raw rate Fitted rate
Motor vehicle assemblers, inspectors and testers 1 -3.68 0.238 -0.021 -3.463 0.0% 3.0%
Metallurgical and materials engineers 2 -3.68 -0.230 -0.039 -3.949 0.0% 1.9%
Machine fitters 2 -3.68 0.083 -0.052 -3.649 0.0% 2.5%
Aircraft assemblers and aircraft assembly inspectors 2 -3.68 0.179 -0.057 -3.558 0.0% 2.8%
Pharmacists 486 -3.68 -1.207 -1.154 -6.041 0.0% 0.2%
Massage therapists 451 -3.68 -0.862 -1.303 -5.845 0.0% 0.3%
Education policy researchers, consultants and program officers 396 -3.68 -0.386 -1.504 -5.570 0.0% 0.4%
Dental hygienists and dental therapists 388 -3.68 -0.862 -1.218 -5.760 0.0% 0.3%

Every occupation in this table has the same raw rate — zero. Under a naive ranking they would be indistinguishable, all tied for lowest.

The four at the top have almost no data. Their own offset is close to zero, so the model declines to move them away from their group: their fitted rate is essentially their major group’s rate, and their zero has bought them nothing. The four at the bottom observed no unemployment across a substantial sample. Their own offset is doing real work, and their fitted rate sits well below their group’s.

Same raw number, different amounts of evidence, different answers.

The main result

Code
main_breaks <- c(0.005, 0.01, 0.02, 0.05, 0.10, 0.20, 0.40)

# Single common position for every zero, below the smallest non-zero raw rate.
# See the group-32 plot above for why this is not 0.5/n.
main_floor <- eb |> filter(urate_raw > 0) |> pull(urate_raw) |> min() / 2

plt_dat <- eb |>
  mutate(
    zero_raw = urate_raw <= 0,
    lo_raw   = qlogis(if_else(zero_raw, main_floor, pmin(urate_raw, 1 - 1e-6))),
    lo_eb    = qlogis(urate_eb)
  )

plt <- ggplot(plt_dat, aes(
  lo_raw,
  lo_eb,
  size = n,
  text = paste(
    class_title,
    "<br>Survey sample: ", round(n,1), " respondent-months",
    "<br>Major group: ", noc_2,
    "<br>Raw Unemployment Rate: ", scales::percent(urate_raw, accuracy = .1),
    "<br>Group Rate: ", scales::percent(urate_eb_noc2, accuracy = .1),
    "<br>Shrunken Unemployment Rate: ", scales::percent(urate_eb, accuracy = .1)
  )
)) +
  geom_abline(slope = 1, intercept = 0, lwd = 2, colour = "white") +
  geom_point(aes(colour = zero_raw), alpha = .25) +
  scale_colour_manual(values = c(`FALSE` = "black", `TRUE` = "firebrick"),
                      guide = "none") +
  scale_x_continuous(breaks = qlogis(main_breaks),
                     labels = scales::percent(main_breaks, accuracy = 0.1)) +
  scale_y_continuous(breaks = qlogis(main_breaks),
                     labels = scales::percent(main_breaks, accuracy = 0.1)) +
  labs(
    size = "Respondent-months",
    x = "Raw Unemployment Rate",
    y = "Shrunken Unemployment Rate",
    title = paste0("Unemployment rates by Occupation ", min(years_use), ":", max(years_use),
                   ":\nRaw vs Shrunken via Hierarchical Empirical Bayes")
  )

plotly::ggplotly(plt, tooltip = "text")

Raw vs. shrunken occupational unemployment rates. The white line is where the two agree; distance from it is the size of the correction. Both axes are log-odds, labelled as rates.

The white line is where the two rates agree. Large points — occupations with a lot of survey data — cluster along it, because the model has no reason to move them. Small points sit further from it: their raw rates have been overruled.

Look at the left edge. The occupations shown in red had a raw unemployment rate of exactly zero, and under a naive ranking every one of them would tie for the lowest rate in British Columbia. They are stacked in a vertical line, because they all share the same horizontal position — but their fitted rates, read off the vertical axis, are spread over a wide range. The model has sorted them by how much their zero was worth.

Note what the pull is toward. It is not toward a single global rate. The occupations with the least data — the red points at the left edge, and the small points furthest from the white line — do not converge on one height. They settle onto their own major groups’ rates, which sit at different levels. That vertical spread among the no-evidence occupations is the middle layer earning its keep.

Note

Zero has no position on a log-odds axis, so the red points are all drawn at a single common position at the left edge, below the lowest non-zero rate in the data. That position is a convention and carries no meaning — every one of these occupations reported the same raw number.

Their fitted values are unaffected: the model handles zero cells correctly, and their position on the vertical axis is real. Only the horizontal position is invented, and it is invented identically for all of them, which is the point.

What to be careful about

  • Proportionality and the panel structure. The survey’s total respondent count is measured from the public microdata, but how those respondents distribute across occupations is assumed to follow employment shares — and the LFS samples geography, not occupation, so occupations with unusual geographic concentration are misjudged. Separately, the panel structure is uncorrected: the effective-sample-size calculation needs quantities RTRA will not return. The panel problem has a known direction — it makes n overstate the information available, so the model shrinks less than it should. Neither is bounded by the sensitivity analysis, which varies the total rather than its distribution.
  • Rounding is not corrected. Published counts are rounded to the nearest 250, which forces many occupations — including large, well-sampled ones — to a recorded zero. Shrinkage addresses sampling noise, not this: a large occupation with a recorded zero has a large n, so the model nudges it only slightly off zero, and toward its group rate rather than toward the low interval the rounding actually implies. Treating those cells as intervals rather than points is the correct repair and is left for future work. Its effect is confined to occupations at or near a recorded zero.
  • Ranks are less robust than rates. The sensitivity analysis found the fitted rates stable across the plausible range but the rankings considerably less so — and the movement is concentrated in the occupations that reported no unemployed respondents. Because those are a large share of the distribution, their movement displaces well-measured occupations that have not moved at all. Downstream use should prefer scaled rates over rank-based or percentile-based transformations of them.
  • Empirical Bayes ignores its own uncertainty. The two variance components are estimated, but the shrinkage treats them as if they were known exactly. The resulting intervals are slightly too narrow. With ~45 major groups and ~500 occupations this is a modest effect, but it is real.
  • Shrunken rates are not raw rates and don’t aggregate like them. They are the best estimate of each occupation’s underlying rate. Summing them back up will not reproduce published totals, and shouldn’t be expected to.
  • The model has no covariates. The only structure it knows about is the NOC hierarchy — nothing about region, age, industry, or education. Two occupations in the same major group are treated as exchangeable, which is a strong claim for some groups.

Extending the method: The ratio of EI claimants to employment.

A rate bounded in [0, 1] with the same hierarchical structure — related occupations face related conditions — so the identical model would apply: shrink each occupation’s ratio toward its group, weighted by cell size.

Two things differ from the unemployment case, and together they suggest the payoff would be smaller.

First, the noise sits in a different place. The EI numerator is administrative — a count of claimants, measured essentially without sampling error — while the denominator, employment, comes from the LFS and carries all the noise. That is the reverse of the unemployment rate, where the numerator (a rare count of unemployed people) was the fragile part. So the zero-floor pathology that motivated this whole exercise would be absent, but shrinkage would still smooth the estimates for small occupations, but it would be smoothing mild noise rather than rescuing a broken statistic.

Second, the cell is employment, not the labour force, so the occupational weights change. The respondent count that governs shrinkage would be employed respondent-months, allocated by each occupation’s employment share — not the labour-force count and shares used above. The two differ most for occupations whose unemployment rate is far from average, which for a claimant-based indicator may be exactly the occupations of interest. It is the same measured total from the microdata, filtered to the employed and re-allocated; not a reuse of the figure derived above.

The honest expectation is a real but modest improvement: the EI ratio is a better-behaved statistic than the raw unemployment rate to begin with, so there is less for shrinkage to fix. The larger unfixed problem for that indicator is not sampling noise at all but whether its numerator and denominator describe the same occupational population — EI claims are coded from administrative records, LFS employment from self-report — which no amount of shrinkage addresses.