---
title: "Bulk Billing Rose After November 2025. Did GP Access Improve?"
description: "Billing changed sharply after the November 2025 GP incentives, especially for working-age adults. The evidence for more care is much weaker."
date: "2026-05-28"
date-modified: "2026-10-20"
categories: [medicare, bulk-billing, health-policy, australia]
author: "Aydin"
---
```{r}
#| label: setup
#| include: false
library(readxl)
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(lubridate)
library(scales)
library(knitr)
source(file.path("..", "..", "R", "data_sources.R"))
theme_set(
theme_minimal(base_size = 12) +
theme(
plot.title.position = "plot",
plot.caption.position = "plot",
panel.grid.minor = element_blank()
)
)
policy_date <- as.Date("2025-11-01")
latest_expected_month <- as.Date("2026-06-01")
latest_expected_quarter <- as.Date("2026-06-01")
raw_dir <- "data-raw"
source_registry <- read_source_registry("sources.csv")
source_files <- lapply(
seq_len(nrow(source_registry)),
function(index) download_registered_source(source_registry[index, , drop = FALSE], raw_dir)
)
source_manifest <- write_source_manifest(source_files, file.path(raw_dir, "source-manifest.csv"))
source_path <- function(source_id) {
index <- match(source_id, source_registry$source_id)
if (is.na(index)) stop("Missing source in sources.csv: ", source_id, call. = FALSE)
source_files[[index]]$path
}
clean_names <- function(x) {
x |>
str_to_lower() |>
str_replace_all("&", "and") |>
str_replace_all("[^a-z0-9]+", "_") |>
str_replace_all("^_|_$", "")
}
assert_true <- function(condition, message) {
if (!isTRUE(condition)) stop(message, call. = FALSE)
}
quarter_to_date <- function(x) {
financial_year_start <- as.integer(str_sub(x, 1, 4))
quarter <- as.integer(str_match(x, "Q([1-4])")[, 2])
month <- c(9L, 12L, 3L, 6L)[quarter]
year <- if_else(quarter <= 2L, financial_year_start, financial_year_start + 1L)
as.Date(sprintf("%04d-%02d-01", year, month))
}
parse_aihw_table <- function(path, sheet, id_names) {
raw <- read_excel(path, sheet = sheet, col_names = FALSE, .name_repair = "minimal")
assert_true(ncol(raw) > 500, "The AIHW workbook has fewer month columns than expected.")
assert_true(
identical(as.character(unlist(raw[2, seq_along(id_names)])), id_names),
paste0("The AIHW ", sheet, " header has changed; review the parser.")
)
first_month_column <- length(id_names) + 1L
month_dates <- as.Date(
suppressWarnings(as.numeric(as.character(unlist(raw[2, first_month_column:ncol(raw)])))),
origin = "1899-12-30"
)
assert_true(!anyNA(month_dates) && !anyDuplicated(month_dates), "AIHW month headers are invalid or duplicated.")
body <- raw[-c(1, 2), ]
names(body) <- c(clean_names(id_names), as.character(month_dates))
body |>
mutate(across(first_month_column:ncol(body), as.character)) |>
pivot_longer(first_month_column:ncol(body), names_to = "month", values_to = "value") |>
mutate(month = as.Date(month), value = suppressWarnings(as.numeric(value)))
}
parse_aihw_table_1 <- function(path) {
parse_aihw_table(path, "Table 1", c("Demography", "Demography value", "Metric")) |>
rename(group = demography_value)
}
parse_aihw_table_2 <- function(path) {
parse_aihw_table(path, "Table 2", c("State", "Demography", "Demography value", "Metric")) |>
rename(group = demography_value)
}
parse_sa3_quarterly <- function(path) {
data <- read_excel(path, sheet = "GP NRA by SA3", skip = 4)
names(data) <- make.unique(clean_names(names(data)))
required <- c(
"quarter", "state", "sa3", "broad_type_of_service", "services",
"bulk_billed_services", "patient_billed_services", "mbs_bulk_billing_rate",
"patient_billed_benefits", "patient_billed_fee_charged"
)
assert_true(all(required %in% names(data)), "The Department SA3 workbook columns have changed.")
data |>
transmute(
quarter,
quarter_date = quarter_to_date(quarter),
state,
sa3,
sa3_id = paste(state, sa3, sep = " / "),
service_type = broad_type_of_service,
services = as.numeric(services),
bulk_billed_services = as.numeric(bulk_billed_services),
patient_billed_services = as.numeric(patient_billed_services),
bulk_billing_rate = as.numeric(mbs_bulk_billing_rate),
patient_billed_benefits = as.numeric(patient_billed_benefits),
patient_billed_fee_charged = as.numeric(patient_billed_fee_charged),
patient_contributions = patient_billed_fee_charged - patient_billed_benefits
)
}
# Bartlett-kernel HAC covariance for the short monthly interrupted time series.
vcov_hac <- function(model, lag = NULL) {
design <- model.matrix(model)
errors <- residuals(model)
observations <- nrow(design)
parameters <- ncol(design)
if (is.null(lag)) lag <- floor(4 * (observations / 100)^(2 / 9))
meat <- crossprod(design * errors)
for (offset in seq_len(lag)) {
covariance <- crossprod(
design[(offset + 1):observations, , drop = FALSE] * errors[(offset + 1):observations],
design[1:(observations - offset), , drop = FALSE] * errors[1:(observations - offset)]
)
meat <- meat + (1 - offset / (lag + 1)) * (covariance + t(covariance))
}
bread <- solve(crossprod(design))
bread %*% meat %*% bread * observations / (observations - parameters)
}
# One-way cluster-robust covariance, with a finite-sample correction.
vcov_cluster <- function(model, cluster) {
design <- model.matrix(model)
errors <- residuals(model)
model_weights <- weights(model)
if (is.null(model_weights)) model_weights <- rep(1, length(errors))
bread <- solve(crossprod(design * sqrt(model_weights)))
scores <- design * (model_weights * errors)
groups <- split(seq_along(cluster), cluster)
meat <- Reduce(`+`, lapply(groups, function(index) {
score_sum <- colSums(scores[index, , drop = FALSE])
tcrossprod(score_sum)
}))
observations <- nrow(design)
parameters <- ncol(design)
group_count <- length(groups)
correction <- (group_count / (group_count - 1)) * ((observations - 1) / (observations - parameters))
bread %*% meat %*% bread * correction
}
aihw_path <- source_path("aihw_monthly_bulk_billing")
aihw <- parse_aihw_table_1(aihw_path)
aihw_state <- parse_aihw_table_2(aihw_path)
sa3_raw <- parse_sa3_quarterly(source_path("health_sa3_quarterly"))
national_monthly <- aihw |>
filter(
demography == "Total",
group == "Total",
metric == "Bulk-billing rate",
month >= as.Date("2023-01-01")
) |>
arrange(month) |>
rename(bulk_billing_rate = value)
assert_true(max(national_monthly$month) == latest_expected_month, "The latest AIHW month is not June 2026.")
assert_true(nrow(national_monthly) == 42, "The recent national AIHW series is not a complete 42-month panel.")
assert_true(!anyDuplicated(national_monthly$month), "The national AIHW series has duplicate months.")
assert_true(
identical(national_monthly$month, seq(as.Date("2023-01-01"), latest_expected_month, by = "month")),
"The national AIHW series has missing months."
)
assert_true(all(is.finite(national_monthly$bulk_billing_rate)), "The national AIHW series contains missing or non-finite rates.")
assert_true(all(between(national_monthly$bulk_billing_rate, 0, 1)), "AIHW rates must lie between zero and one.")
age_rates <- aihw |>
filter(
demography == "Age group", metric == "Bulk-billing rate",
month %in% as.Date(c("2025-06-01", "2026-06-01"))
) |>
select(group, month, bulk_billing_rate = value)
seifa_rates <- aihw |>
filter(
demography == "SEIFA", metric == "Bulk-billing rate",
month %in% as.Date(c("2025-06-01", "2026-06-01"))
) |>
select(group, month, bulk_billing_rate = value)
assert_true(n_distinct(age_rates$group) == 3 && nrow(age_rates) == 6, "AIHW age data are incomplete.")
assert_true(n_distinct(seifa_rates$group) == 5 && nrow(seifa_rates) == 10, "AIHW SEIFA data are incomplete.")
assert_true(
all(is.finite(age_rates$bulk_billing_rate)) && all(between(age_rates$bulk_billing_rate, 0, 1)),
"AIHW age rates contain missing, non-finite or out-of-range values."
)
assert_true(
all(is.finite(seifa_rates$bulk_billing_rate)) && all(between(seifa_rates$bulk_billing_rate, 0, 1)),
"AIHW SEIFA rates contain missing, non-finite or out-of-range values."
)
age_change <- age_rates |>
mutate(period = paste0("rate_", str_replace_all(as.character(month), "-", "_"))) |>
select(-month) |>
pivot_wider(names_from = period, values_from = bulk_billing_rate) |>
transmute(group, change_pp = 100 * (rate_2026_06_01 - rate_2025_06_01))
seifa_change <- seifa_rates |>
mutate(period = paste0("rate_", str_replace_all(as.character(month), "-", "_"))) |>
select(-month) |>
pivot_wider(names_from = period, values_from = bulk_billing_rate) |>
transmute(
group = case_when(
str_detect(group, "quintile 1") ~ "Quintile 1 (most disadvantaged)",
str_detect(group, "quintile 5") ~ "Quintile 5 (most advantaged)",
TRUE ~ str_to_sentence(group)
),
june_2025 = rate_2025_06_01,
june_2026 = rate_2026_06_01,
change_pp = 100 * (rate_2026_06_01 - rate_2025_06_01)
)
age_rate_summary <- age_rates |>
mutate(period = paste0("rate_", str_replace_all(as.character(month), "-", "_"))) |>
select(-month) |>
pivot_wider(names_from = period, values_from = bulk_billing_rate) |>
mutate(change_pp = 100 * (rate_2026_06_01 - rate_2025_06_01))
working_age_monthly <- aihw |>
filter(
demography == "Age group", group == "16-64",
metric %in% c("Bulk-billing rate", "Bulk-billed services"),
month >= as.Date("2023-01-01")
) |>
transmute(month, metric = clean_names(metric), value) |>
pivot_wider(names_from = metric, values_from = value) |>
mutate(visits = bulk_billed_services / bulk_billing_rate) |>
arrange(month)
assert_true(nrow(working_age_monthly) == 42, "The working-age monthly series is incomplete.")
assert_true(
all(is.finite(working_age_monthly$visits)) && all(working_age_monthly$visits > 0),
"Working-age visit counts could not be reconstructed from the published numerator and rate."
)
working_age_rates <- age_rate_summary |>
filter(group == "16-64")
state_working_age_rates <- aihw_state |>
filter(
demography == "Age group", group == "16-64", metric == "Bulk-billing rate",
month %in% as.Date(c("2025-06-01", "2026-06-01"))
) |>
select(state, month, bulk_billing_rate = value) |>
mutate(period = paste0("rate_", str_replace_all(as.character(month), "-", "_"))) |>
select(-month) |>
pivot_wider(names_from = period, values_from = bulk_billing_rate) |>
mutate(change_pp = 100 * (rate_2026_06_01 - rate_2025_06_01))
assert_true(
nrow(state_working_age_rates) == 8 && !anyNA(state_working_age_rates),
"The state-by-working-age rate comparison is incomplete."
)
national_model_data <- national_monthly |>
mutate(
time = row_number(),
post = as.integer(month >= policy_date),
post_time = if_else(post == 1L, interval(policy_date, month) %/% months(1), 0),
calendar_month = factor(month(month))
)
fit_national_model <- function(data) {
model <- lm(bulk_billing_rate ~ time + post + post_time + calendar_month, data = data)
covariance <- vcov_hac(model)
tibble(
term = names(coef(model)),
estimate = unname(coef(model)),
std_error = sqrt(diag(covariance))
) |>
mutate(lower = estimate - 1.96 * std_error, upper = estimate + 1.96 * std_error)
}
national_model <- fit_national_model(national_model_data)
step_estimate <- filter(national_model, term == "post")
latest_rate <- filter(national_monthly, month == latest_expected_month) |> pull(bulk_billing_rate)
june_2025_rate <- filter(national_monthly, month == as.Date("2025-06-01")) |> pull(bulk_billing_rate)
working_age_model_data <- working_age_monthly |>
mutate(
time = row_number(),
post = as.integer(month >= policy_date),
post_time = if_else(post == 1L, interval(policy_date, month) %/% months(1), 0),
calendar_month = factor(month(month))
)
fit_working_age_visit_model <- function(formula, specification) {
model <- lm(formula, data = working_age_model_data)
covariance <- vcov_hac(model)
estimate <- coef(model)[["post"]]
std_error <- sqrt(covariance["post", "post"])
tibble(
specification,
estimate,
std_error,
change_percent = 100 * (exp(estimate) - 1),
lower_percent = 100 * (exp(estimate - 1.96 * std_error) - 1),
upper_percent = 100 * (exp(estimate + 1.96 * std_error) - 1)
)
}
working_age_visit_models <- bind_rows(
fit_working_age_visit_model(
log(visits) ~ time + post + post_time + calendar_month,
"Step and post-policy slope (primary)"
),
fit_working_age_visit_model(
log(visits) ~ time + post + calendar_month,
"Step only (sensitivity)"
)
)
working_age_visit_step <- slice(working_age_visit_models, 1)
working_age_visit_comparison <- working_age_monthly |>
filter(
month %in% as.Date(c(
"2025-01-01", "2025-02-01",
"2025-03-01", "2025-04-01", "2025-05-01", "2025-06-01",
"2026-01-01", "2026-02-01",
"2026-03-01", "2026-04-01", "2026-05-01", "2026-06-01"
))
) |>
mutate(calendar_month = month(month, label = TRUE, abbr = TRUE), year = year(month)) |>
select(calendar_month, year, visits) |>
pivot_wider(names_from = year, values_from = visits, names_prefix = "visits_") |>
mutate(change = visits_2026 / visits_2025 - 1)
working_age_visit_change <- sum(working_age_visit_comparison$visits_2026) /
sum(working_age_visit_comparison$visits_2025) - 1
sa3_numeric <- sa3_raw |>
select(
services, bulk_billed_services, patient_billed_services, bulk_billing_rate,
patient_billed_benefits, patient_billed_fee_charged, patient_contributions
)
assert_true(!anyNA(sa3_raw$quarter_date), "Department quarter labels contain missing or invalid dates.")
assert_true(
!anyNA(sa3_numeric) && all(vapply(sa3_numeric, function(value) all(is.finite(value)), logical(1))),
"Department service fields contain missing or non-finite values."
)
assert_true(max(sa3_raw$quarter_date) == latest_expected_quarter, "The latest Department quarter is not June 2026.")
assert_true(all(sa3_raw$services >= 0), "Department service counts must be non-negative.")
assert_true(all(between(sa3_raw$bulk_billing_rate, 0, 1)), "Department rates must lie between zero and one.")
assert_true(
all(abs(sa3_raw$bulk_billed_services + sa3_raw$patient_billed_services - sa3_raw$services) < 0.01),
"Department billed-service components do not add to total services."
)
assert_true(
all(abs(sa3_raw$bulk_billed_services / sa3_raw$services - sa3_raw$bulk_billing_rate) < 1e-10),
"Department published rates do not reconcile to their numerator and denominator."
)
national_quarterly <- sa3_raw |>
filter(state == "Australia", sa3 == "Australia") |>
arrange(quarter_date)
sa3_panel <- sa3_raw |>
filter(
str_to_upper(state) %in% c("NSW", "VIC", "QLD", "SA", "WA", "TAS", "NT", "ACT"),
str_to_upper(sa3) != str_to_upper(state),
quarter_date >= as.Date("2023-09-01"), quarter_date <= latest_expected_quarter,
services > 0
)
quarter_counts <- sa3_panel |> count(quarter_date)
balanced_ids <- sa3_panel |> count(sa3_id) |> filter(n == 12) |> pull(sa3_id)
assert_true(nrow(quarter_counts) == 12 && n_distinct(quarter_counts$n) == 1, "The SA3 panel is not balanced by quarter.")
assert_true(length(balanced_ids) == 332, "The expected 332 classified SA3 areas are not present in all 12 quarters.")
sa3_reconciliation <- sa3_panel |>
group_by(quarter_date) |>
summarise(sa3_services = sum(services), .groups = "drop") |>
left_join(national_quarterly |> select(quarter_date, national_services = services), by = "quarter_date") |>
mutate(coverage = sa3_services / national_services)
assert_true(all(sa3_reconciliation$coverage > 0.997), "Classified SA3 services do not closely reconcile to the national total.")
baseline <- sa3_panel |>
filter(
sa3_id %in% balanced_ids,
quarter_date >= as.Date("2024-12-01"), quarter_date <= as.Date("2025-09-01")
) |>
group_by(sa3_id) |>
summarise(
baseline_non_bulk_billed_share = 1 - sum(bulk_billed_services) / sum(services),
baseline_services = mean(services), baseline_quarters = n(), .groups = "drop"
)
assert_true(all(baseline$baseline_quarters == 4), "Each SA3 must have four baseline quarters.")
sa3_model_data <- sa3_panel |>
filter(sa3_id %in% balanced_ids, quarter_date != as.Date("2025-12-01")) |>
inner_join(baseline, by = "sa3_id") |>
mutate(
post = as.integer(quarter_date >= as.Date("2026-03-01")),
exposure_10pp = baseline_non_bulk_billed_share / 0.10,
post_exposure = post * exposure_10pp,
log_services = log(services),
model_weight = baseline_services / mean(baseline_services)
)
assert_true(
nrow(sa3_model_data) == 332 * 11 && sum(sa3_model_data$post) == 332 * 2,
"The SA3 model must contain 11 quarters, including two clean post-policy quarters."
)
fit_sa3_model <- function(data, weighted = FALSE) {
model <- if (weighted) {
lm(log_services ~ post_exposure + factor(sa3_id) + factor(quarter_date), data = data, weights = model_weight)
} else {
lm(log_services ~ post_exposure + factor(sa3_id) + factor(quarter_date), data = data)
}
covariance <- vcov_cluster(model, data$sa3_id)
estimate <- coef(model)[["post_exposure"]]
std_error <- sqrt(covariance["post_exposure", "post_exposure"])
tibble(
specification = if_else(weighted, "Baseline-service weighted", "Unweighted (primary)"),
estimate, std_error,
change_percent = 100 * (exp(estimate) - 1),
lower_percent = 100 * (exp(estimate - 1.96 * std_error) - 1),
upper_percent = 100 * (exp(estimate + 1.96 * std_error) - 1)
)
}
sa3_models <- bind_rows(
fit_sa3_model(sa3_model_data, weighted = FALSE),
fit_sa3_model(sa3_model_data, weighted = TRUE)
)
patient_cost_national <- national_quarterly |>
filter(quarter_date %in% as.Date(c("2025-03-01", "2025-06-01", "2026-03-01", "2026-06-01"))) |>
mutate(period = if_else(year(quarter_date) == 2025, "Before: Mar + Jun 2025", "After: Mar + Jun 2026")) |>
group_by(period) |>
summarise(
services = sum(services),
bulk_billed_services = sum(bulk_billed_services),
patient_billed_services = sum(patient_billed_services),
patient_contributions = sum(patient_contributions),
average_per_service = patient_contributions / services,
average_per_patient_billed_service = patient_contributions / patient_billed_services,
.groups = "drop"
)
patient_cost_change <- patient_cost_national |>
summarise(
bulk_billed_service_change = bulk_billed_services[period == "After: Mar + Jun 2026"] -
bulk_billed_services[period == "Before: Mar + Jun 2025"],
contribution_change = patient_contributions[period == "After: Mar + Jun 2026"] -
patient_contributions[period == "Before: Mar + Jun 2025"],
average_change = average_per_service[period == "After: Mar + Jun 2026"] -
average_per_service[period == "Before: Mar + Jun 2025"],
patient_billed_average_change = average_per_patient_billed_service[period == "After: Mar + Jun 2026"] -
average_per_patient_billed_service[period == "Before: Mar + Jun 2025"]
)
cost_area_groups <- baseline |>
transmute(
sa3_id,
quintile = ntile(1 - baseline_non_bulk_billed_share, 5),
group = factor(
quintile,
levels = 1:5,
labels = c(
"Q1 - lowest bulk billing", "Q2", "Q3", "Q4", "Q5 - highest bulk billing"
)
)
)
patient_cost_by_area_group <- sa3_panel |>
filter(
sa3_id %in% balanced_ids,
quarter_date %in% as.Date(c("2025-03-01", "2025-06-01", "2026-03-01", "2026-06-01"))
) |>
inner_join(cost_area_groups, by = "sa3_id") |>
mutate(period = if_else(year(quarter_date) == 2025, "Before", "After")) |>
group_by(group, period) |>
summarise(
services = sum(services),
patient_billed_services = sum(patient_billed_services),
patient_contributions = sum(patient_contributions),
average_per_service = patient_contributions / services,
average_per_patient_billed_service = patient_contributions / patient_billed_services,
.groups = "drop"
) |>
pivot_wider(
names_from = period,
values_from = c(
services, patient_billed_services, patient_contributions,
average_per_service, average_per_patient_billed_service
)
) |>
mutate(
average_change = average_per_service_After - average_per_service_Before,
patient_billed_average_change =
average_per_patient_billed_service_After - average_per_patient_billed_service_Before
)
state_patient_costs <- sa3_raw |>
filter(
str_to_upper(state) %in% c("NSW", "VIC", "QLD", "SA", "WA", "TAS", "NT", "ACT"),
str_to_upper(sa3) == str_to_upper(state),
quarter_date %in% as.Date(c("2025-03-01", "2025-06-01", "2026-03-01", "2026-06-01"))
) |>
mutate(period = if_else(year(quarter_date) == 2025, "Before", "After")) |>
group_by(state, period) |>
summarise(
services = sum(services),
patient_contributions = sum(patient_contributions),
average_per_service = patient_contributions / services,
.groups = "drop"
) |>
pivot_wider(names_from = period, values_from = c(services, patient_contributions, average_per_service)) |>
mutate(
contribution_change = patient_contributions_After - patient_contributions_Before,
average_change = average_per_service_After - average_per_service_Before
)
```
::: {.callout-note title="Executive summary"}
- Australia's GP bulk-billing rate rose from `r percent(june_2025_rate, accuracy = 0.1)` to `r percent(latest_rate, accuracy = 0.1)` in the year to June 2026.
- For 16–64-year-olds it rose from `r percent(working_age_rates$rate_2025_06_01, accuracy = 0.1)` to `r percent(working_age_rates$rate_2026_06_01, accuracy = 0.1)`, while reconstructed visit volumes were `r percent(working_age_visit_change, accuracy = 0.1)` higher across January–June.
- Across all ages, `r number(patient_cost_change$bulk_billed_service_change / 1e6, accuracy = 0.1)` million more GP services were bulk billed and recorded patient contributions were `r dollar(abs(patient_cost_change$contribution_change) / 1e6, suffix = "m", accuracy = 0.1)` lower across the first six months of 2026 than a year earlier.
:::
The policy period produced a clear change in **how GP care was billed**. The evidence for a matching change in **how much care people received** is much weaker.
The comparisons below start in 2023, after the acute COVID-era disruption. January–June 2026 provides the consistent six-month post-policy window used for the year-on-year access and cost comparisons; December 2025 is excluded from the quarterly area model as a transition quarter. [Methods, caveats, detailed tables and the source audit are in the expandable appendix.](#appendix)
## 1. National bulk-billing rates rose from `r percent(june_2025_rate, accuracy = 0.1)` to `r percent(latest_rate, accuracy = 0.1)`
The national service-level rate rose `r number(100 * (latest_rate - june_2025_rate), accuracy = 0.1)` percentage points in the year to June 2026. The break is visible immediately after the 1 November 2025 incentive changes.
```{r}
#| label: fig-national-rate
#| fig-cap: "Source: AIHW monthly GP bulk-billing data. The final three months are preliminary."
#| fig-alt: "Line chart of Australia's monthly GP bulk-billing rate from January 2023 to June 2026. The rate rises sharply after the November 2025 policy line, from 78.8 per cent in June 2025 to 83.6 per cent in June 2026. April to June 2026 are shaded as preliminary."
#| fig-width: 9
#| fig-height: 5.2
ggplot(national_monthly, aes(month, bulk_billing_rate)) +
annotate("rect", xmin = as.Date("2026-04-01"), xmax = as.Date("2026-07-01"), ymin = -Inf, ymax = Inf, fill = "grey80", alpha = 0.35) +
geom_vline(xintercept = policy_date, colour = "#B8472D", linewidth = 0.8) +
geom_line(colour = "#245A73", linewidth = 1) +
geom_point(data = filter(national_monthly, month == latest_expected_month), colour = "#245A73", size = 2.8) +
annotate("text", x = policy_date, y = 0.765, label = "1 Nov 2025\nchanges begin", hjust = 1.05, vjust = 0, colour = "#8F3522", size = 3.4) +
annotate("text", x = latest_expected_month, y = latest_rate, label = percent(latest_rate, accuracy = 0.1), hjust = -0.15, vjust = -0.5, colour = "#245A73", fontface = "bold") +
scale_x_date(date_breaks = "6 months", date_labels = "%b\n%Y", expand = expansion(mult = c(0.01, 0.08))) +
scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0.74, 0.86)) +
labs(
title = paste0(
"National bulk-billing rate rose from ", percent(june_2025_rate, accuracy = 0.1),
" to ", percent(latest_rate, accuracy = 0.1)
),
subtitle = "Monthly share of claimed GP attendances bulk billed, Australia",
x = NULL, y = NULL
)
```
A descriptive interrupted time-series model, controlling for the 2023–26 trend and calendar month, estimates an immediate `r number(100 * step_estimate$estimate, accuracy = 0.1)`-percentage-point step (HAC 95% interval `r number(100 * step_estimate$lower, accuracy = 0.1)` to `r number(100 * step_estimate$upper, accuracy = 0.1)`). There are only eight post-policy monthly observations. The estimate marks a break in timing; it is not proof that the policy was the only cause.
## 2. Bulk billing for ages 16–64 rose from `r percent(working_age_rates$rate_2025_06_01, accuracy = 0.1)` to `r percent(working_age_rates$rate_2026_06_01, accuracy = 0.1)`
The age pattern is the cleanest sign that the expansion reached its target group. The rate for 16–64-year-olds rose `r number(working_age_rates$change_pp, accuracy = 0.1)` percentage points; the changes for children and people aged 65 and over were much smaller.
```{r}
#| label: fig-age-change
#| fig-cap: "Source: AIHW monthly GP bulk-billing data. June 2026 is preliminary. Changes are percentage points."
#| fig-alt: "Horizontal dot plot of the June 2025 to June 2026 change in GP bulk-billing rates by age. Ages 16 to 64 increased by 8.0 percentage points, compared with 1.0 point for ages 0 to 15 and 1.4 points for ages 65 and over."
#| fig-width: 8
#| fig-height: 4.2
age_change |>
mutate(group = factor(group, levels = c("65+", "0-15", "16-64")), emphasis = group == "16-64") |>
ggplot(aes(change_pp, group)) +
geom_segment(aes(x = 0, xend = change_pp, yend = group), colour = "grey75", linewidth = 1.2) +
geom_point(aes(colour = emphasis), size = 4, show.legend = FALSE) +
geom_text(aes(label = paste0(number(change_pp, accuracy = 0.1), " pp")), hjust = -0.2, fontface = "bold", size = 3.8) +
scale_colour_manual(values = c(`FALSE` = "#77838D", `TRUE` = "#B8472D")) +
scale_x_continuous(limits = c(0, 9), breaks = seq(0, 8, 2), expand = expansion(mult = c(0, 0.03))) +
labs(
title = paste0(
"Bulk billing for ages 16-64 rose from ",
percent(working_age_rates$rate_2025_06_01, accuracy = 0.1), " to ",
percent(working_age_rates$rate_2026_06_01, accuracy = 0.1)
),
subtitle = "Change in GP bulk-billing rate, June 2025 to June 2026",
x = "Percentage-point change", y = NULL
)
```
From here, the age-specific analysis focuses on 16–64-year-olds. Children and older Australians were already covered by earlier bulk-billing incentives. The 16–64 band is still broader than the newly eligible population because it includes concession-card holders who were already eligible.
## 3. The rise was nationwide and smallest in the most advantaged areas
Among 16–64-year-olds, the bulk-billing rate rose in every state and territory. The increase ranged from `r number(min(state_working_age_rates$change_pp), accuracy = 0.1)` percentage points in the ACT to `r number(max(state_working_age_rates$change_pp), accuracy = 0.1)` points in the Northern Territory.
```{r}
#| label: fig-state-working-age
#| fig-cap: "Source: AIHW monthly GP bulk-billing data. Ages 16–64; June 2026 is preliminary."
#| fig-alt: "Dumbbell chart comparing June 2025 and June 2026 working-age bulk-billing rates in each state and territory. Every jurisdiction increased, with the largest increase in the Northern Territory and the smallest in the Australian Capital Territory."
#| fig-width: 9
#| fig-height: 5.6
state_working_age_rates |>
mutate(state = reorder(state, change_pp)) |>
ggplot(aes(y = state)) +
geom_segment(aes(x = rate_2025_06_01, xend = rate_2026_06_01, yend = state), colour = "grey72", linewidth = 1.5) +
geom_point(aes(x = rate_2025_06_01), colour = "#77838D", size = 3) +
geom_point(aes(x = rate_2026_06_01), colour = "#B8472D", size = 3.5) +
geom_text(aes(x = rate_2026_06_01, label = paste0("+", number(change_pp, accuracy = 0.1), " pp")), hjust = -0.15, size = 3.2) +
scale_x_continuous(labels = percent_format(accuracy = 1), limits = c(0.4, 0.93), expand = expansion(mult = c(0, 0.04))) +
labs(
title = "Working-age bulk billing rose in every state and territory",
subtitle = "June 2025 (grey) to June 2026 (orange)",
x = NULL, y = NULL
)
```
All-age rates were highest in the most disadvantaged areas before the change and rose across every SEIFA quintile. The changes were not monotonic by disadvantage: quintiles 2 and 3 rose most, while the most advantaged quintile recorded the smallest increase, at `r number(filter(seifa_change, str_detect(group, "most advantaged"))$change_pp, accuracy = 0.1)` points.
```{r}
#| label: fig-seifa-rates
#| fig-cap: "Source: AIHW monthly GP bulk-billing data. These SEIFA aggregates cover all ages because the published workbook does not cross age with SEIFA. June 2026 is preliminary."
#| fig-alt: "Dumbbell chart comparing all-age bulk-billing rates across five socioeconomic-area quintiles in June 2025 and June 2026. Rates rose in every quintile, with the most disadvantaged areas remaining highest and the most advantaged areas recording the smallest increase."
#| fig-width: 9
#| fig-height: 5.4
seifa_change |>
mutate(group = factor(group, levels = rev(group))) |>
ggplot(aes(y = group)) +
geom_segment(aes(x = june_2025, xend = june_2026, yend = group), colour = "grey72", linewidth = 1.5) +
geom_point(aes(x = june_2025), colour = "#77838D", size = 3) +
geom_point(aes(x = june_2026), colour = "#245A73", size = 3.5) +
geom_text(aes(x = june_2026, label = paste0("+", number(change_pp, accuracy = 0.1), " pp")), hjust = -0.15, size = 3.2) +
scale_x_continuous(labels = percent_format(accuracy = 1), limits = c(0.65, 0.96), expand = expansion(mult = c(0, 0.05))) +
labs(
title = "Bulk billing rose across every socioeconomic-area quintile",
subtitle = "All ages; June 2025 (grey) to June 2026 (blue)",
x = NULL, y = NULL
)
```
SEIFA describes the area where a service was delivered, not the circumstances of each patient. It is useful distributional context, not a causal test.
## 4. Working-age visits were `r percent(working_age_visit_change, accuracy = 0.1)` higher across January–June
The working-age billing shift did not translate into a clear, sustained rise in claimed visits. Across January–June 2026, reconstructed 16–64 visit volumes were `r percent(working_age_visit_change, accuracy = 0.1)` above the same months in 2025. The monthly pattern moved between increases and declines.
```{r}
#| label: fig-working-age-visits
#| fig-cap: "Source: AIHW monthly GP bulk-billing data. Total visits for ages 16–64 are reconstructed as bulk-billed services divided by the published bulk-billing rate. The final three months are preliminary."
#| fig-alt: "Bar chart of year-on-year changes in reconstructed working-age GP visits for January to June 2026. The six months together were approximately flat."
#| fig-width: 8
#| fig-height: 4.8
working_age_visit_comparison |>
ggplot(aes(calendar_month, change)) +
geom_hline(yintercept = 0, colour = "grey55", linewidth = 0.7) +
geom_col(aes(fill = change >= 0), width = 0.62, show.legend = FALSE) +
geom_text(
aes(y = if_else(change >= 0, change + 0.004, change - 0.004), label = percent(change, accuracy = 0.1)),
vjust = if_else(working_age_visit_comparison$change >= 0, 0, 1), fontface = "bold", size = 3.6
) +
scale_fill_manual(values = c(`FALSE` = "#77838D", `TRUE` = "#245A73")) +
scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(-0.065, 0.07)) +
labs(
title = paste0("Working-age visits were ", percent(working_age_visit_change, accuracy = 0.1), " higher across January-June"),
subtitle = "Change from the same month in 2025, ages 16-64",
x = NULL, y = NULL
)
```
The post-COVID time-series check is similarly inconclusive. The primary specification estimates an immediate `r percent(working_age_visit_step$change_percent / 100, accuracy = 0.1)` shift in working-age visits, with a 95% interval from `r percent(working_age_visit_step$lower_percent / 100, accuracy = 0.1)` to `r percent(working_age_visit_step$upper_percent / 100, accuracy = 0.1)`. A simpler step-only specification reaches the same broad conclusion, and the area-level check finds no clear relative visit gain in places with more scope to switch billing methods. With only eight post-policy months, these models are supporting checks rather than the headline evidence.
This does not mean access did not improve. Claims cannot show whether waiting times shortened, unmet need fell or continuity improved. It means the available volume evidence does not yet show a decisive access response.
## 5. Recorded patient contributions were lower—but the pattern was uneven
::: {.callout-important title="3.5 million more GP services were bulk billed"}
Across January–June 2026, `r number(patient_cost_change$bulk_billed_service_change / 1e6, accuracy = 0.1)` million more GP services were bulk billed than in the same six months of 2025. Recorded patient contributions were `r dollar(abs(patient_cost_change$contribution_change) / 1e6, suffix = "m", accuracy = 0.1)` lower.
These are observed all-age changes, not estimates of the policy's causal effect. The total number and mix of claimed services also changed.
:::
Across the March and June quarters, recorded patient contributions were `r dollar(abs(patient_cost_change$contribution_change) / 1e6, suffix = "m", accuracy = 0.1)` lower than in the same quarters in 2025. Average recorded contributions across **all** claimed GP services fell from `r dollar(filter(patient_cost_national, str_starts(period, "Before"))$average_per_service, accuracy = 0.01)` to `r dollar(filter(patient_cost_national, str_starts(period, "After"))$average_per_service, accuracy = 0.01)`.
The pattern was not evenly distributed. When SA3s are grouped by their pre-policy bulk-billing rate, the average fell in four of five groups and rose slightly in the group that began with the lowest bulk-billing rates and highest patient costs.
```{r}
#| label: fig-patient-costs-area
#| fig-cap: "Source: Department of Health, Disability and Ageing quarterly Medicare statistics. All ages; classified SA3s; combined March and June quarters. Average patient contribution is total patient-billed fees less Medicare benefits, divided by all claimed services."
#| fig-alt: "Dumbbell chart comparing average patient contributions per GP service before and after the policy across five groups of SA3 areas. Costs fell in four groups but rose slightly in the group with the lowest baseline bulk-billing rate."
#| fig-width: 9
#| fig-height: 5.4
patient_cost_by_area_group |>
ggplot(aes(y = group)) +
geom_segment(aes(x = average_per_service_Before, xend = average_per_service_After, yend = group), colour = "grey72", linewidth = 1.5) +
geom_point(aes(x = average_per_service_Before), colour = "#77838D", size = 3) +
geom_point(aes(x = average_per_service_After), colour = "#B8472D", size = 3.5) +
geom_text(
aes(
x = pmax(average_per_service_Before, average_per_service_After),
label = paste0(if_else(average_change > 0, "+$", "-$"), number(abs(average_change), accuracy = 0.01))
),
hjust = -0.15, size = 3.2
) +
scale_x_continuous(labels = dollar_format(accuracy = 1), limits = c(0, 25), expand = expansion(mult = c(0, 0.04))) +
labs(
title = "Average recorded contributions fell in four of five area groups",
subtitle = "Average contribution per GP service; before (grey) and after (orange)",
x = NULL, y = NULL
)
```
There is another important split: among the visits that were still patient billed, the average contribution rose from `r dollar(filter(patient_cost_national, str_starts(period, "Before"))$average_per_patient_billed_service, accuracy = 0.01)` to `r dollar(filter(patient_cost_national, str_starts(period, "After"))$average_per_patient_billed_service, accuracy = 0.01)`. In other words, many services moved out of patient billing, lowering the average across all care; patients who still paid faced a higher average bill.
## Three-line takeaway
- **Billing:** the policy period coincided with a large, nationwide rise in bulk billing, concentrated among 16–64-year-olds.
- **Access:** working-age visit volumes were only modestly higher and inconsistent month to month; the models do not show a decisive increase.
- **Cost:** `r number(patient_cost_change$bulk_billed_service_change / 1e6, accuracy = 0.1)` million more services were bulk billed and recorded patient contributions were `r dollar(abs(patient_cost_change$contribution_change) / 1e6, suffix = "m", accuracy = 0.1)` lower, but remaining patient-billed visits became more expensive on average.
## Appendix {#appendix}
::: {.callout-warning collapse="true" title="Caveats and interpretation"}
- Bulk billing is measured per service, not per patient. It does not tell us how many people had every visit bulk billed.
- The 16–64 band includes some concession-card holders who were already eligible before November 2025. It is a useful target-age proxy, not a treatment group.
- The published AIHW workbook crosses age with state and LGA, but not with SEIFA. The SEIFA results are therefore all ages.
- The Department's SA3 patient-cost data are also all ages. Patient contributions are patient-billed fees less Medicare benefits; they do not include every cost of obtaining care.
- AIHW marks April–June 2026 as preliminary. Results may be revised.
- Claims measure delivered services. They cannot identify waiting time, unmet need, clinical need, provider supply or continuity of care.
- The policy began on 1 November 2025. December 2025 mixes pre- and post-policy months, so the quarterly model excludes it.
- Working-age visits are approximate. The public monthly workbook does not expose the age-specific denominator directly, so visits are reconstructed as bulk-billed services divided by a published rate rounded to three decimal places.
- These are descriptive and associational comparisons. Other events at the same time may contribute to the observed changes.
:::
::: {.callout-note collapse="true" title="Detailed rates and changes"}
```{r}
#| label: tbl-age-rates
age_rate_summary |>
transmute(
`Age group` = group,
`June 2025` = percent(rate_2025_06_01, accuracy = 0.1),
`June 2026` = percent(rate_2026_06_01, accuracy = 0.1),
`Change (pp)` = number(change_pp, accuracy = 0.1)
) |>
kable(align = c("l", "r", "r", "r"))
```
```{r}
#| label: tbl-state-working-age
state_working_age_rates |>
arrange(desc(change_pp)) |>
transmute(
Jurisdiction = state,
`June 2025` = percent(rate_2025_06_01, accuracy = 0.1),
`June 2026` = percent(rate_2026_06_01, accuracy = 0.1),
`Change (pp)` = number(change_pp, accuracy = 0.1)
) |>
kable(align = c("l", "r", "r", "r"))
```
```{r}
#| label: tbl-seifa
seifa_change |>
transmute(
Quintile = group,
`June 2025` = percent(june_2025, accuracy = 0.1),
`June 2026` = percent(june_2026, accuracy = 0.1),
`Change (pp)` = number(change_pp, accuracy = 0.1)
) |>
kable(align = c("l", "r", "r", "r"))
```
:::
::: {.callout-note collapse="true" title="Visit-volume and patient-cost tables"}
```{r}
#| label: tbl-working-age-visits
working_age_visit_comparison |>
transmute(
Month = calendar_month,
`Visits 2025` = comma(visits_2025, accuracy = 1),
`Visits 2026` = comma(visits_2026, accuracy = 1),
Change = percent(change, accuracy = 0.1)
) |>
kable(align = c("l", "r", "r", "r"))
```
```{r}
#| label: tbl-patient-cost-area
patient_cost_by_area_group |>
transmute(
`Baseline area group` = group,
`Before: average per service` = dollar(average_per_service_Before, accuracy = 0.01),
`After: average per service` = dollar(average_per_service_After, accuracy = 0.01),
`Change per service` = dollar(average_change, accuracy = 0.01),
`Change among patient-billed services` = dollar(patient_billed_average_change, accuracy = 0.01)
) |>
kable(align = c("l", "r", "r", "r", "r"))
```
```{r}
#| label: tbl-patient-cost-state
state_patient_costs |>
arrange(average_change) |>
transmute(
Jurisdiction = state,
`Before: average per service` = dollar(average_per_service_Before, accuracy = 0.01),
`After: average per service` = dollar(average_per_service_After, accuracy = 0.01),
`Change per service` = dollar(average_change, accuracy = 0.01),
`Change in total contributions` = dollar(contribution_change / 1e6, suffix = "m", accuracy = 0.1)
) |>
kable(align = c("l", "r", "r", "r", "r"))
```
:::
::: {.callout-note collapse="true" title="Models and sensitivity checks"}
The national billing-rate and working-age visit models cover January 2023 to June 2026, leaving eight post-policy months. The primary specifications include linear time, a November 2025 level indicator, a post-policy slope and calendar-month fixed effects. The visit sensitivity removes the post-policy slope. Confidence intervals use a Bartlett-kernel heteroskedasticity-and-autocorrelation-consistent covariance matrix with the automatic Newey–West lag rule. The short post window and number of estimated terms make these supporting checks only; they do not identify a causal effect.
```{r}
#| label: tbl-working-age-models
working_age_visit_models |>
transmute(
Specification = specification,
`Immediate change` = percent(change_percent / 100, accuracy = 0.1),
`HAC 95% interval` = paste0(
percent(lower_percent / 100, accuracy = 0.1), " to ",
percent(upper_percent / 100, accuracy = 0.1)
)
) |>
kable(align = c("l", "r", "r"))
```
The area-level model is:
$$
\log(services_{it}) = \alpha_i + \lambda_t + \beta(Post_t \times Exposure_i) + \varepsilon_{it}
$$
Exposure is an SA3's non-bulk-billed share over the four clean quarters from December 2024 to September 2025, scaled so $\beta$ represents a 10-percentage-point difference. The balanced panel contains `r n_distinct(sa3_model_data$sa3_id)` classified SA3s and two clean post-policy quarters. Standard errors are clustered by SA3.
```{r}
#| label: tbl-sa3-models
sa3_models |>
transmute(
Specification = specification,
`Change per 10 pp higher exposure` = percent(change_percent / 100, accuracy = 0.1),
`Clustered 95% interval` = paste0(
percent(lower_percent / 100, accuracy = 0.1), " to ",
percent(upper_percent / 100, accuracy = 0.1)
)
) |>
kable(align = c("l", "r", "r"))
```
:::
::: {.callout-note collapse="true" title="Sources, checks and reproducibility"}
The monthly series is from the [AIHW Medicare bulk-billing dashboard](https://www.aihw.gov.au/reports/medicare/medicare-bulk-billing-gp-attendances-monthly-data/contents/dashboard), with definitions in its [technical notes](https://www.aihw.gov.au/reports/medicare/medicare-bulk-billing-gp-attendances-monthly-data/contents/technical-notes). AIHW [corrected the downloadable workbook on 31 August 2026](https://www.aihw.gov.au/reports/medicare/medicare-bulk-billing-gp-attendances-monthly-data/notes) after rates had been rounded to whole percentages. The quarterly service and patient-cost data are from the Department's [June-quarter SA3 publication](https://www.health.gov.au/resources/publications/medicare-quarterly-statistics-statistical-area-sa3-summary-june-quarter-2025-26?language=en).
The [Department's commencement notice](https://www.health.gov.au/news/bulk-billing-practice-incentive-program-registration-now-open) documents the 1 November 2025 expansion. The Department also publishes a separate [per-patient bulk-billing dashboard](https://www.health.gov.au/resources/publications/medicare-statistics-per-patient-bulk-billing-dashboard-2025-26).
The code validates coverage dates, panel completeness, rate ranges, numerator-denominator identities and SA3-to-national coverage before drawing charts or fitting models. Downloads are pinned to exact official files and SHA-256 checksums in `sources.csv`. Rendering with `BLOG_REFRESH_DATA=true` deliberately fetches the registered files again; a changed upstream file fails checksum validation before replacing the cache.
```{r}
#| label: tbl-sources
source_registry |>
transmute(
Publisher = publisher,
Dataset = paste0("[", dataset, "](", landing_url, ")"),
`Source/revision date` = release_date,
`Coverage end` = coverage_end,
`SHA-256 prefix` = substr(sha256, 1, 12)
) |>
kable(escape = FALSE, align = c("l", "l", "l", "l", "l"))
```
One source oddity is recorded rather than hidden: the Department's June-quarter publication currently links to a file whose URL still says “March quarter”. The workbook itself contains `r n_distinct(sa3_raw$quarter_date)` quarters through June 2026, and the parser stops if that coverage changes.
:::