Impact of COVID-19 pandemic in the Brazilian maternal mortality ratio: A comparative analysis of Neural Networks Autoregression, Holt-Winters exponential smoothing, and Autoregressive Integrated Moving Average models

Background and objectives The acute respiratory infection caused by severe acute respiratory syndrome coronavirus disease (COVID-19) has resulted in increased mortality among pregnant, puerperal, and neonates. Brazil has the highest number of maternal deaths and a distressing fatality rate of 7.2%, more than double the country’s current mortality rate of 2.8%. This study investigates the impact of the COVID-19 pandemic on the Brazilian Maternal Mortality Ratio (BMMR) and forecasts the BMMR up to 2025. Methods To assess the impact of the COVID-19 pandemic on the BMMR, we employed Holt-Winters, Autoregressive Integrated Moving Average (ARIMA), and Neural Networks Autoregression (NNA). We utilized a retrospective time series spanning twenty-five years (1996–2021) to forecast the BMMR under both a COVID-19 pandemic scenario and a controlled COVID-19 scenario. Results Brazil consistently exhibited high maternal mortality values (mean BMMR [1996–2019] = 57.99 ±6.34/100,000 live births) according to World Health Organization criteria. The country experienced its highest mortality peak in the historical BMMR series in the second quarter of 2021 (197.75/100,000 live births), representing a more than 200% increase compared to the previous period. Holt-Winter and ARIMA models demonstrated better agreement with prediction results beyond the sample data, although NNA provided a better fit to previous data. Conclusions Our study revealed an increase in BMMR and its temporal correlation with COVID-19 incidence. Additionally, it showed that Holt-Winter and ARIMA models can be employed for BMMR forecasting with lower errors. This information can assist governments and public health agencies in making timely and informed decisions.


Introduction
The acute respiratory infection caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has led to increased mortality among pregnant, puerperal, and neonates [1].Maternal deaths related to coronavirus disease (COVID- 19) in Brazil have exceeded global figures due to factors such as clinical characteristics, social conditions, and barriers to adequate care [2].According to the COVID-19 Observatory Newsletter, Brazil has the highest number of maternal deaths and a distressing fatality rate of 7.2%, more than double the country's current mortality rate of 2.8% [3].A study in Brazil reported 9,609 cases of pregnant and puerperal individuals with SARS-CoV-2 from December 2019 to August 2020, of which 4,230 (44.0%) were confirmed as COVID-19 positive.Among these cases, 533 (12.60%) of pregnant and puerperal individuals died, with 354 (64.0%) of these deaths attributed to COVID-19 [4].A study conducted in the Democratic Republic of Congo revealed that most COVID-19-related deaths were due to delayed care and advanced age, with common causes, including hemorrhages, uterine ruptures, infections, and dystocias [5].
Access to quality prenatal, birthing, and postpartum care was compromised for pregnant and puerperal women during the SARS-CoV-2 pandemic [4].In 2021, pregnant and puerperal women faced unfavorable outcomes and severe clinical presentations related to COVID-19 due to the presence of high-risk SRAG symptoms and a greater risk of intubation and intensive care unit (ICU) admission compared to other groups [6].Additionally, pregnant women encountered challenges in accessing ventilators and intensive care during the pandemic [3].
Contingency measures focused on maternal health are essential to enhance access to prenatal services and intensive care for pregnant, puerperal women, and neonates [7].The World Health Organization (WHO) emphasizes the need to reduce maternal mortality by 2030, as outlined in the Sustainable Development Goals (SDGs) [8].However, barriers to accessing specialized care services and appropriate monitoring of obstetric complications persist in Brazilian hospitals and primary care settings [1].
Maternal mortality values over the years can be analyzed as a time series, and forecasting methods utilize time series data to predict future trends.Two widely recognized time series models for forecasting COVID-19 deaths and recoveries are the Holt-Winters and Auto-Regressive Integrated Moving Average (ARIMA) models [9].The application of novel forecasting methods, including Neural Network Autoregression (NNA), has gained popularity. NNA mimics the structure and functionality of the human brain, using "neurons" to address complex challenges [10].NNAs have been successfully employed in various fields to predict clinical events, such as peripherally inserted central catheter-related thrombosis in breast cancer patients [11].
Hence, the application of forecasting methods offers a quantitative framework for scientists to explore hypotheses related to potential disease mechanisms or causes of mortality.To the best of our knowledge, this is the first published comparative study, including Holt-Winters, ARIMA, and NNA models, on Brazilian Maternal Mortality Ratios (BMMR).These models are valuable for assessing the impact of interventions, optimizing the effectiveness of control strategies, and generating short-and long-term forecasts.Therefore, this study aims to analyze the relationship between BMMR and the COVID-19 pandemic and forecast BMMR trends up to 2025.

