A joint hierarchical model for the number of cases and deaths due to COVID-19 across the boroughs of Montreal.

As of July 2021, Montreal is the epicentre of the COVID-19 pandemic in Canada with highest number of deaths. We aim to investigate the spatial distribution of the number of cases and deaths due to COVID-19 across the boroughs of Montreal. To this end, we propose that the cumulative numbers of cases and deaths in the 33 boroughs of Montreal are modelled through a bivariate hierarchical Bayesian model using Poisson distributions. The Poisson means are decomposed in the log scale as the sums of fixed effects and latent effects. The areal median age, the educational level, and the number of beds in long-term care homes are included in the fixed effects. To explore the correlation between cases and deaths inside and across areas, three different bivariate models are considered for the latent effects, namely an independent one, a conditional autoregressive model, and one that allows for both spatially structured and unstructured sources of variability. As the inclusion of spatial effects change some of the fixed effects, we extend the Spatial+ approach to a Bayesian areal set up to investigate the presence of spatial confounding. We find that the model which includes independent latent effects across boroughs performs the best among the ones considered, there appears to be spatial confounding with the diploma and median age variables, and the correlation between the cases and deaths across and within boroughs is always negative.

Motivation

As of July 25th, 2021, the coronavirus disease 2019 (COVID-19) pandemic counts a total of 194,248,750 cases and 4,163,599 deaths worldwide (Dong et al., 2020). This disease is caused by an infection with the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The World Health Organisation declared COVID-19 to be a pandemic in March 2020 (Ghebreyesus, 2020) and since then, countries worldwide have instituted infection control measures (e.g., lockdowns, curfews, mask mandates) to control the spread of the disease. We are interested in studying the spatial distribution of the number of cases and deaths due to COVID-19 in the urban agglomeration of Montreal. As of July 25th, 2021, Canada cumulates 1,426,903 cases and 26,547 deaths due to COVID-19 (Government of Canada, 2021) with 9.2% of the cases and 17.9% of the deaths recorded in Montreal (Direction régionale de santé publique, 2021). The agglomeration covers approximately 500 and is the most populous administrative region in the province of Quebec, with over 1.9 million inhabitants which represents, approximately, 6% of the Canadian population (Statistics Canada, 2017a, Statistics Canada, 2017b). The agglomeration incorporates 19 boroughs of the City of Montreal and 15 related cities. From now on, we will refer to these areas as boroughs. These boroughs vary in demography, socio-economic status, and healthcare infrastructure.

This analysis focuses on data on the cumulative number of COVID-19 confirmed cases and deaths in Montreal since the beginning of the pandemic, as of July 25th, 2021. The cases and deaths distributions across the boroughs of Montreal are shown in Fig. 1. This dataset is made publicly available by the regional director of public health (Direction régionale de santé publique, 2021). The data consist of 34 boroughs. However, L’İle-Dorval did not record cumulative data. We exclude this borough from the analysis, yielding a total of 33 boroughs.

Fig. 1

Maps of the COVID-19 cases (left) and deaths (right) across the 33 boroughs of Montreal.

The Ministry of Health and Social Services defines a confirmed case or death due to COVID-19 as follows. Cases are confirmed either by a laboratory or by an epidemiological link. A laboratory confirms a case through the detection of nucleic acids of SARS-CoV-2. To confirm a case by an epidemiological link, the manifestation of clinical symptoms must be compatible with COVID-19 with a high risk of exposure, during the period of contagion, to a case that was confirmed by a laboratory and no other apparent cause. Deaths may also be confirmed by a laboratory or by an epidemiological link. A laboratory confirms a death due to COVID-19 through the manifestation of clinical symptoms compatible with COVID-19 before death and the detection of nucleic acids of SARS-CoV-2. A death is confirmed to be due to COVID-19 by an epidemiological link if there was a manifestation of clinical symptoms compatible with COVID-19 before death with a high risk of exposure, during the period of contagion, to a case that was confirmed by a laboratory and no other apparent cause (Ministère de la Santé et des Services sociaux, 2021a). The death counts of two boroughs are censored for privacy issues. Namely, Baie-d’Urfé and Montréal-Ouest recorded between 1 and 4 deaths. We also obtained data on the population size of each borough from the 2016 Canadian Census (Ville de Montréal, 2016).

