Modelling antimicrobial prescriptions in Scotland: A spatio-temporal clustering approach

TL;DR

This paper develops a novel spatio-temporal clustering method using DBSCAN, enhanced with dynamic time warping and KDE-based spatial weighting, to analyze antimicrobial prescription patterns in Scotland, accounting for diverse population densities.

Contribution

It introduces an extended DBSCAN approach with dynamic time warping and KDE-based spatial weights, improving clustering accuracy in heterogeneous density environments.

Findings

Enhanced clustering performance over traditional methods

Effective identification of GPs' prescription behavior patterns

Actionable insights for antimicrobial stewardship

Abstract

In 2016 the British government acknowledged the importance of reducing antimicrobial prescriptions in order to avoid the long-term harmful effects of over-prescription. Prescription needs are highly dependent on factors that have a spatio-temporal component, such as the presence of a bacterial outbreak and the population density. In this context, density-based clustering algorithms are flexible tools to analyse data by searching for group structures. The case of Scotland presents an additional challenge due to the diversity of population densities under the area of study. We present here a spatio-temporal clustering approach for highlighting the behaviour of general practitioners (GPs) in Scotland. Particularly, we consider the density-based spatial clustering of applications with noise algorithm (DBSCAN) due to its ability to include both spatial and temporal data, as well as its flexibility to be extended with further variables. We extend this approach into two directions. For the temporal analysis, we use dynamic time warping to measure the dissimilarity between warped and shifted time series. For the spatial component, we introduce a new way of weighting spatial distances with continuous weights derived from a KDE-based process. This makes our approach suitable for cases involving spatial clusters with differing densities, which is a well-known issue for the original DBSCAN. We show an improved performance compared to both the latter and the popular k-means algorithm on simulated, as well as empirical data, presenting evidence for the ability to cluster more elements correctly and deliver actionable insights.

Tables

Table 1. Table 1 : Number of outliers and number of correctly and incorrectly clustered data points for each algorithm when tested on the generated data set. The proposed method has a higher number of correctly clustered elements than both the original DBSCAN and k 𝑘 k -means, and a lower number of outliers than the original DBSCAN.

	Original DBSCAN	Proposed method	$k$-means
Number of outliers	256	183	NA
Correctly clustered elements	2092	2175	1189
Incorrectly clustered elements	308	225	1211

Table 2. Table 2 : Number of GPs for each territorial healthboard in Scotland. Each of the fourteen healthboards is listed with their respective GP ‘headcount’. The number of GPs varies greatly between the different healthboards, from as low as 30 to as high as 1076.

Healthboard	Number of GPs	Healthboard	Number of GPs
Ayrshire & Arran	338	Borders	126
Dumfries & Galloway	120	Fife	280
Forth Valley	270	Grampian	533
Greater Glasgow and Clyde	1,076	Highland	398
Lanarkshire	441	Lothian	941
Orkney	36	Shetland	30
Tayside	395	Western Isles	32

Table 3. Table 3 : Variables in the aggregated data matrix after filtering. Each GP practice which prescribed Amoxicillin between October 2015 and September 2017 is represented by a table row using a unique practice code as an identifier. The table lists their postcode and the corresponding latitude and longitude values, as well as their individual time series of summed up Amoxicillin prescriptions per month.

v1	v2	v3		v26	v27	v28
Practice code	Postcode	201510	…	201709	Latitude	Longitude

Equations

$$N_{\epsilon }(p)=\{q\in D|\mathrm{dist}(p,q)\le \epsilon \}.$$

$$\begin{array}{r}{N}_{{\epsilon }_1}(p)=\{q\in D|\mathrm{dist}(p,q)\le {\epsilon }_1\}\text{ and}\\ N_{{\epsilon }_2}(p)=\{q\in D|\mathrm{dist}(p,q)\le {\epsilon }_2\}.\end{array}$$

$$\widehat{f}(\mathbf{y})=\frac{1}{nh^{d}}\sum^{n}_{i=1}K\Big{(}\frac{\mathbf{x}_{i}-\mathbf{y}}{h}\Big{)},$$

$${\int }_{-\infty }^{\infty }K(x)dx=1.$$

$$K(\mathbf{x},\mathbf{x}^{\prime})=\frac{1}{(2\pi h^{2})^{d/2}}e\bigg{(}-\frac{||\mathbf{x}-\mathbf{x}^{\prime}||^{2}}{2h^{2}}\bigg{)},$$

$$h=n^{(-1/(p+4))},$$

$$f(x)=\frac{L}{1+e^{-k\cdot (x-x_0)}},$$

$$f(d_i)=\frac{2}{1+e^{-k\cdot (d_i-\bar{d})}}.$$

$$w_{ij}=s_{ij}\cdot\Big{(}\frac{f(d_{i})+f(d_{j})}{2}\Big{)},$$

Full text

MODELLING ANTIMICROBIAL PRESCRIPTIONS IN SCOTLAND:

A SPATIO-TEMPORAL CLUSTERING APPROACH

Antonia Gieschena*

*Antonia.Gieschen@ed.ac.uk

Jake Ansella

*J.Ansell@ed.ac.uk

Raffaella Calabresea

*Raffaella.Calabrese@ed.ac.uk

Belen Martin-Barragana

Belen.Martin@ed.ac.uk

a University of Edinburgh Business School, 29 Buccleuch Place, Edinburgh, EH8 9JS , United Kingdom

Abstract

In 2016 the British government acknowledged the importance of reducing antimicrobial prescriptions in order to avoid the long-term harmful effects of over-prescription. Prescription needs are highly dependent on factors that have a spatio-temporal component, such as the presence of a bacterial outbreak and the population density. In this context, density-based clustering algorithms are flexible tools to analyse data by searching for group structures. The case of Scotland presents an additional challenge due to the diversity of population densities under the area of study. We present here a spatio-temporal clustering approach for highlighting the behaviour of general practitioners (GPs) in Scotland.

Particularly, we consider the density-based spatial clustering of applications with noise algorithm (DBSCAN) due to its ability to include both spatial and temporal data, as well as its flexibility to be extended with further variables. We extend this approach into two directions. For the temporal analysis, we use dynamic time warping to measure the dissimilarity between warped and shifted time series. For the spatial component, we introduce a new way of weighting spatial distances with continuous weights derived from a KDE-based process. This makes our approach suitable for cases involving spatial clusters with differing densities, which is a well-known issue for the original DBSCAN. We show an improved performance compared to both the latter and the popular k-means algorithm on simulated, as well as empirical data, presenting evidence for the ability to cluster more elements correctly and deliver actionable insights.

