We used a diverse set of features or variables to forecast vaccine utilization. We obtained daily vaccine utilization data from 710 health facilities across three different regions in Tanzania—Arusha, Tanga and Kilimanjaro. The data had 13 features for every observation. We split the datetime feature into 12 features using the fastai add-datepart function (Howard and Gugger, 2020). While these features can partially contribute to predicting vaccine consumption, we believed that analyzing additional features describing individual facilities and their catchment populations would improve predictive performance. Therefore, we examined the Tanzania Health Facility Registry (THFR, see website reference) and the Tanzania National Bureau of Statistics (TNBS, see website reference) to extract several additional dimensions of data to augment our feature set. We also used a web-based tool (GPS visualizer, see website reference) to add elevation to each health facility. The data we used for building predictive models included several new features such as geo-coordinates, distance to nearest facility, type of facility and regional population (See Table 1 for a comprehensive list of features). Additionally, since recent vaccine utilization can serve as a useful feature, we used a rolling 3-month average as a feature. This resulted in a total of 32 features which went into our machine learning model. These features encompassed key intuitive vaccine utilization determinants — 1) Details of the nature of the vaccine, 2) details of each health facility and 3) features related to the catchment population around each facility.

We used fastai modules—train_cats, add_datepart, and proc_df for initial data pre-processing (Howard and Gugger, 2020). Specifically, we carried out the following steps — 1) Assessing fraction null values for each column. No column had more than 4% null values, 2) For categorical variables, null values were treated as a separate level, and imputed cells were recorded in a separate variable, 3) Median values were used to impute missing continuous variables, 4) Date field was split into 12 separate fields such as time elapsed from the start of the year, is date start of month, etc. All date fields are documented elsewhere (fast.ai, see Howard and Gugger, 2020), and 5) Categorical mapping rules and imputation value to feature mappings were stored in a dictionary which was re-used for pre-processing test data.

We aggregated data per facility and vaccine type into biweekly utilization, which we attempted to forecast. The decision to forecast biweekly utilization was based on discussions with healthcare providers in Tanzania (internal communication). Additionally, summarizing the data into monthly or bimonthly rows does not give us sufficient data size to make robustly extrapolatable forecasts. A factor that led to this decision was the frequency of power outages in those regions. We trained, tuned, and evaluated an RFR to forecast biweekly vaccine utilization at a given health facility in Arusha, Tanga, and Kilimanjaro.

Further, approximately 70% of the data was used to train the model. We equally partitioned the remaining data into validation and test sets. The validation data was used for hyper-parameter tuning whereas the test set was used to report final performance scores. Since the goal was to predict future vaccine utilization, train-validation-test split was done in a temporally sorted manner making sure that the hold-out sets contained only data from dates that were in the future relative to the training, or the validation data.

We used two different measures to evaluate model performance, 1) Root Mean Square Error (RMSE) between predicted and actual biweekly utilization values, and 2) Fraction error (F. E), to measure how far each predicted target variable value was from the actual value of that variable. We also used F. E as the optimization function for our model. We define F. E as F .   E = | ( A . U − P . U ) | ÷ A . U Where AU is actual biweekly utilization and P.U is predicted bimonthly utilization.

We implemented a wide range of commonly used machine learning algorithms to select the best model. Among all the algorithms—regularized linear regression, support vector machine, k-nearest neighbors, RFR and autoregression based univariate time series models, random forests had the best performance on the validation set (Table 2). Our exploratory model building also included multilayered neural networks, including recurrent neural nets; however, the models 1) failed to yield comparable performances to some of the other classic machine learning models and, 2) gave largely uninterpretable predictions with no simple way to find feature importances. This bolstered our choice of RFR for vaccine utilization forecasting.

We trained, tuned, and evaluated an RFR to forecast biweekly vaccine utilization at health facilities in Arusha, and Tanga. We chose RFR for further downstream analyses because, 1) As measured by the Root Mean Square Error (RMSE) RFR outperformed all five of the most commonly used machine learning algorithms and, 2) RFR is a powerful and generalizable predictive framework which can also be leveraged to understand the data better.

We used the validation data to set the values for 1) number of estimators, 2) max features, and 3) min samples leaf. Further, we optimized the model to obtain a minimal validation set RMSE. Our optimized values for these hyperparameters were 40, 0.95, and 7, respectively. We also set n_jobs = −1 in scikit-learn, to effectively utilize all available compute cores.

We made two different versions of our RFR (Leo Brieman, 2001) ↓ model I which uses all 32 features and model II which uses all 32 features except the 3-months rolling average. The two models are both based on random forests but differ in the number of features. Importantly, recent vaccine utilization averages are not used for forecasting in model II.

