Markov switching zero-inflated space-time multinomial models for comparing multiple infectious diseases

Comparing the transmission dynamics of several co-circulating infectious diseases is an important problem in epidemiology. In this paper, we are specifically concerned with comparing the transmission dynamics of the vector-borne diseases Zika, chikungunya and dengue. Since 2015 several countries in Latin America have experienced triple epidemics of these diseases, placing an immense burden on the local populations (Rodriguez-Morales et al., 2016; Bisanzio et al., 2018). Despite being transmitted by the same vector, the Aedes aegypti, several studies have found spatio-temporal differences in clusters of the three diseases (Freitas et al., 2019; Kazazian et al., 2020). Investigating how these differences may be related to covariates could help more efficiently distribute medical resources, as health interventions can vary between the diseases (Freitas et al., 2019). Schmidt et al. (2022) attempted this using a multinomial model, but they only investigated how the relative distribution of the diseases varied spatially. There is also an interest in seeing how the transmission of the three diseases may differ by spatio-temporal factors. For instance, laboratory studies have shown that the transmission of Zika may be more sensitive to temperature compared to dengue, which is important as it implies Zika will not be able to spread as effectively into North America and the highlands (Tesla et al., 2018).

As a motivating example, we analyze bi-weekly reported case counts of Zika, chikungunya, and dengue in the 160 neighborhoods of Rio de Janeiro, Brazil, between 2015-2016, during the first triple epidemic there. Since we are mainly interested in comparing the diseases, we condition on the total number of cases in each area and bi-week and use a multinomial distribution to model their relative allocations (Dabney and Wakefield, 2005; Dreassi, 2007; Schmidt et al., 2022). Several authors have proposed modeling multinomial data across time, space, or both by adding correlated random effects to multinomial link functions, such as the log relative odds (Pereira et al., 2018; Sosa et al., 2023). Proposed distributions for the random effects include conditional autoregressive distributions (Schmidt et al., 2022), multivariate logit-beta distributions (Bradley et al., 2019) and multivariate normal random walks (Cargnoni et al., 1997). Autoregressive approaches (Loredo-Osti and Sutradhar, 2012) have also been popular, where the current multinomial probabilities in an area are assumed to depend on past observations in the area (Fokianos and Kedem, 2003) and neighboring areas (Tepe and Guldmann, 2020). Since new infectious disease cases are transmitted from previous cases, conditioning on past observations is often a natural way to model infectious disease counts across time and space (Bauer and Wakefield, 2018). Therefore, we primarily take an autoregressive-driven approach in this paper. We assume that if conditions are more favorable for one disease over the others, then the share of that disease will grow over time. We represent the relative favorability of the diseases through a series of autoregressive parameters that are linked to covariates. This allows us to investigate whether a covariate is associated with conditions that are more favorable toward the transmission of one of the diseases. As we will show, our approach is epidemiologically motivated and can be derived from a series of Reed-Frost chain binomial susceptible-infectious-recovered (SIR) models (Abbey, 1952; Vynnycky, 2010; Bauer and Wakefield, 2018).

A potential issue with applying a multinomial model to multivariate spatio-temporal infectious disease counts is that many infectious diseases, especially vector-borne diseases, can be absent for long periods of time in an area (Bartlett, 1957; Coutinho et al., 2006; Adams and Boots, 2010). This behavior is illustrated in Figure 1 which shows bi-weekly reported case counts of dengue, Zika and chikungunya in the western Rio neighborhood of Praça Seca, between 2015-2016. While dengue cases are reported consistently throughout the two years, very few Zika cases and no chikungunya cases are reported in the first year. These patterns are typical of most neighborhoods. If a disease is absent in an area, then there should be no chance of cases being reported. However, a multinomial model will always assign a positive probability to cases being reported, which may lead to many more zeroes in the counts than can be realistically produced by the fitted model (Koslovsky, 2023). Also, we do not want to compare the transmission dynamics of a present disease and an absent one as it would likely overestimate the favorability of the present disease.

