Capturing spatial dependence of COVID-19 case counts with cellphone mobility data.

Spatial dependence is usually introduced into spatial models using some measure of physical proximity. When analysing COVID-19 case counts, this makes sense as regions that are close together are more likely to have more people moving between them, spreading the disease. However, using the actual number of trips between each region may explain COVID-19 case counts better than physical proximity. In this paper, we investigate the efficacy of using telecommunications-derived mobility data to induce spatial dependence in spatial models applied to two Spanish communities' COVID-19 case counts. We do this by extending Besag York Mollié (BYM) models to include both a physical adjacency effect, alongside a mobility effect. The mobility effect is given a Gaussian Markov random field prior, with the number of trips between regions as edge weights. We leverage modern parametrizations of BYM models to conclude that the number of people moving between regions better explains variation in COVID-19 case counts than physical proximity data. We suggest that this data should be used in conjunction with physical proximity data when developing spatial models for COVID-19 case counts.

Introduction

Spatial analyses of COVID-19 case data were first published as early as March of 2020 (Huang et al., 2020, Arab-Mazar et al., 2020, Giuliani et al., 2020), in an attempt to characterize, predict, and attenuate the severity of the pandemic. Subsequent studies have noted substantial spatial dependence in COVID-19 case counts (Kang et al., 2020, Bilal et al., 2020). This makes sense as regions that are close to each other likely have more people moving between them, spreading the disease to nearby regions.

Many groups have attempted to model COVID-19 case counts as a function of climate (Liu et al., 2020, Shi et al., 2020, Briz-Redón and Serrano-Aroca, 2020), healthcare quality (Sugg et al., 2021), socioeconomic factors (Baum and Henry, 2020) and more. More recently, mobility data has become more abundant and popular for modelling COVID-19 transmission. This makes sense because the disease spreads through human contact, meaning that case counts are likely to be a function of the number of people moving around. Such mobility data has been used to model the evolution of the epidemic in Spain (Aràndiga et al., 2020, Iacus et al., 2020), assess the effectiveness of the Spanish lockdown (Orea and Álvarez, 2020), monitoring the epidemic in Switzerland (Persson et al., 2021), identify at-risk populations in France during a lockdown (Pullano et al., 2020), individual-level infection tracing in China (Kraemer et al., 2020), assess the timing of stay-home orders (Audirac et al., 2020), and evaluate the effectiveness of social distancing in the United States (Badr et al., 2020). This data can be found in many forms, but is commonly found in the form of aggregated areal mobility matrices. If we denote a mobility matrix , corresponds to the number of trips from region to region , and represents the number of trips within region .

These data have been applied in a variety of different models to answer numerous questions, but lack of available methods makes it difficult for researchers to use this data to its full potential. In this paper, we demonstrate a novel method for analysing this data, whereby the mobility data is used as edge weights in a Gaussian Markov random field (network) model. Previous work using network models have been applied to mobility data in the form of a network compartment model (Chang et al., 2021) which was used to conduct inference regarding societal inequities, and inform reopening. This work does not aim to make such claims, but rather demonstrate the efficacy of mobility data in modern parametrizations of Besag, York, and Mollié (BYM) models (Besag et al., 1991) and their extensions.

BYM models have been used frequently in the spatial analysis literature due to their effectiveness and computational efficiency. In these models, the spatial component is comprised of Conditional Autoregressive (CAR) (Besag, 1974) models and conventional random effects. This means that the spatial effect of region depends only on its “neighbours”. Neighbours could be defined by any quantity the analyst has access to, but is most often defined by physical adjacency, i.e. if two regions share a common border, they are considered neighbours. Several ICAR/BYM models have been applied to COVID-19 data with neighbours defined in this way (DiMaggio et al., 2020, Huang and Brown, 2021, Brainard et al., 2020). Although these spatial model components based on physical adjacency are powerful and computationally efficient, it makes more sense to use mobility between regions to induce spatial dependence in COVID-19 models because the disease spreads via person-to-person contact.

In this paper, we build a BYM model where mobility data is used to induce spatial dependence between regions. Using mobility data within two Communities in Spain, Madrid and Castilla-Leon, we demonstrate the value of mobility data for COVID-19 spatial modelling applications. Furthermore, we extend modern parametrizations of BYM models to account for both physical adjacency and mobility simultaneously, and show that mobility data captures spatial variation in COVID-19 case counts much more accurately than physical adjacency alone.