Background Brazilian Maternal Mortality Ratios
The BMMR series (Fig 1) reveals that Brazil consistently had high maternal mortality values, considering WHO parameters (mean BMMR  = 57.99±6.34/100,000 live births).Notably, during the third quarter of 2009, when analyzing the peak of maternal deaths, the BMMR reached 81.00/100,000 live births, marking a 39.67% increase compared to the mean.This increase was associated with the H1N1 pandemic, which is known to have adverse effects on maternal and perinatal outcomes [12].Studies have shown that pregnancy and the postpartum period were statistically associated with deaths from pandemic influenza [13].
In August 2020, the Pan American Health Organization (PAHO) highlighted the elevated risk of severe forms of COVID-19 in pregnant women and recommended that countries intensify efforts related to women's and newborns' health to ensure access to prenatal and disease prevention services [14].However, the Brazilian health system faced challenges in identifying cases of pregnant and postpartum women with COVID-19, resulting in delays in the immunization of pregnant and postpartum women, which only began in July 2021 [15].Consequently, the country experienced the highest mortality ratio in the historical BMMR series during the second quarter of 2021 (169.00/100,000live births), representing an approximately 191% increase compared to the previous period and a notable contrast with global reports on maternal deaths during the COVID-19 pandemic [2].
Ecological spatial analysis research conducted in Brazil revealed that COVID-19 deaths among pregnant and postpartum women exhibited a heterogeneous geographic distribution, with distinct spatial clusters primarily located in rural areas [16].Most maternal deaths occurred during the postpartum period, as investigations using data from the Brazilian Ministry of Health's Epidemiological Surveillance system identified 40,640 women of reproductive age hospitalized for COVID-19.Among them, 3,372 were pregnant, and 794 were postpartum.Postpartum women had a worse prognosis compared to pregnant women, with higher rates of ICU admissions, invasive ventilatory support, and fatalities [17].
A study comparing the first and second waves of COVID-19 in pregnant and postpartum women in Brazil reported 377 maternal deaths in 2020 and 804 in 2021.Furthermore, the COVID-19 mortality rate doubled compared to the previous year.In the obstetric group, the rate increased from 7.7% to 15.4%, while in the non-obstetric group, it increased from 13.9% to 22.9% from 2020 to 2021.Among women with comorbidities, the mortality rate increased by 1.7 times (from 13.3% to 23.3%) and by 1.4 times (from 22.8% to 31.4%) [18].A retrospective time series study using the Holt-Winters exponential smoothing additive model revealed devastating consequences for maternal mortality during the COVID-19 pandemic in the state of Bahia, Brazil [19].

Background forecast methods
Several studies have indicated that NNA, Holt-Winters exponential smoothing, and ARIMA are established approaches for forecasting COVID-19 case counts and deaths.For example, Lynch and Gore (2021) assessed seven forecasting methods with varying look-back and forecast lengths at the country, health district, and state levels [20].Papastefanopoulos et al. (2020) employed statistical and machine learning-inspired time series methods to estimate the percentage of active cases relative to the total population in the ten countries with the highest active cases [21].Barrı ´a-Sandoval et al. (2021) conducted a comparative analysis of ARIMA models, Exponential Smoothing, State Space models, the Bayesian approach, and the GLARMA model to predict confirmed COVID-19 cases and deaths in Chile [22].Gecili et al. (2021) evaluated four forecast models for predicting COVID-19 confirmed cases, deaths, and recoveries in the USA and Italy based on daily reported data [9].

