A Comprehensive Statistical Analysis of COVID-19 Trends: Global and U.S. Insights through ARIMA, Regression, and Spatial Models

3.1 Theoretical Basis of the ARIMA model

The ARIMA model is a widely utilized statistical method for analyzing and forecasting time series data. Its general form for an ARIMA model of order (p, d, q) is given by the following equation:

Where:

• ∇^d = (1–B)^d is the differencing operator, with B representing the backshift operator (Box, Jenkins, Reinsel, & Ljung, 2015).

• ϕ(B) = 1 – ϕ₁B – …– ϕ_pB^p is the autoregressive (AR) coefficient polynomial (Box et al., 2015).

• θ(B) = 1 – θ₁B – …– θ_qB^q is the moving average (MA) coefficient polynomial (Box et al., 2015).

• ε_t denotes the white noise error terms, which satisfy the following properties: E (ε_t) = 0, , and E (ε_tε_s) = 0 for s ≠ t (Box et al., 2015).

• Furthermore, E (x_sε_t) = 0 for all s < t, ensuring that the noise terms are uncorrelated with past values of the series (Hamilton, 1994).

The ARIMA model is composed of three primary components:

1. Autoregressive (AR): The AR part represents the dependence between an observation and several lagged observations (Box et al., 2015).

2. Integrated (I): The integrated component represents the differencing needed to make the time series stationary (Box et al., 2015).

3. Moving Average (MA): The MA part models the dependency between an observation and residual errors from a moving average model applied to lagged observations (Box et al., 2015).

The ARIMA modeling process begins with a stationarity test, commonly performed using the Augmented Dickey-Fuller (ADF) test (Dickey & Fuller, 1979). If the time series is non-stationary, transformations such as differencing or logarithmic scaling are applied to achieve stationarity (Hamilton, 1994). Next, model identification involves determining the order of the model, specifically the values of p and q, which are the autoregressive and moving average terms, respectively. This is usually done by analyzing the autocorrelation function (ACF) and partial autocorrelation function (PACF) plots (Box et al., 2015). The differencing order d is selected based on the transformations applied during the stationarity phase. After the model is identified, the parameters ϕ_i and θ_j are estimated, typically using maximum likelihood estimation (MLE) (Hamilton, 1994). Model validation then follows, where statistical tests such as the Ljung-Box test are used to ensure the residuals behave like white noise, indicating that the model has adequately captured the time series’ structure (Ljung & Box, 1978). Model selection is based on criteria like the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC), with the model having the lowest criterion value generally preferred (Akaike, 1974). Finally, once validated, the model is used for forecasting future values of the time series (Box et al., 2015).

The process is visually represented in Figure 1, which illustrates the steps of fitting an ARIMA model to a time series, starting from stationarity checks to forecasting.

• Download figure
• Open in new tab

FIGURE 1

ARIMA Model Construction Flow Chart

3.2 Rolling Window Cross-Validation and Comparison with auto.arima

In this study, rolling window cross-validation was used to evaluate the performance of ARIMA models for time series forecasting. The primary goal was to identify the optimal ARIMA model parameters by minimizing the Root Mean Squared Error (RMSE) and to compare the results with those obtained from the automated model selection function, auto.arima (Hyndman & Athanasopoulos, 2018).

Rolling window cross-validation is a method specifically designed for time series data, which preserves the temporal order of the observations. In each iteration, a model is trained on a fixed-length window of historical data and then validated on the subsequent observation. This approach ensures that the evaluation reflects real-world forecasting scenarios, where future values are predicted based on past data (Bergmeir & Benítez, 2012). For each ARIMA model evaluated, the one-step-ahead forecast errors were calculated, and RMSE was used as the evaluation metric. RMSE is given by: where x_i are the actual values and are the predicted values. Lower RMSE values indicate better model performance (Chai & Draxler, 2014). A grid search was performed over various combinations of autoregressive (p) and moving average (q) parameters, with the differencing order (d) fixed at 1. This process was parallelized to efficiently explore the parameter space (Hyndman & Athanasopoulos, 2018).

The reason for selecting RMSE as the evaluation metric lies in its ability to quantify the average prediction error while giving more weight to larger errors (Willmott & Matsuura, 2005). RMSE is particularly useful in contexts where significant deviations in the forecast could have a considerable impact, as it highlights larger discrepancies between predicted and actual values more than other metrics, such as Mean Absolute Error (MAE). Additionally, RMSE is measured in the same units as the original data, making the results easier to interpret in practical applications.

