Geographic distribution of healthcare facilities across Canada

Supply, access and inequity: Open Database of Healthcare Facilities linked to the 2021 Census

Objective

To describe how healthcare facilities are distributed across Canada, and to identify which populations have the weakest geographic access to hospitals.

Describe facility supply per 100,000 population by province and facility type.
Measure population-weighted geographic access to hospitals.
Quantify inequality in the distribution of facilities.
Test which community characteristics are associated with facility supply.
Identify spatial clusters of poor hospital access.
Assess residential care supply relative to the older population, and document data completeness.

Data and setup

Sources. The Open Database of Healthcare Facilities (ODHF, version 1.1, Statistics Canada) is downloaded automatically. Population, age structure, household income and community boundaries for 2021 census subdivisions (CSDs) come from the CensusMapper API through the cancensus package.

Before the first run, you need a free CensusMapper API key (create an account at censusmapper.ca, then open the API tab). Save it once on your computer by running this in R, then restart R:

cancensus::set_cancensus_api_key("YOUR_KEY", install = TRUE)

The first run takes several minutes. Results are cached in the _cache folder next to this file. Delete _cache/analysis_inputs.rds to force a full rebuild.

Code

# ---------- helper: find or download the ODHF csv ----------
# If several versions are present, the highest version number is used.
pick_latest <- function(paths) {
  nm <- basename(paths)
  v  <- suppressWarnings(as.numeric(sub(".*v([0-9]+\\.[0-9]+).*", "\\1", nm)))
  v[!grepl("v[0-9]+\\.[0-9]+", nm) | is.na(v)] <- 0
  paths[order(v, decreasing = TRUE)][1]
}

get_odhf <- function() {
  local <- list.files(".", pattern = "^odhf.*\\.csv$", ignore.case = TRUE, full.names = TRUE)
  if (length(local)) return(pick_latest(local))
  zipf <- "_cache/ODHF.zip"
  if (!file.exists(zipf))
    download.file("https://www150.statcan.gc.ca/n1/pub/13-26-0001/2020001/ODHF.zip",
                  zipf, mode = "wb", quiet = TRUE)
  unzip(zipf, exdir = "_cache/odhf", junkpaths = TRUE)
  csv  <- list.files("_cache/odhf", pattern = "\\.csv$", ignore.case = TRUE, full.names = TRUE)
  pref <- csv[grepl("odhf", basename(csv), ignore.case = TRUE)]
  pick_latest(if (length(pref)) pref else csv)
}

clean_csd_id <- function(x) {
  x <- trimws(as.character(x))
  x[x %in% c("", "NA", "NaN")] <- NA_character_
  x <- sub("\\.0+$", "", x)
  x[!is.na(x) & !grepl("^[0-9]{7}$", x)] <- NA_character_
  x
}

