2.1. Data sources

The incidence and mortality data for COVID-19 in the United States were obtained from the Johns Hopkins University Data Repository for COVID-19 (COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, 2025) that collected aggregated data from multiple sources for multiple countries. For the USA specifically, these data were accumulated from health departments nationwide. Data counts were available daily at the county and state levels. At the same time, the reporting relied on the case definition issued by the CDC and was updated twice throughout the course of the pandemic (Coronavirus Disease 2019 (COVID-19), 2020). The dataset included both COVID-19 cumulative incidence and cumulative related mortality. The reported cases were available from January 22, 2020 until March 9, 2023. The available daily COVID-19 incidence and mortality data were aggregated into monthly counts for two main reasons: to reduce numerical instabilities, noise, and outliers caused by irregular reporting (such as sporadic delays or gaps on weekends and holidays, which vary between states); and, more importantly, because Google Trends data covering multiple years are only available at a monthly frequency. It should be noted, however, that the methods presented and utilized can be applied to any equally spaced time intervals, such as days or weeks, provided that reliable data are available for the early stages of future pandemics. The application of these methods remains identical. The only potentially necessary adjustment for shorter time intervals involves modifying the length of the lags used in the cross-correlations, as discussed further.

The Google Trends search queries data are publicly available and can be downloaded from the corresponding website (Google Trends, 2025). These data are anonymized (ensuring the absence of personal data), categorized (assigning topics to search queries), and aggregated (grouped) (FAQ about Google Trends data - Trends Help, 2025). In particular, Google Trends data do not contain complete individual searches but instead contain information about random samples from all search queries within a given period. While the sample sizes used for such summaries and the details of the sampling method are not disclosed, the Google Trends summaries are claimed to be sufficient for a selected geographic region, i.e. states in the presented analysis, along with Washington, DC. For each preselected keyword or key phrase, the resulting data are available as a scaled time series with a time interval length dependent on the length of the requested timeframe chosen to generate the series. In particular, for requested reporting periods longer than twelve months, the reporting time interval duration is one month. The numbers returned by Google Trends for a given search keyword or key phrase are already scaled to a range of 0 to 100, with 100 representing the highest value for the selected timeframe and 0 representing the lowest (Rogers, Reference Rogers2016). In this study, for each of the preselected keywords and key phrases, for each geographic location, the resulting time series is used as an auxiliary predictor variable, which is used to improve the prediction of incidence rates and mortality counts. In particular, the Google Trends data were downloaded monthly between January 2017 and January 2023 to cover 5 years. However, only monthly records from January 2020 corresponded with the dates when COVID-19 reporting in the United States started. Data were downloaded using modified scraper code from Emiliano Mastragostino’s GitHub repository (Scraper for extracting data from Google Trends, 2025). The modified scraper code has been made available on GitHub, as well as the rest of the project code (GitHub repository for the Google Trends Study, 2025). The source Google Trends can be redownloaded from the official website (Google Trends, 2025) using this scraper. The complete lists of the keywords selected for the study for incidence and mortality predictions are summarized in Supplementary Tables S1 and S2.

2.2. Keywords selection and time series models

2.2.2. Granger causality test

The first approach to selecting a keyword pool was based on the Granger causality test. The Granger causality test is a statistical technique used to investigate how helpful a one time series is for predicting another time series (Granger, Reference Granger1969). This test does not test or justify any causal relationships. Instead, it provides the utility of one time series for predicting the other series. The test is based on the idea that if the values of series $ X(t) $ significantly (using some criteria) improve the prediction of the values of series $ Y(t) $ beyond the information contained in the lagged values of series $ Y(t) $ alone, then $ X(t) $ helps predict $ Y(t) $ .

The test has been implemented via the function Granger test from the R package lmtest. The Granger causality test is performed as a two-step procedure. In the first step, autoregressive models are fitted individually for search frequency and COVID-19 incidence or mortality. These models only include lagged values of the respective variable as predictors (Granger, Reference Granger1969):

(1) $$ X(t)=\sum \limits_{i=1}^p{\unicode{x03D5}}_{X,i}X\left(t-i\right)+{\unicode{x025B}}_X(t) $$

(2) $$ Y(t)=\sum \limits_{i=1}^p{\unicode{x03D5}}_{Y,i}Y\left(t-i\right)+{\unicode{x025B}}_Y(t) $$

