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)