Maternal mortality ratio ¼ maternal death number of live births
x 100; 000 The WHO has established parameters to categorize the maternal mortality ratio as low, medium, high, or very high: Low: up to 20/100,000 live births; Average: from 20 to 49/100,000 live births; High: from 50 to 149/100,000 live births; and Very high: > 150/100,000 live births.In this study, the BMMR time series is presented by quarter.

Study design
We conducted an observational time-series cross-sectional data analysis of quarterly BMMR data.To investigate the impact of COVID-19, we used available pre-pandemic data (2020-2019) to forecast the BMMR up to 2025 (COVID-19-free scenario).Subsequently, we compared the forecast values (2020-2023( 1)) with the BMMR reported by the official source.Additionally, we used all historical data (1996-2021) to forecast the BMMR up to 2025 (pandemic scenario).To avoid issues of underfitting or overfitting due to improper training of the dataset, the models underwent cross-validation, with a validation window spanning three years (2016-2019).
Package, data set, and reproducibility.All models were constructed using the "forecast" package version 8.15, which incorporates Forecasting Functions for Time Series and Linear Models [10].All datasets are available in the S1 File.Computer code used for the analyses is provided in the S2 File.All statistical analyses were performed on a computer equipped with two Intel 1 Xeon 1 Silver 4214 Processors (13.75 M Cache, 2.20 GHz), 96 GB of RAM, 24 cores, and 64-bit Windows 10 Enterprise, using the RStudio statistical software, Version 1.4.1717.

Outliers identification
Outliers are observations that significantly differ from the majority of observations in a time series.They may represent errors or unusual events, such as pandemic occurrences.Outliers can alter the dynamics of a time series either temporarily or permanently.These changes are typically non-systematic and cannot be captured by standard time series models.Detecting outliers is crucial because they impact model selection, parameter estimation, and, consequently, forecasts.Outliers were identified using the"tsoutliers()" function in the "forecast" package for R (version 8.15), through an iterative outlier detection and adjustment process using Seasonal Decomposition of Time Series by Loess (with iterate = 2 and lambda = NULL) [24].

Fitting Holt-Winters model
The HW smoothing technique, also known as triple exponential smoothing, was employed to model the BMMR data.This method offers two variations, differing primarily in their treatment of seasonal components: additive and multiplicative.The choice between these variations depends on the nature of seasonal variations in the data.Additive models are preferable when seasonal variations remain roughly constant throughout the series, while multiplicative models are suitable when seasonal variations change proportionally with the level of the series.In this study, the multiplicative method (Eq 2), the seasonal component is expressed in relative terms (percentages), and the series is seasonally adjusted by dividing by the seasonal component [10].The multiplicative models provide better results; thus, they were used in this study.
where ℓ t represents the level; b t is the trend; s t represents the seasonal component; m represents the period of the seasonality (quarterly data), and h represents the forecast.Optimal smoothing weights were computed as follows: α = 0.1571, β* = 0.2434, and γ = 0.1979.The adequacy of the Holt-Winters model was assessed by examining the randomness of the model residuals (S1 Fig) .The ACF plot of the residuals showed that all autocorrelations fell within threshold limits, indicating that the residuals behaved like white noise.The histogram of the residuals displayed a symmetrical distribution.Furthermore, the Ljung-Box test (Q* = 5.1668, df = 8, p-value = 0.7396) confirmed that the residuals did not significantly deviate from a white noise series.

