Region-wide synchrony and traveling waves of dengue across eight countries in Southeast Asia

1.1. Data Collection.

We collected monthly dengue surveillance data at the provincial level for Thailand, Cambodia, Laos, Vietnam, Malaysia, Singapore, the Philippines, and Taiwan. Data for Thailand were provided by the Thai Ministry of Public Health Bureau of Epidemiology (1968–2010). Data for Cambodia were provided by the Cambodia Ministry of Health (MOH) National Dengue Control Program (1998–2009) and supplemented with data from the WHO DengueNet Database (1998–2010) and the WHO Western Pacific Regional Office (1998–2001) where needed. Data for Laos were provided by the Lao People's Democratic Republic MOH National Center for Laboratory and Epidemiology (1998–2010) and supplemented with data from the WHO DengueNet Database (1998–2007) and the WHO Western Pacific Regional Office (1998–2001) where needed. Data for Vietnam were provided by the National Institute of Hygiene and Epidemiology (1994–2010). Data for Malaysia were provided by the Malaysia MOH Disease Control Division (1988–2010). Data for Singapore were provided by the Singapore National Environment Agency (1993–2010). Data for the Philippines were provided by the Philippines MOH National Epidemiology Center (1993–2010), and data for Taiwan were provided by the Taiwan Centers for Disease Control (1998–2010). Dengue surveillance systems were similar across countries and mostly passive, except for the Philippines and Cambodia (sentinel system). In most countries, dengue is predominantly reported among children. Cambodia was the only country for which dengue surveillance was explicitly restricted to children <15 y old. Despite these slight differences, population-level disease patterns were comparable across all countries. For each country, data were collected at the provincial level. We defined a province as the smallest geographical area for which data were available (the first or second administrative level). Where needed, surveillance reports for dengue fever and dengue hemorrhagic fever were combined into the total number of dengue cases. Population counts were obtained for provinces by year for each country. Population counts for Thailand were provided by the Thai Ministry of Public Health Bureau of Epidemiology (1968–2010); population counts for Cambodia were provided by the National Institute of Statistics (1998) and the MOH National Dengue Control Program (2002–2010). Population counts for Laos were obtained from the Laos Statistics Bureau (1995 and 2000–2010). Population counts for Vietnam were provided by the General Statistics Office of Vietnam (1995–2010). Population counts for Malaysia were provided by the Department of Statistics (2000 and 2003–2010). Population counts for Singapore were obtained from the US Census Bureau International Database (25). Population counts for the Philippines were obtained from the Philippines Census Bureau (2000, 2005, and 2010), and population counts for Taiwan were provided by the Taiwan Centers for Disease Control (1998–2010). Administrative boundary files were obtained from the Global Administrative Areas (GADM) Database (29). Monthly gridded (0.5°×0.5°) average air temperature (26) and total precipitation (27) data from 1993 to 2010 were obtained from the National Oceanic and Atmospheric Administration Earth System Research Laboratory Physical Sciences Division for latitudes between 0° and 26.8° and longitudes between 96.2° and 127.7°. Gridded (2.5 × 2.5 km) data on population density (per kilometer²) for the year 2000 were obtained from the Center for International Earth Science Information Network at Columbia University (28). Finally, we obtained monthly data on the Oceanic Niño Index (ONI) from the National Oceanic and Atmospheric Administration Center for Weather and Climate Prediction (www.cpc.ncep.noaa.gov/data/indices). The ONI indicates the strength of the El Niño climate variation based on sea surface temperature in the Niño 3.4 region.

1.2. Data Exclusions.

Province names were standardized based on the administrative level one or two divisions in the GADM Database as of 2013 (29). Provinces that could not be located were excluded (e.g., because of a split after the last update of the boundary definition files in the GADM Database). The study period was determined separately for each province. No data before 1993 were included because of lack of data for most countries. For each province, the first and last years of a study period were defined by the first and last years with at least six monthly observations. Provinces were excluded if 12 or more consecutive monthly observations were missing, two-thirds or more of all observations were missing, or if the study period was less than 4 y (2.5 cycles of the smallest multiannual period of 19 mo). These in- and exclusions resulted in a total of 273 of 289 provinces included (94%).

1.3. Data Management.