build_inputs <- function(odhf_path) {
  # ---- API key and cache for cancensus ----
  key <- Sys.getenv("CM_API_KEY")
  if (!nzchar(key)) key <- getOption("cancensus.api_key", "")
  if (!nzchar(key))
    stop("No CensusMapper API key found. Get a free key at censusmapper.ca, then run ",
         "cancensus::set_cancensus_api_key('YOUR_KEY', install = TRUE) and restart R.")
  options(cancensus.api_key = key)
  dir.create("_cache/cancensus", showWarnings = FALSE, recursive = TRUE)
  options(cancensus.cache_path = normalizePath("_cache/cancensus"))

  # ---- 1. ODHF ----
  odhf <- read_csv(odhf_path, col_types = cols(.default = "c"), show_col_types = FALSE)
  names(odhf) <- tolower(names(odhf))
  need <- c("facility_name", "odhf_facility_type", "provider", "province",
            "csduid", "latitude", "longitude")
  stopifnot(all(need %in% names(odhf)))
  odhf <- odhf |>
    mutate(latitude  = suppressWarnings(as.numeric(latitude)),
           longitude = suppressWarnings(as.numeric(longitude)),
           csduid    = clean_csd_id(csduid),
           province  = toupper(province),
           odhf_facility_type = str_squish(odhf_facility_type),
           row_id     = row_number(),
           has_coords = is.finite(latitude) & is.finite(longitude))

  # ---- 2. Census 2021 by community ----
  pop <- get_census(
    dataset = "CA21", regions = list(C = "01"), level = "CSD",
    vectors = c(pop_65plus = "v_CA21_251", median_income = "v_CA21_906"),
    use_cache = TRUE, quiet = TRUE) |>
    transmute(csduid = as.character(GeoUID),
              pruid  = as.character(PR_UID),
              name   = `Region Name`,
              population = Population,
              area_km2   = `Area (sq km)`,
              pop_65plus, median_income,
              density  = population / area_km2,
              share_65 = pop_65plus / population)

  # ---- 3. Community boundaries ----
  csd_sf <- get_census(dataset = "CA21", regions = list(C = "01"), level = "CSD",
                       geo_format = "sf", use_cache = TRUE, quiet = TRUE) |>
    transmute(csduid_geo = as.character(GeoUID)) |>
    st_transform(3347) |>
    st_make_valid()

  # ---- 4. Place each facility in a community using its coordinates ----
  with_xy <- odhf |>
    filter(has_coords) |>
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove = FALSE) |>
    st_transform(3347) |>
    st_join(csd_sf, join = st_intersects) |>
    st_drop_geometry() |>
    distinct(row_id, .keep_all = TRUE)

  odhf2 <- bind_rows(with_xy, filter(odhf, !has_coords)) |>
    mutate(csduid_final = coalesce(csduid_geo, csduid))

  # ---- 5. Distance from each community to its nearest hospital ----
  hosp <- odhf2 |>
    filter(odhf_facility_type == "Hospitals", has_coords) |>
    st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove = FALSE) |>
    st_transform(3347)

  cent <- st_sf(csduid   = csd_sf$csduid_geo,
                geometry = st_point_on_surface(st_geometry(csd_sf)))
  nearest <- st_nearest_feature(cent, hosp)
  cent$dist_km <- as.numeric(st_distance(cent, hosp[nearest, ], by_element = TRUE)) / 1000

  dat <- cent |>
    st_drop_geometry() |>
    left_join(pop, by = "csduid") |>
    filter(!is.na(population), population > 0)

  list(odhf2 = odhf2, pop = pop, csd_sf = csd_sf, dat = dat)
}

odhf_path  <- get_odhf()
cache_file <- file.path("_cache", paste0("inputs_", tools::file_path_sans_ext(basename(odhf_path)), ".rds"))
if (file.exists(cache_file)) {
  inp <- readRDS(cache_file)
} else {
  inp <- build_inputs(odhf_path)
  saveRDS(inp, cache_file)
}
odhf2 <- inp$odhf2; pop <- inp$pop; csd_sf <- inp$csd_sf; dat <- inp$dat

prov_lab <- c("10"="NL","11"="PE","12"="NS","13"="NB","24"="QC","35"="ON",
              "46"="MB","47"="SK","48"="AB","59"="BC","60"="YT","61"="NT","62"="NU")

# ---- community-level analysis table ----
counts <- odhf2 |>
  filter(!is.na(csduid_final)) |>
  mutate(type = case_when(
    odhf_facility_type == "Hospitals" ~ "hospitals",
    odhf_facility_type == "Ambulatory health care services" ~ "ambulatory",
    odhf_facility_type == "Nursing and residential care facilities" ~ "residential",
    TRUE ~ "other")) |>
  count(csduid_final, type) |>
  pivot_wider(names_from = type, values_from = n, names_prefix = "n_", values_fill = 0)
for (v in c("n_hospitals", "n_ambulatory", "n_residential", "n_other"))
  if (!v %in% names(counts)) counts[[v]] <- 0