Keywords

OR in health services – Multivariate statistics – Cluster analysis – Spatio–temporal data

1 Introduction

The need to identify unwarranted variation and extreme behaviour in drug prescribing was highlighted recently in 2018 by Catherine Calderwood, the Chief Medical Officer for Scotland, in her report titled ‘Practising Realistic Medicine’: “Unwarranted variation […] is variation in healthcare that cannot be explained by need, or by explicit patient or population preferences. Recognising unwarranted variation is of vital importance because it allows the identification of: Underuse of higher value interventions – i.e. under treatment. Over use of interventions which should be used less frequently. Over use of interventions which may result in harm.” Healthcare Quality and Improvement Directorate (\APACyear2018)). Motivated by this, our empirical application leads to a better understanding of prescription behaviour in Scotland regarding its temporal trends, how it differs across regions and differences between GP practices. The areas showing different behaviour patterns can present evidence for distinct policy or communication strategies. Regions with higher demand for specific medications could have a higher rate of underlying health issues. Highlighting extreme behaviour in the form of outliers as well as extreme behaviour within clusters can also offer valuable insights into the behaviour of specific GPs or specific, relatively small-scale areas. In a broader sense, our findings can be used in connection with a disease recognition and, ideally, prevention system. The goal of our research is to inform more effective prescribing by identifying unwarranted variation, to minimise harm, and reduce waste for the National Health Service (NHS) in Scotland. The specific focus on Co-Amoxiclav or Amoxicillin is partially driven by the consideration of the UK 5 Year Antimicrobial Resistance Strategy 2013-2018, which had the aim of: “optimising prescribing practice through implementation of antimicrobial stewardship programmes that promote rational prescribing […]” Department of Health and Social Care (\APACyear2016).

Using data obtained from the NHS Scotland Open Data platform, we can explore the issue of antimicrobial prescriptions. The goal is to identify groups or outliers that reflect positive or negative prescribing behaviours. As the groups are unknown, this suggests the use of clustering. We expect antimicrobial prescriptions to follow a yearly cycle Durkin \BOthers. (\APACyear2018), calling for the need of considering the time element. Spatial aspects reflect the role of health and socio-demographic effects Kjærulff \BOthers. (\APACyear2016); Mölter \BOthers. (\APACyear2018). Our choice is therefore to explore the application of spatio-temporal clustering in this context. Specifically, we analyse the characteristics of our proposal on antimicrobial prescriptions in Scotland, by using a density-based spatio-temporal clustering algorithm to group together (cluster) general practitioners (GPs) with respect to their prescription volume behaviour over time and over space in Scotland. The use of spatial and spatio-temporal clustering techniques to analyse diseases or drug prescriptions has been subject to extensive discussion in the literature C. Anderson \BOthers. (\APACyear2016); Blangiardo \BOthers. (\APACyear2016). \citeAkatz2010 analyse prescriptions on an individual patient level whereas we use data aggregated at GP level. Antimicrobial prescriptions in particular have also been subject of research Petersen \BOthers. (\APACyear2007). When analysing the spatial distribution of diseases, other factors are also considered and can be included in our approach, for example socio-demographic data Kjærulff \BOthers. (\APACyear2016). \citeAmolter2018 analyse antibiotic prescription behaviour in England using a hotspot analysis, and connect their findings to the socio-demographic situation in these areas. Differences between urban and rural areas in the context of health have also been subject of research. \shortciteAcairns2011 use simulation to analyse first-responder schemes for cardiac arrests in Northern Ireland, discussing spatial complexities that come with it for example driven by differences between urban and rural areas.

Our reason to choose a clustering algorithm for this analysis lies in their many advantages when analysing complex data sets to search for groups, for example by eliminating the need for pre-labelled data and the pre-knowledge about classes. Clustering algorithms can be classified depending on the way they assign the data points to their groups, for example by partitioning, building a hierarchy in the form of a dendrogram, or looking for density Han \BOthers. (\APACyear2011). They all group together data points into clusters based on their similarity and the choice of this similarity measure differs depending on the data type. Many clustering methods, especially for spatial clustering, use distance measures to calculate the dissimilarity between objects. Some examples belonging to this group include the k-means algorithm (Han \BOthers., \APACyear2011, pp.451-454), the “Clustering Large Applications based on RANdomized Search” (CLARANS) algorithm Ng \BBA Han (\APACyear2002) and the density-based algorithm DBSCAN Ester \BOthers. (\APACyear1996). Density-based clustering algorithms, such as DBSCAN, use the distance measures to assess the density of points in a region to determine whether a cluster exists in that location. These algorithms, however, suffer from a drawback. Varying densities in spatial dimensions can lead to the algorithm not being able to identify clusters in both very dense and not dense areas. In a practical context, the ability to handle varying densities becomes relevant when comparing urban and rural regions due to different spatial prevalence across areas. This is the case, for example, in Scotland, where highly dense urban areas are concentrated in the South while the centre, the North and the islands remain quite unpopulated. In order to solve this, related research usually focusses on using a nearest neighbour approach to determine density, such as DECODE Pei \BOthers. (\APACyear2009) or the use of shared nearest neighbours Ertöz \BOthers. (\APACyear2003). \citeAertoez2003 propose to define the similarity between two points based on the number of shared nearest neighbours. For this solution, however, one needs to decide the number of nearest neighbours which should be considered by the algorithm. This choice can strongly impact the clustering results.

Clustering of time series often uses transformations of the time series by representing them in a stochastic way for example with Markov Models (model-based approaches) or, alternatively, features (feature-based approaches). Then subsequently the clustering is conducted with the extracted model parameters or features De Angelis \BBA Dias (\APACyear2014). Alternatively, for shape-based approaches, clustering algorithms compare the raw, untransformed, time series by trying to match them by stretching or shifting them Aghabozorgi \BOthers. (\APACyear2015). In all three cases, one needs to determine the measure of similarity between time series. This can be difficult, because effects may not be synchronised, but delayed or shifted. One example is the increase in the number of cases of a disease like the flu, which might occur in different regions at separate times. The challenge of shifted time series has been discussed in the literature and lead to the introduction of methods such as dynamic time warping (DTW) Serra \BBA Arcos (\APACyear2014). DTW calculates the dissimilarity of two time series by minimising the ’cost’ needed to match them. While it is a popular method in time series clustering, to the best of our knowledge it has not been used in connection with spatial data as input for the DBSCAN algorithm.