Fitting ARIMA model
ARIMA is a statistical analysis model that utilizes time series data to enhance understanding of the dataset or predict future trends.ARIMA models assume a linear correlation between timeseries values and aim to use these linear dependencies in observations to extract local patterns while eliminating high-frequency noise from the data [25].
It combines differencing with autoregression and a moving average model, as described by Eq 3.
where y 0 t is the different time series.The "predictors" include both lagged y t values and errors (ε t ).The model is represented by ARIMA (p, d, q), where: p represents the autoregressive process, the influence of the previous value of the variable on the value under consideration; d is related to the number of differentiations to induce stationarity; and q, is related to the influence of the noise generated in the previous values [10].A seasonal ARIMA model is formed by incorporating additional seasonal terms into ARIMA and can be represented as ARIMA (p, d, q) (P, D, Q) m , where uppercase and lowercase notations denote seasonal and non-seasonal components of the model, and m represents the seasonal period.
Although the ARIMA model is valuable and powerful in time series analysis, selecting the appropriate orders for its components can be somewhat challenging.For our purposes focused on the usefulness of the forecast, we chose to determine the orders of the ARIMA model automatically, using the "auto.arima"function in the "forecast" package for R. The function performs a search on the possible model within the constraints of order provided according to the selected metric (such as AIC, AICc, or BIC value).This choice facilitates both the implementation of the model, and its updates as new data becomes available, allowing greater applicability of the model [9].
For optimizing ARIMA models, we employed the Box and Jenkins methodology in a fourstep iterative process, as follows: model identification, parameter estimation, diagnostic checking, and prediction [26].
1) Model identification.The ARIMA model is suitable for stationary time series data, as supported by the literature [27,28].We verified stationarity using a unit root test, specifically the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test, calculated using the "urca" package [29].The obtained p-value of 0.1728 for the time series significantly exceeded the 1% critical value, indicating non-stationarity and the need for differencing [29].After one differencing, the pvalue became 0.0234, confirming that the differenced data are stationary (S2 Fig) .2) Parameter's estimation.To determine the structure of correlation between time lags of the differenced data, we assessed autocorrelation function (ACF) and partial autocorrelation (PACF) plots (S2 Fig) .After evaluating several models, we selected the ARIMA (0,1,3)(2,0,0) [4] model due to its lowest Akaike's Information Criteria (AIC) and Bayesian Information Criteria (BIC), making it suitable for fitting the BMMR (1996-2021).
3) Model diagnostics (goodness of fit).We evaluated the adequacy of the ARIMA (0,1,3) (2,0,0) [4] model by examining the randomness of the residuals (S3 Fig) .The ACF plot of the residuals indicated that all autocorrelations were within threshold limits, suggesting that the residuals behaved like white noise.Additionally, the histogram displayed symmetry.The Ljung-Box test (Q* = 3.7432, df = 3, p-value = 0.2906) confirmed that the residuals were indistinguishable from a white noise series.
Fitting NNA model.NNA is a forecasting method based on mathematical models of the brain, enabling the modeling of complex non-linear relationships between the response variable and its predictors.A key characteristic of NNA is its ability to self-learn without prior knowledge of the intricate non-linear relationships between input and output variables.
In the context of time series data, lagged values of the time series can serve as inputs.The predictors, or inputs, form the bottom layer, while the intermediate layers include "hidden neurons" (optional), and the forecasts or outputs constitute the top layer.In a multilayer feedforward network, each node layer receives inputs from the previous layers and combines them using weighted linear combinations.For instance, the inputs into hidden neuron j are linearly combined, as shown in Eq 4, and this is subsequently modified using a non-linear function like a sigmoid, as depicted in Eq 5.
Where n is the number of neurons in the input layer, parameters b j and w j are "learned" (or estimated) from the data during the training of the model.The notation NNA (p, k, n) indicates there are p-lagged inputs and k nodes in the hidden layer [10].
Various parameters were tested and evaluated (S1 Table ).Considering all combinations of these parameters, 2,400 different NNA models were trained.Among these models, the NNA (8,6,1) [6] demonstrated superior performance by yielding the lowest Root Mean Squared Error (RMSE) [30].In this particular model, the last eight observations are utilized as non-seasonal predictors, six as seasonal lag predictors, and one neuron is placed in the hidden layer.
The suitability of the NNA (8,6,1) [6] model was assessed based on the randomness of model residuals (S4 Fig) .The ACF plot of the residuals indicates that all autocorrelations fall within the threshold limits, suggesting that the residuals exhibit behavior similar to white noise.Additionally, the histogram displays symmetry.The Ljung-Box test (Q* = 6.6877, df = 8, pvalue = 0.5707) confirms that the residuals cannot be distinguished from a white noise series.
Models' comparison.When comparing forecast methods applied to a single time series or to multiple time series with the same units, the mean absolute error (MAE, Eq 6) is the most commonly used metric because it is easy to understand and compute.On the other hand, the root mean squared error (RMSE, Eq 7) minimizes the impact of bias and measures the dispersion of prediction errors, indicating model stability, despite being somewhat challenging to interpret [10].A forecasting method that minimizes MAE tends to produce forecasts closer to the median, while minimizing RMSE results in forecasts closer to the mean.RMSE and MAE were employed as performance indicators for comparing the fitted models: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Additionally, p-values were computed using a two-sample Student's t-Test on data sets from two independent populations (reported versus predicted) with unequal variances.H 0 : true difference between reported and forecast BMMR is equal to zero (significance level of ɑ = 0.02).

