Predicting COVID-19 Hospital Stays with Kolmogorov–Gabor Polynomials: Charting the Future of Care

Abstract

Optimal allocation of ward beds is crucial given the respiratory nature of COVID-19, which necessitates urgent hospitalization for certain patients. Several governments have leveraged technology to mitigate the pandemic’s adverse impacts. Based on clinical and demographic variables assessed upon admission, this study predicts the length of stay (LOS) for COVID-19 patients in hospitals. The Kolmogorov–Gabor polynomial (a.k.a., Volterra functional series) was trained using regularized least squares and validated on a dataset of 1600 COVID-19 patients admitted to Khorshid Hospital in the central province of Iran, and the five-fold internal cross-validated results were presented. The Volterra method provides flexibility, interactions among variables, and robustness. The most important features of the LOS prediction system were inflammatory markers, bicarbonate (HCO₃), and fever—the adj. R² and Concordance Correlation Coefficients were 0.81 [95% CI: 0.79–0.84] and 0.94 [0.93–0.95], respectively. The estimation bias was not statistically significant (p-value = 0.777; paired-sample t-test). The system was further analyzed to predict “normal” LOS ≤ 7 days versus “prolonged” LOS > 7 days groups. It showed excellent balanced diagnostic accuracy and agreement rate. However, temporal and spatial validation must be considered to generalize the model. This contribution is hoped to pave the way for hospitals and healthcare providers to manage their resources better.

1. Introduction

The fast spread of the SARS-CoV-2 coronavirus has placed immense strain on healthcare systems across the globe. As infected individuals surged, the demand for hospital admissions grew accordingly [1]. Past outbreaks have demonstrated that limited bed capacity and hospital resources significantly contribute to higher infectious disease mortality rates [2]. Hence, guidelines for prioritizing patients and determining who should be admitted for essential care are instrumental in addressing resource limitations. Neglecting this could jeopardize the lives of COVID-19 patients [3].

Nine to eleven percent of COVID-19 hospitalizations required enhanced life-support interventions [4]. However, the ICU faced challenges accommodating these needs due to limited beds and shortages in monitoring equipment, life-sustaining machinery, and skilled staff crucial for top-tier care [5]. In a study encompassing 183 nations in 2021, Sen-Crowe et al. [2] reported that high-income areas registered the highest average ICU beds at 12.79 and 402.32 hospital beds for every 100,000 individuals. On the other hand, regions with upper-middle income showed dominance in average acute-care beds, numbering 424.75 per 100,000 inhabitants. This is not the case for low- and middle-income countries, where the number of ICU beds is often insufficient, and the equipment is often old and poorly serviced. This number was five beds per one million people in Africa [6].

Challenges in managing hospital capacity throughout this pandemic spanned various phases, including testing, treatment, and preparation for future patients. As a result, there is a pressing need to accurately predict and prioritize patients based on the likelihood of their condition escalating in severity. It is part of the pandemic preparedness action plan.

In predicting hospital length of stay (LOS), the overarching challenges introduced by COVID-19 cannot be overlooked. The pandemic has significantly strained hospital capacities, potentially altering standard care pathways and discharge protocols. Furthermore, heightened fatigue [7], burnout [8], and stress among healthcare professionals [9], a byproduct of the ongoing crisis, may also have indirect implications for the duration of patient stays. These combined factors elucidate the multifaceted dynamics influencing hospital operations during these unprecedented times.

Numerous studies have explored predicting hospital resource needs for COVID-19 patients. Many of these investigations have leveraged machine learning (ML). ML has established itself as an invaluable tool in the medical realm, adept at sifting through and synthesizing vast amounts of data to discern intricate patterns. Most health-related challenges nowadays rely heavily on ML to disentangle the complexities inherent in large-scale data, facilitating informed healthcare decisions.

During outbreaks like COVID-19, forecasting the imminent demand for medical resources such as beds and nasal oxygen support becomes crucial. In this context, ML methodologies have proven invaluable [10,11]. For instance, researchers from London designed an ML algorithm that outperformed clinical experts in predicting COVID-19 patient mortality [12]. Another ML study successfully predicted which COVID-19 patients would transition into a severe respiratory phase with a 70–80% accuracy rate [13].

Furthermore, an AI-based tool named “ambient warning and response evaluation” has been employed to refine ICU clinical settings. This tool significantly enhanced timely and accurate decision-making, leading to a 37% reduction in LOS [14].

LOS estimation remains crucial for efficient healthcare management, offering insights into patient health trajectories, resource allocation, and the quality-of-care delivery. The state-of-the-art research listed encompasses a myriad of methodologies and priorities, thereby revealing both the advancements and the persisting gaps in LOS prediction.

Nemati et al. (2020) [15] utilized a global dataset and focused on a limited set of five variables, primarily age and sex, to estimate LOS. Their approach, which involved stagewise gradient boosting, did not venture into comprehensive features but mainly centered around symptoms onset date and symptoms. Given the minimalistic input feature set, this focus might limit its applicability in varied clinical settings.

Working in a tertiary care hospital in China, Hong et al. (2020) [16] used logistic regression with a set of 37 variables, including lymphocyte and neutrophil count, heart rate, and procalcitonin levels, D-dimer, and partial thrombin time. Their dataset was also relatively small, including 75 patients considering the number of predictors, and reached an AUC of 0.85 to classify prolonged (>14 days) versus normal (≤14 days) hospital LOS. Their work lacks internal and external validation, indicating potential overfitting risks.

Ebinger et al. (2021) [17] embarked on an extensive exploration of 966 patients with 353 variables of electronic health records (EHRs) to classify patients based on extended stays (i.e., LOS > 8 vs. LOS ≤ 8 days) in the Cedars-Sinai Medical Center. Forty-two machine learning models were used as ensemble models of 12 base classifiers (including Elastic-net and random forest). Such models were trained using the first 1, 2, and 3 days of hospital admission. Advanced Average (AVG) Blender for the day 3 model outperformed the others. Age, Interleukin 6, blood urea nitrogen level, and oxygen flow rate were among the selected features. The best model had an area under the ROC curve (AUC) of 0.82 and a precision of 67%.

Usher et al. (2021) [18] analyzed data from 36 hospitals across Minnesota, Wisconsin, and the Dakotas. Using 20 variables, which included diverse features such as age, critical illness, mechanical ventilator (MV) application, and oxygen requirement, their approach adopted the random forest method, considering it as the best model. The classification output was the LOS ≤ 5 days (reference), LOS between 5 and 10 days, LOS between 10 and 15 days, and LOS > 15 days. With five-fold cross-validation, they achieved an AUC of 0.89, highlighting the potential of integrating diverse input features for LOS category prediction.

Mahboub et al. (2021) [19] at Rashid Hospital in Dubai took a distinct route by incorporating treatments as input features and variables such as urea, platelets, and D-dimer. Utilizing decision trees on a dataset of 2017 patients, they achieved a coefficient of determination (R²) of 0.5, suggesting the relevance of treatment variables in predicting LOS.