d <- dat |>
  left_join(counts, by = c("csduid" = "csduid_final")) |>
  mutate(across(starts_with("n_"), ~ replace_na(.x, 0)),
         n_total = n_hospitals + n_ambulatory + n_residential + n_other,
         prov    = prov_lab[pruid],
         is_iri  = grepl("\\(IRI\\)", name),
         size_class = cut(population, c(0, 1000, 10000, 100000, Inf),
                          labels = c("<1,000", "1,000-9,999", "10,000-99,999", "100,000+"),
                          right = FALSE))

wmedian <- function(x, w) { o <- order(x); x <- x[o]; w <- w[o]
                            x[which(cumsum(w) >= sum(w) / 2)[1]] }
acc <- function(g) g |>
  summarise(communities = n(), pop_total = sum(population),
            median_km      = round(wmedian(dist_km, population), 1),
            pct_over_25km  = round(100 * sum(population[dist_km > 25])  / sum(population), 1),
            pct_over_50km  = round(100 * sum(population[dist_km > 50])  / sum(population), 1),
            pct_over_100km = round(100 * sum(population[dist_km > 100]) / sum(population), 1))
nat <- acc(d)

The ODHF file used was odhf_v1.csv. The analysis covers 9,039 listed facilities and 4,830 census subdivisions with a non-zero population (36,991,981 people).

Data completeness

ODHF is compiled from provincial and other sources that do not report identically, so counts are not strictly comparable across provinces. This table is the basis for reading every province comparison below.

Code

t_comp <- odhf2 |>
  group_by(province) |>
  summarise(facilities = n(),
            data_sources = n_distinct(provider),
            hospitals = sum(odhf_facility_type == "Hospitals", na.rm = TRUE),
            hospitals_with_coordinates = sum(odhf_facility_type == "Hospitals" & has_coords, na.rm = TRUE),
            ambulatory = sum(odhf_facility_type == "Ambulatory health care services", na.rm = TRUE),
            residential = sum(odhf_facility_type == "Nursing and residential care facilities", na.rm = TRUE),
            type_missing = sum(is.na(odhf_facility_type)),
            pct_no_coordinates = round(100 * mean(!has_coords), 1),
            pct_no_community_code = round(100 * mean(is.na(csduid_final)), 1)) |>
  save_tbl("table_completeness")
kable(t_comp)

Table 1: Data completeness by province (ODHF v1.1)

province	facilities	data_sources	hospitals	hospitals_with_coordinates	ambulatory	residential	type_missing	pct_no_coordinates	pct_no_community_code
AB	331	4	207	144	13	111	0	53.5	54.7
BC	1492	3	158	144	226	1107	1	13.7	13.7
MB	214	3	131	104	41	41	1	30.4	30.4
NB	62	4	58	53	2	2	0	8.1	9.7
NL	96	3	44	35	28	23	1	40.6	41.7
NS	195	3	55	46	5	135	0	5.1	5.1
NT	7	3	5	4	2	0	0	14.3	14.3
NU	3	2	3	3	0	0	0	0.0	0.0
ON	3733	3	520	400	1104	2104	5	15.4	15.6
PE	35	3	12	11	12	11	0	2.9	2.9
QC	2359	5	500	344	1114	742	3	12.5	12.6
SK	473	3	143	121	167	161	2	25.2	25.2
YT	39	2	17	13	13	7	0	15.4	15.4

Analysis 1: Facility supply per 100,000 population

Code

t1 <- d |>
  group_by(prov) |>
  summarise(communities = n(), pop_total = sum(population),
            across(c(n_hospitals, n_ambulatory, n_residential, n_total), sum)) |>
  mutate(across(starts_with("n_"), ~ round(.x / pop_total * 1e5, 1), .names = "{.col}_per100k")) |>
  select(prov, pop_total, ends_with("per100k")) |>
  save_tbl("table1_supply_per_100k")
kable(t1)

Table 2: Facilities per 100,000 population by province and type

