Skip to main content
  • Research article
  • Open access
  • Published:

Geographic disparities in Saskatchewan prostate cancer incidence and its association with physician density: analysis using Bayesian models

Abstract

Background

Saskatchewan has one of the highest incidence of prostate cancer (PCa) in Canada. This study assesses if geographic factors in Saskatchewan, including location of where patients live and physician density are affecting the PCa incidence. First, the objective of this study is to estimate the PCa standardized incidence ratio (SIRs) in Saskatchewan stratified by PCa risk-level. Second, this study identifies clusters of higher than and lower than expected PCa SIRs in Saskatchewan. Lastly, this study identifies the association (if any) between family physician density and estimated PCa SIRs in Saskatchewan.

Methods

First, using Global Moran’s I, Local Moran’s I, and the Kuldorff’s Spatial Scan Statistic, the study identifies clusters of PCa stratified by risk-levels. Then this study estimates the SIRs of PCa and its association with family physician density in Saskatchewan using the Besag, York, and Mollie (BYM) Bayesian method.

Results

Higher than expected clusters of crude estimated SIR for metastatic PCa were identified in north-east Saskatchewan and lower than expected clusters were identified in south-east Saskatchewan. Areas in north-west Saskatchewan have lower than expected crude estimated SIRs for both intermediate-risk and low-risk PCa. Family physician density was negatively associated with SIRs of metastatic PCa (IRR: 0.935 [CrI: 0.880 to 0.998]) and SIRs of high-risk PCa (IRR: 0.927 [CrI: 0.880 to 0.975]).

Conclusions

This study identifies the geographical disparities in risk-stratified PCa incidence in Saskatchewan. The study identifies areas with a lower family physician density have a higher-than-expected incidences of metastatic and high-risk PCa. Hence policies to increase the number of physicians should ensure an equitable geographic distribution of primary care physicians to support early detection of diseases, including PCa.

Peer Review reports

Background

Prostate cancer (PCa) accounts for about 20% of all new cancer cases among men in Canada [1]. Within Canada, the incidence rate of PCa varies between provinces. In 2019, Saskatchewan had the third highest projected age-standardized PCa incidence rate (117.8 cases per 100,000 in 2019) when compared to other Canadian provinces [1]. In addition, the Saskatchewan age-standardized PCa incidence rates have remained higher than the national Canadian rates for the majority of the past 10 years [1, 2]. Previous studies have shown geographic factors influence PCa outcomes in Saskatchewan [3, 4], hence this study explores the influence of geographic patterns on PCa incidence rates in Saskatchewan.

Saskatchewan, in terms of geography, has the second lowest population density in Canada (after Newfoundland), with a majority of province sparsely populated and nearly 40% of the Saskatchewan population living in rural areas [5]. Because cancer patient outcomes are worse for rural dwellers compared to urban dwellers [6,7,8,9,10] and Saskatchewan has a relatively large rural population [5], it is possible that the rates for different PCa risk-levels are associated with the geographic distribution of Saskatchewan residents.

The low population density of Saskatchewan results in the geographic factors of remoteness and commute time being healthcare access barriers for Saskatchewan residents [3, 4, 11, 12]. The low population density also impacts the distribution of physicians in the province. Saskatchewan has one of the lowest per capita physician supplies (also known as physician density) compared to the other provinces in Canada (190.3 per 100,000 in 2014 and 204.5 per 100,000 in 2018) [13]. Because one mechanism for improving health outcomes, including reductions in PCa-specific mortalities, is increasing physician supply [14,15,16,17,18,19,20], understanding the association (if any) between physician density and PCa risk-level incidence is crucial to improving PCA outcomes in Saskatchewan.

While the factors leading to such high Saskatchewan PCa incidence rates are unknown, we hypothesize that the unique geography of Saskatchewan may be contributing to the high incidence of PCa in Saskatchewan. In this study we explore the geographic distribution of PCa cases in Saskatchewan. In addition, since the incidence for advanced cancers is known to decrease with the increase in availability of physicians [18, 21, 22], we identify the association (if any) that exists between the family physician density and PCa standardized incidence ratios (SIRs) in Saskatchewan.

The first study objective is to estimate the PCa SIRs in Saskatchewan stratified by PCa risk-levels. The second objective is to identify clusters of higher than and lower than expected PCa SIR stratified by PCa risk-levels in Saskatchewan. The final objective is to identify any associations between family physician density and estimated PCa SIRs in Saskatchewan.

Methods

Data and study area