In the literature, many studies identify the age of individuals as a risk factor of mortality due to COVID-19 (Zhou et al., 2020, Bonanad et al., 2020, Yanez et al., 2020). Some studies also highlight the relationship between socio-economic status and COVID-19 cases (Hawkins et al., 2020, Gangemi et al., 2020) and deaths (Hawkins et al., 2020, Bermudi et al., 2021). Moreover, apart from COVID-19, there exist important health disparities due to socio-economic inequalities in Montreal. These have been shown to lead to a difference in life expectancy of up to 11 years between boroughs (Le Blanc et al., 2011). On the other hand, healthcare infrastructure has been positively associated with reactive responses to COVID-19 (Sharma et al., 2021). In Quebec, there exist five classes of senior residences. Three of these classes of residences are called residential and long-term care centres (Centres d’hébergement et de soins de longue durée, CHSLDs), which may be public, private and subsidised by the government, or private and non-subsidised. Most senior residences are of the two remaining types: intermediary resources and private seniors’ residences. Although all of these residences can provide similar services, CHSLDs specifically serve seniors with an important loss of autonomy (Gouvernement du Québec, 2021). Hospitals in Quebec, under the direction of the Ministry of Health, were generally better prepared to respond to the COVID-19 pandemic than CHSLDs during the first months of the pandemic (Doucet, 2020). The vast majority of COVID-19 deaths in senior residences in Quebec are linked to CHSLDs. In fact, as of July 25th, 2021, 51% of all COVID-19 deaths in Quebec are linked exclusively to CHSLDs (Institut national de santé publique du Québec, 2021).

This literature review motivates our covariates’ selection in order to model the cases and deaths due to COVID-19 in the 33 boroughs of Montreal. To account for the age profile of each borough, we obtain data on the median age of the population by borough from the 2016 Canadian census (Ville de Montréal, 2016). As a proxy for the socio-economic status of each borough, we consider the percentage of the population between 25 and 64 years old with a university diploma. These data are also extracted from the 2016 Canadian census (Ville de Montréal, 2016). Finally, regarding the healthcare infrastructure in Montreal, we include the number of beds in all CHSLDs of each borough. These records are available in public databases (Ministère de la Santé et des Services sociaux, 2021b, Ministère de la Santé et des Services sociaux, 2021c).

In this paper, we analyse the joint distribution of COVID-19 cases and deaths across the 33 boroughs of Montreal through a Bayesian hierarchical joint model. We are interested in measuring the strength of the associations between the three available covariates and the cases and deaths due to COVID-19. Also, we aim at investigating how these two outcomes are correlated within and across boroughs. We examine whether the covariates impact the risk of being a case or the risk of dying from COVID-19 in a similar fashion. Additionally, we investigate if there is any structure left in the data after accounting for the available covariates. We explore different prior specifications for these local latent effects, ranging from independent to spatially structured ones. It is known that the inclusion of spatially structured latent effects can affect the estimation of the fixed effects (Reich et al., 2006, Khan and Calder, 2020, Dupont et al., 2022); this is known in the literature of Spatial Statistics as spatial confounding. We fit a Bayesian version of the Spatial+ approach recently proposed by Dupont et al. (2022) to investigate if there is spatial confounding in the fitted models.

This paper is organised as follows: Section 2 describes the proposed models, the method for spatial confounding adjustment and the inference procedure. In Section 3, the results of the analyses of the COVID-19 cases and deaths are presented. Finally, Section 4 concludes and discusses our findings.

Methods