prov	pop_total	n_hospitals_per100k	n_ambulatory_per100k	n_residential_per100k	n_total_per100k
AB	4262635	3.3	0.1	0.1	3.5
BC	5000879	2.9	4.5	18.3	25.7
MB	1342153	7.7	3.1	0.2	11.1
NB	775610	6.7	0.3	0.4	7.3
NL	510550	6.7	3.7	0.4	11.0
NS	969383	4.7	0.5	13.8	19.1
NT	41070	9.7	4.9	0.0	14.6
NU	36858	8.1	0.0	0.0	8.1
ON	14223942	2.8	7.5	11.7	22.1
PE	154331	7.1	7.8	7.1	22.0
QC	8501833	4.0	12.1	8.0	24.1
SK	1132505	10.7	13.0	7.4	31.3
YT	40232	32.3	32.3	12.4	82.0

Code

t1 |>
  select(prov, ends_with("per100k")) |>
  pivot_longer(-prov) |>
  mutate(name = str_remove_all(name, "^n_|_per100k$")) |>
  ggplot(aes(value, prov)) +
  geom_point(size = 2.5) +
  facet_wrap(~ name, scales = "free_x") +
  labs(x = "Facilities per 100,000 population", y = NULL) +
  theme_minimal()

Figure 1: Facilities per 100,000 population by province and type

Analysis 2: Geographic access to hospitals

Distance is measured in a straight line from a point inside each community to the nearest listed hospital. Percentages are weighted by population, so they describe people rather than communities. Nationally, 1.3% of the population lives more than 50 km from a hospital (median distance 3.9 km).

Code

t2 <- d |> group_by(prov) |> acc() |> arrange(desc(pct_over_50km)) |>
  save_tbl("table2_access_by_province")
kable(t2)

Table 3: Population-weighted access to hospitals by province or territory

prov	communities	pop_total	median_km	pct_over_25km	pct_over_50km	pct_over_100km
NU	26	36858	250.3	67.0	67.0	64.6
NT	37	41070	2.7	34.0	33.8	28.2
YT	31	40232	1.4	15.1	6.0	1.4
NL	365	510550	8.8	22.0	5.2	1.3
SK	899	1132505	1.4	15.6	4.1	0.9
MB	221	1342153	2.1	9.1	3.2	1.3
AB	406	4262635	5.7	7.9	3.0	0.1
NB	265	775610	8.6	18.6	1.7	0.0
BC	643	5000879	4.2	5.0	1.5	0.3
QC	1193	8501833	3.8	5.7	0.7	0.3
ON	554	14223942	3.3	1.8	0.4	0.2
NS	92	969383	18.1	8.5	0.0	0.0
PE	98	154331	5.1	5.7	0.0	0.0

Code

t3 <- d |> group_by(size_class) |> acc() |> save_tbl("table3_access_by_size")
kable(t3)

Table 4: Access to hospitals by community size

size_class	communities	pop_total	median_km	pct_over_25km	pct_over_50km	pct_over_100km
<1,000	2721	993216	24.8	49.3	14.7	5.7
1,000-9,999	1680	5227370	14.5	22.3	4.9	1.5
10,000-99,999	373	10134407	6.2	3.6	0.9	0.0
100,000+	56	20636988	2.6	0.0	0.0	0.0

Code

t4 <- d |> group_by(is_iri) |> acc() |> save_tbl("table4_access_reserves")
kable(t4)

Table 5: Access to hospitals: First Nations reserves (IRI) versus other communities

is_iri	communities	pop_total	median_km	pct_over_25km	pct_over_50km	pct_over_100km
FALSE	4054	36641351	3.8	5.1	1.1	0.2
TRUE	776	350630	21.9	47.0	27.1	12.5

Code

csd_sf |>
  left_join(select(d, csduid, dist_km), by = c("csduid_geo" = "csduid")) |>
  ggplot() +
  geom_sf(aes(fill = dist_km), color = NA) +
  scale_fill_viridis_c(trans = "sqrt", na.value = "grey90", name = "km") +
  theme_void()

Figure 2: Straight-line distance from each community to the nearest hospital (km, square-root scale)

Code