To compare the performance of manual model selection with automated methods, the auto.arima function was used. auto.arima automatically selects the best ARIMA model by optimizing information criteria such as the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) (Hyndman & Khandakar, 2008). While auto.arima quickly identifies a suitable model based on a global fit to the entire dataset, rolling window cross-validation provides a more dynamic evaluation by assessing the model’s performance across different periods in the data (Tashman, 2000). This allows for a comparison of how well the automated model selection aligns with the results of cross-validated manual tuning.

By visualizing the RMSE values across different parameter combinations, we were able to compare the best-performing model identified by rolling window cross-validation with the model chosen by auto.arima. This comparison provided insights into the trade-offs between automated and manual ARIMA model selection in time series forecasting.

3.4 Theoretical Basis of the ARIMAX model

To improve the accuracy of time series forecasting, we utilized the AutoRegressive Integrated Moving Average with eXogenous variables(ARIMAX) model, which integrates external factors into the standard ARIMA framework. This extension allows the model to account for influences beyond the inherent patterns in the target time series (Hyndman & Athanasopoulos, 2018). For this study, we explored whether including an exogenous variable, such as vaccination rates, would improve forecast accuracy compared to the ARIMA model, which relies solely on the historical values of the time series.

The ARIMAX model expands the ARIMA model by introducing exogenous regressors believed to impact the dependent variable. Mathematically, the ARIMAX model can be expressed as: where y_t represents the value of the dependent variable at time t, ϕ_i are the autoregressive coefficients, θ_j are the moving average coefficients, ϵ_t is the error term, and X_t–k represents the exogenous variable lagged by k periods (Pankratz, 1991). The exogenous variables are incorporated to capture additional influences on the time series that are not explained by the series’ own historical values (Box et al., 2015).

The ARIMAX model fitting was conducted using an automated selection of ARIMA parameters (p, d, q), while incorporating the exogenous variable. To evaluate the performance of the ARIMAX model, we compared it with the ARIMA model using standard evaluation metrics such as Akaike Information Criterion (AIC), Root Mean Squared Error (RMSE), and Mean Absolute Error (MAE). Since these metrics were introduced earlier, they will not be repeated here (Chai & Draxler, 2014; Hyndman & Athanasopoulos, 2018).

Model forecasts were generated for a holdout period to assess prediction accuracy. The inclusion of exogenous variables in the ARIMAX model allowed us to examine whether external factors could improve the forecasting performance, providing a more comprehensive understanding of the dynamics affecting the time series. The comparison between ARIMA and ARIMAX models highlighted both the benefits and limitations of incorporating external information into the forecasting process (Pankratz, 1991).

3.5 Analysis of the Impact of Vaccination on New COVID-19 Cases

To analyze the relationship between vaccination rates and the number of new COVID-19 cases, several statistical methods were employed, including Granger causality testing, segmented regression, and regression discontinuity design (RDD). These methods help in understanding both the temporal relationships and potential causal effects of vaccination on the incidence of new cases (Box et al., 2015; Hyndman & Athanasopoulos, 2018).

3.6 Regression Analysis of COVID-19 Infection Rates and Determinants

3.7 Spatial Autocorrelation and Hotspot Analysis of COVID-19 Cases

In this study, spatial analysis techniques were applied to examine the distribution of COVID-19 infection rates across different regions. The analysis involved calculating Moran’s I for global spatial autocorrelation and performing the Getis-Ord Gi* statistic to identify local hotspots and coldspots. Visualization of the results was conducted using traditional red-blue color schemes, effectively highlighting areas with significant spatial clustering of high or low infection rates (Anselin, 1995; Ord & Getis, 1995).

4.1 Short-Term Forecasting with ARIMA Models and Anomaly Detection

To evaluate the short-term predictive performance of ARIMA models on COVID-19 case counts, forecasts were generated for four distinct periods using training data from prior months. The predictive performance was then evaluated against the actual observed data.

The first forecast, covering September 27 to December 27, 2020, was based on data from January 5, 2020, to September 27, 2020. As shown in Figure 2a, the forecast generally follows the actual case trajectory, though some deviations appear towards the end of the period, suggesting the model’s limitations in capturing sudden changes. The ACF and PACF plots (Figures 2b and 2c) indicate some residual autocorrelation, highlighting potential areas for model improvement. The Box-Ljung test returned a p-value of 0.3746, indicating that there is no significant residual autocorrelation.

• Download figure
• Open in new tab

FIGURE 2

ARIMA Model Analysis in the U.S. Over Different Time Periods

