Code
library(tidyverse)
library(lubridate)
library(brms)
library(posterior)Looks like resident doctors are going back on strike next month. This latest round of industrial action is part of a long running campaign for full pay restoration.
Leaving aside the politics and ethics of the industrial action, it has been noted that the strikes seem to fall on periods of school holidays. But is this really the case?
We can start by listing all the periods of industrial action in England since these strikes started in 2023:
strikes <- tribble(
~start, ~end,
"2023-03-13", "2023-03-15",
"2023-04-11", "2023-04-14",
"2023-06-14", "2023-06-16",
"2023-07-13", "2023-07-17",
"2023-08-11", "2023-08-14",
"2023-09-20", "2023-09-22",
"2023-10-02", "2023-10-04",
"2023-12-20", "2023-12-22",
"2024-01-03", "2024-01-08",
"2024-02-24", "2024-02-28",
"2024-06-27", "2024-07-01",
"2025-07-25", "2025-07-29",
"2025-11-14", "2025-11-18",
"2025-12-17", "2025-12-21",
"2026-04-07", "2026-04-12"
) %>%
mutate(
start = ymd(start),
end = ymd(end),
length = as.integer(end - start) + 1L
)We can also pull school holidays over this time period for my local authority:
holidays <- tribble(
~start, ~end, ~label,
"2023-04-03", "2023-04-14", "Spring break 2023",
"2023-05-29", "2023-06-02", "Half-term May 2023",
"2023-07-22", "2023-09-03", "Summer 2023",
"2023-10-23", "2023-11-05", "Half-term Oct 2023",
"2023-12-23", "2024-01-03", "Christmas 2023/24",
"2024-02-12", "2024-02-18", "Half-term Feb 2024",
"2024-03-29", "2024-04-14", "Spring break 2024",
"2024-05-25", "2024-06-02", "Half-term May 2024",
"2024-07-27", "2024-09-01", "Summer 2024",
"2024-10-19", "2024-11-03", "Half-term Oct 2024",
"2024-12-21", "2025-01-05", "Christmas 2024/25",
"2025-02-15", "2025-02-23", "Half-term Feb 2025",
"2025-03-29", "2025-04-13", "Spring break 2025",
"2025-05-24", "2025-06-01", "Half-term May 2025",
"2025-07-28", "2025-08-31", "Summer 2025",
"2025-10-20", "2025-10-31", "Half-term Oct 2025",
"2025-12-20", "2026-01-04", "Christmas 2025/26",
"2026-02-16", "2026-02-20", "Half-term Feb 2026",
"2026-03-28", "2026-04-12", "Spring break 2026"
) %>%
mutate(
start = ymd(start),
end = ymd(end)
)Simple analysis (and why it’s statistically wrong)
The simplest analysis would be to just look at the number of strike days that fall on holidays:
Code
# Return a vector of all dates within each block/holiday table
expand_intervals_to_days <- function(df) {
map2(df$start, df$end, seq.Date, by = "day") %>%
unlist() %>%
as.Date(origin = "1970-01-01") %>%
unique() %>%
sort()
}
# Build a full study calendar with logical strike/holiday indicators
build_calendar <- function(strikes_df, holidays_df) {
study_start <- min(strikes_df$start)
study_end <- max(strikes_df$end)
all_days <- tibble(
date = seq.Date(study_start, study_end, by = "day")
)
strike_days <- expand_intervals_to_days(strikes_df)
holiday_days <- holidays_df %>%
filter(end >= study_start, start <= study_end) %>%
mutate(
start = pmax(start, study_start),
end = pmin(end, study_end)
) %>%
expand_intervals_to_days()
all_days %>%
mutate(
is_strike = date %in% strike_days,
is_holiday = date %in% holiday_days
)
}
# Return start/end/length metadata and calendar vectors for a study period
prepare_study_objects <- function(strikes_df, holidays_df) {
calendar_df <- build_calendar(strikes_df, holidays_df)
study_start <- min(calendar_df$date)
study_end <- max(calendar_df$date)
list(
study_start = study_start,
study_end = study_end,
calendar_df = calendar_df,
holiday_vec = calendar_df$is_holiday,
strike_vec = calendar_df$is_strike,
lengths = strikes_df$length,
n_blocks = nrow(strikes_df),
total_days = nrow(calendar_df),
total_strike_days = sum(strikes_df$length),
free_days = nrow(calendar_df) - sum(strikes_df$length)
)
}# Expand dates into a calendar of each day tagged as a strike and/or school holiday
naive_calendar <- build_calendar(strikes, holidays)
# Basic counts
naive_summary <- naive_calendar %>%
summarise(
study_days = n(),
strike_days = sum(is_strike),
holiday_days = sum(is_holiday),
strike_holiday_days = sum(is_strike & is_holiday),
pct_days_on_strike = 100 * strike_days / study_days,
pct_days_in_holiday = 100 * holiday_days / study_days,
pct_strike_days_in_holiday = 100 * strike_holiday_days / strike_days
)So during the entire 1127 day period that the doctors have been involved in the dispute, they’ve been on strike for 65 days (6%). Of these days, 29% have been during the school holiday… but 27% of the total days since the start have been school holidays anyway, so is this really a surprisng finding?
One approach might be an exact binomial test. Under a naive null, each strike day has probability equal to the background proportion of holiday days in the study window:
binom.test(
x = naive_summary$strike_holiday_days,
n = naive_summary$strike_days,
p = naive_summary$holiday_days / naive_summary$study_days,
alternative = "greater"
)
Exact binomial test
data: naive_summary$strike_holiday_days and naive_summary$strike_days
number of successes = 19, number of trials = 65, p-value = 0.4048
alternative hypothesis: true probability of success is greater than 0.2724046
95 percent confidence interval:
0.2006267 1.0000000
sample estimates:
probability of success
0.2923077 This is a non-significant finding (p = 0.405) - there doesn’t seem to be a bias towards school holidays, just lots of school holidays!
However, this approach (or even the more simple proportion test) is not really the right test here, because it treats each strike day as if it were an independent random observation.
In reality, strike days are clustered into a small number of multi-day strike episodes, so a \(n\)-day strike is only one scheduling decision, not \(n\) independent ones. That breaks the independence assumption and makes the data look more informative than it really is, which can exaggerate statistical significance. The test also ignores the fact that both strikes and school holidays occur in blocks on a calendar, so the key question is not simply whether “some strike days fell in holidays,” but whether the timing of whole strike blocks aligned with holiday periods more than expected by chance.
Block-preserving permutation test
A better approach is a permutation test. The way this works is that we keep the observed strike blocks at their real lengths, then randomly relocate those 15 blocks across the study window many times and ask how often you get as much school-holiday overlap as actually observed.
This tests “more overlap than expected by chance?” without pretending each strike day was independent. It also respects clustering, both of the strikes and the school holidays, and avoids forcing a wrong sampling model (as we don’t actually have a natural parametric model saying strike dates arise from a binomial distribution as we used in our naive analysis, or anything to suggest a Poisson or hypergeometric process for that matter).
The null hypothesis is that there’s no association - so the test is essentially asking if given the actual way strikes are called, is the observed overlap with school holidays more than we would expect by chance?
In other words, where our binomial test asked “out of 65 strike days, is 19 holiday days too many?” the block-preserving permutation test instead asks “out of 15 strike episodes with these exact lengths, is their placement unusually aligned with holidays?”
Full period of industrial action
Lets start by considering this entire round of industrial action:
Code
# Sample a random composition of free days into gaps between/around blocks
sample_gaps <- function(total, parts) {
cuts <- sort(sample.int(total + parts - 1L, parts - 1L, replace = FALSE))
cuts <- c(0L, cuts, total + parts - 1L)
diff(cuts) - 1L
}
# Calculate total holiday overlap from sampled gaps using vector indexing
overlap_from_gaps <- function(gaps, lengths, holiday_vec) {
pos <- gaps[1] + 1L
overlap <- 0L
for (i in seq_along(lengths)) {
idx <- pos:(pos + lengths[i] - 1L)
overlap <- overlap + sum(holiday_vec[idx])
pos <- pos + lengths[i] + gaps[i + 1L]
}
overlap
}
# Calculate holiday overlap within a restricted index window
overlap_from_gaps_window <- function(gaps, lengths, holiday_vec, idx_start, idx_end) {
pos <- gaps[1] + 1L
overlap <- 0L
for (i in seq_along(lengths)) {
block_start <- pos
block_end <- pos + lengths[i] - 1L
overlap_start <- max(block_start, idx_start)
overlap_end <- min(block_end, idx_end)
if (overlap_start <= overlap_end) {
overlap <- overlap + sum(holiday_vec[overlap_start:overlap_end])
}
pos <- pos + lengths[i] + gaps[i + 1L]
}
overlap
}
# Calculate strike days within a restricted index window
strike_days_from_gaps_window <- function(gaps, lengths, idx_start, idx_end) {
pos <- gaps[1] + 1L
strike_days <- 0L
for (i in seq_along(lengths)) {
block_start <- pos
block_end <- pos + lengths[i] - 1L
overlap_start <- max(block_start, idx_start)
overlap_end <- min(block_end, idx_end)
if (overlap_start <= overlap_end) {
strike_days <- strike_days + (overlap_end - overlap_start + 1L)
}
pos <- pos + lengths[i] + gaps[i + 1L]
}
strike_days
}
# Run a fast block-preserving permutation test over a whole study window
run_block_permutation_test <- function(strikes_df,
holidays_df,
n_sim = 100000,
seed = 123) {
prep <- prepare_study_objects(strikes_df, holidays_df)
observed_overlap <- sum(prep$holiday_vec[prep$strike_vec])
observed_prop <- observed_overlap / prep$total_strike_days
set.seed(seed)
sim_overlap <- replicate(
n = n_sim,
expr = {
gaps <- sample_gaps(prep$free_days, prep$n_blocks + 1L)
overlap_from_gaps(gaps, prep$lengths, prep$holiday_vec)
}
)
tibble(
study_start = prep$study_start,
study_end = prep$study_end,
n_blocks = prep$n_blocks,
strike_days = prep$total_strike_days,
observed_overlap = observed_overlap,
observed_prop = observed_prop,
expected_overlap = mean(sim_overlap),
expected_prop = mean(sim_overlap) / prep$total_strike_days,
p_value = mean(sim_overlap >= observed_overlap),
null_lo = quantile(sim_overlap, 0.025),
null_hi = quantile(sim_overlap, 0.975),
sim_overlap = list(sim_overlap)
)
}study_start <- min(strikes$start)
study_end <- max(strikes$end)
# Full period
full_result <- run_block_permutation_test(
strikes_df = strikes,
holidays_df = holidays,
n_sim = 100000,
seed = 123
)Since the start of the industrial action, strikes have overlapped 19 times (or 29%), while we’d expext this to happen 17.7 times (27%). So this isn’t significantly different to what we’d expect if the strikes were selected randomly (p = 0.43625). Or to put it graphically:
tibble(overlap = full_result$sim_overlap[[1]]) %>%
ggplot(aes(x = overlap)) +
geom_histogram(binwidth = 1, boundary = -0.5) +
geom_vline(
xintercept = full_result$observed_overlap,
linetype = "dashed",
linewidth = 1
) +
labs(
subtitle = "Holiday-overlap days under block-preserving random placement",
x = "Holiday-overlap days",
y = "Count"
) +
theme_bw()
One caveat: even this is still a simplified model. It assumes that, under the null, strike blocks could have been placed anywhere in the study period with equal plausibility. Real-world strike timing is influenced by negotiations, notice periods, weekends, exam periods, politics, and NHS operational pressures. So the test is best understood as evidence against random calendar placement, not proof of strategic intent.
Last 18 months
Interestingly, a colleague suggested that the timing of the strikes might have become more strategic in the last 18 months. Let’s repeat this analysis but focusing just on the last year and a half:
recent_start <- study_end %m-% months(18)
# Last 18 months only
strikes_last18m <- strikes %>%
filter(end >= recent_start) %>%
mutate(
start = pmax(start, recent_start),
length = as.integer(end - start) + 1L
)
last18_result <- run_block_permutation_test(
strikes_df = strikes_last18m,
holidays_df = holidays,
n_sim = 100000,
seed = 123
)This is slightly more suggestive of intentional choices, as the strikes have overlapped 10 times (or 48%), while we’d expect this to happen 6.6 times (32%). This still isn’t a statistically significantly value (p = 0.293), but plotting the null distribution shows that because there are only a few recent strike blocks, this subgroup test will have less power than the full-period analysis, making a non-significant p-value harder to interpret strongly.
tibble(overlap = last18_result$sim_overlap[[1]]) %>%
ggplot(aes(x = overlap)) +
geom_histogram(binwidth = 1, boundary = -0.5) +
geom_vline(
xintercept = last18_result$observed_overlap,
linetype = "dashed",
linewidth = 1
) +
labs(
subtitle = "Holiday-overlap days under block-preserving random placement",
x = "Holiday-overlap days",
y = "Count"
) +
theme_bw()
Compare the two periods
The best way around this is not to “force significance,” but to switch emphasis from a binary test to using a model that borrows more structure from the full dataset. Instead of splitting the data into two small datasets and testing each separately, we can test whether the last 18 months’ strikes show more overlap than the earlier period:
Code
# Run a fast recent-vs-earlier contrast permutation test using one shared simulation loop
run_recent_vs_earlier_contrast <- function(strikes_df,
holidays_df,
recent_start,
n_sim = 100000,
seed = 123) {
prep <- prepare_study_objects(strikes_df, holidays_df)
idx_recent_start <- as.integer(recent_start - prep$study_start) + 1L
idx_recent_end <- prep$total_days
idx_earlier_start <- 1L
idx_earlier_end <- idx_recent_start - 1L
observed_recent_overlap <- sum(prep$strike_vec[idx_recent_start:idx_recent_end] &
prep$holiday_vec[idx_recent_start:idx_recent_end])
observed_recent_days <- sum(prep$strike_vec[idx_recent_start:idx_recent_end])
observed_earlier_overlap <- sum(prep$strike_vec[idx_earlier_start:idx_earlier_end] &
prep$holiday_vec[idx_earlier_start:idx_earlier_end])
observed_earlier_days <- sum(prep$strike_vec[idx_earlier_start:idx_earlier_end])
observed_recent_prop <- observed_recent_overlap / observed_recent_days
observed_earlier_prop <- observed_earlier_overlap / observed_earlier_days
observed_prop_diff <- observed_recent_prop - observed_earlier_prop
observed_overlap_diff <- observed_recent_overlap - observed_earlier_overlap
set.seed(seed)
sim_results <- replicate(
n = n_sim,
simplify = "matrix",
expr = {
gaps <- sample_gaps(prep$free_days, prep$n_blocks + 1L)
sim_recent_overlap <- overlap_from_gaps_window(
gaps, prep$lengths, prep$holiday_vec, idx_recent_start, idx_recent_end
)
sim_earlier_overlap <- overlap_from_gaps_window(
gaps, prep$lengths, prep$holiday_vec, idx_earlier_start, idx_earlier_end
)
sim_recent_days <- strike_days_from_gaps_window(
gaps, prep$lengths, idx_recent_start, idx_recent_end
)
sim_earlier_days <- strike_days_from_gaps_window(
gaps, prep$lengths, idx_earlier_start, idx_earlier_end
)
sim_recent_prop <- if (sim_recent_days > 0) {
sim_recent_overlap / sim_recent_days
} else {
NA_real_
}
sim_earlier_prop <- if (sim_earlier_days > 0) {
sim_earlier_overlap / sim_earlier_days
} else {
NA_real_
}
sim_prop_diff <- if (!is.na(sim_recent_prop) && !is.na(sim_earlier_prop)) {
sim_recent_prop - sim_earlier_prop
} else {
NA_real_
}
c(
sim_recent_overlap = sim_recent_overlap,
sim_earlier_overlap = sim_earlier_overlap,
sim_recent_days = sim_recent_days,
sim_earlier_days = sim_earlier_days,
sim_recent_prop = sim_recent_prop,
sim_earlier_prop = sim_earlier_prop,
sim_prop_diff = sim_prop_diff,
sim_overlap_diff = sim_recent_overlap - sim_earlier_overlap
)
}
) %>%
t() %>%
as_tibble()
tibble(
recent_start = recent_start,
recent_end = prep$study_end,
earlier_start = prep$study_start,
earlier_end = recent_start - days(1),
observed_recent_overlap = observed_recent_overlap,
observed_earlier_overlap = observed_earlier_overlap,
observed_recent_days = observed_recent_days,
observed_earlier_days = observed_earlier_days,
observed_recent_prop = observed_recent_prop,
observed_earlier_prop = observed_earlier_prop,
observed_prop_diff = observed_prop_diff,
observed_overlap_diff = observed_overlap_diff,
expected_recent_prop = mean(sim_results$sim_recent_prop, na.rm = TRUE),
expected_earlier_prop = mean(sim_results$sim_earlier_prop, na.rm = TRUE),
expected_prop_diff = mean(sim_results$sim_prop_diff, na.rm = TRUE),
expected_overlap_diff = mean(sim_results$sim_overlap_diff, na.rm = TRUE),
p_value_prop_diff = mean(sim_results$sim_prop_diff >= observed_prop_diff, na.rm = TRUE),
p_value_overlap_diff = mean(sim_results$sim_overlap_diff >= observed_overlap_diff, na.rm = TRUE),
prop_diff_null_lo = quantile(sim_results$sim_prop_diff, 0.025, na.rm = TRUE),
prop_diff_null_hi = quantile(sim_results$sim_prop_diff, 0.975, na.rm = TRUE),
overlap_diff_null_lo = quantile(sim_results$sim_overlap_diff, 0.025, na.rm = TRUE),
overlap_diff_null_hi = quantile(sim_results$sim_overlap_diff, 0.975, na.rm = TRUE),
n_valid_prop_sims = sum(!is.na(sim_results$sim_prop_diff)),
n_invalid_prop_sims = sum(is.na(sim_results$sim_prop_diff)),
sim_results = list(sim_results)
)
}contrast_result <- run_recent_vs_earlier_contrast(
strikes_df = strikes,
holidays_df = holidays,
recent_start = recent_start,
n_sim = 100000,
seed = 123
)This suggests that we went from 20% overlap with the school holidays to 48%. It’s a suspicious increase, but remains non-significant (p = 0.105).
contrast_result$sim_results[[1]] %>%
ggplot(aes(x = sim_prop_diff)) +
geom_histogram(binwidth = 0.02, boundary = 0) +
geom_vline(
xintercept = contrast_result$observed_prop_diff,
linetype = "dashed",
linewidth = 1
) +
labs(
subtitle = "Null distribution under block-preserving random placement",
x = "Recent overlap proportion - earlier overlap proportion",
y = "Count"
) +
theme_bw()
What can we take from this? Recent strikes seem to be much more holiday-aligned than earlier strikes in the observed data. But under the random-placement null, a difference this large still happens about 10.5% of the time. So the data suggests a recent shift towards school holiday timing, but does not provide especially strong evidence against chance under this model.
The main limitation here is that the null model still assumes that apart from block lengths strike timing was otherwise unconstrained, when in reality strike dates were further constrained by notice requirements, timing of negotiations, media and political considerations, and maybe even the timing of post-graduate exams and rotations.
Bayesian approach
When things look suspicious, but the frequentist finding is “not significant”, we can instead use a Bayesian approach to flip things from the probabilty of the finding being surprising given the null hypothesis (\(P(\text{data}|H_0)\)) to what’s the probability that there’s been a delibrate shift, given the data (\(P(H|\text{data})\)).
Code
# Create block-level holiday counts for Bayesian models
make_block_level_df <- function(strikes_df, holidays_df, recent_start) {
study_start <- min(strikes_df$start)
study_end <- max(strikes_df$end)
holiday_days <- holidays_df %>%
filter(end >= study_start, start <= study_end) %>%
mutate(
start = pmax(start, study_start),
end = pmin(end, study_end)
) %>%
expand_intervals_to_days()
strikes_df %>%
mutate(
block_id = row_number(),
recent = if_else(end >= recent_start, 1L, 0L)
) %>%
rowwise() %>%
mutate(
block_days = list(seq.Date(start, end, by = "day")),
n_days = length(block_days),
n_holiday = sum(block_days %in% holiday_days),
prop_holiday = n_holiday / n_days
) %>%
ungroup() %>%
select(block_id, start, end, length, recent, n_days, n_holiday, prop_holiday)
}
# Summarise a fitted brms model into interpretable posterior quantities
summarise_bayes_fit <- function(fit, label = "model") {
draws <- as_draws_df(fit)
draws %>%
transmute(
p_earlier = plogis(b_Intercept),
p_recent = plogis(b_Intercept + b_recent),
risk_diff = p_recent - p_earlier,
odds_ratio = exp(b_recent),
recent_gt_earlier = p_recent > p_earlier
) %>%
summarise(
model = label,
p_earlier_mean = mean(p_earlier),
p_earlier_lo = quantile(p_earlier, 0.025),
p_earlier_hi = quantile(p_earlier, 0.975),
p_recent_mean = mean(p_recent),
p_recent_lo = quantile(p_recent, 0.025),
p_recent_hi = quantile(p_recent, 0.975),
risk_diff_mean = mean(risk_diff),
risk_diff_lo = quantile(risk_diff, 0.025),
risk_diff_hi = quantile(risk_diff, 0.975),
odds_ratio_mean = mean(odds_ratio),
odds_ratio_lo = quantile(odds_ratio, 0.025),
odds_ratio_hi = quantile(odds_ratio, 0.975),
pr_recent_gt_earlier = mean(recent_gt_earlier)
)
}We’ll use a beta-binomial model because our sample is small, the strike blocks themselves are heterogeneous, and overdispersion of the strike blocks is plausible. This model allows for the strike blocks to not all be generated from one single fixed overlap probability, but for each block has its own true overlap probability, and those block-specific probabilities vary around a group mean.
More mathematically speaking, the outcome for each strike block \(i\) is the number of strike days falling within school holidays \(Y_i\), out of the total number of strike days in that block \(n_i\).
The regression model specifies \(Y_i \sim \text{Beta-Binomial}(n_i, \mu_i, \phi)\), where \(\mu_i\) is the mean probability that a strike day in block \(i\) fell during a school holiday and \(\phi\) is a precision parameter governing residual between-block heterogeneity.
The mean probability is modelled as \(\text{logit}(\mu_i) = \alpha + \beta \,\text{Recent}_i\), where \(\text{Recent}_i\) indicates whether the block occurred in the most recent 18 months.
We’ll use weakly informative priors \(\alpha \sim \mathcal{N}(0, 1.5)\), \(\beta \sim \mathcal{N}(0, 1)\), and \(\phi \sim \text{Exponential}(1)\):
# Block-level data used for Bayesian model
block_df <- make_block_level_df(
strikes_df = strikes,
holidays_df = holidays,
recent_start = recent_start
)
# Fit a beta-bionomial model
fit_betabinom <- brm(
formula = n_holiday | trials(n_days) ~ recent,
data = block_df,
family = beta_binomial(),
prior = c(
prior(normal(0, 1.5), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(exponential(1), class = "phi")
),
chains = 4,
iter = 4000,
warmup = 1000,
seed = 123,
refresh = 0
)
print(summary(fit_betabinom)) Family: beta_binomial
Links: mu = logit
Formula: n_holiday | trials(n_days) ~ recent
Data: block_df (Number of observations: 15)
Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
total post-warmup draws = 12000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.10 0.54 -2.20 -0.07 1.00 9892 8567
recent 0.65 0.72 -0.78 2.08 1.00 11412 8532
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 0.57 0.37 0.12 1.54 1.00 10153 7796
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Very nice… but what does this mean?
Code
# Convert posterior draws into intuitive quantities for one fitted model
extract_bayes_quantities <- function(fit, model_label) {
draws <- as_draws_df(fit)
draws %>%
transmute(
model = model_label,
p_earlier = plogis(b_Intercept),
p_recent = plogis(b_Intercept + b_recent),
risk_diff = p_recent - p_earlier,
risk_ratio = p_recent / p_earlier,
odds_ratio = exp(b_recent),
recent_gt_earlier = p_recent > p_earlier
)
}
# Summarise posterior draws into readable quantities
summarise_bayes_human <- function(draws_df) {
draws_df %>%
group_by(model) %>%
summarise(
earlier_pct = 100 * mean(p_earlier),
earlier_pct_lo = 100 * quantile(p_earlier, 0.025),
earlier_pct_hi = 100 * quantile(p_earlier, 0.975),
recent_pct = 100 * mean(p_recent),
recent_pct_lo = 100 * quantile(p_recent, 0.025),
recent_pct_hi = 100 * quantile(p_recent, 0.975),
abs_diff_pct_points = 100 * mean(risk_diff),
abs_diff_lo = 100 * quantile(risk_diff, 0.025),
abs_diff_hi = 100 * quantile(risk_diff, 0.975),
rel_risk = mean(risk_ratio),
rel_risk_lo = quantile(risk_ratio, 0.025),
rel_risk_hi = quantile(risk_ratio, 0.975),
odds_ratio = mean(odds_ratio),
odds_ratio_lo = quantile(odds_ratio, 0.025),
odds_ratio_hi = quantile(odds_ratio, 0.975),
pr_recent_gt_earlier = mean(recent_gt_earlier)
) %>%
ungroup()
}draws_betabinom <- extract_bayes_quantities(fit_betabinom, "Beta-binomial")
bayes_human_summary <- summarise_bayes_human(draws_betabinom)The estimated probability that a strike day within a block fell in a school holiday was 26.1% for earlier strikes (95% credible interval 10 to 48.2%) and 39.8% for the last 18 months (14.2 to 69.4). We can show this graphically:
plot_prob_df <- draws_betabinom %>%
select(p_earlier, p_recent) %>%
pivot_longer(
cols = c(p_earlier, p_recent),
names_to = "period",
values_to = "probability"
) %>%
mutate(
period = recode(period,
p_earlier = "Earlier strikes",
p_recent = "Last 18 months"),
probability_pct = 100 * probability
)
ggplot(plot_prob_df, aes(x = probability_pct, fill = period)) +
geom_density(alpha = 0.4) +
labs(
subtitle = "Estimated probability that a strike day within a block fell in a school holiday",
x = "Probability (%)",
y = "Density",
fill = NULL
) +
theme_bw()
This corresponds to an absolute increase of 13.7% (-15 to 44.1), or a relative risk of a strike falling on a school holiday of 1.71 (0.57 to 3.94). The posterior probability that recent overlap exceeded earlier overlap was 82.1% i.e. there’s an 18% chance that this isn’t a real signal that we’re seeing:
# Build density data and flag whether each x-value is > 0
dens_df <- density(draws_betabinom$risk_diff, n = 2000) %>%
with(tibble(
risk_diff = x,
density = y,
above_zero = x > 0
))
# Plot
ggplot(dens_df, aes(x = 100 * risk_diff, y = density)) +
geom_area(
data = dens_df %>% filter(!above_zero),
alpha = 0.35
) +
geom_area(
data = dens_df %>% filter(above_zero),
alpha = 0.75
) +
geom_vline(xintercept = 0, linetype = "dashed", linewidth = 1) +
annotate(
"text",
x = 100 * mean(draws_betabinom$risk_diff),
y = max(dens_df$density) * 0.9,
label = paste0(
"Pr(recent > earlier) = ",
scales::percent(mean(draws_betabinom$recent_gt_earlier), accuracy = 0.1)
),
vjust = 0
) +
labs(
subtitle = "Shaded area to the right of 0 is the posterior probability",
x = "Recent minus earlier overlap probability (percentage points)",
y = "Posterior density"
) +
theme_bw()
Conclusions
Overall, these analyses suggest that English resident doctor strikes have overlapped school holidays more often than would be expected from a simple day-level comparison, but the strength of that signal weakens once the calendar structure of strike action is modelled properly.
In the more appropriate block-preserving permutation analyses, the full-period overlap only appeared slightly higher than expected by chance, and the recent-versus-earlier comparison suggested that the more recent strike period may have been more holiday-aligned than the earlier one without providing strong proof that this difference was unusual under random block placement.
The Bayesian model pointed in the same direction, estimating higher holiday overlap in recent strike blocks than earlier ones, but with substantial uncertainty once between-block heterogeneity was allowed for.
Taken together, the results are consistent with some evidence of increased holiday alignment, particularly in the more recent period, but not strong enough evidence to support a confident claim of deliberate targeting.
Statistically, these findings speak only to whether the observed timing looks unusual relative to chance-based calendar placement under the chosen models; they do not establish intent, strategy, or political motivation, because real-world strike timing is also shaped by negotiations, notice periods, operational pressures, and other non-random constraints that are not fully captured in the analysis. Having said that, industrial action is intended to be disruptive; targetting school holidays (when consultant activity is already stretched due to childcare and people travelling) wouldn’t be an unreasonable approach to take with that goal in mind.