Liuzzi et al. (2022) [20] from the Fondazione Don Carlo Gnocchi Living COVID-19 Registry in Italy incorporated a comprehensive set of 829 variables, with a focus on 55 primary variables spanning across admission clinical scales, symptoms, and therapies. Their method, employing sequential convolutional neural networks, was validated with repeated five-fold cross-validation, resulting in a median absolute deviation of 2.7 days.

Orooji et al. (2022) [21] in Iran, with data from 1225 patients, utilized 53 variables and emphasized 20 key features such as age, creatinine, and lymphocyte/neutrophil count. They applied statistical feature selection combined with multi-layer perceptron and 12 training algorithms, reaching a root-mean-square error (RMSE) of 1.6213 days.

In 2022, Alabbad et al. [22] from King Fahad University Hospital in Saudi Arabia classified ICU LOS into nine categories using 43 variables. The synthetic minority oversampling technique (SMOTE) was used to balance the class distribution. Their best model employed random forest, and they also explored gradient boosting and extreme gradient boosting. With three-fold cross-validation, their model boasted a positive predictive value (PPV) of 94%, indicating high precision in prediction.

Alam et al. (2023) [23] from Prince Sultan Hospital in Riyadh incorporated 89 variables, including laboratory data, X-ray results, clinical data, and treatments, to classify LOS into seven categories. Their model utilized the Tab Transformer and achieved impressive results, with an F1 score of 93% for discharged patients. The SMOTE-N oversampling technique was also noted to balance the class distribution.

Zhang et al. (2023) [24] analyzed 83 variables, including immunotherapy and heparin, to predict LOS for 384 patients at Zhengzhou University Hospital. Using the least absolute shrinkage and selection operator (LASSO) and linear regression, they explained 30% of LOS variability (R² = 0.30). Missing data were managed with imputations, and results were verified via bootstrap validation.

Overall, while significant strides have been made in predicting LOS through diverse methodologies, ranging from classical regression models to neural networks, gaps in validation and comprehensive feature inclusion, conditioning on future events (e.g., therapies) resulting in selection bias, incorporating time-dependent predictors (e.g., treatments) as time-fixed, leading to immortal-time bias [25], and balancing the dataset, resulting in biased performance indices [26] remained a consistent challenge. Moreover, sample size insufficiency based on the number of input features [27] was the other problem of some methods proposed in the literature. Further research was required to enhance model generalizability across varied clinical settings.

Our research aimed to employ multivariable analysis and the Kolmogorov–Gabor polynomial to craft a predictive model. This model aimed to precisely forecast the LOS of COVID-19 patients in a nationally representative sample of the pediatric population in the Middle East and North Africa (MENA) based on their demographic and clinical data upon hospital admission.

Our primary model was designed to predict the continuous LOS. We evaluated and presented performance metrics for this continuous prediction. Additionally, we derived a binary representation of LOS from the predicted and actual data, categorizing it as either prolonged (LOS > 7 days) or normal (LOS ≤ 7 days). This binary classification’s performance was also examined. We adopted this approach to accommodate the existing literature’s categorical and continuous LOS representations. Our primary focus remained the continuous prediction model, which can seamlessly be converted to a binary prediction through straightforward post-processing.

2. Materials and Methods

2.1. Data Source

In this retrospective study, we examined the clinical records of N = 1600 confirmed COVID-19 cases with complete information from Isfahan, situated in the center of Iran, from 6 March to 7 May 2020. These patients were admitted to Khorshid Hospital, which caters to the vast metropolitan area of Isfahan, home to over 15 million residents. Given that this hospital functioned as the primary referral center for critical COVID-19 cases during this period, our study exclusively focused on the patients admitted to the hospital. Patients with a positive RT-PCR test confirming SARS-CoV-2 infection or confirmed chest computed tomography (CT) results were enrolled in this study.

All participants’ LOS was calculated from their initial hospital ward or ICU admission until discharge. It is noteworthy to mention that this LOS represents the first recorded admission. Comprehensive information regarding the study design and the methods used to register variables can be found in our Khorshid COVID Cohort (KCC) study [28]. The data gathered included demographic details such as age and sex, pertinent dates including COVID-19 diagnosis and hospital or ICU admission, and the patient’s most recent known clinical status.

2.2. Data Description and Pre-Processing

This study extracted and used patients’ records, including non-clinical, clinical, and symptom data. Non-clinical data included sex, age, occupation, education, body mass index, family size, number of family members infected, house area, travel history, duration of symptoms before admission, and history of influenza vaccination. Clinical patient data included principal diagnosis, admission unit, medical history, and comorbidities. Laboratory data included the results of all blood tests performed at patient admission. The latest available laboratory tests included were CBC results, sodium (Na+), potassium (K+), urea, creatinine, alkaline phosphatase (ALP), aspartate transaminase (AST), alanine aminotransferase (ALT), bilirubin, international normalized ratio (INR), lactate dehydrogenase (LDH), C-reactive protein (CRP), ferritin, hemoglobin A1c (HbA1c), D-dimer, erythrocyte sedimentation rate (ESR), and vitamin D. To assess patient health status and identify the required level of care, parameters such as blood pressure, heart rate, and respiratory rate were recorded. Comorbidity categories were evaluated by the Charlson comorbidity index (CCI), which is one of the most commonly used methods to evaluate comorbid factors and predict mortality [29]. It was calculated based on age category, history of myocardial infarction (MI), congestive heart failure (CHF), peripheral vascular disease, history of a cerebrovascular accident or transient ischemic attacks, dementia, chronic obstructive pulmonary disease (COPD), connective tissue disease, peptic ulcer disease, liver disease, diabetes mellitus, hemiplegia, moderate to severe chronic kidney disease (CKD), presence of solid tumor, leukemia, lymphoma, and AIDS, ranging from 0 to 37. These medical conditions were classified by the International Classification of Diseases, Tenth Revision, Clinical Modification (ICD-10-CM) codes that are available in Appendix I Table SI-1 in [30]. The CCI was categorized into five groups: CCI score 0, CCI score 1–2, CCI score 3–4, CCI score 5–6, and CCI score ≥7 [31].

In addition to fever (body degree “up to 39.4 °C”), other symptoms, including fatigue, cough, sore throat, headache, nasal congestion, shortness of breath, severe chest pain, severe muscle pain, vomiting, dry cough, nausea, diarrhea, abdominal pain, muscle and joint pain, general weakness, smell-taste disorder, and dyspnea were identified by the medical interview [32]. Primary composite endpoints (PCEP) were defined as death, the use of mechanical ventilation, or admission to intensive care [33].

2.3. Statistical Data Analysis