To account for long periods of disease absence, we will place our proposed autoregressive multinomial model within a zero-inflated multinomial modeling framework (Zeng et al., 2022; Koslovsky, 2023). Zero-inflated multinomial models allow categories, in our case diseases, to be absent from observations. If a category is absent, its multinomial probability is fixed at 0. The presence of a category is usually treated as random and generated through an independent Bernoulli process (Koslovsky, 2023). While a zero-inflated multinomial model would allow diseases to be absent in areas, existing approaches have made two independence assumptions that are inappropriate for spatio-temporal infectious disease counts. Firstly, they have assumed that the multivariate counts are independent conditional on the present categories of each observation. Infectious disease counts usually exhibit a high amount of correlation across space and time due to disease transmission (Bauer and Wakefield, 2018). Secondly, past approaches assumed that if a category is present in one observation it does not affect whether that category or other categories are present in other observations (independence of categorical presence). If an infectious disease is present in an area, then it is more likely to be present in the next time period, as disease presence and absence are usually persistent (Coutinho et al., 2006). Also, the presence of one disease may affect whether other diseases are present due to potential disease interactions (Paul et al., 2008). Finally, due to between-area disease spread (Grenfell et al., 2001), the status of the disease in an area will affect its presence in neighboring areas.

To address the above issues, we first choose a baseline disease that is well-established in the region and assume it is always present. This avoids all diseases being absent in an area at once which would be impossible if the total number of cases were not 0. We assume all other diseases switch back and forth between periods of presence and absence in each area through a series of Markov chains. The Markov chains account for long periods of disease absence, disease interactions and covariates, including the prevalence of the diseases in neighboring areas. For the diseases present in an area together, we assume that the relative distribution of their cases follows our proposed space-time autoregressive multinomial model. The absent diseases always report 0 cases. The rest of this paper is organized as follows. Section 2 introduces our proposed autoregressive multinomial model for comparing the transmission dynamics of multiple infectious diseases in a space-time setting. Section 2.1 extends the model to incorporate zero inflation and deal with long periods of disease absence. Section 3 details our Bayesian inferential procedure which utilizes the forward filtering backward sampling (FFBS) algorithm of Chib (1996) for efficient inference. In Section 4, we apply our model to cases of dengue, Zika and chikungunya in Rio de Janeiro during the first triple epidemic there. We close with a general discussion in Section 5.

Assume we have reported case counts for $K$ different infectious diseases across $i=1,\dots,N$ areas and $t=1,\dots,T$ time periods. Let $\bm{y}_{it}=(y_{1it},\dots,y_{Kit})^{T}$ represent the vector of case counts for all diseases in area $i$ during the time interval $(t-1,t]$, where $y_{kit}$ represents the reported case counts for disease $k$. We will first present our model without zero inflation and then extend to the zero-inflated case in Section 2.1. Since we are mainly interested in comparing the diseases, we condition on the total number of cases and model $\bm{y}_{it}$ using a multinomial distribution,

$\displaystyle\bm{y}_{it}|\text{total}_{it},\bm{y}_{t-1}\sim\text{Multinom}(\bm{\pi}_{it},\text{total}_{it}),$

(1)

where $\text{total}_{it}=\sum_{k=1}^{K}y_{kit}$ and $\bm{y}_{t-1}=(\bm{y}_{1(t-1)},\dots,\bm{y}_{N(t-1)})^{T}$ is the vector of all case counts reported in the previous time, giving the model an autoregressive structure. In (1), $\bm{\pi}_{it}=(\pi_{1it},\dots,\pi_{Kit})^{T}$ represents the expected relative distribution of the disease counts at time $t$ in area $i$ given the previous observed counts $\bm{y}_{t-1}$. If conditions in $(t-1,t]$ favor some diseases over others then the share of the favored diseases should grow relative to those less favored. Therefore, we model the relative odds as,

$\displaystyle\frac{\pi_{kit}}{\pi_{1it}}$