Our proposed method is able to handle both of these identified challenges, namely clusters of varying densities and shifted time series. In contrast to the nearest neighbour solutions in the literature Kriegel \BOthers. (\APACyear2011); Ertöz \BOthers. (\APACyear2003), our approach is continuous, covering points with a smooth function, without considering a specified number of neighbours. The main contribution of this work can be found in its methodological novelty of using continuous density estimates around spatial data points to weight their distances to one another by deriving weights through a logistic function. This idea of using the density measures derived via kernel density estimation (KDE) Silverman (\APACyear1986), which we introduce here in a DBSCAN algorithm, is very flexible. It can be implemented for a variety of density-based clustering algorithms which consider spatial information in the form of a distance matrix. Examples for this include, but are not limited to, DBSCAN Ester \BOthers. (\APACyear1996) and many of its extensions. In this paper, we implement our approach using ST-DBSCAN, a spatio-temporal version of the density-based clustering algorithm DBSCAN Birant \BBA Kut (\APACyear2007). By weighting the spatial distance matrix with continuous distance weights derived from a KDE-based process, we are able to make our approach suitable for spatial clusters with varying densities. The algorithm makes use of KDE in order to formulate an area’s density Terrell \BOthers. (\APACyear1992); Silverman (\APACyear1986). While previous work formulates a region’s density based on a point’s nearest neighbours, this approach provides us with a smooth density estimate in every location of our spatial map. The specific advantage of our approach lies in its increased suitability for real-life applications in which spatial clusters exhibit varying densities and the global neighbourhood radius used by DBSCAN to identify clusters is not accurate, as can be seen in our application. Furthermore, the methodological contribution of our work is usage of KDE for scaling which contributes to the area of local scaling and dealing with varying densities for clustering. This is of special interest for organisations and companies which are also looking for ways of graphically presenting their results, making them more accessible. Furthermore, regarding the temporal aspect, this approach is combined with a flexible way of measuring temporal dissimilarity by employing DTW. This results in a spatio-temporal clustering algorithm suitable for the analysis of time series which are spatially unequally distributed.

We test our approach against k-means as the most commonly used clustering algorithm, as well as against the original DBSCAN algorithm. Our results show an improved performance regarding the number of correctly clustered elements when compared with both these algorithms. When tested on 3000 simulated spatial data points, our proposed method reduces the number of incorrectly clustered elements from 1211/3000 (k-means) and 308/3000 (DBSCAN) to 225/3000. At the same time, the number of outliers is reduced from 256/3000 (DBSCAN) to 183/3000. Applied on spatio-temporal empirical data, the proposed method is able to cluster both a dense and a sparse area simultaneously, while the original ST-DBSCAN algorithm either over-smoothed them or identified a very high number of small clusters.

We provide background information on both clustering and KDE in Section 2, while Section 3 describes the proposed methodology in detail. Our approach is tested against the original DBSCAN to show its improvements, as well as against k-means as a common baseline algorithm, in a simulation in Section 4. In Section 5 our method is applied in a case study of GP prescription behaviour in Scotland to demonstrate its performance for spatio-temporal clusters with varying densities and shifted time series. The paper concludes and offers suggestions for further research in Section 6.

2 Literature review

Clustering is an unsupervised machine learning technique for grouping data points. In contrast to supervised classification, it finds and groups data structures according to their (dis)similarity instead of requiring a training set of pre-labelled examples Han \BOthers. (\APACyear2011); Kaufman \BBA Rousseeuw (\APACyear2009). Clustering applications are broad and range from biology Kiselev \BOthers. (\APACyear2019) to economics Dupor \BBA McCrory (\APACyear2017), and pattern recognition Schroff \BOthers. (\APACyear2015) to medicine Dougherty \BOthers. (\APACyear2002). The choice of algorithm depends on the type and amount of data, as well as on the objective. Well-known algorithms include partitioning methods such as k-means, which optimises an objective function calculating the intra-cluster and inter-cluster distances. There are numerous adaptations which, for example, respect upper and lower bounds for the number of clusters Borgwardt \BOthers. (\APACyear2017). Instead of finding one outcome, hierarchical approaches produce tree-like structures called dendrograms Han \BOthers. (\APACyear2011). Probabilistic approaches, including mixture models, assume that each cluster is one of a mixture of subpopulations Mai \BOthers. (\APACyear2018); Liao (\APACyear2005). Graph clustering looks for clusters of vertices in graphs in such a way that there are many edges within a cluster and relatively few between them Benati \BOthers. (\APACyear2017); Schaeffer (\APACyear2007). One group of algorithms in this area are spectral methods, where an eigenvector or a combination of eigenvectors is used for computing the similarity of clusters Nascimento \BBA De Carvalho (\APACyear2011); Schaeffer (\APACyear2007), while another approach is to formulate a min-sum problem which minimises the sum of distances between points in a cluster Hassin \BBA Or (\APACyear2010).

Spatial clustering describes using cluster analysis to cluster data points in geographic regions. For example, \shortciteAcoll2014 use hierarchical clustering to investigate patterns of vehicle collisions within an identified hot spot area in Northern Ireland. In the literature spatial clustering is sometimes treated as a special case due to the assumptions made in the wider area of, for example, spatial statistics, namely interdependency structure between locations Menafoglio \BBA Secchi (\APACyear2017); Cressie (\APACyear1992). This means that we assume that locations which are close influence each other more strongly than those further away from each other. There have been efforts to capture this dependency structure through copulas based clustering, which can be interpreted as a sub-category of model-based clustering in which copula information is used to derive the cluster composition Di Lascio \BOthers. (\APACyear2017). \shortciteAdisegna2017 use copulas to capture spatial dependency structures in their dissimilarity matrix, which they then feed into a fuzzy clustering algorithm. The latter does not allocate each object to one and only one cluster, but instead allows for membership in multiple clusters by assigning membership degrees (Everitt, \APACyear2011, pp.242-245). \shortciteAmarbac2017 use a mixture of Gaussian copulas for clustering mixed data, commenting on the advantage of being able to interpret the dependency structure through the copula.

