The past several years have brought substantial advances in analytical methods for predictive modeling, even amidst worldwide impacts of Covid-19 vaccines. Gradient boosted trees and deep neural networks have emerged as state-of-the-art methods for time series forecasting. You can review the recent special issue of the International Journal of Forecasting, including a paper by David Lander and me on our winning solution for the M5 Uncertainty competition, in which we exclusively use gradient boosted trees. This note demonstrates an easy way to use the XGBoost Add-in in JMP Pro for basic time series forecasting on a current dataset from the US Center for Disease Control (CDC) at the center of a controversy surrounding Covid-19 vaccine safety.
In a series of three long blog posts, The Ethical Skeptic (TES) raises alarm bells regarding substantial increases in excess cancer and other mortalities since the rollouts of Covid-19 vaccines in early 2021. He uses publicly available data from the CDC, and one of his primary graphs is the following:
In rebuttal, Jeffrey Morris (JM) reanalyzes the data and concludes that the signal does not appear compelling through March 12, 2022. Two of his graphs follow:
JM writes:
The key to assessing his claims is of course to see the 2022 data as they become more mature and we have more complete data through the summer. If the increase in cancer deaths from the nadir in Spring 2021 continues throughout 2022, then this would provide more support for his hypothesis, but if the decrease that appears to have started in early 2022 continues, then this would weaken the claim of inherently higher excess cancer deaths initiated by some factor introduced in early 2021.
We now have the summer 2022 data to analyze in JMP Pro to see where we stand. The data are available in two .csv files covering 2014-2019 and 2020-2022. The latter are deemed “provisional” since new data are continually being added, and the final few weeks tend to have abnormally low numbers since they have not yet been fully tabulated. The subsequent analyses are based on data downloaded on November 22, 2022.
To investigate questions of excess mortality, we need a decent estimate of baseline mortality over time that can be subtracted from current data to estimate the discrepancies. TES describes a fairly complex workflow to accomplish this, including a 6-week smoothing function and a “pull forward effect” to accommodate individuals whose death has been attributed earlier to Covid-19 instead of cancer. Unfortunately, TES does not appear to offer any code to reproduce his results exactly, although the steps in his approach appear sensible. JM estimates baseline as a basic linear trend plus weekly seasonal means (the yellow and blue lines in the left graph above, respectively). He admirably provides R code, which I was able to successfully run and reproduce his results. The baselines of both TES and JM are effectively forecasts of 2020-2022 based on the 2014-2019 data, a reasonable time split given the pandemic began around January 2020.
Given the strong recent successes of gradient boosted trees for time series forecasting, an appealing approach for reanalysis is to use XGBoost to forecast 2020-2022 mortality with the 2014-2019 data and use this as a baseline for assessing excess deaths during 2020-2022. The key to forecasting with XGBoost is to engineer appropriate time-based features that enable generalizability into the future. Depending on the structure of the data and problem, many features are possible; see our IJF paper for some details in this regard. To facilitate comparisons with TES and JM, I use a simple approach with only two numeric features: week of year and 156-week lag values. More sophisticated features are certainly possible, although these two appear to do a reasonably good job and provide a nice initial illustration of the power of XGBoost for problems like this. Results are shown below in Graph Builder:
The top panel shows the raw data in blue and the XGBoost predictions in red. The simple XGBoost forecasts pass the “eyeball” test, although they could very likely be improved with more careful modeling. The middle panel has the residuals from the fit, and results from 2020 forward can be considered estimates of excess cancer deaths for each week. We see the observed data exceed the forecasts for most weeks past mid-2021. The bottom panel has the cumulative sums of the residuals, and the final number of 15,000 is the approximate estimate of total excess deaths for October 1, 2022. An initial conclusion from this analysis provides support for the TES claim of excessive excess cancer deaths from mid-2021 onward, suggesting the vaccines may be playing a causal role.
Estimating excess deaths is a complex topic, and the preceding analysis is designed to be a straightforward starting point that is more illustrative than definitive. The CDC and UK OHID have pages with a lot more depth, and there are numerous analyses and disputes online. For example, see the disagreements between JM and Mathew Crawford on some surprising insurance data exhibiting a regional Simpson’s Paradox, a retracted paper by Kostoff et al., and provocative articles by Norman Fenton and colleagues and Steve Kirsch on data from the UK. It will be interesting to see what data and debates the winter months bring, and thankfully we have a reliable tool like JMP Pro to efficiently and directly analyze data ourselves.
Steps to reproduce the XGBoost forecasting analysis in JMP Pro are as follows, which you can use to perform your own further explorations now or at a later time when more data become available:
- Download, open, and concatenate the two CDC tables (2014-2019 and 2020-2022). Subset the data to Jurisdiction of Occurrence == “United States”. Select the final five rows, and then click Rows > Hide and Exclude. I excluded the final five rows and chose October 1, 2022 as a cutoff date based on an initial Graph Builder display of the raw data. Always plot your data!
- Select rows for years 2020-2022 and then Rows > Row Selection > Name Selection in Column to create a new variable named Val20 that has value 0 for years 2014-2019 (unselected) and 1 for years 2020-2022 (selected). Change the Analysis Type of this variable to Continuous by right-clicking on the red histogram next to it in the far left panel of Columns. Val20 is used as the validation column in XGBoost to hold out 2020-2022 during model fitting.
- The numeric variable MMWR Week is already in the table and is our first feature. To engineer the second feature, select the Malignant neoplasms (C00-C97) column, which is the primary cancer response. In the table, right-click on this variable > New Formula Column > Row > Lag Multiple. Select 156 as both the start and end values then click OK. This creates a new column which is a 3-year lag variable. We use a 3-year lag in order to completely cover the 2020-2022 forecast horizon with observed values from 2014-2019 without any leakage.
- Install the XGBoost add-in. Click Add-Ins > XGBoost > XGBoost. Specify Malignant neoplasms (C00-C97) as the Y Response, MMWR Week and the new 3-year lag feature as X, and Val20 as Validation. Click OK, check Autotune, then Go. Results should be similar to the following.
- In the XGBoost output window, click the red triangle beside Model 1 and click Save Predicteds. This saves predictions for all rows to the JMP table. Those for years 2020-2022 are forecasts since we held out that data during training. Note the plot on the right above also shows that observed values exceed forecasts for 2020-2022, as evidenced by the majority of points lying above the 45-degree reference line.
- In the table, select both Malignant neoplasms (C00-C97) and the new prediction column. Right-click > New Formula Column > Combine > Difference. This creates a new column of residuals. Select and right-click on this new column > New Formula Column > Row > Cumulative Sum.
- In Graph Builder, drag Week Ending Date onto X and the previously created variables onto Y to create the display shown previously.
A JMP table containing all of the data (Steps 1-3) and embedded XGBoost and Graph Builder scripts (Steps 4-7) is attached. Please feel free to contact me regarding the topics here.
