Predicting dengue importation into Europe, using machine learning and model-agnostic methods

Dengue data

We obtained monthly data for imported cases of dengue in Europe, for 2010 – 2015, from the European Centre for Disease Prevention and Control (ECDC)²⁴. Here, we utilized confirmed dengue cases (as defined the European Union generic case definition for viral haemorrhagic fevers) with known travel history²⁵. A total of 21 European Union/European Economic Area (EU/EEA) countries reported data on imported dengue, from a total of 98 different source countries between 2010 and 2015 (inclusive of zero reporting). The monthly level case counts were aggregated by country of infection (as source country) and the reporting country in Europe (as destination country). We transform the absolute count data into a binary response variable, that indicates whether there was an imported case of dengue (1) or not (0) in a destination country, in a month.

Air passenger’s data

Comprehensive air passengers travel data for 2010 – 2015, was obtained from the International Air Transport Association (IATA)²⁶. The data included true origin, connecting points and final-destination airports for all routes in the world and their corresponding passengers’ volume. The data contains over 11,996 airports in 229 different countries and their territorial dependencies. The passengers’ travel volume for each route worldwide was available at the country level and at a monthly timescale. This data was used to construct a monthly passenger flow from all countries worldwide with a final destination in Europe (accounting for all connecting flights) between 2010 and 2015. The data also included the passengers’ flow between European countries.

Connectivity indices between a source and destination country

Drawing on the underlining concept of spatial interaction modelling, that inflow between two locations is a function of the attributes of the source and destination and their corresponding interaction²⁷. Connectivity indices between a source country and a destination country in Europe were previously developed¹³, using different factors that potentially mediate dengue importation risk. The indices were decomposed into components describing the source strength’ (the risk of dengue infection) and the transport or importation potential (the connection between a source country and a potential destination country in Europe). Source strength for all indices was modelled to represent the endemicity of dengue in a source country. While transport and importation potential were modelled to characterize seasonal dengue activity, incidence rates, geographical proximity, epidemic vulnerability, air passenger volume, population size and wealth of a source country as mediating risk factors. The connectivity indices and their descriptions are listed in Table 1.

Centrality measures of the air transport network

Using the monthly air passenger’s data, we constructed a weighted directed network. The network for each month was denoted by G_m = (V, E), where V_G is a set containing all the nodes (or vertices), while E_G contains all the edges, with m indicating the month (m = 1, 2, 3… 72, covering the years of 2010 – 2015). Nodes represented all countries worldwide, while edges represent the flow of passengers from a source country to a destination country in Europe. Four different centrality measures were used to analyses the network and quantify the capacity of a source node to influence transportation of dengue or a destination node’s propensity to receive an imported case of dengue, by virtue of their connection topology within the network. Centrality measures and their descriptions are listed in Table 1.

Feature (variable) engineering

One fundamental step in building machine learning models is the process of feature engineering, i.e. using domain-specific knowledge to create new features (i.e. variables) or transform and encode existing original data into a more informative format^28,29. For this analysis, we created an additional variable based on our a priori knowledge of the data. The reporting rate for dengue data was heterogeneous across the EU/EEA member countries, as some countries were not consistent in monitoring and reporting to the ECDC (we considered “zero reporting”, which designate that no imported case was recorded during the reporting period, but a report was submitted to the ECDC). Hence to account for this variability in reporting rates, we created a time effect variable – coded as a set of dummy variables, 1 for a given month and 0 for all other months. This variable not only controls for reporting rates, but for other potential time-specific effects that might be restricted to a given period. Time events that might increase or decrease passengers’ traffic to a specific country, and in turn affect dengue importation (e.g. the introduction and/or discontinuation of an airline carrier or route).

Data pre-processing and splitting

The dataset used to build our machine learning models consists of the connectivity indices and centrality measures of the air transport network. The single unit of analysis is a source-destination country pair, at a monthly timescale, and a binary response variable coded to indicate an imported dengue case (1) or not (0).

Before model training, we performed the following pre-processing analyses to the full dataset. First, we examine the correlation between our predictor variables, by using Spearman’s correlation to rank the statistical dependencies. Several pairs of continuous variables displayed moderate-to-high pairwise correlations. Figure 1 shows the correlation matrix between the continuous predictors in our dataset. Most of the centrality measures (for source and destination countries) were highly correlated, example the betweenness and eigenvector centrality for source countries had a Spearman’s ρ = 0.99 This is not unusual, as on average, centrality measures are highly correlated in a network³⁵. The highly correlated pairs are practically redundant in conventional regression modelling and the heuristic approach to dealing with this is to exclude one. However, in our case, we did not manually exclude variables, as we focus on the predictive power of the entire bundle of variables as oppose the estimated coefficients of individual variables. Likewise, the suite of algorithms compared, each utilizes a combination of inbuilt feature selection and penalization functions to exclude redundant variables in their ensembling and mitigate the effect of multicollinearity.