$\displaystyle=\lambda_{kit}^{*}=\lambda_{kit}\frac{\mathopen{}\mathclose{{}\left(y_{ki(t-1)}+1}\right)^{\zeta_{k}}}{\mathopen{}\mathclose{{}\left(y_{1i(t-1)}+1}\right)^{\zeta_{1}}},$

(2)

for $k=2,\dots,K$, where we assume disease 1 is the baseline disease and we add 1 to the top and bottom to avoid dividing by zero, adding such constants is common in autoregressive count models (Liboschik et al., 2017; Fritz and Kauermann, 2022). The parameters $\zeta_{k}$ and $\zeta_{1}$, assumed to be between 0 and 1, are meant to account for nonhomogeneous mixing (Wakefield et al., 2019). That is, as cases increase there tends to be a dampening effect on transmission as individuals avoid infection and governments implement interventions. In (2), $\lambda_{kit}>0$, which we link to covariates below, is a parameter that represents the favorability of conditions for the transmission of disease $k$ relative to the baseline disease during $(t-1,t]$. For instance, if $\lambda_{kit}>1$ ($\lambda_{kit}<1$) then we would expect the share of disease $k$ to grow (shrink) relative to the baseline disease in $(t-1,t]$, after adjusting for nonhomogeneous mixing with $\zeta_{k}$ and $\zeta_{1}$. If $\zeta_{k}\approx\zeta_{1}<<1$, like for our motivating example in Section 4, then this means that if $\lambda_{kit}>1$ the share of disease $k$ will grow initially from an equal share. However, as disease $k$ becomes more dominant the share will be pulled back towards equality as individuals avoid infection from the dominant disease.

We link $\lambda_{kit}$ to covariates and random effects using a log-linear model,

$\displaystyle\log(\lambda_{kit})=\alpha_{0ki}+\bm{x}_{kit}^{T}\bm{\alpha}_{k}+\phi_{kit},$

(3)

where $\alpha_{0ki}\sim N(\alpha_{0k},\sigma_{k}^{2})$ is a normal random intercept meant to account for between area differences and $\bm{x}_{kit}$ is a vector of covariates that might affect the favorability of disease $k$ relative to the baseline disease. We model the random effects $\phi_{kit}$ using a multivariate normal distribution (Xia et al., 2013; Zeng et al., 2022),

$\displaystyle\bm{\phi}_{it}$

$\displaystyle=(\phi_{2it},\dots,\phi_{Kit})^{T}\sim\text{MVN}_{K-1}(\bm{0},\bm{\Sigma}),$

(4)

where $\bm{\Sigma}$ is a $K-1$ by $K-1$ variance-covariance matrix. These random effects are mainly meant to account for overdispersion, which is a very common issue when modeling multinomial counts (Xia et al., 2013). They also allow for a more flexible correlation structure between the individual disease counts, see the supplementary materials (SM) Section 2.

There is an epidemiological interpretation of the parameters in Equations (1)-(4). As we show in SM Section 1, if all diseases follow a Reed-Frost chain binomial SIR model (Abbey, 1952; Vynnycky, 2010; Bauer and Wakefield, 2018) with roughly the same serial interval, the time between successive generations of cases, then we can derive Equations (1)-(4). Under this derivation, $\lambda_{kit}$ in (2) represents the ratio of the effective reproduction number of disease $k$ and the baseline disease. The effective reproduction number is an important measure of transmission in epidemiology that gives the average number of new infections produced by a single infectious individual in the current population before they recover (Vynnycky, 2010). The covariate effect $\alpha_{kj}$ then represents the difference between the effect of covariate $x_{kitj}$ on the effective reproduction number of disease k and the baseline disease, see SM Section 1 for more details.

