This work was funded by the Belgian Science Policy Office (BELSPO), grant number: BR/143/A3/GRESP-HEALTH. The authors acknowledge the Intermutualistisch Agentschap – Agence Inter-Mutualiste (IMA-AIM) and Mariska Bauwelinck and Professor Patrick Deboosere (Vrije Universiteit Brussel) for providing the data and for their contribution to the early version of this paper.

1. Introduction

1Living environment and health are closely associated. Among others, salutogenic effects of the proximity to green spaces have been largely demonstrated on physical and mental health (see e.g. Ulrich, 1984 ; Sacker and Cable, 2005 ; Nielsen and Hansen, 2007 ; Maas et al., 2009 ; Van Herzele and de Vries, 2012 ; Gascon et al., 2015 ; Pietilä et al., 2015 ; Mukherjee et al., 2017), while high levels of air pollutants such as PM₁₀ or NO₂ are associated with respiratory and cardiovascular problems and consequent medication intakes (see e.g. Casas et al., 2016 ; Sofianopoulou et al., 2013). Recent studies insist on additional adverse effects of air pollution on the brain, birth outcomes (e.g. low birth weight) or metabolic disorders like diabetes mellitus (Landrigan et al., 2017 ; Thurston et al., 2017).

2The association between environment and health is most often investigated at the local level and within selected urban areas (see e.g. Cheng et al., 2011), in particular when tackling issues of increasing urban population, for targeted policies and investments on green infrastructures (Conway et al., 2010 ; Coutts and Hahn, 2015 ; Adli et al., 2017). City dwellers are generally considered as living in less healthy environments, while enjoying the proximity of a large variety of facilities. Non-urban areas, experiencing lower built-up densities are generally considered as healthier environments, although they also suffer from some types of pollution due, for instance, to industrial or agricultural practices (Paulot and Jacob, 2014), including intensive animal breeding or the use of pesticides. Studies analyzing spatial variations at the level of an entire country are, however, very scarce (Henau et al., 2015).

3Health status is generally proxied by different types of data such as individual perceived health obtained from surveys (Maas et al., 2006 ; Triguero-Mas et al., 2015), hospital records (Dadvand et al., 2012 ; Pereira et al., 2012), and more recently by medication prescriptions (Halasa et al., 2002 ; Naughton, Bennett and Feely, 2006 ; Bianchi et al., 2009), or data derived from time-series and cohort studies, in the case of mortality (Samet and Krewski, 2007). Although medication prescriptions have been used in the literature for specific medications, mapping their spatial variations is quite novel and rare (Cheng et al., 2011). Here, werefer to Wangia and Shireman (2013) for a meta-analysis on the application of GIS tools for analyzing spatial variations of medication prescriptions. The vast majority of studies have investigated spatial variations of health by looking at individual morbidity or mortality causes, rather than looking at their covariation in space (see e.g. Waldhoer et al., 2008).

4In a sister paper, we described the spatial variation of asthma medication prescriptions for preschoolers (0-5 years old) in Belgium (Trabelsi et al., 2019); the paper uncovered clear spatial structures for the prevalence of asthma, but it also revealed that the spatial structure of the medication cost per case differed substantially from that of the prevalence. The present paper deals with a larger number and variety of health disorders (via ATC medication codes) and considers all age groups. We a priori presumed that the history of the Belgian economic and social development (former coal mines, etc.), the urbanization level, the social realities, have generated concentrations of specific diseases, and hence over/under consumptions of specific types of medications. Here, we conduct exploratory spatial data analyses (E.S.D.A.) to describe the spatial distribution and clustering of age standardized prevalence of medication purchases in Belgium for five groups of medications and their correlation with environmental and socio-economic variables.

2. Data and method

2.1 Medication purchases