For multiple countries, the numbers and names of provinces changed over time because of merging or splitting events. We used province boundaries as in the GADM Database as of 2013; we merged data for provinces that merged and only used data after a splitting event for provinces that had split. Geographical distances between all province pairs were computed using latitude/longitude of province centroids according to the World Geodetic System 1984 Revision Coordinate System. We computed the average air temperature and precipitation per latitude for each month. We centered and reduced the time series for each latitude to z scores around the latitude mean. We also extracted the average value of air temperature, precipitation, and population density for province polygons from the gridded data. Population estimates were not available for the entire study period of each province. We applied country-specific, nationwide annual population growth rates obtained from the US Census Bureau International Database (25) to estimate the provincial-level population for missing years. The growth rate rc,y in country c and year y was computed asrc,y=Nc,yNc,y−1,[S1]where Nc,y is the population size in country c in year y. The population size in province p of country c in year y was then computed fromNp,c,y=rc,y×Np,c,y−1.[S2]

We used dengue case data and population estimates to compute monthly IRs per 100,000 people for each province.

1.4. Data Transformation.

IRs were log-transformed aslog10(105×I+1N),[S3]where I is the incidence (i.e., number of cases), and N is the population size. The log-transformed IRs were then detrended by subtracting fitted values of a linear model from the observed values. The fitted linear model readsIRp,m^=C+βp×m,[S4]where β_p is the linear regression coefficient, and m is the study month. The detrended IRs for each province (p) and month (m) are thusIRp,m˜=IRp,m−IRp,m^,[S5]where m refers to the month and IR_p,m refers to the reported IRs for each province and month. After detrending, missing IRs were imputed separately for each province by randomly selecting an observation for the same month but from a different year. Finally, IRs were centered (discounted the mean) and reduced (divided by the SD) to z scores to increase cross-province comparability.

1.5. Reconstructing Annual and Multiannual Cycles.

We used wavelet analysis to isolate cycles with multiannual and annual periodicities from the reported IRs. Wavelet analysis is well-suited to characterize epidemiological time series containing cyclical variability with periodicities that change over time (nonstationarity). Wavelet analysis has been used previously to study a wide range of infectious diseases, such as measles, pertussis, and dengue (2, 5, 8, 9, 11, 13). Detailed methods for wavelet analysis have been described previously (30, 31). We computed Morlet wave transforms for each province using a nondimensional frequency ω0 = 6 and a periodicity step size δj of 0.25 on a linear scale. The Morlet wavelet enables a high resolution of the periodicity scale and is complex, which makes it possible to extract phase angles to represent epidemic timing (7, 30). The Morlet wavelet has been widely used previously for the study of infectious disease IRs (2, 8, 9, 13). The value for ω0 of 6 has been used previously (8, 30), and a δj of 0.25 gives a detailed resolution of the periodicity scale (1-mo intervals). We conducted a sensitivity analysis to assess the effect of alternative values for ω0 and δj on the wavelet transforms and resulting synchrony between province pairs. Only extreme values of ω0 reduced synchrony, and changes in δj had no noticeable effect (Fig. S6).

For each scale s and time interval δt, the continuous wavelet transform of a time series x_k is defined as in the work by Torrence and Campo (30):Wn(s)=∑k=0N−1xkψ*(k−n)δts[S6]where ψ represents the Morlet wavelet function and * the complex conjugate. n represents the time index, ranging from zero to the total number of time points N. We computed wavelet transforms for periodicities (or scales) ranging from 2 to 60 mo. This range contains all periodic cycles previously observed for dengue IRs. Statistical significance of wavelet transforms was tested by comparing the wave signal with a red noise background signal (30, 31).

For each province, we computed the average of wavelet transforms for each periodic intervals over all time points n (global wavelet spectrum) according to the work by Torrence and Campo (30):W¯2(s)=1N∑n=0N|Wn(s)|2,[S7]where N represents the number of observations per province. For each province, we only explored periodicity s below the maximum periodicity smax,p, to reduce edge effects. We determined the maximum periodicity for each province p, s_max,p, as the number of monthly observations divided by 2.5.