For the second period, December 27, 2020, to March 28, 2021, the model was extended to include data up to December 27, 2020. Figure 2d illustrates a closer alignment between the predicted and actual cases, with only minor deviations. The ACF and PACF plots (Figures 2e and 2f) further support the model’s adequacy, though some residual correlations persisted. The Box-Ljung test for this period yielded a p-value of 0.6327, further indicating that residual autocorrelation is not a concern.

The third forecast, covering March 28 to June 27, 2021, used data up to March 28, 2021. As depicted in Figure 2g, the model closely tracks the actual case counts, demonstrating robust predictive capability. The corresponding ACF and PACF plots (Figures 2h and 2i) show that the model effectively captures the data’s temporal structure, though the Box-Ljung test resulted in a p-value of 0.07281, suggesting that there may still be some minor residual autocorrelation.

Finally, the model was applied to forecast cases from September 26 to December 26, 2021, using training data from January 3, 2021, to September 26, 2021. As seen in Figure 2j, the model continues to perform well, closely aligning with the observed case counts during the forecast period. The ACF and PACF plots (Figures 2k and 2l) indicate that the model has successfully captured the underlying patterns, with the Box-Ljung test returning a p-value of 0.2876, indicating minimal residual autocorrelation.

Throughout these periods, the ARIMA model demonstrated consistent predictive accuracy, although the residual autocorrelation evident in the ACF and PACF plots suggests areas where further refinement could enhance the model’s performance.

Among the four distinct periods analyzed using ARIMA models, the forecast for the interval from September 27 to December 27, 2020, stands out as the least accurate. Several factors may contribute to this discrepancy between the forecasted and actual data during this period.

One possible reason for the inaccuracy could be the inherent limitations of the ARIMA model, which is a linear model designed to predict future values based on past data. This model may struggle to capture sudden, nonlinear changes or external shocks that occur during the forecast period. ARIMA models assume a level of stationarity in the data, and if the underlying time series experiences structural breaks or sudden shifts, the model’s predictions can become less reliable.

Another contributing factor could be the presence of significant outliers or unexpected spikes in COVID-19 cases during the forecast period. Such anomalies might arise due to the emergence of new virus variants, changes in public health policies, or sudden shifts in public behavior, all of which can lead to rapid increases in case numbers that the model, trained on historical data, might not adequately predict.

To investigate this hypothesis, an outlier detection analysis was conducted on the data from January 5, 2020, to December 27, 2020. The detected outliers, shown in Table 2 and illustrated in Figure C2, highlight several key dates where significant anomalies were observed. These anomalies correspond to periods of sharp increases in case counts, suggesting that the forecast discrepancies could be attributed to these sudden and unexpected changes.

As shown in Table 2, significant outliers were detected on November 8, November 15, and December 13, 2020, where the cumulative cases sharply increased. These dates likely correspond to specific events or conditions that triggered a surge in cases, such as the spread of more transmissible variants or changes in testing or reporting practices.

To explore potential anomalies in COVID-19 case trends, an outlier detection analysis was performed on the cumulative COVID-19 case data for both the United States and globally. The goal was to identify points in time where the actual case numbers significantly deviated from the expected trend, potentially indicating periods when new variants emerged and spread.

The results of the outlier detection for the United States from January 5, 2020, to December 27, 2020, are shown in Figure 3a, with the detected outliers summarized in Table C1a. Notably, several of these dates align with the emergence of significant COVID-19 variants, such as the Omicron variant (B.1.1.529), first identified in November 2021 in South Africa and Botswana (Organization, 2021). Other variants, like BQ.1 and BQ.1.1, began spreading rapidly in late 2022, contributing to the increased number of cases that may have led to forecast inaccuracies (Organization, 2022a). The figure 3b shows Time Series Plot with Detected Outliers for USA.

• Download figure
• Open in new tab

FIGURE 3

Side-by-side comparison of outliers detected and time series plot for USA

Further analysis was conducted on a global scale, with the results depicted in Figure C1. The corresponding dates and case numbers for these global outliers are listed in Table C1b. Similar to the U.S. data, these global outliers coincide with key dates when variants such as XBB, CH.1.1, and BF.7 were identified and began spreading across multiple regions, leading to significant surges in cases (Organization, 2022b). These variants, first noted in late 2022 and early 2023, had a notable impact on regions such as Asia and Europe, leading to substantial deviations from the predicted trends (Organization, 2023).

