Measles surveillance data

In Vietnam, measles is one of the 20 infectious diseases that are required to be notified online to the Electronic Communicable Disease Surveillance System (ECDS) within 24 h post clinical diagnosis. In the case report form, several fields including demographic characteristics, date of illness onset, date of hospitalisation or medical examination of in- and out-patients are obligated to be reported. More information about the disease surveillance and reporting system in Vietnam can be found in [18, 19]. In this study, we analysed measles cases reported between 1 January 2018 and 30 June 2020 in the Southern region. On 11 October 2020, 26 047 individual cases were extracted from the ECDS, accessed by the Pasteur Institute in Ho Chi Minh City, the public health institute that manages the disease surveillance system in the South. For the age-structured spatio-temporal analysis of the outbreak, we categorised individuals into four age groups: 0–4 (children), 5–14 (school children), 15–24 (adolescent and young adults) and 25+ (adults) years of age. Overall, the median age at disease onset was 3 years, ranging from <1 year to 84 years. The most affected age group was 0–4 years, which accounted for 61.4% of measles cases while the 5–14 years, 15–24 years and 25+ years groups represented 22.0%, 4.5% and 12.1% of the total cases, respectively. Based on the date of onset, we aggregated the daily number of cases by age group in each province in which the cases resided. Figure 1 depicts the evolution of daily counts of measles infections and monthly incidence per 100 000 population by age group and Figure 2 presents maps of the age-specific cumulative incidence per 100 000 individuals across all provinces.

Fig. 1. Evolution of age-specific measles cases by date of symptom onset (left axes) and incidence per 100 000 population by month (blue lines, right axes) in Southern Vietnam, 1 January 2018 to 30 June 2020. Lunar New Year in 2019 and 2020 are highlighted in red.

Fig. 2. Maps of age-specific cumulative incidence per 100 000 population in Southern Vietnam, 1 January 2018 to 30 June 2020.

Population data

Age- and province-specific population data were obtained from the census data in 2019 [20]. The Southern population was approximately 36 million in 2019, and in each age group, the population fraction was 6.7% for 0–4 years, 14.9% for 5–14 years, 14.3% for 15–24 years and 64.2% for 25+ years. We assumed that the total population was constant over the period 2018–2020.

Social contact data