While the Reed-Frost model is commonly used (Wakefield et al., 2019), it makes assumptions that are not appropriate for many diseases (Abbey, 1952), and the model can be sensitive to underreporting and reporting delay when it comes to estimating the effective reproduction number (Bracher and Held, 2021; Quick et al., 2021). Regardless, the Reed-Frost derivation does reveal an important source of potential confounding due to differences in the susceptible populations, i.e., the population of individuals who are not immune to the disease. From the Reed-Frost derivation, see SM Equation (8), ideally, we would add $\log(\delta_{ki(t-1)}/\delta_{1i(t-1)})$ to $\log(\lambda_{kit})$ in Equation (3), where $\delta_{ki(t-1)}$ is the size of the susceptible population for disease $k$ in area $i$ at time $t-1$. This makes sense intuitively, if a disease has a larger susceptible population compared to another then transmission for that disease will naturally be favored. The effect of any covariates correlated with the ratio of the susceptible populations could then be confounded. We conduct a sensitivity analysis to attempt to control for this for our motivating example in SM Section 8.

Finally, we note that the number of cases in neighboring areas may affect the favorability of the diseases as individuals will mix with those in neighboring areas (Grenfell et al., 2001). Therefore, we allow functions of disease cases in neighboring areas, such as the prevalence across neighboring areas, to be added to the covariate vector $\bm{x}_{kit}$. Also, we allow previous counts of other diseases, e.g., $\log(y_{ji(t-1)}+1)$ for $j\neq k$, to be added to $\bm{x}_{kit}$ to account for potential disease interactions (Freitas et al., 2019), see Section 4 for an example. We will call the model defined by (1)-(4) the autoregressive multinomial (ARMN) model.

Here we assume the number of diseases $K$ is small enough for matrix multiplication with a $2^{K}$ by $2^{K}$ matrix to be computationally feasible. This should be the case in most epidemiological applications, for instance, in our motivating example we compare $K=3$ diseases. Our inferential strategy revolves around reparametrizing the MS-ZIARMN model as a Markov switching model (Frühwirth-Schnatter, 2006), which allows for very efficient Bayesian inference using the forward filtering backward sampling (FFBS) algorithm (Chib, 1996). To simplify the explanations we will assume, without loss of generality, that $K=3$.

In the case of $K=3$, we have $\bm{S}_{it}=(S_{2it},S_{3it})^{T}$ as $S_{1it}$ is always equal to 1. Let $S_{it}^{*}$ be an indicator for the possible values of $\bm{S}_{it}$, so that $S_{it}^{*}=1$ if $\bm{S}_{it}=(1,1)^{T}$, $S_{it}^{*}=2$ if $\bm{S}_{it}=(0,1)^{T}$, $S_{it}^{*}=3$ if $\bm{S}_{it}=(1,0)^{T}$ and $S_{it}^{*}=4$ if $\bm{S}_{it}=(0,0)^{T}$. Then, from Equation (6),

$\displaystyle\bm{y}_{it}|\bm{y}_{t-1},S_{it}^{*},\text{total}_{it}\sim\begin{cases}\text{Multinom}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\frac{1}{1+\lambda_{2it}^{*}+\lambda_{3it}^{*}},\frac{\lambda_{2it}^{*}}{1+\lambda_{2it}^{*}+\lambda_{3it}^{*}},\frac{\lambda_{3it}^{*}}{1+\lambda_{2it}^{*}+\lambda_{3it}^{*}}}\right)^{T},\text{total}_{it}}\right),&\text{ if $S_{it}^{*}=1$}\\[5.0pt] \text{Multinom}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\frac{1}{1+\lambda_{3it}^{*}},0,\frac{\lambda_{3it}^{*}}{1+\lambda_{3it}^{*}}}\right)^{T},\text{total}_{it}}\right),&\text{ if $S_{it}^{*}=2$}\\[5.0pt] \text{Multinom}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\frac{1}{1+\lambda_{2it}^{*}},\frac{\lambda_{2it}^{*}}{1+\lambda_{2it}^{*}},0}\right)^{T},\text{total}_{it}}\right),&\text{ if $S_{it}^{*}=3$}\\[5.0pt] \text{Multinom}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(1,0,0}\right)^{T},\text{total}_{it}}\right),&\text{ if $S_{it}^{*}=4$},\end{cases}$

