Force of infection as a catalytic model

The force of infection is the key quantity governing disease transmission within a given population and is defined as the instantaneous per capita rate at which susceptible individuals acquire infection (i.e. the hazard of infection). It reflects both the degree of contact between susceptible and infected individuals and the transmissibility of the pathogen per contact. In cases where contact is age dependent, the force of infection is itself a function of age analogous to the hazard being a function of time in a survival model [Reference Farrington21].

Formally, let P(a) denote the probability that an individual susceptible at birth is still susceptible (i.e. uninfected) at age a; this is identical to the cumulative survival function in a survival analysis. Let λ(a) denote the force of infection at age a, which represents the rate of infection of the susceptible population P(a). In terms of a survival model, the force of infection λ(a) represents the instantaneous infection hazard. Therefore, for lifelong infections (or infections conferring lifelong immunity) that do not significantly affect mortality, λ(a) can be defined as:

(1)$$\lambda (a) = - \displaystyle{1 \over {P(a)}}\displaystyle{{dP(a)} \over {da}}$$

under the assumption that the population is in dynamic equilibrium (meaning that P(a) is constant over calendar time and not changing by birth cohort).

Equation (1) specifies a so-called catalytic epidemic model, first defined by Muench [Reference Muench22]. The term catalytic model derives from its origins in chemistry, where these models were used to study chemical reaction kinetics. Catalytic epidemic models have been widely used to estimate the force of infection of different infectious disease such as measles, rubella, mumps, hepatitis A and yellow fever using seroprevalence epidemiological data from settings where these diseases do not significantly affect mortality [Reference Hens16].

Catalytic models can easily be written in terms of the cumulative probability of infection by age a, F(a) = 1 − P(a), which more directly corresponds to how seroprevalence data are collected [Reference Farrington21, Reference Muench22]. In terms of F(a), the force of infection, λ(a), becomes:

(2)$$\lambda (a) = \displaystyle{1 \over {1 - F(a)}}\displaystyle{{dF(a)} \over {da}}.$$

Note that rearranging Equation (2) we can solve for F(a) as a function of λ(a), which is equivalent to a survival model where the probability of infection by age a is the cumulative distribution function of the time to infection:

(3)$$F(a) = 1 - {\rm exp}\left[ { - \mathop \int \limits_0^a \lambda (s)ds} \right].$$

Seroprevalence data of H. pylori in Mexico

In 1987–1988, the National Seroepidemiological Survey was conducted in Mexico as a nationally representative seroprevalence study of H. pylori infection [Reference Torres15, Reference Tapia, Gutiérrez and Sepúlveda23]. The survey was designed to represent the country by including all 32 states of Mexico, all ages, all socio-economic levels and both sexes [Reference Tapia, Gutiérrez and Sepúlveda23]. In total, 32 200 households were surveyed and more than 70 000 serum samples were collected. Of these samples, 11 605 individuals were used for H. pylori testing, and out of these, 7 720 (66%) were seropositive for H. pylori. H. pylori infection was determined by ELISA detection of IgG antibodies to specific H. pylori antigens. H. pylori prevalence in Mexico varies by state and can be as low as 48% in the state of Colima and as high as 83% in the state of Baja California Sur (Fig. 1).

Fig. 1. Map of Mexico with observed prevalence of Helicobacter pylori by state in 1987–88.

Functional form of H. pylori infection in Mexico

The National Seroepidemiological Survey was collected prior to widespread antibiotic treatment of H. pylori in 1991, so we can assume that the population is in dynamic equilibrium. The age-specific prevalence measured nationally is shown in Figure 2 and it is clear that the trend is monotonically non-decreasing.

Fig. 2. Empirical prevalence of Helicobacter pylori in Mexico by age nationally and in three different states with different amount of data.

Based on the shape of these data, we assume that the state-specific cumulative probability of H. pylori infection, F _i(a), for each state i follows a modified exponential functional form:

(4)$$F_i(a) = \alpha _i(1 - e^{ - \gamma _ia}),$$

where 0 ⩽ α _i ⩽ 1 and γ _i ⩾ 0 are the state-specific asymptote and rate parameters of the catalytic model. The asymptote parameter α is often interpreted as the proportion of the population that can be infected (i.e. the susceptible population), if it is assumed that it is the exhaustion of susceptible individuals that leads to the eventual plateau in the proportion infected at older ages [Reference Hens16]. With this interpretation, the rate parameter γ is then the constant rate at which these susceptible individuals become infected. We allow α and γ to vary by state to reflect the heterogeneity in H. pylori prevalence trends resulting from differences in socio-economic and environmental factors. Note that for all states, we assume that no children are infected at birth (F(0) = 0) since there is no biological evidence for vertical transmission of H. pylori from mother to child during pregnancy, labour or delivery [Reference Espinoza Monje and Marin Vega24].

Based on this functional form for the cumulative probability of infection, we used Equation (3) to solve for the state-specific force of infection, λ _i(a), in terms of α _i and γ _i:

(5)$$\lambda _i(a) = \displaystyle{{\alpha _i\gamma _ie^{ - \gamma _ia}} \over {1 - \alpha _i(1 - e^{ - \gamma _ia})}}.$$