5Health insurance is compulsory in Belgium and most residents (98%) are enrolled in the Health Security System. In other words, a very small part (2%) of the population living in Belgium is not covered by the social security system and therefore, not included in this study. These correspond to undocumented immigrants, people with complex social problems (such as homelessness), or employees of international bodies (EU, NATO, etc.). Most health care expenditure including medication prescribed by physicians is reimbursed (at least in part) by one of the seven Belgian health insurance funds. The Intermutualistic Agency (IMA-AIM) is an umbrella organization that centralizes all information on health care reimbursements collected by all Belgian health insurance funds. Information on reimbursed medication is coded using the ATC (Anatomical Therapeutic Chemical) code (see WHO Collaborating Centre for Drug Statistics Methodology, 2018). After medication prescription by a medical doctor and the purchase of that medication in a pharmacy located in Belgium (including hospital pharmacies), the ATC codes are registered in the health insurance system. The database does not include non-reimbursable medications and medications that do not require a prescription.

6In this analysis, medication reimbursement information is provided at the municipality level for all 589 municipalities of Belgium and concerns all individuals having purchased one of the studied prescribed drugs at least once in a calendar year (here named “a case”). Information on its actual use, frequency of use, or about co-morbidities, meaning here the presence of the same individual in more than one medication group, cannot be deducted from the data. In addition, this database does not provide information on the frequency of purchase of the medication for one year. We know that the intensity of use of a medication could lead to other spatial structures (see Trabelsi et al. 2019). 

7Prevalence is here defined as the number of cases divided by the total population registered in each municipality. Age-standardized prevalence is expressed per 100 inhabitants (Age-Standardized Medication reimbursement Prevalence, noted here ASMP). It is provided separately for men and women (see Costa, Bauwelinck and Deboosere, (2019) for details about computation) as it is presumed that health disorders are often gender specific (Glaeske et al., 2012 ; Che et al., 2014 ; Skoog et al., 2014). In order to minimize the effects of small numbers generating statistical and privacy issues, an average prevalence is computed over 3 years (2012-2013-2014).

8We here limit our investigations to five medication groups corresponding to common health disorders, managed at the level of primary care, and that are a priori related to environmental conditions, according to the epidemiological literature (see e.g. Thurston et al., 2017 ; Fong, Hart and James, 2018). Table 1 presents the target system, the ATC codes, the name of the codes and the abbreviations used in this paper for each medication group. RS includes medications for respiratory obstructive diseases such as asthma or Chronic Obstructive Pulmonary Disease (COPD). CVS includes medications used to reduce blood pressure or treat other cardiovascular problems such as ischemic diseases. LIP includes medications prescribed to reduce hyperlipidemia, while IS includes medication for treating allergy symptoms (local and systemic). Finally, CNS includes drugs acting on the central nervous system mainly indicated to treat symptoms of anxiety and depression. Nevertheless, it is important to highlight that some of the listed medications may be used for different disorders [e.g. diuretics (C03) can be used to reduce blood pressure as well as for some renal disorders] and hence their specificity is not absolute.

Table 1 : Five groups of medications used in this study.

	Target system	ATC code	Name
RS	Respiratory system	R03	Drugs for obstructive airway diseases
CVS	Cardiovascular system	B01A	Antithrombotic agents
		C01	Cardiac therapy
		C02	Antihypertensives
		C03	Diuretics
		C07	Betablocking agents
		C08	Calcium channel blockers
		C09	Agents acting on renin-angiotensin system
LIP	Lipids	C10	Lipid modifying agents
IS	Immune system (allergy)	D04A	Antipruritics, incl. antihistamines, anesthetics, etc.
		R06	Antihistamines for systemic use.
CNS	Central Nervous system	N05	Psycholeptics
		N06	Psychoanaleptics

2.2 Study area

9Belgium is a 11 Mio inhabitants, densely populated, European country (374 inhab./km²). It is divided into three administrative regions: Flanders in the North (tight urban network, densely inhabited), Wallonia in the South (greener, less densely populated, looser city network) and the Brussels Capital Region (BCR) in the center. The BCR corresponds to the urban core of the city of Brussels, while the urban agglomeration sprawls into both Flanders (in the Province of Vlaams-Brabant) and Wallonia (in the Province of Brabant Wallon). Belgium is divided into provinces (10 provinces + BCR) as illustrated in Figure S1 in the Supplementary Material, and further into 589 municipalities at the date of the analysis.  The municipality of Herstappe was removed from the analyses because of its small size (88 inhabitants in 2018) (StatBel, 2018). 