Descriptive statistics, including means, frequencies, and proportions, are summarized for the collected data. The disease severity level stratifies summaries. Chi-squared and Fisher exact tests were used whenever appropriate to examine differences among categorical predictors. The endpoint of this study was LOS, which was calculated according to the number of days of hospitalization. The paired-sample t-test was used to identify if the LOS bias was statistically significant [34]. The Bland–Altman plot, (also known as the Tukey mean-difference plot) [35] was provided to analyze the LOS error. Patients were divided into two groups for descriptive analysis, according to the quartile LOS value: ≤7 days as normal and > 7 days as prolonged LOS [36]. Such a cutoff was used in terms of healthcare utilization. We considered p < 0.05 as statistically significant. Predictive modeling was performed offline using MATLAB version 9.6 R2019a (Natick, MA, USA: The MathWorks Inc.), while statistical analysis was performed using IBM SPSS Statistics for Windows, Version 29.0 (Armonk, NY, USA: IBM Corp).

2.4. Predictive Modeling

Volterra functional series, also known as Kolmogorov–Gabor polynomials [37], were used in our study for prediction. The level of interaction was limited to two to reduce the computational complexity and overfitting. The proposed model is provided in Equation (1).

$$\mathrm{y}={\mathrm{a}}_0+\sum _{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{\mathrm{j}=1}^{\mathrm{m}}{\mathrm{a}}_{\mathrm{ij}}{\mathrm{x}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{j}}$$

(1)

where y is the output of the model (LOS), x_i is the i^th input feature (i = 1, …, m), m is the number of features, and the model parameters are a₀ (the offset) and a_ij (two-way interaction coefficients; i, j = 1, …, m). Prior to estimation, the output variable was detrended by subtracting its average. After the model was constructed, this offset was subsequently added back. Since some input features were categorical, one-hot and ordinal encoding were used for nominal and ordinal features, respectively [38], allowing capturing the system’s response for each of multiple generated binary features. Prior to estimating the coefficients, highly correlated (i.e., with an absolute correlation coefficient higher than or equal to 0.8) features and two-way interactions were identified, and some of those were selected to avoid collinearity and multicollinearity [39]. Note that multicollinearity was further reduced by dropping one of the one-hot encoded columns, also known as “dropping one level”. Since Equation (1) contains all two combinations of the input coded features, we used regularized least squares (RLS) [40], with the Euclidean norm penalization, also known as the ridge regression, to estimate the coefficients in the under-determined system:

$${\mathrm{A}}_{\mathrm{R}\mathrm{L}\mathrm{S}}={\left({\mathrm{X}}^{\mathrm{T}}\times \mathrm{X}+\mathsf{\lambda }\mathrm{I}\right)}^{-1}\times {\mathrm{X}}^{\mathrm{T}}\times \mathrm{Y}$$

(2)

where A_RLS are the estimated coefficients in Equation (1), Y is the target LOS vector, X is the data matrix for the selected input features of the training set, and T is the transpose operator. The regularization parameter (λ) was estimated during the cross-validation on the training set [40]. The ridge regression can help reduce the model’s variance and improve its generalization to unseen data and more stable estimates, mitigate the risk of overfitting, manage model complexity, and provide some feature stability [41].

2.5. Model Validation

Five-fold cross-validation was used in our study, and the cross-validated results were provided. The goodness-of-fit of the LOS estimation algorithm was assessed using root-mean-squared-error (RMSE), mean and median absolute error, as well as the coefficient of determination (R²) [42], adjusted R² (adj. R²), and the concordance correlation coefficient (${\mathsf{\rho }}_{\mathrm{c}}$) [43]. For the LoS_i and y_i pairs (i = 1, …, N), such indices were calculated as follows:

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left(\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left(\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathrm{i}}-\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathsf{\mu }}\right)}^2}$$

(3)

where,

$$\mathrm{a}\mathrm{d}\mathrm{j}.\stackrel{}{} {\mathrm{R}}^2=1-\left(\frac{\mathrm{N}-1}{\mathrm{N}-\mathrm{p}-1}\right)\times \left(1-{\mathrm{R}}^2\right)$$

(4)

where p is the number of selected input features of the model.

$${\mathsf{\rho }}_{\mathrm{c}}=\frac{2\times \mathrm{C}\mathrm{o}\mathrm{V}\left(\mathrm{L}\mathrm{o}\mathrm{S},\mathrm{y}\right)}{{\mathsf{\sigma }}_{\mathrm{L}\mathrm{o}\mathrm{S}}^2+{\mathsf{\sigma }}_{\mathrm{y}}^2+{\left(\overline{\mathrm{y}}-\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathsf{\mu }}\right)}^2}$$

(5)

where CoV is the covariance, ${\mathsf{\sigma }}_{\mathrm{L}\mathrm{o}\mathrm{S}}^2=\left(\frac{1}{\mathrm{N}}\right)\times {\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left(\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathrm{i}}-\mathrm{L}\mathrm{o}{\mathrm{S}}_{\mathrm{u}}\right)}^2$ is the variance of the LOS, LOS_u is the mean of LOS, ${\mathsf{\sigma }}_{\mathrm{y}}^2$ is the variance of the predicted LOS, and $\overline{\mathrm{y}}$ is the mean of the predicted LOS.

We further analyzed the binary outcome of the prediction system for the normal (LOS ≤ 7 days) and prolonged LOS (LOS > 7 days) [44]. The performance indices were calculated based on the cross-validated confusion matrix:

• TP (True Positives) = The number of accurately identified prolonged LOS
• TN (True Negatives) = The number of accurately identified normal LOS
• FP (True Positives) = The number of inaccurately identified prolonged LOS
• FN (True Positives) = The number of inaccurately identified normal LOS

The following performance indices were then calculated:

$$\mathrm{S}\mathrm{e}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

(6)

$$\mathrm{S}\mathrm{p}=\frac{\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}}$$

(7)

$$\mathrm{P}\mathrm{P}\mathrm{V}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}}$$

(8)

$$\mathrm{D}\mathrm{O}\mathrm{R}=\frac{\mathrm{T}\mathrm{P}\times \mathrm{T}\mathrm{N}}{\mathrm{F}\mathrm{P}\times \mathrm{F}\mathrm{N}}$$

(9)

$$\mathrm{A}\mathrm{U}\mathrm{C}=\frac{\mathrm{S}\mathrm{e}+\mathrm{S}\mathrm{p}}{2}$$

(10)

$${\mathrm{F}}_1=\frac{2\times \mathrm{T}\mathrm{P}}{2\times \mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}+\mathrm{F}\mathrm{P}}$$

(11)

$$\mathrm{M}\mathrm{C}\mathrm{C}=\frac{\mathrm{T}\mathrm{P}\times \mathrm{T}\mathrm{N}-\mathrm{F}\mathrm{P}\times \mathrm{F}\mathrm{N}}{\sqrt{\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}\right)\times \left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}\right)\times \left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}\right)\times \left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}\right)}}$$

(12)

$$\mathrm{K}\left(\mathrm{C}\right)=\frac{2\times \left(\mathrm{T}\mathrm{P}\times \mathrm{T}\mathrm{N}-\mathrm{F}\mathrm{P}\times \mathrm{F}\mathrm{N}\right)}{\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}\right)\times \left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}\right)\times \left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}\right)\times \left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}\right)}$$

(13)