• Download figure
• Open in new tab

Figure 1. Spearman correlation matrix of continuous variables.

Correlation is computed from the full dataset and coloured according to magnitude. Blue colours indicate strong positive correlations, red indicates strong negative correlations, and white implies no empirical relationship between the variables. Figure was generated using R programming language version 3.6.1³⁹.

Next, we randomly split the dataset into two sets, 70% into a training subset and 30% into a testing subset. This split was done based on the distribution of our outcome variable (i.e. binary response of an imported case of dengue (1) or not (0)), with sampling occurring within each category, thereby preserving the overall class distribution of our data. Our full dataset contained a total of 2, 055 unique country pairs (i.e. source -destination country pair), with a corresponding total of 147,960 monthly observations. An imported case was recorded in 1,937 observations, i.e. 1% of the total observations had an imported case. These observations were split into 1366 for the training subset and 573 for the testing subset, at a 70:30 ratio (a similar split was done for the observations with no imported cases). The training subset is used to build and tune the various machine learning models, while the test subset is used to evaluate the predictive performance of the models.

Due to the imbalance class distribution of the outcome variable we model this as a rare event, using a post hoc sampling approach to attenuate the effects of the imbalance during model training³⁶. The synthetic minority over-sampling technique (SMOTE), was used to subsample the training subset data to create a roughly equal distribution within the classes. SMOTE utilizes a hybrid of either up-sampling, to synthesize new data points in the minority class or down-sampling, to down-size the majority class³⁷. Our training dataset was balanced to a 3:4 ratio of an imported case to no imported case. The testing subset was maintained to reflect the original imbalance as a quality assurance of the predictive model performance. Lastly, before training, we apply a data transformation on all continuous variables in the dataset, by centring and scaling them. This transformation ensures that the variables have a zero mean and a common standard deviation of one, thereby improving normality and the numerical stability of the model calculations³⁸.

Model selection

Our training model is a classification-based model, to predict the probability that a case of dengue is imported into a country in Europe. There are several classification techniques (or classifiers) employed in machine learning models. The choice of the suite of algorithms we tested, was a trade-off between, meta-algorithm that fits our classification problem and those with built-in feature selection. Other algorithmic and systematic features that were considered include regularization (to handle the effects of multicollinearity), hyperparameter optimization (model tuning capabilities), and efficient computation time. To build our predictive model, we compare four widely used classifiers algorithms in machine learning, as listed below:

Partial least squares (pls) implements a supervised version of principal component analysis, using a dimension reduction technique. This technique first summarizes the original variables into a few new variables called principal components (PCs), as supervised by their relationship to the outcome variable. These components are then used to fit a linear regression model⁴⁰. For classification problems, the partial least squares discriminant analysis variant is fitted. This method has an embedded feature selection and regularization⁴¹.

Lasso and elastic-net regularized generalized linear models (glmnet) implement a logistic generalized linear model via penalized maximum likelihood. The addition of a penalty shrinks the coefficients of the less contributive variables toward zero (L2 ridge penalty) or absolute zero (L1-Lasso penalty)⁴². The glmnet implements a combination of both L1 & L2 penalties (otherwise called elastic net penalty), for its regularization and simultaneous feature selection.

Random forest (randomForest) is a bootstrap aggregated (or bagged) decision tree-based ensemble technique. The algorithm constructs multiple decision trees by repeat resampling of the training dataset and outputs the mode of the classes as a consensus prediction. The trees are created independently from a random vector distribution; hence each individual tree is heterogeneous with high variance and casts a unit vote for the most popular class⁴³. By averaging several decision trees, it intuitively avoids overfitting and performs an embedded feature selection.

Extreme gradient boosting (xgboost) implementation of a gradient boosted decision trees ensemble technique. The gradient boosting framework iteratively refines its model, to create a strong classifier by combining multiple weak classifiers in a stage-wise manner to minimize the loss function^44,45. The xgboost algorithm is a commonly preferred classifier, because it utilizes parallelization and distributed computation for implementation, thereby ensuring high efficiency in computation time and resources^45-47.

Model tuning and validation

Machine learning models can be prone to overfitting, to mitigate this we implemented a model building approach that encompasses model tuning and repeated evaluation during training. We use a methodological resampling technique of the training dataset, i.e. five repeats of 10-fold cross-validation (CV). The 10-fold CV randomly partitions the training dataset into 10 sets of roughly equal size, one set retained, and the others used to fit a model. The retained set is used to estimate model performance. The first set is then returned to the training set and the procedure iterated until each set has been used for validation. This whole process is repeated five times before results are aggregated and summarized. This procedure automatically chooses tuning parameters associated with optimal model performance.