10The three Regions have their own political administration, language (Dutch in Flanders, French in Wallonia, German in the German-speaking community, and both Dutch and French in the BCR), and universities delivering medical degrees. Reimbursement criteria are defined at the national level, but regional factors may induce differences in practices.

2.3 Environmental and socio-economic variables

11Our aim is not to statistically explain the spatial variation of medication consumption by a set of variables, but to look at their co-occurrence in a purely exploratory (descriptive) perspective. Hence, we gathered information on several standard environmental (ENV) and socio-economic (SE) characteristics at the level of the municipalities. 

12Data on air pollution were provided by the Belgian Interregional Environment Agency (Irceline, 2018) : data used in this study include average values of Particulate Matter 2.5 (PM_2.5), Particulate Matter 10 (PM10), Nitrogen Dioxide (NO2), and Black Carbon (BC) in 2015. These are measured in fixed monitoring stations and further modelled by 4 ⋅ 4 km cells. To link the air pollution data with the medication data, the pollutant concentrations were estimated at the municipality level. 

13Variables on green spaces were obtained from the Corine Land Cover dataset (European land-use dataset) developed by the European Environment Agency (Larsen, 2002). Green space is here defined as the sum of urban green, forests, pastures and agricultural land, and semi-natural areas. For this study we used the percentage of green spaces (hereafter noted “Total GS”) in each municipality, and the average distance from the centroid of the built-up surface to the nearest green space (hereafter noted “Ave Near GS”). 

14Standard socio-economic variables based on existing international literature (Steele et al., 2007 ; Reiss, 2013 ; Mayer and Österle, 2015 ; Vanthomme and Gadeyne, 2019) were used as proxies for education, poverty, type of employment. For the purpose of our study, we used the share of unemployed persons in the 18–64-year-old population (noted Unemployed), the share of individuals aged 25-64 with primary education only (noted Primary ed. 25-64y), the share of individuals aged 25-64 with a higher education degree (noted Higher ed. 25-64y), the share of manual workers (noted Manual workers) and the share of foreign-born inhabitants coming from low income countries (noted Foreign-born). The Belgian Statistical Office (StatBel) provided data on median income and population density in 2014, while socio-economic variables were derived from the latest Belgian Census (2001) as no other census was conducted since. The descriptive statistics about these variables are reported in Table S1 in the Supplementary Material.

2.4 Methodology

15The five age-standardized medication reimbursement prevalence (ASMP) were first mapped at the municipality level. For all these maps, Jenks discretization was applied in order to minimize the intra-class and maximize the inter-class differences. The visual output of each map is accompanied by the computation of a Moran’s I index for quantifying spatial autocorrelation. Spatial autocorrelation has been assessed considering all municipalities sharing at least one border (direct neighbors).

16Correlations between all five ASMP variables were computed using Pearson coefficients (Table 5). A hierarchical clustering of the five ASMP variables was then computed using Ward’s method to identify clusters of municipalities with similar medication reimbursement spatial patterns (Figure 3).

17Pearson correlations were computed between all ASMP and ENV / SE variables (not reported here) and further summarized by a Principal Component Analysis (PCA) (see Table S2) as well as on the ENV and SE variables only (see Table S3). The PCA is here used as an exploratory tool aiming at simply reducing the size of the dataset in terms of variables : components are a linear combination of initial variables that co-vary spatially. PCA was therefore used in an exploratory way, to present which variables co-vary best and if they do. A Varimax rotation was applied. We limited ourselves to factors with eigenvalues >1. A Ward clustering is then performed following the PCA to present a spatial partition of Belgium in terms of the variables considered.

18A summary of the data used, and the analysis carried out in the paper can be found in Figure 1. 

3. Analysis and results

3.1 Mapping the Age-Standardised Medication reimbursement Prevalence 

19Basic descriptive statistics about the ASMP variables are reported in Table 2. Results for both men and women are presented separately. Since the gender difference is not significant (almost identical descriptive statistics and comparable medication reimbursement patterns - see correlation coefficients for each medication group between men and women noted r in Table 2), we further limit our report to men only (women data and results are not illustrated, but available upon request : statistical differences are very limited).

