Role of Area-Level Access to Primary Care on the Geographic Variation of Cardiometabolic Risk Factor Distribution: A Multilevel Analysis of the Adult Residents in the Illawarra—Shoalhaven Region of NSW, Australia

Background: Access to primary care is important for the identification, control and management of cardiometabolic risk factors (CMRFs). This study investigated whether differences in geographic access to primary care explained area-level variation in CMRFs. Methods: Multilevel logistic regression models were used to derive the association between area-level access to primary care and seven discrete CMRFs after adjusting for individual and area-level co-variates. Two-step floating catchment area method was used to calculate the geographic access to primary care for the small areas within the study region. Results: Geographic access to primary care was inversely associated with low high density lipoprotein (OR 0.94, CI 0.91–0.96) and obesity (OR 0.91, CI 0.88–0.93), after adjusting for age, sex and area-level disadvantage. The intra-cluster correlation coefficient (ICCs) of all the fully adjusted models ranged between 0.4–1.8%, indicating low general contextual effects of the areas on CMRF distribution. The area-level variation in CMRFs explained by primary care access was ≤10.5%. Conclusion: The findings of the study support proportionate universal interventions for the prevention and control of CMRFs, rather than any area specific interventions based on their primary care access, as the contextual influence of areas on all the analysed CMRFs were found to be minimal. The findings also call for future research that includes other aspects of primary care access, such as road-network access, financial affordability and individual-level acceptance of the services in order to gain an overall picture of the area-level contributing role of primary care on CMRFs in the study region.


Data
The study used three different databases: (1) the CMRF pathology test data; (2) primary care provider data; and (3) the estimated resident populations from the 2011 Australian Census of Population and Housing. The CMRF test data were extracted from the Southern IML Research (SIMLR) Study database. The SIMLR Study database comprises de-identified and internally linked pathology results from a major pathology provider in the study region and provides near-census coverage of the study population [8]. The CMRF test data were extracted for multiple risk factors on the most recent test results, of nonpregnant adults aged 18 years and over, undergoing a laboratory test between 1 January 2012 and 31 December 2017.
The primary care provider data were manually extracted in 2016 from publicly available data sources, including Yellow Pages, White Pages, online general practitioner (GP) appointment booking services and Google search results. The 2011 Australian Census of Population and Housing data were accessed to extract the population denominator data of the study region at SA1 level [34].

Dependent Variable
Dichotomised results of the individual CMRF tests were the dependent variables in this study. The CMRF test results included: fasting blood sugar level (FBSL); glycated haemoglobin (HbA1c); total cholesterol (TC); high density lipoprotein (HDL); urinary albumin creatinine ratio (ACR); estimated glomerular filtration rate (eGFR); and objectively measured body mass index (BMI). The pathology service routinely collects BMI for each of the remaining CMRF tests and thus became available for

Dependent Variable
Dichotomised results of the individual CMRF tests were the dependent variables in this study. The CMRF test results included: fasting blood sugar level (FBSL); glycated haemoglobin (HbA1c); total cholesterol (TC); high density lipoprotein (HDL); urinary albumin creatinine ratio (ACR); estimated glomerular filtration rate (eGFR); and objectively measured body mass index (BMI). The pathology service routinely collects BMI for each of the remaining CMRF tests and thus became available for analyses in this study [12]. However, it should be noted that the data do not include blood pressure readings. Although blood pressure is an important CMRF it is not routinely collected with any of the pathology test samples and thus not available for analyses in this study. During analyses, all the retrieved CMRF test results were dichotomised into higher risk and lower risk values based on established risk classification guideline values. Table 1 shows the CMRF definitions adopted in this study to dichotomies the test results. High HbA1c HbA1c > 7.5% [38]. 3.

Independent Variable
Primary care access index calculated for the small areas at SA1-level was the independent study variable. An access index score was calculated for each SA1 using a two-step floating catchment area (2SFCA) method, which balanced both supply and demand of primary care services in the study region.
The 2SFCA method was developed by Luo and Wang in 2003 to measure geographic accessibility of health care services [43]. The method has undergone several enhancements since its inception but essentially consists of two steps underpinned by a gravity model [43,44]. The first step computes a population-to-provider ratio for each primary care service location by aggregating the population size of the SA1s whose centroids (i.e., the geometric centers) are located within a defined spatial buffer distance [45]. The total number of general practitioners working in the primary care service locations within this buffer distance were the numerators for the provider-to-population ratio calculations.
Thus, step 1: where Sj is the number of general practitioners at location j, pi is the number of adult residents in the SA1s (those SA1 geographic centroids are located within the spatial buffer distance of the primary care locations) and Rj is the population-to-provider ratio for service j [45].
In step 2, a population-to-provider ratio (access score) is computed for each geographic centroid of the SA1s by aggregating all primary care service population-to provider ratios of the primary care services that are located within the same spatial buffer distance [45].
Thus, step 2: where A i is the access index for population location i. The resulting access indices were retained as a continuous variable for the analyses. A higher score indicated better geographic access of the SA1s to primary care services.
A spatial radial buffer distance of 30 km was chosen to compute primary care access for SA1s in the study region. In the preliminary stage, sensitivity analyses were performed using 1 km, 16 km and 30 km spatial radial buffer distances. In step 1 2SFCA analyses, the 1 km distance covered only 545 (56%) SA1 centroids in the study region in relation to the primary care provider locations, whereas a 16 km radial buffer distance covered 973 (~99%) and a 30 km radial buffer distance covered 978 (~100%) SA1s' geographic centroids. Therefore, a radial buffer distance of 30 km was chosen to determine the access which was observed to cover the mixed rural, semi-rural and urban distribution of the population in the study region well.