This is a short focused paper with the following plan. Section 2 presents the data and the modelling strategy based on particular parametrizations of BYM models. The results come in Section 3, and the paper ends with a final discussion in Section 4.

Methods

Data

This paper is focused on two regions in Spain. Castilla-Leon is the largest Community in Spain by area and is located in the northwest part of Spain, with a population of 2.5 million. The Community of Madrid is located in the central part of Spain and has a population of around 6.8 million, and it is home of the capital of the country, Madrid City, with 3.3 million inhabitants.

The human mobility data was obtained from Barcelona Supercomputing Center Flow-map dashboard (Valencia, 2021). Trips within Madrid and Castilla-Leon were extracted from over 13 million phone records provided by a Spanish cellphone company. Both passive (GPS) and active (text messages, calls etc.) data were aggregated to construct daily movement matrices in each of the Communities, prior to the authors acquisition of the data. Given that trips were only recorded from one cellphone company, adjustment was made to estimate the number of total trips between each region. As a result, the entries of the mobility matrices are non-integer values.

Figs. 1(a) and 1(b) show the total daily movement between regions in Madrid, and Castilla-Leon, respectively. There is a sharp drop in the number of trips around March 14th 2020, which corresponds to a nation-wide lockdown. Lockdown restrictions began to ease around May 11th, where the number of trips slowly began to rise. Fig. 1, Fig. 1 show the number of cases of COVID-19 cases in both Communities. COVID-19 daily cases data were retrieved from the open data portal of Castilla-Leon (General Directorate of Information Systems, Quality and Pharmaceutical Provision, 2021) and from the Epidemiological Surveillance Network of Madrid (Anon, 2021). Notice that the movement drops as cases rise, because a lockdown was implemented in response to the increasing severity of the epidemic. In order to avoid this potential “reverse causality” problem, we will only use movement data in the first week of March. Our justification for this is that there is a time lag between when the virus spreads and the resulting COVID cases are confirmed. That is, the “first wave” of the epidemic was likely influenced mostly by the movement that occurred prior to the peak in cases, and less by the movement that occurred during it.

Fig. 1

Number of trips greater than 500 metres (a and b) and daily case counts (c and d) in the two Communities of Spain from March to June 2020.

Fig. 2 shows the spatial distribution of the COVID-19 case rates up until May 31, 2020. The cases per thousand people range from (approximately) in Madrid, and for Castilla-Leon. We can see that there is substantial variation in the case rates within each of these Communities. Note that the extreme values in these plots are mostly small regions, which makes sense since the variance of case rates is higher when population is small. In the north of Madrid, there is a cluster of municipalities that have very low case rates. In Castilla-Leon, case rates are highest near the southeast border, which is the border to Madrid.

Fig. 2

COVID-19 cases per thousand, up to May 31 2020 for two communities in Spain. Background map ©Stamen Design.

Fig. 3 shows the number of trips to, from, and within each Municipality of Madrid (there are 179 of these small regions), and Castilla-Leon (there are 245 health zones). Madrid and Castilla-Leon are considered separately throughout this paper. Although they are adjacent, data on movements between the two communities are not available. In Madrid, there is a lot of movement in and around Madrid City, and less movement in the more rural areas. Castilla-Leon shows a less predictable movement pattern, as there is not a single capital city that accounts for most of the movement. This movement data will be used to induce spatial correlation between regions, as described in Section 2.3.

Fig. 3

Number of trips (incoming, outgoing, and within) the 179 regions of Madrid, and 245 health zones of Castilla-Leon, for the period March 1 to March 7 2020.

Spatial autoregressive models