Table 2 : Descriptive statistics of the prevalence of the 5 groups of reimbursed medications (ASMP) in Belgium, by gender 

Medication Group	Gender	Mean	Std dev	Min	Max	Median	r
RS	men women	14.23 14.01	2.33 2.52	5.44 5.34	20.4 21.99	14.09 13.84	0.95
CVS	men women	27.64 27.48	2.87 2.64	14.47 15.05	34.09 33.01	27.6 27.57	0.95
LIP	men women	15.26 15.21	1.85 2.09	6.51 7.23	20.57 21.16	15.43 15.26	0.86
IS	men women	10.78 10.67	2.33 2.12	3.09 2.84	19.44 19.61	10.73 10.51	0.97
CNS	men women	16.7 16.83	2.88 2.75	7.33 8.01	25.54 23.59	16.72 16.66	0.96

n= 588   Std dev : standard deviation, Min : minimum value, Max : maximum value,  r : Pearson correlation coefficient between ASMP for men and ASMP for women for each medication group.

20The 5 choropleth maps of the prevalence of medication reimbursement are reported in Figure 2. Contrary to our expectations, spatial structures look alike, except for IS which concerns mainly medications associated with allergies. Medication reimbursements for the immune system IS (Figure 2) are more prevalent in the North-East of the country and more particularly in the Province of Antwerp, and partially in the South-East near the border with France (South in the Province of Hainaut, and in the western part of the Province of Namur). The Province of Hainaut and Namur are both characterized by a high prevalence of medication reimbursements for the respiratory system (RS - Figure 2), the cardiovascular system (CVS – Figure 2), and the central nervous system (CNS - Figure 2), all expanding towards the East. Wallonia, with the exception of Walloon Brabant, is generally characterized by high prevalence of medication reimbursements, regardless of the type of ASMP. Brussels (BCR) tends to have lower values of most ASMP variables, whereas the eastern part of the province of Antwerp seems to concentrate medication purchase for allergic symptoms.

21A Moran’s I index is associated with each map. All five ASMP measures are characterized by a positive and significant autocorrelation, all varying between +0.50 and +0.67 : municipalities that look alike tend to be closely located.

3.2 Comparing ASMPs 

22The Pearson correlation coefficients between the five medication groups (Table 3) confirm the visual analyses: coefficients are positive and statistically significant. The strongest correlations are observed between the cardiovascular system medication group (CVS) and all other medication groups with the exception of the immune system (IS) targeted medications. From Figure 2 it is clear that IS has a different spatial structure, hardly explainable by common knowledge. The high concentration in the East of Antwerp may be due to the presence of a natural allergen or to polluted air emitted in and around Antwerp and pushed by western dominant winds. These hypotheses should be investigated and verified by further individual medical analyses and diagnoses.

23High correlations between ASMPs can also be explained by multimorbidity and polypharmacy as recently discussed by Van Den Akker et al., (2019) for Flanders.

Table 3 : Pearson correlation coefficients between the five groups of medications in Belgium

	RS	CVS	LIP	IS	CNS
CVS	0.82	1
LIP	0.54	0.71	1
IS	0.41	0.44	0.40	1
CNS	0.68	0.81	0.57	0.39	1

All r are significant at 0.0001. n = 588 municipalities.

3.3 Clustering Belgian municipalities with the 5 ASMP variables

24A Ward clustering method is applied on the 5 ASMP variables to obtain groups of municipalities that are similar in terms of medication reimbursement practices. Three clusters are optimally identified. Contrary to expectations, fragmentation is small: each cluster consists of contiguous municipalities, meaning that close municipalities have more similar ASMP profiles than distant ones (Figure 3). Cluster 1 (in yellow) extends principally in the North of the country (Flanders), but also in the South-East in the province of Liège (including, among others, the German speaking Community of Belgium) and partly in the province of Luxembourg. Despite its location within Wallonia, most of the municipalities of the province of Brabant Wallon are also grouped into Cluster 1. Cluster 2 (green) mainly includes (1) municipalities in and around the Eastern and South-eastern part of Brussels, and (2) municipalities close to the border of Belgium and/or in the vicinity of an important international workplace (i.e. Luxembourg). Cluster 2 is characterized by the lowest values of the means (Table 4), suggesting low prescription prevalence in those municipalities and/or the purchase of medications in the neighboring country where they may be cheaper or more easily available. Cluster 3 (in red) extends mainly in Wallonia, although it also includes a few municipalities scattered in Flanders. Cluster 3 is characterized by high means for each variable, pointing to higher prescription rates of medications, whatever the type. 