d |> group_by(prov) |> arrange(dist_km) |>
  mutate(cum_pop = cumsum(population) / sum(population)) |>
  ggplot(aes(dist_km + 1, cum_pop, colour = prov)) +
  geom_step() + scale_x_log10() +
  labs(x = "Distance to nearest hospital, km (log scale)",
       y = "Cumulative share of population", colour = NULL) +
  theme_minimal()

Figure 3: Cumulative share of population by distance to the nearest hospital

Analysis 3: Inequality in the distribution of facilities

Communities are ranked by facilities per person. The Gini coefficient measures how unevenly facilities are spread relative to population (0 is perfectly even, 1 is maximal inequality). Because hospitals are rare and sit in larger places, high values are expected for them.

Code

lorenz <- function(n, pop) { o <- order(n / pop)
  tibble(p = c(0, cumsum(pop[o]) / sum(pop)), f = c(0, cumsum(n[o]) / sum(n))) }
gini <- function(n, pop) { if (sum(n) == 0) return(NA_real_)
  l <- lorenz(n, pop); 1 - sum(diff(l$p) * (l$f[-1] + l$f[-nrow(l)])) }

t5 <- d |> group_by(prov) |>
  summarise(across(c(n_hospitals, n_ambulatory, n_residential, n_total),
                   ~ round(gini(.x, population), 3), .names = "gini_{.col}")) |>
  save_tbl("table5_gini")
kable(t5)

Table 6: Gini coefficient of facility distribution across communities, by province

prov	gini_n_hospitals	gini_n_ambulatory	gini_n_residential	gini_n_total
AB	0.803	0.996	0.874	0.800
BC	0.737	0.405	0.356	0.336
MB	0.824	0.927	0.624	0.811
NB	0.811	0.940	0.740	0.786
NL	0.843	0.954	0.784	0.811
NS	0.618	0.686	0.458	0.448
NT	0.522	0.932	NA	0.550
NU	0.773	NA	NA	0.773
ON	0.572	0.479	0.364	0.381
PE	0.774	0.836	0.719	0.759
QC	0.640	0.618	0.582	0.514
SK	0.811	0.766	0.805	0.691
YT	0.769	0.283	0.369	0.399

Code

map_dfr(c("n_hospitals", "n_ambulatory", "n_residential"),
        ~ lorenz(d[[.x]], d$population) |> mutate(type = str_remove(.x, "^n_"))) |>
  ggplot(aes(p, f, colour = type)) +
  geom_line(linewidth = 1) + geom_abline(linetype = "dashed") +
  labs(x = "Cumulative share of population (ranked by facilities per capita)",
       y = "Cumulative share of facilities", colour = NULL) +
  theme_minimal()

Figure 4: Lorenz curves of facility distribution across Canadian communities

Analysis 4: Community characteristics and facility supply

Negative binomial regression models the number of facilities in each community, with population as the offset, so coefficients are rate ratios for facilities per person. Province is adjusted for (Ontario is the reference). Communities with suppressed census values are excluded.

Code

dm <- d |>
  filter(!is.na(share_65), density > 0) |>
  mutate(prov = relevel(factor(prov), ref = "ON"),
         share_10pp   = share_65 * 10,
         log2_density = log2(density))

fit_nb <- function(y) MASS::glm.nb(
  as.formula(paste(y, "~ share_10pp + log2_density + is_iri + prov + offset(log(population))")),
  data = dm)

term_lab <- c(share_10pp = "Seniors share (per 10 percentage points)",
              log2_density = "Population density (per doubling)",
              is_iriTRUE = "First Nations reserve (vs other)")

t6 <- map_dfr(c("n_hospitals", "n_ambulatory", "n_residential"),
              ~ tidy(fit_nb(.x), exponentiate = TRUE) |>
                mutate(outcome = str_remove(.x, "^n_"))) |>
  filter(term %in% names(term_lab)) |>
  mutate(term = term_lab[term],
         lo = estimate * exp(-1.96 * std.error),
         hi = estimate * exp( 1.96 * std.error)) |>
  select(outcome, term, rate_ratio = estimate, lo, hi, p.value) |>
  mutate(across(c(rate_ratio, lo, hi), ~ round(.x, 3)),
         p.value = signif(p.value, 2)) |>
  save_tbl("table6_rate_ratios")