To reflect the amount of mixing between age groups, we used social contact data, adapted from an empirical contact matrix in a survey in the Red River Delta region of Northern Vietnam in 2007 [Reference Horby11]. The matrix was transformed because of the difference in demographic structures between Northern and Southern Vietnam and hence, directly using the original contact matrix would not be valid in our study. First, we extracted the social contact patterns aggregated to the age groups of interest from the Social Contact Rates (SOCRATES) Data Tool (http://www.socialcontactdata.org/socrates/) [Reference Willem21]. This age-structured contact matrix ${\boldsymbol C} = ( {{\boldsymbol c}_{{g}^{\prime}g}} )$ provided the average non-negative number of contacts of a person in age group g ^′ (rows) with a contact in age group g (columns) in 1 day (Fig. 3a) aggregated over weekdays or weekends, contact duration, physical or non-physical contacts and gender. Next, we projected the social contact matrix for Southern Vietnam ${\boldsymbol C}_{( {\boldsymbol P} ) } = ( {{\boldsymbol c}_{( {\boldsymbol P} ) {g}^{\prime}g}} )$ using the density correction method proposed by Arregui et al. [Reference Arregui22]. The projected contact matrix is a product of an intrinsic connectivity matrix C(N/N _g) and the fraction of individuals in the age group of the contact $N_g^{\prime} /{N}^{\prime}$, where N _g and $N_g^{\prime}$ are the demographic structures in 2009 (Red River Delta region) and 2019 (Southern region), respectively. Note that because the age-structured population numbers of Red River Delta region are not available for 2007, we used population of Red River Delta region from the 2009 census [20]. The obtained contact matrix is shown in Figure 3b.

Fig. 3. (a) Original age-structured contact matrix C estimated in Northern Vietnam anno 2007 aggregated to the age groups of interest and (b) the age-structured contact matrix projected for Southern Vietnam C_(P) based on (a). The entries contain the mean number of contacts made by one participant per day. (c), (d), (e) and (f) refer to the power transformation of row-normalised contact matrix C_(P) for different values of κ.

Age-structured spatio-temporal analysis

In general, we leveraged an endemic–epidemic modelling framework for multivariate infectious disease counts first introduced by Held et al. [Reference Held, Höhle and Hofmann16] and extended in a series of publications [Reference Herzog, Paul and Held5, Reference Meyer and Held17, Reference Paul and Held23–Reference Held, Meyer and Bracher26, Reference Geilhufe38, Reference Bracher and Held39]. The framework subsequently incorporated the age-structured contact matrix (possibly adjusted) to better understand disease spread in the scenario of heterogeneous mixing [Reference Meyer and Held17].

Formally, let Y _grt denote the number of cases in age group g = 1, …, G in province r = 1, …, R at time t = 1, …, T. Conditional on the number of cases at the previous time point t − 1, the counts are assumed to follow a negative binomial distribution with conditional mean μ _grt:

(1)$$\mu _{grt} = e_{gr}\nu _{grt} + \phi _{grt}\mathop \sum \limits_{{g}^{\prime}, r^{\prime}} \lfloor{{\boldsymbol c}_{( {\boldsymbol P} ) g^{\prime}g}{\cal W}_{{r}^{\prime}r}\;} \rfloor Y_{{g}^{\prime}, {r}^{\prime}, t-1}$$

and a variance μ _grt(1 + μ _grtψ _g) with a group-specific overdispersion parameter ψ _g > 0 [Reference Meyer and Held17, Reference Held, Meyer and Bracher26]. Note that if ψ _g = 0, the distribution simplifies to the Poisson distribution. The mean μ _grt is decomposed into endemic and epidemic components. The former component exhibits baseline patterns. The latter component involves an autoregressive effect that links cases at time point t in unit r with observations at the previous time point t − 1 and in units r ^′ = 1, …, R. Specifically, the non-negative parameters ν _grt and ϕ _grt are modelled as log-linear predictors:

(2)$$\log ( {\nu_{grt}} ) = \alpha _g^{( \nu ) } + \alpha _{mekong}^{( \nu ) } + \beta _{lunar}^{( \nu ) } x_t + \beta _{trend}^{( \nu ) } t + \beta _{{\rm sin}}\sin ( {\omega t} ) + \beta _{{\rm cos}}{\rm cos}( {\omega t} ) $$

(3)$$\log ( {\phi_{grt}} ) = \alpha _g^{( \phi ) } + \alpha _r^{( \phi ) } + \beta _{lunar}^{( \phi ) } x_t + {\rm \tau }\log ( {e_{gr}} ) .$$

The above two equations contain age-specific fixed effects ($\alpha _g^{( \nu ) }$, $\alpha _g^{( \phi ) }$). Because we deemed that fewer cases were reported during the Lunar New Year, we included an indicator for the holiday in 2019 and 2020 with coefficient β _lunar (x _t = 1 for dates from 2 to 10 February 2019 and from 23 to 29 January 2020, otherwise x _t = 0). In the endemic component, to adjust for the possibly different number of individuals at risk in each province and age group, we included the population size e _gr as an offset [Reference Meyer and Held25]. To allow for the differences amongst the two main administrative regions, $\alpha _{mekong}^{( \nu ) }$ is included to resemble the effect of the Mekong River Delta region. We also assumed that the disease incidence varies with a linear time effect (β _trend) and an overall seasonal sine–cosine term where the sinusoidal wave of frequency ω identified as 2π/365 for daily continuous measurement [Reference Held and Paul24]. In the epidemic component, we allowed for province-specific effects $\alpha _r^{( \phi ) }$ and accounted for population size e _gr to quantify how ‘attraction’ to a province r scales with population size in group g, in which the strength of population scaling factor τ is to be estimated [Reference Meyer and Held25, Reference Xia, Bjørnstad and Grenfell27].

To determine transmission weights from age stratum g ^′ to age stratum g (i.e. ${\boldsymbol c}_{( {\boldsymbol P} ) g^{\prime}g}$), and from area r ^′ to r (i.e. ${\cal W}_{{r}^{\prime}r}$), the product ${\boldsymbol c}_{( {\boldsymbol P} ) g^{\prime}g}{\cal W}_{{r}^{\prime}r}$, which is row-normalised, i.e. $\sum\nolimits_{g, r} {\lfloor{{\boldsymbol c}_{( {\boldsymbol P} ) g^{\prime}g}{\cal W}_{{r}^{\prime}r}} \rfloor } = 1$ was introduced in the epidemic component. In ideal circumstances, to best reflect the transmission between strata, the matrices of contact and mobility should be displayed by age and province. Nevertheless, such data sources are not easily available as collecting contact and movement patterns is cumbersome. In our study, we took the overall estimate of contact data for age-group weights and used the power law approximation for the spatial weights. First, the age-group weights ${\boldsymbol c}_{( {\boldsymbol P} ) g^{\prime}g}$ are row-normalised and then raise it to the power κ ≥ 0, i.e. ${\boldsymbol C}_{( {\boldsymbol P} ) }^\kappa$ [Reference Meyer and Held17]. In an easy interpretation, the limit κ = 0 corresponds to no mixing between different age groups, i.e. the diagonal contact matrix C_(P) = I (Fig. 3c). When κ = 1, the contact matrix represents the given projected contact matrix (Fig. 3d). As κ → ∞, the transmission from an infected person to any individual of any age group has the same distribution with other groups regardless of the group they are in [Reference Meyer and Held17]. We also consider homogenous mixing scenario in the epidemic component. Second, the non-negative weight ${\cal W}_{{r}^{\prime}r}$ in the epidemic component describes the strength of transmission between geographical units. In the absence of mobility data, it can be estimated using a power law formulation in terms of adjacency order $o_{{r}^{\prime}r}$, which is a discrete distance measure of neighbourhood order between unit r ^′ and r [Reference Meyer and Held25]. The power law weights ${\cal W}_{{r}^{\prime}r} = ( {o_{{r}^{\prime}r} + 1} ) ^{{-}d}$, where d > 0 is the decay parameter to be estimated, thus give unit weight to local transmission when r ^′ =  r and then decay to promote the spatial transmission from unit r ^′ to unit r. The power law weights can be age-dependent (replacing d by $d_{{g}^{\prime}}$) [Reference Meyer and Held17]. In this study, $o_{{r}^{\prime}r}$ ranges from 0 to 7.

All procedures were performed using R software version 4.0.5, packages surveillance version 1.19.1 [Reference Meyer, Held and Höhle28] and hhh4contacts version 0.13.1 [Reference Meyer and Held17]. In each model, maximum likelihood estimates of parameters and 95% confidence intervals (95% CIs) were obtained numerically. Model selection is performed according to the smallest Akaike information criterion (AIC) value.

Sensitivity analysis

We performed a sensitivity analysis using weekly aggregation of the surveillance data. We also ran another sensitivity analysis to assess the impact of different forms of contact matrix on our results, including the original contact matrix, and the per capita contact rates (i.e. dividing the mean number of contacts per day per participant in group g ^′ to the Vietnamese population size in 2009 and in 2019 [20] in contact group g).

Ethical consideration

As part of public health surveillance system in Vietnam, case-based data of measles were routinely collected for disease control purposes. Anonymised data, i.e. without identification of patient information, were provided for use in this study. Therefore, this study did not require ethical approval.