The data for PCa were from the Saskatchewan Cancer Registry (SCR) and consisted of demographic, clinical, and geographic information for 3526 patients diagnosed with PCa between 2010 and 2014. Based on the demographic information, all PCa patients were age 35 years or over. The study area contained 82 geographic areas (GAs) in central and southern Saskatchewan categorized (for privacy reasons) by SCR using residence codes (Figs. 1 and 2) [23]. From 2010 to 2014, the study area contained 3289 PCa patients, after excluding those living out-of-province at the time of diagnosis (194 patients) and those (43 patients) living in the three northern regions (Mamawetan Churchill River, Keewatin Yatthe, and Athabasca). The northern regions were excluded because these regions could not be subdivided due to privacy reasons. Of these 3289 PCa patients, the analysis further excluded 298 patients because their PCa risk levels were unknown. Therefore, the final sample had 2991 patients, each categorized per the GA in which the patient lived at the time of diagnosis.

Fig. 1
figure 1

Crude estimated SIRs for metastatic, high-risk, intermediate-risk, and low-risk PCa cases in Saskatchewan (2010–2014)

Fig. 2
figure 2

Metastatic PCa crude estimated SIR clustering analysis: (A) Local Moran’s I; (B) Kuldorff’s Spatial Scan Statistic

To calculate SIRs (described in the following sections), population counts for 2012 from the Saskatchewan Covered Population (SCP) were used in the denominator in the formula to calculate PCa SIR (see Definitions section) [23]. The SCP is a count of residents with provincial health insurance in Saskatchewan in a given year and is maintained by Government of Saskatchewan [23]. Because all PCa patients in the SCR dataset were age 35 or over, the SCP data used was for men over the age of 35. The overall SCP for men age 35 or over deviated less than 3% each year [24], therefore we chose to use statistics from the midpoint year 2012 (between 2010 and 2014) for the denominators in the calculations [25].

To calculate physician density (described in the following sections) for the period 2010 to 2014, the required data from the Canadian Medical Association was only available for 2011 [26]. Hence our estimated physician densities are based on the year 2011. Canadian Medical Association data consists of family physicians and general practitioners licensed to practice medicine in Saskatchewan.

The relative variability of PCa incidence between the time period (2010 to 2014) was assessed using the coefficient of variation (ratio of the standard deviation to the average).

The University of Saskatchewan BioMedical Research Ethics Board provided ethics approval (Bio-REB certificate #15–34).

Definitions

The risk levels (low, intermediate, high) for PCa were based on the Genitourinary Radiation Oncologists of Canada (GUROC) definitions [27] and a fourth risk level (“metastatic”) was added to include patients diagnosed with metastatic cancer. For each risk level, the expected number of PCa cases in the ith GA (Ei) was calculated as follows [28]:

$$ {E}_i={n}_i\left(\frac{\sum_i{O}_i}{\sum_i{n}_i}\right) $$

where ni and Oi respectively denote the population count of men age 35 or over and the observed number of PCa cases in the ith GA. For each PCa risk level, the standardized incidence ratio (SIR) (which will be referred to as the crude estimated SIR) was estimated by dividing the number of observed cases in by the number of expected cases in each GA [28].

Independent variables

The study variables of interest were physician density, remoteness index, and closest PCa assessment centre. For each GA, the family physician density was calculated using the 2011 Canadian Medical Association data and the same population denominator used for the expected count of PCa cases. A remoteness index for a GA was calculated using the average of the Statistics Canada remote indices for regions forming the GA [29]. For each GA, the closest PCa assessment centre was categorized as Regina or Saskatoon, based on the shortest Euclidean distance between the centroid of the GA and the centroids of Saskatoon and Regina. Further details regarding remoteness index and closest PCa assessment centre variables used in this study can be found in the literature [3, 4].

Statistical methods

Clustering analysis was conducted to identify spatial clusters of PCa SIR by each risk level. Second, for each PCa risk level, a null model was built where the crude estimated SIRs were smoothed using the method proposed by Besag, York and Mollie (BYM model) [29]. The estimated SIRs from the BYM models will be referred to as the smoothed estimated SIRs. Third, ecological analyses were conducted to assess associations between the independent variables and the smoothed estimated SIRs for each of four PCa risk levels.

Clustering analysis

For each risk level, Global Moran’s I was calculated using the crude estimated SIR value for each GA [30]. The statistical significance for Global Moran’s I statistic was calculated using 999 permutations [30], which, if significant, demonstrates that the GAs sharing common boundaries have similar SIRs instead of having random geographically-distributed SIRs [31]. GAs within a 120-km radius of a GA were identified as neighbours of the GA. The corresponding weight matrix for the analysis was then computed using the inverse of the Euclidean distances between the centroids of a GA and its neighbours. This weight matrix was chosen to reflect the suspected correlation structure of the data [32].

For each risk level with statistically significant Global Moran’s I values, the crude estimated SIRs were studied further using the Local Moran’s I and Kuldorff’s Spatial Scan statistics [33,34,35].

BYM modeling

The SIRs were estimated using a Bayesian model-based approach to ensure, if spatial correlation exists, the estimated SIRs (i.e the smoothed estimated SIRs) were corrected for any spatial dependence between the GAs.

First, for each PCa risk level, a null model was built where the smoothed estimated SIRs were computed using the Bayesian BYM method [29]. Due to the count nature of the data, we assume our observed data Oi follows a Poisson distribution [36] with mean Eiθi where Ei and θi respectively denote the expected number of PCa cases and the “true” SIR in the ith GA [37, 38]. The BYM method models the log of the SIR as follows:

$$ \mathrm{Log}\ \left({\uptheta}_{\mathrm{i}}\right)=\mathrm{c}+{\mathrm{u}}_{\mathrm{i}}+{\mathrm{v}}_{\mathrm{i}} $$

where intercept c is the mean, and the terms ui and vi respectively denote the spatially structured and unstructured random effects [37, 38].

The parameters used in this model are based on literature [37,38,39,40,41]. The random effects and the intercept are assigned prior distributions. The intercept was assigned a uniform prior that extends over the whole real line [37, 38]. The structured random effect ui was assumed to follow a conditional auto-regressive distribution and the unstructured random effect vi was assumed to follow a normal distribution with mean zero [37, 38]. The variability for both random effects were controlled by a precision parameter. The precision parameter for the random effects were assigned a Gamma distribution with hyper-prior specification of (0.5, 0.0005) [39, 41].

The simulation for each model consisted of three chains [42, 43]. Each chain consisted of 200,000,000 iterations to obtain 50,000 data points: one for each 4000 time steps taken. A burn-in period of 8,000,000 iterations was selected based on the characteristics of the Brooks-Gelman-Rubin plots [38, 42, 43]. To determine whether the generated estimates for each parameter were from the correct distribution, the following diagnostic tests were performed: potential scale reduction factor, [42] stationarity and half-width tests, [44] Z-score for equality of the means, [45] and run length control [46, 47]).