The detected outliers in both the U.S. and global datasets underscore the significant impact of COVID-19 variant emergence on the spread of the virus. Although the Alpha variant (B.1.1.7) and Gamma variant (P.1) were not explicitly captured by the outlier detection process due to their emergence towards the end of 2020, the general trend observed in Figure 3a (the U.S. outlier detection graph) does indicate a notable spike in cases during this period (Organization, 2021). This spike aligns with the timeline of Alpha and Gamma variants becoming prevalent, suggesting that the rapid spread of these variants, coupled with their increased transmissibility and potential for immune escape, likely contributed to the observed surge in case numbers. Consequently, almost all significant spikes in the data, as identified through outlier detection, correspond closely with the emergence of a new COVID-19 variant.

The dates of these outliers, as presented in Tables 1 and 2, closely match the timelines for variant identification and global spread. For example, the spikes in U.S. case counts from November 2021 to February 2022 correspond with the rise of Omicron and its subvariants (Organization, 2022a). Similarly, the global spikes identified in late 2021 and throughout 2022 align with the spread of Omicron and its subvariants, reinforcing the idea that these variants significantly impacted the accuracy of forecasted versus actual case numbers.

4.2 Regional COVID-19 Forecasting Across Continents

The ARIMA model was employed to forecast COVID-19 cases across different continents, including Asia, Europe, Africa, the Americas, and South America. The forecasting approach for South America involved excluding Canada, the United States, and Mexico from the Americas dataset to better approximate the COVID-19 trajectory specific to South America. The forecast results for each continent are illustrated in the figure D3, with predictions covering the period from January 2020 to early 2021.

In Asia, the ARIMA model’s predictions aligned closely with the observed data, capturing the overall upward trend in COVID-19 cases. The prediction intervals encompassed the actual case numbers, indicating the model’s robustness in this region.

In Europe, the ARIMA model’s predictions were less accurate, as the predicted cases significantly deviated from the actual observed data. This discrepancy is particularly evident towards the end of 2020, where a sharp and sudden increase in cases occurred, which the ARIMA model failed to predict effectively. Upon reviewing the forecast and considering our earlier anomaly detection for the U.S. and global data, alongside reports from the WHO regarding the emergence of variants, it is plausible to attribute this rapid rise to the Alpha variant (B.1.1.7). This variant, first detected in September 2020 in the United Kingdom, was the first to be classified by the WHO as a ”variant of concern.” The Alpha variant’s high transmissibility likely contributed to the significant increase in cases, which was not fully captured by the ARIMA model, highlighting the challenges of forecasting during periods of rapid epidemiological change.

For Africa, the model demonstrated a good fit with the actual data, though, like in Europe, the rapid rise in cases near the year’s end pushed the limits of the prediction interval. In the Americas, the ARIMA model performed strongly, with predictions closely matching the steep increase in case numbers. This region, which experienced one of the most significant surges in cases, was effectively modeled, with predictions falling within the confidence intervals.

In South America, after excluding the northern countries, the ARIMA model continued to perform well. The predicted cases remained within reasonable bounds compared to the observed data, similar to the other continents.

Across all regions, the Box-Ljung test p-values were significantly above the 0.05 threshold, indicating no significant auto-correlation in the residuals. This suggests that the ARIMA models successfully captured the temporal patterns of COVID-19 case progression in each region. The occasional underestimations, particularly during periods of rapid increases in cases, underscore the challenges posed by the pandemic’s dynamic nature and the emergence of new variants, which may not have been fully accounted for in models trained on earlier data. Nonetheless, the overall performance of the ARIMA models across these diverse regions was robust, providing valuable insights into the spread of COVID-19 during the forecasted periods.

4.3 Rolling Window Cross-Validation and Comparison with auto.arima

In our previous ARIMA forecasting efforts, we relied on the auto.arima function to automatically select the ARIMA model parameters p, d, and q. The auto.arima function is based on minimizing the Akaike Information Criterion (AIC), which balances model fit and complexity by penalizing excessive parameters. This approach offers several advantages, including speed, automation, and often, reasonably good results. However, its reliance on AIC may not always yield the most accurate forecasts, especially when dealing with nonstationary time series or when the model is used for long-term predictions.

To explore whether other parameter selection methods could improve forecast accuracy, we applied a rolling window cross-validation technique to optimize the p and q parameters, while keeping d fixed as determined by auto.arima. The underlying reason for fixing d lies in its role in differencing, which is generally determined by the data’s trend and seasonality structure-a topic well-supported by statistical theory. For example, once the time series is made stationary through differencing, the differencing order d should remain constant to maintain this stationarity, regardless of changes in p and q.

