Computing the final epidemic size distributions of a multi-type Galton--Watson process

The gist The present study proposes a framework for computing the final size distribution of multi-type Galton–Watson processes using a method based on Cauchy integral contours, providing both theoretical justification and practical numerical procedures for statistical inference in epidemiological modeling.

Introduction and Motivation

Stochastic epidemic models are applied to analyze real-time transmission dynamics, demonstration of transmission heterogeneity, and determination of the eventual size of epidemics 1. When used with surveillance-based cluster size data, the Galton–Watson process (GWP) can elicit information about the extent of transmission heterogeneity even when each transmission process is not directly observable 1. For multi-type GWPs, real-world applications remain limited possibly because of the absence of conceptually and practically straightforward approaches for deriving the closed-form solution of the final size distribution 1. The paper proposes a framework for computing the final size distribution using a method for the choice of Cauchy integral contour 1.

Theoretical Framework and Final Size Equation

The application of single-type GWPs to epidemic cluster size data is facilitated by the fact that the explicit form of the final size distribution for a negative binomial offspring is known and can be directly used for likelihood-based statistical inference 2. For multi-type GWPs, the criticality of multi-type GWPs is defined by the spectral radius of the Jacobian matrix of G(z) at z = 1, or K = JG(1) 2. The final size equation for the PGF of a GWP originating from a “parent” of type i is (2.4) H(z) = diag(z)G(H(z)) 4.

Cauchy Integral Calculation and Contour Selection

The coefficients of H(i) in Eq. (2.3) can be written as a multivariate Cauchy integral representation 4. The critical step enabled by Proposition 3.1 is the principled choice of the contour radii r = (r1,..., rn) to guarantee the existence and uniqueness of the integrand 4. For subcritical or critical processes, r = (0.95, 0.95) is used 4. For supercritical processes, a vector r is found such that ρ(JG(r)) = 0.975 4.

Numerical Procedure and Verification

The numerical framework combines standard methods for calculating Cauchy integrals and for optimization 4. The calculation of PGF coefficients is based on a multivariate Cauchy integral evaluated on multidimensional grids using the discrete Fourier transform (DFT) 5. To solve the final size equation, the system Fz(h) = h − diag(z)G(h) = 0 is solved as a minimization problem of ψz(h) =∥Fz(h)∥22 or the product of Fz(h) and its conjugate transpose using Newton’s method with backtracking 5. Proposition 3.1 ensures both the existence of a unique solution and the regularity of the Jacobian matrix JFz(h) within the region h ≤ r 5.

Case Study and Parameter Inference

The framework is applied to Middle East respiratory syndrome (MERS) cluster size data, modeling a two-type GWP based on age groups 8. The analysis estimated parameters using the Markov Chain Monte Carlo (MCMC) method, yielding estimates such as R = 0.77 and k = 0.30 12. The results show that the likelihood surface is smoothly connected across the criticality threshold log(R) = 0 8, and a trade-off between mean transmission and dispersion is observed in the supercritical region 8.

Conclusion and Advantages

The study provides a straightforward framework for evaluating Cauchy integral representations of final size PGF coefficients for the multi-type GWP 13. The approach does not involve the change of variables required in the Lagrange–Good method, and it allows seamless handling of subcritical and supercritical processes within the same framework 13. It also provides a firm basis for obtaining the numerator of the integrand in Eq. (4.1) without explicit discussion on convergence or uniqueness of solutions 13. The framework is applicable to multi-type GWPs with general offspring PGFs not limited to the negative binomial or negative multinomial family 13. It is also encouraging that we could actually infer epidemiological parameters without imposing an assumption on the criticality of the underlying GWP from both our simulated data and the MERS cluster size data 13. The study has revealed that extinction probabilities for GWPs with negative multinomial offspring are monotonically decreasing functions of both the elements of the mean offspring matrix and the overdispersion parameters for all types of hosts 13.

Limitations

Several limitations merit further investigation, including the need for efficient computation in higher-dimensional complex spaces for greater numbers of host types 13, and the lack of focus on the identifiability of criticality itself 13. The study has focused on estimating the spectral radius of the mean offspring matrix and the overdispersion parameter from cluster size data 13. Future work should explore extending the framework to real-time forecasting or to joint inference on additional epidemiological parameters 13.

How it works

The application of single-type GWPs to epidemic cluster size data is facilitated by the fact that the explicit form of the final size distribution for a negative binomial offspring is known and can be directly used for likelihood-based statistical inference 2. The framework relies on calculating PGF coefficients using a multivariate Cauchy integral (4.1) 4, which is then evaluated on multidimensional grids along the contour Γz using the discrete Fourier transform (DFT) 5. The critical step enabled by Proposition 3.1 is the principled choice of the contour radii r = (r1,..., rn) to guarantee the existence and uniqueness of the integrand 4.

Numerical Procedure

To solve for H(i) on each grid point via DFT, one must find the unique solution h to Eq. (3.1), which is equivalent to finding roots of the system Fz(h) = h − diag(z)G(h) = 0 5. This is solved as a minimization problem of ψz(h) =∥Fz(h)∥22 or the product of Fz(h) and its conjugate transpose using Newton’s method with backtracking 5. Proposition 3.1 ensures both the existence of a unique solution and the regularity of the Jacobian matrix JFz(h) within the region h ≤ r 4.

Parameter Inference

The framework enables likelihood-based inference by calculating marginal likelihoods LD using formulas such as (6.4) LD = X2i=1 p(ID = i)pi(d1, d2;θ) 9. Parameters are estimated using the Markov Chain Monte Carlo (MCMC) method, and the posterior estimates for R, k, S, and π are presented in Table 2 14. The correlation plot suggests a ridge toward the high log R - log k direction that extends into the supercritical domain 15.

1. This AI system can perform parameter inference for multi-type Galton–Watson processes using an unconditional likelihood framework to seamlessly handle the sub- and supercritical regimes, allowing it to infer epidemiological parameters without imposing an assumption on criticality from both our simulated data and the MERS cluster size data.

2. The system can compute the final size distribution coefficients for multi-type GWPs using a novel method based on the Cauchy integral representation, providing a straightforward framework that is applicable to general offspring PGFs not limited to the negative binomial or negative multinomial family.

3. The AI can conduct likelihood-based inference on real-world epidemic data, such as MERS cluster size distributions, by estimating parameters like the basic reproduction number (R) and dispersion parameter (k) through Markov Chain Monte Carlo (MCMC) methods.

4. The system can identify potential parameter identifiability challenges in supercritical domains by analyzing the "negative correlation between R and k for posterior draws in the supercritical (R > 1) region, which suggests a trade-off that stems from the extinction probabilities."