Besag, York, and Mollié (BYM) models (Besag et al., 1991) are widely used in spatial epidemiology and disease mapping due to their simplicity and computational efficiency. They assume the incidence of disease in region follows a Poisson distribution where is the number of infected cases in region , and is some form of expected count or offset, which could be the at-risk population, exposure time, etc. The log-relative risk, , is often modelled as where is the overall intercept, is the effect of spatial covariates, is the structured spatial random effect, and is the unstructured spatial random effect which allows for overdispersion in the response. In the spatial formulation of the BYM model, = 1 when regions and share a common border, and 0 otherwise. That is, region ’s structured spatial effect is only conditionally dependent on its neighbours, given all other regions. The distributions are known as the full conditionals, where is short hand for the set . We can see from (1) that is a weighted average of its neighbours, resulting in spatial smoothing. These full conditionals correspond to the joint distribution of the ’s being a Gaussian Markov random field (GMRF) (Rue and Held, 2005), with where is a matrix of weights such that for and , and is a variance parameter to be estimated. is a diagonal matrix such that . This definition ensures that the precision matrix, , is both symmetric and positive definite. In addition to the 0–1 weights based on regions being adjacent, other weighting schemes, such as inverse of Euclidean distance between regions, have been used. For a comparison of common weighting schemes, see Duncan et al. (2017). When we specify in this way, we refer to this as an Intrinsic Autoregressive (ICAR) model for . The joint density function has a computationally convenient form with which is sometimes referred to as the pairwise difference formula. Notice that this density is invariant to the addition of a constant to each , leaving the spatial random effects unidentifiable up to a constant. This is typically remedied by imposing the constraint (Duncan et al., 2017). We will now modify this BYM model to account for movement between regions, in addition to physical adjacency.

Movement augmented BYM model

In order to extend the BYM model to allow for spatial correlation based on movement data, a second ICAR term, , with dependence structure governed by the movement data is added to the model. We also retain an adjacency-determined spatial effect in order to infer the relative importance of mobility-based and adjacency-based spatial dependence in determining COVID-19 case counts. The resulting model is where and are the spatial random effects with priors based on the physical data and movement data respectively. The geographically-defined process has weights if regions and share a common border and are 0 otherwise, while the movement-defined process has weights representing the number of trips between regions and . Using mobility as edge weights in network models has shown to be effective in the context of infectious diseases (Schrödle et al., 2012, Volkova et al., 2010, Geilhufe et al., 2014). Schrödle et al. (2012) used mobility weights in an autoregressive term, which allowed the weights matrices to be asymmetric. However, given that our mobility data is being used in a Gaussian prior for a random effect, the precision matrices of and , and , must be symmetric. Therefore we require and . While the first equality will always be true, the mobility matrices are not perfectly symmetric, thus symmetry was induced by defining as the sum of the numbers of trips from to and from to . The GRMF does not account for the movement within a region, so the movement within a region was included in the model as a spatial covariate (fixed effect). That is, was computed as where is the number of trips per person within a region, and and are the mean and standard deviations of the trips per person in all other regions. This model was run on both the Madrid and Castilla-Leon data.

There are two main drawbacks with the formulations of BYM models presented thus far. Firstly, the interpretation of the parameters and depend on the average number of neighbours and the total number of trips for each region, and hence their magnitudes are not comparable (Sørbye and Rue, 2014). Secondly, , and are hard to estimate without very careful choices of hyperpriors (Leroux et al., 2000). We will now address these shortcoming via reparametrizations.

Reparametrizations and priors

In order to solve issues with comparability, interpretability, and estimation, we apply a reparameterization of our model that is inspired by Riebler et al. (2016) with where and . The priors for and are Note that

Here, is the combined variance of the spatial effects, and the ’s are mixing parameters, interpreted as the proportion of the combined spatial variance explained by each model component. Note that reduces the spatial component to purely overdispersion, reduces the spatial component of the model to an adjacency ICAR model for the spatial effects, and reduces the spatial component to a mobility ICAR model. Most importantly, if then this means that the mobility data better explains variation in COVID-19 case counts than the adjacency data. As long as the spatial weights matrix and the mobility weights matrix are linearly independent, then having both spatial and mobility terms in our model present no issues with identifiability (Rodrigues and Assunçao, 2012). Finally, and are scaling factors, such that the geometric means of and are both 1 for each , meaning that and are the log relative risk contributions from the movement data and physical data respectively (Sørbye and Rue, 2014). Scaling is absolutely necessary in order to conduct inference on the ’s. We compute the scaling factors as follows where is the generalized inverse of the precision matrix (Freni-Sterrantino et al., 2018). In order to scale the precision matrices of the spatial effects, the generalized inverse for sparse matrices from (Rue et al., 2017) was used. The diagonal elements, , of are referred to as the marginal variances of the structured spatial effects, i.e and .