where Se is the sensitivity, Sp is the specificity, PPV is the positive predictive value, DOR is the diagnostic odds ratio, AUC is the balance diagnostic accuracy (area under the ROC curve), F₁ is the F₁ score, MCC is the Matthews’s correlation coefficient [45,46], and K(C) is the Cohen’s Kappa agreement rate.

Also, the unbiased PPV was calculated based on the sensitivity and specificity of the developed dichotomous LOS model using different prevalence (P) measures of the prolonged LOS in the hospital. PPV is the probability that a patient has prolonged LOS when the dichotomous LOS model results are positive. The related formula was presented in Equation (14). It was estimated using the Bayes’ theorem [26]:

$$\mathrm{u}\mathrm{n}\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{d}\stackrel{}{} \mathrm{P}\mathrm{P}\mathrm{V}=\frac{\mathrm{S}\mathrm{e}\times \mathrm{P}}{\mathrm{S}\mathrm{e}\times \mathrm{P}+\left(1-\mathrm{S}\mathrm{p}\right)\times \left(1-\mathrm{P}\right)}$$

(14)

Following the Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD) guideline [47], a CI of 95% of the performance indices was reported.

2.6. Ethical Considerations

The study protocol was reviewed and approved by the Isfahan University of Medical Sciences Research Ethical Committee (IUMSREC), with the following approvals: Modeling of incidence and outcomes of COVID-19: IR.MUI.RESEARCH.REC.1399.479 and Longitudinal epidemiologic investigation of patients’ characteristics with coronavirus infection referring to Isfahan Khorshid Hospital: IR.MUI.MED.REC.1399.029, conforming to the Declaration of Helsinki. Patient informed consent was obtained before admission to the current study. All data were kept confidential and had no personal identifiers. No minors participated in our study.

3. Results

Descriptive statistics were used to summarize the baseline characteristics of the study population. In our setting, 1600 COVID-19 patients were included in the study. Patients were categorized according to their LOS (≤7 days (n = 1165) as “normal”, >7 days (n = 435) as “prolonged”) in univariate comparison analysis. The median length of stay during the study period was 7.2 (IQR 4−9) days.

Table 1. Characteristics of hospitalized patients with COVID-19 in the Khorshid Cohort Study.

Parameters	Total (n = 1600)	Length of Stay (LOS)		p-Value b
		≤7 Days “Normal” (n = 1165)	>7 Days “Prolonged” (n = 435)
LOS, days a	6.01 (4.85)	3.76 (1.94)	12.11 (5.09)	<0.001
Age (>65 years)	562 (56.10%)	507 (43.50%)	55 (12.60%)	<0.001
Gender (% Female)	670 (48.80%)	464 (39.80%)	206 (47.30%)	0.001
Charlson Comorbidity Index (CCI) a	2.67 (2.13)	2.49 (2.11)	3.13 (2.13)	<0.001
Temperature maximum (≥38 degrees Celsius)	412 (25.75%)	322 (23.64%)	90 (20.68%)	0.745
Heart rate, beats per minute (<60 or >100)	478 (53.98%)	388 (33.3%)	90 (20.68%)	0.028
Respiratory rate, breaths per minute a	22.41 (5.67)	22.02 (5.27)	23.49 (6.56)	0.006
Systolic blood pressure (≥120 mmHg)	574 (35.80%)	247 (21.01%)	277 (63.70%)	<0.001
Diastolic blood pressure (≥90 mmHg)	218 (13.60%)	113 (9.60%)	105 (24.10%)	0.046
% O2 saturation minimum (<90)	754 (47.10%)	606 (52.01%)	148 (34.02%)	0.001
Neutrophils (<4 × 109/L)	956 (59.75%)	620 (53.22%)	336 (77.20%)	0.028
Lymphocytes (<1 × 109/L)	900 (96.40%)	621 (53.30%)	279 (64.10%)	0.028
Hemoglobin (<12 g/dL)	356 (22.30%)	293 (20.50%)	63 (14.40%)	0.085
Platelets (<150 × 109/L)	678 (59.75%)	480 (41.20%)	198 (45.51%)	0.142
Ferritin (>500 ng/mL)	94 (5.80%)	72 (6.01%)	22 (5.05%)	0.298
CRP (>30 mg/L)	685 (42.80%)	542 (46.52%)	143 (32.87%)	0.017
ESR (>60 mm/h)	420 (26.30%)	245 (21.03%)	175 (40.20%)	0.027
LDH (>222 U/L)	672 (42.00%)	416 (35.70%)	256 (58.80%)	0.046
D-dimer (>0.5 mg/L)	381 (23.80%)	95 (8.20%)	286 (65.70%)	0.036
AST (>35 IU/L)	1156 (72.30%)	749 (64.30%)	407 (93.50%)	0.330
HCO3 (mEq/L)	23.65 (3.67)	17.25 (3.76)	20.45 (2.78)	0.0123
ALT (>45 IU/L)	401 (25.10%)	305 (26.18%)	96 (22.06%)	0.204
Creatinine (>1 mg/dL)	822 (51.40%)	591 (45.40%)	231 (53.10%)	<0.001
Phosphorus (mg/dL) a	3.06 (0.85)	2.97 (0.85)	3.24 (0.81)	<0.001
Magnesium (mg/dL) a	1.96 (0.51)	1.95 (0.27)	1.99 (0.74)	0.335
Sodium (mEq/L) a	136.30 (4.13)	136.42 (3.94)	136.09 (4.46)	0.054
Potassium (mEq/L) a	4.02 (0.56)	3.99 (0.54)	4.08 (0.60)	0.055
BUN (mg/dL) a	19.79 (13.47)	18.67 (12.37)	21.92 (15.13)	<0.001
Total bilirubin (mg/dL) a	1.03 (2.17)	1.06 (2.61)	0.98 (0.61)	0.361

^a The percentage of the high-risk group (i.e., exposure) was provided in parentheses in total, “normal” or “prolonged” LOS subgroups when the high-risk cutoff was mentioned for parameters, and the standard deviation (SD) was provided (with “a” superscript) otherwise for the variables with an interval measurement scale. Such cutoffs were taken from the literature, and their citations were provided in the manuscript. For the predictor gender, the percentage of female subjects was provided in parentheses as the reference group. Statistical tests were selected based on the data’s nature and the variables’ distribution. ^b An independent-sample t-test was used for interval variables if the data were normally distributed; otherwise, the Mann–Whitney U test was employed. The Chi-square test was utilized for binary variables to compare proportions between the two independent groups. Note that the expected frequencies in any of the cells of the contingency table were more than five. ESR: erythrocyte sedimentation rate, LDH: lactate dehydrogenase, AST: aspartate transferase, HCO3: bicarbonate, ALT: alanine transaminase, BUN: blood urea nitrogen.

and

Table 2. Symptoms distribution between patients with normal and prolonged LOS.

Symptoms	Total	Length of Stay		p-Value a
		≤7 Days “Normal” (n = 1165)	>7 Days “Prolonged” (n = 435)