(8)

meaning, for $K=3$, the MS-ZIARMN model is a mixture of 4 different multinomial distributions where the multinomial probabilities are fixed at 0 for the diseases absent.

Note that, from Equation (7), as $\bm{S}_{it}$ only depends on $\bm{S}_{i(t-1)}$, $\bm{y}_{t-1}$ and covariates, $S_{it}^{*}$ follows a first order nonhomogeneous Markov chain. Let $\Gamma(S_{it}^{*}|\bm{y}_{t-1})$ be the transition matrix of $S_{it}^{*}$, where $\Gamma(S_{it}^{*}|\bm{y}_{t-1})_{jk}=P(S_{it}^{*}=k|S_{i(t-1)}^{*}=j,\bm{y}_{t-1})$ for $j,k=1,2,3,4$, $i=1,\dots,N$ and $t=2,\dots,T$. As $S_{it}^{*}$ is an indicator for $\bm{S}_{it}$, the transition matrix can be derived through

$\displaystyle P(\bm{S}_{it}=\bm{s}_{it}|\bm{S}_{i(t-1)}=\bm{s}_{i(t-1)},\bm{y}_{t-1})=\prod_{k=2}^{3}P(S_{kit}=s_{kit}|S_{2i(t-1)}=s_{2i(t-1)},S_{3i(t-1)}=s_{3i(t-1)},\bm{y}_{t-1}),$

for $\bm{s}_{i(t-1)},\bm{s}_{it}=(1,1),(0,1),(1,0),(0,0)$, using Equation (7). Finally, as $S_{i2}^{*}$ depends on $S_{i1}^{*}$ we require an initial distribution for the Markov chain, that is, $P(S_{i1}^{*}=j)$ for $j=1,2,3,4$ and $i=1,\dots,N$. The modeler can set the initial distribution based on how likely they believe the diseases to be present at the beginning of the study period.

Given $p(\bm{y}_{it}|\bm{y}_{t-1},S_{it}^{*},\text{total}_{it})$ and $\Gamma(S_{it}^{*}|\bm{y}_{t-1})$ for $i=1,\dots N$ and $t=2,\dots T$, and $p(S_{i1}^{*})$ for $i=1,\dots,N$, we completely define a Markov switching model (Frühwirth-Schnatter, 2006). A Markov switching model assumes a time series can be described by several submodels, often called states or regimes, where switching between submodels is governed by a first-order Markov chain (Hamilton, 1989). In our case, the submodels are given by the 4 multinomial distributions in (8) and we switch between these submodels through the transition matrix $\Gamma(S_{it}^{*}|\bm{y}_{t-1})$. Typically, $S_{it}^{*}$ is called the state indicator.

Let $\bm{S}^{*}=(\bm{S}_{1}^{*},\dots,\bm{S}_{T}^{*})^{T}$ be the vector of all state indicators, where $\bm{S}_{t}^{*}=(S_{1t}^{*},\dots,S_{Nt}^{*})^{T}$, $\bm{y}=(\bm{y}_{1},\dots,\bm{y}_{T})^{T}$ be the vector of all observations, and $\bm{\beta}$ the vector of all parameters in the multinomial part of the model (see (2) and (4)); further, let $\bm{\theta}$ be the vector of all parameters in the Markov chain part of the model (see (7)), and, finally, let $\bm{v}=(\bm{\beta},\bm{\theta})^{T}$ be the vector of all model parameters. The likelihood of $\bm{v}$ given $\bm{S}^{*}$ and $\bm{y}$ is given by,