Approaches commonly used for spatial data are density-based methods such as the DBSCAN algorithm Ester \BOthers. (\APACyear1996), which is explained in more detail in Section 3. This algorithm looks for high-density neighbourhoods of points which are similar or close to each other, and segregated by areas of lower density. Another method using the concept of neighbours is CLARANS Ng \BBA Han (\APACyear2002), which is able to cluster not only points but also polygonal objects. There are two main advantages of density-based methods. The first is their ability to find clusters of arbitrary shape. The second is the fact that many density-based methods do not require the number of clusters as an input as does, for example, k-means.

There are many adaptations to DBSCAN. These include the constraint-based C-DBSCAN Ruiz \BOthers. (\APACyear2007), which uses pre-knowledge about group memberships, sampling-based IDBSCAN Borah \BBA Bhattacharyya (\APACyear2004), which is especially suitable for very large spatial databases, and the spatio-temporal ST-DBSCAN Birant \BBA Kut (\APACyear2007). Another extension on the original DBSCAN algorithm is “Ordering Points to Identify the Clustering Structure” (OPTICS) Ankerst \BOthers. (\APACyear1999), which produces not a clustering but an ordering of the database that represents the underlying clustering structure. It is ordered in such a way that spatially close points are neighbours in the ordering, and objects in dense clusters are processed first. Besides being a good tool for visualisation purposes, OPTICS’s main advantage is its ability to handle varying densities by analysing the neighbours for each point. Similarly, “Locally Scaled Density Based Clustering” (LSDBC) Biçici \BBA Yuret (\APACyear2007) also uses neighbouring points, in this case to determine the local density maxima in regions which act as cluster centres. It aims to identify clusters within background noise and is thought to be more robust towards parameter changes. LSDBC connects density regions until the overall density falls below a pre-defined threshold. \citeArodriguez2014 propose “Clustering by Fast Search and Find of Density Peaks” (CFSFDP), which identifies cluster centres by defining them as points with a high local density and a high distance to other such centres. Subsequent work proposes approaches to identify the optimal values for these parameters Mehmood \BOthers. (\APACyear2016) and fuzzy implementations of the algorithm Bie \BOthers. (\APACyear2016).

Density-based clustering methods are applied to a broad range of research areas. \citeAgomide2011 use Twitter data to analyse a Dengue fever epidemic. The authors use ST-DBSCAN to cluster cities that show a similar incidence rate of Dengue fever occurring during the same period of time, exploring whether Twitter can be used to predict the number of Dengue fever cases. For this purpose they compare the number of cases reported by official statistics and the number of tweets about the disease during the same time period. This enables them to develop a Dengue fever surveillance tool based on four dimensions, namely volume, location, time, and public perception, showcasing the real life applicability of their approach. In addition, \shortciteAanbaroglu2014 comment on the use of spatio-temporal clustering to detect non-recurrent traffic congestion in London, as clusters of slowly-moving vehicles can cause congestion in an urban road network.

DBSCAN and adaptations like ST-DBSCAN face problems when trying to detect clusters with varying point densities due to the use of global parameters when looking for close neighbouring points that indicate existing clusters. An example for a situation like this can be seen in Scotland, where the population density varies across the region. This issue has been addressed in the literature, for example by \shortciteAkriegel2011. They explain how density-based clustering methods often label points in less dense areas as noise or outliers when there are regions of high density separated by areas of lower density. As potential solutions to this problem, they propose hierarchical approaches of ordering points such as OPTICS Ankerst \BOthers. (\APACyear1999), as well as nearest neighbour approaches Ertöz \BOthers. (\APACyear2003); Pei \BOthers. (\APACyear2009).

\citeA

birant2007 discuss varying densities in their paper on ST-DBSCAN and outline the problem of identifying actual noise in a scenario of high-density areas surrounded by regions of lower density. They propose assigning each cluster a density factor, calculated by looking at the maximum and minimum distances within a respective cluster. This single factor might, however, not reflect the various different density regions which are present across each cluster, especially if the area described by one factor is relatively large.

Similarly, \citeAzelnik2005 discuss the idea of local scaling in the context of their self-tuning spectral clustering algorithm. Based on this, \citeAbiccici2007 combine this local scaling with an adapted density-based clustering algorithm. Their density estimation, however, is based on a $k$-nearest neighbour approach. There are two aspects to this approach which can be considered drawbacks. First, one has to decide the number of considered neighbours $k$, which is not the case for our approach. Secondly, in contrast to our method, the approach by \citeAbiccici2007 does not produce a continuous density estimation map.

Clustering can also be employed in time series analysis De Angelis \BBA Dias (\APACyear2014); Serra \BBA Arcos (\APACyear2014). For clustering such data, emphasis in the literature has been put on the selection of the most suitable dissimilarity measure. In general, there are four ways of measuring the dissimilarity between time series, as reviewed by \citeAserra2014. Lock-step methods such as the Euclidean distance measure similarity by comparing values at the same time steps in each series, their main advantage being computational simplicity Xing \BOthers. (\APACyear2012). Feature-based methods employ the same approach, but work on a transformation of the time series. This can, under some circumstances, improve the accuracy of the dissimilarity compared to the raw Euclidean distance and is, therefore, commonly used Serra \BBA Arcos (\APACyear2014); Montani \BOthers. (\APACyear2006). \citeAinniss2006 introduces an approach for seasonal clustering of time series applied to weather and aviation data. They form clusters of contiguous months with similar weather conditions, using the mean value of a month, as measured over several years, or the difference between distributions for each month, over a number of years. Model-based approaches assume an underlying model generating the time series and measure dissimilarity based on the model parameters. For example, mixture models are used by \shortciteApovinelli2004, who comment that the method is able to perform well on a range of application areas with minimal input tuning. \shortciteAdias2015 present an approach which is able to account for time-constant unobserved heterogeneity and hidden regimes within time series by using hidden Markov models. Elastic methods compare the time series by their overall shape over time instead of their absolute values measured at each time step. Dynamic Time Warping (DTW) aims to optimally align the two time series by minimising a warping cost. It has been used in various contexts, including for incomplete medical time series Tormene \BOthers. (\APACyear2009) and to measure the dissimilarity between time series in combination with fuzzy clustering Izakian \BOthers. (\APACyear2015). “Edit Distance on Real Sequence” (EDR) calculates dissimilarity as the number of edit operations necessary for transforming one time series into the other Chen \BOthers. (\APACyear2005). $k$-Shape and $k$-MS are iterative time series clustering algorithms which use shape-based distance (SBD) as their dissimilarity measure Paparrizos \BBA Gravano (\APACyear2017). The authors comment on the fact that their algorithm is able to outperform DTW, especially regarding computational time and the fact that their approach requires no parameter tuning.