Fever	1118 (69.9%)	721 (61.9%)	397 (91.3%)	< 0.001
Cough	1125 (70.3%)	990 (85.0%)	135 (31.0%)	< 0.001
Myalgia	838 (52.4%)	562 (48.2%)	276 (63.4%)	< 0.001
Throat pain	255 (15.9%)	168 (14.4%)	87 (20.0%)	0.058
Weight Loss	259 (16.2%)	164 (14.1%)	95 (21.8%)	0.018
Chest pain	394 (24.6%)	279 (23.9%)	115 (26.4%)	0.365
Dizziness	97 (6.1%)	64 (5.5%)	33 (7.6%)	0.540
Headache	515 (32.2%)	372 (31.9%)	143 (32.9%)	0.112
Loss of smell and taste	186 (11.6%)	134 (11.5%)	52 (12.0%)	0.260
Diarrhea	377 (23.6%)	247 (21.2%)	130 (29.9%)	0.113
Vomiting	352 (22.0%)	233 (20.0%)	119 (27.4%)	0.478
Nausea	543 (33.9%)	373 (32.0%)	170 (39.1%)	0.518
Shortness of breath	995 (62.2%)	646 (55.5%)	349 (80.2%)	0.032
Stomachache	243 (15.2%)	166 (14.2%)	77 (17.7%)	0.393

^a The Chi-square test was utilized for binary variables to compare proportions between the two independent groups. Note that the expected frequencies in any of the cells of the contingency table were more than five.

summarize the descriptive statistics and the characteristics and symptoms of the patients considered in the study according to the length of stay categories.

For an example, an 86-year-old male COVID-19 patient with fever, cough, myalgia, sore throat, dizziness, diarrhea, stomachache and weight loss symptoms, but without chest pain, headache, loss of smell, vomiting, nausea, and short breath with CCI of 4, maximum body temperature of 36°, heart rate of 84 (beats per minute), respiratory rate of 16 (breaths per minute), systolic blood pressure of 105 (mmHg), diastolic blood pressure of 66 (mmHg), %O₂ saturation minimum of 90, neutrophils of 715 (×10⁹/L), lymphocytes of 264 (×10⁹/L), hemoglobin of 12.10 (g/dL), platelet of 126.00 (×10⁹/L), ferritin of 255.50 (ng/mL), CRP of 14.00 (mg/L), ESR of 26.00 (mm/h), LDH of 487.00 (U/L), D-dimer of 98.30 (mg/L), AST of 39.00 (IU/L), HCO3 of 31.00 (mEq/L), ALT of 16.00 (IU/L), creatinine of 0.73 (mg/dL), phosphorus of 2.56 (mg/dL), magnesium of 1.90 (mg/dL), sodium of 135.00 (mEq/L), potassium of 4.00 (mEq/L), BUN of 24.80 (mg/dL), and total bilirubin of 0.96 (mg/dL) had an LOS of 5 days in Khorshid Hospital.

Figure 1 shows the frequency distribution of LOS, which was right-skewed. The median age of patients was 59 (IQR 47–79) years (range 5–91), and 58% were male. Comorbidities were present in more than half of the patients, with hypertension being the most common comorbidity, followed by diabetes. The Charlson comorbidity index is presented in Table 1 for patients admitted for ≤7 and >7 days in hospital. The comorbidities score was significantly higher in patients with longer LOS (p-value < 0.05, Table 1). The cutoffs used for high-risk exposures were provided by the following: age [48], body temperature [49], heart rate [50,51], blood pressure [52], oxygen saturation [53], neutrophils, lymphocytes, hemoglobin, platelets, D-dimer [54], ferritin [55], CRP [56], ESR [57], LDH [58], AST [59], ALT [60], and creatinine [61].

Figure 2 shows the rate of different PCEP events among patients for both short and prolonged LOS. ICU admission is the most prevalent (55%) status among patients with prolonged LOS. There was a significant association between LOS and the PCEP binary variable (p-value < 0.001).

Fever (p-value < 0.001), cough (p-value < 0.001), myalgia (p-value < 0.001), weight loss (p-value 0.018), and shortness of breath (p-value 0.032) were significantly different in the LOS groups (Table 2). The cross-validated results of the proposed algorithm for the estimation of LOS as well as its “normal” and “prolonged” categories is provided in

Table 3. The cross-validated results of the proposed prediction algorithm in percent.

Indices	RMSE	MAE1	MAE2	R2	adj. R2	${\mathsf{\rho }}_{\mathbf{c}}$	Se	Sp	PPV	DOR	AUC	F1	MCC	K(C)
Value	1.58	1.22	0.98	89	81	94	92	91	79	112	91	80	79	79
95% CI-Lower	1.51	1.16	0.92	88	79	93	89	89	75	71	89	76	77	75
95% CI-Upper	1.64	1.28	1.05	91	84	95	95	93	83	179	94	85	81	83

MAE₁: mean absolute deviation; MAE₂: median absolute deviation.

. The most important features of the LOS prediction system were inflammatory markers, HCO₃, and fever.

The Bland–Altman plot of residual analysis is provided in Figure 3. Although the bias was not statistically significant (p-value = 0.777; paired-sample t-test), the residual error was higher in higher target LOS than lower LOS values. The residual error was further analyzed in “normal” and “prolonged” LOS groups. The estimation error values of 3.9% and 3.2% subjects surpassed the lower or higher 1.96 limits in the “normal” and “prolonged” LOS groups, respectively.

The ROC curve was then provided for the predicted LOS versus the binary ground truth (Figure 4). The best cutoff was calculated using the Youden index (J = Se + Sp − 1), estimated as >6.95, almost identical to our a-priori threshold.

We further predicted the importance of significant factors in the Kolmogorov–Gabor polynomials. The main predictors were only analyzed based on their normalized coefficients in the model. The seven most important factors are provided in Figure 5.

The unbiased PPV plot is provided in Figure 6 based on the prevalence of the prolonged LOS and Equation (14). The required parameters of Bayes’ theorem were assessed from Table 3.

4. Discussion

4.1. Implications

Medical researchers have recently been striving to enhance the quality and efficiency of healthcare systems and services. A significant aspect of this endeavor pertains to the LOS in the context of future outbreaks. Given the emergence of various variants of the virus responsible for COVID-19, accurately assessing or predicting LOS is becoming increasingly vital. An extended LOS not only impacts hospital capacity [62] but also escalates costs associated with outbreak management [63]. Hence, nations must plan for even the worst-case scenarios. This research delved into the risk factors of hospital admissions that influence the LOS among COVID-19 patients in Isfahan, Iran. We utilized a novel nonlinear artificial intelligence method for continuous data, focusing on comprehensive predictors. Our findings indicate that patients with prolonged hospital stays typically exhibited higher inflammatory markers, increased HCO₃, and more prevalent fever. These insights can guide clinicians in pinpointing specific risk factors linked to extended LOS. Moreover, our results serve as a benchmark for various models that could be applied in similar analyses, allowing healthcare professionals to narrow down critical variables for predicting LOS from the multitude recorded in hospital systems.

