Prediction of Ross River virus incidence in Queensland, Australia: building and comparing models

Introduction

Ross River virus (RRV) is a mosquito-transmitted Alphavirus that causes arthritis, rash, and constitutional symptoms of fever, fatigue, and myalgia (Liu, Tharmarajah & Taylor, 2017; Harley, Sleigh & Ritchie, 2001). It is the most common human arboviral infection in Australia. Between 1993 and 2020, almost half of the nationally notified cases (49.0%, 63,880/130,271) occurred in Queensland with a total of 43,699 cases reported from 2001 to 2020 (Qian et al., 2021; Australian Department of Health, 2021). Outbreaks of RRV have occurred regularly in Queensland over the past two decades (Qian et al., 2021). Timely prediction of outbreaks based on epidemiological models is necessary for disease prevention and control. Transmission of RRV is influenced by mosquito abundance, reservoir host populations, climate and weather, geographical exposures, and socio-economic indices (Qian et al., 2020; Murphy et al., 2020; Tong et al., 2008). Weather patterns may influence RRV transmission differently at different times and in different parts of Queensland (Qian et al., 2020). The dominant vector and reservoir host species in coastal and inland areas differ (e.g., Culex annulirostris breeds in freshwater habitats, Aedes vigilax breeds in mangroves and salt marshes, and Aedes notoscriptus breeds in artificial containers) (Stephenson et al., 2018; Harley et al., 2000). In Queensland, temperature, rainfall, relative humidity, and high tide are possible predictors of RRV transmission in coastal regions, while temperature, or relative humidity and rainfall are associated with inland RRV activity (Tong & Hu, 2002; Bi & Parton, 2003). The number of exposures and the complex relationships between these exposures (Ng et al., 2014; Gatton, Kay & Ryan, 2005) make RRV modelling challenging.

Many studies have applied epidemiological models to predict RRV cases, incidence, or outbreaks, with exposures including mosquito abundance, rainfall, temperature, tidal height, humidity, and river flow (Qian et al., 2020). Half of these studies employed generalised linear models, which are straightforward and simple, and hence commonly used in predicting infectious diseases. Other approaches such as time-series models, spatial and temporal models, hurdle models, and Generalised Additive Models (GAMs) were also used. A few studies have published data related to model performance. Performance accuracies varied from 63.0% to 100.0%, sensitivities and specificities varied from 0.0 to 1.0, overall agreements varied from 75.8% to 88.5%, and true positives varied from 0.0% to 100.0% (Hu et al., 2006; Woodruff et al., 2002; Woodruff et al., 2006; Jacups, Whelan & Harley, 2011; Pelecanos, Ryan & Gatton, 2010; Gatton, Kay & Ryan, 2005; Koolhof, Bettiol & Carver, 2017; Ng et al., 2014). Most weather and climatic exposures used in previous studies (e.g., rainfall, temperature, vapour pressure, and evaporation) have been found to have positive effects on disease transmission, which means that increases of these exposure values are related to an increase in disease incidence. High relative humidity generally reduces RRV incidence. Other factors include mosquito abundance, tidal height, and river flow (Qian et al., 2020). Stepwise variable selection was often used in developing RRV predictive models (Koolhof et al., 2021). Regional analysis was suggested but rarely used in previous studies (Ng et al., 2014). Studies assessing predictive model comparison are lacking.

Selecting important predictors can improve our understanding of disease transmission. A model building process with detailed considerations of lagged predictors, variable selection, and validation, and a head-to-head comparison of the models will allow assessment of suitable model frameworks for future prediction of RRV. We applied 10 models and compared their model fit and performance in predicting RRV incidence. The Poisson and negative binomial generalised linear models were implemented to fit the distribution of RRV cases. Zero-inflated Poisson and Zero-inflated negative binomial models were developed due to the excess zeros caused by relatively small spatial and temporal resolution. Generalised additive models with and without zero-inflated Poisson or zero-inflated negative binomial models were also used. A variable selection strategy to screen the most important variables for prediction was developed and regional analysis was used.

By building and comparing these models, we aimed to (1) identify the most robust approach for predicting RRV incidence with good model fit among generalised linear models, zero-inflated models, and GAMs; and (2) develop a model building process that includes an effective variable selection strategy to build models with both excellent predictive performance and clear epidemiological logic in different regions. The results and the strategy provide important guidance for mosquito-borne disease modelling.