where $ p $ represents the number of lagged terms included in the model, $ {\unicode{x03D5}}_{X,i} $ and $ {\unicode{x03D5}}_{Y,i} $ are the estimated autoregressive coefficients, and $ {\unicode{x025B}}_X(t) $ and $ {\unicode{x025B}}_Y(t) $ are the error terms.

In the second step, the corresponding autoregressive models are extended by including the lagged values of the other trend, which is either the search trend or the incidence or mortality series, as predictors (Granger, Reference Granger1969):

(3) $$ X(t)=\sum \limits_{i=1}^p{\unicode{x03D5}}_{X,i}X\left(t-i\right)+\sum \limits_{i=1}^q{\unicode{x03B8}}_{X,i}Y\left(t-i\right)+{\unicode{x025B}}_X(t) $$

(4) $$ Y(t)=\sum \limits_{i=1}^p{\unicode{x03D5}}_{Y,i}Y\left(t-i\right)+\sum \limits_{i=1}^q{\unicode{x03B8}}_{Y,i}X\left(t-i\right)+{\unicode{x025B}}_Y(t) $$

where $ q $ represents the number of lagged terms of the other series included in the model, $ {\theta}_{X,i} $ and $ {\theta}_{Y,i} $ are the coefficients representing the influence of the lagged values of search frequency and incidence or mortality, respectively.

The significance of the Granger causality test is evaluated via the corresponding hypothesis testing (Granger, Reference Granger1969; Hamilton, Reference Hamilton2020). More specifically, the null hypothesis ( $ {H}_0 $ ) states that Google Trend frequency does not “Granger cause” the corresponding COVID-19 incidence or mortality, that is, it does not help predict COVID-19 incidence and mortality, implying that the additional information provided by the trend does not improve the prediction of incidence (or mortality) beyond the lagged values of incidence (or mortality) alone. The alternative hypothesis suggests the presence of Granger causality, that is, that the corresponding Google Trend frequency helps predict COVID-19 incidence or mortality. The significance is evaluated using the F-test. The F-statistic is computed by comparing the squared residuals from the restricted model (without the additional predictors from the second series) to the sum of squared residuals from the extended model (with the additional predictors from the second series). The corresponding test result has been considered statistically significant, that is, the presence of Granger causality, if the corresponding test p-value was below $ \alpha =0.05 $ to determine the presence of statistically significant Granger causality. A p-value below this threshold indicated that search frequency had a significant causal influence on COVID-19 incidence. No multiple testing adjustment was performed since the goal was to include the most predictors possible in the available set of potentially useful ones. This work compares the response series $ Y(t) $ with multiple Google Trends series $ X(t) $ using pairwise Granger causality tests.

2.2.3. Cross-correlation coefficient

The second approach for selecting a pool of keywords involved calculating cross-correlation coefficient values to assess the correlation between the corresponding time series.

The cross-correlation coefficient is a statistic used to assess the strength of the relationship between two time series for a specific lag value (Box et al., Reference Box, Jenkins, Reinsel and Ljung2015). More precisely, the cross-correlation function at an integer-valued time lag of $ \unicode{x03C4} $ , denoted as $ {\unicode{x03C1}}_{XY}\left(\unicode{x03C4} \right) $ , represents the correlation coefficient between the series $ X(t) $ and $ Y(t) $ . The estimate of the cross-correlation coefficient $ {\hat{\unicode{x03C1}}}_{XY}\left(\unicode{x03C4} \right) $ is calculated as (Box et al., Reference Box, Jenkins, Reinsel and Ljung2015):

(5) $$ {\hat{\unicode{x03C1}}}_{XY}\left(\unicode{x03C4} \right)=\frac{1}{n-\mid \unicode{x03C4} \mid}\sum \limits_{t=1}^{n-\mid \unicode{x03C4} \mid}\left[X(t)-\overline{X}(t)\right]\cdot \left[Y\left(t+\unicode{x03C4} \right)-\overline{Y}\left(t+\unicode{x03C4} \right)\right] $$

where $ n $ is the number of time slots in the original series, $ X(t) $ and $ Y(t) $ are the values of Google Trend and COVID-19 incidence at time $ t $ , $ \unicode{x03C4} $ is the integer time lag, and $ \overline{X}(t) $ and $ \overline{Y}\left(t+\unicode{x03C4} \right) $ are the means of the original time series and the version shifted by the value of $ \unicode{x03C4} $ . For a lag value of $ \unicode{x03C4} >0 $ , the number of time slots in the original series used to compute the correlation coefficient decreases by $ \unicode{x03C4} $ , respectively.