Let be the recorded number of COVID-19 cases in borough , for , where is the number of non-overlapping boroughs in Montreal. Let be the number of deaths due to COVID-19 in borough . We model the deaths and cases through a joint hierarchical Bayesian model in order to accommodate the natural relationship one expects between cases and deaths of a disease. First, we define a marginal distribution for the cases, and then, conditional on the number of cases, we define a distribution for the number of recorded deaths. Sahu and Böhning (2021) modelled the cases and deaths due to COVID-19 in England in a similar fashion, using a 2-stage Bayesian hierarchical model. However, their interest lied in the temporal structure of the data and did not consider joint latent effects. Our approach differs because we do not conduct a spatio-temporal analysis, and we are interested in measuring the correlation between cases and deaths inside and across the boroughs of Montreal. To that end, we include a correlation parameter that helps borrow strength from the cases and deaths inside and across boroughs. More specifically, we assume, where stands for the Poisson distribution, and denote, respectively, the relative risk of the cases and the deaths in borough , whereas and are offsets. Note that if , then with probability one, since the number of deaths is necessarily smaller or equal to the number of cases. For the analysis performed here, all the counts are strictly positive. If zeros are observed, the codes to fit the proposed model and made available here need to be adapted to accommodate this case. The offset for the number of cases is computed based on the population size of each borough, while the offset for the number of deaths given the number of cases is computed based on the number of cases observed in borough . More specifically, and . To model the number of deaths conditionally on the number of cases, the Poisson distribution is used as an approximation to a Binomial distribution with size and probability . This approximation is reasonable because of high values of ’s and small ’s as estimated in an exploratory data analysis (Wakefield, 2013 section 7.3.2). In the next step, we model the log relative risks as follows: where and denote, respectively, the overall mean log risks of cases and deaths, is a vector of covariates used to model both the cases and deaths associated with the -dimensional vectors of coefficients, and , and and are latent random effects that accommodate whatever is left after accounting for the available covariates. Further, the inclusion of the latent random effects in both log risks’ decompositions allows for the accommodation of overdispersion which is observed in the exploratory data analysis.

We explore bivariate models for the latent random effects, . The models that we consider are special cases of the general formulation , where the vector of independent random effects, , , is independent of , the vector of spatially structured effects following independent CAR models (Besag, 1974), and where is the -dimensional identity matrix and Let the CAR distributed effects , , where , and is the matrix of weights that defines the neighbourhood structure. Commonly, we define if boroughs and share a border, denoted by , and , otherwise. This neighbourhood structure is used throughout this paper. Let and , we may write , which is not properly defined since is not positive definite (Banerjee et al., 2014).

The first special case of the model for that we consider is one with independent random effects across the boroughs, denoted the IID model: . Hence, denotes the intra-borough correlation between the latent effects for cases and deaths. When , the IID joint model results in independent random effects for the cases and deaths inside each borough. Further, if one assumes prior independence between the fixed effects in both models, then fitting the joint model provides the same results as those obtained by fitting the two models, for cases and deaths, separately. In other words, the number of cases and deaths do not borrow strength from each other within boroughs. If , this bivariate IID model allows for dependence between the deaths and cases within a particular borough.

The IID model, however, does not allow for spatial autocorrelation. Yet, one may expect the cases from neighbouring boroughs to be more correlated than cases from further apart boroughs, and similarly for the numbers of deaths. To adjust for this possible spatial autocorrelation, the second special case that we consider for is a multivariate CAR model (Lawson, 2020): . One may write the joint distribution for the long vector of latent effects, , as  (Jin et al., 2007). From this joint formulation, we see that the multivariate CAR model allows for dependence between the cases and deaths of a particular borough as well as between the counts of neighbouring boroughs.

The multivariate CAR model assumes a priori that the latent effects are necessarily distributed according to a spatial structure. To relax this assumption, we consider the multivariate extension of the BYM model (Besag et al., 1991) that allows for both unstructured and spatially structured sources of variability. This corresponds to the general formulation, . The long vector of latent effects is distributed as , for

Let . Although the marginal distribution of with respect to the latent effects is not available in closed form, it is possible to obtain the marginal moments using the properties of conditional expectations and the law of total covariance. Table 1 shows the marginal moments of the cases and deaths obtained from the models discussed above. The detailed computation of these marginal moments is available in section 1 of the Supplementary Material.