Table 4 : Cluster composition for the 5 ASMP variables (mean and standard deviation)

	Cluster 1 (n=361)	Cluster 2 (n =56)	Cluster 3 (n =171)
RS	13.68 (1.44)	10.27 (1.60)	16.68 (1.30)
CVS	27.05 (1.55)	22.15 (2.16)	30.67 (1.34)
LIP	15.24 (1.50)	12.21 (1.57)	16.28 (1.45)
IS	10.89 (2.35)	7.46 (1.66)	11.64 (1.32)
CNS	15.87 (1.97)	12.53 (1.84)	19.82 (1.48)

n refers to the number of municipalities included in cluster

3.4 Co-variation of ASMP and environmental / socio-economic variables

25The thirteen environmental and socio-economic variables defined in Section 2.3 are now also taken into consideration. These variables are potentially related to the five ASMP variables as they are a priori assumed to play a (major) role in health issues (and hence medication consumption). Correlations are significant but not extremely high (see Table 5).  

26Most correlations between ENV variables and ASMP are negative with the exception of (1) the percentage of total green (TotalGS), meaning that the ASMP increases with greenness, which may appear as counter-intuitive but consistent with the North South divisions illustrated in the preceding maps, and is also consistent with recent research about childhood asthma (Aerts et al., 2020), and (2) between allergies and PM_2.5. In other words, medication reimbursements are poorly and negatively correlated to environmental variables: the higher the pollution, the smaller the ASMP values. This does not mean that pollution and health are not associated but that the spatial variations of the measures are little associated, at least if no control variable is considered. 

27Looking at SE variables, we observe that most ASMP variables are positively correlated to the proportion of unemployed, while central nervous system (CNS) and respiratory (RS) medication reimbursements are positively correlated to the share of population (25-64y) with primary education only. Strong negative correlations are to be found between the proportion of people with higher education (25-64y), and foreign born and most groups of medications (except for the immune system (IS), and central nervous system (CNS), respectively). Also, median income and population density both show a negative sign with asthma and chronic obstructive pulmonary disease (RS), and cardiovascular problems (CVS), and median income is also negatively correlated with depression and anxiety (CNS). 

Table 5 : Pearson correlation coefficients between the five ASMP and the ENV (top rows) and SE variables (bottom rows)

	RS	CVS	LIP	IS	CNS
PM2.5	-0.29***	-0.29***	n.s.	0.15***	-0.30***
PM10	-0.20***	-0.19***	n.s.	0.12**	-0.22***
NO2	-0.42***	-0.43***	-0.22***	n.s.	-0.29***
BC	-0.33***	-0.32***	-0.18***	n.s.	-0.23***
Total GS	0.32***	0.30***	0.16***	n.s.	0.14***
Ave near GS	-0.33***	-0.30***	-0.26***	-0.14***	-0.15***
Manual workers	n.s.	n.s.	0.18***	0.17***	n.s.
Unemployed	0.37***	0.41***	0.16***	n.s.	0.52***
Foreign-born	-0.28***	-0.29***	-0.28***	-0.19***	-0.08*
Primary ed. 25-64y	0.14***	n.s.	n.s.	-0.14***	0.11***
Higher ed. 25-64y	-0.34***	-0.35***	-0.29***	-0.16**	-0.32***
Median income	-0.26***	-0.50***	-0.10**	n.s.	-0.50***
Population density	-0.30***	-0.31***	-0.26**	-0.13**	-0.13**

*** : 0.001, ** : 0.01, n.s. : non significant 