This paper makes the following two contributions. The first is related to DBSCAN’s shortcomings when confronted with varying point densities. Existing literature has tackled this problem with nearest neighbour approaches Kriegel \BOthers. (\APACyear2011); Pei \BOthers. (\APACyear2009); Biçici \BBA Yuret (\APACyear2007); Zelnik-Manor \BBA Perona (\APACyear2005); Ertöz \BOthers. (\APACyear2003) or by assigning density factors to clusters Birant \BBA Kut (\APACyear2007). To the best of our knowledge, no continuous approach using KDE has been used for this purpose so far. The second gap is the usage of a time series dissimilarity measure which is able to handle shifted or warped times series. Here, we tackle this by employing DTW with the spatio-temporal clustering algorithm ST-DBSCAN.

3 Methodology

The objective of this paper is to enable the spatio-temporal clustering of a large amount of data points forming clusters with varying densities, each of which is associated with a univariate time series. Our approach can be used in connection with a variety of clustering algorithms but will be demonstrated by ST-DBSCAN allowing spatio-temporal clustering Birant \BBA Kut (\APACyear2007). This algorithm is chosen due to its computational efficiency and its flexibility with regard to the addition of more variables alongside the temporal and spatial information.

In order to tackle ST-DBSCAN’s inability to handle unevenly distributed spatial points, our approach manipulates the input spatial distance matrix using weights based on the location density derived through KDE. This makes the spatial distances of locations with varying point densities comparable to each other. In doing so, the use of a global spatial radius parameter, which identifies other cluster members belonging to the same group, becomes possible.

3.1 Clustering using DBSCAN and ST-DBSCAN

DBSCAN operates by visiting a random data point and looking for neighbouring points within a pre-defined maximum spatial radius $\varepsilon_{1}$. \citeAester1996 define the $\varepsilon$-neighbourhood of a point $p$ as $N_{\varepsilon}(p)$, with

[TABLE]

If $N_{\varepsilon}(p)\geq\eta$, that means the number of neighbouring points in the neighbourhood described by Equation (1) is more than a pre-defined number $\eta$, the respective area is considered a dense area, and a cluster is formed. In that case $p$ is called a core point. Points which are subsequently assigned as members of this cluster without fulfilling the condition are called border points. Points outside of clusters are labelled as noise.

In order to form the clusters, \citeAester1996 define a number of concepts regarding the reachability of points. To capture both border and core points, these concepts are directly density-reachable, density-reachable and density-connected.

A point $q$ is directly density-reachable (DDR) from a point $p$ if two conditions hold true, namely

$q\in N_{\varepsilon}(p)$ and 2. 2.

$|N_{\varepsilon}(p)|\geq\eta$,

with the second condition being called the core point condition. This means that $q$ must be in the $\varepsilon$-neighbourhood of $p$, and $p$ must be a core point. This is visualised in Figure 1.

A point $o$ is density-reachable (DR) from a point $p$ if there exists a chain of points which are directly density-reachable from each other. This condition makes DBSCAN flexible for clusters of arbitrary shape due to a chaining effect which is, however, constricted to points which fulfil the two conditions for direct density-reachability. This chaining effect is shown in Figure 2.

A point $t$ is density-connected (DC) to a point $p$ if there is a point $o$ which is density-reachable from both $t$ and $p$. This concept is shown in Figure 3.

Based on these definitions, \citeAester1996 define a cluster $C$ as a non-empty subset of a database $D$ that fulfils the two conditions

$\forall p,q$: if $p\in C$ and $q$ is density-reachable from $p$ then $q\in C$, and 2. 2.

$\forall p,q\in C$: $q$ is density-connected to $p$.

ST-DBSCAN operates in a similar manner as DBSCAN, but in addition to the spatial parameter for determining the $\varepsilon$-neighbourhood, $\varepsilon_{1}$, the algorithm takes an additional parameter $\varepsilon_{2}$, which determines the temporal neighbourhood. In order to be considered a member of a cluster, a data point $q$ has to be in both the spatial and the temporal neighbourhood of point $p$. This extends Equation (1) as illustrated in Equation (2).

[TABLE]

DBSCAN and ST-DBSCAN both look for neighbouring points for each data point separately, which means that the distance matrices are calculated within the main algorithm body. If one can assume that the spatial distances are stable over time, and that the time series are not updated very frequently, it makes sense to take this matrix calculation outside of the neighbourhood calculation. The two matrices are calculated and stored before running DBSCAN or ST-DBSCAN, after which the algorithm then just has to access the stored information relevant for the respective visited data point.

Both DBSCAN and ST-DBSCAN suffer from the drawback of choosing $\varepsilon_{1}$ for $N_{\varepsilon_{1}}(p)$. Note that $\varepsilon_{1}$ is global, in the sense that it is constant across the whole region of the study. If the spatial points are distributed with varying densities, as it is, for example, the case in our application in Scotland, then one $\varepsilon_{1}$ is not be able to capture structures in both dense and sparse areas. If $\varepsilon_{1}$ is chosen too small, it will not be able to identify areas of relatively high density in less dense areas. If it is chosen too large, it will not be able to distinguish between structures in high-density areas. This drawback is the primary motivation for the development of our method. We consider two possible directions to solve this issue. The straightforward approach is to use a distance that is modified according to the density of the area. Alternatively, the global $\varepsilon_{1}$ can be replaced by a local parameter that varies according to the area’s density.

3.2 Kernel density estimation

In order to make more and less dense areas comparable and, therefore, a global parameter suitable, we employ KDE as a non-parametric way of formulating an area’s density. KDE is a density estimation technique which approximates the probability density function from samples. It takes the following form Terrell \BOthers. (\APACyear1992); Matioli \BOthers. (\APACyear2018):

[TABLE]

for a random sample $\mathbf{x}$ of size $n$, with $\mathbf{x}\in\mathbb{R}^{d}$, $h$ as a smoothing parameter, and $K$ as the kernel function which satisfies the condition

[TABLE]

KDE is used to estimate the density function from observed data. As a non-parametric approach, it assumes that the data originates from a distribution with the probability density $f$, and then directly uses the data to estimate it. For this, a suitable kernel function is selected, which models the influences of the data points in the data space, with the Gaussian kernel as a commonly used option Hinneburg \BBA Gabriel (\APACyear2007). The sum of all kernels at a location then provides an estimate of the density function at this location. A Gaussian kernel function is described as