In our analysis, we selected a period where the ARIMA predictions notably diverged from the actual data, such as the case with COVID-19 cases in the USA and Europe from January 5, 2020, to December 27, 2020. The divergence was primarily due to the unforeseen surge in cases caused by the emergence of new variants, highlighting the limitations of traditional ARIMA models in capturing such abrupt changes.

Using the rolling window cross-validation approach, we evaluated different combinations of p and q based on the Root Mean Squared Error (RMSE) metric. This method, which assesses the out-of-sample performance of models across multiple training windows, is particularly valuable when forecasting non-stationary time series with varying patterns over time. Table 3 summarizes the RMSE values for the USA’s ARIMA model with parameters selected via cross-validation versus those obtained using auto.arima, while Figure 4a provides a heatmap visualizing the RMSE across different p and q combinations.

• Download figure
• Open in new tab

FIGURE 4

Heatmap of RMSE and forecast comparison of USA

Note: The blue points represent the observed data, and the red line represents the fitted regression discontinuity model. The vertical black line indicates the intervention date (December 13, 2020).

As depicted in Figure 4a, the RMSE heatmap clearly shows that the cross-validated ARIMA parameters (p = 2, q = 2) outperform the auto.arima parameters (p = 1, q = 0) in terms of RMSE. The heatmap provides a comprehensive view of how different combinations of p and q impact the forecast accuracy, with lower RMSE values indicating better performance.

Furthermore, Figure 4b compares the forecasted COVID-19 cases in the USA using the auto.arima model and the cross-validated ARIMA model. Although both models exhibit significant deviations from the actual data due to the unexpected surge in cases, the cross-validated model’s predictions are slightly closer to the actual values than those of auto.arima, suggesting that the cross-validation approach can yield more accurate forecasts under certain conditions.

For the ARIMA model applied to Europe, a similar approach was employed as described in the previous. Table E2 presents the RMSE values for the European ARIMA model, comparing the parameters selected via cross-validation with those determined by auto.arima. The RMSE heatmap shown in Figure E4a visualizes the performance across different combinations of p and q. Figure E4b illustrates the comparison of the forecasted COVID-19 cases in Europe using the ARIMA model with parameters selected by cross-validation against those obtained by auto.arima. As observed, the forecast line derived from the cross-validated model is significantly closer to the actual data than the one produced by auto.arima, although both models still show notable deviations from the actual trajectory. These findings closely mirror the conclusions drawn from the cross-validation results for the ARIMA model applied to the USA, further reinforcing the potential advantages of cross-validation in parameter selection for ARIMA models in the context of highly volatile and non-stationary time series data.

4.4 Effect of Vaccination on New COVID-19 Cases

From December 2020, the global effort to vaccinate against COVID-19 began, raising the critical question of whether the vaccination campaign effectively reduced the number of new COVID-19 cases. To address this question, several statistical methods were employed, including the Granger Causality Test, segmented regression analysis, the Chow Test, and Regression Discontinuity Design (RDD).

The analysis began with the Granger causality test, aimed at determining whether the number of people vaccinated could serve as a predictor for future new COVID-19 cases, while controlling for past cases. Two models were compared: the first model included lags of both new cases and vaccinations, while the second model included only lags of new cases. The results did not indicate a significant causal relationship between vaccination and a reduction in new cases, as the F-statistic was 0.24 with a p-value of 0.9746, suggesting that the inclusion of vaccination data did not improve the predictive power of the model within the lags tested. Specifically, with a lag of 7 (equivalent to 49 days), the Granger causality test showed no significant effect of vaccination on new cases within this period. As shown in table4.

To further investigate the potential impact of vaccination on the trend in COVID-19 cases, segmented regression analysis was employed, introducing a breakpoint at the onset of the vaccination campaign. The regression model accounted for time, an indicator for the postintervention period, and the interaction between time and the post-intervention phase. The analysis revealed that while the overall trend in new cases showed a positive slope (β = 18987, p = 0.02136), the interaction term (time_post_intervention) was negative and significant (β = –24115, p = 0.00445), indicating a deceleration in the growth rate of new cases following the intervention. Despite this, the post-intervention indicator itself was not statistically significant (p = 0.31000), which aligns with the results of the Granger causality test, further suggesting that the immediate effect of vaccination on reducing new cases was not significant. However, the negative and significant interaction term implies that vaccination had a significant long-term impact on reducing new cases, indicating a beneficial effect over time. Figure 5a illustrates the segmented regression results, showing how the predicted number of cases diverges from the actual cases over time. As shown in Table F3, the segmented regression results demonstrate these trends clearly.

• Download figure
• Open in new tab

FIGURE 5

Comparison of Segmented Regression of COVID-19 Cases and Regression Discontinuity Plot

