Spatio-temporal modeling of co-dynamics of smallpox, measles, and pertussis in pre-healthcare Finland

Data

During the study period, 1820–1850, the parishes in Finland kept track, among others, of baptisms, burials and causes of deaths, according to common and long held principles (Pitkänen, 1977). Even though the death diagnostics were based on symptoms, some infections can be considered to be diagnosed rather accurately due to their characteristic features. These diseases include pertussis (whooping cough), measles, and smallpox, which we consider. These infections were the main reason for child mortality, and, overall, they were responsible for approximately 5, 3, and 3 percent of total deaths, respectively, according to our data. Based also on the available records in our data, the median ages of deaths in complete years were 0 (sd = 3.6) for pertussis, 2 (sd = 3.9) for measles, and 2 (sd = 7.5) years for smallpox. Smallpox vaccinations were started in Finland in 1802 and were slowly progressing during the study period (Briga, Ketola & Lummaa, 2022). However, general healthcare was almost non-existent as in 1820 there were only 373 hospital sickbeds for 1.2 million inhabitants (Saarivirta, Consoli & Dhondt, 2012).

Our data consist of the daily numbers of deaths, classified by the cause of death, between January 1820 and December 1850 from N = 387 different regions (towns) in mainland Finland with the exclusion of northern areas. The time window is chosen such that there were no major famines, wars, border changes, or other potentially confounding events, which could have altered the geographical partition or the dynamics of the epidemics. The general stability achieved is advantageous in the modeling.

Although using a daily time scale would, theoretically, be ideal for modeling disease dynamics, the infrequency of deaths implies that our data lack sufficient information to study temporal dependencies at such a detailed level. This issue is likely further pronounced by the spatial heterogeneity of the data and the potentially complex lagged auto- and cross-dependencies between the infection dynamics of the diseases. Therefore, the daily counts of deaths are aggregated over time into a monthly level, decreasing the number of zero observations yet maintaining a reasonable time resolution for observing the spread of the diseases on our geographical scale. This yields a total of T = 372 time points. The counts of the observed numbers of deaths by disease, considering the aggregated data, are visualized in Fig. 1. In each case of the three infections, about 93–96% of the death counts are zero or one, and less than 2% of the counts are more than three, despite the aggregation. The reliability of the actual death counts varies considerably both temporally and spatially owing to the heterogeneous quality of the parish records and the cause of death classifications. Moreover, these deficiencies are not necessarily independent of the true number of deaths. It is also noteworthy that despite the records of baptisms and burials, there are no reliable estimates of the population sizes at the town level (Voutilainen, Helske & Högmander, 2020). Hence the intuitive idea of using local relative mortality is unfortunately beyond reach. Because of this and the aforementioned reliability issues—and since most of the counts are still either zero or one even after the aggregation to a monthly level—only the dichotomous knowledge of the death occurrence is used in the main analysis. Even so, the count data are considered in the model comparison section.

About 24% of the data concerning death occurrences in a town and a month are missing, and from 57 out of the overall 387 regions there are no observations at all. The missingness pattern is common to all the diseases, as for a particular town and month, we have observed the death counts for all the diseases or for none of them. This is because the missing data can be attributed to the absence of parish records that document all deaths. Such unavailability of data may be a result of incomplete digitization of the parish records, loss of the documents for example due to a fire in the local rectory or church, or simply because there were no deaths in a given month. Therefore, the missing data could potentially depend on the unknown regional population size, but not on the specific cause of death, given that the proportions of deaths attributable to the particular cause of death are relatively small. Nevertheless, our proposed model provides estimates of the monthly probabilities of observing at least one death also for the towns with missing observations.

Model

We construct a general model to describe the spatial and temporal dependencies both within and between the infections under study. We also want to enable exploiting other relevant information as explanatory variables. In epidemiological context, there typically occur spatial or temporal trends or seasonal effects, which can be included as separate components in the model. Due to the nature of our data, we model the probability of observing at least one death caused by a disease in a certain town in a certain month. The model consists of a trend, a seasonal effect, and a regression part reflecting the local effects of the previous state of infection in the focal town and its neighboring towns.