Covariates
The individual-level variables adjusted at SA1-level were: sex (male and female) and age group (18-29, 30-39, 40-49, 50-59, 60-69, 70-79 and 80+ years). The area-level covariate adjusted at SA1-level was the area-level socioeconomic disadvantage. The Index of Relative Socioeconomic Disadvantage (IRSD) score of the SA1s in the study region based on the 2011 Australian Bureau of Statistics conducted census of population and housing was used as the measurement variable for the area-level socioeconomic disadvantage of the SA1s [37]. The IRSD summarises a range of measures of relative socioeconomic disadvantage of people and/or households within SA1s and includes: level of income; education; employment; family structure; disability; housing; transportation; and internet connection [37]. A higher IRSD score indicated lower levels of disadvantage [37]. The Illawarra-Shoalhaven region has a diverse socio-economic profile, making this landscape useful for area-level population health studies [46].

Statistical Analyses
Multilevel logistic regression models were fitted to individual CMRF test data at the SA1 level. For each of the seven CMRFs analysed in this study, five nested models were fit that included fixed effects for access index after adjusting for sex, age and IRSD score; and random effect intercepts for SA1s. In the nested models, Model 1 (M1) was a null model of CMRF at SA1-level; Model 2 (M2) included the area-level study variable (access index) only; Model 3 (M3) included individual-level factors at SA1-level (age and sex) only; Model 4 (M4) included individual and area-level factors (age, sex and IRSD score) at SA1-level; and Model 5 (M5) included M4 variables plus access index. Thus, the final model (M5) estimated the effect of primary care access after adjusting for individual and area-level factors. Odds ratios (ORs) were derived from the exponentials of regression coefficients from fitted models. As the IRSD scores and access index of the SA1s were fitted as mean-centred continuous variables, ORs were expressed per standard deviation unit change in these variables. Statistical significance of the models was evaluated using likelihood ratio tests and a type I error rate of 0.05.

Model Comparison
Model fit was compared using the Akaike Information Criterion (AIC). The models were also evaluated for: area-level variance (τ 2 ); proportional change in variance (PCV) in comparison with the null model; intra-cluster correlation coefficient (ICC) of the model; and the median odds ratios (MORs). The ICC and MOR of the models were used to index the between-area variability. A latent variable approach was used to derive the ICC of models [47]. The MOR translates the area-level variance into an easily interpretable OR and is assumed to be statistically independent of the test specific prevalence of the CMRFs [48]. The unique contribution of the primary care access of the SA1s to the area-level variance of CMRF was estimated through the reduction in PCV between M4 and M5.

Statistical Package
All mapping and geospatial measurements were performed using ArcGIS version 10.4.1 (ESRI Inc. Redlands, CA, USA) [49]. All statistical analyses were performed using R version 3.4.4. (R Foundation for Statistical Computing, Vienna, Austria) [50]. Multilevel models were fit using the glmer function in the lme4 package [51]; and likelihood ratio tests were calculated using the lrtest function in the lmtest package [52]. The glmer function fit the generalized linear mixed model, which incorporates both fixed-effects parameters and random effects in a linear predictor, via maximum likelihood [53,54].