Ecological analysis

Using the BYM models, unconditional analyses were conducted to identify any associations between the independent variables and the SIRs for each risk level. The statistical significance of an independent variable was determined via its 95% credible interval (CrI).

Global and Local Moran’s I statistics were computed using Geoda 1.12 [48]. Kuldorff’s Spatial Scan Statistic was computed using SatScan™ v9.4 [49]. SIR maps were built using quantum Geographical Analysis System (QGIS.org) Version 3.12 [50]. BYM models were built in OpenBUGS version 3.2.3 [51]. Convergence diagnostics for the BYM models were conducted in R using the package ‘coda’ [52].

Results

The study sample consisted of an average count of 598 PCa cases per year between 2010 and 2014, and the coefficient of variation for PCa incidence between 2010 and 2014 was 7.9%. During the five-year period, the coefficient of variation for remoteness index and physician density were 1.4 and 5.0%, respectively. Based on the age demographic information of all cases, a majority of the PCa cases were 70 years or older, followed by those who were 60 to 69 years old (Table 1). However, the distribution of the age demographics varied by risk level. Low and intermediate risk PCa cases had higher proportions of cases in the younger age groups. In contrast, high risk and metastatic cases had higher proportions of cases in the older age groups. Among all cases, each year (between 2010 and 2014) the proportion of cases diagnosed was about 20% with deviations of less than 2%. See Table 1 for details.

Table 1 Demographic information of the PCa cases stratified by GUROC risk levels (n = 2991)

The highest proportion of cases were high-risk PCa (36.4%) followed by intermediate-risk (35.3%), low-risk (14.8%) and metastatic cases (13.4%). In nearly a third of GAs (32.9%), the observed incidence of metastatic PCa was more than 50% than the expected incidence. In 28, 18 and 24% of GAs, the observed incidences of high-risk, intermediate-risk and low-risk PCa, respectively, were more than 50% than the expected incidence. See Table 2 for details.

Table 2 Crude estimated SIRs for PCa cases diagnosed within geographic areas by GUROC risk level

Clustering analysis

The pattern of crude estimated SIRs for each PCa risk level in Saskatchewan is visualized in Fig. 1. Spatial patterns within Fig. 1 are identified using clustering analysis. The Global Moran’s I statistics for the crude estimated SIRs for each PCa risk level (except for high-risk) show evidence of positive spatial autocorrelation (Table 3). Hence, there was evidence that some geographical areas in Saskatchewan sharing boundaries have similar crude estimated SIRs for metastatic, intermediate and low risk PCa, instead of a random distribution of incidence patterns.

Table 3 Result of the Global and Local Moran’s I for each of GUROC risk level