28A PCA is first performed on all environmental (ENV), socio-economic (SE) and ASMP variables in order to see how they are related to each other at the municipal level within Belgium (see Table S2 in Supplementary Material). Factors selected were limited to those with an eigenvalue greater than 1. Environmental, socio-economic and medication reimbursements are independent from each other. Environmental variables (PM_2.5, PM₁₀, NO₂, BC) have high positive loadings on Factor 1 as well as the percentage of green space, which has a negative loading. Factor 2 gathers socio-economic variables with high positive loadings on the share of unemployed, the share of foreigners born in low-income countries, population density and distance to the closest green space, while median income has a negative sign. Factor 3 is composed of the five ASMP variables, confirming their specific and independent spatial structure. Factor 4 gathers SE variables related to education (high proportion of low education, manual jobs; low proportion of high education). Several PCAs were also performed separately for the two Belgian Regions (Flanders and Wallonia). Given the geography of the country, Brussels Capital Region was included once in the Flemish Region, and once in the Walloon Region. Results (not illustrated here), confirm the independence of the three groups of variables (ASMP, ENV and SE) discussed above.

29As previously mentioned, and as represented in Figure 1, an additional PCA was conducted on the thirteen ENV and SE variables (no ASMP variable included) in order to synthetize the environmental and socio-economic structure of the country into a few new variables and to further explore their link with the reimbursement prevalence. Three factors appear (Table S3 in Supplementary Material). Factor 1 is characterized by positive loadings on air pollution (PM_2.5, PM₁₀, NO₂, and BC), population density, distance to green spaces, and the proportion of individuals born in low-income countries and a negative loading on percentage of green areas. Factor 2 is characterized by lower socio-economic characteristics: positive loadings on share of unemployed, presence of foreigners born in low-income countries, population density, and negative loadings with median income. Factor 3 is composed of a high share of manual workers and share of population (25-64y) with a low education, and consequently a low share of population (25-64y) with higher education.  

30A second cluster analysis was then performed on the scores of these three factors as well as the five ASMP medication groups (see Figure 4 and Table 6). The BCR with its better-off South/Southeastern periphery clearly pops up in Cluster 3 (green), alongside with some border municipalities, close to Luxembourg and Germany. This cluster is characterized by low prevalence of reimbursement for all groups of medications, as well as a higher education level (Factor 3). Cluster 4 (in red) mainly groups Walloon municipalities, with the exception of the Province of Brabant Wallon that clusters with the Northern part of the country. This red cluster is predominantly characterized by high reimbursement rates for medication associated with the respiratory system, the cardiovascular system, and the central nervous system. The Province of Brabant Wallon is part of Cluster 1 (light blue). This cluster covers mainly the North of the country (Flemish Region), with some exceptions in the provinces of Liège and Luxembourg (around the green cluster). It is characterized by low prevalence of reimbursements for medications for the respiratory system and the central nervous system, a higher income, and low unemployment. The last cluster (Cluster 2, in orange) spreads primarily into two areas in the Flemish Region (West Vlaanderen and South-West Limburg/East of Antwerp). Cluster 2 is characterized by high medication prescriptions for the immune system and lipid-lowering medication, a lower share of higher educated people and more manual workers. Of note, South-West of Limburg was deeply hit by the Covid-19 (see Sciensano, 2019).    

31This last clustering confirms the strong and clear spatial structure observed in Figure 3. With some exceptions, clusters seem to follow, again, the linguistic/regional border between Flanders and Wallonia. Adding environmental and socio-economic variables does not change the trend.

Table 6 : Cluster composition : mean and standard deviation of the five ASMP variables and three PCA factors