Materials & Methods

Variable selection

Before building the model, we developed a variable selection process based on the purposeful selection of covariates strategy proposed by Hosmer, Lemeshow & May (1999) to obtain a parsimonious variable set which fitted the hypotheses of the models and had strong epidemiological and statistical correlations with RRV activity. In this study, exposures were the potential influencing factors of RRV transmission regardless of lags. The impact of exposures on RRV transmission may lag by several weeks. We referred to the exposures at specific lags, including lag zero which represents the current value of an exposure, as predictors (or lagged variables). For example, temperature was an exposure, and maximum weekly temperature at lag 5, which was at 5 weeks earlier, was a specific model predictor.

Optimal weekly lags of temperature, range of temperature, relative humidity, rainfall, vapour pressure, evaporation, evapotranspiration, SOI, El Niño events, La Niña events, NDVI, bushfires and flood events were selected based on the maximum positive and minimum negative values of the cross-correlation function, which showed the correlation between RRV cases and exposures at different lags (Curriero, Shone & Glass, 2005). Due to the possible time periods of mosquito breeding, transmission to reservoir hosts, the incubation period in humans and lags for Ross River virus reporting, we considered that lags greater than one year would not plausibly have meaning for associations of exposures and RRV incidence.

Then, a univariable analysis was implemented to exclude predictors with p value greater than 0.1 or standardised regression coefficients less than 0.1. These predictors were considered as having a weak impact on RRV incidence. Spearman correlation evaluates the monotonic relationship between two variables. Predictors with absolute value of Spearman correlations greater than 0.9 were highly correlated with each other (Akoglu, 2018). One of the two correlated predictors was considered as redundant and was excluded. The Variance Inflation Factor (VIF) was used to calculate the inflation of coefficient variance and address the magnitude of multicollinearity in models. Predictors with high VIF values were removed to reduce redundancy and to avoid possible problem caused by multicollinearity (e.g., distorted statistical significance or predictor coefficients which could misinterpret the effect of predictors on RRV transmission). A repeated backward stepwise screen and reassessment process was used to screen predictors improving model fit in predicting RRV incidence with tolerable multicollinearity. Predictors with VIF greater than five were removed (Menard, 2002). The predictor with lowest Bayesian Information Criterion (BIC) increment was removed in each iteration until the total BIC increment reached 0.5%. The predictors were then considered in the model one at a time, and those predictors with VIF values less than ten and BIC decrements greater than 0.1% were included. We continued this process until no predictor was excluded or included, two consecutive iterations returned the same variable set, or it reached 11 iterations. Initially, this process was applied to the climatic variable set, and to the geographical and socio-economic variable sets separately, then, we applied to all variables selected from the two sets. Finally, predictors with VIF greater than five were excluded. The remaining predictors were selected for model building. Further details of this process are shown in Fig. 1.

Results

Regions of Queensland have different ecological and climatic patterns, vector species, and host populations, which influence RRV transmission. As shown in Fig. 2, 14 SA3 areas have a hot humid summer (abbreviated as “Hot” region); 61 have a warm humid summer, or warm summer and cool winter (“Warm” region); and the remaining seven areas have a hot dry summer, and mild or cold winter (“Dry” region; Fig. 2).

A total of 43,699 cases were reported in Queensland, Australia from 2001 to 2020. Among these cases, 13,422 were recorded in the Hot region, 25,949 were recorded in the Warm region, and 4,325 were reported in the Dry region. The mean incidence rate was five per 10,000 people per annum during the study period. Some SA3 areas on the coast and in Western Queensland reported high cumulative incidence rates across the 20-year study period (Fig. 2). The distribution of weekly RRV notifications in Queensland is shown in Fig. S1.

Ross River virus cases had strong seasonal variations with annual peaks in late summer and autumn (Fig. 3). Since 2006, RRV cases have increasingly been reported throughout the year. The highest notification number and highest incidence rate were recorded in 2015, while the lowest were reported in 2002. The weekly average values for maximum temperature, rainfall, and relative humidity at minimum temperatures are displayed in Fig. 3.