Results
A total of 1,132,029 CMRF test results for 256,525 individual residents in the Illawarra-Shoalhaven region between 2012 and 2017 were extracted for analysis. The mean number of tests undertaken per person was 4.4 (SD = 1.8, range = 1-7). After excluding 1162 (1.0%) test results with incomplete details, a total of 1,130,894 tests were retained in the final data set. IRSD scores of the SA1s were the most frequent missing variable, as this was not available for some SA1s in the study region [55]. Available IRSD scores ranged between 446.7 and 1143.7 (mean = 976.7, SD = 98.6) for SA1s, with a higher score indicating lower area-level disadvantage. Table 2 details the individual-level CMRFs risk proportions of the final data set.
For primary care access, a total of 165 primary care service locations with 611 general practitioners were identified in the study area in 2016. The primary care access index of the SA1s in the study region ranged between 0 and 5.41 general practitioners per 1000 people (mean = 2.1, SD = 0.77). Figure 2 illustrates the distribution of the primary care access index within the study region. total of 1,130,894 tests were retained in the final data set. IRSD scores of the SA1s were the most frequent missing variable, as this was not available for some SA1s in the study region [55]. Available IRSD scores ranged between 446.7 and 1143.7 (mean = 976.7, SD = 98.6) for SA1s, with a higher score indicating lower area-level disadvantage. Table 2 details the individual-level CMRFs risk proportions of the final data set. For primary care  Multilevel logistic regression models for each CMRF are presented in Tables 3-9. The null models indicated significant geographic variation in the distribution of all CMRFs at the SA1 level. Model 2s showed inverse associations between access index and all CMRFs except TC, which displayed a positive association with the access index. Model 3s adjusted CMRF models for individual-level age and sex, which accounted for 1.5% (obesity) to 87.3% (eGFR) of unexplained variation in the null model. The general contextual effect of areas over and above their individual composition, such as age and sex, is obtained by a measure of clustering (i.e., ICCs) in the model 3s, which ranged between 0.6-3.4% in the CMRFs models presented. Model 4s demonstrated significant inverse associations between area-level socioeconomic disadvantage and all CMRFs except for TC after adjusting for individual-level factors. Total cholesterol again showed a positive association with area-level socioeconomic disadvantage. In the final models (M5s), the access index was found to be inversely associated with low HDL (HDL < 1 mmol/L) and obesity (BMI ≥ 30 kg/m 2 ), after adjusting for individual and area-level factors. Including the access index in the final models did not attenuate associations between area-level disadvantage and CMRFs observed in M4s.
Model fit of the nested models of each CMRF were compared using the Akaike Information Criterion (AIC). The AIC estimated the out-of-sample prediction error rates of individual models and thus the relative quality of individual models for a given set of nested models [56]. Reductions in the AIC values were observed for all CMRFs models, except in TC and eGFR, from the null model (M1) to the final model (M5), indicating a better fit of the final models. The AIC for TC and eGFR models indicated M4 was the best fitting model for these CMRFs.
In the null models (M1s), low eGFR demonstrated the most area-level variance and high TC showed the least. The access only models (M2s) showed a reduction in the residual variance of all CMRFs from those of null models. In Model 3s, adjusting for age and sex initially increased the residual variance of FBSL (PCV = +1.9%), HbA1c (PCV = +3.0%), HDL (PCV = +15.3%) and BMI (PCV = +1.5%). In Model 4s, adjusting the CMRFs for individual-level age and sex and area-level disadvantage resulted in major reductions of variance from −33.1% (in TC) to −93.3% (in eGFR). In the final models (M5s), including access index in the models after adjusting for the covariates extended the reduction in variance in all CMRFs, except for TC and eGFR. Including the access index had been observed to increase the variance in the TC and eGFR final models, compared with the lower level model.
Similarly, in the unadjusted models, the MORs, which indicate the odds of having a higher risk CMRF test result for a person from the most, compared to the least, area-level disadvantage, were the highest among eGFR (τ 2 = 0.189; ICC = 5.4%; MOR = 1.51) and the least among TC (τ 2 = 0.025; ICC = 0.8%; MOR = 1.16). The ICCs of CMRFs in all the models were comparatively small (Table 4) in all the models, indicating minimal contextual effect of areas on any of the CMRFs. In the fully adjusted models, the ICCs further reduced and ranged between 0.4% and 1.8% in low eGFR and BMI respectively. Table 10 presents a summary and comparison of the model fit.