Using the Local Moran’s I statistic, clusters of crude estimated SIRs for each PCa risk level were identified. In Fig. 2, “high-high” clusters of metastatic PCa were identified in the north-east part of the study area. Hence areas in north-east Saskatchewan have higher-than-average crude estimated SIRs for metastatic PCa. For intermediate-risk and low-risk PCa, “low-low” clusters are identified in the north-west part of the study area. Therefore, areas in north-west Saskatchewan have lower-than-average crude estimated SIRs for both intermediate-risk and low-risk PCa (Fig. 3 and Fig. 4).

Fig. 3
figure 3

Intermediate-risk PCa crude estimated SIR clustering analysis: (A) Local Moran’s I; (B) Spatial Scan Statistic

Fig. 4
figure 4

Low-risk PCa crude estimated SIR clustering analysis: (A) Local Moran’s I; (B) Spatial Scan Statistic

For Kuldorff’s Spatial Scan Statistic, the maximum spatial window size for metastatic, intermediate and low risk PCa were equal to or less than 30, 25 and 25%, respectively, of the total population. Kuldorff’s Spatial Scan Statistic identified a higher-than-the-average cluster of crude estimated SIRs for metastatic PCa in north-east Saskatchewan and lower-than-the-average cluster in south-east Saskatchewan, analogous to the clusters identified using the Local Moran’s I statistics (Fig. 2). Similarly, the spatial scan statistic results for intermediate-risk and low-risk PCa were comparable to the clusters identified using the Local Moran’s I statistics described earlier (Fig. 3 and Fig. 4).

BYM modeling

The crude and smoothed estimated SIRs for a GA are illustrated in Fig. 5. For both metastatic and high-risk PCa, the smoothed BYM estimates highlight areas of elevated incidence in north-east part of Saskatchewan (Fig. 5). Also in Fig. 5, the smoothed estimated SIRs for intermediate-risk and low-risk PCa identify areas of low incidence in north-west part of Saskatchewan. Table 4 illustrates how the crude estimated minimum and maximum SIR values are adjusted by the BYM modelling.

Fig. 5
figure 5

Quantile distribution of PCa crude estimated SIR and smoothed estimated SIR by GUROC risk levels

Table 4 Comparison of minimum/maximum values of crude and smooth estimated SIRs for each GUROC risk level

Ecological analysis

Family physician density was negatively associated with the smoothed estimated SIRs for metastatic PCa (IRR: 0.935 [CrI: 0.880 to 0.998]) and for high-risk PCa (IRR: 0.927 [CrI: 0.880 to 0.975]). Based on the mean coefficient of family physician density for metastatic PCa (Table 5), one unit increase (or increase of 1 physician per 1000 population) would be equal to an average decrease in metastatic PCa SIR by 6.5%. Similarly, an average increase of 1 physician per 1000 population would be equal to an average decrease in high-risk SIR by 7.3%. Figure 6 provides geographic pattern of family physician density in Saskatchewan and comparison with Fig. 5 visually compliments the negative correlation with metastatic PCa and high-risk PCa. For intermediate-risk and low-risk PCa, based on the credible intervals, there was no evidence of association with family physician density (Table 5).

Table 5 Result of the ecological analysis using Bayesian BYM analysis for each of GUROC risk level
Fig. 6
figure 6

Quantile distribution of family physician density in Saskatchewan in 2011

There was no evidence of any association between the SIR of each PCa risk levels and the two remaining independent variables (closest PCa assessment centre and remoteness index).

Discussion

This study estimated risk stratified PCa SIRs in Saskatchewan to identify if any geographic patterns and disparities. The geographic patterns of the risk stratified SIRs identified areas of concern (higher than expected SIRs) in Saskatchewan using Bayesian models and traditional clustering analysis methods. This study found clustering of higher than expected incidence for metastatic PCa in north-east part of Saskatchewan, and lower than expected incidence in south-east part of Saskatchewan. This study also identified lower than expected incidence of intermediate-risk and low-risk PCa in the north-west part of Saskatchewan. The estimation of SIRs using the BYM method led to adjustment of the crude estimated SIRs to facilitate identification of spatial trends [53].

Our study also shows that areas with lower density of family physician have higher than expected incidence of metastatic and high-risk PCa. A similar trend has been observed in the United States where increases in primary care physician density were associated with a decrease in late-stage diagnosis of cancers including PCa [21, 54]. The findings of this study highlight the effect that increasing physician supply may have on improving health outcomes, as identified in previous studies including PCa [14,15,16,17,18,19,20]. The results also highlight the wide-ranging distribution of family physician within Saskatchewan, acknowledging Saskatchewan also has one of the lowest per capita physician supply compared to the other provinces in Canada [13]. Hence, policies to increase physician supply should ensure equitable geographic distribution of primary care physicians to support early detection of diseases including PCa.