For each training dataset in each region, predictors (exposures at different lags) were selected to fit the standard negative binomial generalised linear model and then used in the training and validation of all the other models using the same dataset of this region. The coefficients and 95% confidence intervals of chosen predictors in the NB model using the third training dataset for three regions are listed in Table 1. For each training dataset, nine to thirteen predictors were selected. Disparity in the predictor sets selected for the three regions indicates that regional analysis is necessary for detecting important predictors influencing RRV transmission under different environmental and weather conditions. The performance of models in preliminary analyses, which used the same model in selecting predictors and building models, are provided in Table S1. The results show that the NB model is appropriate for variable selection.

Evapotranspiration, NDVI, and vapour pressure with appropriate lags are selected for modelling RRV incidence in all regions, and a lagged El Niño variable was selected in three regions. Different lags for a single exposure may have both positive and negative effects on RRV incidence. In the Dry region, maximum temperature, vapour pressure and the index of education and occupation had positive effects on disease transmission, while evapotranspiration, El Niño events, proportion of agriculture land and NDVI were negatively related to RRV incidence. In the Hot region, rainfall, vapour pressure, evaporation, evapotranspiration, elevation and NDVI were positively correlated with RRV incidence. Temperature and types of land use were negatively associated with disease activity in the Hot region. In the Warm region, high values of weather exposures (relative humidity, vapour pressure and evapotranspiration) accelerated RRV transmission, while El Niño events and high values of SOI and NDVI decreased RRV incidence. A high proportion of nature land and scores indicating less access to service centres also indicated areas with high disease activity.

The average AIC, BIC, and HQIC of the training datasets were calculated for each region and the whole of Queensland. The NB-based models had lower BICs than the Poisson-based counterparts, which indicated that over-dispersion exists in all datasets (Fig. 4). A lower variation in BICs indicated NB-based models were more stable than Poisson-based models. All NB-based models had similar BIC values regardless of model complexity across all regions. So, in general, NB-based models are robust, stable, and performed well in the epidemiological modelling of RRV surveillance data. The AIC and HQIC of all models are in Figs. S2–S3.

Within the Poisson-based models, Zero-Inflated models (ZI models) fitted the data better than the others, suggesting that the inclusion of zero-inflation in the model is needed if data dispersion is not considered. Considering the linear and non-linear models have similar BIC values, there is no evidence that non-linear models fit the data substantially better. Over the three model components considered here, the over-dispersion effect is most important, followed by zero inflation, with non-linearity being the least important in fitting these data.

The MSE, RMSE and MAE of the validation data in all regions are provided in Fig. S7 and Table S2 to show the unpenalized relative residuals of the models. These relative residuals were lower for the generalised linear models and ZI models. The GAM models with zero-inflation tended to have the greatest relative residuals in the validation sets.

The predicted values of the NB and the NBGAM models, which were stable and had lower BICs in general, together with real RRV cases in Queensland are displayed in Fig. 5. The results for the three regions are provided in Figs. S4–S6. These two methods were considered as the representative linear and non-linear models for this study. Both models fitted the seasonal patterns and timing of peaks well but were unable to fit the size of the peaks, especially when the number of cases were high. Both models predicted one to two years of data comparatively well but failed to predict the latter two years. The non-linear models were less stable than linear models in prediction, especially in the final year, and over-fitting may have occurred.

Discussion

We modelled RRV incidence using climatic, geographical, and socio-economic exposures in different regions of Queensland. Using a substantial amount of data in a wide geographical area across 20 years with numerous exposures provided a reliable basis for model comparison. A purposeful selection strategy was developed to identify important predictors among numerous lagged variables. Ten models including generalised linear models (Poisson and negative binomial models), zero-inflated models (ZIP, ZINB, ZIPe, and ZINBe), and GAM models (PGAM, NBGAM, ZIPGAM, ZINBGAM) were implemented, where ZI models were applied to RRV data for the first time. We first compared model performance and evaluated the effects of over-dispersion, zero-inflation, and non-linearity of the epidemiological surveillance data on model fit. The standard negative binomial generalised linear model fit the data well, and it is simple, robust, and efficient. Over-dispersion is the most important effect to be considered when analysing this type of data.