kable(t6)

Table 7: Rate ratios for facilities per capita (negative binomial regression, province-adjusted)

outcome	term	rate_ratio	lo	hi	p.value
hospitals	Seniors share (per 10 percentage points)	2.071	1.889	2.270	0.0e+00
hospitals	Population density (per doubling)	1.044	1.021	1.069	2.2e-04
hospitals	First Nations reserve (vs other)	0.239	0.098	0.581	1.6e-03
ambulatory	Seniors share (per 10 percentage points)	1.782	1.650	1.924	0.0e+00
ambulatory	Population density (per doubling)	0.978	0.959	0.998	3.0e-02
ambulatory	First Nations reserve (vs other)	3.923	2.932	5.250	0.0e+00
residential	Seniors share (per 10 percentage points)	1.806	1.672	1.952	0.0e+00
residential	Population density (per doubling)	1.087	1.065	1.109	0.0e+00
residential	First Nations reserve (vs other)	0.223	0.108	0.461	5.2e-05

Code

ggplot(t6, aes(rate_ratio, term, xmin = lo, xmax = hi, colour = outcome)) +
  geom_pointrange(position = position_dodge(.5)) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_log10() +
  labs(x = "Rate ratio (log scale)", y = NULL, colour = NULL) +
  theme_minimal()

Figure 5: Rate ratios (95% CI) for facilities per capita

Analysis 5: Spatial clustering of poor hospital access

Global Moran’s I tests whether communities with long distances sit next to each other more than chance would predict. Local indicators (LISA) locate the clusters. Distance varies smoothly across space, so a high global value is expected; the local map is the informative output.

Code