Similarly, for each province, we also computed the average of statistically significant wavelet transforms per time interval δt across periodicities in the annual (6–18 mo) and multiannual (19–60 mo) band using the scale-averaged wavelet power according to Torrence and Campo (30):Wn¯2=δjδtCδ∑j=j1j2|Wn(sj)|2sj,[S8]where s represents the periodicity included, δj represents the scale interval size, and C_δ is a constant for the Morlet wavelet (0.776). We only explored periodicities below the maximum smax,p supported by the length of province time series. We used j ranging from 6 to 18 mo for the average for annual cycles and from 19 to 60 mo for the average for multiannual cycles.

We reconstructed annual epidemic cycles defined as the average of 6–18 mo periodic cycles and multiannual cycles defined as the average of 19–60 mo periodic cycles. We reconstructed these cycles using the filter described in the work by Torrence and Campo (30):x′n=δjδtCδΨ0∑j=j1j2ℜ{Wn(sj)}sj,[S9]where δj was 0.25, Cδ was 0.776 as previously used (30), and values for j1 and j2 were 6 and 18 mo, respectively, to reconstruct annual cycles and 19 and 60 mo, respectively, for multiannual cycles. We only reconstructed waves for provinces that had statistically significant wavelet transforms for at least 50 time points within the annual or multiannual periodicity bands.

1.6. Estimating Synchrony.

We defined synchrony between provinces as the pairwise Pearson correlation coefficient, taking into account both timing and amplitude of the signals. We computed synchrony for both annual and multiannual cycles. To assess the association between synchrony and geographical distance between provinces, we used a previously described nonparametric covariance function that estimates the underlying correlation without assuming any particular shape (2, 33). It uses a smoothing spline (bandwidth of 300 data points) to estimate the curve and bootstrapping (n = 1,000) to estimate confidence limits:ρ˜(δ)=∑i=1n∑i=1nK(δijh)(ρ^ij)∑i=1n∑i=1nK(δijh),[S10]where K is a kernel function, h is the bandwidth, and ρij^ is the autocorrelation matrix.

We measured changes in synchrony over time by computing synchrony separately for parts of the time series within moving and overlapping 5-y time windows. For each time window, we computed the average synchrony per province weighted by the number of province pairs included. We did this for annual and multiannual cycles and the unfiltered IRs (Fig. 3). We used a linear model to estimate the association between synchrony and the pairwise average log₁₀ population density, precipitation and temperature for province pairs, adjusted for geographical distance between pairs:ρp^=C+β1dp+β2t¯p+β3κ¯p+β4log10(λ)¯p,[S11]where ρ is the synchrony of multiannual or annual cycles or the unfiltered IRs between unique province pair p, d is geographical distance in kilometers, t is temperature in degrees Celsius, κ is precipitation in millimeters, and λ is population density per kilometer².

1.7. Wavelet Coherency.

As an alternative metric for synchrony over time and for comparison with the wave correlation analysis, we measured wavelet coherency between province pairs. As previously described, wavelet coherency uses wave transforms of two time series to indicate their localized phase relationship in a time–frequency spectrum (7, 32). Wavelet coherency ranges from zero to one, and high wavelet coherency requires that statistically significant cycles of a particular periodicity are detected in both time series and that these cycles are phase-dependent (positively or negatively):Rn2(s)=|S(s−1WnXY(s))|2S(s−1|WnX(s)|2)×S(s−1|WnY(s)|2),[S12]where S is a smoothing operator. Statistical significance of wavelet coherency was tested using Monte Carlo methods (n = 600) (32). For the wavelet transform in the wave coherency function, we used the same parameters as specified for the wave correlation analysis (i.e., a value for ω0 = 6 and a δj resulting in a periodicity step size of 1 mo along a linear scale). We only measured wavelet coherency for periodicities below the maximum for each province smax,p.

We studied changes over time in wavelet coherency for each province by computing the percentages of other provinces that were statistically significantly coherent with this province at each time point (Fig. S1 C and D). In addition, we used wavelet coherency to estimate the association between the ONI and the multiannual dengue cycle for each province. For each province, we averaged all statistically significant values of this wavelet coherency matrix within the multiannual periodicity band (19–60 mo).

1.8. Phase Differences.

As previously described, we computed phase angles for annual and multiannual cycles as indicators of epidemic timing (2, 30, 31). For each province and for annual and multiannual periodicities separately,ϕn(s)=tan−1(ℑ{WnXY(s)}ℜ{WnXY(s)}).[S13]