Formally, let $y_{i,t}^d$ denote a dichotomous variable of an event where at least one death occurs due to a disease d in a region i = 1, …, N at a time point t = 1, …, T, where N is the number of regions and T the number of time points. Let $K_x^d$ indicate the number of explanatory variables x, based on the features of the region i. Accordingly, z denotes the explanatory variables, and $K_z^d$ their number, related to the neighborhood of the region i. Thus the model for observing at least one death caused by the disease d in the region i at the time point t can be written as follows: (1)${\displaystyle y_{i,t}^d\sim \text{Bernoulli}\left({\text{logit}}^{-1}\left({\eta }_{i,t}^d\right)\right),}$ where (2)${\displaystyle {\eta }_{i,t}^d={\lambda }_i^d{\tau }_t^d+s_t^d+a_i^d+b_i^d\sum _{k=1}^{K_x^d}{\beta }_k^d{x}_{i,k}^d+c_i^d\sum _{k=1}^{K_z^d}{\gamma }_k^d{z}_{i,k}^d.}$ Here Bernoulli distribution with a logit link is chosen due to our dichotomous consideration of the occurrences of death, but other distributions with appropriate link functions can be applied for different types of response variables.

The first three terms being summed in Eq. (2) form a base level, in our case, for the probability of observing at least one death caused by the disease d at each town and month. More specifically, the first term consists of the time dependent latent factor τ_t, describing the nationwide incidence (on log-odds scale), or trend, of the disease d, and the regional adjustments or loadings λ_i, with respect to the mean level. As in general dynamic factor models, the products λ_iτ_t are not identifiable without constraints (Bai & Wang, 2015). Instead of the common approach of fixing one of the loadings λ_i to 1, we constrain the mean of the loadings to 1, enabling the interpretation of the factor τ_t as the nationwide incidence level. Due to the nature of the other terms in Eq. (2), this incidence level gives the national average log-odds of observing at least one death in an “average” town in a given month if no deaths were observed in the previous month in the focal town or in its neighboring towns. The second term s_t is a monthly seasonal effect, which is the average deviation from the nationwide incidence level, summing up to 0 over the months. The third term a_i is a regional, zero-mean constant reflecting local deviations from the nationwide incidence level τ_t due to unobserved local demographic, geographic, social or other characteristics associated with mortality.

The last two sum terms in Eq. (2) form the regression part of the model. The first sum includes the covariates x_i,k related to the focal region i, and the second sum the covariates z_i,k related to the neighboring regions. These variables have both nationwide coefficients β and γ, and their local adjustments b_i and c_i amplifying or diminishing the nationwide level. As the regional constants a_i, also the multiplicative local adjustments may account for any unobserved heterogeneity between the regions, such as the local population sizes or densities. The adjustment parameter b_i reflects the features of the focal region i, and c_i those of the neighborhood of the region i (possibly relative to i). We assume that b_i is the same for all covariates x_k, and, accordingly, c_i for all z_k, since the underlying regional characteristics modifying the nationwide mortality effects β and γ do not depend on the covariates.

In our study, there are three covariates assigned to the town i and another three to its neighbors for all diseases d = p, m, s, where p stands for pertussis, m measles, and s smallpox. The local explanatory variables x_i,k are the presences of deaths caused by the three different diseases in the previous month in the focal town, whereas the regional neighborhood predictors z_i,k are the averages of the same presences of deaths over the local neighborhood. We define two regions being neighbors when they share a border. By this definition, all the towns have at least one neighbor. Other definitions of neighborhood could be used as well, for example based on the transportation network or distance, leading to weighted averages of death occurrences.

As noted earlier, our data contain a large number of missing observations. We assume that the probability of having a missing response or predictor variable is independent of the value of the response variable. This presumption may be considered valid, given that the lack of parish records on deaths is unlikely to depend on the causes of deaths in a particular month. Under this assumption, we can model the observed data analogously to a complete-case analysis in an unbiased manner, assuming our model is correctly specified (van Buuren, 2018). This eliminates the need for multiple imputation or sampling missing observations using MCMC algorithms which do not use gradient information (i.e., algorithms capable of sampling discrete variables), both of which would be computationally unfeasible in our Bayesian spatio-temporal context.

In practice, the complete-case analysis in our context means that to use an observation as a response, we require that both the current response variable as well as all of the related covariates are observed. If any of them is missing, we omit the particular combination of town and month as a response. On the other hand, when calculating the neighborhood covariates z_i,k, which in our case are the averages over the observations within the neighborhood, we omit the neighbors with missing observations so that they are not included even as a denominator in the evaluation of the mean. If all neighbor observations are missing, the corresponding covariate is defined as missing.