Table 1

Marginal moments of for each model, where .

Model	Marginal moments
IID	E(Yiℓ)=Eiℓexpβ0ℓ+x_iβ^ℓ+σℓ,v2/2≡μiℓ
	V(Yiℓ)=μiℓ1+μiℓexp(σℓ,v2)−1
	Cov(C_i,D_j)=μiCμiDexp(ρ_vσ_C,vσ_D,v)−1,ifi=j,0,otherwise.
CAR	E(Yiℓ)=Eiℓexpβ0ℓ+x_iβ^ℓ+σℓ,u2[Q^−]_ii/2≡μiℓ
	V(Yiℓ)=μiℓ1+μiℓexpσℓ,u2[Q^−]_ii−1
	Cov(C_i,D_j)=μiCμiDexpρ_uσ_C,uσ_D,u[Q^−]_ii−1,ifi=j,μiCμjDexpρ_uσ_C,uσ_D,u[Q^−]_ij−1,ifi∼j,0,otherwise.
BYM	E(Yiℓ)=Eiℓexpβ0ℓ+x_iβ^ℓ+σℓ,v2/2+σℓ,u2[Q^−]_ii/2≡μiℓ
	V(Yiℓ)=μiℓ1+μiℓexp(σℓ,v2)−1expσℓ,u2[Q^−]_ii−1
	Cov(C_i,D_j)=μiCμiDexpρ_uσ_C,uσ_D,u[Q^−]_ii+ρ_vσ_C,vσ_D,v−1,ifi=j,μiCμjDexpρ_uσ_C,uσ_D,u[Q^−]_ij−1,ifi∼j,0,otherwise,

From Table 1 it is clear that the IID model is only able to capture correlation between cases and deaths within a borough. The correlation will be negative if . On the other hand, the CAR and BYM models are able to accommodate correlations both within and among neighbouring boroughs. In the CAR model, if the correlation within and among neighbouring boroughs will be negative, and positive if . In the BYM model, there are two parameters capturing the correlation within a borough, and a negative correlation results if . On the other hand, for neighbouring boroughs the correlation is captured only by the parameter so that if a negative correlation results in this case.

Investigating spatial confounding: a Bayesian alternative to Spatial+

When modelling disease risks through a Poisson model, Clayton et al. (1993) noted that covariates effects may change in the presence of spatially structured latent effects. This issue is termed spatial confounding (Reich et al., 2006). Spatial confounding corresponds to the situation where spatially structured latent effects are correlated with the covariates, resulting in biased estimates of the fixed effects (Reich et al., 2006, Page et al., 2017). Reich et al. (2006) proposed a restricted spatial regression (RSR) model to overcome this issue, by removing the collinearity between the covariates and the spatial effects. However, Khan and Calder (2020) point out that RSR models are challenging to use. Additionally to leading to a loss of computational efficiency, these models cannot be naturally extended to other spatial models than the conditional autoregressive ones. Khan and Calder (2020) also note that their use may entail an increased type-S error for both correctly and incorrectly specified models. A type-S error occurs when a posterior credible interval for a particular coefficient does not include 0 even though the true value is 0. Recently, Dupont et al. (2022) proposed an approach to accommodate spatial confounding called Spatial+. The proposal is to “regress away” the spatial structure from the covariates and only include the residuals from such regression in the modelling of the variable of interest. Dupont et al. (2022) suggest to use Spatial+ as a tool to investigate for the presence of spatial confounding. If after fitting a model with and without spatial effects the fixed effects estimates are similar, this might be an indication that there is no spatial confounding. The Spatial+ approach of Dupont et al. (2022) is based on generalised additive models (GAM). As pointed out by Dupont et al. (2022), because of the relationship between splines and Gaussian Markov random fields (GMRF), Spatial+ can also be used when one captures the spatial structure through a GMRF. Next, we follow Schmidt (2022) and describe how to adapt the Spatial+ approach when the latent effect follows a CAR prior distribution.