Ethical aspects
All data utilized in this study were obtained from publicly available sources through the Office of Public Health.Consequently, the study did not undergo ethical review.

Pandemic events
The sudden increase in BMMR during the second quarter of 2009 (81.00/100,000 live births) is associated with the H1N1-Influenza (Flu) pandemic, while the peak (197.75/100,000live births) in the second quarter of 2021 is related to the COVID-19 pandemic (2020-2021).Predicting the occurrence of the next pandemic is unfeasible due to the random nature of such events.Therefore, for modeling purposes, these events were treated as outliers [31].Four outliers were identified in the BMMR time series for the years 1996 to 2021 (Fig 1).

Holt-Winters, ARIMA, and NNA forecasts
The Holt-Winters forecasts for controlled and pandemic COVID-19 scenarios are presented in Fig 2 .In a controlled COVID-19 scenario, the BMMR exhibits a decreasing trend, reaching a ratio of 49.90/100,000 live births, with an 80% confidence interval (CI) of [28.07, 71.74] in 2025.However, in a pandemic scenario, the BMMR shows an increasing trend, reaching a ratio of 119.47/100,000 live births, with an 80% CI of [102.15, 136.80] in 2025.
The ARIMA forecast for controlled and pandemic COVID-19 scenarios is presented in Fig 3 .In a controlled COVID-19 scenario, the BMMR exhibits a stable trend, reaching a ratio of 58.62/100,000 live births, with an 80% CI of [44.78, 72.46] in 2025.However, in a pandemic scenario, the BMMR shows an increasing trend, reaching a ratio of 96.29/100,000 live births, with an 80% CI of [85. 25, 107.34] in 2025.
The NNA forecast for controlled and pandemic COVID-19 scenarios is presented in Fig 4 .In a controlled COVID-19 scenario, the BMMR exhibits a stable trend, reaching a ratio of 65.47/100,000 live births, with an 80% CI of [61.08, 68.96] in 2025.However, in a pandemic scenario, the BMMR shows a slightly increasing trend, reaching a ratio of 85.60/100,000 live births, with an 80% CI of [81.10, 92.75] in 2025.
Models performance.The performance of the models was evaluated based on a test set with twelve samples.The RMSE for NNA (4.5701/100,000 live births) was lower than that for Holt-Winters (7.1026/100,000 live births) and ARIMA (4.7613/100,000 live births).In contrast, the MAE was the lowest for ARIMA (3.6750/100,000 live births), compared to Holt-Winters (6.0426/100,000 live births) and NNA (3.7579/100,000 live births).Additionally, NNA provided a small value for Cross-Validation (CV) and a better-adjusted R 2 , while ARIMA provided small values for Akaike's Information Criterion (AIC) and Schwarz's Bayesian Information Criterion (BIC), as shown in Table 1.The ARIMA and NNA models appear to be more suitable at this stage.