We model the latent factor ${\tau }_t=\left({\tau }_t^p,{\tau }_t^m,{\tau }_t^s\right)$, the temporal process describing the baseline of the nationwide incidence rates, as an intercorrelated random walk, τ_t+1 ∼ N(τ_t, Σ). Here Σ is an unconstrained 3 × 3 covariance matrix parametrized using the standard deviations ${\sigma }_{\tau }^d$ and the correlation matrix R. In our application this latent process, together with the regional constant and the monthly effect, can be interpreted as the probability to observe at least one new death in a particular town when no deaths caused by any of the three diseases were observed in the previous month in the focal town or in its neighborhood.

In the Bayesian modeling framework, we need prior distributions for all the parameters to be estimated. The incidence factors τ^d follow a N(−2, 2²) prior at the first time point and form a random walk at the later time points. For the correlation matrix R we use an LKJ(1) prior with a Cholesky parametrization (Lewandowski, Kurowicka & Joe, 2009), i.e., a uniform prior over valid 3 × 3 correlation matrices. For the regional parameters we apply Gaussian priors: ${\lambda }_i^d\sim \text{N}\left(1,{\sigma }_{{\lambda }^d}^2\right),a_i^d\sim \text{N}\left(0,{\sigma }_{a^d}^2\right),b_i^d\sim \text{N}\left(1,{\sigma }_{b^d}^2\right)$, and $c_i^d\sim \text{N}\left(1,{\sigma }_{c^d}^2\right)$. The prior means of the local adjustments $b_i^d$ and $c_i^d$ are set to 1, which enables us to interpret ${\beta }_k^d$ and ${\gamma }_k^d$ as nationwide effects as is the case with additive multilevel models having population-level and group-level effects. While we use a hard equality constraint for the mean of the ${\lambda }_i^d$s to ensure the identifiability and more efficient estimation of the model, the hierarchical priors for $b_i^d$ and $c_i^d$ are sufficient for their identifiability. For the unknown deviations σ_τ^d, σ_λ^d, σ_a^d, σ_b^d, and σ_c^d we assign Gamma priors with shape parameter 2 and rate parameter 1. The nationwide coefficients ${\beta }_k^d$ and ${\gamma }_k^d$ have N(0, 2²) priors. The seasonal effects $s_t^d$ follow a standard normal prior with the aforementioned sum-to-zero constraint. The priors can be seen as weakly informative, and they are chosen primarily to enhance the computational efficiency (Banner, Irvine & Rodhouse, 2020).

Naturally, any of the components in the model could be excluded by setting the corresponding coefficients or standard deviations to zero. Our Bayesian model encompasses all such simplified alternatives, with the corresponding model and parameter uncertainty reflected by the estimated posterior distributions, leading to more truthful uncertainty estimates compared to merely imposing prior constraints on certain effects to be zero.

Model comparison

While our main interest was studying the past spatio-temporal dynamics of infections and disease associations within and between the diseases, we also examined the necessity and reasonableness of modeling the disease interdependencies and the response aggregation. We compared our model with a corresponding one without the dependencies between the infections by excluding the other diseases as explanatory variables and omitting the correlation between the incidence factors τ_t in the simpler model. Since the original data contained the numbers of deaths instead of the dichotomous aggregates we used as a response, we also estimated corresponding models with the difference of using the counts as a response and a negative binomial distribution to model them. Additionally, the briefly aforementioned model without a seasonal component was included in the comparison in the case of both types of responses. This resulted in six different model versions for comparison: dependent diseases, independent diseases, and dependent diseases without a seasonal effect, each for both Bernoulli and negative binomial distributions.

The negative binomial model can be formally written as (3)${\displaystyle y_{i,t}^d\sim \text{NB}\left(\text{exp}\left({\eta }_{i,t}^d\right),\text{exp}\left({\alpha }_{{\varphi }^d}+{\varphi }_i^d\right)\right),}$ where the mean parameter ${\eta }_{i,t}^d$ is defined as in Eq. (2), and the nationwide dispersion parameters α_ϕ^d and the local dispersion parameters ${\varphi }_i^d$ depend on the response disease. The priors are the same as with the Bernoulli model, with the addition of α_ϕ^d ∼ N(0, 1²), ${\varphi }^d\sim \text{N}\left(0,{\sigma }_{{\varphi }^d}^2\right)$, and σ_ϕ^d ∼ Gamma(2, 1). From the negative binomial model, we could then compute our quantity of interest, the probability of observing at least one death in a specific town and month, which could be compared with the corresponding estimates of the Bernoulli model.