Let be the vector for the th covariate. For each covariate , we model the covariates through Gaussian distributions with latent spatial effects that follow a CAR distribution a priori. More specifically, let where, a priori, . After assigning prior distributions to the parameters for the model of each we follow the Bayesian paradigm and obtain a sample from the resultant posterior distribution. Then, we define as the posterior mean of the residual . That is, , where and are, respectively, the point estimates (posterior means) of and . Finally, each vector of potentially spatially confounded covariate, is replaced by in the models for the cases and deaths that include spatial effects.

Inference procedure

In the fully Bayesian framework that we consider, regardless of the prior specification of the parameters, none of the Poisson models fitted to the data result in closed form posterior distributions. Hence, we approximate these posterior distributions using computational methods, specifically, Markov chain Monte Carlo (MCMC) methods. The MCMC procedures are run in R using the Nimble package, version 0.11.0 (de Valpine et al., 2017, de Valpine et al., 2021). The advantage of Nimble is that customary samplers are automatically available making the implementation of the MCMC simpler. Also, censored outcome values are naturally accounted for, as is necessary for the boroughs Baie-d’Urfé and Montréal-Ouest.

The CAR prior is not a proper distribution and is invariant to the addition of a constant (Rue and Held, 2005). This implies that there is an identifiability issue between the spatially structured effects and the global mean, the intercept. To ensure that the posterior distribution is proper, we impose a sum-to-zero constraint for each of the CAR distributed effects. For each of the MCMC iterations, we impose . Lawson (2020) discusses how to implement the multivariate CAR model in Nimble.

The analysis presented in this paper is conducted using R, version 4.0.5. All the Nimble codes and datasets used are available at https://github.com/vicmic13/Covid_Joint_CasesDeaths (Michal et al., 2022).

Results

For each of the 33 boroughs of Montreal, we have available the case and death counts due to COVID-19 that were recorded from the beginning of the pandemic, until July 25th, 2021. Let be the case count in borough , and , the death count. The spatial distributions of the cases and deaths are shown in Fig. 1 in Section 1. To model the cases and deaths, we consider as covariates the number of CHSLD beds, the percentage of the population with a university diploma, and the median age in each borough. Fig. 2 shows the maps of the variables as included in the models. The median age is scaled and the log of the number of beds is computed. Note that the log scale is used to assume a normal approximation of the number of beds. This is necessary, when fitting the Spatial+ models, as we assume that the covariates are normally distributed when adjusting for potential spatial confounding following Dupont et al. (2022). The resulting covariates from the spatial confounding adjustment are shown in Figure 1 in the Supplementary Material.

Fig. 2

Maps with the distribution of the three covariates: log beds, diploma and age, included in the fitted models.

We are interested in assessing the association between the selected covariates and the number of cases and deaths. We are also interested in measuring the correlation between the cases and deaths inside a particular borough or between neighbouring boroughs. To that end, we allow for a spatial autocorrelation between the cases and deaths, as described in Section 2. Hence, we fit six different Poisson models, as described in Section 2. Table 2 presents the models characteristics, namely the form of the latent effects, , and the covariates included. For all these models, the same priors are defined for the parameters involved in the fixed effects: . Additionally, the relevant standard deviations that appear in the latent effects’ prior distributions are given a half-Cauchy prior, , as suggested by Gelman (2006). Finally, when appropriate, the latent effects’ correlation parameters are assumed to follow a noninformative uniform prior distribution: .

Table 2

List of the models fitted to the number of cases and deaths due to COVID-19 across the boroughs of Montreal. The symbol ✓denotes which components were included in the respective model.

	Latent effects		Covariates included in the fixed effects
	Unstructured	Structured	Original	Adjusted for spatial confounding
	Γ_vv_i	Γ_uu_i	x_⋅k	r_⋅k
Simple			✓
IID	✓		✓
CAR		✓	✓
BYM	✓	✓	✓	–
CAR+		✓		✓
BYM+	✓	✓		✓

