Team:Lambert GA/CALM SUPPLEMENTARY
C A L M
S U P P L E M E N T A R Y
S U P P L E M E N T A R Y
Uniqueness of Approach
The authors have attempted to be as comprehensive as possible in representing the literature on existing models for the Yemeni outbreak. As of the writing of this page, various models, some incorporating machine learning, most not, have been constructed by others. Many of these models are accurate in their specific use cases-such as the one constructed by Jutla, Akanda, and Islam (2010)- but are applied in areas where cholera is seasonal and non-sporadic, such as Bangladesh (Jutla, Akanda, Unnikrishnan, Huq, & Colwell, 2015), and are thus fairly simple (often using various kinds of regression(s) or logistic models and modeling linear relationships). Cholera, in general, is seasonal, but is subject to non-seasonal influences (Emch et al., 2008). In fact, the Yemeni outbreak has been especially subject to many non-seasonal, sporadic influences, namely the Yemeni civil war, necessitating a more complex model that can capture these nonlinear, nonseasonal relations (Camacho et al., 2018). Our extreme gradient boosting approach provides this, offering a robust, principled approach used widely by data scientists to achieve state-of-the-art results on many machine learning challenges (Chen & Guestrin, 2016). The use of machine learning beyond regression is key, as by deriving a deeper understanding of the breadth of data available CALM is able to deliver a more useful forecast. While more complex machine learning algorithms like XGBoost can come at the cost of overfitting, viable complex models are possible without overfitting, as Pezeshki et al. (2016) have demonstrated by predicting cholera in Chabahar City, Iran, using an artificial neural network.
Additionally, forecasts produced by other models often undersupply comprehensiveness, lacking details on when an outbreak might strike and exactly how many will be impacted (for example, Jutla et al. developed a model predicting cholera risk and not cases (Cole, 2018)). In contrast, CALM forecasts the exact number of cholera cases any given Yemeni governorate will experience in 2-week time intervals ranging from 2 weeks to 2 months, providing fundamentally different information than a risk indicator or a broad cumulative incidence count ) to an aid organization or government official.
Finally, existing models often do not make use of the full breadth of cholera-predictive data available, usually making use of only seasonal environmental factors or only cholera incidence. Given that Yemen is currently in a civil war, we propose the incorporation of civil war fatality data along with environmental and epidemiological data to span the entire range of factors that can affect cholera. When paired with extensive feature engineering, CALM’s use of rainfall, past cholera cases and deaths, and civil war fatalities allows it to find key patterns in cholera incidence in Yemen to create a model capable of strongly modeling the nonlinear trends of cholera.
Additionally, forecasts produced by other models often undersupply comprehensiveness, lacking details on when an outbreak might strike and exactly how many will be impacted (for example, Jutla et al. developed a model predicting cholera risk and not cases (Cole, 2018)). In contrast, CALM forecasts the exact number of cholera cases any given Yemeni governorate will experience in 2-week time intervals ranging from 2 weeks to 2 months, providing fundamentally different information than a risk indicator or a broad cumulative incidence count ) to an aid organization or government official.
Finally, existing models often do not make use of the full breadth of cholera-predictive data available, usually making use of only seasonal environmental factors or only cholera incidence. Given that Yemen is currently in a civil war, we propose the incorporation of civil war fatality data along with environmental and epidemiological data to span the entire range of factors that can affect cholera. When paired with extensive feature engineering, CALM’s use of rainfall, past cholera cases and deaths, and civil war fatalities allows it to find key patterns in cholera incidence in Yemen to create a model capable of strongly modeling the nonlinear trends of cholera.
Extra Details
Data
Cholera Epidemiological Data
Cholera case and death statistics are reported by the World Health Organization (WHO) where health experts and researchers work directly with Yemeni health authorities at both the country and local level. Through this direct connection, the WHO is able to record all reported cholera cases and deaths caused by cholera (WHO presence in Yemen, 2018). The data, collected by the WHO, was accessed through the Humanitarian Data Exchange (https://data.humdata.org /group/yem). It provided reports of accumulated new cholera cases and deaths per governorate from up to May 22, 2017, to February 18, 2018. Past cholera cases and deaths were included with the simple assumption that they would be predictive of future cases. Vibrio cholerae requires aquatic environments and can transfer between humans through the transfer of bodily fluids. Thus, the incidence of cholera in one region can indicate the contamination of several food and water sources and therefore indicate a further spread of cholera.
Rainfall
As Vibrio cholerae is indigenous to aquatic environments, rainfall is a significant predictor of the transmission of cholera. In areas exposed to heavy rainfall, through the collapse of sanitary and health infrastructure, interaction between contaminated water and human activities accelerates, resulting in an outbreak (Jutla et al., 2013). This scenario can be found in Yemen, where when exposed to heavy rainfall and deterioration of health facilities, there was a surge in cholera cases (Camacho et al., 2018). Global Lancet Researchers analyzing surveillance date for the Yemen Cholera Outbreak from 2016 to 2018 have found a positive and nonlinear association between weekly rainfall and suspected cholera incidence: the relative risk of cholera 10 days after a weekly rainfall of 25 mm is 42% higher than compared with a week without rain (Camacho et al., 2018). In addition, despite the inability to establish that rainfall is directly causal to the increase in cholera outbreaks, the use of unsafe water sources during the drought season, contamination of water sources during the rainy season, and changing levels of zooplankton and iron in water (which help cholera bacteria survive), may contribute to the increasing levels of cholera during the rainy season (Camacho et al., 2018). These correlations demonstrate the need to measure rainfall in the machine learning model, as rainfall is a predictor for possible climate changes and the corresponding human response and subsequently indicates the spread of cholera in Yemen. Daily rainfall data for Yemen from January 1st ,2017 to March 30th, 2018 was accessed through NASA’s Goddard Earth Sciences Data and Information Services Center (GES DISC), which provides Global Precipitation Measurement data through the Simple Subset Wizard (SSW) database. The Global Precipitation Measurement mission (GPM), launched on February 27th, 2014, is an international network of satellites that use microwave imagers and precipitation radars to measure the volume of rainfall in several regions of the world (Global Precipitation Measurement, 2011). The rainfall data was initially in a netcdf4 format. The 452 files were then parsed and converted to comma-separated-values (CSV). As there were individual data points for every .25 degrees of both latitude and longitude, Reverse geolocation was performed to match coordinates with corresponding Yemeni governorates.
Conflict Data (Yemeni Civil War)
Yemen is currently in the grip of a devastating civil war, which is heavily impacting the cholera crisis in Yemen. While cholera is preventable and treatable under normal circumstances, the collapse of Yemen’s health, water, and sanitation sectors amidst the ongoing armed conflict have fueled the spread of cholera across the country. With direct attacks against hospitals and the bombing of water supplies, the conflict has dissolved 55% of the country's medical, wastewater, and solid waste management infrastructure, making access to clean water and healthcare difficult and expensive (Camacho et al, 2018; Yemen’s Cholera Crisis: Fighting Disease During Armed Conflict, 2017; Yemen: The Forgotten War, 2018). This has led to 15 million Yemenis in need of water and sanitation assistance. Information regarding the status of ongoing conflicts, namely the severity in terms of death toll, was collected with the hope of it being predictive of the region’s infrastructure ability to provide treatments in cholera in the following weeks. Data gathered by the Armed Conflict Location and Event Data Project (ACLED) was retrieved from the Humanitarian Data Exchange (https://data.humdata.org/group/yem). ACLED reported the type of conflict, agents, locations, dates, and other characteristics of the politically charged conflict from January 1, 2016, to June 6, 2018 (Raleigh and Dowd, 2017). The number of daily casualties due to conflict in each Yemeni governorate was used as a metric for civil war related violence.
Methods
Dataset Preparation
With the objective of predicting new cholera cases in any given governorate in Yemen from week to week, there were a number of steps taken to prepare the data. In order to produce models that did not simply rely on seasonal trends and were able to predict spikes in cholera cases, the case and death report time series were made stationary through temporal differencing. It should be noted that the country of Yemen encompasses 21 governorates or administrative divisions. While the CALM models were trained on data from all 21 governorates, data preparation on each governorate was performed separately to preserve each governorate’s unique time series. As the interval between each WHO cholera case/death report was not standard, the data was linearly interpolated into a daily time series. The Yemeni Cholera outbreak is seasonal and endemic as outbreaks spike during the rainy season (April-August) - however, the outbreaks rely on non seasonal factors such as conflict and damage to health and sanitation (Camacho et al., 2018a). Parsing data required finding the number of new cholera cases in a single day, given the total number of cases in the previous day. The values were then normalized by the population of each governorate (e.g new cases per 10,000 people). Finally, we calculated our four target variables: the number of new cholera cases 0-2 weeks from the present day, 2-4 weeks from the present, 4-6 weeks from the present, and 6-8 weeks from the present.
Our dataset was split into three portions: training, cross-validation, and a hold-out test set. The hold-out set was left untouched until the completion of our methods to provide an accurate real-world simulation of our models’ performance. Our base training set was defined from July 1 to August 15th. While WHO reports extended back as far as May 22, we chose to start on July 1 in order to have enough prior data for feature calculation. Our cross-validation dataset was defined from August 15 to November 10. Finally, our hold-out set started from November 11 and extended to a final date in January/February, which varied for each defined target variable depending on the respective range: a 6-8 week forecast implies a larger time frame between current and forecast date than a 2-4 week forecast, and so the 6-8 week forecast holdout set would end prior to the 2-4 week forecast. It may seem that the cross-validation set outweighs the training set significantly, but this was mitigated with the use of a rolling window forecast - a gold standard for cross-validation in time series forecasting. Rolling window cross-validation is easiest understood with the following example. Given a dataset spanning four weeks, a rolling window forecast would dictate that we train on the first week, predict on the second week, then train on the first two weeks, predict on the third, and finally train on the first three weeks and predict on the fourth. In this example, the first week would be the base-training set (as it was never predicted on and was included in the training set of each fold), the second and third weeks the cross-validation set (as they varied between prediction and training sets), and the fourth the week the hold-out set (as it was never trained on). Our five cross-validation sets were defined as follows: August 16 to August 31, August 31 to September 15, September 15 to September 30, September 30 to October 15, and Finally October 15 to October 30 (it should be noted that the final fold included data from October 30 to November 10 as a prediction set, though this does not cross into the hold-out set). The cross-validation sets were used to select features and find optimal hyperparameters for our model, and the hold-out set was used to simulate real-world performance of our model.
Feature Engineering and Tuning
Feature engineering is the crux of applied machine learning, and so we went through an exhaustive feature extraction and selection process in order to arrive at our final features. First, we extracted 45,000 potentially relevant features using the tsFresh package, which calculates an expansive array of time series features on our data (Christ et al., 2018). The objective of calculating these many features was the hope to capture ideal representations of our data: while the majority of these features would not be used in the final model, our coverage of this expansive set allowed us to ensure the best features would be found. We also calculated features over a series of overlapping time frames in order to provide varying frames of reference and lags: 8 weeks prior, 6 weeks prior, 4 weeks prior, 2 weeks prior, and 1 week prior. Features describing geographically neighboring governorates (through taking the mean) were also calculated. While having more data is usually beneficial, in this case, our number of training examples was far outnumbered by the number of features. Therefore, a demanding feature selection process was required. Using tsFresh’s scalable hypothesis tests with a false discovery rate of 0.001, we were able to calculate features statistically relevant to each time-range prediction, providing us with four sets of features ~15,000 in number for each time-frame prediction. Next, we removed collinear features, or those that were 97% correlated with each other, as these features would be redundant to our model. This provided us with sets of ~10,000 features to further narrow. We trained and tuned an extreme gradient boosting model, XGBoost, to rank the features in order of importance for each time-range prediction. Utilizing the ranking produce, we recursively added features based on if they added to our cross-validation loss (the root mean square error across all five cross-validation folds). This allowed us to arrive at the best 30-50 features for each time-range. All in all, we were able to remove ~99.9% of our original features.
Model
We utilized XGBoost, a random forest-based, extreme gradient boosting algorithm, to construct each of our models. Through bootstrap aggregation, the construction of multiple (often hundreds) of decision trees that are trained on random subsets of the data and then collectively vote for the final prediction, XGBoost is able to address variance-related error (overfitting). XGBoost also addresses the converse, bias-related error (underfitting), through gradient boosting: the process by which each decision tree is constructed with a greater focus on the samples the prior trees had difficulties with (Chen and Guestrin, 2016). As opposed to simpler regression techniques utilized by previous models (refer to the background), XGBoost is able to gain a far deeper understanding of the data through nonlinear relations (while being able to distinguish from noise), making it an ultimately more robust choice of algorithm.
Tuning
We utilized Bayesian Optimization to find optimal hyperparameters for our model. In contrast with a brute-force search over a defined set of hyperparameters, Bayesian Optimization tracks prior evaluations to form probabilistic assumptions on an objective function given a set of hyperparameters, allowing informed choices to be made on which hyperparameters to try (Snoek et al., 2012). This allowed us to converge at optimal hyperparameters with far greater efficiency.
Feature Engineering Results
Feature engineering is integral to the application of machine learning by transforming raw data into feature vectors that cumulatively affect the model’s output. Feature engineering allows the conversion of variables to numeric values in order to forecast future incidence. For instance, instead of using individual cases of cholera, the model uses multiple numerical representations of this data to gain further insights into the case data. Lambert iGEM utilized this technique to build the machine learning model CALM. In order to make a prediction, CALM takes into account multiple features of varying importance. Lambert iGEM used multiple unique timeseries in the design of CALM. Features were calculated for each governorate and for governorates neighboring the respective governorate. Features were also calculated over multiple time frames: 8 weeks prior, 6 weeks prior, 4 weeks prior, 2 weeks prior, and 1 week prior. Data from which features were calculated includes conflict fatalities, rainfall, past cases, and past deaths.
Cholera Epidemiological Data
Cholera case and death statistics are reported by the World Health Organization (WHO) where health experts and researchers work directly with Yemeni health authorities at both the country and local level. Through this direct connection, the WHO is able to record all reported cholera cases and deaths caused by cholera (WHO presence in Yemen, 2018). The data, collected by the WHO, was accessed through the Humanitarian Data Exchange (https://data.humdata.org /group/yem). It provided reports of accumulated new cholera cases and deaths per governorate from up to May 22, 2017, to February 18, 2018. Past cholera cases and deaths were included with the simple assumption that they would be predictive of future cases. Vibrio cholerae requires aquatic environments and can transfer between humans through the transfer of bodily fluids. Thus, the incidence of cholera in one region can indicate the contamination of several food and water sources and therefore indicate a further spread of cholera.
Rainfall
As Vibrio cholerae is indigenous to aquatic environments, rainfall is a significant predictor of the transmission of cholera. In areas exposed to heavy rainfall, through the collapse of sanitary and health infrastructure, interaction between contaminated water and human activities accelerates, resulting in an outbreak (Jutla et al., 2013). This scenario can be found in Yemen, where when exposed to heavy rainfall and deterioration of health facilities, there was a surge in cholera cases (Camacho et al., 2018). Global Lancet Researchers analyzing surveillance date for the Yemen Cholera Outbreak from 2016 to 2018 have found a positive and nonlinear association between weekly rainfall and suspected cholera incidence: the relative risk of cholera 10 days after a weekly rainfall of 25 mm is 42% higher than compared with a week without rain (Camacho et al., 2018). In addition, despite the inability to establish that rainfall is directly causal to the increase in cholera outbreaks, the use of unsafe water sources during the drought season, contamination of water sources during the rainy season, and changing levels of zooplankton and iron in water (which help cholera bacteria survive), may contribute to the increasing levels of cholera during the rainy season (Camacho et al., 2018). These correlations demonstrate the need to measure rainfall in the machine learning model, as rainfall is a predictor for possible climate changes and the corresponding human response and subsequently indicates the spread of cholera in Yemen. Daily rainfall data for Yemen from January 1st ,2017 to March 30th, 2018 was accessed through NASA’s Goddard Earth Sciences Data and Information Services Center (GES DISC), which provides Global Precipitation Measurement data through the Simple Subset Wizard (SSW) database. The Global Precipitation Measurement mission (GPM), launched on February 27th, 2014, is an international network of satellites that use microwave imagers and precipitation radars to measure the volume of rainfall in several regions of the world (Global Precipitation Measurement, 2011). The rainfall data was initially in a netcdf4 format. The 452 files were then parsed and converted to comma-separated-values (CSV). As there were individual data points for every .25 degrees of both latitude and longitude, Reverse geolocation was performed to match coordinates with corresponding Yemeni governorates.
Conflict Data (Yemeni Civil War)
Yemen is currently in the grip of a devastating civil war, which is heavily impacting the cholera crisis in Yemen. While cholera is preventable and treatable under normal circumstances, the collapse of Yemen’s health, water, and sanitation sectors amidst the ongoing armed conflict have fueled the spread of cholera across the country. With direct attacks against hospitals and the bombing of water supplies, the conflict has dissolved 55% of the country's medical, wastewater, and solid waste management infrastructure, making access to clean water and healthcare difficult and expensive (Camacho et al, 2018; Yemen’s Cholera Crisis: Fighting Disease During Armed Conflict, 2017; Yemen: The Forgotten War, 2018). This has led to 15 million Yemenis in need of water and sanitation assistance. Information regarding the status of ongoing conflicts, namely the severity in terms of death toll, was collected with the hope of it being predictive of the region’s infrastructure ability to provide treatments in cholera in the following weeks. Data gathered by the Armed Conflict Location and Event Data Project (ACLED) was retrieved from the Humanitarian Data Exchange (https://data.humdata.org/group/yem). ACLED reported the type of conflict, agents, locations, dates, and other characteristics of the politically charged conflict from January 1, 2016, to June 6, 2018 (Raleigh and Dowd, 2017). The number of daily casualties due to conflict in each Yemeni governorate was used as a metric for civil war related violence.
Methods
Dataset Preparation
With the objective of predicting new cholera cases in any given governorate in Yemen from week to week, there were a number of steps taken to prepare the data. In order to produce models that did not simply rely on seasonal trends and were able to predict spikes in cholera cases, the case and death report time series were made stationary through temporal differencing. It should be noted that the country of Yemen encompasses 21 governorates or administrative divisions. While the CALM models were trained on data from all 21 governorates, data preparation on each governorate was performed separately to preserve each governorate’s unique time series. As the interval between each WHO cholera case/death report was not standard, the data was linearly interpolated into a daily time series. The Yemeni Cholera outbreak is seasonal and endemic as outbreaks spike during the rainy season (April-August) - however, the outbreaks rely on non seasonal factors such as conflict and damage to health and sanitation (Camacho et al., 2018a). Parsing data required finding the number of new cholera cases in a single day, given the total number of cases in the previous day. The values were then normalized by the population of each governorate (e.g new cases per 10,000 people). Finally, we calculated our four target variables: the number of new cholera cases 0-2 weeks from the present day, 2-4 weeks from the present, 4-6 weeks from the present, and 6-8 weeks from the present.
Our dataset was split into three portions: training, cross-validation, and a hold-out test set. The hold-out set was left untouched until the completion of our methods to provide an accurate real-world simulation of our models’ performance. Our base training set was defined from July 1 to August 15th. While WHO reports extended back as far as May 22, we chose to start on July 1 in order to have enough prior data for feature calculation. Our cross-validation dataset was defined from August 15 to November 10. Finally, our hold-out set started from November 11 and extended to a final date in January/February, which varied for each defined target variable depending on the respective range: a 6-8 week forecast implies a larger time frame between current and forecast date than a 2-4 week forecast, and so the 6-8 week forecast holdout set would end prior to the 2-4 week forecast. It may seem that the cross-validation set outweighs the training set significantly, but this was mitigated with the use of a rolling window forecast - a gold standard for cross-validation in time series forecasting. Rolling window cross-validation is easiest understood with the following example. Given a dataset spanning four weeks, a rolling window forecast would dictate that we train on the first week, predict on the second week, then train on the first two weeks, predict on the third, and finally train on the first three weeks and predict on the fourth. In this example, the first week would be the base-training set (as it was never predicted on and was included in the training set of each fold), the second and third weeks the cross-validation set (as they varied between prediction and training sets), and the fourth the week the hold-out set (as it was never trained on). Our five cross-validation sets were defined as follows: August 16 to August 31, August 31 to September 15, September 15 to September 30, September 30 to October 15, and Finally October 15 to October 30 (it should be noted that the final fold included data from October 30 to November 10 as a prediction set, though this does not cross into the hold-out set). The cross-validation sets were used to select features and find optimal hyperparameters for our model, and the hold-out set was used to simulate real-world performance of our model.
Feature Engineering and Tuning
Feature engineering is the crux of applied machine learning, and so we went through an exhaustive feature extraction and selection process in order to arrive at our final features. First, we extracted 45,000 potentially relevant features using the tsFresh package, which calculates an expansive array of time series features on our data (Christ et al., 2018). The objective of calculating these many features was the hope to capture ideal representations of our data: while the majority of these features would not be used in the final model, our coverage of this expansive set allowed us to ensure the best features would be found. We also calculated features over a series of overlapping time frames in order to provide varying frames of reference and lags: 8 weeks prior, 6 weeks prior, 4 weeks prior, 2 weeks prior, and 1 week prior. Features describing geographically neighboring governorates (through taking the mean) were also calculated. While having more data is usually beneficial, in this case, our number of training examples was far outnumbered by the number of features. Therefore, a demanding feature selection process was required. Using tsFresh’s scalable hypothesis tests with a false discovery rate of 0.001, we were able to calculate features statistically relevant to each time-range prediction, providing us with four sets of features ~15,000 in number for each time-frame prediction. Next, we removed collinear features, or those that were 97% correlated with each other, as these features would be redundant to our model. This provided us with sets of ~10,000 features to further narrow. We trained and tuned an extreme gradient boosting model, XGBoost, to rank the features in order of importance for each time-range prediction. Utilizing the ranking produce, we recursively added features based on if they added to our cross-validation loss (the root mean square error across all five cross-validation folds). This allowed us to arrive at the best 30-50 features for each time-range. All in all, we were able to remove ~99.9% of our original features.
Model
We utilized XGBoost, a random forest-based, extreme gradient boosting algorithm, to construct each of our models. Through bootstrap aggregation, the construction of multiple (often hundreds) of decision trees that are trained on random subsets of the data and then collectively vote for the final prediction, XGBoost is able to address variance-related error (overfitting). XGBoost also addresses the converse, bias-related error (underfitting), through gradient boosting: the process by which each decision tree is constructed with a greater focus on the samples the prior trees had difficulties with (Chen and Guestrin, 2016). As opposed to simpler regression techniques utilized by previous models (refer to the background), XGBoost is able to gain a far deeper understanding of the data through nonlinear relations (while being able to distinguish from noise), making it an ultimately more robust choice of algorithm.
Tuning
We utilized Bayesian Optimization to find optimal hyperparameters for our model. In contrast with a brute-force search over a defined set of hyperparameters, Bayesian Optimization tracks prior evaluations to form probabilistic assumptions on an objective function given a set of hyperparameters, allowing informed choices to be made on which hyperparameters to try (Snoek et al., 2012). This allowed us to converge at optimal hyperparameters with far greater efficiency.
Feature Engineering Results
Feature engineering is integral to the application of machine learning by transforming raw data into feature vectors that cumulatively affect the model’s output. Feature engineering allows the conversion of variables to numeric values in order to forecast future incidence. For instance, instead of using individual cases of cholera, the model uses multiple numerical representations of this data to gain further insights into the case data. Lambert iGEM utilized this technique to build the machine learning model CALM. In order to make a prediction, CALM takes into account multiple features of varying importance. Lambert iGEM used multiple unique timeseries in the design of CALM. Features were calculated for each governorate and for governorates neighboring the respective governorate. Features were also calculated over multiple time frames: 8 weeks prior, 6 weeks prior, 4 weeks prior, 2 weeks prior, and 1 week prior. Data from which features were calculated includes conflict fatalities, rainfall, past cases, and past deaths.