$\displaystyle p(\bm{S}^{*},\bm{y}|\bm{v})=\prod_{i=1}^{N}\prod_{t=2}^{T}p(\bm{y}_{it}|\bm{y}_{t-1},S_{it}^{*},\text{total}_{it},\bm{\beta})\prod_{i=1}^{N}p(S_{i1}^{*})\prod_{t=2}^{T}p(S_{it}^{*}|S_{i(t-1)}^{*},\bm{y}_{t-1},\bm{\theta}).$

(9)

Note, from Equation (8), $S_{it}^{*}$ is only fully known when $y_{1it},y_{2it},y_{3it}>0$, in which case we must have $S_{it}^{*}=1$. It is possible, assuming $K$ is not too large, to use the forward filter (Hamilton, 1989) to completely marginalize $\bm{S}^{*}$ from (9) and calculate $p(\bm{y}|\bm{v})$, see SM Section 3. However, we want to make inferences about $\bm{S}^{*}$ to investigate when the model believes the diseases were present or absent for model checking, see SM Section 5. Therefore, $\bm{S}^{*}$ is estimated along with $\bm{v}$ by sampling both from their joint posterior distribution $p(\bm{v},\bm{S}^{*}|\bm{y})\propto p(\bm{S}^{*},\bm{y}|\bm{v})p(\bm{v})$, where $p(\bm{v})$ is the prior distribution of $\bm{v}$. We specified wide independent priors for most lower-level elements of $\bm{v}$. We specified a conjugate inverse-Wishart prior for $\bm{\Sigma}$ with $K$ degrees of freedom and identity matrix, as it implies the correlations in $\bm{\Sigma}$ have marginal uniform prior distributions (Gelman et al., 2013).

As the joint posterior is not available in closed form, we resorted to Markov chain Monte Carlo (MCMC) methods, in particular, we used a hybrid Gibbs sampling algorithm with some steps of the Metropolis-Hastings algorithm to sample from it. The full details of the Gibbs sampler are given in SM Section 3. We sampled all of $\bm{S}^{*}$ jointly from $p(\bm{S}^{*}|\bm{v},\bm{y})$ using the FFBS algorithm for Markov switching models (Chib, 1996). We found our MCMC algorithm mixed much faster when jointly sampling $\bm{S}^{*}$ compared to sampling each presence indicator individually. This is typically the case for Markov switching models (Frühwirth-Schnatter, 2006) and was one of our main motivations for reparametrizing the MS-ZIARMN model as a Markov switching model.

Our MCMC algorithm was implemented using the R package Nimble (de Valpine et al., 2017). We implemented the FFBS algorithm using Nimble’s custom sampler feature. All other MCMC samplers mentioned in SM Section 3 are built into Nimble. All code and data are available on GitHub https://github.com/Dirk-Douwes-Schultz/MS_ZIARMN_code.

Between 2015 and 2016 Rio de Janeiro experienced its first triple epidemic of dengue, Zika and chikungunya. Using space-time cluster analysis, Freitas et al. (2019) found that Zika and dengue clusters were much more likely in the western region of Rio (there were no chikungunya clusters found in the west), while other regions contained many clusters of each disease though they often differed in time. However, they did not investigate how these differences in clustering might depend on covariates or other factors. Schmidt et al. (2022), using a spatial multinomial regression model, found that cases of chikungunya were more likely in urban areas and Zika cases were less likely in population-dense areas, compared to dengue cases. A spatial analysis cannot investigate how space-time factors might have affected differences in the transmission intensity of the diseases, which may be important in this example, see below. Therefore, we apply our MS-ZIARMN model to cases of the three diseases, $K=3$, across the 160 neighborhoods of Rio de Janeiro, $N=160$, for each bi-week between 2015-2016, $T=52$. Our model seems appropriate for this example as there appear to be long periods of disease absence in both the time series of chikungunya and Zika cases, see Figure 1 for example (overall 65% of Zika cases are equal to 0 and 75% of chikungunya cases are equal to 0).