The models are fitted through the R package Nimble (de Valpine et al., 2021). The MCMC procedure consists of 2 chains of 800,000 iterations each, with a burn-in period of 400,000 iterations and a thinning factor of 160. An elliptical slice sampler is used for the regression coefficients, as these are normally distributed a priori. This sampler was chosen to efficiently sample from these posterior full conditionals. Metropolis–Hastings algorithms are used for the independent random effects in the IID, BYM and BYM+ models. For the spatial effects included in the CAR, BYM, CAR+, and BYM+ models, the default CAR_normal sampler from Nimble is used. The chains have mixed well, as assessed by the trace plots, effective sample sizes and the statistic proposed by Gelman and Rubin (1992). The reason for the need of these long chains seems to be related to the possible collinearity between the variable diploma and the models that include a spatial effect. The IID model and the ones that accommodate spatial confounding showed convergence before reaching the burn in of 400,000 iterations. For consistency, all models were fitted with the same number of iterations.

The posterior summaries for the fixed effects’ coefficients are presented in Fig. 3 for each model. The solid circles correspond to the estimated posterior means, while the segments are the limits of the 95% posterior credible intervals. The first row corresponds to the fixed effects in the log relative risk of cases and the second row, to the ones in the log relative risk of recorded deaths due to COVID-19. Clearly, the inclusion of random effects in each of the equations change the importance of each of the covariates. The log number of beds is negatively associated with the log relative risk of cases under the Simple model and as a random effect is included in the model, zero falls within the 95% posterior credible interval, suggesting that the number of beds is not associated with the cases. For the log relative risk of deaths, on the other hand, the association with the number of beds is strictly positive. This is expected as the majority of deaths in Quebec are connected to CHSLDs.

Fig. 3

Posterior summaries for the model coefficients across all the fitted models. Solid circles: posterior means; Vertical lines: 95% posterior credible intervals; Dashed line: indicates no association.

The diploma variable is negatively associated with the log relative risk of cases. The range of the 95% posterior credible interval is much wider under the CAR+ and BYM+ models than the ones obtained under the IID and CAR models. Although diploma results in a positive association with the log relative risk of death under the Simple model, once a random effect is included, zero is within the 95% posterior credible of the coefficient associated with diploma. This suggests that the average educational level of the borough is not associated with the risk of death.

For the median age of the borough, the Simple and IID models suggest a negative association with the log relative risk of cases and a positive association with the log relative risk of deaths. For the log relative risk of cases once a spatial structure is included in the model, zero falls within the 95% posterior credible interval. This is also true for the CAR+ and BYM+ models which account for spatial confounding.

Another goal of this analysis, other than to investigate the association of the log relative risks of cases and deaths with the available covariates, is to estimate the correlation between the cases and deaths within and across boroughs. We model the cases and deaths jointly in order to borrow strength from the different recordings both within each borough (e.g., IID model) and between boroughs (e.g., CAR model). Fig. 4 shows the posterior summaries (posterior mean and posterior 95% credible interval) for the correlation coefficients, and , respectively included in the latent effects’ IID, CAR, CAR+, BYM or BYM+ models. The BYM and BYM+ yield posterior 95% credible intervals that include 0 for and . On the other hand, the IID, CAR and CAR+ result in negative correlations between the latent effects for the cases and for the deaths (, and , respectively). This suggests that modelling the cases and deaths jointly is adequate. Note that for the BYM and BYM+ models, if the posterior 95% credible limits for and include zero, this does not necessarily imply that the marginal correlations will also include zero. From Table 1, it is clear that under models BYM and BYM+ the correlation between cases and deaths within and across boroughs also depends on the latent effects’ standard deviations as well as the spatial structure. Figure 2 of the Supplementary Material shows the posterior summaries for the intra-borough correlation between cases and deaths for each model. Although and have posterior credible intervals that include zero, the BYM model yields strictly negative intervals in 19 boroughs.

Fig. 4