sp <- csd_sf |> inner_join(select(d, csduid, dist_km), by = c("csduid_geo" = "csduid"))
y  <- log1p(sp$dist_km)
nb <- poly2nb(sp, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

mi <- moran.test(y, lw, zero.policy = TRUE)
kable(tibble(morans_I = round(unname(mi$estimate[1]), 3),
             p_value  = signif(mi$p.value, 2),
             communities = length(y)))

Table 8: Global Moran’s I for log distance to nearest hospital

morans_I	p_value	communities
0.451	0	4830

Code

lisa <- localmoran(y, lw, zero.policy = TRUE)
z    <- as.numeric(scale(y))
lagz <- lag.listw(lw, z, zero.policy = TRUE)
sig  <- lisa[, 5] < 0.05
sp$cluster <- case_when(!sig ~ "Not significant",
                        z > 0 & lagz > 0 ~ "High-High (remote cluster)",
                        z < 0 & lagz < 0 ~ "Low-Low (well-served cluster)",
                        z > 0 & lagz < 0 ~ "High-Low",
                        TRUE ~ "Low-High")

ggplot(sp) + geom_sf(aes(fill = cluster), color = NA) +
  scale_fill_manual(values = c("Not significant" = "grey90",
                               "High-High (remote cluster)" = "firebrick",
                               "Low-Low (well-served cluster)" = "steelblue",
                               "High-Low" = "orange", "Low-High" = "lightblue")) +
  theme_void() + labs(fill = NULL)

Figure 6: Local spatial clusters of hospital access (LISA, p < 0.05)

Code

kable(count(st_drop_geometry(sp), cluster, name = "communities"))

Table 9: Number of communities in each cluster type

cluster	communities
High-High (remote cluster)	258
High-Low	7
Low-High	106
Low-Low (well-served cluster)	161
Not significant	4298

Analysis 6: Residential care supply relative to the older population

Counts of residential care facilities per 10,000 people aged 65 and over. Only communities where the census reports the senior population are included. ODHF has no bed counts, so this is a rough proxy for capacity.

Code

t7 <- d |>
  filter(!is.na(pop_65plus)) |>
  group_by(prov) |>
  summarise(seniors = sum(pop_65plus), residential_facilities = sum(n_residential)) |>
  mutate(per_10k_seniors = round(residential_facilities / seniors * 1e4, 2)) |>
  save_tbl("table7_residential_care")
kable(t7)

Table 10: Residential care facilities per 10,000 seniors by province

prov	seniors	residential_facilities	per_10k_seniors
AB	629075	4	0.06
BC	1016085	916	9.01
MB	228970	3	0.13
NB	177145	3	0.17
NL	120530	2	0.17
NS	215315	134	6.22
NT	4085	0	0.00
NU	1610	0	0.00
ON	2637615	1666	6.32
PE	32715	11	3.36
QC	1753290	680	3.88
SK	197525	84	4.25
YT	6005	5	8.33

Code

ggplot(t7, aes(seniors, per_10k_seniors, label = prov)) +
  geom_text() + scale_x_log10() +
  labs(x = "Population aged 65+ (log scale)", y = "Residential care facilities per 10,000 seniors") +
  theme_minimal()

Figure 7: Residential care facilities per 10,000 seniors against provincial senior population

Limitations

ODHF records the location and type of facilities, not beds, staff, services or wait times, so access here means geographic proximity only.
Distance is straight-line from a point in each community, so it understates travel in remote areas and where there are no roads. Road-network travel time would strengthen Analysis 2.
Provinces contribute to ODHF differently (see the completeness table), so province comparisons of counts can reflect reporting differences. Hospitals are the most comparable facility type.
The analysis is ecological: communities are the unit, so findings do not apply to individuals.
Census values are suppressed for small communities, mainly seniors’ share and income, which affects Analyses 4 and 6.
A community’s characteristics are linked to its facilities by census subdivision, which is large in rural areas and small in cities.

Reproducibility

Code

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Toronto
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] knitr_1.50      broom_1.0.8     spdep_1.3-8     spData_2.3.4   
 [5] cancensus_0.5.7 sf_1.0-21       lubridate_1.9.4 forcats_1.0.0  
 [9] stringr_1.5.1   dplyr_1.1.4     purrr_1.1.0     readr_2.1.5    
[13] tidyr_1.3.1     tibble_3.3.0    ggplot2_3.5.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.52          htmlwidgets_1.6.4  lattice_0.22-6    
 [5] tzdb_0.5.0         vctrs_0.6.5        tools_4.3.2        generics_0.1.4    
 [9] parallel_4.3.2     proxy_0.4-27       pkgconfig_2.0.3    KernSmooth_2.23-26
[13] RColorBrewer_1.1-3 lifecycle_1.0.4    compiler_4.3.2     farver_2.1.2      
[17] deldir_2.0-4       htmltools_0.5.8.1  class_7.3-23       yaml_2.3.10       
[21] crayon_1.5.3       pillar_1.11.0      MASS_7.3-60        classInt_0.4-11   
[25] wk_0.9.4           boot_1.3-31        tidyselect_1.2.1   digest_0.6.37     
[29] stringi_1.8.7      labeling_0.4.3     fastmap_1.2.0      grid_4.3.2        
[33] archive_1.1.10     cli_3.6.5          magrittr_2.0.3     dichromat_2.0-0.1 
[37] e1071_1.7-16       withr_3.0.2        scales_1.4.0       backports_1.5.0   
[41] bit64_4.6.0-1      sp_2.2-0           timechange_0.3.0   rmarkdown_2.29    
[45] igraph_2.1.4       bit_4.6.0          hms_1.1.3          evaluate_1.0.4    
[49] viridisLite_0.4.2  s2_1.1.7           rlang_1.1.6        Rcpp_1.1.0        
[53] glue_1.8.0         DBI_1.2.3          geojsonsf_2.0.3    vroom_1.6.5       
[57] rstudioapi_0.17.1  jsonlite_2.0.0     R6_2.6.1           units_0.8-7