Model I. We picked this optimized RFR and achieved an F. E of less than 0.2 for a significant subset of facilities, that is, an error of ±2 doses where actual vaccine utilization was 10 doses (Figure 1).

Further, we found that for the majority of health facilities where our predictive performance passed a given threshold for one vaccine, it also passed the same threshold for the other vaccine as well. We were thus able to predict bi-weekly vaccine utilization within an F. E < 0.1 for 22–27% of the facilities depending on the vaccine (Figure 2). Our model performed best for rotavirus vaccine utilization, with almost 50% of facilities approaching a F. E < 0.2.

Model II. This version of the model is more generalizable since all its features can be derived from publicly available resources. It has an F. E < 0.2 for 15–20% of facilities (Supplementary Figure S1). Vaccine type was identified as the most important feature for this model (Supplementary Figure S2).

We used the “mean decrease accuracy” method to calculate importance scores for all features in our model. This involved randomly permuting each column of data, and then calculating the decrease in R ² on the out of the box datasets (Brieman and Friedman. 1984). This feature importance scoring scheme, as implemented in scikit-learn, outputs relative feature importance scores.

Model I. By far, the three-months utilization rolling average has the greatest impact on the model prediction. This is followed by the time, relative to year beginning, to vaccination. A number of features that relate to the health facility come next in our model—public or private, GIS coordinates, district, and ward. We hypothesize that features related to the facility implicitly encode demand characteristics of the catchment population (Figure 3).

Model II. Vaccine type was identified as the most important feature for this model (Supplementary Figure S3). Again, a number of features that relate to the health facility showed up next.

We made sure our model performance evaluation included completely non-overlapping “future” data. We included actual geo-coordinates and altitude to remove subjective levels like “High altitude” or “Low altitude”. Removing the geo-coordinates still gives us a meaningful model with only a small (0.7) decrease in forecasting accuracy. This is hardly surprising since the relative feature importance based on the Random Forest Regressor falls off rather sharply after feature #3 “Births”.

We hypothesize that features related to the facility implicitly encode demand characteristics of the catchment population.

In order to get a measure of vaccine utilization forecasts based on the existing system, we mined stock addition data. We did not find evidence of any statistical forecasting model. The amount of new stock additions immediately prior to our validation and test dates, were treated as “forecasts” made by the baseline, existing system. Since we had actual utilization data for the validation and test period, we were able to calculate F. E and RMSE for the existing system. Evaluation of RFR model performance against this baseline was made on the basis of F. E and RMSE (see Table 3).

The Global Alliance for Vaccines and Immunization anticipates a six-fold increase from 2010 to 2020 in the number of vaccine doses given to complete immunization. Global coverage for basic childhood vaccines has reached a record 86%, but there has been a parallel increase in vaccine wastage, decreasing resource efficiency. Vaccine stock-outs compound the problem by wasting opportunities for immunization. There is no good reason why the correlation between higher rates of immunization and supply chain waste should persist. Further, wastage of all kinds is largely a result of inaccurate, univariate or static models of vaccine demand assumptions; in summation the wrong quantity and type of vaccines at the wrong time. There is thus an enormous economic incentive to reduce vaccine wastage and stock-outs, without sacrificing high immunization coverage rates.

Random Forests (Leo Breiman, 2001) is an algorithm that uses bootstrap to sample multiple data observations or rows from the original data, builds decision trees for each bootstrap sample, then integrates predictions of multiple decision trees, and finally uses majority vote or averaging to arrive at final predictions. The RFR) is conceptualized as a strong predictor combining a bunch of weak predictors.

To our knowledge there has only been one paper where machine learning has been applied to predict vaccine utilization and/or demand (Fruggiero et al., 2012). In that study, the researchers used a combination of autoregressive integrated moving average model and neural networks to forecast annual MMR demand in Taipei County, Taiwan. Specifically, the authors use 10 features to build a decision model, using data related to vaccine demand relative variables and population growth relative variables. Our study differs from theirs in several respects— 1) their goal is to forecast annual demand whereas we aim at forecasting bi-weekly demand, 2) they aim at forecasting demand using variables related to population or stock. No features related to health facility location or the facility itself were included. Our model is significantly more comprehensive, as it includes many granular details of the health facilities including their geo-location, altitude, facility details and 3) they aim at forecasting county wide demand whereas we forecast health facility level demand. Our model is therefore, significantly more fine grained.

Vaccine campaigns sometimes deliver vaccines directly or serve to increase demand. Here, we did not have information on vaccine campaign data. In the next round of data gathering, we may be able to access that data and build an increasingly multidimensional, and more accurate model.