We have proposed a zero-inflated spatio-temporal multinomial model for comparing the transmission dynamics of multiple co-circulating infectious diseases. Our approach can account for long periods of disease absence while investigating how certain factors are related to differences in the transmission intensity of the diseases. Past zero-inflated multinomial models (Tang and Chen, 2019; Zeng et al., 2022; Koslovsky, 2023) have made independence assumptions that are inappropriate for spatio-temporal infectious disease counts, such as assuming categorical presence was independent between observations. We accounted for correlations across space, time and disease, in both the multinomial and zero-inflated components, through a combination of autoregression, regressing on past observations, and assuming the presence of the diseases followed a series of coupled Markov chains. This approach allowed for efficient and computationally feasible Bayesian inference using the FFBS algorithm.

We assume that the cases of the diseases present in an area jointly follow a multinomial model, conditioning on their total. It is more standard in the disease mapping literature to jointly model multiple disease counts with some form of a multivariate Poisson distribution (Jack et al., 2019). There are advantages to both approaches and we give a more detailed comparison of the two (without model fitting) in SM Section 1.2. If one is only interested in comparing the transmission of the diseases, like here, then a multinomial model eliminates many nuisance parameters and any shared space-time factors. A multivariate Poisson approach could estimate the overall effects of a covariate on transmission, not just the relative effects, and does not need to assume one of the diseases is always present. To the best of our knowledge, only Pavani and Moraga (2022) and Rotejanaprasert et al. (2021) have considered zero-inflated spatio-temporal multivariate Poisson models. Both relied on random effects to capture correlations in the data. However, these proposals have some important limitations, including not considering any space-time or space-time-disease interactions and not allowing the probability of disease presence to change over time. These modeling assumptions are inappropriate for our motivating example, and likely others, see Figure 4 for instance. An approach based on random effects could become intractable as one would need to consider random effects correlated across three dimensions, space, time and disease, in both the count and zero-inflated components of the model. Our approach seems much more computationally feasible and could be easily extended to the multivariate Poisson case, see SM Section 1.2. We will focus on these extensions and comparisons with the multinomial model in future work.

In our application, we found that Zika generally had more intense transmission compared to dengue and chikungunya, but was also not able to transmit as well at a lower temperature. This was not a major factor in tropical Rio, see Figure 3, but could be important in North America and Europe. In contrast, chikungunya transmission was relatively higher at lower temperatures meaning it might fare better in colder regions compared to Zika and dengue. Laboratory studies have come to similar conclusions (Tesla et al., 2018; Mercier et al., 2022). Alternative models that did not account for zero inflation or did not model correlations in disease presence (Zeng et al., 2022) did not fit as well and produced estimates inconsistent with our knowledge about the epidemiology of the diseases.

There are also some important limitations to our work. Being an observational study, unmeasured confounding could be an important issue, especially concerning differences in the susceptible populations which are unobservable. We tried using differences in cumulative incidence as a proxy for differences in the susceptible populations, in a sensitivity analysis in SM Section 8, and our results did not change. However, differences in cumulative incidence could be an imperfect proxy due to, for example, changes in reporting rates across space and time. Another limitation is that climate covariates, such as temperature, can have lagged and non-linear effects on the transmission of arboviral diseases (Lowe et al., 2018). As our multinomial model estimates differences in effects on transmission, our approach might be somewhat robust to non-linearity if the non-linear patterns are similar between diseases. It would be challenging to incorporate non-linear and lagged effects into our already complex modeling scheme. Shrinkage methods could help deal with the many parameters that would need to be added to the model by eliminating unimportant effects (Wang et al., 2023).

This work is part of the PhD thesis of D. Douwes-Schultz under the supervision of A. M. Schmidt in the Graduate Program of Biostatistics at McGill University, Canada. Douwes-Schultz is grateful for financial support from IVADO and the Canada First Research Excellence Fund/Apogée (PhD Excellence Scholarship 2021-9070375349). Schmidt is grateful for financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada (Discovery Grant RGPIN-2017-04999). This research was enabled in part by support provided by Calcul Québec (www.calculquebec.ca) and Compute Canada (www.computecanada.ca).