[TABLE]

where $h$ is the bandwidth, $d$ the dimensionality of $\mathbf{x}$, and $||\mathbf{x}||=\sqrt{\mathbf{x}^{T}\mathbf{x}}$ denotes the Euclidean norm Sugiyama (\APACyear2015).

Besides selecting the kernel function, there is also the need to select the bandwidth, or smoothing factor, $h$, which affects the smoothness of the surface laid over the points. In other words, it is the influence a data point has on more distant regions , and can be held constant or made variable across the data space Hinneburg \BBA Gabriel (\APACyear2007); Terrell \BOthers. (\APACyear1992). $h$ is often chosen by cross validation or ”Scott’s Rule’ Scott (\APACyear1992).

An advantage of KDE is its continuous nature, as it approximates the density function at all locations of the data space, not just at the locations with existing data points. This makes our approach also suitable for situations in which new data points are added in between existing ones. KDE is often used to identify hotspots which are areas of high point density. \shortciteAanderson2009 identifies hotspot areas of road accidents with a combination of KDE and $k$-means. They create a density map which is subsequently divided into hotspot cells, with the latter then being clustered using $k$-means. \shortciteAsmith2018 use KDE to create a map of drug-resistant tuberculosis cases in South Africa. \citeAgerber2014 develops a combined approach of spatio-temporally tagged Twitter data and KDE in order to predict crime in a US city, whereas \citeAeslinger2017 analyse areas of high numbers of vehicle-bear collisions in Florida. \citeAmatioli2018 introduce ClusterKDE, an algorithm based on univariate KDE combined with optimisation techniques, to iteratively obtain new clusters by minimising a smooth kernel function.

3.3 ST-DBSCAN with KDE-based local scaling

Let $\mathbf{x}$ be a set of observations in a two-dimensional space. Using these observations, we estimate the density function with KDE as described in Section 3.2. For the selection of the bandwidth $h$, there are a number of optimisation techniques available. We use “Scott’s Rule” Scott (\APACyear1992), which calculates $h$ as

[TABLE]

with $n$ as the number of data points and $p$ the number of dimensions. Scott’s Rule has previously been shown to perform well compared to other approaches Scott (\APACyear2009), and is frequently used as a way of determining the optimal bin width, for example by \citeAbernacchia2011. Let $d_{i}$ be the density of point $i$, calculated using Equation (5). Let $\bar{d}$ be the average density of all points in the dataset.

Our aim is to re-scale distances in such a way that the result is a distance matrix in which entries for more and less dense areas are comparable. Let $S$ be the spatial distance matrix, then the elements of $S$ are $s_{ij}$ for the distance between points $i$ and $j$, calculated for example using a great-circle distance, which is the shortest distance between two points on a sphere such as the Earth Williams (\APACyear2011). We then re-scale the distance matrix with weights derived from the density estimates. This means that the resulting scaled distance matrix includes information about the density in the location surrounding an object.

A logistic function is chosen for the calculation of the weights from the density estimates. The reasons for this is that the slope of the function can be adjusted by changing a parameter, which makes it possible to change the impact the density estimate has on the weights. Other reasons for choosing a logistic function are its continuous nature and the fact that it can be symmetrically centred around one. The function takes the general form

[TABLE]

with $x_{0}$ as the horizontal axis value of the sigmoid’s midpoint where the curve changes from convexity to concavity, $L$ as the sigmoid’s maximum value, and $k$ as the steepness of the curve.

Our aim is to use the function shown in Equation (7) as a way of “translating” density estimates into weights. In order to stretch distances in dense areas and compress them in sparse areas, we use a function that increases or decreases weights while converging to a chosen maximum value. Therefore, in our approach, the curve’s midpoint $x_{0}$ is chosen to be the mean of the density estimates $\bar{d}$ and the curve’s maximum $L$ is chosen to be two, thus centring the curve around one on the vertical axis, which gives us for the density $d$ of point $i$

[TABLE]

This means that if the density estimate for an area is equal to the overall average of all areas, these distances are multiplied with a weight of one, resulting in no modification. The two extreme ends represent the most dense and least dense areas, where the weights are two and zero, respectively. We suggest to weight the distance between two points with a factor of two, the maximum weight, in the highest-density area of the space. We choose the value two as a maximum for symmetry reasons, to account for zero being the minimum and one being the mean weight. The weight of exactly zero is only used to scale the distance between spatially identical points, in which case the distance is already zero. The steepness $k$ of the curve can be adjusted depending on the application.

The distance matrix is now modified with

[TABLE]

for a point pair $i$ and $j$, and their respective distance $s_{ij}$ in the spatial distance matrix $S$. Equation (9) calculates the weights for both $i$ and $j$ using Equation (8), takes the mean of these two weights, and multiplies this combined weight with the respective distance $s_{ij}$.

4 Simulation

In order to test the properties of the proposed method, we first apply it to a simulated dataset. The main purpose of this simulation study is to determine whether our approach is able to handle data with varying point densities between different areas, with an improvement in the results when compared to the original algorithm.

4.1 Set up of the simulations

For a more straightforward implementation, as we want to highlight the extent to which our proposed distance weighting improves the outcome, we only consider spatial variables. For this purpose, we implement our approach in combination with the original DBSCAN algorithm. An implementation with both spatial and temporal dimensions using ST-DBSCAN can be found in Section 5.

We generate 2400 data points in order to have sufficiently large individual clusters to be detected with our method while being sufficiently small to allow for fast runtimes. They consist of four clusters with 600 randomly Gaussian-distributed points per cluster. More specifically, each cluster consists of 300 points that are randomly sampled from a Gaussian with a higher variance, and 300 points randomly sampled from a Gaussian with a smaller variance and the same mean. The clusters can be seen in Figure 4.

All four datasets are clustered using the original DBSCAN approach, the proposed method of DBSCAN with distance weighting, and $k$-means as a benchmark approach.

The chosen methods allow us to compare our approach with the original DBSCAN algorithm proposed by \citeAester1996 to demonstrate its improved performance under certain conditions. Furthermore, the comparison with $k$-means acts as a benchmark, as $k$-means is considered one of the most common clustering algorithms that also uses distance measures to form clusters (Han et al., 2011, pp.451-454).