Additionally, the Chow Test was conducted to formally assess the presence of a structural break at the intervention point. The test provided strong evidence of a structural change, with a p-value of 6.437e-06, indicating that the introduction of the vaccination campaign significantly altered the underlying relationship between time and new cases. This result corroborates the findings from the segmented regression analysis, suggesting that vaccination led to a structural shift in the trend of new cases. Finally, a Regression Discontinuity Design (RDD) was applied to further validate these findings. The RDD analysis, focusing on the sharp change in the trend of new cases at the point of intervention, yielded a conventional coefficient estimate of 76662.154 with a non-significant p-value (p = 0.636). This suggests that while there may have been an observable shift in the trend of new cases at the intervention point, it was not statistically significant at conventional levels. As shown in Table F4, the regression discontinuity results further support the conclusion that the immediate impact of vaccination was not statistically significant. Figure 5b provides a visualization of the RDD results, highlighting the discontinuity at the intervention point. The non-significant result from the RDD analysis is consistent with the findings from both the Granger causality test and the segmented regression, indicating that the immediate impact of vaccination was not significant.

4.5 Forecast on COVID-19 cases using ARIMAX with vaccine as an exogenous variable

Building upon the results in the previous section, where it was demonstrated that vaccination had a significant long-term impact on the reduction of new COVID-19 cases, a logical extension was to incorporate the number of vaccinations as an exogenous variable in ARIMAX models. The hypothesis was that the inclusion of this variable might improve the accuracy of forecasts compared to a standard ARIMA model, which does not account for such external factors.

To test this hypothesis, we utilized our dataset where vaccination began on December 13, 2020. Given the evidence that the impact of vaccination is more pronounced over the long term, the first training period selected for the model was from January 5, 2020, to June 27, 2021, approximately six months after the start of vaccination. This period was used to predict the cumulative number of cases for the subsequent three months. Subsequently, the training data span was gradually extended in two further scenarios:

1. January 5, 2020 – June 27, 2021

2. January 5, 2020 – December 26, 2021

3. January 5, 2020 – September 25, 2022

Figures 6, G5, and G6 respectively compare the ARIMA and ARIMAX model predictions for these three periods, and tables 5, G5, G6 compare the metrics of ARIMA and ARIMAX models across three time periods, including AIC, RMSE, and MAE.

• Download figure
• Open in new tab

FIGURE 6

ARIMAX forecast for period 1

From these results, it can be observed that the ARIMAX model sometimes produces forecasts closer to the actual data than the ARIMA model, as indicated by lower RMSE and MAE values in certain scenarios. However, there are also cases where the ARIMAX model deviates more from the actual data, resulting in higher RMSE values. Interestingly, the Akaike Information Criterion (AIC) does not always correlate directly with RMSE and MAE improvements. For example, in the third period, despite the ARIMAX model providing predictions closer to the actual data (as reflected by lower RMSE and MAE), its AIC is higher than that of the ARIMA model. This higher AIC reflects the trade-off between model complexity and goodness-of-fit inherent in the AIC calculation.

4.6 Multivariate Regression Analysis of Global COVID-19 Infection Rates

To investigate the factors influencing COVID-19 infection rates across different countries, we began by hypothesizing that more developed countries, with their advanced medical infrastructure and higher availability of healthcare resources, might exhibit lower infection rates. However, an examination of the top 10 countries by infection rate as of December 31, 2023 (Figure H7), revealed that many of the countries with the highest infection rates are, in fact, highly developed. For example, countries such as Luxembourg, Denmark, and Austria, all of which are considered highly developed, show among the highest infection rates, challenging the initial hypothesis.

To explore this relationship further, we conducted a linear regression analysis between the COVID-19 infection rate and GDP per capita as a proxy for a country’s level of development. The scatterplot with the regression line is displayed in Figure 7. The regression output showed a significant positive relationship between GDP per capita and infection rate, with the coefficient for GDP per capita being positive and highly significant (p < 2e-16), indicating that higher GDP per capita is associated with higher infection rates. The regression model, however, had a relatively low R-squared value of 0.4763, suggesting that while GDP per capita is a significant predictor, it explains only about 47.63% of the variance in infection rates.

• Download figure
• Open in new tab

FIGURE 7

The scatterplot with regression line