As was the case with the ICAR model in (1), we can derive the full conditionals of the combined spatial effect, , for the model described in Section 2.3 These full conditionals can help provide some intuition as to the mechanism by which this model provides spatial smoothing. As , is simply the weighted sum of the other regions, where the weights are the proportion of region ’s total movement between each other region. If , the conditional mean of reduces to the arithmetic average of the spatial effects of its neighbours. If , then the conditional mean shrinks to 0 (remember that ). Given that is positive, the conditional mean is always shrunk towards 0, resulting in spatial smoothing. In practice, the conditional mean will be a weighted average of the estimates smoothed by the movement GMRF, the physical GMRF, and 0. It is important to note here that the and are relative measures due to the scaling factors. That is, doubling the total amount of movement has no effect on the conditional mean or variance of . This is in contrast to the combined spatial effects in the commonly used Leroux model (Leroux et al., 2000). Additionally, the variance of is lower when region has a lot of movement or many neighbours, relative to the other regions.

Inference, computation, and validation

Four chains each with 3000 iterations of No U-Turn Sampling were used for parameter estimation within Stan (Stan Development Team, 2021). The first 1500 iterations were used as a warm-up, the 1500 remaining iterations from each chain were thinned by a factor of 10, leaving 600 total posterior samples to perform inference. As mentioned in Section 2.2, we require . In practice, we use the soft constraint for computation purposes (as recommended by the Stan team Morris et al., 2019). To complete the model, priors for and were . To ensure the robustness of our results, we also ran BYM models using the adjacency data and the movement data separately. That is, for both Madrid and Castilla-Leon, we ran a model where we assumed , and a separate model where . The results of these four models are presented in Section 3.2.

Our code and posterior samples are posted at https://github.com/cghr-toronto/public/tree/master/covid/spain_public_code.

Results

Joint model

Table 1 shows posterior medians and credible intervals for the mixing parameters for the model with both movement and adjacency spatial effects. For both Madrid and Castilla-Leon, the proportion of spatial variation explained by is much higher than that of and . The posterior probability that was 0.997 for Madrid, and 0.998 for Castilla-Leon. However, does seem to account for a non-trivial amount of spatial variation in both Madrid and Castilla-Leon. This means that although movement data is likely more explanatory, adjacency data can help with explaining variation in COVID-19 cases. Additionally, there is a substantial amount of spatial variation explained by the unstructured spatial effect for Madrid. This is not the case for Castilla-Leon, as most of the mass of the posterior of is near 0. This makes sense given that Madrid has a large metropolitan centre surrounded by a mix of suburbs and rural areas, so there are probably spatial confounders that our model is missing. For a plot of the posterior densities of , see Appendix A.

Table 1

Posterior medians, and 95% credible intervals for in BYM models using movement and physical (adjacency) data in the same model.

	Parameter	Madrid	Castilla-Leon
		Est (95% CrI)	Est (95% CrI)
ρ	Movement	0.76 (0.54, 0.89)	0.88 (0.66, 0.98)
	Neighbour	0.13 (0.01, 0.39)	0.09 (0.01, 0.30)
	Independent	0.10 (0.02, 0.25)	0.02 (0.00, 0.09)
	μ	−5.36 (−5.51, −5.24)	−3.75 (−3.78, −3.73)
	β	0.12 ( 0.05, 0.20)	−0.01 (−0.04, 0.02)
	σ	0.65 ( 0.55, 0.78)	0.72 ( 0.63, 0.83)

Figs. 4(a) through 4(d) show the spatial distribution and , plotted using the same colour scale for comparability. We can see that ’s log-relative risks have a lot more spatial variation in both Communities. The log-relative risks for tend to have smooth spatial gradients, while tends to identify clusters of regions as high-risk areas. As seen in Eq. (2), the expectation of the combined spatial effects are a weighted average of these spatial effects, and 0 (notice that the numerator can be rewritten as where ). Fig. 4, Fig. 4 show the predicted cases per 1000 people per region, showing highly similar patterns to the observed values in Fig. 2.

Fig. 4

Log-relative risk contributions (a–d) from the movement effects () and spatial effects (). The predicted cases per thousand people are also presented (e–f).