The results, as displayed in Figure 5, are compared based on their ability to identify the more dense clusters while still accounting for any clusters in the less dense areas surrounding them. For assessment purposes, the results are compared regarding their number of outliers and the number of correctly classified points111All three methods are implemented using R version 3.5.1 and RStudio version 0.99.896. For the original DBSCAN method, we use the R package ‘dbscan’ version 1.1-2, for $k$-means we use the R package ‘stats’ version 3.4.4. Our proposed method uses KDE from the package ‘ks’ version 1.11.0. The code is available from the authors upon request..

4.2 Results of the simulations

We aim to choose parameters for all three methods using common approaches, comparing them under individually realistic conditions. This results in DBSCAN with $\varepsilon$ = 0.3, chosen according to the nearest neighbour plot as proposed by \citeAester1996, and $\eta$ = 6. For the proposed adapted DBSCAN, we choose $\varepsilon$ = 0.3 and $\eta$ = 6 based on the same principle. For $k$-means, we choose $k=4$ based on both the known real number of clusters and the within-clusters sum of squares (‘elbow method’).

The initial assessment is a general visual comparison of the results as shown in Figure 5. As the data is generated to represent four clusters, we are able to assess whether the three methods correctly identify these four groups.

Next, we compare the number of outliers, as well as the number of correctly clustered objects as shown in Table 1. As expected, $k$-means does not perform well on the simulated data set, combining the two lower clusters into one and splitting the top right cluster into two. The reason for this can likely be found in the random initial seeds, which strongly influence the clustering outcome. The number of outliers is not applicable (NA) for this method, as $k$-means assigns each data point to a cluster without marking outliers.

The results for DBSCAN and the proposed method show an improvement for the latter regarding the number of correctly and incorrectly clustered objects. Looking closely at the graphical outputs shows that the proposed method is able to capture more of the less densely-spread data points surrounding the denser areas. Due to its ability to adapt to changing environments, the algorithm’s result has a lower number of outliers because the search radius parameter varies with changing point densities. This makes the method especially valuable when outliers are driven due to environmental influences such as more rural or remote areas, and the researcher wishes to account for this fact.

5 Empirical application

Our proposed method is tested for its applicability by analysing prescription data on GPs from the Scottish NHS. The chosen data was collected over a two-year period from October 2015 to September 2017. There are four main objectives for this application. The first is to improve the general understanding of GP behaviour over time as well as differences between GPs in different locations across Scotland. As a further development, the usage of clustering also enables the organisation to generate representative samples of GPs. Sampling from the generated clusters means that the existing underlying structures can be replicated. Additional analysis of GPs can then be conducted using the samples instead of the whole dataset, improving the speed of analysis significantly and allowing for real-time analysis and visualisation tools. Lastly, the developed software tool, which was built based on our methodology, makes our approach accessible for users from different areas to improve their data understanding in an interactive and real-time manner.

5.1 Data source and description

The data sources are NHS Scotland publications and the Open Data Platform of their Information Services Division (ISD)222https://www.opendata.nhs.scot/. Using open data has multiple advantages, such as being readily accessible via the ISD’s website which also makes our research easy to replicate.

NHS Scotland covers 14 territorial healthboards of Scotland as shown in Figure 6 and 4,994 GPs as of September 2018 NHS ISD (\APACyear2018). The number of GPs varies across the different healthboards, as seen in Table 2.

The prescription data from NHS ISD is merged with additional open data containing the GPs’ locations given by their postcode and corresponding latitude and longitude. We filter the data table for one specific drug, Amoxicillin, which is an antibiotic that is used to treat bacterial infections. The time series for each GP is then calculated by aggregating the prescription instances for different types of the drug on a monthly basis. Based on a two-year period from October 2015 to September 2017 with monthly data, this gives us a time series of length 24 for each GP. Additional variables are a unique identifier (practice code), post code, latitude, and longitude, which account for four more columns. There are 1081 rows in total, which includes all GPs who prescribed any type of Amoxicillin over the period of two years. Due to missing values in their coordinates, a number of rows had to be removed, resulting in a 980$\times$28 matrix.

When plotting the locations of the 980 analysed GPs using their latitude and longitude, as shown in Figure 7, it becomes clear the distribution of their locations is not even across the country. This is likely to be the case due to the uneven population distribution in Scotland. The “belt” of more densely populated areas towards the South is clearly visible, as are the much less populated northern regions.

Figure 8 shows the mean sum of prescriptions over all GPs per month. As can be seen, more Amoxicillin is prescribed during the late autumn and winter months (October to January) compared to the summer months (May to August).

5.2 Clustering process and results

We use two algorithms to cluster our data. In a first step, we use the original ST-DBSCAN algorithm. The spatial distance matrix is calculated using the ‘Meeus’ great-circle distance from the R package ‘geosphere’ Hijmans (\APACyear2017). The temporal dissimilarity between the time series is calculated using the DTW function ‘DTWDistance’ from the R package ‘TSDist’ Mori \BOthers. (\APACyear2016). These two dissimilarity matrices are then used by the algorithm to define $\varepsilon$-neighbourhoods as described in Section 3.1.

The time series are clustered by calculating the temporal dissimilarity matrices using DTW. Due to the expected seasonal variations in the prescription levels of Amoxicillin, DTW was chosen as it is able to detect similarity between time series where peaks in prescription volume occur shifted. This seasonality is shown in Figure 8. The ST-DBSCAN algorithm calculates the dissimilarity matrix outside of the main algorithm body and refers back to it during the $N_{\varepsilon_{2}}$ search. $\varepsilon_{2}=1000$ is chosen using the nearest neighbour distance plot as proposed by \citeAester1996.

The uneven spatial data distribution does, however, leads to a problem. The global spatial parameter $\varepsilon_{1}$ is not able to capture small-scale clusters in the South without oversmoothing them. If chosen too small, the algorithm will form many (usually around eleven) very small and dense clusters, with many points within less dense areas marked as outliers. Alternatively, the algorithm will select all points in the denser areas to belong to the same cluster when the chosen parameter is too large, as seen in Figure 9.

The overall goal of our proposed approach is for the algorithm to not overemphasise clusters in dense regions such as cities. This thought leads to the idea of scaling the distances based on their surrounding point density. For this, we use the process described in Section 3. We calculate the density estimates for all regions based on the amount of data points present in each area. Using these estimates as a way of determining how dense or sparse a geographic region is, we derive weights with which we multiply our spatial distance matrix. The updated distance matrix is then used as an input for the original ST-DBSCAN algorithm.