As a scoring rule for the model comparison, we used the expected log predictive density (ELPD), which measures the goodness of the entire predictive distribution (Vehtari, Gelman & Gabry, 2017). The ELPD was estimated via an approximate leave-one-out cross-validation using the R package loo (Vehtari et al., 2023). We left out one month and town from all the diseases at a time to estimate the ELPD. Models with higher values of ELPD correspond to greater posterior predictive accuracy for predicting new data points compared to models with lower ELPD values.

According to the differences in the ELPDs in Table 2, the best performing model is the one utilizing Bernoulli distribution and considering the diseases dependent. Omitting the dependencies results in the second-best model, with the difference in ELPD over three times the standard error. As could be expected, omitting the seasonal effect further impairs the model. When it comes to the negative binomial models with counts as responses, the order of the dependent, independent, and seasonless models is the same. The Bernoulli models outperform the negative binomial ones in all cases. Overall, we see that directly using the dichotomized data versus modeling the count data has a greater impact than considering the infection dependencies or seasonality in our model. However, even though in terms of predictive performance the differences between different Bernoulli models are relatively small, we used the most complex model in our main analysis. This is in line with the common Bayesian paradigm of incorporating the uncertainty of the model structure in the model (Vehtari & Ojanen, 2012).

For the dependent and independent models, we performed additional prediction checks by discarding the last two years of the data and estimating the probabilities for those years. We also calculated the ELPDs considering the removed years, see Table 2. The modifications without the seasonal effect were not included in this comparison due to their already evident poor performance and the fact that they were originally fitted merely to investigate the importance of the obvious seasonal variation. The posterior means were quite similar in all cases, but the posterior intervals were wider for measles and smallpox in the case of the negative binomial model, as can be seen from Fig. 7.

Overall, the results of the Bernoulli, as well as the negative binomial, models seem to indicate similar interdependencies between the diseases. In the case of the negative binomial model, the estimates of all the time and town invariant parameters with their prior and posterior intervals are shown in Supplementary Table 2 on GitHub. Also, figures corresponding to the ones representing the results of the Bernoulli model are available in GitHub (Supplementary Figures 3-9).

Furthermore, to inspect the effect of the earlier disease history, a Bernoulli model according to Eqs. (1) and (2) was fitted with additional lags of two months for the three focal and three neighbor covariates. The approach resulted in having six scalar regression coefficients (three βs and three γs) more for each response disease than in our main model. The results are well aligned with the one-month lag model. The amount of available observations decreases by about three percentage points when introducing the two-month lag since we must know not only the previous observation but also the one preceding that. Thus this model is not completely comparable to our main model. Using the same data for both one-month and two-month lag models, the model with two-month lag performs better, measured with the ELPD: one-month lag model results in ELPD_diff =  − 989.87 and SE_diff = 53.89, compared with the two-month lag model. Nevertheless, there are some convergence and efficacy issues. Out of the 20, 000 iterations, there are 28 diverging ones which can potentially bias the results, thus they are not completely reliable (Betancourt, 2018). The outperformance of the two-month lag model in the sense of the ELPD might be also related to smaller variance regardless of the possible bias, which is consistent with the seemingly better fitting predictions gained from our main model than the two-month lag model (see Supplementary Figure 10). The full results are in GitHub in Supplementary Table 3 and Supplementary Figures 10-16.

Additionally, we fitted a model using our data aggregated on a weekly level, which increases the amount of missing data from 24% in the monthly data to 48% in the weekly data. We used lags from 1 to 4 weeks to cover as much delayed effect as with our main model. Due to the unequal number of weeks per calendar month, incorporating the monthly effect is not straightforward, so we omitted the seasonal effect. Unfortunately, this model did not converge, potentially because of the increased amount of missing data, or the complex dependency structures due to varying sub-monthly incubation and time-to-fatality times. The weekly data and model code are available in the supplementary material.