Note that if the entire population is eventually infected in a given state (i.e. α _i = 1), λ _i(a) = γ _i and Equation (4) reduces to a traditional exponential model in survival analysis.

Estimation of the force of H. pylori infection

We used statistical methods to estimate the parameters α and γ that best fit the cumulative probability of infection in each state and then calculated the corresponding force of infection. Due to the high degree of variation in sample composition for each state in the National Seroepidemiological Survey, traditional non-linear least-squares estimation yielded imprecise, unstable or unrealistic estimates of α and γ. This could translate into estimates of force of infection that are implausibly high or low with wide confidence intervals (CIs) when only a small number of people are sampled (and, for instance, none are infected). For example, some states had individuals sampled at all ages while others had limited data for certain ages. In Figure 3, we show the observed prevalence in four different states with varying amounts of data and the predicted prevalence for these states using the estimates from non-linear least-squares estimation (in blue solid lines). Nuevo Leon and Mexico City have more data and less variability in terms of disease status compared with the states of Baja California Sur and Nayarit. The state of Baja California Sur is a state where all adults sampled older than 20 years old (n = 24) are positive for H. pylori. The resulting parameter estimates are implausibly high based on our current knowledge of H. pylori.

Fig. 3. Empirical prevalence of Helicobacter pylori in Mexico by age in four different states with different amount of data, together with model-predicted prevalence using non-linear least-squares (NLS) and hierarchical Bayesian estimation methods. The upper two reflect states with data for all ages and the lower two reflect states with limited amount of data for some ages. The red solid line denotes the model-predicted posterior mean and the red dashed lines denote the 95% credible bounds. The black dotted horizontal line shows the average national H. pylori prevalence.

To solve these issues, we instead implemented a hierarchical Bayesian non-linear model that ‘borrows’ information from all states when estimating state-specific effects. Such models have been described in detail elsewhere [Reference Gelman and Hill25–Reference Witte27] and are known to reduce mean squared error of model estimates in some settings. We are unaware of these models having been applied to catalytic models to estimate force of infection. We first used hierarchical Bayesian non-linear models to estimate the posterior distribution of the α and γ parameters of model (4) for each of the 32 states. We then estimated the aggregated and state-specific force of infection of H. pylori in Mexico.

Hierarchical non-linear Bayesian model

We assume that for each state i = 1, …, 32 the number of infected individuals y _i(a) at age a follows a binomial distribution where n _i(a) is the sample size and p _i(a) is the cumulative probability of infection (i.e. F _i(a)) in state i:

(6)$$\matrix{ {y_i(a)} & { \sim {\rm Binomial(}n_i(a),p_i(a))} \cr {\,p_i(a)} & { = F_i(a) = \alpha _i(1 - e^{ - \gamma _ia})} \cr {\left[ {\matrix{ {{\rm logit(}\alpha _i{\rm )}} \cr {{\rm log(}\gamma _i{\rm )}} } } \right]} &\!\!\!\!\!\! \sim { {\rm MVN}\left( {\left[ {\matrix{ {\beta _1\; } \cr {\beta _2} } } \right],\Sigma } \right)} \cr \Sigma & { \sim {\rm IW}\left( {\left[ {\matrix{ {{10}^3} & 0 \cr 0 & {{10}^3}} } \right],2} \right)}} .$$

Note that each of the 32 states has its own set of parameters α _i and γ _i, but these state-specific parameters ‘borrow’ information from each other by assuming that they collectively follow a multivariate normal (MVN) distribution centred at β ₁ and β ₂. To compute the average national-level proportion of susceptibles and infection rate, α ⁰ and γ ⁰, we applied a logit and log transformation to the β ₁ and β ₂ parameters, respectively.

For the covariance matrix, we used an inverse Wishart (IW) distribution of two degrees of freedom. The IW is a widely used prior distribution for covariance matrices; the large diagonal entries and the zero off-diagonal entries in the matrix represent a weakly informative prior on the variances and no prior correlation between α and γ parameters [Reference Gelman and Hill25]. Weakly informative normal prior distributions (μ = 0 and σ = 10³) were specified for β ₁ and β ₂ parameters.

Parameters of model (6) were estimated through Gibbs sampling, using Just Another Gibbs Sampler (JAGS) [Reference Plummer, Hornik, Leish and Zeileis28]. We ran two Markov chain Monte Carlo (MCMC) chains. For each chain, 50 000 iterations were used to generate a posterior distribution after a 10 000 iterations burn-in period. The mean estimates and 95% credible intervals (CR) were derived from posterior distributions. For each chain, we used starting values computed from the solution of model (4) using non-linear least-squares. We checked each parameter's chains to assess convergence. To confirm whether the model provides a reasonable fit of the data, we plotted the observed prevalence by age with 95% CIs, together with the model-predicted prevalence.

The force of H. pylori infection was estimated within the MCMC model by evaluating Equation (5) at each parameter set drawn from the posterior distribution, respectively. We then calculated their posterior predicted mean and 95% CR. All statistical modelling was conducted in JAGS [Reference Plummer, Hornik, Leish and Zeileis28] and R version 3.2.4 (http://www.r-project.org).