Increase in mortality associated with COVID-19
The BMMR reported during the COVID-19 pandemic (2020-2021) was significantly higher than previous (1996-2019; p < 0.01) and forecasted (2022-2025; p < 0.01) values, as shown in Table 2.This underscores the impact of the pandemic on maternal mortality.Notably, the BMMR returned to historical levels after the second quarter of 2022, coinciding with the Ministry of Health of Brazil declaring the end of the Public Health Emergency of National Importance due to Covid-19 on April 22, 2022.Considering the forecast results beyond the sample period (2022 and beyond), the Holt-Winters model demonstrated its superiority in forecasting the BMMR time series.Although ARIMA and NNA models provided an excellent fit to the The results indicate that even if the COVID-19 pandemic stabilizes, Brazil is unlikely to achieve the maternal mortality target outlined in the SDGs in the coming years.Therefore, forecast studies play a crucial role in informing more effective actions for preventing maternal deaths.Improvements in living conditions, access to quality healthcare, and comprehensive care are essential to preventing maternal deaths related to pregnancy, childbirth, and the postpartum period.The historical data also highlight that events such as H1N1 and COVID-19 expose existing weaknesses in maternal and child health in the country.It is important to note that accurate data preprocessing for time series models can be challenging, often requiring a trial-and-error approach.Additionally, NNA methods are computationally intensive and demand powerful computing resources and time.Lastly, a limitation of the study is the use of secondary data published in Brazil, and the data on maternal deaths in 2022 being preliminary.
Limitations.Threats to internal validity relate to factors that affect the dependent variables without the researchers' knowledge [20].It is possible that some implementation flaws affected subsequent modeling or data analysis.To minimize this risk, the algorithms in our source code were (1) built in established libraries, (2) underwent internal reviews, and (3) are publicly accessible, along with the data and results.Threats to external validity pertain to the generalizability of the analysis results.The results obtained are specific to Brazil, limited to BMMR reported by DATASUS (official database), and based on preliminary data for 2022.These results may not be immediately applicable to (1) different causes of death, (2) other datasets, (3) different time periods, or (4) different geographic regions.Additionally, the selected forecasting methods assume that public health policies remain unchanged during the forecast periods.To enhance model performance, forecast analyses should be regularly updated as new data becomes available, thereby improving the accuracy of these models.

Conclusion
This study highlights the increase in BMMR and its temporal association with the incidence of Influenza A H1N1 and the COVID-19 pandemic in Brazil.Moreover, it demonstrates that ARIMA and Holt-Winters models were the most effective forecasting methods for modeling maternal mortality when compared to neural networks.These forecasts contribute to the formulation of more effective public policies for the maternal and child population and the expansion of access to healthcare services.Furthermore, these models are readily accessible, as they are implemented in an R package, and can be easily applied by analysts.The fitting process for each model only takes a few minutes, making them conducive to regular updates as new data becomes available.

Fig 2 .
Fig 2. Forecast of the quarterly Brazil Maternal Mortality Ratio up to 2025 using the Holt-Winters multiplicative method.: Black line represents the quarterly Brazil Maternal Mortality Ratio (1996-2021) reported by the official source; red line represents values provided by the Holt-Winters model; blue line represents forecast values in a controlled COVID-19 scenario; purple line represents forecast values in a pandemic COVID-19 scenario, and the shaded area indicates the 80% confidence interval limits.https://doi.org/10.1371/journal.pone.0296064.g002

Fig 3 .
Fig 3. Forecast of the quarterly Brazil Maternal Mortality Ratio up to 2025 using ARIMA (0,1,3) (2, 0, 0) [4].The black line represents the quarterly Brazil Maternal Mortality Ratio (1996-2021) reported by the official source; the red line represents values provided by the ARIMA model; the blue line represents forecast values in a controlled COVID-19 scenario; the purple line represents forecast values in a pandemic COVID-19 scenario, and the shaded area indicates the 80% confidence interval limits.https://doi.org/10.1371/journal.pone.0296064.g003

Fig 4 .
Fig 4. Forecast of the quarterly Brazil Maternal Mortality Ratio up to 2025 using the NNA (8,6,1) [6] model.The black line represents the quarterly Brazil Maternal Mortality Ratio (1996-2021) reported by the official source; the red line represents values provided by the NNA model; the blue line represents forecast values in a controlled COVID-19 scenario; the purple line represents forecast values in a pandemic COVID-19 scenario, and the shaded area indicates the 80% confidence interval limits.https://doi.org/10.1371/journal.pone.0296064.g004