Variable	Cluster 1 (n=180)	Cluster 2 (n =144)	Cluster 3 (n (n=55)	Cluster 4 (n =209)
RS	13.18  (2.24)	13.89 (2.14)	10.17  (2.24)	16.43 (2.27)
CVS	26.18  (2.78)	27.58 (2.65)	22.07  (2.78)	30.39 (2.72)
LIP	14.45  (1.87)	15.96 (1.82)	12.18  (1.87)	16.27 (1.77)
IS	9.63  (2.38)	12.54 (2.31)	7.56  (2.38)	11.41 (2.08)
CNS	14.79  (2.86)	16.61 (2.77)	12.52  (2.86)	19.52 (2.92)
Factor 1 (high loadings on PM2.5, PM10   NO2 and BC, Ave near green ; low on Total GS)	0.24  (0.98)	0.31  (0.88)	0.68  (0.99)	-0.60 (0.76)
Factor 2 (high loadings on Unemployment, Foreigners (LIC), Population density ; low on Median income)	-0.33  (0.99)	-0.62  (1.00)	0.68  (1.00)	0.54  (0.82)
Factor 3 (high Manual workers, Primary ed. 25-64y ; low Higher ed. 25-64y)	-0.09  (1.02)	0.40  (1.05)	-0.86  (1.00)	0.03  (0.87)

4. Discussion and conclusion

32We initially intended to use medication reimbursement data as a proxy for health. The Belgian sickness fund provided us data about medication prescriptions geocoded at the place of residence at municipality level. The exploratory spatial analyses conducted here show a very consistent spatial pattern, with a strong positive spatial covariation of the five types of medications, as well as their independence from environmental and socio-economic characteristics. But this pattern questions the real meaning of the data provided.

33The data used are based on medications reimbursed by the Belgian sickness funds. This information has intrinsic limitations: first, it does not include people not covered by the Belgian social security system. Nevertheless, the proportion of residents that are not covered is estimated to be very low (less than 2%), therefore the potential selection bias related to this limitation can be considered minimal. Second, non-reimbursable or not prescribed medications are not included, and the purchase and reimbursement of a medication does not equate its consumption. However, the information on medication purchase used here can be considered as a good indicator of the use of the Belgian health care system, as all purchased and reimbursed medication in Belgium are prescribed by a doctor. Medication is prescribed because of a diagnosis; hence, we can confidently assume that the practitioner of the patient buying the medication does consider that the person needs it. Last, the data does not provide information about the consumer (precise geo-location, and characteristics of the person), nor about the quantity of medication prescribed and bought yearly, nor about the exact pathology the drug was prescribed for, or its severity or chronicity. These biases are unfortunately not possible to estimate, nor to quantify. This database is originally built by and for health insurance companies; it is here diverted for potential health geography purposes in a specific multidisciplinary research project (GrespHealth), which aims at analyzing the relationship between health and green/blue environments. 

34We here used age-standardized data at the municipality level. This level of spatial aggregation has the advantage of having a very low number of missing data, it avoids privacy issues and allows comparison with other contextual variables, but it does not allow any in-depth analysis at the individual level, which would have been desirable for causal analysis (see e.g. Kwan, 2009 ; Park and Kwan, 2017). Aggregating data entails indeed possible modifiable areal unit problems (see e.g. Openshaw and Taylor, 1979 ; Dusek, 2004) and/or ecological fallacy errors, that is, inferring about individuals on the basis of spatially aggregated data (see e.g. Openshaw, 1984 ; Wakefield and Lyons, 2010). In statistical analysis, spatial is special. Here, analyzing spatial variations at the level of all municipalities of one country allows identification of regional structures and to question the role of the linguistic/political border, and is relevant in terms of public health in Belgium. However, it is impossible to disentangle the different components such as medical doctors’ education and practice, medication purchase, pharmaceutical sales promotion, role of commercials, for which individual data would be essential. To our knowledge, few studies have analyzed these variations in space as such. Waldhoer et al. (2008) looked at the spatial distribution of infant mortality in Austria for different causes of death. They found a clear non-random spatial distribution of mortality causes, and a non-relevant effect of socio-economic status of the mother. Moreover, this paper also reminds how representing a country (Belgium) by some selected observations in one part of the country (Tyrstrup et al., 2017) may lead to very questionable and false results. Tyrstrup et al. (2017) in fact take one region of Belgium (Flanders) as representative of the whole country on an international scale comparison with the Netherlands and Sweden.

35Last but not least, environmental and socio-economic variables are not straightforward to measure and far from being perfect here. Averaging pollutants such as BC and NO₂ at the municipality level introduces estimation biases (they are measured punctually, interpolated on a grid and aggregated by municipalities that are irregular in shape and size). 

36The sign and value of the correlation coefficients between the five medication groups fit with medical knowledge. The cardiovascular medication group included in our study is not only prescribed to treat diseases but, like the lipid modifying agents, it is also prescribed to prevent the onset of chronic cardiovascular diseases or acute events. The respiratory medication is mainly used to treat symptoms of asthma (childhood and early adulthood) and chronic obstructive pulmonary disease (COPD) (adulthood and elderly). Since COPD shares many risk factors with cardiovascular diseases, it was expected that COPD medication sales (and reimbursement) would be strongly correlated with those of cardiovascular medication. Finally, the strong correlation observed between medication reimbursement for the central nervous system (mainly used among patients with mood disorders) and cardiovascular system medication sales, is consistent with the increasing evidence that mental health problems may be risk factors for cardiovascular disease (Cohen, Edmondson and Kronish, 2015).  

37Clustering the medication reimbursements for the five medication groups shows a clear spatial pattern. The regional (and linguistic) division appears clearly, as well as the effect of the international borders (with Luxembourg in particular), together with the possible effect of the presence of expats with special status in the Region of Brussels. In addition, medication reimbursement prevalence is generally higher in the southern part of the country, a region characterized by lower pollution levels, and a high share of green spaces. Low prevalence of medication reimbursement can also be found in the big cities, hence presenting an urban-rural difference in purchase and reimbursement. All medication reimbursement types are highly autocorrelated (Moran’s I ≥ 0.5). Consistently, this structure cannot be explained by environmental factors, which are not affected by the administrative partition, but it is suggestive of an effect of prescription practices. We included in this study a selection of environmental and socio-economic variables in addition to the medication reimbursement data to better understand the spatial variations in medication prescription. We computed Pearson correlations among ASMP variables and the environmental and socio-economic ones. This yielded unexpected inverse correlations with air pollution (PM_2.5, PM₁₀, NO₂, and BC) and distance to green spaces, while the correlation was positive with the share of green spaces. These results are contrary to the existing literature based on individual studies indicating the adverse effect of air pollution on health, and the beneficial effects of a natural environment (Pereira et al., 2012 ; Sofianopoulou et al., 2014 ; Triguero-Mas et al., 2015 ; Landrigan et al., 2017 ; Thurston et al., 2017), suggesting potential ecological fallacy in our study when looking at the correlation with environmental factors. On the contrary, our results regarding socio-economic position are in line with results of previous research (see e.g. Bossuyt et al., 2004 ; Renard et al., 2019). We observe that low socio-economic position (share of unemployed, manual workers, young people with primary education only) is positively correlated with prevalence of medication reimbursement, with the exception of individuals coming from low-income countries. 

38Principal Component Analysis including medication reimbursement data, environmental and socio-economic variables indicates that there is no relationship between the prevalence of medication reimbursements and all other collected variables, confirming again the strong independency of such variables. The clustering analysis, including the five medication reimbursement groups and the three PCA factors computed on the environmental and socio-economic variables, presents a strong spatial structure. The addition of environmental and socio-economic variables does not have a relevant impact on the spatial pattern of medication reimbursements, as ASMP appear to be highly correlated between each other and independent from external factors. 

39The patterns of medication reimbursement in Belgium highlighted in this paper indicate that the lack of information at the individual level does not really allow us to measure health, and to draw conclusions on the effects of factors such as environmental characteristics or socio-economic position on wellbeing. Our data rather hints at medical and pharmaceutical practices and reminds us to keep in mind the existence of geographical differences in medical practice when using prescribed medication as a proxy for disease. Information on the education and training of the medical doctors and their current working location might be useful to further understand medication reimbursements distribution in Belgium. Moreover, “commercial” practices of pharmaceutical firms would be an interesting issue to understand further, especially in a country like Belgium divided by a linguistic border.   

40As a conclusion, this paper describes the spatial intra-national variations of medication prescriptions in Belgium. The use of data for an entire country and for several medication groups allows to generate hypotheses on the factors leading to medications purchase in Belgium for all age groups. The study suggests that medication prescription and reimbursement is not (only) an environmentally driven problem, at least at the level of the municipalities. Variation in medication prescription and reimbursement may be linked to curative healthcare matters and language/education of practitioners. Therefore, it is important to be careful when using and interpreting such data as health proxies: this paper reminds us how important it is to describe, control and understand the data before moving into more explanatory analyses.