Session info
sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
## [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.utf8
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] jsonlite_2.0.0 servr_0.32 reader_1.1.0 NCmisc_1.3.1
## [5] patchwork_1.3.2 scales_1.4.0 kableExtra_1.4.0 knitr_1.51
## [9] janitor_2.2.1 lubridate_1.9.4 forcats_1.0.1 stringr_1.6.0
## [13] dplyr_1.1.4 purrr_1.2.0 readr_2.2.0 tidyr_1.3.2
## [17] tibble_3.3.0 ggplot2_4.0.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 xml2_1.5.1 stringi_1.8.7
## [5] hms_1.1.4 digest_0.6.37 magrittr_2.0.4 evaluate_1.0.5
## [9] grid_4.4.2 timechange_0.3.0 RColorBrewer_1.1-3 fastmap_1.2.0
## [13] promises_1.5.0 viridisLite_0.4.2 textshaping_1.0.4 jquerylib_0.1.4
## [17] cli_3.6.5 crayon_1.5.3 rlang_1.1.6 bit64_4.6.0-1
## [21] withr_3.0.2 cachem_1.1.0 yaml_2.3.12 otel_0.2.0
## [25] parallel_4.4.2 tools_4.4.2 tzdb_0.5.0 httpuv_1.6.17
## [29] vctrs_0.6.5 R6_2.6.1 lifecycle_1.0.4 snakecase_0.11.1
## [33] bit_4.6.0 vroom_1.7.1 pkgconfig_2.0.3 later_1.4.8
## [37] pillar_1.11.1 bslib_0.9.0 gtable_0.3.6 Rcpp_1.1.1
## [41] glue_1.8.0 systemfonts_1.3.1 xfun_0.55 tidyselect_1.2.1
## [45] rstudioapi_0.18.0 farver_2.1.2 htmltools_0.5.9 labeling_0.4.3
## [49] rmarkdown_2.30 svglite_2.2.2 compiler_4.4.2 S7_0.2.0
- Use rates (2025 baseline) - NH residents divided by
population by age band per municipality. These are the age-specific
admission probabilities.
# ── National 2025 elderly population from projection data ───────────────────
pop_2025 <- projections_banded |>
filter(year == 2025, scenario == "Strong ageing") |>
# Both scenarios have identical 2025 values — it's the base year
group_by(age_band) |>
summarise(population_2025 = sum(population, na.rm = TRUE), .groups = "drop")
cat("National elderly population 2025:\n")
## National elderly population 2025:
print(pop_2025)
## # A tibble: 3 × 2
## age_band population_2025
## <chr> <dbl>
## 1 67-79 667683
## 2 80-89 225329
## 3 90+ 45319
# ── National long-term institutional demand 2025 ─────────────────────────────
national_lt_demand <- users |>
summarise(
total_lt_stay = sum(long_term_stay_in_institutions, na.rm = TRUE),
total_beds = sum(institution_all_available_beds, na.rm = TRUE),
total_dementia = sum(dementia_beds, na.rm = TRUE)
)
cat("\nNational long-term institutional demand 2025:\n")
##
## National long-term institutional demand 2025:
print(national_lt_demand)
## # A tibble: 1 × 3
## total_lt_stay total_beds total_dementia
## <dbl> <dbl> <dbl>
## 1 30622 41650. 11704
# ── National use rates by age band ───────────────────────────────────────────
# We cannot split long_term_stay by age band from this data directly.
# Instead we use the known national age distribution of NH residents
# from KOSTRA 2025 (from your Residence in institutions sheet in Excel):
# 67-79: ~9,719 residents
# 80-89: ~15,144 residents
# 90+: ~10,380 residents
# Adjust these numbers if your Excel data shows different figures.
lt_by_age <- tibble(
age_band = c("67-79", "80-89", "90+"),
lt_residents_2025 = c(9719, 15144, 10380) # <-- update from your Excel if needed
)
use_rates <- lt_by_age |>
left_join(pop_2025, by = "age_band") |>
mutate(
use_rate = lt_residents_2025 / population_2025
)
cat("\nAge-specific use rates (probability of being in LT institutional care):\n")
##
## Age-specific use rates (probability of being in LT institutional care):
print(use_rates)
## # A tibble: 3 × 4
## age_band lt_residents_2025 population_2025 use_rate
## <chr> <dbl> <dbl> <dbl>
## 1 67-79 9719 667683 0.0146
## 2 80-89 15144 225329 0.0672
## 3 90+ 10380 45319 0.229
# ── Corrected LT institutional use rates ─────────────────────────────────────
# Scale the age-band breakdown to match the actual total_lt_stay from the data
actual_total_lt <- national_lt_demand$total_lt_stay # 30,622
assumed_total <- 9719 + 15144 + 10380 # 35,243
scale_factor <- actual_total_lt / assumed_total
lt_by_age <- tibble(
age_band = c("67-79", "80-89", "90+"),
lt_residents_2025 = round(c(9719, 15144, 10380) * scale_factor)
)
use_rates_institutional <- lt_by_age |>
left_join(pop_2025, by = "age_band") |>
mutate(
use_rate_institutional = lt_residents_2025 / population_2025
)
cat("Corrected institutional use rates:\n")
## Corrected institutional use rates:
print(use_rates_institutional)
## # A tibble: 3 × 4
## age_band lt_residents_2025 population_2025 use_rate_institutional
## <chr> <dbl> <dbl> <dbl>
## 1 67-79 8445 667683 0.0126
## 2 80-89 13158 225329 0.0584
## 3 90+ 9019 45319 0.199
# ── Home care use rates ───────────────────────────────────────────────────────
# We need elderly home care users by age band.
# From users dataset: home_care_elderly = total - under 67
# But we need the 67-79 / 80-89 / 90+ split.
# We can get this from the projection files indirectly via the services data,
# but we do not have it split by age in services_raw.
# Best available proxy: use the age distribution from SSB table 11645
# which shows home nursing users by age (from your Users sheet in the Excel).
# Known approximate national split for home-based services elderly users:
# 67-79: ~36% of elderly home care users
# 80-89: ~44% of elderly home care users
# 90+: ~20% of elderly home care users
# These come from your DifferentDatas sheet analysis earlier in this project.
total_home_care_elderly <- sum(users$home_care_elderly, na.rm = TRUE)
cat("\nTotal elderly home care users (national):", comma(total_home_care_elderly), "\n")
##
## Total elderly home care users (national): 113,343
# Also compute total home nursing and practical assistance from services data
total_home_nursing <- sum(users$home_nursing, na.rm = TRUE)
total_practical <- sum(users$practical_assistance_daily_activities, na.rm = TRUE)
cat("Total home nursing users: ", comma(total_home_nursing), "\n")
## Total home nursing users: 168,586
cat("Total practical assistance users: ", comma(total_practical), "\n")
## Total practical assistance users: 75,068
# Age-band split for home care
hc_by_age <- tibble(
age_band = c("67-79", "80-89", "90+"),
hc_share = c(0.367, 0.298, 0.247), # from your DifferentDatas analysis
hc_users_2025 = round(total_home_care_elderly * hc_share / sum(c(0.367, 0.298, 0.247)))
)
use_rates_homecare <- hc_by_age |>
left_join(pop_2025, by = "age_band") |>
mutate(
use_rate_homecare = hc_users_2025 / population_2025
)
cat("\nHome care use rates by age band:\n")
##
## Home care use rates by age band:
print(use_rates_homecare)
## # A tibble: 3 × 5
## age_band hc_share hc_users_2025 population_2025 use_rate_homecare
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 67-79 0.367 45611 667683 0.0683
## 2 80-89 0.298 37035 225329 0.164
## 3 90+ 0.247 30697 45319 0.677
# ── Combined use rate summary ─────────────────────────────────────────────────
use_rates_combined <- use_rates_institutional |>
select(age_band, lt_residents_2025, use_rate_institutional) |>
left_join(
use_rates_homecare |> select(age_band, hc_users_2025, use_rate_homecare),
by = "age_band"
) |>
left_join(pop_2025, by = "age_band") |>
mutate(
total_care_users = lt_residents_2025 + hc_users_2025,
use_rate_any_care = total_care_users / population_2025,
pct_institutional = lt_residents_2025 / total_care_users
)
cat("\nCombined use rate summary:\n")
##
## Combined use rate summary:
use_rates_combined |>
mutate(across(starts_with("use_rate"), ~ percent(.x, accuracy = 0.1)),
across(where(is.numeric), ~ comma(.x))) |>
kable(caption = "2025 care service use rates by elderly age band",
col.names = c("Age band", "LT institutional users", "Use rate (institutional)",
"Home care users", "Use rate (home care)",
"Population 2025", "Total care users",
"Use rate (any care)", "Share in institutions")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
2025 care service use rates by elderly age band
|
Age band
|
LT institutional users
|
Use rate (institutional)
|
Home care users
|
Use rate (home care)
|
Population 2025
|
Total care users
|
Use rate (any care)
|
Share in institutions
|
|
67-79
|
8,445
|
1.3%
|
45,611
|
6.8%
|
667,683
|
54,056
|
8.1%
|
0.156
|
|
80-89
|
13,158
|
5.8%
|
37,035
|
16.4%
|
225,329
|
50,193
|
22.3%
|
0.262
|
|
90+
|
9,019
|
19.9%
|
30,697
|
67.7%
|
45,319
|
39,716
|
87.6%
|
0.227
|
- Demand projection - apply use rates to projected
populations year by year under both ageing scenarios.
# ════════════════════════════════════════════════════════════════════════════
# DEMAND PROJECTION
# ════════════════════════════════════════════════════════════════════════════
# ── Parameters (adjust these to test different scenarios) ───────────────────
SUPPLY_GROWTH_RATE <- 0.004 # historical +0.4% per year
EXPANSION_ANNUAL <- 500 # extra beds added per year in expansion scenario
MIN_POP_THRESHOLD <- 50 # minimum elderly pop per band for stable use rate
# ── Step 1: National use rates (already computed above) ─────────────────────
# use_rates_institutional has: age_band, use_rate_institutional
national_use_rates <- use_rates_institutional |>
select(age_band, use_rate = use_rate_institutional)
# ── Step 2: National demand projection ──────────────────────────────────────
national_demand <- national_proj |>
left_join(national_use_rates, by = "age_band") |>
mutate(projected_lt_users = population * use_rate) |>
group_by(scenario, year) |>
summarise(
projected_demand = sum(projected_lt_users, na.rm = TRUE),
.groups = "drop"
)
# National supply scenarios
national_supply_base <- sum(users$institution_all_available_beds, na.rm = TRUE)
national_supply <- tibble(year = 2025:2036) |>
mutate(
flat = national_supply_base,
historical = national_supply_base * (1 + SUPPLY_GROWTH_RATE)^(year - 2025),
expansion = national_supply_base + EXPANSION_ANNUAL * (year - 2025)
) |>
pivot_longer(-year, names_to = "supply_scenario", values_to = "supply")
# ── Combine demand and supply ────────────────────────────────────────────────
national_combined <- national_demand |>
left_join(national_supply, by = "year",
relationship= "many-to-many") |>
mutate(
gap = projected_demand - supply,
over_capacity = gap > 0
)
cat("National supply baseline (2025):", comma(national_supply_base), "beds\n")
## National supply baseline (2025): 41,650 beds
cat("National LT demand (2025): ",
comma(round(
national_combined |>
filter(year == 2025,
scenario == "Strong ageing",
supply_scenario == "flat") |>
pull(projected_demand)
)), "\n")
## National LT demand (2025): 30,622
# Clean join
national_combined <- national_demand |>
left_join(national_supply, by = "year",
relationship = "many-to-many") |>
mutate(
gap = projected_demand - supply,
over_capacity = gap > 0
)
cat("National supply baseline (2025):", comma(national_supply_base), "beds\n")
## National supply baseline (2025): 41,650 beds
cat("National LT demand (2025): ",
comma(round(filter(national_combined, year == 2025,
scenario == "Strong ageing",
supply_scenario == "flat")$projected_demand)), "\n")
## National LT demand (2025): 30,622
# ── Step 3: National crossover year ─────────────────────────────────────────
national_crossover <- national_combined |>
filter(over_capacity) |>
group_by(scenario, supply_scenario) |>
summarise(crossover_year = min(year), .groups = "drop")
cat("\nNational crossover years (demand exceeds supply):\n")
##
## National crossover years (demand exceeds supply):
national_crossover |>
pivot_wider(names_from = supply_scenario,
values_from = crossover_year,
values_fill = NA) |>
kable(caption = "First year projected LT demand exceeds supply — national") |>
kable_styling(bootstrap_options = c("striped", "hover"))
First year projected LT demand exceeds supply — national
|
scenario
|
flat
|
historical
|
|
Strong ageing
|
2034
|
2034
|
# ── Step 4: Municipality use rates ──────────────────────────────────────────
# For each municipality: use rate = lt_residents_2025 / population_2025
# We do not have lt_residents split by age band per municipality,
# so we apply the NATIONAL age-band use rates scaled by the municipality's
# overall lt intensity (their total LT stay / national average LT stay per capita)
national_lt_per_capita <- sum(users$long_term_stay_in_institutions, na.rm = TRUE) /
sum(pop_2025$population_2025)
muni_use_rates <- users |>
select(muni_code, muni_name, long_term_stay_in_institutions,
institution_all_available_beds) |>
mutate(
# Municipality-specific scaling factor relative to national average
lt_per_capita = long_term_stay_in_institutions /
sum(long_term_stay_in_institutions, na.rm = TRUE) *
sum(pop_2025$population_2025),
muni_scale = (long_term_stay_in_institutions /
replace_na(long_term_stay_in_institutions, 0)) /
(national_lt_per_capita * 1),
) |>
# Simpler and more robust: just compute LT users per elderly resident
# using the municipality's own total elderly population from projections
select(muni_code, muni_name, long_term_stay_in_institutions,
supply_2025 = institution_all_available_beds)
# Get 2025 elderly population per municipality from projections
muni_pop_2025 <- projections_banded |>
filter(year == 2025, scenario == "Strong ageing") |>
group_by(muni_code, age_band) |>
summarise(population = sum(population, na.rm = TRUE), .groups = "drop")
muni_pop_2025_total <- muni_pop_2025 |>
group_by(muni_code) |>
summarise(elderly_pop_2025 = sum(population), .groups = "drop")
# Join and compute municipality use rate (aggregate, not age-banded)
muni_use_rates <- muni_use_rates |>
left_join(muni_pop_2025_total, by = "muni_code") |>
mutate(
use_rate_muni = long_term_stay_in_institutions / elderly_pop_2025,
# Flag municipalities with too few elderly for stable rate
stable = elderly_pop_2025 >= MIN_POP_THRESHOLD
) |>
filter(!is.na(use_rate_muni))
cat("\nMunicipality use rates — summary:\n")
##
## Municipality use rates — summary:
summary(muni_use_rates$use_rate_muni)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.005995 0.027099 0.032022 0.033312 0.038778 0.085938
# ── Step 5: Municipality demand projection ───────────────────────────────────
# For each municipality: apply their own use rate to projected elderly pop
# For small municipalities (stable == FALSE): use national age-band rates instead
muni_proj_pop <- projections_banded |>
group_by(muni_code, scenario, year) |>
summarise(elderly_pop = sum(population, na.rm = TRUE), .groups = "drop")
muni_demand <- muni_proj_pop |>
left_join(muni_use_rates |> select(muni_code, use_rate_muni,
supply_2025, stable, muni_name),
by = "muni_code") |>
# Fall back to national aggregate use rate for unstable municipalities
mutate(
use_rate_applied = if_else(
stable,
use_rate_muni,
sum(users$long_term_stay_in_institutions, na.rm = TRUE) /
sum(muni_pop_2025_total$elderly_pop_2025, na.rm = TRUE)
),
projected_demand = elderly_pop * use_rate_applied
)
# ── Step 6: Municipality supply scenarios ────────────────────────────────────
muni_supply <- muni_demand |>
distinct(muni_code, supply_2025) |>
crossing(year = 2025:2036) |>
mutate(
flat = supply_2025,
historical = supply_2025 * (1 + SUPPLY_GROWTH_RATE)^(year - 2025),
expansion = supply_2025 + (EXPANSION_ANNUAL *
(supply_2025 / national_supply_base)) * (year - 2025)
) |>
pivot_longer(c(flat, historical, expansion),
names_to = "supply_scenario", values_to = "supply")
muni_combined <- muni_demand |>
left_join(muni_supply, by = c("muni_code", "year")) |>
mutate(
gap = projected_demand - supply,
over_capacity = gap > 0
)
# ── Step 7: Municipality crossover years ─────────────────────────────────────
muni_crossover <- muni_combined |>
filter(over_capacity) |>
group_by(muni_code, muni_name, scenario, supply_scenario) |>
summarise(crossover_year = min(year), .groups = "drop")
# Summary table: how many municipalities cross over by year
cat("\nMunicipalities crossing over — national summary (flat supply, strong ageing):\n")
##
## Municipalities crossing over — national summary (flat supply, strong ageing):
muni_crossover |>
filter(scenario == "Strong ageing", supply_scenario == "flat") |>
count(crossover_year) |>
arrange(crossover_year) |>
mutate(cumulative = cumsum(n)) |>
kable(caption = "Number of municipalities where demand exceeds supply (flat supply, strong ageing)",
col.names = c("Year", "New crossovers", "Cumulative")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Number of municipalities where demand exceeds supply (flat supply,
strong ageing)
|
Year
|
New crossovers
|
Cumulative
|
|
2025
|
1
|
1
|
|
2026
|
3
|
4
|
|
2027
|
1
|
5
|
|
2028
|
1
|
6
|
|
2029
|
2
|
8
|
|
2030
|
4
|
12
|
|
2031
|
7
|
19
|
|
2032
|
4
|
23
|
|
2033
|
9
|
32
|
|
2034
|
8
|
40
|
# Top 20 most urgent municipalities
cat("\nMost urgent municipalities (earliest crossover, flat supply, strong ageing):\n")
##
## Most urgent municipalities (earliest crossover, flat supply, strong ageing):
muni_crossover |>
filter(scenario == "Strong ageing", supply_scenario == "flat") |>
arrange(crossover_year) |>
head(20) |>
left_join(users |> select(muni_code, institution_all_available_beds,
long_term_stay_in_institutions),
by = "muni_code") |>
kable(caption = "20 most urgent municipalities — earliest demand-supply crossover",
format.args = list(big.mark = ",")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
20 most urgent municipalities — earliest demand-supply crossover
|
muni_code
|
muni_name
|
scenario
|
supply_scenario
|
crossover_year
|
institution_all_available_beds
|
long_term_stay_in_institutions
|
|
5,536
|
Lyngen - Ivgu - Yyke?
|
Strong ageing
|
flat
|
2,025
|
24.0
|
29
|
|
1,145
|
Bokn
|
Strong ageing
|
flat
|
2,026
|
9.0
|
9
|
|
1,840
|
Saltdal
|
Strong ageing
|
flat
|
2,026
|
41.0
|
41
|
|
5,620
|
Nordkapp
|
Strong ageing
|
flat
|
2,026
|
28.0
|
28
|
|
5,622
|
Porsanger - Pors?ngu - Porsanki
|
Strong ageing
|
flat
|
2,027
|
36.0
|
34
|
|
5,610
|
K?r?sjohka - Karasjok
|
Strong ageing
|
flat
|
2,028
|
16.0
|
15
|
|
301
|
Oslo - Oslove
|
Strong ageing
|
flat
|
2,029
|
4,121.1
|
3,759
|
|
1,828
|
Nesna
|
Strong ageing
|
flat
|
2,029
|
15.0
|
14
|
|
1,532
|
Giske
|
Strong ageing
|
flat
|
2,030
|
24.0
|
21
|
|
3,222
|
L?renskog
|
Strong ageing
|
flat
|
2,030
|
240.0
|
206
|
|
3,336
|
Rollag
|
Strong ageing
|
flat
|
2,030
|
16.0
|
15
|
|
5,544
|
Nordreisa - R?isa - Raisi
|
Strong ageing
|
flat
|
2,030
|
52.0
|
48
|
|
1,160
|
Vindafjord
|
Strong ageing
|
flat
|
2,031
|
81.0
|
70
|
|
1,838
|
Gildesk?l
|
Strong ageing
|
flat
|
2,031
|
32.0
|
29
|
|
3,112
|
R?de
|
Strong ageing
|
flat
|
2,031
|
51.0
|
46
|
|
3,238
|
Nannestad
|
Strong ageing
|
flat
|
2,031
|
68.0
|
56
|
|
4,614
|
Stord
|
Strong ageing
|
flat
|
2,031
|
109.0
|
95
|
|
5,029
|
Skaun
|
Strong ageing
|
flat
|
2,031
|
50.0
|
42
|
|
5,534
|
Karls?y
|
Strong ageing
|
flat
|
2,031
|
42.0
|
38
|
|
3,220
|
Enebakk
|
Strong ageing
|
flat
|
2,032
|
71.0
|
57
|
- Back-calibration - what reduction in use rate keeps
demand below the supply ceiling? Expressed as N patients per year that
must stay in home care.
# ════════════════════════════════════════════════════════════════════════════
# BACK-CALIBRATION
# How many institutional admissions need to be prevented each year?
# ════════════════════════════════════════════════════════════════════════════
# ── Step 1: National back-calibration ───────────────────────────────────────
# ── Step 1: Total projected elderly population by scenario and year ──────────
national_elderly_proj <- national_proj |>
group_by(scenario, year) |>
summarise(total_elderly = sum(population, na.rm = TRUE), .groups = "drop")
# ── Step 2: National back-calibration table ──────────────────────────────────
national_backcalib <- national_combined |>
filter(supply_scenario == "flat") |>
left_join(national_elderly_proj, by = c("scenario", "year")) |>
mutate(
# Use rate that would keep demand exactly at supply
required_use_rate = supply / total_elderly,
# Current use rate (what we actually observe)
current_use_rate = projected_demand / total_elderly,
# How many admissions need to be prevented (positive = action needed)
admissions_to_prevent = projected_demand - supply,
# What percentage reduction in admissions does that represent
pct_reduction_needed = admissions_to_prevent / projected_demand * 100,
# Use rate reduction needed (percentage points)
use_rate_reduction_pp = (current_use_rate - required_use_rate) * 100
)
# Print national summary
cat("National back-calibration (flat supply, both scenarios):\n")
## National back-calibration (flat supply, both scenarios):
national_backcalib |>
filter(year >= 2028) |>
select(scenario, year, projected_demand, supply,
admissions_to_prevent, pct_reduction_needed) |>
mutate(
projected_demand = comma(round(projected_demand)),
supply = comma(round(supply)),
admissions_to_prevent = comma(round(admissions_to_prevent)),
pct_reduction_needed = round(pct_reduction_needed, 1)
) |>
kable(
caption = "National back-calibration — admissions to prevent per year (flat supply)",
col.names = c("Scenario", "Year", "Projected demand",
"Supply", "Admissions to prevent", "% reduction needed")
) |>
kable_styling(bootstrap_options = c("striped", "hover")) |>
row_spec(which(as.numeric(
str_extract(national_backcalib |>
filter(year >= 2028) |>
pull(admissions_to_prevent), "-?\\d+"
)) > 0), background = "#FCE4D6")
National back-calibration — admissions to prevent per year (flat supply)
|
Scenario
|
Year
|
Projected demand
|
Supply
|
Admissions to prevent
|
% reduction needed
|
|
Strong ageing
|
2028
|
34,446
|
41,650
|
-7,205
|
-20.9
|
|
Strong ageing
|
2029
|
35,861
|
41,650
|
-5,790
|
-16.1
|
|
Strong ageing
|
2030
|
37,267
|
41,650
|
-4,384
|
-11.8
|
|
Strong ageing
|
2031
|
38,704
|
41,650
|
-2,946
|
-7.6
|
|
Strong ageing
|
2032
|
40,015
|
41,650
|
-1,636
|
-4.1
|
|
Strong ageing
|
2033
|
41,636
|
41,650
|
-14
|
0.0
|
|
Strong ageing
|
2034
|
43,405
|
41,650
|
1,755
|
4.0
|
|
Weak ageing
|
2028
|
33,292
|
41,650
|
-8,358
|
-25.1
|
|
Weak ageing
|
2029
|
34,298
|
41,650
|
-7,352
|
-21.4
|
|
Weak ageing
|
2030
|
35,280
|
41,650
|
-6,370
|
-18.1
|
|
Weak ageing
|
2031
|
36,242
|
41,650
|
-5,409
|
-14.9
|
|
Weak ageing
|
2032
|
37,050
|
41,650
|
-4,600
|
-12.4
|
|
Weak ageing
|
2033
|
38,106
|
41,650
|
-3,544
|
-9.3
|
|
Weak ageing
|
2034
|
39,236
|
41,650
|
-2,415
|
-6.2
|
# ── Step 3: Municipality back-calibration ────────────────────────────────────
muni_elderly_proj <- projections_banded |>
group_by(muni_code, scenario, year) |>
summarise(total_elderly = sum(population, na.rm = TRUE), .groups = "drop")
muni_backcalib <- muni_combined |>
filter(supply_scenario == "flat") |>
left_join(
muni_elderly_proj,
by = c("muni_code", "scenario", "year")
) |>
left_join(muni_lookup, by = "muni_code") |>
# Clean up duplicate columns from the joins
rename(scenario = scenario,
muni_name = muni_name.y) |>
select(-any_of(c("muni_name.x", "scenario.x", "scenario.y",
"muni_name", "supply_2025.x", "supply_2025.y"))) |>
mutate(
required_use_rate = supply / total_elderly,
current_use_rate = projected_demand / total_elderly,
admissions_to_prevent = projected_demand - supply,
pct_reduction_needed = admissions_to_prevent / projected_demand * 100,
use_rate_reduction_pp = (current_use_rate - required_use_rate) * 100
)
muni_backcalib <- muni_combined |>
filter(supply_scenario == "flat") |>
left_join(
muni_elderly_proj,
by = c("muni_code", "scenario", "year")
) |>
select(-any_of(c("muni_name", "muni_name.x", "muni_name.y",
"supply_2025.x", "supply_2025.y",
"scenario.x", "scenario.y"))) |>
mutate(
required_use_rate = supply / total_elderly,
current_use_rate = projected_demand / total_elderly,
admissions_to_prevent = projected_demand - supply,
pct_reduction_needed = admissions_to_prevent / projected_demand * 100,
use_rate_reduction_pp = (current_use_rate - required_use_rate) * 100
) |>
left_join(muni_lookup, by = "muni_code") # add name at the very end
# ── Step 4: Municipality summary table ───────────────────────────────────────
# Show each municipality's situation in their crossover year
# (the first year action is needed)
muni_action_needed <- muni_backcalib |>
filter(admissions_to_prevent > 0, scenario == "Strong ageing") |>
group_by(muni_code, muni_name) |>
slice_min(year, n = 1) |> # first year action is needed
ungroup() |>
arrange(year, desc(admissions_to_prevent)) |>
select(muni_code, muni_name, year,
total_elderly, projected_demand, supply,
admissions_to_prevent, pct_reduction_needed,
use_rate_reduction_pp)
cat("\nMunicipalities requiring action — first year and scale of intervention:\n")
##
## Municipalities requiring action — first year and scale of intervention:
muni_action_needed |>
mutate(
total_elderly = comma(round(total_elderly)),
projected_demand = comma(round(projected_demand)),
supply = comma(round(supply)),
admissions_to_prevent = comma(round(admissions_to_prevent)),
pct_reduction_needed = paste0(round(pct_reduction_needed, 1), "%"),
use_rate_reduction_pp = paste0(round(use_rate_reduction_pp, 3), " pp")
) |>
select(-muni_code) |>
kable(
caption = "Municipality back-calibration — action required (strong ageing, flat supply)",
col.names = c("Municipality", "Action needed from", "Elderly population",
"Projected demand", "Supply", "Admissions to prevent",
"% reduction needed", "Use rate reduction (pp)")
) |>
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = TRUE, font_size = 11)
Municipality back-calibration — action required (strong ageing, flat
supply)
|
Municipality
|
Action needed from
|
Elderly population
|
Projected demand
|
Supply
|
Admissions to prevent
|
% reduction needed
|
Use rate reduction (pp)
|
|
Lyngen - Ivgu - Yyke?
|
2025
|
748
|
29
|
24
|
5
|
17.2%
|
0.668 pp
|
|
Bokn
|
2026
|
201
|
10
|
9
|
1
|
6%
|
0.284 pp
|
|
Saltdal
|
2026
|
1,133
|
42
|
41
|
1
|
1.3%
|
0.049 pp
|
|
Nordkapp
|
2026
|
619
|
28
|
28
|
0
|
1.5%
|
0.067 pp
|
|
Porsanger - Pors?ngu - Porsanki
|
2027
|
911
|
36
|
36
|
0
|
0.2%
|
0.006 pp
|
|
K?r?sjohka - Karasjok
|
2028
|
563
|
17
|
16
|
1
|
3.9%
|
0.117 pp
|
|
Oslo - Oslove
|
2029
|
95,450
|
4,217
|
4,121
|
95
|
2.3%
|
0.1 pp
|
|
Nesna
|
2029
|
406
|
15
|
15
|
0
|
2.6%
|
0.099 pp
|
|
L?renskog
|
2030
|
8,071
|
244
|
240
|
4
|
1.7%
|
0.052 pp
|
|
Giske
|
2030
|
1,625
|
24
|
24
|
0
|
0.8%
|
0.012 pp
|
|
Rollag
|
2030
|
408
|
16
|
16
|
0
|
0.9%
|
0.036 pp
|
|
Nordreisa - R?isa - Raisi
|
2030
|
1,092
|
52
|
52
|
0
|
0.1%
|
0.005 pp
|
|
Nannestad
|
2031
|
2,476
|
69
|
68
|
1
|
1.1%
|
0.03 pp
|
|
Stord
|
2031
|
3,716
|
110
|
109
|
1
|
0.7%
|
0.02 pp
|
|
R?de
|
2031
|
1,759
|
52
|
51
|
1
|
1.4%
|
0.04 pp
|
|
Vindafjord
|
2031
|
1,855
|
81
|
81
|
0
|
0.4%
|
0.019 pp
|
|
Karls?y
|
2031
|
679
|
42
|
42
|
0
|
0.4%
|
0.024 pp
|
|
Gildesk?l
|
2031
|
591
|
32
|
32
|
0
|
0.3%
|
0.016 pp
|
|
Skaun
|
2031
|
1,490
|
50
|
50
|
0
|
0.1%
|
0.004 pp
|
|
Nittedal
|
2032
|
4,420
|
125
|
124
|
2
|
1.3%
|
0.037 pp
|
|
Enebakk
|
2032
|
2,135
|
72
|
71
|
1
|
1.9%
|
0.063 pp
|
|
B?tsfjord
|
2032
|
428
|
17
|
16
|
1
|
3.6%
|
0.14 pp
|
|
Bardu
|
2032
|
832
|
37
|
37
|
0
|
0%
|
0 pp
|
|
Lillestr?m
|
2033
|
16,892
|
490
|
473
|
17
|
3.4%
|
0.099 pp
|
|
Stavanger
|
2033
|
26,841
|
1,136
|
1,119
|
17
|
1.5%
|
0.062 pp
|
|
Sandnes
|
2033
|
14,180
|
523
|
518
|
5
|
0.9%
|
0.034 pp
|
|
Nes
|
2033
|
5,067
|
129
|
125
|
4
|
2.8%
|
0.071 pp
|
|
Trondheim - Tr?ante
|
2033
|
37,366
|
1,545
|
1,542
|
3
|
0.2%
|
0.009 pp
|
|
R?lingen
|
2033
|
3,340
|
89
|
86
|
3
|
3.1%
|
0.081 pp
|
|
Sola
|
2033
|
4,904
|
188
|
187
|
1
|
0.4%
|
0.014 pp
|
|
Gausdal
|
2033
|
1,599
|
63
|
63
|
0
|
0.8%
|
0.031 pp
|
|
Hol
|
2033
|
1,199
|
39
|
39
|
0
|
0.5%
|
0.018 pp
|
|
Kristiansand
|
2034
|
22,932
|
735
|
726
|
9
|
1.2%
|
0.04 pp
|
|
Bergen
|
2034
|
54,563
|
2,615
|
2,609
|
6
|
0.2%
|
0.011 pp
|
|
Troms?
|
2034
|
14,272
|
478
|
474
|
4
|
0.8%
|
0.027 pp
|
|
Vennesla
|
2034
|
3,043
|
101
|
100
|
2
|
1.6%
|
0.052 pp
|
|
Halden
|
2034
|
7,807
|
219
|
218
|
1
|
0.5%
|
0.013 pp
|
|
Alstahaug
|
2034
|
1,839
|
57
|
56
|
1
|
1.5%
|
0.048 pp
|
|
Sunndal
|
2034
|
1,851
|
78
|
78
|
0
|
0.2%
|
0.008 pp
|
|
Lur?y
|
2034
|
502
|
13
|
13
|
0
|
0.1%
|
0.002 pp
|
# ── Step 4b: Municipality trajectory — 5 years beyond crossover ──────────────
# Get each municipality's crossover year
muni_crossover_year <- muni_backcalib |>
filter(admissions_to_prevent > 0, scenario == "Strong ageing") |>
group_by(muni_code) |>
summarise(crossover_year = min(year), .groups = "drop")
# Filter backcalib to crossover year through min(crossover + 5, 2036)
muni_trajectory <- muni_backcalib |>
filter(scenario == "Strong ageing") |>
left_join(muni_crossover_year, by = "muni_code") |>
filter(
!is.na(crossover_year), # only municipalities that cross over
year >= crossover_year, # from crossover year onwards
year <= pmin(crossover_year + 5, 2036) # up to 5 years after or 2036
) |>
arrange(crossover_year, muni_code, year) |>
select(muni_code, muni_name, crossover_year, year,
total_elderly, projected_demand, supply,
admissions_to_prevent, pct_reduction_needed,
use_rate_reduction_pp)
cat("\nMunicipality trajectories from crossover year (strong ageing, flat supply):\n")
##
## Municipality trajectories from crossover year (strong ageing, flat supply):
muni_trajectory |>
mutate(
total_elderly = comma(round(total_elderly)),
projected_demand = comma(round(projected_demand)),
supply = comma(round(supply)),
admissions_to_prevent = comma(round(admissions_to_prevent)),
pct_reduction_needed = paste0(round(pct_reduction_needed, 1), "%"),
use_rate_reduction_pp = paste0(round(use_rate_reduction_pp, 3), " pp")
) |>
select(-muni_code, -crossover_year) |>
kable(
caption = "Municipality trajectories from crossover year + 5 years (strong ageing, flat supply)",
col.names = c("Municipality", "Year", "Elderly population",
"Projected demand", "Supply", "Admissions to prevent",
"% reduction needed", "Use rate reduction (pp)")
) |>
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = TRUE, font_size = 11) |>
collapse_rows(columns = 1, valign = "top")
Municipality trajectories from crossover year + 5 years (strong ageing,
flat supply)
|
Municipality
|
Year
|
Elderly population
|
Projected demand
|
Supply
|
Admissions to prevent
|
% reduction needed
|
Use rate reduction (pp)
|
|
Lyngen - Ivgu - Yyke?
|
2025
|
748
|
29
|
24
|
5
|
17.2%
|
0.668 pp
|
|
2026
|
753
|
29
|
24
|
5
|
17.8%
|
0.69 pp
|
|
2027
|
768
|
30
|
24
|
6
|
19.4%
|
0.752 pp
|
|
2028
|
777
|
30
|
24
|
6
|
20.3%
|
0.788 pp
|
|
2029
|
784
|
30
|
24
|
6
|
21%
|
0.816 pp
|
|
2030
|
786
|
30
|
24
|
6
|
21.2%
|
0.824 pp
|
|
Bokn
|
2026
|
201
|
10
|
9
|
1
|
6%
|
0.284 pp
|
|
2027
|
203
|
10
|
9
|
1
|
6.9%
|
0.328 pp
|
|
2028
|
208
|
10
|
9
|
1
|
9.1%
|
0.435 pp
|
|
2029
|
210
|
10
|
9
|
1
|
10%
|
0.476 pp
|
|
2030
|
209
|
10
|
9
|
1
|
9.6%
|
0.456 pp
|
|
2031
|
211
|
10
|
9
|
1
|
10.4%
|
0.497 pp
|
|
Saltdal
|
2026
|
1,133
|
42
|
41
|
1
|
1.3%
|
0.049 pp
|
|
2027
|
1,139
|
42
|
41
|
1
|
1.8%
|
0.068 pp
|
|
2028
|
1,148
|
42
|
41
|
1
|
2.6%
|
0.096 pp
|
|
2029
|
1,170
|
43
|
41
|
2
|
4.4%
|
0.163 pp
|
|
2030
|
1,205
|
44
|
41
|
3
|
7.2%
|
0.265 pp
|
|
2031
|
1,216
|
45
|
41
|
4
|
8.1%
|
0.296 pp
|
|
Nordkapp
|
2026
|
619
|
28
|
28
|
0
|
1.5%
|
0.067 pp
|
|
2027
|
638
|
29
|
28
|
1
|
4.4%
|
0.201 pp
|
|
2028
|
657
|
30
|
28
|
2
|
7.2%
|
0.328 pp
|
|
2029
|
671
|
31
|
28
|
3
|
9.1%
|
0.417 pp
|
|
2030
|
682
|
31
|
28
|
3
|
10.6%
|
0.485 pp
|
|
2031
|
699
|
32
|
28
|
4
|
12.7%
|
0.584 pp
|
|
Porsanger - Pors?ngu - Porsanki
|
2027
|
911
|
36
|
36
|
0
|
0.2%
|
0.006 pp
|
|
2028
|
926
|
37
|
36
|
1
|
1.8%
|
0.07 pp
|
|
2029
|
945
|
37
|
36
|
1
|
3.8%
|
0.149 pp
|
|
2030
|
965
|
38
|
36
|
2
|
5.7%
|
0.228 pp
|
|
2031
|
980
|
39
|
36
|
3
|
7.2%
|
0.285 pp
|
|
2032
|
1,002
|
40
|
36
|
4
|
9.2%
|
0.365 pp
|
|
K?r?sjohka - Karasjok
|
2028
|
563
|
17
|
16
|
1
|
3.9%
|
0.117 pp
|
|
2029
|
578
|
17
|
16
|
1
|
6.4%
|
0.19 pp
|
|
2030
|
589
|
17
|
16
|
1
|
8.2%
|
0.242 pp
|
|
2031
|
612
|
18
|
16
|
2
|
11.6%
|
0.344 pp
|
|
2032
|
626
|
19
|
16
|
3
|
13.6%
|
0.403 pp
|
|
2033
|
642
|
19
|
16
|
3
|
15.8%
|
0.466 pp
|
|
Oslo - Oslove
|
2029
|
95,450
|
4,217
|
4,121
|
95
|
2.3%
|
0.1 pp
|
|
2030
|
98,133
|
4,335
|
4,121
|
214
|
4.9%
|
0.218 pp
|
|
2031
|
100,949
|
4,459
|
4,121
|
338
|
7.6%
|
0.335 pp
|
|
2032
|
103,916
|
4,591
|
4,121
|
469
|
10.2%
|
0.452 pp
|
|
2033
|
106,889
|
4,722
|
4,121
|
601
|
12.7%
|
0.562 pp
|
|
2034
|
109,868
|
4,853
|
4,121
|
732
|
15.1%
|
0.667 pp
|
|
Nesna
|
2029
|
406
|
15
|
15
|
0
|
2.6%
|
0.099 pp
|
|
2030
|
415
|
16
|
15
|
1
|
4.7%
|
0.18 pp
|
|
2031
|
424
|
16
|
15
|
1
|
6.8%
|
0.256 pp
|
|
2032
|
434
|
16
|
15
|
1
|
8.9%
|
0.338 pp
|
|
2033
|
444
|
17
|
15
|
2
|
11%
|
0.416 pp
|
|
2034
|
452
|
17
|
15
|
2
|
12.5%
|
0.475 pp
|
|
Giske
|
2030
|
1,625
|
24
|
24
|
0
|
0.8%
|
0.012 pp
|
|
2031
|
1,661
|
25
|
24
|
1
|
3%
|
0.044 pp
|
|
2032
|
1,699
|
25
|
24
|
1
|
5.2%
|
0.077 pp
|
|
2033
|
1,729
|
26
|
24
|
2
|
6.8%
|
0.101 pp
|
|
2034
|
1,760
|
26
|
24
|
2
|
8.4%
|
0.126 pp
|
|
L?renskog
|
2030
|
8,071
|
244
|
240
|
4
|
1.7%
|
0.052 pp
|
|
2031
|
8,350
|
253
|
240
|
13
|
5%
|
0.152 pp
|
|
2032
|
8,670
|
262
|
240
|
22
|
8.5%
|
0.258 pp
|
|
2033
|
8,964
|
271
|
240
|
31
|
11.5%
|
0.348 pp
|
|
2034
|
9,269
|
280
|
240
|
40
|
14.4%
|
0.437 pp
|
|
Rollag
|
2030
|
408
|
16
|
16
|
0
|
0.9%
|
0.036 pp
|
|
2031
|
415
|
16
|
16
|
0
|
2.6%
|
0.102 pp
|
|
2032
|
419
|
17
|
16
|
1
|
3.5%
|
0.139 pp
|
|
2033
|
429
|
17
|
16
|
1
|
5.8%
|
0.228 pp
|
|
2034
|
429
|
17
|
16
|
1
|
5.8%
|
0.228 pp
|
|
Nordreisa - R?isa - Raisi
|
2030
|
1,092
|
52
|
52
|
0
|
0.1%
|
0.005 pp
|
|
2031
|
1,097
|
52
|
52
|
0
|
0.6%
|
0.026 pp
|
|
2032
|
1,130
|
54
|
52
|
2
|
3.5%
|
0.165 pp
|
|
2033
|
1,148
|
55
|
52
|
3
|
5%
|
0.237 pp
|
|
2034
|
1,161
|
55
|
52
|
3
|
6%
|
0.288 pp
|
|
Vindafjord
|
2031
|
1,855
|
81
|
81
|
0
|
0.4%
|
0.019 pp
|
|
2032
|
1,904
|
84
|
81
|
3
|
3%
|
0.132 pp
|
|
2033
|
1,955
|
86
|
81
|
5
|
5.5%
|
0.243 pp
|
|
2034
|
1,975
|
87
|
81
|
6
|
6.5%
|
0.285 pp
|
|
Gildesk?l
|
2031
|
591
|
32
|
32
|
0
|
0.3%
|
0.016 pp
|
|
2032
|
600
|
33
|
32
|
1
|
1.8%
|
0.097 pp
|
|
2033
|
605
|
33
|
32
|
1
|
2.6%
|
0.141 pp
|
|
2034
|
615
|
33
|
32
|
1
|
4.2%
|
0.227 pp
|
|
R?de
|
2031
|
1,759
|
52
|
51
|
1
|
1.4%
|
0.04 pp
|
|
2032
|
1,793
|
53
|
51
|
2
|
3.2%
|
0.095 pp
|
|
2033
|
1,836
|
54
|
51
|
3
|
5.5%
|
0.162 pp
|
|
2034
|
1,878
|
55
|
51
|
4
|
7.6%
|
0.224 pp
|
|
Nannestad
|
2031
|
2,476
|
69
|
68
|
1
|
1.1%
|
0.03 pp
|
|
2032
|
2,596
|
72
|
68
|
4
|
5.7%
|
0.157 pp
|
|
2033
|
2,699
|
75
|
68
|
7
|
9.3%
|
0.257 pp
|
|
2034
|
2,820
|
78
|
68
|
10
|
13.1%
|
0.365 pp
|
|
Stord
|
2031
|
3,716
|
110
|
109
|
1
|
0.7%
|
0.02 pp
|
|
2032
|
3,810
|
113
|
109
|
4
|
3.1%
|
0.092 pp
|
|
2033
|
3,914
|
116
|
109
|
7
|
5.7%
|
0.168 pp
|
|
2034
|
4,044
|
119
|
109
|
10
|
8.7%
|
0.258 pp
|
|
Skaun
|
2031
|
1,490
|
50
|
50
|
0
|
0.1%
|
0.004 pp
|
|
2032
|
1,542
|
52
|
50
|
2
|
3.5%
|
0.117 pp
|
|
2033
|
1,599
|
54
|
50
|
4
|
6.9%
|
0.233 pp
|
|
2034
|
1,657
|
56
|
50
|
6
|
10.2%
|
0.342 pp
|
|
Karls?y
|
2031
|
679
|
42
|
42
|
0
|
0.4%
|
0.024 pp
|
|
2032
|
697
|
43
|
42
|
1
|
3%
|
0.183 pp
|
|
2033
|
708
|
44
|
42
|
2
|
4.5%
|
0.277 pp
|
|
2034
|
715
|
44
|
42
|
2
|
5.4%
|
0.335 pp
|
|
Enebakk
|
2032
|
2,135
|
72
|
71
|
1
|
1.9%
|
0.063 pp
|
|
2033
|
2,207
|
75
|
71
|
4
|
5.1%
|
0.172 pp
|
|
2034
|
2,276
|
77
|
71
|
6
|
7.9%
|
0.269 pp
|
|
Nittedal
|
2032
|
4,420
|
125
|
124
|
2
|
1.3%
|
0.037 pp
|
|
2033
|
4,572
|
130
|
124
|
6
|
4.6%
|
0.13 pp
|
|
2034
|
4,757
|
135
|
124
|
11
|
8.3%
|
0.235 pp
|
|
Bardu
|
2032
|
832
|
37
|
37
|
0
|
0%
|
0 pp
|
|
2033
|
849
|
38
|
37
|
1
|
2%
|
0.089 pp
|
|
2034
|
856
|
38
|
37
|
1
|
2.8%
|
0.125 pp
|
|
B?tsfjord
|
2032
|
428
|
17
|
16
|
1
|
3.6%
|
0.14 pp
|
|
2033
|
429
|
17
|
16
|
1
|
3.8%
|
0.149 pp
|
|
2034
|
439
|
17
|
16
|
1
|
6%
|
0.233 pp
|
|
Stavanger
|
2033
|
26,841
|
1,136
|
1,119
|
17
|
1.5%
|
0.062 pp
|
|
2034
|
27,585
|
1,167
|
1,119
|
48
|
4.1%
|
0.175 pp
|
|
Sandnes
|
2033
|
14,180
|
523
|
518
|
5
|
0.9%
|
0.034 pp
|
|
2034
|
14,620
|
539
|
518
|
21
|
3.9%
|
0.144 pp
|
|
Sola
|
2033
|
4,904
|
188
|
187
|
1
|
0.4%
|
0.014 pp
|
|
2034
|
5,062
|
194
|
187
|
7
|
3.5%
|
0.133 pp
|
|
Lillestr?m
|
2033
|
16,892
|
490
|
473
|
17
|
3.4%
|
0.099 pp
|
|
2034
|
17,540
|
509
|
473
|
36
|
7%
|
0.202 pp
|
|
R?lingen
|
2033
|
3,340
|
89
|
86
|
3
|
3.1%
|
0.081 pp
|
|
2034
|
3,456
|
92
|
86
|
6
|
6.3%
|
0.168 pp
|
|
Nes
|
2033
|
5,067
|
129
|
125
|
4
|
2.8%
|
0.071 pp
|
|
2034
|
5,244
|
133
|
125
|
8
|
6.1%
|
0.154 pp
|
|
Hol
|
2033
|
1,199
|
39
|
39
|
0
|
0.5%
|
0.018 pp
|
|
2034
|
1,221
|
40
|
39
|
1
|
2.3%
|
0.076 pp
|
|
Gausdal
|
2033
|
1,599
|
63
|
63
|
0
|
0.8%
|
0.031 pp
|
|
2034
|
1,633
|
65
|
63
|
2
|
2.8%
|
0.113 pp
|
|
Trondheim - Tr?ante
|
2033
|
37,366
|
1,545
|
1,542
|
3
|
0.2%
|
0.009 pp
|
|
2034
|
38,378
|
1,587
|
1,542
|
45
|
2.8%
|
0.118 pp
|
|
Sunndal
|
2034
|
1,851
|
78
|
78
|
0
|
0.2%
|
0.008 pp
|
|
Alstahaug
|
2034
|
1,839
|
57
|
56
|
1
|
1.5%
|
0.048 pp
|
|
Lur?y
|
2034
|
502
|
13
|
13
|
0
|
0.1%
|
0.002 pp
|
|
Halden
|
2034
|
7,807
|
219
|
218
|
1
|
0.5%
|
0.013 pp
|
|
Kristiansand
|
2034
|
22,932
|
735
|
726
|
9
|
1.2%
|
0.04 pp
|
|
Vennesla
|
2034
|
3,043
|
101
|
100
|
2
|
1.6%
|
0.052 pp
|
|
Bergen
|
2034
|
54,563
|
2,615
|
2,609
|
6
|
0.2%
|
0.011 pp
|
|
Troms?
|
2034
|
14,272
|
478
|
474
|
4
|
0.8%
|
0.027 pp
|
# ── Step 5: Trajectory plot for selected municipalities ──────────────────────
# Show demand vs supply over time for the 6 most urgent municipalities
top6_munis <- muni_crossover |>
filter(scenario == "Strong ageing", supply_scenario == "flat") |>
arrange(crossover_year) |>
slice_head(n = 10) |>
pull(muni_code)
muni_backcalib |>
filter(muni_code %in% top6_munis, scenario == "Strong ageing") |>
left_join(muni_lookup, by = "muni_code") |>
ggplot(aes(x = year)) +
geom_line(aes(y = projected_demand, colour = "Projected demand"),
linewidth = 1.1) +
geom_line(aes(y = supply, colour = "Supply (flat)"),
linewidth = 1.1, linetype = "dashed") +
geom_ribbon(aes(ymin = supply, ymax = pmax(projected_demand, supply),
fill = projected_demand > supply),
alpha = 0.15) +
scale_colour_manual(values = c("Projected demand" = "#C00000",
"Supply (flat)" = "#4472C4")) +
scale_fill_manual(values = c("FALSE" = "#4472C4", "TRUE" = "#C00000"),
guide = "none") +
scale_x_continuous(breaks = seq(2025, 2036, 2)) +
facet_wrap(~ muni_name.y, scales = "free_y", ncol = 3) +
labs(
title = "Demand vs supply trajectory — 6 most urgent municipalities",
subtitle = "Strong ageing scenario, flat supply. Red shading = demand exceeds supply.",
x = NULL, y = "Long-term institutional users / beds",
colour = NULL
) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank(),
strip.text = element_text(face = "bold"))
Markov Model
!!!!!!!!
# ════════════════════════════════════════════════════════════════════════════
# PART A — MARKOV MODEL ENGINE
# ════════════════════════════════════════════════════════════════════════════
# ── A1: Define the transition matrix ─────────────────────────────────────────
# Rows = state you are IN, Columns = state you move TO
# Each row must sum to 1.0 (you must go somewhere each year)
#
# States:
# 1 = Home care light
# 2 = Home care moderate
# 3 = Home care intensive
# 4 = Nursing home
# 5 = Dead
build_transition_matrix <- function(
# Progression probabilities (moving to next home care level)
p_light_to_moderate = 0.15,
p_moderate_to_intense = 0.12,
# Nursing home admission probabilities (from each home care level)
p_light_to_nh = 0.02,
p_moderate_to_nh = 0.08,
p_intense_to_nh = 0.20,
# Annual mortality probabilities (by state)
p_die_light = 0.03,
p_die_moderate = 0.05,
p_die_intense = 0.08,
p_die_nh = 0.25
) {
# Build the 5x5 matrix row by row
# Each row: probabilities of moving to states 1, 2, 3, 4, 5
m <- matrix(0, nrow = 5, ncol = 5,
dimnames = list(
from = c("Light", "Moderate", "Intensive", "NH", "Dead"),
to = c("Light", "Moderate", "Intensive", "NH", "Dead")
))
# Row 1: Light home care
m[1, 2] <- p_light_to_moderate # progress to moderate
m[1, 4] <- p_light_to_nh # admit to NH directly
m[1, 5] <- p_die_light # die
m[1, 1] <- 1 - sum(m[1, ]) # stay (remainder)
# Row 2: Moderate home care
m[2, 3] <- p_moderate_to_intense # progress to intensive
m[2, 4] <- p_moderate_to_nh # admit to NH
m[2, 5] <- p_die_moderate # die
m[2, 2] <- 1 - sum(m[2, ]) # stay
# Row 3: Intensive home care
m[3, 4] <- p_intense_to_nh # admit to NH
m[3, 5] <- p_die_intense # die
m[3, 3] <- 1 - sum(m[3, ]) # stay
# Row 4: Nursing home
m[4, 5] <- p_die_nh # die
m[4, 4] <- 1 - sum(m[4, ]) # stay
# Row 5: Dead (absorbing — once dead, stay dead)
m[5, 5] <- 1
# Validate: every row must sum to exactly 1
row_sums <- rowSums(m)
if (any(abs(row_sums - 1) > 1e-9)) {
stop("Transition matrix rows do not sum to 1. Check your probabilities.\n",
"Row sums: ", paste(round(row_sums, 6), collapse = ", "))
}
# Validate: no negative probabilities
if (any(m < 0)) {
stop("Negative probability detected. Your input probabilities sum to more than 1 ",
"in at least one row. Reduce them.")
}
return(m)
}
# ── A2: The model engine ──────────────────────────────────────────────────────
# Takes a starting cohort and runs it forward year by year
run_markov <- function(
transition_matrix, # 5x5 matrix from build_transition_matrix()
cohort_start, # named vector: how many people start in each state
# e.g. c(Light=1000, Moderate=500, Intensive=200, NH=300, Dead=0)
n_years = 11, # how many years to run (2025 to 2036 = 11 steps)
start_year = 2025
) {
# Storage: one row per year, one column per state
n_states <- nrow(transition_matrix)
state_names <- rownames(transition_matrix)
years <- start_year:(start_year + n_years - 1)
results <- matrix(
NA,
nrow = n_years,
ncol = n_states,
dimnames = list(year = years, state = state_names)
)
# Year 1 = starting cohort
results[1, ] <- cohort_start
# Run forward: multiply current state vector by transition matrix
for (t in 2:n_years) {
results[t, ] <- results[t - 1, ] %*% transition_matrix
}
# Convert to tidy data frame for easy plotting and joining
results_df <- as.data.frame(results) |>
tibble::rownames_to_column("year") |>
dplyr::mutate(year = as.integer(year)) |>
tidyr::pivot_longer(
cols = -year,
names_to = "state",
values_to = "n_people"
) |>
dplyr::mutate(
state = factor(state,
levels = c("Light", "Moderate", "Intensive", "NH", "Dead"))
)
return(results_df)
}
# ── A3: Helper — extract annual NH admissions ────────────────────────────────
# NH admissions in year t = people who were NOT in NH in year t-1
# but ARE in NH in year t. This is different from total NH residents.
get_nh_admissions <- function(markov_results) {
markov_results |>
dplyr::filter(state == "NH") |>
dplyr::arrange(year) |>
dplyr::mutate(
nh_admissions = n_people - dplyr::lag(n_people, default = n_people[1])
) |>
dplyr::select(year, nh_residents = n_people, nh_admissions)
}
# ── A4: Quick test run ───────────────────────────────────────────────────────
# Build a transition matrix with default probabilities
baseline_matrix <- build_transition_matrix()
cat("Transition matrix (baseline):\n")
## Transition matrix (baseline):
print(round(baseline_matrix, 3))
## to
## from Light Moderate Intensive NH Dead
## Light 0.8 0.15 0.00 0.02 0.03
## Moderate 0.0 0.75 0.12 0.08 0.05
## Intensive 0.0 0.00 0.72 0.20 0.08
## NH 0.0 0.00 0.00 0.75 0.25
## Dead 0.0 0.00 0.00 0.00 1.00
cat("\nRow sums (all must be 1.0):\n")
##
## Row sums (all must be 1.0):
print(rowSums(baseline_matrix))
## Light Moderate Intensive NH Dead
## 1 1 1 1 1
# Define a starting cohort of 10,000 home care users
# Distribution matches the national need-level shares from your data:
# Light 37%, Moderate 30%, Intensive 25%, NH 8% (already there), Dead 0%
cohort <- c(Light = 3700, Moderate = 3000, Intensive = 2500,
NH = 800, Dead = 0)
cat("\nStarting cohort:\n")
##
## Starting cohort:
print(cohort)
## Light Moderate Intensive NH Dead
## 3700 3000 2500 800 0
cat("Total:", sum(cohort), "people\n\n")
## Total: 10000 people
# Run the model
baseline_run <- run_markov(
transition_matrix = baseline_matrix,
cohort_start = cohort,
n_years = 12,
start_year = 2025
)
# Show results
cat("State populations by year:\n")
## State populations by year:
baseline_run |>
tidyr::pivot_wider(names_from = state, values_from = n_people) |>
dplyr::mutate(across(where(is.numeric), ~ round(.x))) |>
kable(caption = "Markov model — cohort state populations 2025-2036 (baseline)") |>
kable_styling(bootstrap_options = c("striped", "hover"))
Markov model — cohort state populations 2025-2036 (baseline)
|
year
|
Light
|
Moderate
|
Intensive
|
NH
|
Dead
|
|
2025
|
3700
|
3000
|
2500
|
800
|
0
|
|
2026
|
2960
|
2805
|
2160
|
1414
|
661
|
|
2027
|
2368
|
2548
|
1892
|
1776
|
1416
|
|
2028
|
1894
|
2266
|
1668
|
1962
|
2210
|
|
2029
|
1516
|
1984
|
1473
|
2024
|
3004
|
|
2030
|
1212
|
1715
|
1298
|
2002
|
3773
|
|
2031
|
970
|
1468
|
1141
|
1922
|
4499
|
|
2032
|
776
|
1247
|
997
|
1807
|
5173
|
|
2033
|
621
|
1051
|
868
|
1670
|
5790
|
|
2034
|
497
|
882
|
751
|
1522
|
6348
|
|
2035
|
397
|
736
|
646
|
1372
|
6848
|
|
2036
|
318
|
611
|
554
|
1225
|
7292
|
# NH admissions
cat("\nNursing home admissions per year:\n")
##
## Nursing home admissions per year:
get_nh_admissions(baseline_run) |>
mutate(across(where(is.numeric), ~ comma(round(.x)))) |>
kable(caption = "Annual NH admissions and total NH residents (baseline cohort)",
col.names = c("Year", "NH residents", "New admissions")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Annual NH admissions and total NH residents (baseline cohort)
|
Year
|
NH residents
|
New admissions
|
|
2,025
|
800
|
0
|
|
2,026
|
1,414
|
614
|
|
2,027
|
1,776
|
362
|
|
2,028
|
1,962
|
186
|
|
2,029
|
2,024
|
62
|
|
2,030
|
2,002
|
-22
|
|
2,031
|
1,922
|
-79
|
|
2,032
|
1,807
|
-116
|
|
2,033
|
1,670
|
-137
|
|
2,034
|
1,522
|
-147
|
|
2,035
|
1,372
|
-150
|
|
2,036
|
1,225
|
-147
|
# Plot
baseline_run |>
filter(state != "Dead") |>
ggplot(aes(x = year, y = n_people, fill = state)) +
geom_area(alpha = 0.8, colour = "white", linewidth = 0.3) +
scale_fill_manual(values = c(
"Light" = "#4472C4",
"Moderate" = "#ED7D31",
"Intensive" = "#FFC000",
"NH" = "#C00000"
)) +
scale_x_continuous(breaks = seq(2025, 2036, 2)) +
scale_y_continuous(labels = comma) +
labs(
title = "Cohort flow through care states 2025-2036",
subtitle = "Baseline transition probabilities — cohort of 10,000 home care users",
x = NULL, y = "Number of people", fill = "State"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank())
Pat B: Using real life
data:
# ════════════════════════════════════════════════════════════════════════════
# PART B — CALIBRATION + OPEN COHORT ENGINE
# ════════════════════════════════════════════════════════════════════════════
# ── B1: Calibration targets (unchanged from before) ──────────────────────────
total_hc_elderly <- sum(users$home_care_elderly, na.rm = TRUE)
observed_light <- round(total_hc_elderly * 0.367)
observed_moderate <- round(total_hc_elderly * 0.298)
observed_intensive <- round(total_hc_elderly * 0.247)
observed_nh <- sum(users$long_term_stay_in_institutions, na.rm = TRUE)
avg_nh_stay_years <- 3
observed_nh_admissions_annual <- round(observed_nh / avg_nh_stay_years)
p_die_light_obs <- 0.030
p_die_moderate_obs <- 0.050
p_die_intensive_obs <- 0.080
p_die_nh_obs <- 0.250
cat("=== 2025 CALIBRATION TARGETS ===\n")
## === 2025 CALIBRATION TARGETS ===
cat("Home care — light: ", comma(observed_light), "\n")
## Home care — light: 41,597
cat("Home care — moderate: ", comma(observed_moderate), "\n")
## Home care — moderate: 33,776
cat("Home care — intensive: ", comma(observed_intensive), "\n")
## Home care — intensive: 27,996
cat("NH long-term residents: ", comma(observed_nh), "\n")
## NH long-term residents: 30,622
cat("Estimated annual NH admissions:", comma(observed_nh_admissions_annual), "\n")
## Estimated annual NH admissions: 10,207
# ── B2: Calibration loss function (unchanged) ────────────────────────────────
calibration_loss <- function(params) {
p_l2m <- params[1]; p_m2i <- params[2]
p_l2nh <- params[3]; p_m2nh <- params[4]; p_i2nh <- params[5]
if (any(params < 0) || any(params > 1)) return(1e9)
if (p_l2m + p_l2nh + p_die_light_obs > 1) return(1e9)
if (p_m2i + p_m2nh + p_die_moderate_obs > 1) return(1e9)
if (p_i2nh + p_die_intensive_obs > 1) return(1e9)
m <- tryCatch(
build_transition_matrix(
p_light_to_moderate = p_l2m,
p_moderate_to_intense = p_m2i,
p_light_to_nh = p_l2nh,
p_moderate_to_nh = p_m2nh,
p_intense_to_nh = p_i2nh,
p_die_light = p_die_light_obs,
p_die_moderate = p_die_moderate_obs,
p_die_intense = p_die_intensive_obs,
p_die_nh = p_die_nh_obs
),
error = function(e) NULL
)
if (is.null(m)) return(1e9)
cohort_2025 <- c(Light = observed_light, Moderate = observed_moderate,
Intensive = observed_intensive, NH = observed_nh, Dead = 0)
result <- run_markov(m, cohort_2025, n_years = 2, start_year = 2025)
year2 <- result |> filter(year == 2026) |>
select(state, n_people) |> tibble::deframe()
nh_deaths_pred <- observed_nh * p_die_nh_obs
nh_admissions_pred <- (year2["NH"] - observed_nh) + nh_deaths_pred
loss <-
((year2["Light"] - observed_light) / observed_light)^2 * 1 +
((year2["Moderate"] - observed_moderate) / observed_moderate)^2 * 1 +
((year2["Intensive"]- observed_intensive)/ observed_intensive)^2* 1 +
((nh_admissions_pred - observed_nh_admissions_annual) /
observed_nh_admissions_annual)^2 * 3
return(loss)
}
# ── B3: Run optimiser and build calibrated matrix ────────────────────────────
cat("\nRunning calibration...\n")
##
## Running calibration...
start_params <- c(p_light_to_moderate = 0.15,
p_moderate_to_intense = 0.12,
p_light_to_nh = 0.02,
p_moderate_to_nh = 0.08,
p_intense_to_nh = 0.20)
calibration_result <- optim(
par = start_params,
fn = calibration_loss,
method = "Nelder-Mead",
control = list(maxit = 10000, reltol = 1e-10)
)
calibrated_params <- calibration_result$par
names(calibrated_params) <- names(start_params)
cat("Converged:", calibration_result$convergence == 0,
"| Loss:", round(calibration_result$value, 8), "\n\n")
## Converged: TRUE | Loss: 0.06413883
calibrated_matrix <- build_transition_matrix(
p_light_to_moderate = calibrated_params["p_light_to_moderate"],
p_moderate_to_intense = calibrated_params["p_moderate_to_intense"],
p_light_to_nh = calibrated_params["p_light_to_nh"],
p_moderate_to_nh = calibrated_params["p_moderate_to_nh"],
p_intense_to_nh = calibrated_params["p_intense_to_nh"],
p_die_light = p_die_light_obs,
p_die_moderate = p_die_moderate_obs,
p_die_intense = p_die_intensive_obs,
p_die_nh = p_die_nh_obs
)
cat("Calibrated transition matrix:\n")
## Calibrated transition matrix:
print(round(calibrated_matrix, 4))
## to
## from Light Moderate Intensive NH Dead
## Light 0.8266 0.1154 0.0000 0.0281 0.03
## Moderate 0.0000 0.7171 0.1318 0.1011 0.05
## Intensive 0.0000 0.0000 0.7242 0.1958 0.08
## NH 0.0000 0.0000 0.0000 0.7500 0.25
## Dead 0.0000 0.0000 0.0000 0.0000 1.00
# ── B4: New entrant projection ───────────────────────────────────────────────
# New entrants = elderly people entering home care for the first time each year.
# Estimated as: growth in elderly population × home care entry rate.
#
# Home care entry rate = share of elderly population entering HC in a year.
# We estimate this from the current stock and average HC duration.
# Average time in home care before NH or death ≈ 4 years (literature estimate).
# Entry rate ≈ current HC users / (elderly population × avg duration)
avg_hc_duration_years <- 4.0
# 2025 total elderly population (from projections, both scenarios identical in 2025)
elderly_pop_2025 <- national_proj |>
filter(year == 2025, scenario == "Strong ageing") |>
summarise(total = sum(population)) |>
pull(total)
# Annual home care entry rate
hc_entry_rate <- total_hc_elderly / (elderly_pop_2025 * avg_hc_duration_years)
cat("\nHome care entry rate (annual probability of entering HC):",
round(hc_entry_rate, 4), "\n")
##
## Home care entry rate (annual probability of entering HC): 0.0302
# Project new entrants for each year and scenario
new_entrants_proj <- national_elderly_proj |>
mutate(
new_entrants = round(total_elderly * hc_entry_rate),
# New entrants start in light home care state
entrants_light = new_entrants,
entrants_moderate = 0,
entrants_intensive = 0
)
cat("\nProjected new HC entrants by year (strong ageing):\n")
##
## Projected new HC entrants by year (strong ageing):
new_entrants_proj |>
filter(scenario == "Strong ageing") |>
select(year, total_elderly, new_entrants) |>
mutate(across(where(is.numeric), comma)) |>
kable(caption = "Projected annual new home care entrants",
col.names = c("Year", "Total elderly", "New HC entrants")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Projected annual new home care entrants
|
Year
|
Total elderly
|
New HC entrants
|
|
2,025
|
938,331
|
28,336
|
|
2,026
|
961,814
|
29,045
|
|
2,027
|
985,854
|
29,771
|
|
2,028
|
1,009,922
|
30,498
|
|
2,029
|
1,034,215
|
31,231
|
|
2,030
|
1,058,794
|
31,973
|
|
2,031
|
1,084,221
|
32,741
|
|
2,032
|
1,111,425
|
33,563
|
|
2,033
|
1,139,078
|
34,398
|
|
2,034
|
1,167,117
|
35,245
|
# ── B5: Open cohort Markov engine ────────────────────────────────────────────
# Runs the Markov model with new entrants added each year.
run_markov_open <- function(
transition_matrix,
cohort_start,
new_entrants_df,
scenario_name,
start_year = 2025
) {
years <- sort(unique(new_entrants_df$year))
n_years <- length(years)
state_names <- rownames(transition_matrix)
results <- matrix(NA, nrow = n_years, ncol = length(state_names),
dimnames = list(year = as.character(years),
state = state_names))
results[1, ] <- cohort_start
entrants <- new_entrants_df |>
filter(scenario == scenario_name) |>
arrange(year)
for (t in 2:n_years) {
# %*% returns a matrix — drop to a named vector with as.vector()
transitioned <- as.vector(results[t - 1, ] %*% transition_matrix)
names(transitioned) <- state_names
# Add new entrants for this year
yr <- years[t]
ent_row <- entrants |> filter(year == yr)
if (nrow(ent_row) > 0) {
transitioned["Light"] <- transitioned["Light"] + ent_row$entrants_light
transitioned["Moderate"] <- transitioned["Moderate"] + ent_row$entrants_moderate
transitioned["Intensive"] <- transitioned["Intensive"] + ent_row$entrants_intensive
}
results[t, ] <- transitioned
}
as.data.frame(results) |>
tibble::rownames_to_column("year") |>
mutate(year = as.integer(year), scenario = scenario_name) |>
pivot_longer(cols = -c(year, scenario),
names_to = "state",
values_to = "n_people") |>
mutate(state = factor(state,
levels = c("Light","Moderate","Intensive","NH","Dead")))
}
# ── B6: Run open cohort model for both scenarios ─────────────────────────────
cohort_2025 <- c(Light = observed_light,
Moderate = observed_moderate,
Intensive = observed_intensive,
NH = observed_nh,
Dead = 0)
open_strong <- run_markov_open(
transition_matrix = calibrated_matrix,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
)
open_weak <- run_markov_open(
transition_matrix = calibrated_matrix,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Weak ageing",
start_year = 2025
)
open_combined <- bind_rows(open_strong, open_weak)
# ── B7: NH trajectory with supply comparison ─────────────────────────────────
national_supply_flat <- sum(users$institution_all_available_beds, na.rm = TRUE)
nh_open_trajectory <- open_combined |>
filter(state == "NH") |>
arrange(scenario, year) |>
group_by(scenario) |>
mutate(
nh_admissions = n_people - lag(n_people, default = first(n_people)),
supply = national_supply_flat,
gap = n_people - supply,
over_capacity = gap > 0
) |>
ungroup()
cat("\nNH residents vs supply — open cohort model:\n")
##
## NH residents vs supply — open cohort model:
nh_open_trajectory |>
mutate(
n_people = comma(round(n_people)),
nh_admissions = comma(round(nh_admissions)),
supply = comma(supply),
gap = comma(round(gap))
) |>
select(scenario, year, n_people, nh_admissions, supply, gap, over_capacity) |>
kable(
caption = "Open cohort — NH residents vs supply 2025-2036",
col.names = c("Scenario", "Year", "NH residents", "Net change",
"Supply", "Gap", "Over capacity")
) |>
kable_styling(bootstrap_options = c("striped", "hover")) |>
row_spec(which(nh_open_trajectory$over_capacity), background = "#FCE4D6")
Open cohort — NH residents vs supply 2025-2036
|
Scenario
|
Year
|
NH residents
|
Net change
|
Supply
|
Gap
|
Over capacity
|
|
Strong ageing
|
2025
|
30,622
|
0
|
41,650
|
-11,028
|
FALSE
|
|
Strong ageing
|
2026
|
33,029
|
2,407
|
41,650
|
-8,622
|
FALSE
|
|
Strong ageing
|
2027
|
34,326
|
1,297
|
41,650
|
-7,324
|
FALSE
|
|
Strong ageing
|
2028
|
35,149
|
823
|
41,650
|
-6,501
|
FALSE
|
|
Strong ageing
|
2029
|
35,929
|
780
|
41,650
|
-5,721
|
FALSE
|
|
Strong ageing
|
2030
|
36,927
|
998
|
41,650
|
-4,723
|
FALSE
|
|
Strong ageing
|
2031
|
38,275
|
1,348
|
41,650
|
-3,376
|
FALSE
|
|
Strong ageing
|
2032
|
40,017
|
1,742
|
41,650
|
-1,633
|
FALSE
|
|
Strong ageing
|
2033
|
42,141
|
2,124
|
41,650
|
490
|
TRUE
|
|
Strong ageing
|
2034
|
44,599
|
2,459
|
41,650
|
2,949
|
TRUE
|
|
Weak ageing
|
2025
|
30,622
|
0
|
41,650
|
-11,028
|
FALSE
|
|
Weak ageing
|
2026
|
33,029
|
2,407
|
41,650
|
-8,622
|
FALSE
|
|
Weak ageing
|
2027
|
34,321
|
1,292
|
41,650
|
-7,329
|
FALSE
|
|
Weak ageing
|
2028
|
35,132
|
810
|
41,650
|
-6,519
|
FALSE
|
|
Weak ageing
|
2029
|
35,887
|
755
|
41,650
|
-5,763
|
FALSE
|
|
Weak ageing
|
2030
|
36,844
|
958
|
41,650
|
-4,806
|
FALSE
|
|
Weak ageing
|
2031
|
38,134
|
1,290
|
41,650
|
-3,516
|
FALSE
|
|
Weak ageing
|
2032
|
39,797
|
1,663
|
41,650
|
-1,854
|
FALSE
|
|
Weak ageing
|
2033
|
41,817
|
2,020
|
41,650
|
167
|
TRUE
|
|
Weak ageing
|
2034
|
44,147
|
2,330
|
41,650
|
2,497
|
TRUE
|
# ── B8: Plot ─────────────────────────────────────────────────────────────────
nh_open_trajectory |>
ggplot(aes(x = year)) +
geom_line(aes(y = n_people, colour = scenario), linewidth = 1.2) +
geom_hline(aes(yintercept = national_supply_flat),
linetype = "dashed", colour = "#4472C4", linewidth = 1) +
geom_ribbon(
data = nh_open_trajectory |> filter(over_capacity),
aes(ymin = supply, ymax = n_people, fill = scenario),
alpha = 0.2
) +
annotate("text", x = 2026, y = national_supply_flat + 500,
label = "Supply ceiling (flat)", hjust = 0,
colour = "#4472C4", size = 3.5) +
scale_colour_manual(values = c("Strong ageing" = "#C00000",
"Weak ageing" = "#ED7D31")) +
scale_fill_manual(values = c("Strong ageing" = "#C00000",
"Weak ageing" = "#ED7D31"),
guide = "none") +
scale_x_continuous(breaks = 2025:2036) +
scale_y_continuous(labels = comma) +
labs(
title = "Projected NH residents vs supply — open cohort Markov model",
subtitle = "New entrants added annually from population projections. Dashed = flat supply.",
x = NULL, y = "NH residents / beds", colour = "Ageing scenario"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank())
Part C: What if we
stabalise?
# ════════════════════════════════════════════════════════════════════════════
# PART C — SCENARIO COMPARISON
# Effect of early-stage stabilisation on nursing home demand
# ════════════════════════════════════════════════════════════════════════════
# ── C1: Intervention logic ───────────────────────────────────────────────────
# The intervention reduces progression probabilities at light and moderate stages.
# This slows the pipeline feeding into intensive care and eventually into NH.
# Effect is parameterised as a percentage reduction in progression probabilities.
# A secondary smaller effect is allowed at moderate-to-intensive and
# intensive-to-NH to reflect partial downstream benefit.
#
# intervention_strength: 0 = no effect, 1 = 100% reduction (impossible in practice)
# Realistic range: 0.10 to 0.40 (10% to 40% reduction in progression rates)
build_intervention_matrix <- function(
base_matrix,
intervention_strength, # primary effect at light stage
secondary_multiplier = 0.4 # secondary effect = 40% of primary, at moderate stage
) {
# Extract baseline probabilities from the calibrated matrix
p_l2m_base <- base_matrix["Light", "Moderate"]
p_m2i_base <- base_matrix["Moderate", "Intensive"]
p_l2nh_base <- base_matrix["Light", "NH"]
p_m2nh_base <- base_matrix["Moderate", "NH"]
p_i2nh_base <- base_matrix["Intensive","NH"]
# Apply reductions
# Primary effect: slow light → moderate progression
p_l2m_new <- p_l2m_base * (1 - intervention_strength)
# Secondary effect: slow moderate → intensive progression
p_m2i_new <- p_m2i_base * (1 - intervention_strength * secondary_multiplier)
# Small tertiary effect at intensive → NH (40% of secondary = 16% of primary)
p_i2nh_new <- p_i2nh_base * (1 - intervention_strength *
secondary_multiplier * secondary_multiplier)
# NH admission from light and moderate: minor effect from stabilisation
p_l2nh_new <- p_l2nh_base * (1 - intervention_strength * 0.2)
p_m2nh_new <- p_m2nh_base * (1 - intervention_strength * secondary_multiplier * 0.5)
# Build new matrix with intervention probabilities
build_transition_matrix(
p_light_to_moderate = p_l2m_new,
p_moderate_to_intense = p_m2i_new,
p_light_to_nh = p_l2nh_new,
p_moderate_to_nh = p_m2nh_new,
p_intense_to_nh = p_i2nh_new,
p_die_light = p_die_light_obs,
p_die_moderate = p_die_moderate_obs,
p_die_intense = p_die_intensive_obs,
p_die_nh = p_die_nh_obs
)
}
# ── C2: Show what the intervention does to transition probabilities ───────────
cat("=== EFFECT OF INTERVENTION ON TRANSITION PROBABILITIES ===\n\n")
## === EFFECT OF INTERVENTION ON TRANSITION PROBABILITIES ===
intervention_levels <- c(0, 0.10, 0.20, 0.30, 0.40)
prob_comparison <- map_dfr(intervention_levels, function(strength) {
m <- build_intervention_matrix(calibrated_matrix, strength)
tibble(
intervention_strength = paste0(strength * 100, "%"),
p_light_to_moderate = round(m["Light", "Moderate"], 4),
p_moderate_to_intensive = round(m["Moderate", "Intensive"], 4),
p_intensive_to_nh = round(m["Intensive","NH"], 4),
p_light_to_nh = round(m["Light", "NH"], 4),
p_moderate_to_nh = round(m["Moderate", "NH"], 4)
)
})
prob_comparison |>
kable(
caption = "Transition probabilities under each intervention strength",
col.names = c("Intervention", "Light→Moderate", "Moderate→Intensive",
"Intensive→NH", "Light→NH", "Moderate→NH")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Transition probabilities under each intervention strength
|
Intervention
|
Light→Moderate
|
Moderate→Intensive
|
Intensive→NH
|
Light→NH
|
Moderate→NH
|
|
0%
|
0.1154
|
0.1318
|
0.1958
|
0.0281
|
0.1011
|
|
10%
|
0.1038
|
0.1266
|
0.1927
|
0.0275
|
0.0990
|
|
20%
|
0.0923
|
0.1213
|
0.1895
|
0.0269
|
0.0970
|
|
30%
|
0.0808
|
0.1160
|
0.1864
|
0.0264
|
0.0950
|
|
40%
|
0.0692
|
0.1108
|
0.1833
|
0.0258
|
0.0930
|
# ── C3: Run open cohort model for each intervention level ────────────────────
cat("\nRunning scenarios...\n")
##
## Running scenarios...
scenarios_all <- map_dfr(intervention_levels, function(strength) {
int_matrix <- build_intervention_matrix(calibrated_matrix, strength)
run_markov_open(
transition_matrix = int_matrix,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
) |>
mutate(intervention = paste0(strength * 100, "% reduction"))
})
cat("Done.\n")
## Done.
# ── C4: Extract NH trajectory for each scenario ──────────────────────────────
nh_scenarios <- scenarios_all |>
filter(state == "NH") |>
arrange(intervention, year) |>
group_by(intervention) |>
mutate(
nh_admissions = n_people - lag(n_people, default = first(n_people)),
supply = national_supply_flat,
gap = n_people - supply,
over_capacity = gap > 0
) |>
ungroup() |>
mutate(
intervention = factor(intervention,
levels = paste0(intervention_levels * 100, "% reduction"))
)
# ── C5: Crossover year by intervention strength ──────────────────────────────
crossover_by_intervention <- nh_scenarios |>
filter(over_capacity) |>
group_by(intervention) |>
summarise(crossover_year = min(year), .groups = "drop")
# Add scenarios with no crossover
all_interventions <- tibble(
intervention = factor(paste0(intervention_levels * 100, "% reduction"),
levels = paste0(intervention_levels * 100, "% reduction"))
)
crossover_by_intervention <- all_interventions |>
left_join(crossover_by_intervention, by = "intervention") |>
mutate(crossover_year = if_else(is.na(crossover_year),
"No crossover by 2036", as.character(crossover_year)))
cat("\n=== CROSSOVER YEAR BY INTERVENTION STRENGTH ===\n")
##
## === CROSSOVER YEAR BY INTERVENTION STRENGTH ===
crossover_by_intervention |>
kable(
caption = "First year NH demand exceeds supply — by intervention strength",
col.names = c("Intervention strength", "Crossover year")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
First year NH demand exceeds supply — by intervention strength
|
Intervention strength
|
Crossover year
|
|
0% reduction
|
2033
|
|
10% reduction
|
2034
|
|
20% reduction
|
No crossover by 2036
|
|
30% reduction
|
No crossover by 2036
|
|
40% reduction
|
No crossover by 2036
|
# ── C6: People kept out of nursing homes ─────────────────────────────────────
# Compared to baseline (0% intervention)
baseline_nh <- nh_scenarios |>
filter(intervention == "0% reduction") |>
select(year, baseline_residents = n_people)
avoided_admissions <- nh_scenarios |>
left_join(baseline_nh, by = "year") |>
mutate(
residents_avoided = baseline_residents - n_people,
cumulative_avoided = cumsum(pmax(residents_avoided, 0))
) |>
filter(intervention != "0% reduction") |>
select(intervention, year, n_people, baseline_residents,
residents_avoided, cumulative_avoided, over_capacity)
cat("\nPeople kept out of nursing homes vs baseline:\n")
##
## People kept out of nursing homes vs baseline:
avoided_admissions |>
filter(year >= 2029) |>
mutate(
n_people = comma(round(n_people)),
baseline_residents = comma(round(baseline_residents)),
residents_avoided = comma(round(residents_avoided)),
cumulative_avoided = comma(round(cumulative_avoided))
) |>
select(intervention, year, baseline_residents,
n_people, residents_avoided, cumulative_avoided) |>
kable(
caption = "NH residents avoided vs baseline — from crossover year onwards",
col.names = c("Intervention", "Year", "Baseline NH residents",
"Intervention NH residents", "Avoided per year",
"Cumulative avoided")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), font_size = 11)
NH residents avoided vs baseline — from crossover year onwards
|
Intervention
|
Year
|
Baseline NH residents
|
Intervention NH residents
|
Avoided per year
|
Cumulative avoided
|
|
10% reduction
|
2029
|
35,929
|
35,283
|
646
|
1,643
|
|
10% reduction
|
2030
|
36,927
|
36,102
|
824
|
2,467
|
|
10% reduction
|
2031
|
38,275
|
37,259
|
1,016
|
3,483
|
|
10% reduction
|
2032
|
40,017
|
38,801
|
1,216
|
4,699
|
|
10% reduction
|
2033
|
42,141
|
40,724
|
1,416
|
6,116
|
|
10% reduction
|
2034
|
44,599
|
42,989
|
1,611
|
7,726
|
|
20% reduction
|
2029
|
35,929
|
34,640
|
1,290
|
11,005
|
|
20% reduction
|
2030
|
36,927
|
35,281
|
1,646
|
12,651
|
|
20% reduction
|
2031
|
38,275
|
36,242
|
2,032
|
14,684
|
|
20% reduction
|
2032
|
40,017
|
37,580
|
2,437
|
17,121
|
|
20% reduction
|
2033
|
42,141
|
39,293
|
2,847
|
19,968
|
|
20% reduction
|
2034
|
44,599
|
41,352
|
3,247
|
23,216
|
|
30% reduction
|
2029
|
35,929
|
34,000
|
1,930
|
28,125
|
|
30% reduction
|
2030
|
36,927
|
34,462
|
2,465
|
30,589
|
|
30% reduction
|
2031
|
38,275
|
35,227
|
3,048
|
33,637
|
|
30% reduction
|
2032
|
40,017
|
36,353
|
3,663
|
37,301
|
|
30% reduction
|
2033
|
42,141
|
37,850
|
4,291
|
41,592
|
|
30% reduction
|
2034
|
44,599
|
39,690
|
4,909
|
46,501
|
|
40% reduction
|
2029
|
35,929
|
33,363
|
2,566
|
53,033
|
|
40% reduction
|
2030
|
36,927
|
33,647
|
3,280
|
56,313
|
|
40% reduction
|
2031
|
38,275
|
34,213
|
4,062
|
60,375
|
|
40% reduction
|
2032
|
40,017
|
35,124
|
4,893
|
65,268
|
|
40% reduction
|
2033
|
42,141
|
36,394
|
5,746
|
71,014
|
|
40% reduction
|
2034
|
44,599
|
38,005
|
6,594
|
77,608
|
# ── C7: Main trajectory plot ─────────────────────────────────────────────────
nh_scenarios |>
ggplot(aes(x = year, y = n_people,
colour = intervention, group = intervention)) +
geom_line(linewidth = 1.1) +
geom_hline(yintercept = national_supply_flat,
linetype = "dashed", colour = "black", linewidth = 0.8) +
annotate("text", x = 2025.2, y = national_supply_flat + 600,
label = "Supply ceiling", hjust = 0, size = 3.5, colour = "black") +
scale_colour_manual(values = c(
"0% reduction" = "#C00000",
"10% reduction" = "#ED7D31",
"20% reduction" = "#FFC000",
"30% reduction" = "#70AD47",
"40% reduction" = "#4472C4"
)) +
scale_x_continuous(breaks = 2025:2036) +
scale_y_continuous(labels = comma,
limits = c(NA, max(nh_scenarios$n_people) * 1.05)) +
labs(
title = "NH residents under different intervention strengths",
subtitle = "Early-stage stabilisation in home care — strong ageing scenario",
x = NULL, y = "NH residents", colour = "Intervention"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank())
# ── C8: Crossover postponement bar chart ─────────────────────────────────────
crossover_plot_data <- crossover_by_intervention |>
mutate(
crossover_numeric = suppressWarnings(as.integer(crossover_year)),
label = if_else(
is.na(crossover_numeric),
"No crossover\nby 2036",
as.character(crossover_numeric)
),
# For plotting: use 2037 as a placeholder for "no crossover"
plot_year = if_else(is.na(crossover_numeric), 2037L, crossover_numeric),
baseline_year = min(crossover_numeric, na.rm = TRUE),
years_postponed = plot_year - baseline_year,
postponed_label = case_when(
years_postponed == 0 ~ "Baseline",
plot_year == 2037 ~ "Prevented",
TRUE ~ paste0("+", years_postponed, " yrs")
),
full_label = paste0(label, "\n(", postponed_label, ")")
)
crossover_plot_data |>
ggplot(aes(x = intervention, y = plot_year, fill = intervention)) +
geom_col(width = 0.6, alpha = 0.85) +
geom_text(aes(label = full_label),
vjust = -0.3, size = 3.5, fontface = "bold") +
geom_hline(yintercept = 2036, linetype = "dashed",
colour = "black", linewidth = 0.7) +
annotate("text", x = 0.5, y = 2036.3,
label = "End of projection window (2036)",
hjust = 0, size = 3, colour = "black") +
scale_fill_manual(values = c(
"0% reduction" = "#C00000",
"10% reduction" = "#ED7D31",
"20% reduction" = "#FFC000",
"30% reduction" = "#70AD47",
"40% reduction" = "#4472C4"
)) +
scale_y_continuous(limits = c(2020, 2040),
breaks = seq(2020, 2040, 2)) +
labs(
title = "Crossover year postponed by early-stage intervention",
subtitle = "Bars reaching 2037 indicate no crossover within the projection window",
x = "Intervention strength", y = "Crossover year"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
panel.grid.minor = element_blank())
# ════════════════════════════════════════════════════════════════════════════
# C9: BACK-CALIBRATION — fine-grained intervention search
# ════════════════════════════════════════════════════════════════════════════
# ── Step 1: Create crossover_fine ────────────────────────────────────────────
cat("Running fine-grained intervention search...\n")
## Running fine-grained intervention search...
fine_strengths <- seq(0, 1, by = 0.01)
crossover_fine <- map_dfr(fine_strengths, function(strength) {
m <- tryCatch(
build_intervention_matrix(calibrated_matrix, strength),
error = function(e) NULL
)
if (is.null(m)) return(tibble(strength = strength,
max_nh = NA_real_,
crossover = NA))
result <- run_markov_open(
transition_matrix = m,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
)
nh_max <- result |>
filter(state == "NH") |>
summarise(max_nh = max(n_people)) |>
pull(max_nh)
tibble(
strength = strength,
max_nh = round(nh_max),
crossover = nh_max > national_supply_flat
)
})
cat("Search complete.", nrow(crossover_fine), "scenarios evaluated.\n\n")
## Search complete. 101 scenarios evaluated.
# ── Step 2: Minimum intervention check ───────────────────────────────────────
if (nrow(crossover_fine |> filter(!crossover, !is.na(crossover))) == 0) {
cat("No intervention strength tested fully prevents crossover through 2036.\n")
cat("The new entrant flow alone is sufficient to exceed supply.\n\n")
best_case <- crossover_fine |>
filter(!is.na(crossover)) |>
slice_min(max_nh, n = 1)
cat("Best case (strength =", paste0(best_case$strength * 100, "%):\n"))
cat(" Max NH residents:", comma(best_case$max_nh), "\n")
cat(" Supply ceiling: ", comma(national_supply_flat), "\n")
cat(" Remaining gap: ", comma(best_case$max_nh - national_supply_flat), "\n\n")
} else {
min_intervention <- crossover_fine |>
filter(!crossover) |>
slice_min(strength, n = 1)
cat("Minimum intervention to prevent crossover through 2036:\n")
cat(" Strength: ", paste0(min_intervention$strength * 100, "%"), "\n")
cat(" Max NH residents:", comma(min_intervention$max_nh), "\n")
cat(" Supply ceiling: ", comma(national_supply_flat), "\n\n")
m_min <- build_intervention_matrix(calibrated_matrix, min_intervention$strength)
cat("In plain language:\n")
cat(" Light→Moderate: from",
round(calibrated_matrix["Light", "Moderate"], 4), "to",
round(m_min["Light", "Moderate"], 4), "\n")
cat(" Moderate→Intensive: from",
round(calibrated_matrix["Moderate", "Intensive"], 4), "to",
round(m_min["Moderate", "Intensive"], 4), "\n")
cat(" Intensive→NH: from",
round(calibrated_matrix["Intensive","NH"], 4), "to",
round(m_min["Intensive","NH"], 4), "\n\n")
}
## Minimum intervention to prevent crossover through 2036:
## Strength: 19%
## Max NH residents: 41,517
## Supply ceiling: 41,650
##
## In plain language:
## Light→Moderate: from 0.1154 to 0.0934
## Moderate→Intensive: from 0.1318 to 0.1218
## Intensive→NH: from 0.1958 to 0.1898
# ── Step 3: Gap closure table ─────────────────────────────────────────────────
cat("Gap between max NH residents and supply ceiling by intervention strength:\n")
## Gap between max NH residents and supply ceiling by intervention strength:
crossover_fine |>
filter(strength %in% seq(0, 1, by = 0.10)) |>
select(strength, max_nh) |>
mutate(
gap = max_nh - national_supply_flat,
gap_closed_pct = (1 - gap / first(gap)) * 100
) |>
mutate(
strength = paste0(strength * 100, "%"),
max_nh = comma(max_nh),
gap = comma(gap),
gap_closed_pct = paste0(round(gap_closed_pct, 1), "%")
) |>
kable(
caption = "Gap reduction by intervention strength (strong ageing, flat supply)",
col.names = c("Intervention", "Max NH residents",
"Gap vs supply", "Gap closed vs baseline")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Gap reduction by intervention strength (strong ageing, flat supply)
|
Intervention
|
Max NH residents
|
Gap vs supply
|
Gap closed vs baseline
|
|
0%
|
44,599
|
2,949
|
0%
|
|
10%
|
42,989
|
1,339
|
54.6%
|
|
20%
|
41,352
|
-298
|
110.1%
|
|
40%
|
38,005
|
-3,645
|
223.6%
|
|
50%
|
36,298
|
-5,352
|
281.5%
|
|
70%
|
32,823
|
-8,827
|
399.4%
|
|
80%
|
31,689
|
-9,961
|
437.8%
|
|
90%
|
31,415
|
-10,235
|
447.1%
|
|
100%
|
31,236
|
-10,414
|
453.2%
|
# ── Step 4: Gap closure plot ──────────────────────────────────────────────────
crossover_fine |>
mutate(
gap = max_nh - national_supply_flat,
baseline_gap = first(gap)
) |>
ggplot(aes(x = strength * 100, y = gap)) +
geom_line(colour = "#4472C4", linewidth = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed",
colour = "#C00000", linewidth = 0.8) +
annotate("text", x = 5, y = 200,
label = "Supply ceiling", hjust = 0,
colour = "#C00000", size = 3.5) +
scale_x_continuous(breaks = seq(0, 100, 10),
labels = paste0(seq(0, 100, 10), "%")) +
scale_y_continuous(labels = comma) +
labs(
title = "Gap between NH demand and supply — by intervention strength",
subtitle = "Even at 100% intervention the gap is nearly closed but not fully eliminated",
x = "Intervention strength (reduction in progression rates)",
y = "Remaining gap (NH residents above supply ceiling)"
) +
theme_minimal(base_size = 13) +
theme(panel.grid.minor = element_blank())
# ── Step 5: Combined intervention + supply expansion ─────────────────────────
cat("\n=== COMBINED INTERVENTION + SUPPLY EXPANSION ===\n")
##
## === COMBINED INTERVENTION + SUPPLY EXPANSION ===
crossover_fine |>
filter(strength %in% seq(0, 1, by = 0.10)) |>
select(strength, max_nh) |>
mutate(
beds_needed = pmax(max_nh - national_supply_flat, 0),
pct_of_current = beds_needed / national_supply_flat * 100,
annual_beds_needed = round(beds_needed / 11)
) |>
mutate(
strength = paste0(strength * 100, "%"),
max_nh = comma(round(max_nh)),
beds_needed = comma(beds_needed),
pct_of_current = paste0(round(pct_of_current, 1), "%"),
annual_beds_needed = comma(annual_beds_needed)
) |>
kable(
caption = "Additional beds needed alongside each intervention level",
col.names = c("Intervention", "Max NH residents", "Additional beds needed",
"% of current supply", "Avg beds/year needed")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Additional beds needed alongside each intervention level
|
Intervention
|
Max NH residents
|
Additional beds needed
|
% of current supply
|
Avg beds/year needed
|
|
0%
|
44,599
|
2,949
|
7.1%
|
268
|
|
10%
|
42,989
|
1,339
|
3.2%
|
122
|
|
20%
|
41,352
|
0
|
0%
|
0
|
|
40%
|
38,005
|
0
|
0%
|
0
|
|
50%
|
36,298
|
0
|
0%
|
0
|
|
70%
|
32,823
|
0
|
0%
|
0
|
|
80%
|
31,689
|
0
|
0%
|
0
|
|
90%
|
31,415
|
0
|
0%
|
0
|
|
100%
|
31,236
|
0
|
0%
|
0
|
# ── Step 6: Commercial argument ───────────────────────────────────────────────
cat("\n=== THE COMMERCIAL ARGUMENT IN PLAIN NUMBERS ===\n\n")
##
## === THE COMMERCIAL ARGUMENT IN PLAIN NUMBERS ===
for (pct in c(0, 10, 20, 30, 40)) {
row <- crossover_fine |> filter(strength == pct / 100)
if (nrow(row) == 0) next
gap <- row$max_nh - national_supply_flat
cat(paste0(" ", pct, "% intervention: max NH = ", comma(row$max_nh),
", gap = ", comma(round(gap)),
" beds (",
round((1 - gap / (crossover_fine$max_nh[1] - national_supply_flat)) * 100, 1),
"% of gap closed)\n"))
}
## 0% intervention: max NH = 44,599, gap = 2,949 beds (0% of gap closed)
## 10% intervention: max NH = 42,989, gap = 1,339 beds (54.6% of gap closed)
## 20% intervention: max NH = 41,352, gap = -298 beds (110.1% of gap closed)
## 30% intervention: max NH = 39,690, gap = -1,960 beds (166.5% of gap closed)
## 40% intervention: max NH = 38,005, gap = -3,645 beds (223.6% of gap closed)
Adding costs
# ════════════════════════════════════════════════════════════════════════════
# PART D — COST INTEGRATION
# ════════════════════════════════════════════════════════════════════════════
# ── D1: Define unit costs ─────────────────────────────────────────────────────
unit_costs <- c(
Light = 67921, # Small need of assistance
Moderate = 180340, # Average-to-much need
Intensive = 657041, # Extensive need
NH = 1585357, # Nursing home placement
Dead = 0 # No cost
)
# Co-located housing cost (between intensive home care and NH)
# This is not a Markov state in our model but used as a reference benchmark
cost_collocated <- 1238172
cat("=== UNIT COSTS PER STATE (NOK/person/year) ===\n")
## === UNIT COSTS PER STATE (NOK/person/year) ===
tibble(
State = names(unit_costs),
Cost_NOK = unit_costs,
Cost_formatted = scales::comma(unit_costs)
) |>
kable(caption = "Unit costs by care state",
col.names = c("State", "NOK/person/year", "Formatted")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Unit costs by care state
|
State
|
NOK/person/year
|
Formatted
|
|
Light
|
67921
|
67,921
|
|
Moderate
|
180340
|
180,340
|
|
Intensive
|
657041
|
657,041
|
|
NH
|
1585357
|
1,585,357
|
|
Dead
|
0
|
0
|
# ── D2: Cost function — apply unit costs to Markov output ────────────────────
apply_costs <- function(markov_df) {
markov_df |>
mutate(
unit_cost = unit_costs[as.character(state)],
total_cost = n_people * unit_cost
)
}
# ── D3: Total annual system cost — baseline vs intervention scenarios ─────────
# Run cost calculation for all intervention scenarios from Part C
costed_scenarios <- scenarios_all |>
apply_costs() |>
group_by(intervention, year) |>
summarise(
total_cost_NOK = sum(total_cost, na.rm = TRUE),
total_people = sum(n_people[state != "Dead"], na.rm = TRUE),
nh_residents = n_people[state == "NH"],
.groups = "drop"
) |>
mutate(
total_cost_MNOK = total_cost_NOK / 1e6,
intervention = factor(intervention,
levels = paste0(c(0,10,20,30,40), "% reduction"))
)
cat("\nTotal annual system cost by scenario (NOK million):\n")
##
## Total annual system cost by scenario (NOK million):
costed_scenarios |>
select(intervention, year, total_cost_MNOK) |>
pivot_wider(names_from = intervention, values_from = total_cost_MNOK) |>
mutate(across(where(is.numeric), ~ round(.x, 0))) |>
mutate(across(where(is.numeric), comma)) |>
kable(caption = "Total annual system cost — all states (NOK million)",
col.names = c("Year", "0% reduction", "10% reduction",
"20% reduction", "30% reduction", "40% reduction")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Total annual system cost — all states (NOK million)
|
Year
|
0% reduction
|
10% reduction
|
20% reduction
|
30% reduction
|
40% reduction
|
|
2,025
|
75,858
|
75,858
|
75,858
|
75,858
|
75,858
|
|
2,026
|
78,151
|
77,800
|
77,448
|
77,096
|
76,745
|
|
2,027
|
79,355
|
78,681
|
78,007
|
77,335
|
76,663
|
|
2,028
|
80,537
|
79,522
|
78,509
|
77,497
|
76,487
|
|
2,029
|
82,315
|
80,926
|
79,535
|
78,145
|
76,754
|
|
2,030
|
84,985
|
83,190
|
81,386
|
79,576
|
77,760
|
|
2,031
|
88,629
|
86,409
|
84,168
|
81,907
|
79,627
|
|
2,032
|
93,198
|
90,552
|
87,864
|
85,135
|
82,368
|
|
2,033
|
98,568
|
95,511
|
92,384
|
89,188
|
85,926
|
|
2,034
|
104,584
|
101,144
|
97,599
|
93,952
|
90,201
|
# ── D4: Cost savings vs baseline ─────────────────────────────────────────────
baseline_cost <- costed_scenarios |>
filter(intervention == "0% reduction") |>
select(year, baseline_cost = total_cost_NOK)
cost_savings <- costed_scenarios |>
left_join(baseline_cost, by = "year") |>
mutate(
annual_saving_NOK = baseline_cost - total_cost_NOK,
annual_saving_MNOK = annual_saving_NOK / 1e6
) |>
filter(intervention != "0% reduction")
cat("\nAnnual cost savings vs baseline (NOK million):\n")
##
## Annual cost savings vs baseline (NOK million):
cost_savings |>
select(intervention, year, annual_saving_MNOK) |>
pivot_wider(names_from = intervention, values_from = annual_saving_MNOK) |>
mutate(across(where(is.numeric), ~ round(.x, 1))) |>
mutate(across(where(is.numeric), comma)) |>
kable(caption = "Annual cost savings vs baseline (NOK million)",
col.names = c("Year", "10% reduction", "20% reduction",
"30% reduction", "40% reduction")) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Annual cost savings vs baseline (NOK million)
|
Year
|
10% reduction
|
20% reduction
|
30% reduction
|
40% reduction
|
|
2,025
|
0
|
0
|
0
|
0
|
|
2,026
|
352
|
703
|
1,055
|
1,406
|
|
2,027
|
674
|
1,348
|
2,020
|
2,692
|
|
2,028
|
1,014
|
2,028
|
3,039
|
4,050
|
|
2,029
|
1,389
|
2,779
|
4,170
|
5,560
|
|
2,030
|
1,795
|
3,598
|
5,409
|
7,225
|
|
2,031
|
2,220
|
4,461
|
6,722
|
9,002
|
|
2,032
|
2,646
|
5,334
|
8,063
|
10,830
|
|
2,033
|
3,058
|
6,185
|
9,380
|
12,642
|
|
2,034
|
3,441
|
6,985
|
10,632
|
14,383
|
# ── D5: Cumulative cost savings 2025-2034 ────────────────────────────────────
cumulative_savings <- cost_savings |>
group_by(intervention) |>
summarise(
cumulative_saving_BNOK = sum(annual_saving_NOK, na.rm = TRUE) / 1e9,
avg_annual_saving_MNOK = mean(annual_saving_NOK, na.rm = TRUE) / 1e6,
saving_by_2034_MNOK = annual_saving_NOK[year == 2034] / 1e6,
.groups = "drop"
)
cat("\nCumulative cost savings 2025-2034:\n")
##
## Cumulative cost savings 2025-2034:
cumulative_savings |>
mutate(across(where(is.numeric), ~ round(.x, 1))) |>
kable(
caption = "Cumulative cost savings 2025-2034 by intervention strength",
col.names = c("Intervention", "Cumulative saving (NOK billion)",
"Average annual saving (NOK million)",
"Annual saving by 2034 (NOK million)")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Cumulative cost savings 2025-2034 by intervention strength
|
Intervention
|
Cumulative saving (NOK billion)
|
Average annual saving (NOK million)
|
Annual saving by 2034 (NOK million)
|
|
10% reduction
|
16.6
|
1658.8
|
3440.6
|
|
20% reduction
|
33.4
|
3342.0
|
6984.8
|
|
30% reduction
|
50.5
|
5049.0
|
10632.4
|
|
40% reduction
|
67.8
|
6779.1
|
14382.9
|
# ── D6: Cost per state breakdown — where does the money go? ──────────────────
cost_by_state <- scenarios_all |>
apply_costs() |>
filter(intervention %in% c("0% reduction", "30% reduction"),
state != "Dead") |>
group_by(intervention, year, state) |>
summarise(total_cost_MNOK = sum(total_cost, na.rm = TRUE) / 1e6,
.groups = "drop")
cost_by_state |>
filter(year %in% c(2025, 2029, 2034)) |>
mutate(total_cost_MNOK = round(total_cost_MNOK, 0)) |>
pivot_wider(names_from = year, values_from = total_cost_MNOK) |>
arrange(intervention, state) |>
kable(
caption = "Cost breakdown by state — baseline vs 30% intervention (NOK million)",
col.names = c("Intervention", "State", "2025", "2029", "2034")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Cost breakdown by state — baseline vs 30% intervention (NOK million)
|
Intervention
|
State
|
2025
|
2029
|
2034
|
|
0% reduction
|
Light
|
2825
|
7648
|
11103
|
|
0% reduction
|
Moderate
|
6091
|
5883
|
9557
|
|
0% reduction
|
Intensive
|
18395
|
11823
|
13218
|
|
0% reduction
|
NH
|
48547
|
56961
|
70706
|
|
30% reduction
|
Light
|
2825
|
8248
|
12684
|
|
30% reduction
|
Moderate
|
6091
|
5023
|
7955
|
|
30% reduction
|
Intensive
|
18395
|
10972
|
10389
|
|
30% reduction
|
NH
|
48547
|
53902
|
62923
|
# ── D7: NH cost specifically — the avoided bed-cost argument ─────────────────
nh_cost_comparison <- costed_scenarios |>
mutate(nh_cost_MNOK = nh_residents * unit_costs["NH"] / 1e6) |>
select(intervention, year, nh_residents, nh_cost_MNOK)
baseline_nh_cost <- nh_cost_comparison |>
filter(intervention == "0% reduction") |>
select(year, baseline_nh_cost = nh_cost_MNOK,
baseline_residents = nh_residents)
nh_savings <- nh_cost_comparison |>
left_join(baseline_nh_cost, by = "year") |>
mutate(
nh_cost_saving_MNOK = baseline_nh_cost - nh_cost_MNOK,
residents_avoided = baseline_residents - nh_residents
) |>
filter(intervention != "0% reduction", year >= 2029)
cat("\nNH cost savings specifically (NOK million) — from crossover year:\n")
##
## NH cost savings specifically (NOK million) — from crossover year:
nh_savings |>
select(intervention, year, residents_avoided, nh_cost_saving_MNOK) |>
mutate(nh_cost_saving_MNOK = round(nh_cost_saving_MNOK, 0),
residents_avoided = round(residents_avoided)) |>
kable(
caption = "NH cost savings vs baseline — from 2029 onwards",
col.names = c("Intervention", "Year", "Residents avoided", "NH cost saving (M NOK)")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), font_size = 11)
NH cost savings vs baseline — from 2029 onwards
|
Intervention
|
Year
|
Residents avoided
|
NH cost saving (M NOK)
|
|
10% reduction
|
2029
|
646
|
1025
|
|
10% reduction
|
2030
|
824
|
1307
|
|
10% reduction
|
2031
|
1016
|
1611
|
|
10% reduction
|
2032
|
1216
|
1928
|
|
10% reduction
|
2033
|
1416
|
2246
|
|
10% reduction
|
2034
|
1611
|
2553
|
|
20% reduction
|
2029
|
1290
|
2045
|
|
20% reduction
|
2030
|
1646
|
2610
|
|
20% reduction
|
2031
|
2032
|
3222
|
|
20% reduction
|
2032
|
2437
|
3864
|
|
20% reduction
|
2033
|
2847
|
4514
|
|
20% reduction
|
2034
|
3247
|
5148
|
|
30% reduction
|
2029
|
1930
|
3059
|
|
30% reduction
|
2030
|
2465
|
3908
|
|
30% reduction
|
2031
|
3048
|
4832
|
|
30% reduction
|
2032
|
3663
|
5808
|
|
30% reduction
|
2033
|
4291
|
6803
|
|
30% reduction
|
2034
|
4909
|
7783
|
|
40% reduction
|
2029
|
2566
|
4068
|
|
40% reduction
|
2030
|
3280
|
5200
|
|
40% reduction
|
2031
|
4062
|
6440
|
|
40% reduction
|
2032
|
4893
|
7757
|
|
40% reduction
|
2033
|
5746
|
9110
|
|
40% reduction
|
2034
|
6594
|
10455
|
# ── D8: The commercial number — cost per avoided NH admission ─────────────────
# If an intervention achieves 30% stabilisation, what is the cost saving
# per person kept out of a nursing home per year?
# This is the number to put in front of a municipality.
cat("\n=== THE COMMERCIAL ARGUMENT IN COST TERMS ===\n\n")
##
## === THE COMMERCIAL ARGUMENT IN COST TERMS ===
for (int_level in c("10% reduction", "20% reduction", "30% reduction")) {
total_avoided <- nh_savings |>
filter(intervention == int_level) |>
summarise(
total_residents_avoided = sum(residents_avoided),
total_nh_saving_MNOK = sum(nh_cost_saving_MNOK)
)
total_system_saving <- cost_savings |>
filter(intervention == int_level) |>
summarise(total = sum(annual_saving_NOK)) |>
pull(total)
cat("Intervention:", int_level, "\n")
cat(" NH residents avoided (2029-2034): ",
comma(round(total_avoided$total_residents_avoided)), "\n")
cat(" NH cost saving (2029-2034): NOK",
comma(round(total_avoided$total_nh_saving_MNOK)), "million\n")
cat(" Total system cost saving (all years): NOK",
comma(round(total_system_saving / 1e6)), "million\n")
cat(" Cost saving per avoided NH resident: NOK",
comma(round(total_avoided$total_nh_saving_MNOK * 1e6 /
max(total_avoided$total_residents_avoided, 1))), "\n\n")
}
## Intervention: 10% reduction
## NH residents avoided (2029-2034): 6,730
## NH cost saving (2029-2034): NOK 10,669 million
## Total system cost saving (all years): NOK 16,588 million
## Cost saving per avoided NH resident: NOK 1,585,357
##
## Intervention: 20% reduction
## NH residents avoided (2029-2034): 13,500
## NH cost saving (2029-2034): NOK 21,403 million
## Total system cost saving (all years): NOK 33,420 million
## Cost saving per avoided NH resident: NOK 1,585,357
##
## Intervention: 30% reduction
## NH residents avoided (2029-2034): 20,306
## NH cost saving (2029-2034): NOK 32,192 million
## Total system cost saving (all years): NOK 50,490 million
## Cost saving per avoided NH resident: NOK 1,585,357
# ── D9: Visualisation — total system cost trajectories ───────────────────────
costed_scenarios |>
ggplot(aes(x = year, y = total_cost_MNOK,
colour = intervention, group = intervention)) +
geom_line(linewidth = 1.1) +
scale_colour_manual(values = c(
"0% reduction" = "#C00000",
"10% reduction" = "#ED7D31",
"20% reduction" = "#FFC000",
"30% reduction" = "#70AD47",
"40% reduction" = "#4472C4"
)) +
scale_x_continuous(breaks = 2025:2034) +
scale_y_continuous(labels = ~ paste0(comma(.x / 1000), "B")) +
labs(
title = "Total system care cost — by intervention strength",
subtitle = "All home care states + nursing home. Strong ageing scenario.",
x = NULL, y = "Total cost (NOK billion)", colour = "Intervention"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank())
# ── D10: Visualisation — annual cost saving vs baseline ──────────────────────
cost_savings |>
ggplot(aes(x = year, y = annual_saving_MNOK,
colour = intervention, group = intervention)) +
geom_line(linewidth = 1.1) +
geom_point(size = 2) +
scale_colour_manual(values = c(
"10% reduction" = "#ED7D31",
"20% reduction" = "#FFC000",
"30% reduction" = "#70AD47",
"40% reduction" = "#4472C4"
)) +
scale_x_continuous(breaks = 2025:2034) +
scale_y_continuous(labels = comma) +
labs(
title = "Annual cost saving vs baseline — by intervention strength",
subtitle = "Saving grows over time as intervention effect compounds through the pipeline",
x = NULL, y = "Annual saving (NOK million)", colour = "Intervention"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom",
panel.grid.minor = element_blank())
Sensitivity
Analysis:
# ════════════════════════════════════════════════════════════════════════════
# PART E — SENSITIVITY ANALYSIS
# How much do results change if transition probabilities are higher or lower?
# ════════════════════════════════════════════════════════════════════════════
# ── E1: Define sensitivity scenarios ─────────────────────────────────────────
# We vary all transition probabilities simultaneously by a percentage
# Low = -20% (optimistic — slower progression than calibrated)
# Base = 0% (calibrated values)
# High = +20% (pessimistic — faster progression than calibrated)
#
# We also test individual parameter uncertainty for the three most
# important transitions: light→moderate, moderate→intensive, intensive→NH
sensitivity_multipliers <- c(
"Low (-20%)" = 0.80,
"Base (0%)" = 1.00,
"High (+20%)" = 1.20
)
# ── E2: Build sensitivity matrices ───────────────────────────────────────────
build_sensitivity_matrix <- function(multiplier) {
# Scale all progression and NH admission probabilities
# Keep mortality fixed — these come from life tables not calibration
p_l2m_new <- calibrated_matrix["Light", "Moderate"] * multiplier
p_m2i_new <- calibrated_matrix["Moderate", "Intensive"] * multiplier
p_l2nh_new <- calibrated_matrix["Light", "NH"] * multiplier
p_m2nh_new <- calibrated_matrix["Moderate", "NH"] * multiplier
p_i2nh_new <- calibrated_matrix["Intensive","NH"] * multiplier
# Validate: probabilities cannot exceed row constraints
# If scaling up pushes a row over 1, cap at 95% of the headroom
cap <- function(p, max_val) min(p, max_val * 0.95)
headroom_light <- 1 - p_die_light_obs
headroom_moderate <- 1 - p_die_moderate_obs
headroom_intense <- 1 - p_die_intensive_obs
p_l2m_new <- cap(p_l2m_new, headroom_light - p_l2nh_new)
p_m2i_new <- cap(p_m2i_new, headroom_moderate - p_m2nh_new)
p_i2nh_new <- cap(p_i2nh_new, headroom_intense)
build_transition_matrix(
p_light_to_moderate = p_l2m_new,
p_moderate_to_intense = p_m2i_new,
p_light_to_nh = p_l2nh_new,
p_moderate_to_nh = p_m2nh_new,
p_intense_to_nh = p_i2nh_new,
p_die_light = p_die_light_obs,
p_die_moderate = p_die_moderate_obs,
p_die_intense = p_die_intensive_obs,
p_die_nh = p_die_nh_obs
)
}
# Show what the sensitivity scenarios do to key probabilities
cat("=== SENSITIVITY SCENARIO TRANSITION PROBABILITIES ===\n")
## === SENSITIVITY SCENARIO TRANSITION PROBABILITIES ===
map_dfr(names(sensitivity_multipliers), function(label) {
mult <- sensitivity_multipliers[label]
m <- build_sensitivity_matrix(mult)
tibble(
scenario = label,
light_to_moderate = round(m["Light", "Moderate"], 4),
moderate_to_intensive = round(m["Moderate", "Intensive"], 4),
intensive_to_nh = round(m["Intensive","NH"], 4)
)
}) |>
kable(
caption = "Key transition probabilities under each sensitivity scenario",
col.names = c("Scenario", "Light→Moderate",
"Moderate→Intensive", "Intensive→NH")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Key transition probabilities under each sensitivity scenario
|
Scenario
|
Light→Moderate
|
Moderate→Intensive
|
Intensive→NH
|
|
Low (-20%)
|
0.0923
|
0.1055
|
0.1566
|
|
Base (0%)
|
0.1154
|
0.1318
|
0.1958
|
|
High (+20%)
|
0.1384
|
0.1582
|
0.2350
|
# ── E3: Run open cohort model for each sensitivity scenario ──────────────────
cat("\nRunning sensitivity scenarios...\n")
##
## Running sensitivity scenarios...
sensitivity_runs <- map_dfr(names(sensitivity_multipliers), function(label) {
mult <- sensitivity_multipliers[label]
m <- build_sensitivity_matrix(mult)
run_markov_open(
transition_matrix = m,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
) |>
mutate(sensitivity = label)
})
cat("Done.\n")
## Done.
# ── E4: NH trajectory under each sensitivity scenario ────────────────────────
nh_sensitivity <- sensitivity_runs |>
filter(state == "NH") |>
arrange(sensitivity, year) |>
group_by(sensitivity) |>
mutate(
nh_admissions = n_people - lag(n_people, default = first(n_people)),
supply = national_supply_flat,
gap = n_people - supply,
over_capacity = gap > 0
) |>
ungroup() |>
mutate(sensitivity = factor(sensitivity,
levels = names(sensitivity_multipliers)))
# Crossover year by sensitivity scenario
crossover_sensitivity <- nh_sensitivity |>
filter(over_capacity) |>
group_by(sensitivity) |>
summarise(crossover_year = min(year), .groups = "drop")
cat("\nCrossover year by sensitivity scenario:\n")
##
## Crossover year by sensitivity scenario:
crossover_sensitivity |>
kable(
caption = "Crossover year — sensitivity analysis",
col.names = c("Scenario", "Crossover year")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Crossover year — sensitivity analysis
|
Scenario
|
Crossover year
|
|
Base (0%)
|
2033
|
|
High (+20%)
|
2030
|
# ── E5: NH residents table — sensitivity scenarios ───────────────────────────
nh_sensitivity |>
select(sensitivity, year, n_people, gap) |>
mutate(
n_people = comma(round(n_people)),
gap = comma(round(gap))
) |>
pivot_wider(names_from = sensitivity,
values_from = c(n_people, gap),
names_glue = "{sensitivity}_{.value}") |>
select(year,
starts_with("Low"),
starts_with("Base"),
starts_with("High")) |>
kable(
caption = "NH residents and gap vs supply — sensitivity scenarios",
col.names = c("Year",
"Low: residents", "Low: gap",
"Base: residents", "Base: gap",
"High: residents", "High: gap")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), font_size = 11)
NH residents and gap vs supply — sensitivity scenarios
|
Year
|
Low: residents
|
Low: gap
|
Base: residents
|
Base: gap
|
High: residents
|
High: gap
|
|
2025
|
30,622
|
-11,028
|
30,622
|
-11,028
|
30,622
|
-11,028
|
|
2026
|
31,016
|
-10,634
|
33,029
|
-8,622
|
35,041
|
-6,609
|
|
2027
|
31,014
|
-10,636
|
34,326
|
-7,324
|
37,583
|
-4,067
|
|
2028
|
30,953
|
-10,698
|
35,149
|
-6,501
|
39,220
|
-2,430
|
|
2029
|
31,054
|
-10,596
|
35,929
|
-5,721
|
40,615
|
-1,036
|
|
2030
|
31,451
|
-10,200
|
36,927
|
-4,723
|
42,164
|
514
|
|
2031
|
32,206
|
-9,444
|
38,275
|
-3,376
|
44,064
|
2,413
|
|
2032
|
33,336
|
-8,314
|
40,017
|
-1,633
|
46,374
|
4,723
|
|
2033
|
34,824
|
-6,826
|
42,141
|
490
|
49,073
|
7,422
|
|
2034
|
36,635
|
-5,015
|
44,599
|
2,949
|
52,096
|
10,446
|
# ── E6: Cost sensitivity ──────────────────────────────────────────────────────
cost_sensitivity <- sensitivity_runs |>
mutate(
unit_cost = unit_costs[as.character(state)],
total_cost = n_people * unit_cost
) |>
group_by(sensitivity, year) |>
summarise(
total_cost_BNOK = sum(total_cost, na.rm = TRUE) / 1e9,
nh_residents = sum(n_people[state == "NH"], na.rm = TRUE),
.groups = "drop"
) |>
mutate(sensitivity = factor(sensitivity,
levels = names(sensitivity_multipliers)))
# Baseline cost for reference
base_cost <- cost_sensitivity |>
filter(sensitivity == "Base (0%)") |>
select(year, base_cost = total_cost_BNOK)
cost_sensitivity <- cost_sensitivity |>
left_join(base_cost, by = "year") |>
mutate(cost_vs_base_MNOK = (total_cost_BNOK - base_cost) * 1000)
cat("\nTotal system cost by sensitivity scenario (NOK billion):\n")
##
## Total system cost by sensitivity scenario (NOK billion):
cost_sensitivity |>
select(sensitivity, year, total_cost_BNOK) |>
mutate(total_cost_BNOK = round(total_cost_BNOK, 2)) |>
pivot_wider(names_from = sensitivity, values_from = total_cost_BNOK) |>
kable(
caption = "Total system cost — sensitivity scenarios (NOK billion)",
col.names = c("Year", "Low (-20%)", "Base (0%)", "High (+20%)")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Total system cost — sensitivity scenarios (NOK billion)
|
Year
|
Low (-20%)
|
Base (0%)
|
High (+20%)
|
|
2025
|
75.86
|
75.86
|
75.86
|
|
2026
|
78.15
|
81.01
|
75.29
|
|
2027
|
79.35
|
83.96
|
74.66
|
|
2028
|
80.54
|
86.34
|
74.53
|
|
2029
|
82.31
|
89.12
|
75.21
|
|
2030
|
84.98
|
92.75
|
76.83
|
|
2031
|
88.63
|
97.35
|
79.41
|
|
2032
|
93.20
|
102.87
|
82.89
|
|
2033
|
98.57
|
109.14
|
87.18
|
|
2034
|
104.58
|
115.99
|
92.14
|
# ── E7: Individual parameter sensitivity — tornado analysis ──────────────────
# Test each transition probability individually:
# what happens if just one parameter is 20% higher/lower?
cat("\nRunning tornado analysis (individual parameter sensitivity)...\n")
##
## Running tornado analysis (individual parameter sensitivity)...
params_to_test <- list(
"Light→Moderate" = "p_light_to_moderate",
"Moderate→Intensive"= "p_moderate_to_intense",
"Intensive→NH" = "p_intense_to_nh",
"Light→NH" = "p_light_to_nh",
"Moderate→NH" = "p_moderate_to_nh"
)
base_values <- list(
p_light_to_moderate = calibrated_matrix["Light", "Moderate"],
p_moderate_to_intense = calibrated_matrix["Moderate", "Intensive"],
p_intense_to_nh = calibrated_matrix["Intensive","NH"],
p_light_to_nh = calibrated_matrix["Light", "NH"],
p_moderate_to_nh = calibrated_matrix["Moderate", "NH"]
)
tornado_results <- map_dfr(names(params_to_test), function(param_label) {
param_name <- params_to_test[[param_label]]
map_dfr(c(low = 0.80, high = 1.20), function(mult) {
# Build a modified parameter list
p <- base_values
p[[param_name]] <- p[[param_name]] * mult
m <- tryCatch(
build_transition_matrix(
p_light_to_moderate = p$p_light_to_moderate,
p_moderate_to_intense = p$p_moderate_to_intense,
p_light_to_nh = p$p_light_to_nh,
p_moderate_to_nh = p$p_moderate_to_nh,
p_intense_to_nh = p$p_intense_to_nh,
p_die_light = p_die_light_obs,
p_die_moderate = p_die_moderate_obs,
p_die_intense = p_die_intensive_obs,
p_die_nh = p_die_nh_obs
),
error = function(e) NULL
)
if (is.null(m)) return(tibble())
result <- run_markov_open(
transition_matrix = m,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
)
max_nh <- result |>
filter(state == "NH") |>
summarise(max_nh = max(n_people)) |>
pull(max_nh)
tibble(
parameter = param_label,
direction = if_else(mult < 1, "Low (-20%)", "High (+20%)"),
max_nh = round(max_nh),
mult = mult
)
}, .id = "direction_id")
})
# Add baseline for comparison
base_max_nh <- sensitivity_runs |>
filter(sensitivity == "Base (0%)", state == "NH") |>
summarise(max_nh = max(n_people)) |>
pull(max_nh)
tornado_results <- tornado_results |>
mutate(
base_max_nh = round(base_max_nh),
diff_from_base = max_nh - base_max_nh
)
cat("\nTornado analysis — effect of individual parameter variation on max NH residents:\n")
##
## Tornado analysis — effect of individual parameter variation on max NH residents:
tornado_results |>
select(parameter, direction, max_nh, diff_from_base) |>
arrange(parameter, direction) |>
mutate(
max_nh = comma(max_nh),
diff_from_base = paste0(if_else(diff_from_base > 0, "+", ""),
comma(diff_from_base))
) |>
kable(
caption = paste("Tornado analysis — baseline max NH residents:",
comma(round(base_max_nh))),
col.names = c("Parameter", "Direction",
"Max NH residents", "Difference from base")
) |>
kable_styling(bootstrap_options = c("striped", "hover"))
Tornado analysis — baseline max NH residents: 44,599
|
Parameter
|
Direction
|
Max NH residents
|
Difference from base
|
|
Intensive→NH
|
High (+20%)
|
45,549
|
+950
|
|
Intensive→NH
|
Low (-20%)
|
43,305
|
-1,294
|
|
Light→Moderate
|
High (+20%)
|
46,512
|
+1,913
|
|
Light→Moderate
|
Low (-20%)
|
42,490
|
-2,109
|
|
Light→NH
|
High (+20%)
|
46,849
|
+2,250
|
|
Light→NH
|
Low (-20%)
|
42,274
|
-2,325
|
|
Moderate→Intensive
|
High (+20%)
|
44,992
|
+393
|
|
Moderate→Intensive
|
Low (-20%)
|
44,152
|
-447
|
|
Moderate→NH
|
High (+20%)
|
46,429
|
+1,830
|
|
Moderate→NH
|
Low (-20%)
|
42,563
|
-2,036
|
# ── E8: Tornado plot ──────────────────────────────────────────────────────────
tornado_plot_data <- tornado_results |>
select(parameter, direction, diff_from_base) |>
pivot_wider(names_from = direction, values_from = diff_from_base) |>
rename(low = `Low (-20%)`,
high = `High (+20%)`) |>
mutate(
range = abs(high) + abs(low),
parameter = fct_reorder(parameter, range)
)
tornado_plot_data |>
ggplot() +
geom_segment(
aes(x = low, xend = high,
y = parameter, yend = parameter),
colour = "grey70", linewidth = 6
) +
geom_segment(
aes(x = pmin(low, 0), xend = pmax(high, 0),
y = parameter, yend = parameter,
colour = high > 0),
linewidth = 6, alpha = 0.7
) +
geom_vline(xintercept = 0, linewidth = 0.8, colour = "black") +
scale_colour_manual(values = c("TRUE" = "#C00000", "FALSE" = "#4472C4"),
guide = "none") +
scale_x_continuous(labels = comma) +
labs(
title = "Tornado diagram — sensitivity of max NH residents",
subtitle = paste("Each bar shows impact of ±20% change in one parameter.",
"Baseline max NH residents:", comma(round(base_max_nh))),
x = "Change in maximum NH residents vs baseline",
y = NULL
) +
theme_minimal(base_size = 13) +
theme(panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank())
# ── E9: Sensitivity of crossover year to intervention under each scenario ─────
# Does the 30% intervention still work under pessimistic assumptions?
cat("\nSensitivity of intervention effect across uncertainty scenarios:\n")
##
## Sensitivity of intervention effect across uncertainty scenarios:
sensitivity_intervention_runs <- map_dfr(
names(sensitivity_multipliers),
function(sens_label) {
sens_mult <- sensitivity_multipliers[sens_label]
sens_base <- build_sensitivity_matrix(sens_mult)
map_dfr(c(0, 0.20, 0.30, 0.40), function(int_strength) {
# Apply intervention on top of sensitivity scenario
m <- tryCatch(
build_intervention_matrix(sens_base, int_strength),
error = function(e) NULL
)
if (is.null(m)) return(tibble())
result <- run_markov_open(
transition_matrix = m,
cohort_start = cohort_2025,
new_entrants_df = new_entrants_proj,
scenario_name = "Strong ageing",
start_year = 2025
)
nh_traj <- result |>
filter(state == "NH") |>
mutate(over_capacity = n_people > national_supply_flat)
crossover <- nh_traj |>
filter(over_capacity) |>
summarise(crossover = min(year)) |>
pull(crossover)
tibble(
sensitivity = sens_label,
intervention = paste0(int_strength * 100, "% reduction"),
crossover_year = if_else(length(crossover) == 0,
"No crossover", as.character(crossover)),
max_nh = round(max(nh_traj$n_people))
)
})
}
)
sensitivity_intervention_runs |>
select(sensitivity, intervention, crossover_year, max_nh) |>
mutate(max_nh = comma(max_nh)) |>
pivot_wider(names_from = sensitivity,
values_from = c(crossover_year, max_nh),
names_glue = "{sensitivity}_{.value}") |>
select(intervention,
starts_with("Low"),
starts_with("Base"),
starts_with("High")) |>
kable(
caption = "Crossover year and max NH residents — intervention × sensitivity",
col.names = c("Intervention",
"Low: crossover", "Low: max NH",
"Base: crossover", "Base: max NH",
"High: crossover", "High: max NH")
) |>
kable_styling(bootstrap_options = c("striped", "hover"), font_size = 11)
Crossover year and max NH residents — intervention × sensitivity
|
Intervention
|
Low: crossover
|
Low: max NH
|
Base: crossover
|
Base: max NH
|
High: crossover
|
High: max NH
|
|
0% reduction
|
Inf
|
36,635
|
2033
|
44,599
|
2030
|
52,096
|
|
20% reduction
|
Inf
|
34,092
|
Inf
|
41,352
|
2032
|
48,228
|
|
30% reduction
|
Inf
|
32,810
|
Inf
|
39,690
|
2032
|
46,217
|
|
40% reduction
|
Inf
|
31,520
|
Inf
|
38,005
|
2033
|
44,155
|
# ════════════════════════════════════════════════════════════════════════════
# PART F — EXPORT DATA FOR INTERACTIVE DASHBOARD
# ════════════════════════════════════════════════════════════════════════════
library(jsonlite)
# Create data folder if it does not exist
dir.create("data", showWarnings = FALSE)
# ── F1: Model parameters ──────────────────────────────────────────────────────
model_params <- list(
national_supply_flat = national_supply_flat,
supply_growth_rate = SUPPLY_GROWTH_RATE,
expansion_annual = EXPANSION_ANNUAL,
hc_entry_rate = hc_entry_rate,
avg_hc_duration = avg_hc_duration_years,
avg_nh_stay = avg_nh_stay_years,
unit_costs = as.list(unit_costs),
calibrated_probs = list(
light_to_moderate = round(calibrated_matrix["Light", "Moderate"], 4),
moderate_to_intensive= round(calibrated_matrix["Moderate", "Intensive"], 4),
intensive_to_nh = round(calibrated_matrix["Intensive","NH"], 4),
light_to_nh = round(calibrated_matrix["Light", "NH"], 4),
moderate_to_nh = round(calibrated_matrix["Moderate", "NH"], 4),
die_light = p_die_light_obs,
die_moderate = p_die_moderate_obs,
die_intensive = p_die_intensive_obs,
die_nh = p_die_nh_obs
),
calibration_targets = list(
observed_light = observed_light,
observed_moderate = observed_moderate,
observed_intensive = observed_intensive,
observed_nh = observed_nh,
observed_nh_admissions = observed_nh_admissions_annual
),
cohort_start = as.list(cohort_2025)
)
write_json(model_params, "data/model_params.json", pretty = TRUE, auto_unbox = TRUE)
cat("Exported model_params.json\n")
## Exported model_params.json
# ── F2: NH trajectory — all intervention scenarios ────────────────────────────
nh_trajectories <- nh_scenarios |>
select(intervention, year, nh_residents = n_people,
supply, gap, over_capacity) |>
mutate(
nh_residents = round(nh_residents),
gap = round(gap)
)
write_json(nh_trajectories, "data/nh_trajectories.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported nh_trajectories.json\n")
## Exported nh_trajectories.json
# ── F3: Open cohort full state populations (baseline + 30% intervention) ──────
state_populations <- open_combined |>
filter(scenario == "Strong ageing") |>
bind_rows(
scenarios_all |>
filter(intervention == "30% reduction") |>
mutate(scenario = "Strong ageing — 30% intervention")
) |>
mutate(
n_people = round(n_people),
source = if_else(is.na(intervention) | intervention == "0% reduction",
"Baseline", "30% intervention")
) |>
select(source, scenario, year, state, n_people)
write_json(state_populations, "data/state_populations.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported state_populations.json\n")
## Exported state_populations.json
# ── F4: Crossover years ───────────────────────────────────────────────────────
crossover_export <- bind_rows(
# By intervention level
nh_scenarios |>
filter(over_capacity) |>
group_by(intervention) |>
summarise(crossover_year = min(year), type = "intervention", .groups = "drop"),
# By sensitivity scenario (baseline intervention only)
nh_sensitivity |>
filter(over_capacity) |>
group_by(sensitivity) |>
summarise(crossover_year = min(year), .groups = "drop") |>
rename(intervention = sensitivity) |>
mutate(type = "sensitivity")
)
write_json(crossover_export, "data/crossover_years.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported crossover_years.json\n")
## Exported crossover_years.json
# ── F5: Cost savings ──────────────────────────────────────────────────────────
cost_savings_export <- cost_savings |>
select(intervention, year,
total_cost_NOK, baseline_cost,
annual_saving_NOK, annual_saving_MNOK) |>
mutate(across(where(is.numeric), round))
write_json(cost_savings_export, "data/cost_savings.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported cost_savings.json\n")
## Exported cost_savings.json
# ── F6: Gap closure (fine-grained) ───────────────────────────────────────────
gap_closure_export <- crossover_fine |>
filter(!is.na(max_nh)) |>
mutate(
gap = max_nh - national_supply_flat,
baseline_gap = max_nh[strength == 0] - national_supply_flat,
gap_closed_pct = round((1 - gap / baseline_gap) * 100, 1),
gap = round(gap)
) |>
select(strength, max_nh, gap, gap_closed_pct, crossover)
write_json(gap_closure_export, "data/gap_closure.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported gap_closure.json\n")
## Exported gap_closure.json
# ── F7: Sensitivity results ───────────────────────────────────────────────────
sensitivity_export <- nh_sensitivity |>
select(sensitivity, year, nh_residents = n_people,
supply, gap, over_capacity) |>
mutate(across(where(is.numeric), round))
write_json(sensitivity_export, "data/sensitivity.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported sensitivity.json\n")
## Exported sensitivity.json
# ── F8: National population projections ──────────────────────────────────────
pop_export <- national_proj |>
mutate(population = round(population))
write_json(pop_export, "data/population_projections.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported population_projections.json\n")
## Exported population_projections.json
# ── F9: Municipality crossover ────────────────────────────────────────────────
# Check what columns exist after the join
muni_joined <- muni_action_needed |>
left_join(muni_lookup, by = "muni_code")
# Use whatever name column exists
name_col <- intersect(names(muni_joined), c("muni_name", "muni_name.y", "muni_name.x"))
muni_export <- muni_joined |>
mutate(across(where(is.numeric), ~ round(.x, 2))) |>
rename(muni_name = any_of(name_col)[1]) |>
select(muni_code, muni_name, year,
total_elderly, projected_demand, supply,
admissions_to_prevent, pct_reduction_needed,
use_rate_reduction_pp)
write_json(muni_export, "data/municipality_crossover.json",
pretty = TRUE, auto_unbox = TRUE)
cat("Exported municipality_crossover.json\n")
## Exported municipality_crossover.json
# ── F10: Summary confirmation ─────────────────────────────────────────────────
exported_files <- list.files("data", pattern = "\\.json$")
cat("\n=== EXPORT COMPLETE ===\n")
##
## === EXPORT COMPLETE ===
cat("Files written to data/ folder:\n")
## Files written to data/ folder:
for (f in exported_files) {
size <- file.size(file.path("data", f))
cat(sprintf(" %-40s %s KB\n", f, round(size / 1024, 1)))
}
## cost_savings.json 8.1 KB
## crossover_years.json 0.4 KB
## gap_closure.json 12.7 KB
## model_params.json 0.9 KB
## municipality_crossover.json 10.5 KB
## nh_trajectories.json 8.1 KB
## population_projections.json 6.7 KB
## sensitivity.json 4.7 KB
## state_populations.json 14.9 KB
cat("\nNext step: run 'quarto render dashboard.qmd' in the terminal.\n")
##
## Next step: run 'quarto render dashboard.qmd' in the terminal.