Discussion
This study aimed to inform area-specific interventions for the prevention and control of CMRFs, based on the primary care access status of the small areas within Illawarra-Sholhaven region of NSW, Australia. After adjusting for the covariates, we found that: a) greater access to primary care was associated with a reduction in the odds of low HDL and obesity but was not associated with high FBSL, high HbA1C, high TC, high ACR and low eGFR; b) the general contextual effect of areas on each of the CMRFs were minimal; and c) the geographic variation of CMRFs specifically explained by primary care access was small and did not demonstrate any attenuating effect on the contribution of area-level disadvantage on the variation of CMRFs in the study region. The results demonstrate that though the probability of low HDL and obesity decreases with increasing primary care access, the low general contextual effects of the areas on each of the CMRFs (i.e., low ICCs of Model 3s, ranges 0.6-3.4%) indicate minimal difference between the small areas after controlling for the study variables. Thus, the findings suggest that preventive interventions should not only be focused on areas with lower primary care access. Rather, interventions should be universal but proportional to the need and risk level of the people for the prevention and control of CMRFs. Primary care access was associated with all CRMFs in unadjusted models but only with low HDL and obesity in models fully adjusted for individual-and area-level covariates. These findings support the arguments of the possible role of confounders and reverse causality in ecological models [57], which question the previously established associations between primary care access and improved health [58,59]. The study suggests higher odds of being identified with low HDL and obesity with reduced access to primary care. In previous studies, when the relationship between health care service outcomes and travel time was modelled using multilevel logistic regression, it was found that GP consultations were less likely to happen when the travel time was longer, which is more common in rural areas [21]. The current study outcomes are consistent with those findings. However, it should also be noted that the current findings pertain only to the geographical/spatial accessibility of the primary health care services within 30 km distance of an SA1 centroid, rather than their road network access, actual usage and affordability.
The primary care access index, derived from the study region, ranged from 0 to 5.41 general practitioners per 1000 people (mean = 2.1, SD = 0.77). Multiple previous studies have reported inequalities in the geographic access to primary care services, using different enhanced versions of the 2SFCA method [45,[60][61][62][63][64][65][66][67][68]. For example, the spatial accessibility index derived from rural Otago in New Zealand, using the travel time distance, ranged between 1 to 10, where a higher score indicated better access [62]. The accessibility index reported from Thimphu district in Bhutan ranged between 0 and 1, where 1 was the maximum access [69]. The spatial accessibility index of GP accessibility in England has been reported to range between 7.2 (South of England) and 13.3 (in London) [69]. The access map of the study region (Figure 2) clearly shows a polarisation of the higher access indices along the northern and southern ends of the study region, thus a visible inequality in their distribution. The WHO recommends universal access to primary care for all populations, where geographic access is one part of physical access to primary care [70].
Area-level disadvantage explained more geographic variation in CMRFs than area-level access to primary care. Inclusion of the access index in the final model did not demonstrate any reduction in the variance explained by area-level disadvantage on the geographic variation of CMRFs. This finding supports the importance of overall socioeconomic development of areas to reduce CMRF risk. Moreover, the ICC values of the final models were too small to suggest any meaningful area-level difference in the modelled CMRF variables. This would support the call for universal approaches for the prevention and control of CMRFs rather than any targeted area-level approaches, but with a proportional priority to disadvantaged populations in the study region [24,28,71].
This study has to be considered within its limitations. First, the cross-sectional nature of the study does not support causal inference. Second, the CMRF data used in this study are from people already utilising health care service in the study area, so care should be taken in generalising the results to the overall population. The SIMLR database does not include hospital or emergency service based tests. Therefore we believe that the database has a reasonable representation of community dwelling adults in the study region. However, it should be noted that the study sample includes only people who have accessed health care and pathology services, and the omission of those who have not accessed care may have biased our results. Given our population coverage this seem unlikely. Third, the study used a radial buffer distance of 30km for access calculations rather than travel time/distance because proprietary road network data were unavailable for this study. Thus, the patients' actual experiences of seeking physical access in daily life need not exactly reflect the compound measure of access index adopted in this study. Even though the 30km buffer distance helped to include a maximum coverage of the population in relation to the geographic location of the primary care providers, this distance might have also influenced the discriminatory accuracy of the SA1s in the multilevel analyses. In addition, it should also be noted that the access index described in the study pertains only to the geographical reachability, but not to the affordability and acceptability of the available services. Forth, the study did not include blood pressure as a variable, although it is a major CMRF, due to non-availability of data. We were also unable to adjust for ethnicity for the same reason.
The main strength of this study is the use of a large population-derived database comprising a wide range of CMRFs. The research adds to the very few studies which consider multiple CMRF variables from the same region [18][19][20][72][73][74] and is unusual for its hierarchical analysis of the associations between a range of CMRFs and primary care access in a widely dispersed population.
Future research is required to investigate other area-level attributes contributing to the geographic variation of CMRFs in the study region. Our previous research has reported that area-level disadvantage contributes 14.7-57.8% of the geographic variation in CMRFs. The current study extended the previous findings by identifying the specific contribution of area-level primary care access, ranging between 0.0-10.5%. Further area-level analyses are required to identify other factors contributing to the geographic inequality of the CMRFs in the study region.

Conclusions
The findings of the study suggest that adults residing in areas that have a poor primary care access are more likely to be identified with low HDL and obesity. However, the specific contribution of area-level primary care access was small when compared to the contribution of area-level disadvantage. The finding supports the importance of overall socioeconomic development of areas to reduce CMRF risk, while supporting universal approaches for the prevention and control of CMRFs which are proportional to the need and disadvantage level of the individuals. Future research including other aspects of primary care access such as road-network access, financial affordability and acceptance of the services might help to provide an overall picture of the contributing role of primary care in the study region.