In our research, the median LOS was 7.2 days, with an interquartile range (IQR) of 4–9 days. It aligns closely with findings from a Chinese study [64], wherein the median hospital LOS fluctuated between 4 to 53 days over 45 domestic studies and 4 to 21 days across eight international studies. In contrast, a comprehensive report, drawing from data across 25 countries, recorded a median LOS of just 4 days and an IQR of 1–9 days [65]—substantially shorter than our recorded observations. Notably, our results surpassed the median LOS of 6 days documented in Saudi Arabia [66]. However, it fell short of the 16.4 days indicated in Indiana [67], aligning with the 8.5 days reported in the Mediterranean. These regional variances in LOS can be ascribed to many factors, including the infrastructure of healthcare facilities, the severity of treated cases, diverse admission and discharge protocols, and varying treatment approaches. Additionally, sociodemographic variables, especially age, are pivotal in influencing the observed international disparities in hospital stays.

4.2. Risk Factors

The consistency in hospital bed occupancy duration across various demographic groups in our research contrasts starkly with findings from a significant US study by Nguyen et al. [68]. Their research indicated that males typically had a more extended LOS than females. Due to the limited sample size in our cohort, we could not investigate the influence of gender on the likelihood of ICU admission. Furthermore, while our findings showed a consistent LOS regardless of demographic distinctions, European studies suggest a pronounced variability in hospital stays based on both gender and age [69]. In our data, the correlation between age and LOS was relatively weak (r = 0.134; p-value < 0.001).

In our research, the predominant symptoms upon admission, such as cough, fever, and shortness of breath, align with many earlier studies [70,71]. A systematic review and meta-analysis spanning 54 studies identified the most frequent symptoms in COVID-19 patients as follows: fever at 81.2% (95% CI: 77.9–84.4), cough at 58.5% (95% CI: 54.2–62.8), fatigue at 38.5% (95% CI: 30.6–45.3), dyspnea at 26.1% (95% CI: 20.4–31.8), and sputum production at 25.8% (95% CI: 21.1–30.4) [72]. Our findings concur with these percentages concerning fever and cough. However, the prevalence of shortness of breath in our study diverged.

Disturbances in total white blood cells, particularly lymphocytes, are often seen as the immune system’s response to inflammation. There is growing evidence that lymphopenia, characterized by a reduced lymphocyte count, significantly influences the trajectory of COVID-19, right from its onset to the eventual development of viral sepsis. This decrease in lymphocytes has been identified as a symptom of acute COVID-19, potentially resulting from direct damage inflicted by the virus [73]. Our findings regarding lymphopenia echo those of previous studies. Earlier research has outlined prognostic models that gauge the severity of SARS-CoV-2 infection by monitoring the lymphocyte-to-leukocyte ratio [74,75].

Recent studies have illustrated that lymphocyte counts below 5% were predominantly observed in patients exhibiting severe symptoms upon follow-up. There also appears to be a trend wherein lymphopenia is more pronounced and persistent among the gravely affected patients [76,77]. These studies also highlighted that patients with extended hospital stays typically had increased circulating lymphocytes, whereas their neutrophil counts were marginally diminished. This surge in lymphocytes might be attributed to rejuvenated production, given their rise both as a percentage of total blood cells and in absolute terms. Notably, the lymphocyte count was elevated across all patient severity subgroups [78], suggesting its potential role in extended LOS or heightened mortality risk [79].

Our study underscores the substantial influence of D-dimer on the length of hospital stay, aligning with the conclusions of other meta-analyses. These analyses indicate that D-dimer correlates with factors such as comorbidities, demographics, specific laboratory tests, radiological findings, the duration of hospitalization, complications, and ultimate outcomes. Such findings propose that D-dimer is a distinct biomarker, interfacing with other inflammatory cytokine markers indicative of organ or tissue damage. Furthermore, the interaction of acute-phase proteins with D-dimer implies that infection-driven inflammation (comprising cytokines and chemokines) instigates a state of hyper-fibrinolysis, a notion reinforced by D-dimer’s disconnect from the comprehensive coagulation panel [80,81].

Our study showed that patients with prolonged LOS among COVID-19 cases exhibited a significantly higher ESR. Many studies have assessed acute-phase responses to COVID-19 since the pandemic’s onset, and these frequently included ESR data [82,83]. A meta-analysis [84] further highlighted that elevated ESR levels were particularly pronounced in severe and fatal cases of COVID-19. Another comprehensive meta-analysis by Zhang et al. [85], which analyzed 28 studies encompassing 4663 cases, discovered that 61.2% of cases with increased ESR had a longer length of stay and were at a heightened risk for severe disease. Notably, variations in sedimentation rates between the groups were not explored.

There is a noticeable gap in the literature regarding using HCO₃ values as predictors for LOS. It might be attributed to the understanding that abnormal HCO₃ levels already indicate extended hospitalization [86]. It is plausible that these levels act more as process variables than predictors, a sentiment echoed by our findings. The serum HCO₃ level indicates the acid–base balance within the human body and is commonly assessed in routine biochemical tests, particularly as renal diseases advance [87]. Certain clinical studies have posited a potential role for serum HCO₃ levels in forecasting mortality from ailments beyond progressive renal disease. For instance, diminished HCO₃ levels have been linked to mortality from malignancies, while elevated HCO₃ levels have been associated with cardiovascular disease complications and related mortalities [88].

Low HCO₃ serum levels upon ICU admission significantly predict both short-term and long-term mortality. Additionally, a reduced serum HCO₃ serves as an indicator of acidosis. Past research confirms that acidosis can diminish systemic vascular resistance, exacerbating conditions like circulatory shock, impaired myocardial contraction, and tissue malperfusion. This cascade of complications can ultimately precipitate end-organ failure, including acute kidney injury, which might primarily contribute to the grim prognosis observed in critically ill patients.

4.3. The Properties of Kolmogorov–Gabor Polynomials

Among the diverse techniques employed for continuous prediction, we utilized Kolmogorov–Gabor polynomials. These are more popularly recognized as the Volterra series. They serve as instrumental methodologies in identifying and modeling nonlinear systems. They can adeptly capture a broad spectrum of nonlinear behaviors by executing a series expansion based on system input. Within the context of a hospital setting, when compared with other prediction algorithms, the Volterra series boasts several advantages [89]:

Flexibility: The Volterra series can depict many nonlinear systems, endorsing its versatility across diverse modeling landscapes.

Interpretability: A standout feature of the Volterra series is its capacity to demystify the system’s structure. It delineates the input–output relationship across linear, quadratic, and cubic terms. It facilitates a deeper comprehension of the system’s nonlinearity and subsequent impact.

Theoretical Foundations: The mathematical underpinnings of Kolmogorov–Gabor polynomials are well-established and rigorously studied, ensuring a robust theoretical base for their application.

However, it is essential to note that while these advantages are compelling, particular challenges and potential drawbacks also emerge [37]:

Computational Complexity: The expansion order’s escalation leads to an exponential growth in computational demands. It can hinder the feasibility of deploying high-order models, particularly when grappling with multitudinous inputs.