Candidate models were evaluated using the following performance metrics: area under the receiving operating characteristic curve (AUC), sensitivity (true positive rate), and specificity (false positive rate). Our final candidate model was selected based on the receiving operating characteristic curve (ROC) threshold, which maximizes the trade-off between sensitivity and specificity⁴⁸. The ROC curve evaluates the class probabilities across a continuum of thresholds, with an arbitrary (algorithmically set) “optimal” cut point for determining what percentage of probability is accepted in classifying an imported case of dengue.

Model interpretability

We utilized the model-agnostic approach to provide interpretation of our optimal model. Model-agnostic methods work by extracting post-hoc explanations from an original machine learning model²³. This involves training an interpretable model on the predictions of the original model⁴⁹ and/or by changing the inputs of the original model and measure the changes in the prediction output⁵⁰. We employ the use of recent model-agnostic tools^51,52with both global and local scale interpretability functions. Global interpretation helps to understand the modelled relationship and distribution of the predicted target outcome (i.e. dengue importation) based on the input variables, while local interpretation zooms in, to help understand model predictions for a single instance (i.e. a single unit of observation or analysis).

We obtained global interpretations of our final candidate model through the following, variable importance, and partial dependence plots^49,53. Variable importance measures the contribution of each input variable, by calculating the increase in the model’s prediction error after permuting the variable⁵⁴. While the Partial dependence plots (PDP) are graphical renderings of the prediction function that helps visualize the relationship between the variables and predicted outcome^46,55,56. The relative importance of each variable is normalized to have a maximum value of 100, with higher scores indicating the most influential variable. We note that it may not be feasible to explore in detail the relationship of all variables in our model. Hence, we set an arbitrary cut-off on the variable importance measures at a value >50, to determine a subset of variables to focus on.