We computed the average phase angle-ϕn for annual and multiannual cycles as the circular average of periodicity-specific phase angle-ϕn(s) across the annual or multiannual periodicity bandwidth (6–18 and 19–60 mo, respectively) as described previously (32). Because we only reconstructed annual and multiannual cycles for provinces with statistically significant cycles in these bandwidths, we only computed phase angles for provinces with statistically significant cycles.

We used these phase angles to compute the phase angle difference θ between provinces for annual and multiannual cycles. The phase angle difference θ was constrained between −π and π. We expressed θ in months and assumed a cycle length of 12 mo for annual and 39 mo for multiannual cycles. We used a linear model with a linear spline to measure the relation between θ and geographical distance. We used unique absolute values of θ to only include unique province pairs and included a linear spline to detect the distance after which the relation between θ and distance, d_p, would become nonsignificant:|θp^|=C+β1dp+βp,s=500S((dp=s)×(dp−s)),[S14]where we varied s with increments of 100 km from a minimum of 500 km up to a maximum S of 2,500 km and selected the best fitting model based on Akake’s Information Criterion. The association between θ and distance was only statistically significant before the spline point s. We used the spline point of the best fitting model as the distance radius within which we detected local traveling waves. For multiannual cycles, the s of the best fitting model was at 1,000 km, and for annual cycles, it was at 1,500 km.

1.9. Traveling Waves.

We defined a local traveling wave for a province as a statistically significant positive linear association between θ and geographical distances between that province and the other provinces. A province could have a negative or positive θ with other provinces. A negative θ (lag time) indicated that a province was timed later compared to the other province (lagging behind), and a positive θ (lead time) indicated that the province was timed earlier compared to the other province (leading ahead). Outgoing traveling waves can emerge from provinces that are timed ahead of others. Incoming traveling waves can occur for provinces that are timed behind others. For provinces that were lagging behind (negative θ), we defined an incoming traveling wave as a decreasing θ with decreasing distance. For provinces that were leading ahead (positive θ), we defined an outgoing traveling wave as an increasing θ with increasing distance. For each province and separately for annual and multiannual cycles, we used a linear model to detect traveling waves:θp,q={C−βpdp,qθp,q<0; incoming waves for province pC+βpdp,qθp,q>0; outgoing waves for province p,[S15]where p is a province, and q indicates all other provinces. We only included values of d below 1,000 (multiannual cycles) or 1,500 km (annual cycles). We inversed the sign of distance dp,q for θp,q values below zero for more intuitive displays of incoming waves (Fig. 5). We tested statistical significance of βp for incoming and outgoing waves separately for each province. We used Bonferroni-corrected significance levels of 0.0002 per province for a combined significance level of 0.05 for each group of tests: incoming waves of multiannual cycles, outgoing waves of multiannual cycles, incoming waves of annual cycles, and outgoing waves of annual cycles. We measured the association of having multiannual or annual traveling waves (binary variable) for provinces with their average population size, temperature, and precipitation using a logistic regression model. For incoming and outgoing waves of annual and multiannual cycles separately:τp=C+β1sp+β2tp+β3κp,[S16]where τp is a binary variable specifying the presence or not of a traveling wave for a province, sp is the average population size of a province across study years, tp is the average monthly temperature in degrees Celsius for a province, and κp is the average monthly precipitation in millimeters. Given that we used this model separately for four different waves, we used a Bonferroni-corrected significance level for a combined level of 0.05.

1.10. Sensitivity Analysis.

We conducted a sensitivity analysis on the effect of the value for the wavelet parameters ω0 and δj on synchrony. We computed wavelet transforms as described above for annual and multiannual cycles while varying the value of ω0 from 2 to 10. For each value of ω0, we computed the average synchrony between each province and all of the other provinces weighted by the number of pairwise observations (Fig. S6 A and B). We did the same for values of δj ranging from 0.05 to 0.5 (Fig. S6 C and D).

1.11. Computing Environment.

All analyses were conducted using the R System, version 3.2.1. We used the package dplR for the wavelet analysis with a modified Morlet function with a linear periodicity scale. We used the package biwavelet for wavelet coherency with modified wtc and wt functions with a linear periodicity scale. We used the function Sncf for the nonparametric analysis of the correlation between synchrony and geographical distance.