Most researchers applied a stepwise method for variable selection in developing RRV prediction models (Koolhof et al., 2021; Gatton, Kay & Ryan, 2005; Bi et al., 2009; Ng et al., 2014). Stepwise selection, an automated approach using either predictor significance or overall fit, does not account for a complex setting involving the modelling using correlated exposures at various lags. Limitations of stepwise selection include parameter estimation bias, inconsistencies among model selection algorithms, and an inappropriate focus or reliance on a single best model (Whittingham et al., 2006). Appropriate selection methods for climatic, environmental, and socio-economic predictors with lagged effects, taking account of complex correlations among predictors, have not been well studied. To choose important predictors, our strategy screens exposures at different lags based on both univariable and multivariable (adjusted) associations (the repeated backward stepwise screen and reassessment process). The significance and standardised coefficients of the predictors, correlations between predictors, and model fit are considered to select predictors. The stepwise process followed by a reassessment could avoid removing predictors based on a single criterion. The variable selection approach was able to identify the most important predictors by considering statistical significance, standardised regression coefficients, multicollinearity, overall fit and correlation between predictors in a comprehensive way. It was found to be effective and able to generate consistent and reasonable numbers of predictors across all datasets in this study, and we believe it has great potential in many epidemiological predictive studies including surveillance data arising from ecological designs. Different predictors were selected for each region and their effects on RRV activity varied across regions. Regional analysis is beneficial for selecting important predictors relevant to the transmission cycle in each region. However, adequate sample size is required to ensure a robust model fit.

Several exposures are proved to play relatively important roles in predicting RRV incidences, such as evapotranspiration, Normalized Difference Vegetation Index (NDVI), vapour pressure, and El Niño events. We found vapour pressure has a positive association with RRV transmission in all regions regardless of lags, as demonstrated in other studies (Cutcher et al., 2017). Evapotranspiration has a positive effect on RRV incidence in the Hot and Warm regions and a negative effect in the Dry region. As mentioned by Koolhof et al. (2020), evapotranspiration may be a stronger predictor than other climatic factors including rainfall and temperature. The impact of El Niño events on RRV incidence has been less studied. We found a negative association between El Niño events and RRV incidence in the Warm and Dry regions. The decreased rainfall and higher temperature during El Niño periods may lower RRV incidence (Australian Government, 2020). Vegetation cover favours mosquito breeding and provide habitats for reservoir hosts. However, in forests with dense vegetation, mosquitoes may be less likely to survive and breed because of the low temperatures in the shade. The average NDVI value is high after rainy days in late autumn and early winter in Queensland when the disease incidence is low (Australian Government, 2019b). These might partly explain our findings that high NDVI at three months’ lag in Dry and Warm areas was associated with higher RRV activity, while high NDVI in the Hot regions around 10 months earlier may lead to vegetation growth in forests and cause decreased RRV incidence (Ha et al., 2021). A high air temperature may prevent mosquito breeding, thus may decrease disease transmission in the Hot region, whereas high temperatures in the Dry region may help mosquitoes to survive through a cold winter and subsequently increase RRV incidence (Jacups, Whelan & Harley, 2011; Werner et al., 2012). Relative humidity lagged at 1–2 months may accelerate mosquito breeding and vegetation growth, thus indirectly contribute to disease transmission (Ng et al., 2014; Tong & Hu, 2002). The selected predictors, such as vapour pressure, NDVI, and El Niño events, can provide information about conditions that could accelerate mosquito breeding or indicate high-risk areas that require disease prevention. Targeted mosquito surveillance and control programs (e.g., mosquito monitoring, reduction of breeding sites, and control of larval and adult mosquitoes) can reduce transmission of the disease. Interaction between exposures was not considered here due to the large number of exposures; however, there is potential to develop models that capture all three regions simultaneously, with interaction terms used to capture differences between regions.