The standard deviation was slightly larger for Castilla-Leon than it was for Madrid. Fig. B.2 shows the spatial distribution of the standard deviation of the cases per thousand people in both communities. Here, we can see that the standard deviation is pretty small in and around Madrid-city, because the movement to and from Madrid-city is causing a high-degree of spatial smoothing in the surrounding area. The effect of movement within regions, , is associated with larger case counts in Madrid, but this is not the case for Castilla-Leon. This small covariate effect could result in more variance being attributable to the random effects, potentially contributing to the larger in Castilla-Leon.

Fig. B.2

Standard deviations of predicted cases per thousand people.

Model validation - individual models

Table 2 shows posterior medians and credible intervals for the parameter from the movement and physical BYM models described in Section 2.5, fit separately to Madrid and Castilla-Leon (four models total). In both regions, the model where spatial smoothing is induced by population movement explains a higher proportion of the variation in the outcome, indicated by the posterior density of having more mass near 1. Additionally, the BYM model that used physical adjacency as a spatial smoother had a much wider credible interval for , indicating more model uncertainty. Both models show more uncertainty in the region of Madrid than for Castilla-Leon, likely due to the fact that Madrid is more heterogeneous in terms of population density and other factors. For full posterior densities of the parameter, see Fig. A.2.

Table 2

Posterior medians, and 95% credible intervals for in BYM models using movement and physical (adjacency) data in separate models.

	Parameter	Madrid	Castilla-Leon
		Est (95% CrI)	Est (95% CrI)
ρ	Movement	0.82 ( 0.66, 0.91)	0.95 ( 0.89, 0.98)
	Neighbour	0.56 ( 0.22, 0.83)	0.77 ( 0.58, 0.91)
μ	Movement	−5.34 (−5.48, −5.23)	−3.75 (−3.78, −3.73)
	Neighbour	−5.18 (−5.30, −5.09)	−3.74 (−3.78, −3.70)
β	Movement	0.12 ( 0.05, 0.18)	−0.02 (−0.05, 0.02)
	Neighbour	0.13 ( 0.01, 0.24)	−0.01 (−0.05, 0.04)
σ	Movement	0.63 ( 0.55, 0.76)	0.74 ( 0.65, 0.83)
	Neighbour	0.66 ( 0.56, 0.83)	0.58 ( 0.51, 0.66)

Fig. A.2

Posterior Density of the proportion of variance explained by spatial components when adjacency and movement data are used in separate models (model validation).

Discussion

In this paper, we have demonstrated that there is much value in using mobility data in combination with geographical proximity for defining correlation structures COVID-19 incidence data. We showed that even while using only one week of movement data, we were able to explain the spatial variation in COVID-19 counts better than using the classic BYM model. Additionally, we showed that the model can be re-parametrized so that the means by which smoothing occurs in these mobility models is intuitive.

A key limitation of this work is that the models presented in this paper do not serve as individual-level infectious disease models, as correlation is induced by a latent effect rather than direct dependence between the counts. However, this will be a natural extension of this work and would require the addition of many more parameters, including multiple mobility network components at various time points. This will ultimately pose a computational challenge as well.

An additional limitation of this work is that the availability and structure of mobility data will vary across data sources, and may only be available in higher income countries. Furthermore, there is selection bias in the movement data, as it only tracks those who actually have a cellphone, which may tend to be younger and more economically advantaged individuals. Given potential differences in quality of these data, its efficacy in spatial models may need to be assessed on a case by case basis.

Furthermore, the models presented in this paper may suffer from overfitting. A potential remedy for this would be to put a penalized complexity prior (Simpson et al., 2017) on the mixing parameters, which may improve inference by shrinking (and perhaps ) towards 0. An interesting area for future work would be to combine Dirichlet and penalized complexity priors to specify a joint prior for the mixing parameters as described in Fuglstad et al. (2020), which can be implemented using the makemyprior R package (Hem et al., 2021). This was deemed unnecessary for this work, as we were mainly interested in comparing to , and felt that our prior should not favour either one of these terms.

Despite these limitations, this work demonstrates the value of mobility data and provides the foundation for various extensions and future work. This data is only becoming more abundant as time passes, and methods that allow for efficient use of this data are essential to model the current epidemic, and any spatial epidemiological application where population movement is likely a predictor of disease.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.