In the presented study, the cross-correlation coefficients were computed between pairs of series for different values of the lag $ \unicode{x03C4} $ to investigate the strength of the association between selected Google Trend keywords and COVID-19 incidence and mortality. Multiple lag values $ \unicode{x03C4} $ were considered based on the hypothesis that changes in search frequency might precede or follow corresponding changes in incidence and mortality. Since the correlation values range from $ -1 $ to $ 1 $ , the goal is to focus on those whose trends’ absolute values are close to $ 1 $ , and subsequently use those keywords as predictors.

2.2.4. Data rescaling for comparisons

Since the time series were on different scales (i.e., actual counts for incidence and mortality data versus re-scaled Google Trends values on a scale from 0 to 100), all the series used in the analysis were rescaled and normalized. Suppose exact predictions are needed for practical purposes. In that case, the one-to-one backward transformation can be applied to the data resulting from predictions made by the models fitted to the rescaled and normalized data. In particular, the following normalization was used:

(6) $$ {Y}^{\prime }(t)=\frac{Y(t)}{\overline{Y}} $$

where $ {Y}^{\prime }(t) $ is the rescaled value at time point $ t $ , $ Y(t) $ is the original value at time point $ t $ , and $ \overline{Y} $ is the average value of that time series across the available time slots. Since both considered time series (incidence or mortality vs. Google Trends) are non-negative, the presented standardization method allows for the scaling of the series without causing any values to become negative.

2.3. Prediction models and performance metrics

The processed and standardized data were used for prediction models with and without Google Trends as external predictors. In particular, time series methods were used, including Autoregressive Integrated Moving Average with Exogenous Variables (sometimes referred to as either ARIMAX to emphasize the inclusion of predictors or, just as the generic name ARIMA) and Prophet predictive models. In addition, the eXtreme Gradient Boosting method (XGBoost model), based on decision trees, has been adopted for time series analysis. Although numerous methods are available for estimation and prediction, the three models evaluated here enable an assessment of predictor robustness across distinct established modeling frameworks, including likelihood-based (ARIMA), Bayesian (Prophet), and machine learning (XGBoost) approaches.

To evaluate and compare the models’ prediction performances, they were fitted to specific calibration periods, and predictions were considered for the 2 months following these calibration periods. The predictions for those two months were compared with the actual historical data for those two months, which were not used for model calibration. The predictions were compared with the data using the mean absolute error (MAE). The length of the calibration period gradually increased, starting from 12 months and going up to 34 months (i.e., 36 months of data – 2 months of prediction), where 36 months corresponds to the total number of months in the dataset. The increased calibration period mirrors real-life settings, where the calibration window expands as more data about the ongoing epidemic typically become available. The smallest number of calibration points was specifically selected to be 12 to allow the models to capture the inherent annual seasonality of COVID-19.

More specifically, to generate predictions, the number of data points (i.e., months) for each calibration period gradually increased for each model and state combination. The dataset contained 36 months of data. Initially, a 2-month forecast was made using 12 months of historical data to capture a year of possible seasonality. An additional month of data were incorporated for each subsequent prediction, with the final prediction using 34 data points. This process yielded 23 two-month forecasts for each model and state combination. Prediction accuracy was then assessed by comparing the predicted values to the actual observed data for each length of data used for model fitting. Finally, median values were calculated across all states for each model and calibration length, and summaries were produced.

In summary, the analysis directly incorporates evaluations of the robustness of the presented approach both temporally (using multiple training sets ranging from 12 to 34 months with 1-month increments) and geographically (across all 50 states and the District of Columbia, each with distinct COVID-19 dynamics), leveraging a rich dataset.

Specifically, to address temporal robustness, models have been trained on progressively expanding monthly time series and then used to predict the subsequent months. This approach allowed thorough evaluation of model performance over multiple forecast horizons and different phases of the pandemic, including the early stages, which the authors consider among the most crucial from an epidemiological perspective.

To address geographic robustness, each expanding dataset has been evaluated across all 50 US states to account for geographic heterogeneity and to examine the consistency of trends across regions.