Posterior summaries for the correlation coefficient of the latent effects. Solid circles: posterior means; Vertical lines: 95% posterior credible intervals. Dashed line: and .

Additionally to the posterior summaries for the intra-borough correlations between cases and deaths obtained for each model, which are available in Figure 2, section 2, of the Supplementary Material, the posterior means for the correlation between boroughs estimated under the CAR, BYM, CAR+, and BYM+ models are available in section 2 of the Supplementary Material, in Figure 3, Figure 4, Figure 5 and Figure 6, respectively. We note that all the estimated correlations are negative. Fig. 5 maps the posterior means of the correlations between cases and deaths inside each borough for each model that includes latent effects. The IID model does not allow for correlation across boroughs but does accommodate an intra-borough correlation (see Table 1). All the models estimate negative correlations between the cases and deaths on average a posteriori. For example, the correlation’s posterior mean in Senneville is for the CAR model and for the BYM model. In Côte-des-Neiges-Notre-Dame-de-Grâce, the posterior mean is for the CAR+ model and for the BYM. For the IID, CAR and CAR+ models, the correlations’ posterior 95% credible intervals (see Figure 2 in the Supplementary Material) are always negative. In Fig. 5, there does not seem to be a spatial structure in the negative posterior means of the intra-borough correlations between the cases and deaths due to COVID-19. Additionally, Figure 7 in the Supplementary Material maps the relative risks of cases and deaths as estimated by each model. The negative correlations between cases and deaths can be visualised in this Figure 7 in the Supplementary Material. For example, the southwest region of Montreal has low estimated risks of cases but high estimated risks of deaths for each model.

Fig. 5

Maps of the posterior means of the intra-borough correlation between cases and deaths for each model. These correlations were estimated based on the equations for the covariance under each model shown in Table 1.

Finally, we compare the models’ performances using the WAIC (Watanabe and Opper, 2010) in Table 3. In terms of WAIC, for which smaller values are preferred, the Simple model is the least adequate among the fitted ones, with a value of 3320 and 265 effective number of parameters. This suggests that there is a need for latent effects in the modelling of cases and deaths due to COVID-19 in Montreal. However, there does not seem to be a need for a spatial structure in these random effects as the CAR, BYM, CAR+, and BYM+ all perform worse than the IID model. The IID model yields a WAIC of 205 whereas the others yield a WAIC of 212, 225, 243 and 227, respectively. The CAR model with covariates adjusted for potential spatial confounding does not perform as well as the CAR model with the original covariates. However, there is no significant difference in terms of WAIC between the BYM and BYM+ model.

Table 3

Values of WAIC and effective number of parameters for each fitted model. The smallest value indicates the best model among fitted ones (in italics).

	Simple	IID	CAR	BYM	CAR+	BYM+
WAIC	3320.1	205.0	212.2	224.9	242.7	227.1
p_W	265.2	57.2	60.7	67.2	75.3	68.0

Discussion

In this paper, we analyse the case and death counts due to COVID-19 across the 33 boroughs of Montreal, as of July 25th, 2021. We fit six different models with conditional Poisson distributions for the cases and for the deaths in each borough. The areal means are decomposed in the log scale as the sum of an offset, fixed effects, and latent random effects. One model is fitted without latent effects. We then allow for correlation between the cases’ latent effects and the deaths’ latent effects within a borough. In one model, we assume that the latent effects are independent across the boroughs, whereas in the other four models, we accommodate a potential spatial autocorrelation between the latent effects of neighbouring boroughs. In two of these models, we set multivariate CAR and BYM priors for the latent effects. In the final two models, we first adjust the covariates for potential spatial confounding following Dupont et al. (2022) and then consider again multivariate CAR and BYM priors for the latent effects.