Local interpretation of our model was implemented via the use of local surrogate models, otherwise called-Local interpretable model-agnostic explanations (LIME)^23,50. The underlining assumption of LIME is that complex black box models are linear on a local scale, hence a simple (surrogate) model can be fitted for an individual observation that mimics the behaviour of the global model at this locality. The simple model and its variable weights are then used to explain the individual predictions locally. To demonstrate the LIME technique, we selected 10 single observations from our initial testing subset. These observations were sampled methodologically to include both classes (i.e. imported case [1] or not [0]) and representative of countries with a high and low frequency of dengue importation. We set the number of variables to best describe the predicted outcome, as the 5 most influential. The resulting weights for these variables are plotted to explain the local behavior of the model. The plots delineate if a variable increase or decreases the predicted probability of an imported case of dengue (detailed vignette for the LIME techniques can be found here^57,58.

Statistical software

All statistical analyses were performed with R Programming Language version 3.6.1³⁹. For uniformity in our model build, we utilized the classification and regression training (caret) R package, this is an interface to a vast amount of available machine learning algorithms³⁸. The package streamlines the process of building and validating predictive models by using a set of intuitive call functions. Supporting packages for specific functions includes: pls⁴¹, glmnet⁴², randomForest⁵⁹, xgboost⁶⁰, plyr⁶¹, doSNOW⁶², DMwR⁶³, pROC⁶⁴, pdp⁵⁶, iml⁵¹, lime⁵² and their various dependencies.

Model prediction performance

We compared the prediction performance of the different classifiers’ algorithms, in their ability to predict an imported case of dengue. Models were evaluated with the testing dataset, utilizing the area under the receiver operating characteristic curve (AUC) as the quantitative measure for performance comparisons. All four models performed comparably well, with AUC scores above 0.80 (Table 2). AUC score using pls was 0.88 (95% CI, 0.86 to 0.90); glmnet was 0.89 (95% CI, 0.87 to 0.91); randomForest was 0.97 (95% CI, 0.96 to 0.98); and xgboost was 0.97 (95% CI, 0.96 to 0.98). Performance metrics for each model are depicted in Table 2. Figure 2 also shows the ROC curve plots for the different models.

• Download figure
• Open in new tab

Figure 2. Comparison of the receiver operator characteristic (ROC) curve for the different models.

Curves characterize the trade-off between the sensitivity (true positive rate) and 1minus specificity (false positive rate). The y-axis = sensitivity and the x-axis = 1 minus specificity.

The AUC score indicates that predictions from the randomForest and xgboost models were better fitted to the dataset, outperforming the pls and glmnet models (with the pls being the least fitted). The randomForest and xgboost models had similar performance across the metrics with nearly negligible differences (Table 2). However, they had a distinction in their ROC curves effective threshold (Figure 3). The ROC threshold appropriately maximizes the trade-off between sensitivity and specificity. With the best threshold cut-off at 68% (i.e., only probabilities greater than 0.68 were classified as an imported case of dengue, Figure 3a) the xgboost model outperforms the randomForest, (cut-off at 0.64, Figure 3b) in a competitive comparison of prediction accuracy. With a true positive rate of 0.88 and a false positive rate of 0.12, the xgboost model was selected as the optimal model for our dataset. Hence, our final predictive model was able to predict 88% of dengue importation cases in our test dataset accurately (Figure 3a).

• Download figure
• Open in new tab

Figure 3. Comparison of the receiver operator characteristic (ROC) curves for extreme gradient boosting and random forest models.

The dot on both plots indicates the value corresponding to the “best” cut-off point threshold for each model that appropriately maximizes the trade-off between sensitivity and specificity. The numbers in parentheses are (specificity, sensitivity). Extreme gradient boosting (a) cut-off was at 68% (i.e., probabilities greater than 0.68 are classified as an imported case of dengue), delivering a specificity of 0.883, sensitivity of 0.880, while random forest (b) cut-off was at 64%.

Model interpretability

The best performing model (i.e. xgboost), was examined further for interpretation. Figure 4 illustrates the variable importance scores for the10 most influential variables, from our optimal (i.e. xgboost) model. With an arbitrary cut-off at >0.50, the subset of our ‘most important’ variables included the following: Source country’s dengue incidence rate; population; the number of arriving passengers; betweenness, closeness and degree centrality measures of the destination country. Figure 5 illustrates a visual representation of the relationship between this subset of variables and the predicted response while accounting for the average effect of the other predictors in the model. These plots demonstrate that the probability of an imported case of dengue increases on average for source countries with higher incidence rates, large population size and higher passenger traffic. Likewise, the probability increases for destination countries with higher betweenness, closeness and degree centrality measures.

• Download figure
• Open in new tab

Figure 4. Variable importance plots.

Top 10 most influential variables from the extreme gradient boosting model. The relative importance of each variable is normalized to have a maximum value of 100, with higher scores indicating the most influential variable.

• Download figure
• Open in new tab

Figure 5. Partial dependence plots for a sub-set of the most influential variables in the optimal model predicting the probability of an imported case of dengue.

The optimal model is the xgboost model. Sub-set variables represent variables with a variable importance ranking score >50. Y-axis is set on a probability scale since our model was a classification model; Blue rug marks at the inside bottom of plots indicate the min/max and deciles of the variable distribution. Top plots show the 3 most influential connectivity indices, while bottom plots show the 3 most influential network centrality measures.

Note: to guard against over-interpreting the partial responses, we added rug displays to the plots, i.e. tick marks indicating the minimum, maximum and deciles of the variable distribution. Our interpretation of the partial response was limited to the regions within the minimum and maximum.

In addition to providing global explanations of our optimal model, we provide a local explanation for individual predictions given for 10 single observations, i.e. a single unit of analyses, source– destination-month combination (Table 3). Figure 6 gives a visual representation of the first four single observations in local subset data, each plot shows the predicted probability of each observation, being an imported case of dengue. Likewise, it shows the five most influential variables that best explain the model’s prediction at the local region of the single observation, and whether these variables increase (supports) or decrease (contradicts) the probability of an imported case of dengue. With these results, we can infer that for case 1 (i.e. Indonesia-to-Germany, for February 2010), the local model-predicted probability of being an imported case, was 94%. The top five variables influencing this probability were: closeness and betweenness centrality measures of Germany, the incidence rate of dengue and closeness centrality measures of Indonesia and the geographical distance between both countries. Conversely, case 3 (i.e. Tanzania to the United Kingdom, for December 2014), had a similar set of variables as most influential, however, the incidence rate of Tanzania decreases the probability of having an imported case. This demonstrates how variables influencing predictions for a single observation can differ at the local scale. Finally, to ensure the trustworthiness of the local model, we compared the predicted probability of the local model to that of the global optimal model for each observation. There was no difference in the predicted probabilities of the global and local model (analytical comparison not shown). Overall, the local interpretation provides insights into the variations of the individual predictions and provides an important aspect to assuring trust of the model.

• Download figure
• Open in new tab

Figure 6. LIME model plots explaining individual predictions.

Cases 1-4, as presented in Table 3. Each plot shows the 5 most influential variables that best explain the prediction in the local region. The blue bars represent variables that increase the predicted probability (supports), the red bars represent variables that decrease the probability (contradicts).