Modelling RRV is challenging because of the complicated transmission cycle, the characteristics of aggregated notification data, the wide geographical area of Queensland, and the 20-years of data. To avoid the limitations of a single method, we considered the effect of over-dispersion, zero-inflation, and non-linearity, and added corresponding components into the models. Over-dispersion of the data may be due to excess zeros, due to the distribution of RRV cases in SA3 areas, and models that incorporate consideration of over-dispersion were found to fit the data better. Indeed, we found that the standard negative binomial generalised linear model fit the data well suggesting that it’s ability to model over-dispersion may have been sufficient to solve to problems of excess zeros, which often manifests as over-dispersion. Non-linearity can be a good addition to NB models, but the GAM models did little to improve model fit without considering dispersion or zero-inflation of the data. Generalised linear models and ZI models performed better than GAM models, with or without zero-inflation, based on relative residuals in the validation sets. Our aim is to build and obtain models with excellent predictive performance and clear epidemiological logic in different regions, which is to balance a good model fit in the training set with selected explainable exposures and a good model predictive performance in the validation set. Generally, among the models implemented in this study, NB-based models were the best choice and are recommended for similar epidemiological studies that consider vector-borne diseases at different temporal and spatial resolutions.

Two noteworthy peaks of RRV notifications were reported in 2015 and in 2020. The 2015 outbreak were possibly due to a combination of ecological factors including above-average rainfall and consequent high mosquito abundance (Jansen et al., 2019). The high cases reported in 2020 may be related to altered RRV transmission risk and reporting. Changed behaviour due to COVID-19 may have resulted in more time gardening or exercising outdoors, while awareness of symptoms of fever and fatigue may have affected healthcare seeking (Webb, 2020; Jansen et al., 2021). The 2015 outbreak was included in our training set and the 2020 peak was in the validation set, allowing a test of model performance when RRV incidence was extremely high. The predictors selected could explain the common trends of RRV incidence but were unable to predict the rise of disease transmission in 2015 and in 2020. Historical RRV cases, modelled with an autoregressive term, may help to detect incidence within a certain range and over a short period. Some important exposures, such as mosquito abundance and host populations, have been suggested as greatly improving model fit (Ng et al., 2014), but are not included in this study as they are not available across Queensland for our study period. Adding autoregressive terms and using different formats of variables (e.g., using the number of rainy days instead of average rainfall in a week) may improve model fit and predictive performance.

Data quality has an impact on modelling performance. The definition of RRV was revised in 2013 and 2016 to remove possible false-positive diagnostic test results (Selvey et al., 2014; Knope et al., 2019). We are unable to evaluate the impact of these revisions on predicting RRV incidence. As a passive surveillance system, RRV surveillance may be biased for the following reasons: underestimation caused by asymptomatic infections, false-positive diagnostic tests (Selvey et al., 2014), greater efforts consequent upon greater awareness by doctors to ascertain cases during outbreak years (Werner et al., 2012), and lags in reporting due to the incubation period (Queensland Government, 2017). Flooding and bushfires were coded as binary in each SA3 area and do not capture areas influenced by these events or the river height of related river catchments. The SEIFA and ARIA data from censuses in 2011 and 2016 had to be extrapolated to the entire study period. Data on land use and elevation were only available in one year, thus were considered as constant values in this study, with spatial estimates of elevation interpolated to areas with missing data. The weekly population data required linear interpolation and extrapolation when converted from yearly data. Data available with adequate quality were considered and data were processed with reliable procedures.

Potential methodological limitations might exist when we determined data units for analysis and chose models for variable selection. The weekly SA3 area data used in this study contain details but have many zeros. However, using data at a lower spatial or temporal resolution (e.g., monthly data) have fewer zeros but may lose information, especially for daily weather data. The standard negative binomial generalised linear model was used for selecting predictors for all models which may have given it an advantage in model fitting. However, using the same set of variables makes head-to-head model comparison possible. We chose here to select the same set of predictors for all models to enable a comparison of different models, rather than a comparison of the model building process, with our additional analyses suggesting that any advantage to the negative binomial model is relatively small.

Conclusions

This is the first head-to-head comparative modelling study of prediction of RRV incidence in Queensland. We developed a comprehensive strategy to purposefully select variables using multiple criteria followed by reassessments and to allow the choice of a stable and valuable variable set for model fit. The comparison among different types of models provides an example of determining the key element affecting model fit. We found the standard negative binomial generalised linear model to be the best choice in predicting RRV incidence and the statistical over-dispersion effect possibly caused by excess zeros is a principal issue to be considered. We advocate a model selection algorithm, such as what we developed, that accounts for some issues arising while modelling infectious disease incidence including various exposures, lags, data distribution, multicollinearity, and model fit. The model building process introduced is promising for epidemiological surveillance studies of mosquito-borne disease with temporal and spatial data.

Supplemental Information