CRC cases that occurred between 2012 and 2020 in Norway were included in the study (n = 38,968), with available data on staging (approximately 90% of the total data). CRC incidence rates were obtained from the Cancer Registry of Norway and corresponded to the yearly number of new cases of colorectal tumours observed at the municipality level. The municipality (kommune) is the lowest level administrative division in Norway. As the municipalities in Norway are undergoing continuous consolidation, data were adjusted based on the latest administrative changes from 2020, according to which Norway is divided into 356 administrative regions. Cases were drawn using the third edition of the International Classification of Diseases for Oncology codes for colon and rectal cancers (C18.0–18.9, C19.9, C26.0, C20.9). According to Surveillance, Epidemiology, and End Result Program stages were classified as localised (referred to as the early stage in this paper) and regional, and distant (referred to as the late stage).

To evaluate access to healthcare, we chose municipality-level measures of the availability of physicians’ resources for 2011 obtained from Statistics Norway (Municipality Healthcare Services). All physicians involved in the primary diagnostic and preoperative investigations were included in the study such as GPs, gastroenterologists, surgeons, diagnostic radiologists, and pathologists.

There were two primary outcome variables of interest in this study. These variables included: (1) CRC stage at diagnosis based on the SEER summary stage classified as late-stage (if regional or distant) and early stage (if localised) and (2) relative risk or standardised incidence ratio.

Municipality-level relative risk rates for late-stage CRC were calculated as the ratio between the number of observed and expected cases. Expected cases per municipality were equal to the expected counts of cases in the observed area after applying rates that are specific to Norway for the particular age and gender groups. Specifically, we calculated country-level disease rates for four age-sex strata (females 0–74 years, females 75 years and older, males 0–74 years, males 75 years and older) and applied them to each of the observed areas based on their age-sex specific population structure in the reported period. A relative risk of 1 indicated that the number of observed cancer cases in that municipality was equal to the number of cancer cases expected for the entire population. Relative risks greater than 1 indicated that more cases occurred than expected, and relative risks lower than 1 indicated that fewer cases occurred than expected.

Key exposure was the density of municipality-level physicians per 10,000 people. We expressed physician density as physician-per-population-ratio or number of physicians per 10,000 people at the municipality level. The physicians included GP counts per municipality, as well as specialist counts that were closely related to the CRC diagnostic and preoperative investigations, such as gastroenterologists, surgeons, diagnostic radiologists, and pathologists.

Other covariate data included socioeconomic, lifestyle, and urban–rural population information and additional available information, including colonoscopy rates. Municipality-level socioeconomic status was assessed based on the following variables: percentages of the population living in urban areas, the population living in households with persistent low income, and the population not employed or studying. Municipality-level lifestyle data included variables on smoking status and exercise levels in the population. Most of the data were obtained from Statistics Norway (2022) or the Norwegian Institute for Public Health (2022). In addition, number of colonoscopies per 10,000 people was obtained from the Gastronet project (a network for quality assurance of gastroenterology in Norway) at the Cancer Registry of Norway.

The main goal of our analysis was to estimate the causal effect of physician density on late-stage CRC distribution. We considered two analyses to investigate the relationship of interest in the presence of unmeasured spatial confounders [12]. We investigated this relationship by applying a novel causal inference method for spatial data, that is, a neighbourhood adjustment model via spatial smoothing, which we referred to as the NA approach. We compared results from the NA approach to results from the more established Bayesian hierarchical BYM model for count data, which we referred to as the BYM approach. This approach only estimates correlations [14].In the BYM approach, we modelled late-stage CRC incidence as a Poisson-distributed outcome using a log link with exposure (physician density) and other available covariates (Table 1). This approach included two random effects: a spatially structured random effect, which smooths locally towards the values of nearby areas, and an unstructured random effect that smooths globally, towards the overall average. We assumed that spatially correlated random effects follow the CAR distribution. The CAR prior was defined by an adjacency matrix, where geographically adjacent neighbours were given a value of 1, and all other pairs, a value of 0.

However, Paciorek [15] found that simply accounting for spatial correlation may not fully resolve spatial confounding and that commonly used estimators, including the BYM approach, can be biased in the presence of unmeasured spatial confounding. As a response to these findings, Schnell and Papadogeorgou [12] proposed the NA method that adjusts for unobserved spatial confounders by blocking their statistical dependence on either the treatment or the exposure variable. Specifically, in the NA approach, we addressed confounding bias by joint spatial modelling of exposure and unobserved spatial confounders. The bias, denoted as \(B(X)\), is a result of attributing the effect of the confounder on the outcome to the exposure variable when exposure and spatial confounders (U) are correlated [10]. The bias was modelled by specifying a joint distribution of spatial confounders and exposure that allowed for different ranges of spatial correlation and permitted a correlation between U and exposure. To identify the unmeasured confounding bias, Schnell and Papadogeorgou (12) provided a set of assumptions, including the main assumption that the spatial scale of exposure is larger or about the same as the unmeasured confounder (see Additional file 1 for more details).

Both approaches were fit within the Bayesian paradigm, and Gaussian priors were used for regression coefficients with a mean of 0 and standard deviation of 10. Within the Bayesian paradigm, the spatial confounder U was viewed as a missing variable that was iteratively imputed using a Gibbs sampler. Posterior distributions were drawn using 10,000 Gibbs sampler iterations after 1000 burn-in iterations.

The main aim of our statistical analysis was to quantify the causal effect of a one percentage point increase in physician density on late-stage colorectal cancer distribution. We used the most common estimand for continuous treatments, the population average exposure–response curve. The population average exposure–response curve can be interpreted as the posterior distribution of exponentiated exposure coefficient, that is, the relative expected risk of late-stage CRC due to a one percentage point increase in physician density in a randomly chosen municipality. A posterior probability of relative risk lower than 1 indicates a negative association, that is, a protective effect of higher physician density and CRC outcome. The posterior probability of relative risk higher than 1 indicates a positive association, that is, a harmful effect of higher physician density on CRC outcome.

We used 95% credible intervals to assess statistical significance. The statistical analyses were performed in R software, and the codes for the BYM and NA approaches available from Schnell & Papadogeorgou’s GitHub site (https://github.com/schnellp/causal-spatial) were used and adjusted accordingly Choropleth maps (Fig. 1a and b) were made in Python software (mapclassify package) by using the Fisher-Jenks classification method, to determine optimal classes for visual data presentation.

a. Municipality-level relative risk rates for late-stage CRC, 2012–2020, Norway (The risk rate is defined as the ratio between the number of observed and expected CRC cases in a municipality). b. Municipality-level physician density per 100,000 people, 2011, Norway (Physician density is defined as the number of physicians per 100,000 inhabitants in a municipality)