The results for our implementation show a more realistic clustering result, which is more useful for organisations working on such problems, as seen in Figure 10. For both the original ST-DSBCAN algorithm and our proposed method, we choose $\varepsilon_{1}$ = 8, $\varepsilon_{2}$ = 1000, and $\eta$ = 7 according to nearest neighbour distance plots as proposed by \citeAester1996. With the same parameters, our adapted ST-DBSCAN method typically captures around six clusters in the belt area compared to the eleven small ones with the original method. Overall, our approach identifies a smaller number of nine clusters, compared to 17 for the original ST-DBSCAN algorithm, as it is able to combine some smaller dense areas with their surrounding points.

Our proposed approach is able to highlight extreme behaviour such as particularly high or low prescription volumes of single or groups of GPs. This makes the results especially interesting for policy makers and the NHS as a whole, as unusually high prescription levels of a specific drug could be an indicator for an underlying problem in this region. Figure 10 shows, however, that many points are still classified as outliers. This can be explained due to the chosen variables, meaning prescription volume over time and location, not being informative enough to fully cluster GPs. Additional variables such as socio-economic factors can potentially add more useful information to this problem. Figure 11, which shows the average prescription amounts in each cluster and their seasonal changes over the analysed time period, illustrates this issue. Each cluster has a very similar time series, which mostly differs in total prescription volume and not, as one would expect, in the point in time where peaks occur. This exemplifies the need for further research, for example by analysing a different drug or adding additional factors.

6 Conclusion

The well-known DBSCAN algorithm Ester \BOthers. (\APACyear1996) and its spatio-temporal adaptation ST-DBSCAN Birant \BBA Kut (\APACyear2007) are flexible and efficient ways of clustering spatial and spatio-temporal data Catlett \BOthers. (\APACyear2019); Chimwayi \BBA Anuradha (\APACyear2018); Gomide \BOthers. (\APACyear2011). Due to its nature, DBSCAN suffers from a weakness when presented with clusters of varying densities. This is frequently the case when using DBSCAN for real-world spatial data. Previous research has emphasised this issue and tried to solve it using nearest neighbour approaches Pei \BOthers. (\APACyear2009); Biçici \BBA Yuret (\APACyear2007); Zelnik-Manor \BBA Perona (\APACyear2005); Ertöz \BOthers. (\APACyear2003). These approaches do, however, not result in a continuous representation and handling of the varying densities.

When clustering time series, the choice of dissimilarity measure is crucial. In order to cluster shifted time series, elastic measures such as DTW offer a suitable solution and perform better than step-wise methods. To the best of our knowledge, ST-DSBCAN has not been used in combination with DTW before.

In this paper, we present a method of scaling distances using continuous density-based distance weights. Our objective is to introduce a method of manipulating a spatial dissimilarity or distance matrix in such a way that it reflects the density of points in a given area. For this purpose, we propose and implement a way of deriving weights from the density estimates acquired through KDE, using a logistic function. Our proposed method first estimates the density in each spatial region, which is then used to weight the distances between data points depending on how dense or sparse their respective areas are. In denser areas, a higher weighting factor leads to the stretching of a distance, whereas lower weights for less dense areas lead to the compression of a distance.

We demonstrate the applicability of this approach for clustering with density-based methods such as DBSCAN, especially in cases of uneven data point densities in which a global distance parameter might not be able to capture group structures in regions of varying density. For this purpose, we show how our approach works well in both experiments with simulated data and an application case of NHS Scotland prescription data. The second application case also establishes that our approach works well in combination with using DTW for the temporal dimension. In addition, our method is very flexible and can be used in connection with any clustering algorithm that uses a distance matrix as an input.

We see the methodological contribution of our work in the use of distance scaling based on local densities using KDE. The inability of DBSCAN to handle varying point densities has been acknowledged by some of the algorithm’s original authors in \citeAkriegel2011. Previous attempts to tackle this issue focus on nearest neighbour approaches Pei \BOthers. (\APACyear2009); Biçici \BBA Yuret (\APACyear2007); Zelnik-Manor \BBA Perona (\APACyear2005); Ertöz \BOthers. (\APACyear2003) or assigning density factors to clusters Birant \BBA Kut (\APACyear2007). As for the reason behind these developments, \citeAbiccici2007 comment on their choice of a neighbour-based density estimate instead of a Parzen window approach by highlighting the difficulty of choosing the right window size. This issue has, however, been addressed extensively in the literature covering KDE. Our approach also offers added benefits due to the use of a continuous density map, which is not only useful for measuring the density for any location or new point, but also has the advantage of good visualisation possibilities. In contrast to nearest neighbour approaches, our proposed method does not require a specification of the number of neighbours considered to assess an area’s density. We thus see our algorithm as an improvement, especially in terms of interpretability and visualisation. Its flexibility means that it can be applied to any case in which data points are spatially clustered with an algorithm that uses a distance matrix as an input.

Furthermore, with regard to the application case and practical contribution, our research offers actionable insights into drug prescription behaviour. Our findings demonstrate how clustering can be used to group GPs according to their spatial location and behaviour over time. Our proposed method is able to highlight extreme behaviour such as GPs deviating from peers in their group. This has been acknowledged as an important issue by both the NHS and, in the special case of antibiotics, the British government Healthcare Quality and Improvement Directorate (\APACyear2018); Department of Health and Social Care (\APACyear2016), showcasing the increasing awareness of policy makers.

We see the potential for future research in exploring how the smoothing factor in KDE affects the results and can be adjusted depending on the application case. When interpreted in the context of our approach, the smoothing factor accounts for the amount of impact a group of data points has on their surroundings, or how far their “density-impact” reaches spatially. From an application point of view, this should be considered by researchers on a case-by-case basis.

In addition, our application to prescription data has implications for the analysis of prescription behaviour in different areas and over time. While promising, our findings also provide an example for the issue of trying to explain complex dynamics using only two variables. The fact that a relatively common antibiotic has been chosen as the drug to be analysed might also contribute to the issue of relatively little spatial fluctuation. Further research could focus on additional explanatory variables which can shed more light on this behaviour, for example socioeconomic factors, as additional dimensions besides time and space. Research similar to this proposal is, for example, introduced by \citeAkjaerulff2016. This will help to further expand the understanding of prescription behaviour for different medications and in different regions, offering valuable insights for organisation such as the NHS and policy makers alike.

Acknowledgements

The authors thank Wallscope for the support provided with the data access.

Funding: This work was supported by The Data Lab [project registration number Reg-17669].