In addition to the regression analysis, we calculated three correlation metrics to further understand the relationship between GDP per capita and infection rate. The Pearson correlation coefficient was 0.690161, indicating a moderately strong positive linear relationship between the two variables. The Spearman rank correlation coefficient was higher, at 0.8593242, suggesting a strong monotonic relationship. Finally, the Maximal Information Coefficient (MIC) was 0.7256325, pointing to a strong association that may capture nonlinear relationships between GDP per capita and infection rate. These correlations further support the conclusion that higher GDP per capita is associated with increased infection rates, although other factors likely contribute to the remaining unexplained variance.

Given the relatively low R-squared, it was clear that other factors beyond GDP per capita might influence infection rates. Therefore, we expanded the model to include additional variables that could plausibly affect infection rates: Human Development Index (HDI), Gini coefficient (a measure of income inequality), per capita health expenditure, the number of hospital beds per 1,000 people, and population density. The resulting multivariate regression model incorporated both main effects and interaction terms between these variables.

The multivariate regression results, which included interaction terms, revealed a more complex relationship between the predictors and the infection rate. While GDP per capita continued to have a significant effect (p = 0.006544), other variables like health expenditure and certain interaction terms also emerged as significant predictors. For example, the interaction between GDP per capita and HDI (p = 0.006381) and between GDP per capita and the Gini coefficient (p = 0.029676) were both significant, indicating that the effect of GDP per capita on infection rates is moderated by a country’s HDI and income inequality. Additionally, the interaction between HDI and health expenditure (p = 0.000724) was significant, suggesting that the combined effect of these two factors significantly influences infection rates. The detailed results of the regression analysis, including estimates, standard errors, t-values, and p-values, are provided in Table I7.

Despite these findings, the model’s R-squared increased substantially to 0.8179, indicating that approximately 81.79% of the variance in infection rates can be explained by the expanded set of predictors and their interactions. However, the residual plots (Figure 8a) reveal some potential issues with model fit, including non-constant variance (heteroscedasticity) and some deviation from normality in the residuals, as indicated by the Q-Q plot.

• Download figure
• Open in new tab

FIGURE 8

Comparison of Residual Plots and Scatterplot Matrix

The scatterplot matrix (Figure 8b) and coefficient plot (Figure J8) further underscore the complexity of the relationships among the predictors. The scatterplot matrix highlights the correlations between variables, with some expected relationships, such as the positive correlation between GDP per capita and HDI (corr = 0.729) and a negative correlation between GDP per capita and the Gini coefficient (corr = -0.330). The coefficient plot shows the magnitude and direction of the effects, with GDP per capita, health expenditure, and certain interaction terms having the most pronounced impacts on infection rates.

4.6.1 Addressing Multicollinearity in the Regression Model

The initial multivariate regression model that included interaction terms significantly improved the model’s explanatory power, as indicated by a substantial increase in the R² value. However, this complexity introduced severe multicollinearity into the model, as evidenced by the extremely high Variance Inflation Factor (VIF) values. For instance, variables such as GDP per capita, HDI, and health expenditure, along with their interaction terms, exhibited VIF values that were in the tens of thousands, indicating that multicollinearity was indeed a significant problem. Multicollinearity can render the regression coefficients unstable and their interpretation challenging, which necessitated a more rigorous approach to model simplification and stabilization.

To mitigate these issues, we employed stepwise regression as an initial step to reduce model complexity by removing less significant predictors. The simplified model retained key variables and interactions while excluding those that contributed less to the model’s overall explanatory power. The resulting model, while more manageable, continued to exhibit notable multicollinearity, with several VIF values remaining high, albeit reduced from their initial levels.

Given the persistence of multicollinearity, we explored alternative methods to further address this issue. Both Partial Least Squares (PLS) and Principal Component Regression (PCR) were considered, as these techniques are specifically designed to handle multicollinearity by transforming the predictor variables into a set of uncorrelated components.

We applied PLS and PCR to the dataset, each method aiming to reduce the dimensionality of the predictor variables while maximizing the explained variance in the response variable, infection rate. The PLS analysis, as shown in Table 7, was particularly effective in this context. The model explained 67.22% of the variance in infection rates using five components, which was identified as the optimal number of components through cross-validation (Figure 9). Beyond five components, the model’s Mean Squared Error of Prediction (MSEP) began to increase, suggesting that additional components might introduce noise rather than improve predictive accuracy.

View this table:

• View inline
• View popup
• Download powerpoint

TABLE 6 Linear Regression Results: Infection Rate vs. GDP per Capita.

View this table:

• View inline
• View popup
• Download powerpoint

TABLE 7 PLS Results: Variance Explained by Number of Components

• Download figure
• Open in new tab

FIGURE 9

Cross-Validation Mean Squared Error Plot of the PLS model