Nearly half of the patients in the sample were high risk (36.4%) or metastatic (13.4%), which could be indicative of physician practices including receptiveness towards PCa screening policies. Literature shows that physician beliefs regarding PCa screening/diagnosis procedures can influence their practice (physicians who are uncertain about PCa screening/diagnosis procedures are less receptive towards offering PCa screening/diagnosis to their patients). Because such practice variations exist among Saskatchewan physicians [55] and the literature reports that an increase in advanced PCa may be due to a decrease in PCa screening [56], physician beliefs might possibly explain the geographic variations in PCa diagnosis rates.

Although family physician density was not associated with diagnostic pattern for low-risk and intermediate-risk PCa, further research is needed if these regional trends are related to physician practices given the controversy of screening tests for early detection of PCa [57, 58]. Given recent research showing PCa screening and detection of early-stage PCa decreasing, potentially due to mixed PCa screening guidelines, further studies assessing the role of PCa screening guidelines on geographic disparities in early-stage PCa incidence may provide further explanation [59].

This study identifies the geographical disparities in risk-stratified PCa incidence in Saskatchewan. This study suggests that healthcare access factors [60], including availability of physicians and the geographic location of individuals, may affect health outcomes for PCa. This study further highlights the possibility that enhancing health delivery in rural areas may improve health outcomes. A recent report by the Rural Road Map Implementation Committee in Canada shows there are continued challenges regarding healthcare access in rural parts of Canada including difficulties of attracting and retaining physicians [61].

The limitations of this study include the use of aggregate data for the ecological study design due to lack of information on individual-level data on family physician availability to the patient. However, the study uses widely developed Bayesian and conventional spatial analysis methods to identify inherent patterns in the study area. Other limitations include the unavailability of socioeconomic/demographic data and PCa screening data for the study period.

Conclusions

This study identified geographic disparities in PCa incidence in Saskatchewan. There were higher than expected incidence of metastatic PCa in north-east parts of Saskatchewan, and lower than expected incidence of intermediate-risk and low-risk PCa in the north-west part of Saskatchewan. In addition, areas with lower density of family physician had higher than expected incidence of metastatic and high-risk PCa. This study shows that availability of community level healthcare providers and geographic location of patients affects cancer care in Saskatchewan. This highlights the need for adequate availability of primary care physicians in rural and urban areas to improve cancer care in Saskatchewan.

Availability of data and materials

This study was conducted secondary analysis of existing data, which is only available at Saskatchewan Cancer Registry. For further information on accessing data, please contact Saskatchewan Cancer Registry at info@saskcancer.ca.

Abbreviations

PCa:

prostate cancer

SIRs:

standardized incidence ratios

BYM:

Besag, York and Mollie

GA:

Geographic Areas

IRR:

Incidence Rate Ratio

CrI:

Credible Interval

SCR:

Saskatchewan Cancer Registry

SCP:

Saskatchewan Covered Population

GUROC:

Genitourinary Radiation Oncologists of Canada