Overfitting: As with many adaptive models, there is an inherent overfitting risk when complexity overshadows the data’s intricacy. Such scenarios necessitate a meticulous model selection process to safeguard against over-optimization and ensure genuine applicability to new datasets.

In our research, we have employed regularization techniques. Additionally, by capping the interaction level at two, we have strategically mitigated computational demands and curtailed the risk of overfitting.

4.4. Performance Indices

Guarding against testing hypotheses suggested by the data (Type III errors) was guaranteed by cross-validation. The LOS prediction method showed strong agreement with the measured LOS (${\mathsf{\rho }}_{\mathrm{c}}$ = 0.94), and strong goodness-of-fit (R² = 0.8), and did not show a significant bias (p-value = 0.777; paired-sample t-test). However, the Bland–Altman error regression showed higher errors for lower LOS values.

The binary classification algorithm, on the other hand, showed a statistical power of 92%, a Type I error of 0.09, and a precision of 79%. It also had an excellent balanced diagnosis accuracy (AUC = 0.91), a high correlation between predicted and observed class labels (MCC = 0.79), and an excellent class labeling agreement rate (K(C) = 0.79). However, it is not entirely clinically reliable, as Type I errors must be less than 0.05, and the precision must be higher than or equal to 95% [34].

4.5. Comparison with the State-of-the-Art

We searched “Embase” for journal papers with the key words “(‘length of stay’/exp OR ‘length of stay’) AND (‘hospital’/exp OR ‘hospital’) AND (‘prediction’/exp OR ‘prediction’) AND (‘machine learning’/exp OR ‘machine learning’) AND (‘COVID’/exp OR ‘COVID’OR ‘coronavirus’/exp OR ‘coronavirus’)” without publication date condition. Among 64 screened, 45 records were excluded after analyzing their abstracts since they did not predict LOS. Nineteen records were thus assessed for eligibility. Journal papers with at least one prediction performance index and a sound ML methodology were included in

Table 4. The state-of-the-art to predict hospital LOS in COVID-19 patients.

Reference	Center/Region	Sample Size	Inputs	Important Features	Outputs	Models	Validation	Indices (the Best Method)	Important Characteristics
Ebinger et al., 2021 [17]	Cedars-Sinai Medical Center (Los Angeles), USA	966	353 variables	Age, respiratory rate, oxygen flow rate	LOS > 8 days vs. LOS ≤ 8 days	42 models	20% Hold-out	Se = 93% Sp = 63% F1 = 78% PPV = 67% AUC = 0.82	Missing imputation; cumulative day three information was used.
Hong et al., 2020 [16]	A tertiary care hospital in Zhejiang province, China	75	37 variables	Lymphocyte count, heart rate, cough, Epidermis, procalcitonin;	LOS > 14 days vs. LOS ≤ 14 days	Stepwise multivariable regression	No internal or external validation	AUC = 0.85 [CI 95: 0.75–0.94]	Missing imputation;
Orooji et al., 2022 [21]	Ayatollah Taleghani Hospital, Abadan, Iran	1225	53 variables	20 variables: Age, creatinine, WBC, lymphocyte/neutrophil count, BUN, ASP, ALT, LDH, activated PTT, coughing, hypertension, CVD, diabetes, dyspnea, oxygen therapy, pneumonia, GI complications, ESR, and CRP.	LOS	Statistical feature selection (correlation coefficient)+ MLP+ 12 training algorithms	10% Hold-out	RMSE = 1.6213 (days)	Patients who died within three days of admission were excluded (n = 128); selection bias. Missing data imputation.
Zhang et al., 2023 [24]	Zhengzhou University Hospital (Henan), China	384	83 variables	Immunotherapy, heparin, familial cluster, rhinorrhea (runny nose), and APTT	LOS	LASSO+ linear regression	Bootstrap validation (N = 2000)	R2 = 0.30	Missing data imputation (10 imputations);
Alabbad et al., 2022 [22]	King Fahad University hospital, Saudi Arabia	895	43 variables	Age, C-reactive protein (CRP), nasal oxygen support days	9-class ICU LOS	Random forest (RF) (the best classifier), gradient boosting (GB), extreme gradient boosting (XGBoost), and ensemble models	3-fold cross-validation	PPV = 94% Se = 94% F1 = 94%	Missing data imputation; SMOTE was used to balance nine classes to have 144 records each, biased performance indices. The original samples ranged from 12 to 144 for the classes; no admission date was provided.
Nemati et al., 2020 [15]	Global dataset	1182	Five variables	Age, sex	LOS	Stagewise GB (the best method), IPCRidge, CoxPH, Coxnet, Componentwise GB, Fast SVM, Fast Kernel SVM	No internal or external validation	C-index = 0.71	No comprehensive features except symptoms onset date, symptoms, and chronic disease binary variable
Usher et al., 2021 [18]	36 hospitals (Minnesota, Wisconsin, and the Dakotas)	2665	20 variables	Various variables, including age, critical illness, oxygen requirement, weight loss, and nursing home admission	LOS at >5, >10 and >15 days	GLM, RF (the best model)	5-fold cross-validation	AUC = 0.89	ICU admission, mechanical ventilation, and mortality risk are among the input features; selection and immortal-time bias.
Liuzzi et al., 2022 [20]	28 centers (Fondazione Don Carlo Gnocchi (FDG) Living COVID-19 Registry), Italy	222	829	55 variables: anagraphical data, admission clinical scales, admission signs and symptoms, admission supports, COVID-19 therapy, therapy prior to COVID-19, hematochemics	LOS	Sequential convolutional neural network	Repeated (N = 10) 5-fold cross-validation	MAE2 = 2.7 days (IQR = 3.0 days)	17 COVID-19 therapies were included in the input data; selection and immortal-time bias.
Mahboub et al., 2021 [19]	Rashid Hospital (Dubai), UAE	2017	22 variables	Urea, PLT, D-dimer, K+, anti-inflammatory medicine, antiviral medicine, mechanical ventilation, hemoglobin, azithromycin medicine, vitamin C medicine, painkiller medicine	LOS	Decision Tree	25% Hold-out	R2 = 0.5	In addition to mechanical ventilation, treatments were used as input features; selection and immortal-time bias.
Alam et al., 2023 [23]	Prince Sultan Hospital (Riyadh), Saudi Arabia	308	89 variables	Laboratory, X-ray, clinical data, and treatments, including LDH and D-dimer levels, lymphocyte count, and comorbidities such as hypertension and diabetes	Seven-class LOS	Tab Transformer	30% stratified hold-out	Pr = 83%, Se = 93%, F1 = 93% (discharged) Pr = 75%, Se = 98%, F1 = 84% (dead)	SMOTE-N oversampling technique was used to balance the classes and biased performance indices. Treatments, including anticoagulants, antibiotics, antivirals, an immunomodulators, were used as the inputs; selection and immortal-time bias.
This study	Khorshid Hospital (Isfahan), Iran	1600	42	Inflammatory markers (ESR, D-dimer, lymphocyte counts), HCO3, and fever	LOS and also LOS≤ 7 days vs. LOS > 7 days	The Kolmogorov–Gabor polynomial plus regularized least squares	Three-fold cross-validation	LOS: R2 = 0.89 [0.88–0.91], ${\mathsf{\rho }}_{\mathrm{c}}$ = 0.94 [0.93–0.95], RMSE = 1.58 [1.64–1.51] days MAE1 = 1.22 [1.28–1.16] days, MAE2 = 0.98 [0.92–1.05] days LOS categories: Se = 92% [89–95], Sp = 91% [89–93], PPV = 79% [75–83], AUC = 0.87 [84–89], F1 = 80% [76–85]	No class balancing was used. ICU admission, mechanical ventilation, and treatments were not used as the input features.