In addition to the error plot, Table K9 provides detailed cross-validation results for each model. This table presents the MSEP for different numbers of components, highlighting how the error decreases as the number of components increases up to five, and then rises with the inclusion of additional components. These results reinforce the findings illustrated in Figure 9, where five components were found to be optimal.

The component loadings from the PLS model, visualized in the heatmap (Figure M9), highlight the contribution of each variable to the principal components. Variables such as GDP per capita, HDI, and health expenditure had substantial loadings on the first few components, indicating their significant influence on the model. However, complex interactions, such as those between GDP per capita and population density or health expenditure and population density, also played critical roles in later components.

In addition to the heatmap, the detailed PLS loadings are presented in Table L10. This table provides specific loading values for each variable across the first five components, further illustrating the contributions and interactions among variables in shaping the principal components.

The PLS model, through component reduction, provided a more stable set of coefficients, as shown by the reduced VIF values and improved interpretability of the regression coefficients. The final regression coefficients derived from the PLS model (Table L11) illustrate the direct and interaction effects of the predictor variables on infection rates, offering a clearer understanding of the underlying relationships.

PCR, as shown in Table K8, on the other hand, showed similar results but with slightly lower explained variance for the same number of components. While PCR effectively reduced multicollinearity, the trade-off was in the form of lower predictive power compared to PLS. Specifically, PCR explained 59.43% of the variance in infection rates with six components, which increased to 75.81% with 17 components, but did not surpass the PLS model in terms of overall performance.

4.7 Spatial Autocorrelation and Hotspot Analysis of COVID-19 Cases in U.S. States

We performed a spatial autocorrelation and hotspot analysis of COVID-19 cases in U.S. states. Spatial autocorrelation measures the degree to which COVID-19 cases are spatially clustered, while hotspot analysis identifies regions with significantly high or low case counts.

Using the most recent COVID-19 case data across U.S. states, we calculated the number of COVID-19 cases per 100,000 people to normalize for population differences. Figure 10a illustrates the spatial distribution of these normalized case counts across the United States. States with higher case rates per 100,000 people are represented in shades of red, while those with lower rates are in shades of blue. Notably, Alaska and several states in the South exhibit particularly high case rates.

• Download figure
• Open in new tab

FIGURE 10

(a) COVID-19 Cases per 100,000 People by State, (b) Getis-Ord Gi* Hotspot Analysis of COVID-19 Cases per 100,000 People

In addition to the spatial visualization, Table N12 provides the detailed rankings of the top 10 and bottom 10 states in terms of COVID-19 cases per 100,000 people. Alaska ranks highest with 40,576.16 cases per 100,000 people, followed by Rhode Island and Kentucky. Conversely, New York, Maryland, and Oregon report the lowest case rates, with New York at the bottom with 18,251.51 cases per 100,000 people.

To further understand the spatial pattern, we conducted a Moran’s I test, a measure of spatial autocorrelation. The results revealed a Moran’s I statistic of 0.1578 with a p-value of 0.03169, indicating a significant positive spatial autocorrelation. This suggests that states with high COVID-19 case rates tend to be geographically clustered rather than randomly distributed.

To identify specific clusters of high or low case rates, we applied the Getis-Ord Gi* statistic, which provides a measure of local spatial clustering. As shown in Figure 10b, the hotspot analysis identified several ”hot” and ”cold” spots. Southern states such as Arkansas, Georgia, and Mississippi, as well as parts of Texas and Ohio, were identified as hotspots (with high Gi* values), indicating significant clustering of high case rates. Conversely, states like Alaska, Delaware, New Hampshire, and Vermont were identified as cold spots, suggesting significant clustering of low case rates.

In addition, Table N13 lists the top 10 hotspot and coldspot states along with their corresponding Gi* values. For instance, Arkansas has the highest Gi* value of 1.000414, marking it as a significant hotspot, while Alaska has the lowest Gi* value of -1.237323, making it a prominent coldspot. These detailed Gi* values provide a quantitative basis for understanding the spatial clustering patterns observed in the map.

For Alaska, given that it does not share borders with any other state, we designated Washington as its sole neighbor. Despite this adjustment, Alaska, which has the highest infection rate among all states, was classified as a cold spot in the Getis-Ord Gi* Hotspot analysis. This counterintuitive result may be due to the isolation of Alaska, where the lack of adjacent states diminishes the influence of its high case rate on surrounding areas. Consequently, its high infection rate is not part of a broader regional trend, leading the Gi* statistic to categorize it as a cold spot rather than a hotspot. This highlights the importance of considering geographic and relational context in spatial analyses, particularly for isolated regions.