References

  1. Committee CCSA. Canadian Cancer Statistics 2019. Toronto; 2019 Sept 2019.

  2. Table 13-10-0762-01 Number of new cases and age-standardized rates of primary cancer, by stage at diagnosis, selected cancer type and sex: Statistics Canada; 2020.

  3. Andkhoie M, Szafron M. The impact of geographic location on Saskatchewan prostate Cancer patient treatment choices: a multilevel and spatial analysis. J Rural Health 2020 DOI: https://doi.org/10.1111/jrh.12471, 36, 4, 564, 576.

  4. Andkhoie M, Szafron M. Geographic factors associated with time-to-treatment outcomes for radiation therapy among localized prostate cancer patients in Saskatchewan. J Cancer Policy. 2020;26:100259. https://doi.org/10.1016/j.jcpo.2020.100259.

    Article  Google Scholar 

  5. Moazzami B. Fewer & older: population and demographic crossroads in rural Saskatchewan. Canada: Strengthening Rural Canada; 2015.

  6. Laing KA, Bramwell SP, McNeill A, Corr BD, Lam TB. Prostate cancer in Scotland: does geography matter? An analysis of incidence, disease characteristics and survival between urban and rural areas. J Clin Urol. 2014;7(3):176–84. https://doi.org/10.1177/2051415813512303.

    Article  Google Scholar 

  7. Carriere R, Adam R, Fielding S, Barlas R, Ong Y, Murchie P. Rural dwellers are less likely to survive cancer – an international review and meta-analysis. Health & lace. 2018;53:219–27. https://doi.org/10.1016/j.healthplace.2018.08.010.

    Article  Google Scholar 

  8. Unger JM, Moseley A, Symington B, Chavez-MacGregor M, Ramsey SD, Hershman DL. Geographic Distribution and Survival Outcomes for Rural Patients With Cancer Treated in Clinical Trials. JAMA Network Open. 2018;1(4):e181235-e DOI: https://doi.org/10.1001/jamanetworkopen.2018.1235.

  9. Papa N, Lawrentschuk N, Muller D, MacInnis R, Ta A, Severi G, et al. Rural residency and prostate cancer specific mortality: results from the Victorian Radical Prostatectomy Register. Australian and New Zealand Journal of Public Health. 2014;38(5):449–54 DOI: https://doi.org/10.1111/1753-6405.12210.

  10. Afshar N, English DR, Milne RL. Rural-urban residence and cancer survival in high-income countries: a systematic review. Cancer. 2019;125(13):2172–84. https://doi.org/10.1002/cncr.32073.

    Article  PubMed  Google Scholar 

  11. Nair BV, Schuler R, Stewart S, Taylor-Gjevre RM. Self-reported barriers to healthcare access for rheumatoid arthritis patients in rural and northern Saskatchewan: a mixed methods study. Musculoskeletal Care. 2016;14(4):243–51. https://doi.org/10.1002/msc.1146.

    Article  CAS  PubMed  Google Scholar 

  12. Karunanayake CP, Rennie DC, Hagel L, Lawson J, Janzen B, Pickett W, et al. Access to specialist Care in Rural Saskatchewan: the Saskatchewan rural health study. Healthcare (Basel). 2015;3(1):84–99. https://doi.org/10.3390/healthcare3010084.

    Article  Google Scholar 

  13. Physician in Canada, 2018. Ottawa: Canadian Institute for Health Information.; 2019.

  14. Macinko J, Starfield B, Shi L. Quantifying the health benefits of primary care physician supply in the United States. Int J Health Serv. 2007;37(1):111–26. https://doi.org/10.2190/3431-G6T7-37M8-P224.

    Article  PubMed  Google Scholar 

  15. Coughlin SS, Leadbetter S, Richards T, Sabatino SA. Contextual analysis of breast and cervical cancer screening and factors associated with health care access among United States women, 2002. Soc Sci Med. 2008;66(2):260–75. https://doi.org/10.1016/j.socscimed.2007.09.009.

    Article  PubMed  Google Scholar 

  16. Ferrante JM, Gonzalez EC, Pal N, Roetzheim RG. Effects of physician supply on early detection of breast cancer. J Am Board Fam Pract. 2000;13(6):408–14. https://doi.org/10.3122/15572625-13-6-408.

    Article  CAS  PubMed  Google Scholar 

  17. Fleisher JM, Lou JQ, Farrell M. Relationship between physician supply and breast cancer survival: a geographic approach. J Community Health. 2008;33(4):179–82. https://doi.org/10.1007/s10900-008-9090-z.

    Article  PubMed  Google Scholar 

  18. Gorey KM, Luginaah IN, Holowaty EJ, Fung KY, Hamm C. Associations of physician supplies with breast cancer stage at diagnosis and survival in Ontario, 1988 to 2006. Cancer. 2009;115(15):3563–70. https://doi.org/10.1002/cncr.24401.

    Article  PubMed  Google Scholar 

  19. Yao N, Foltz SM, Odisho AY, Wheeler DC. Geographic analysis of urologist density and prostate Cancer mortality in the United States. PLoS One. 2015;10(6):e0131578. https://doi.org/10.1371/journal.pone.0131578.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Aneja S, Yu JB. The impact of county-level radiation oncologist density on prostate cancer mortality in the United States. Prostate Cancer Prostatic Dis. 2012;15(4):391–6. https://doi.org/10.1038/pcan.2012.28.

    Article  CAS  PubMed  Google Scholar 

  21. Nguyen KD, Hyder ZZ, Shaw MD, Maness SB, Cookson MS, Patel SG, et al. Effects of primary care physician density, urologist presence, and insurance status on stage of diagnosis for urologic malignancies. Cancer Epidemiol. 2018;52:10–4. https://doi.org/10.1016/j.canep.2017.10.012.

    Article  PubMed  Google Scholar 

  22. Ananthakrishnan AN, Hoffmann RG, Saeian K. Higher physician density is associated with lower incidence of late-stage colorectal cancer. J Gen Intern Med. 2010;25(11):1164–71. https://doi.org/10.1007/s11606-010-1457-z.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Saskatchewan e. Covered Population 2015. Saskatchewan: eHealth Saskatchewan; 2015 June 30, 2015.

  24. Saskatchewan e. Saskatchewan Health Coverage Reports Saskatchewan: eHealth Saskatchewan; 2020 [Available from: https://opendata.ehealthsask.ca/MicroStrategyPublic/asp/Main.aspx.

  25. Crombie IK, Cramer N. Five-year age-specific incidence rates. I. Their nature and limitations. J Epidemiol Community Health. 1980;34(3):223–8. https://doi.org/10.1136/jech.34.3.223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Canadian Physician Data Canada: Canadian Medical Association; 2020 [Available from: https://www.cma.ca/canadian-physician-data.

  27. Lukka H, Warde P, Pickles T, Morton G, Brundage M, Souhami L, et al. Controversies in prostate cancer radiotherapy: consensus development. Can J Urol. 2001;8(4):1314–22.

    CAS  PubMed  Google Scholar 

  28. Ayubi E, Mansournia MA, Motlagh AG, Mosavi-Jarrahi A, Hosseini A, Yazdani K. Exploring neighborhood inequality in female breast cancer incidence in Tehran using Bayesian spatial models and a spatial scan statistic. Epidemiol Health. 2017;39:e2017021. https://doi.org/10.4178/epih.e2017021.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Ann Inst Stat Math. 1991;43(1):1–20. https://doi.org/10.1007/BF00116466.

    Article  Google Scholar 

  30. Anselin L. Global Spatial Autocorrelation. Github; 2018.

    Google Scholar 

  31. Chaix B, Merlo J, Chauvin P. Comparison of a spatial approach with the multilevel approach for investigating place effects on health: the example of healthcare utilisation in France. J Epidemiol Community Health. 2005;59(6):517–26. https://doi.org/10.1136/jech.2004.025478.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Anselin LF, R. J. G. M.; Rey, S.J. Advances in Spatial Econometrics. Anselin LF, M. M.; Hewings, G. J. D.; Nijkamp, P.; Snickars, F., editor. New York: Springer; 2004, DOI: https://doi.org/10.1007/978-3-662-05617-2.

  33. Anselin L. Local Indicators of Spatial Association - LISA. Geographic Analysis. 1995;27(2):93–115 DOI: j.1538–4632.1995.tb00338.x.

  34. Kulldorff M. A spatial scan statistic. Communications in Statistics - Theory and Methods. 2007;26(6):1481–96. https://doi.org/10.1080/03610929708831995.

    Article  Google Scholar 

  35. Han J, Zhu L, Kulldorff M, Hostovich S, Stinchcomb DG, Tatalovich Z, et al. Using Gini coefficient to determining optimal cluster reporting sizes for spatial scan statistics. Int J Health Geogr. 2016;15(1):27. https://doi.org/10.1186/s12942-016-0056-6.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Dohoo I, Martin W, Stryhn H. Methods in epidemiologic research. Charlotte Town: VER Inc.; 2012.

    Google Scholar 

  37. Goovaerts P, Gebreab S. How does Poisson kriging compare to the popular BYM model for mapping disease risks? Int J Health Geogr. 2008;7(1):6. https://doi.org/10.1186/1476-072X-7-6.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Lawson AB, Browne WJ, Vidal Rodeiro CL. Disease Mapping with WinBUGS and MLwiN2003.

  39. Johnson GD. Small area mapping of prostate cancer incidence in New York state (USA) using fully Bayesian hierarchical modelling. Int J Health Geogr. 2004;3(1):29. https://doi.org/10.1186/1476-072X-3-29.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Thomas A, Best N, Lunn D, Arnold R, Spiegelhalter D. GeoBUGS User Manual2014.

  41. Lawson AB. Bayesian disease mapping: hierarchical modeling in spatial epidemiology: Taylor & Francis group; 2009.

    Google Scholar 

  42. Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992;7(4):457–511.

    Google Scholar 

  43. Brooks SP, Gelman A. General methods for monitoring convergence of iterative simulations. J Comput Graph Stat. 1998;7(4):434–55. https://doi.org/10.1080/10618600.1998.10474787.

    Article  Google Scholar 

  44. Heidelberger P, Welch PD. Simulation run length control in the presence of an initial transient. Oper Res. 1983;31(6):1109–44. https://doi.org/10.1287/opre.31.6.1109.

    Article  Google Scholar 

  45. Geweke JF. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. Federal Reserve Bank of Minneapolis; University of Minnesota; 1991.

  46. Raftery AE, Lewis SM. [practical Markov chain Monte Carlo]: comment: one long run with diagnostics: implementation strategies for Markov chain Monte Carlo. Stat Sci. 1992;7(4):493–7.

    Google Scholar 

  47. Raftery AE, Lewis SM. The number of iterations, convergence diagnostics and generic Metropolis algorithms. Practical Markov chain Monte Carlo. London: Chapman and Hill; 1995.

    Google Scholar 

  48. Anselin L, Syabri I, Kho Y. GeoDa: an introduction to spatial data analysis. Geographic Analysis. 2006;38(1):5–22. https://doi.org/10.1111/j.0016-7363.2005.00671.x.

    Article  Google Scholar 

  49. Kulldorff M, Information Management Services I. SaTScanTM v9.4: Software for the spatial and space-time scan statistics. 2015.

  50. QGIS.org. QGIS Geographic Information System. Open Source Geospatial Foundation 2018.

  51. Lunn D, Spiegelhalter D, Thomas A, Best N. The BUGS project: Evolution, critique, and future directions. Stat Med. 2009;28(25):3049–67. https://doi.org/10.1002/sim.3680.

    Article  PubMed  Google Scholar 

  52. Plummer M, Best N, Cowles K, Vines K, Sarkar D, Bates D, et al. Package ‘coda’. CRAN; 2019.

    Google Scholar 

  53. Catelan D, Lagazio C, Biggeri A. A hierarchical Bayesian approach to multiple testing in disease mapping. Biom J. 2010;52(6):784–97. https://doi.org/10.1002/bimj.200900209.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Plascak JJ, Fisher JL, Paskett ED. Primary care physician supply, insurance type, and late-stage cancer diagnosis. Am J Prev Med. 2015;48(2):174–8. https://doi.org/10.1016/j.amepre.2014.08.014.

    Article  PubMed  Google Scholar 

  55. Ross LE, Hall IJ, Howard DL, Rim SH, Richardson LC. Primary care physicians beliefs about prostate-specific antigen evidence uncertainty, screening efficacy, and test use. J Natl Med Assoc. 2018;110(5):491–500. https://doi.org/10.1016/j.jnma.2017.12.007.

    Article  PubMed  Google Scholar 

  56. Nyame YA, Gulati R, Tsodikov A, Gore JL, Etzioni R. Prostate-Specific Antigen Screening and Recent Increases in Advanced Prostate Cancer. JNCI Cancer Spectr. 2021;5(1):pkaa098 DOI: https://doi.org/10.1093/jncics/pkaa098.

  57. Ilic D, Neuberger MM, Djulbegovic M, Dahm P. Screening for prostate cancer. Cochrane Database Syst Rev. 2013;1:CD004720. https://doi.org/10.1002/14651858.CD004720.pub3.

    Article  Google Scholar 

  58. Shoag JE, Nyame YA, Gulati R, Etzioni R, Hu JC. Reconsidering the trade-offs of prostate Cancer screening. N Engl J Med. 2020;382(25):2465–8. https://doi.org/10.1056/NEJMsb2000250.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Jemal A, Fedewa SA, Ma J, Siegel R, Lin CC, Brawley O, et al. Prostate Cancer incidence and PSA testing patterns in relation to USPSTF screening recommendations. JAMA. 2015;314(19):2054–61. https://doi.org/10.1001/jama.2015.14905.

    Article  CAS  PubMed  Google Scholar 

  60. Penchansky R, Thomas JW. The concept of access: definition and relationship to consumer satisfaction. Med Care. 1981;19(2):127–40. https://doi.org/10.1097/00005650-198102000-00001.

    Article  CAS  PubMed  Google Scholar 

  61. Wilson CR, Rourke J, Oandasan IF, Bosco C. Progress made on access to rural health care in Canada. Can Fam Physician. 2020;66(1):31–6.

    PubMed Central  Google Scholar 

Download references

Acknowledgements

This research was enabled in part by support provided by the Information and Communication Technology group of University of Saskatchewan and its Advanced Research Computing resources.

Funding

This work was supported through funding from the Prostate Cancer Fight Foundation and the TELUS Ride for Dad. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

MA contributed to the conceptualization, methodology, formal analysis, data curation, writing original draft of the manuscript, reviewing & editing the manuscript, and developing the visualization components of the study. MS contributed to the conceptualization, methodology, validation, supervision, reviewing & editing the manuscript, project administration, and funding acquisition components of the study. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Michael Szafron.

Ethics declarations

Consent for publications

Not applicable.

Ethics approval and consent to participate

The University of Saskatchewan’s BioMedical Research Ethics Board provided ethics approval (Bio-REB certificate #15–34). Because this project was a secondary analysis of data, the University of Saskatchewan’s BioMedical Research Ethics Board waived the need for consent.

Competing interests

The authors declare they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Andkhoie, M., Szafron, M. Geographic disparities in Saskatchewan prostate cancer incidence and its association with physician density: analysis using Bayesian models. BMC Cancer 21, 948 (2021). https://doi.org/10.1186/s12885-021-08646-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-021-08646-2

Keywords