MAE₁: mean absolute deviation; MAE₂: median absolute deviation.

(10 methods as the state-of-the-art, besides the proposed algorithm “this study”).

Among the studies in Table 4, only Hong et al., 2020 [16] and our study followed the TRIPOD guideline [47] to report the 95% CI of the performance indices. In addition to transparency in reporting, it quantifies precision, uncertainty, reproducibility, and generalization. Only Alam et al. [23], Mahboub et al. [19], Liuzzi et al. [20], Usher et al. [18] and our study did not use missing imputations. The others used missing imputations. However, no analysis was performed to identify the reasons for missing data, i.e., missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR), which is critical in missing data analysis [90]. In our study, we not only did not include ICU admission, mechanical ventilation, or treatments as the inputs, but we also only used baseline information at hospital admission, which was not the case for Ebinger et al. [17], Usher et al. [18], Liuzzi et al. [20], Mahboub et al. [19], and Alam et al. [23].

Similarly to Alabbad et al. [22], Usher et al. [18], and Liuzzi et al. [20], we used cross-validation. Zhang et al. [24] used bootstrapped validation, though the 0.632+ bootstrap method is preferred in the literature [91]. Hold-out validation used by Ebinger et al. [17], Orooji et al. [21], Mahboub et al. [19], and Alam et al. [23] might introduce Type III error, and the repeated hold-out validation method is preferred. Also, Hong et al. [16] and Nemati et al. [15] did not use validation. Among the studies included in Table 4, only Nemati et al. [15], Usher et al. [18], and Liuzzi et al. [20] were multi-center. Orooji et al. [21] excluded subjects who died within 3 days of hospital admission, resulting in sampling bias. Our study is ranked in the top third based on the sample size. Moreover, Alabbad et al. [22] and Alam et al. [23] balanced the unbalanced training and test datasets, resulting in biased evaluation metrics, potential misleading improvement, and overfitting to the minority class. However, they had a better goodness-of-fit R² = 0.80 compared to other studies. Our study is the only one that reported the Bland–Altman plot critical to analyzing the residual error [35].

Like most studies in Table 4 [15,16,17,18,19,20], our study only focused on the first COVID-19 wave. However, Orooji et al. [21] considered the first, second, and third waves. Alam et al. [23] analyzed the first and second waves, and Zhang et al. [24] considered the Omicron variant. Thus, a direct comparison of the results of the proposed method and the other three methods [21,23,24] is not entirely rigorous.

4.6. Dichotomous LOS Definition

When the median LOS in our dataset was 7 days, then using a 7-day cutoff to dichotomize hospital LOS was statistically motivated: (1) Using the median as a cutoff point ensures that approximately half of the patients are categorized as “short stay” and the other half as “long stay”; (2) the median represents the robust central tendency of the data; patients with “prolonged” LOS are staying longer than the majority of patients, suggesting they might have different clinical characteristics, needs, or outcomes; (3) the median is robust to outliers and is not affected by very short or very long LOS values; and (4) the binary outcome can be directly tied to the dataset’s inherent structure, making the results more interpretable in the context of the data. However, it might make direct comparisons between the binary LOS model and other datasets or studies more challenging unless they also use a median LOS of 7 days.

4.7. Limitations and Future Research

Our study has several limitations. Firstly, given its single-center design with 1600 COVID-19 patients at a major academic hospital following specific institutional treatment protocols, the findings might not directly apply to other hospitals throughout Asian countries. More samples are required to improve the statistical power of the proposed method. Secondly, while we tried to control for disease severity in our analysis, we could not account for more subjective factors, including the nuances of treatment that might influence endpoint decisions. To comprehensively evaluate the potential impact of treatment on LOS, a prospective randomized trial is imperative. Thirdly, the Bland–Altman analysis of residual error highlighted a non-uniform error across measured LOS. Integrating the Bland–Altman parameters into the cost function will be a focal point of our future endeavors. While our initial findings demonstrate promising results, expanding the validation scope will provide a more holistic understanding of the model’s capabilities. Addressing these gaps in temporal and spatial validation will be instrumental in fostering confidence in our approach and ensuring its relevance across broader contexts, which is the focus of our future activity. Moreover, using multimodal image-processing prediction methods could, in principle, improve the reliability of the proposed algorithm [92,93], which is a focus of our future studies.

While the current model has been calibrated based on the original SARS-CoV-2 strain, the underlying framework holds potential for adaptation to newer strains. We can ensure its sustained relevance and accuracy in predicting LOS by continually updating and retraining the model with data from emerging variants. Integrating this dynamic model within the hospital health information system will facilitate real-time adaptability, making it a versatile tool for clinicians across different pandemic phases.

5. Conclusions

In this research, the utilization of machine learning models, notably the Volterra functional series, demonstrated a promising approach to predicting the length of stay (LOS) of COVID-19 patients. Validated on a significant dataset from Khorshid Hospital in Iran, the model showed strong performance metrics, including an R² of 0.8 and a concordance correlation coefficient of 0.94, indicating a good fit and a high agreement with the measured LOS. As noted in multiple studies, key features that played a vital role in LOS prediction were inflammatory markers, bicarbonate, and fever, aligning with the commonly observed symptoms in COVID-19 patients. The binary classification algorithm further provided insights into differentiating between “normal” and “prolonged” LOS groups. While the results present a substantial basis, there is room for improvement in the clinical reliability of the binary classification algorithm, especially concerning its Type I error and precision rate.

However, some limitations and considerations remain in the study. The Bland–Altman error regression indicated a higher error rate for patients with a lower LOS, suggesting potential areas for refinement in the model for this patient subgroup. Moreover, while our findings regarding the most prevalent symptoms upon admission were consistent with several other studies, there were notable discrepancies in the observed prevalence of shortness of breath. As healthcare providers and hospitals globally grapple with the challenges posed by the COVID-19 pandemic, findings from this research could pave the way for better resource management. Nonetheless, further temporal and spatial validation is imperative before generalized application. Future research endeavors could delve deeper into optimizing the model’s clinical reliability and expanding the model’s scope to other pertinent clinical outcomes. Further studies and medical regulations are essential to establish a dependable clinical prediction model suitable for smart hospitals.