In total, the number of evaluations was 51 (locations) $ \times $ (34 – 12 + 1) (time intervals) = 1173 for each keyword combination considered in the analysis for each incidence and mortality. These evaluations were repeated for 15 selected keyword combinations, resulting in 1173 $ \times $ 15 = 17,595 model fits performed separately for incidence and mortality outcomes, as well as for each of the three models considered. For the PCA-transformed predictors, the number of comparisons was also 1173 model fits for each model type separately for incidence and mortality outcomes.

2.4. Autoregressive integrated moving average with exogenous variables (ARIMAX) model

The autoregressive integrated moving average (ARIMA) model is a likelihood-based model for analyzing time series data, used for predictions by incorporating past observations and moving averages. More specifically, the ARIMA model is characterized by three primary components: autoregression (AR), differencing or integration (I), and moving averages (MA) (Shumway and Stoffer, Reference Shumway and Stoffer2017). The ARIMA model is formulated as:

(7) $$ y(t)=c+\sum \limits_{i=1}^p{\unicode{x03D5}}_i\hskip0.1em y\left(t-i\right)+\sum \limits_{j=1}^q\;{\unicode{x03B8}}_j\hskip0.1em {\unicode{x025B}}_{t-j}+{\unicode{x025B}}_t. $$

In formula (7), values $ y(t) $ denote the series values at time $ t $ (possibly after differencing transformations to ensure stationarity, if needed), and $ {\unicode{x03D5}}_i $ are the corresponding autoregressive coefficients, which are estimated. The coefficients $ {\unicode{x03B8}}_j $ are the moving average parameters, which are also estimated. At the same time, $ {\unicode{x025B}}_t $ are the residual errors from the moving average model at time $ t $ . The lowercase values $ y(t) $ emphasize the potential preprocessing differencing transformations of the original values $ Y(t) $ , where for $ d=0 $ there are no transformations, that is, $ y(t)=Y(t) $ ; for $ d=1 $ , the differenced observations are $ y(t)=Y(t)-Y\left(t-1\right) $ , and so on, with higher orders applying the same transformation to the values obtained in the previously applied transformation.

When additional explanatory variables (exogenous variables, which are Google Trends time series in this work) are included in the ARIMA framework, the model is extended to an ARIMAX (autoregressive integrated moving average with exogenous variables) model. The model has the form:

(8) $$ y(t)=c+\sum \limits_{i=1}^p{\unicode{x03D5}}_i\hskip0.1em y\left(t-i\right)+\sum \limits_{j=1}^q{\unicode{x03B8}}_j\hskip0.1em {\unicode{x025B}}_{t-j}+\sum \limits_{k=1}^m{\unicode{x03B2}}_k\hskip0.1em {X}_k\left(t-\unicode{x03C4} \right)+{\unicode{x025B}}_t, $$

where $ {X}_k\left(t-\unicode{x03C4} \right) $ denotes the value of the $ k $ -th predictors series value at time $ t $ with the integer lag value $ \unicode{x03C4} \ge 0 $ , and $ {\unicode{x03B2}}_k $ for $ k=1,2,\dots, m $ is the set of corresponding estimated coefficients for those predictors.

The presented formulation allows the ARIMAX model (8) to account for external predictors, improving the prediction of the studied time series. For this study, the auto.arima function from the R package forecast was used. This function automatically selects the best-fitting parameters of the model based on the Akaike information criterion (AIC). It can incorporate external predictors to fit the ARIMAX model. The numbers of parameters $ p $ and $ q $ were searched with the maximum number corresponding to the length of the calibration period, which was twelve points. Box–Cox transformations were also used to improve the stability of the model input data.

On top of that, to evaluate the robustness of the ARIMA-type models, the autoregressive fractionally integrated moving average (ARFIMA) model was also considered using the arfima function to assess potential differences in the model results.

2.5 Prophet model

The Prophet model was developed by Meta-Facebook’s Core Data Science team for business forecasting (Taylor and Letham, Reference Taylor and Letham2018). It is a GAM variant with nonlinear smoothers applied to the regressors. This tool was designed for the user without expert knowledge about time series and their specific methods. It can handle multiple seasonalities with high magnitude, irregular influences like holidays, missing data, outliers, and changes in general trends. While it was created with business tasks in mind, it proved its usefulness in public health research as well (Navratil and Kolkova, Reference Navratil and Kolkova2019; Samal et al., Reference Samal, Babu, Das and Acharaya2019; Satrio et al., Reference Satrio, Darmawan, Nadia and Hanafiah2021; Kirpich et al., Reference Kirpich, Shishkin, Weppelmann, Tchernov, Skums and Gankin2022; Shishkin et al., Reference Shishkin, Lhewa, Yang, Gankin, Chowell, Norris, Skums and Kirpich2023; Bleichrodt et al., Reference Bleichrodt, Luo, Kirpich and Chowell2024). Prophet model is formulated as the sum of the following components:

(9) $$ Y(t)=g(t)+s(t)+h(t)+\unicode{x025B} (t), $$

where $ Y(t) $ is the COVID cases or death at time point $ t $ , $ g(t) $ is the trend component, $ s(t) $ is the periodical seasonal component, $ h(t) $ is the irregular component (like holidays, for example), and $ \unicode{x025B} (t) $ is modeling error at time $ t $ . Trend component $ g(t) $ is a logistic growth model that incorporates change points in which growth is allowed to change. Those change points could be specified by the model user or selected automatically. The periodical seasonal component $ s(t) $ is a Fourier series, with parameters selected using AIC procedure. Irregular component $ h(t) $ directly adds value at a particular time. The model also has embedded uncertainty estimation.

In this article, Prophet function from the package with the same name for R programming language was used. This implementation can use exogenous variables, so we incorporated Google Trends data in the model. The following model parameters were set: $ yearly. seasonality $ to 4, and $ weekly. seasonality $ as well as $ daily. seasonality $ disabled. Other parameters were not set, so the function used default values described in the package documentation.

2.6. XGBoost model

XGBoost (eXtreme Gradient Boosting) is an implementation of gradient-boosted decision trees (GBDT), widely used for multiple tasks, including classification, regression, and (most relevant for this project) time series forecasting, including epidemiologic forecasting (Chen and Guestrin, Reference Chen and Guestrin2016; Wang and Guo, Reference Wang and Guo2020; Rahman et al., Reference Rahman, Chowdhury and Amrin2022). The model sequentially builds weak learners, typically decision trees (Gehrke and Ye, Reference Gehrke and Ye2003), to correct the errors of previous weak learners and to minimize a specified loss function.

In brief, for a set of $ n $ records and $ m $ predictors per each record represented by the set:

(10) $$ \mathcal{D}=\left\{\left({\boldsymbol{x}}_i,{y}_i\right)\hskip0.62em \mathrm{where}\hskip0.22em |\mathcal{D}|=n,{\boldsymbol{x}}_i\in {\mathrm{\mathbb{R}}}^m,{y}_i\in \mathrm{\mathbb{R}}\right\}, $$

a tree ensemble model that predicts an output as a sum of $ K $ functions has the form:

(11) $$ {\hat{y}}_i=\unicode{x03D5} \left({\boldsymbol{x}}_i\right)=\sum \limits_{k=1}^K{f}_k\left({\boldsymbol{x}}_i\right),\mathrm{where}\hskip0.22em {f}_k\in \mathcal{F}. $$

The set $ \mathcal{F}=\left\{f(x)={w}_{q(x)}\hskip0.42em \mathrm{where}\hskip0.22em q:{\mathrm{\mathbb{R}}}^m\to T,w\in {\mathrm{\mathbb{R}}}^T\right\} $ represents the set of regression trees, where $ q $ is the structure of each tree that maps the input $ x $ to the corresponding leaf index, and $ T $ is the number of leaves in the tree. Each function $ {f}_k $ corresponds to a specific tree structure $ q $ and leaf weights $ w $ . To identify the set of functions for the model, the following regularized objective function is minimized:

(12) $$ \mathcal{L}=\sum \limits_il\left({\hat{y}}_i,{y}_i\right)+\sum \limits_k\Omega \left({f}_k\right), $$

where $ \Omega (f)=\gamma T+\frac{1}{2}\lambda {\left\Vert w\right\Vert}^2 $ is the regularization term described in detail later, which represents a penalty for the complexity of the model. In equation (12), $ l $ represents a differentiable convex function that measures the disagreement between the predictions $ {\hat{y}}_i $ and the data points $ {y}_i $ . It is important to note that the tree model in equation (12) includes functions as parameters and cannot be minimized directly in Euclidean space using traditional optimization methods. Instead, the model is trained in an additive form, using multiple iterations.

More formally, if $ {\hat{y}}_i^{(t)} $ is the prediction of the $ i $ -th function of the model at iteration $ t $ , this approach called boosting, ensures that each subsequent tree in the ensemble (called a forest) corrects the classification error of the previous tree. This algorithm minimizes the following regularized objective function:

(13) $$ {\mathcal{L}}^{(t)}=\sum \limits_{i=1}^n\ell \left({y}_i,{\hat{y}}_i^{\left(t-1\right)}+{f}_t\left({\boldsymbol{x}}_i\right)\right)+\Omega \left({f}_t\right), $$

where $ \ell $ is the differentiable convex loss function that represents the difference between the prediction $ {\hat{y}}_i^{(t)} $ and the real data $ {y}_i $ , and $ \Omega \left({f}_t\right) $ is the regularization term defined in the same way as in (12), having the form:

(14) $$ \Omega \left({f}_t\right)=\gamma T+\frac{1}{2}\lambda {\left\Vert w\right\Vert}^2=\gamma T+\frac{1}{2}\lambda \sum \limits_{j=1}^T{w}_j^2. $$

In equation (14), $ T $ represents the number of leaves in the tree, while $ {w}_j $ denotes the weight assigned to the $ j $ -th leaf. The parameter $ \gamma $ controls the complexity penalty for each leaf, and $ \lambda $ is the $ {L}^2 $ regularization term that helps prevent overfitting by constraining the leaf weights. This objective function promotes the creation of simpler decision trees or, in the context of this model, weaker learners.

For optimization, the XGBoost model utilizes both first-order and second-order gradient information. The model updates the objective function by approximating it with a second-order Taylor expansion:

(15) $$ {\mathcal{L}}^{(t)}\approx \sum \limits_{i=1}^n\left[{g}_i{f}_t\left({x}_i\right)+\frac{1}{2}{h}_i{f}_t{\left({x}_i\right)}^2\right]+\Omega \left({f}_t\right), $$

where $ {g}_i=\frac{\partial \ell \left({y}_i,{\hat{y}}_i\right)}{\partial {\hat{y}}_i} $ and $ {h}_i=\frac{\partial^2\ell \left({y}_i,{\hat{y}}_i\right)}{\partial {\hat{y}}_i^2} $ are the first and second derivatives of the loss function concerning the prediction $ {\hat{y}}_i $ .

The XGBoost model is typically used for classification or regression tasks, where each row in the dataset is treated independently of the others. However, time series data have a fundamentally different structure, where the value at any given time point typically depends on values from both earlier and future time points. Therefore, to apply XGBoost to time series forecasting, the original time series must be preprocessed into a new dataset with a specific format that captures these temporal dependencies and makes the XGBoost framework applicable.

In the presented analysis, this has been implemented by creating “lagged features,” that is, dataset columns that contain the values of the target variable from previous time points (Sharma, Reference Sharma2024). For example, to predict the number of cases for the given month, the model uses the number of cases from one month before, two months before, and so on, up to a fixed, predefined number of previous months for a given month as input features. Each row in this new dataset then corresponds to one point in time and includes its own set of lagged values of a fixed length, that is, a fixed number of previous points. This transformation allows the model to learn patterns over time, such as trends or cycles, based on recent history for each time point. In the presented model, the lag value was set to six previous time points (months). This choice was based on the available length of the time series, which was limited to 2 years, and the assumption that a typical wave of incidence develops within this timeframe, providing the model with sufficient context to identify meaningful patterns and make accurate forecasts.

The other limitation of decision tree-based models such as XGBoost is that they can only predict values within the range of the training data. This can cause prediction instabilities for nonstationary time series data, that is, those where the average level or variability changes over time, because the model may fail to capture upward or downward trends and instead produce flat (plateaued) predictions. To address this issue and to avoid such instabilities, relative changes in values, that is, how much a value increased or decreased compared to the previous time point, were modeled instead of using the absolute values (Sharma, Reference Sharma2024). This helps to stabilize the data, and such an approach is very similar to how ARIMA-type models utilize input data differencing transformations to ensure the stationarity assumption for the modeled time series. More specifically, for the analyzed data, the model was trained to predict whether the value would increase or decrease and by how much, rather than predict the exact absolute count.

This lag calculation approach was applied to both the target variable (e.g., incidence and mortality increment counts) and the external predictors from Google Trends. In this setup, for each keyword considered as a predictor, six additional columns corresponding to its lagged values were added to the resulting data frame to create the desired set of predictive features.

This study used the R implementation of the XGBoost model from the xgboost package. The training parameter values were: nrounds set to 100, max_depth set to 10, and eta set to 0.1.