Interestingly, the posterior summaries for the fixed effects parameters are reversed between the cases and deaths components of the models. Regarding the log number of beds, the coefficient has posterior 95% credible intervals that include zero for the cases (through the IID, CAR, BYM, CAR+, and BYM+ models) but the posterior credible intervals are strictly positive for the deaths. We used the number of beds in CHSLDs as a proxy for the healthcare infrastructure in Montreal. This result suggests that this covariate has no link with the recorded number of cases but that, given the number of cases, an increased number of beds corresponds to a higher risk of death due to COVID-19. This result is consistent with the COVID-19 situation in Montreal and in the province of Quebec. In particular, CHSLDs were the epicentre of the COVID-19 crisis in Quebec (Hsu et al., 2020, Institut national de santé publique du Québec, 2021). The associations between the percentage of the population with a university diploma and the cases and deaths due to COVID-19 in Montreal are also different for each outcome considered. There does not seem to be a relationship between the percentage of people with diploma and the risk of death due to COVID-19, given the cases, for the IID, CAR, BYM, CAR+, and BYM+ models. On the other hand, each model results in negative posterior credible intervals for the effects of the diploma on the number of cases. The diploma variable is used as a proxy to the socio-economic situation of each borough in Montreal. This result seems to be aligned with other studies that also found that the level of education is negatively associated with the log risk of cases (see e.g., Hawkins et al., 2020). Finally, the risk of COVID-19 cases does not seem to be associated with age for the CAR, BYM, CAR+, and BYM+ models whereas older populations seem correlated with higher risks of death due to COVID-19. Once again, this result agrees with the current knowledge on COVID-19: all may contract COVID-19, but younger people are less likely to die from this disease (Williamson et al., 2020).

As the importance of the covariates changed when a latent spatially structured random effect was considered, we also investigated if there is spatial confounding by fitting a Bayesian version of the method proposed by Dupont et al. (2022). We find that the log number of beds does not seem to be spatially confounded whereas the diploma and median age do. As suggested by Dupont et al. (2022), we believe there is an interest in adjusting for spatial confounding even when there is no intuition of spatial confounding. Their method can be used as a tool to investigate the presence of such confounding. However, in terms of WAIC, the model with only IID latent effects seems to perform the best among the fitted ones. This suggests that there is no need to account for spatially structured latent effects. Note that the adjustment for spatial confounding is unnecessary in this IID model.

We believe that there is an interest to the data analysis conducted in this paper, as we are able to provide an estimate of the correlations between the cases and deaths due to COVID-19 within, and across, boroughs. Fitting separate models for cases and deaths, or even joint models that assume parameters to be independent a priori, such as in the Simple model, impose independence between cases and deaths. Our proposed approach allows the data to drive the inference procedure, and if there is correlation between the two outcomes, this can be estimated through the marginal correlations computed from the results shown in Table 1. Note that in terms of WAIC, we find that the Simple model performs worse than the five other models that include latent effects which allow for marginal correlations between cases and deaths both within and across boroughs. Finally, we find that the correlation between the cases and deaths due to COVID-19 within and between boroughs is always negative. This may be due to the aggregation of cases and deaths over time. We only have available aggregated data, which may yield different correlations than what would be observed at a certain time or through a temporal analysis. Looking at aggregated provincial time series of cases and deaths due to COVID-19, the number of cases and deaths were positively correlated at the beginning of the pandemic (Institut national de santé publique du Québec, 2021). We hypothesise that the estimated negative correlation between cases and deaths are related to different policies that were instituted throughout the pandemic and, moreover, due to the fact that people started getting massively vaccinated from May 30th, 2021, on. For instance, in the province of Quebec, the number of cases increased from 573 to 1683 between March 15th, 2021, and April 15th, 2021, whereas the vaccination coverage for the population aged 80 and older increased from 57.9% to 91.1% in the same period (Institut national de santé publique du Québec, 2021). It can be observed in the aggregated time series that, in this period, the daily number of deaths somewhat stabilised at low numbers.

In this paper, we used the median age of the boroughs as a proxy to the age population profile. If data for different age groups are available, one could consider models for case and death counts for each age stratum, allowing for different effects of the available variables and correlation parameters across the different age groups. In this framework, one would need to carefully consider the inclusion of the number of CHSLD beds in the various strata.

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.