2.1 Study population and data description

Of the 1200 patients initially screened, 45 were excluded from the study. Reasons for exclusion included missing critical data, failure to meet the inclusion criteria (such as lack of radiological confirmation of pneumonia or not requiring supplemental oxygen), or withdrawal from the study. Therefore, our study included 1155 adult patients hospitalized with radiologically confirmed pneumonia caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), all of whom required supplemental oxygen (≥24%). Among these patients, 843 were treated with RDV, while 312 received SOC, providing a diverse dataset for analysis. Fig. 1 shows the flow of participants in the study.

Fig. 1.Flowchart of participant selection and allocation. RDV: remdesivir; SOC: standard of care.

The allocation of RDV within our cohort was based on clinical criteria and the patient condition, following evolving medical guidelines. Criteria for RDV administration centered around the severity of the patient condition, indicated by factors such as the extent of hypoxemia and the immediate need for supplemental oxygen. Other considerations included the overall health status of the patient, using underlying conditions that could impact the effectiveness of the treatment or introduce potential risks. The timing of RDV therapy initiation was carefully aligned with national health authority recommendations and the latest clinical research findings. This strategic decision-making process reflects the agility of the healthcare system in integrating new therapeutic insights into patient care protocols, an important aspect of our analysis. Stratifying patients into RDV and SOC groups is integral to assess the impact of RDV alongside other non-invasive predictors on patient mortality. Specifically, our models concentrate on patients who received RDV treatment. The timing of RDV initiation—a key variable in our analysis—is relevant only within this patient group. Such timing allows us to isolate and evaluate the impact of RDV treatment on patient outcomes, especially when analyzed alongside key variables that are readily accessible to healthcare professionals and vital for clinical decision-making. Table 1 provides a comprehensive summary of the clinical characteristics of the patients included in the study, stratified by treatment group (RDV and SOC).

Our preliminary data analysis provides a transparent and detailed view of the patient cohort under study. In the RDV group, 97.1% of patients also received dexamethasone, with a mean age of 69.7 years and a standard deviation of 14.4 years, with 61.7% of them being male. In contrast, the SOC group had a mean age of 73.9 years and a standard deviation of 14.5 years, with 50.0% of them being male. Despite similarities in the prevalence of diabetes, hypertension, and chronic lung diseases between the groups, the RDV group had a higher proportion of overweight individuals, whereas the SOC group included more individuals with immunosuppressive conditions and smokers. Notably, patients in the RDV group had an average hospital stay of 4.25 days, shorter than the 5-day average of the SOC group, suggesting a potential benefit of RDV in reducing hospitalization duration. Furthermore, the relative risk of death during the hospital stay for the RDV group was 0.47 compared to the SOC group, indicating that patients treated with RDV had a 53% lower risk of mortality during hospitalization compared to those receiving SOC. This reduction in relative risk suggests a potential survival benefit associated with RDV treatment as part of a multifaceted approach to treating severely ill COVID-19 patients. To build on these findings, we are actively working on a more detailed analysis, which will be presented in an upcoming publication. This future publication will delve into the clinical implications of RDV treatment, specifically addressing how relative risk and other critical metrics impact patient outcomes. We aim to provide a comprehensive perspective that will further enhance understanding in the area under study.

2.2 Machine learning methods

We recall that our primary objective is to investigate the relationship between the predictor variables and the binary outcome of survival or death in hospitalized COVID-19 patients requiring ventilation, presented in Table 1. As previously mentioned in Subsection 2.1, our analysis is focused exclusively on the group treated with RDV, particularly assessing the impact of RDV treatment timing on patient outcomes. To effectively analyze the mentioned variables, we utilize GLMs due to their aptitude for describing binary outcomes like patient survival or death.

GLMs are advantageous as they describe the probability of an outcome (ranging from zero to one) as a function of various predictor variables. This modeling allows us to consider a range of relevant clinical and demographic factors within our patient cohort. In the GLM framework, the binary outcome, denoted as Y, indicates the survival (Y = 1) or non-survival (Y = 0) of the patients. The predictor variables, collectively referred to as X and presented in Table 1, encompass diverse predictors relevant to patient outcomes.

We use logistic regression, a type of GLM, to analyze the relationship between predictors and patient survival probability p = P(Y = 1| X = x). Logistic regression calculates the odds of survival as an odds ratio (OR), and mathematically is defined as p/(1 − p) providing a direct measure of the probability that a patient survives based on the specific characteristics and treatments. The OR derived from the logistic regression is a crucial metric in our analysis. The OR indicates the strength and direction of the association between each predictor and patient survival. An OR greater than one indicates an increased probability of survival, while an OR less than one suggests a decreased probability. In our logistic regression analysis, we model the logarithm of the OR (log-OR), denoted as Z and its observed value as z, establishing a linear combination of predictors defined by Z = β₀ + β₁x₁ + ^... + β_kx_k, where β₀, β₁, …, β_k are the model parameters, each reflecting the impact of the predictors X₁, …, X_k based on their observed values x₁, …, x_k, on patient survival. To convert log-ORs into the probability of survival, p say, we use the sigmoid function stated as:

Which converts a real number into a value at [0, 1]. This conversion translates the linear combination of predictors into a probability value. Such a conversion provides a more intuitive understanding of the probability of a patient survival based on the characteristics and treatments. To train the model, we optimize the parameters β_j, for j = 0, 1, …, n, to minimize a cost function formulated as:

Where y_i is the observed class (survival or non-survival) for instance i in the dataset, $\hat{p}$_i is the corresponding estimated/predicted probability of survival for that instance, and n is the number of instances in the dataset.

The cost function stated in Eqn. 2 is convex, which ensures that an optimization algorithm, such as gradient descent, can find the global minimum, leading to the best set of parameters for the model.

Once trained, the model uses the estimated regression parameters to calculate the survival probability using the function stated in Eqn. 1. Instance i belongs to the survival class if $\hat{p}$_i > 0.5, that is, $\hat{y}$_i = 1 when $\hat{p}$_i > 0.5; or $\hat{y}$_i = 0, otherwise, where $\hat{y}$_i is the estimated value of y_i. Despite being called “regression”, logistic regression is mainly used for classification. It is effective for linearly separable classes but its performance decreases with non-linear classes. High correlation among predictors can lead to overfitting, destabilizing the coefficients. Overfitting occurs when a model learns of the training data closely, including its noise, leading to poor performance on unobserved data.

For the case of non-linear data or when a more flexible decision boundary is needed, SVM is preferred over logistic regression. SVMs belong to supervised machine learning designed for classification, but also cater to regression tasks [43, 44]. SVM creates a hyperplane that best separates classes in the data. In a two-dimensional space, this hyperplane is a line. For data that are linearly separable, an SVM finds the line that maximizes the margin between the nearest points (support vectors) of each class. This line is termed the decision boundary, and the area delineated by the margin is known as the “street”. For non-linearly separable data, SVM uses a kernel trick, which transforms the data into a high-dimensional space to be partitioned by a hyperplane. Once the hyperplane is determined in this space, it may be mapped back to the original space, resulting in a non-linear decision boundary.

Data often contain noise and classes may overlap, which is challenging to identify a hyperplane that separates the classes. To address such a challenge, SVM introduces a soft margin allowing some points to fall on the wrong side of the margin or even the decision boundary. This wrong classification is avoided balancing maximization of the margin and minimizing misclassification. While SVMs provide classification and regression, we must acknowledge its limitations when dealing with large datasets due to computational demands. However, the concept of optimizing the margin of fit, in classification or regression, remains consistent. SVM regression, addressing non-linearities using kernel trick and its margin-based approach, accommodates as instances as possible within a margin, whose width is regulated by a parameter, balancing complexity of the model and its fit to the data.

In scenarios dealing with large datasets, complex interactions, or the necessity for stability against overfitting, ensemble learning methods, like RF, emerge as powerful solutions. RF is particularly beneficial when the limitations of single-model approaches, such as logistic regression, become evident in handling complex data patterns.

RF creates multiple decision trees during training and aggregate their outputs, following a “wisdom of the crowd” approach. This is represented by the mode of the classes from individual trees [45, 46]. The RF algorithm involves the following steps:

• Step 1—Begin by selecting random samples from the dataset with replacement, known bootstrap samples [47], where reach sample typically includes B instances, with B being the same as the total number of instances n in the original dataset or a different specified number.

• Step 2—Create a decision tree for each bootstrap sample and at each node of the tree considering:

- Step 2.1—Select a subset of variables (features), denoted by m, at random from the total number of variables, k say, in the dataset, where typically m is chosen to be the square root of k.

- Step 2.2—Determine the optimal split on these m variables to divide the node, with the split maximizing the reduction in variance (for regression trees) or minimizing the Gini impurity (for classification trees), and then the node is divided into left and right child nodes. The process of finding the optimal split is mathematically represented as $\mathrm{\Delta }I(S,A)=I(S)-\sum _{v\in values(A)}(|S_v|/|S|)I(S_v)$, where I(S) is the impurity measure of the original set, values(A) represents the set of all possible values for attribute x, S_v is the subset of S for which attribute x has value v, and ∣S_v∣/∣S∣ is the proportion of instances that have value v for attribute x. The attribute that provides the highest information gain is chosen for the split.

• Step 3—Pass a new object through all the trees in the forest, where each tree casts a vote for a class, and the forest selects the class with the majority of votes from all the trees as its prediction.

A crucial parameter of the RF is the number of trees to generate. A larger number of trees reduces the possibility of overfitting but increases the computational complexity. RFs are robust to outliers, scalable, and capable of modeling non-linear decision boundaries due to their hierarchical structure. Moreover, RF automatically selects the variables that are important in the classification.

RFs can model complex non-linear decision boundaries and do not require input features to be scaled, that is, to have the same range, unlike methods such as SVMs or logistic regression. Despite the power and robustness of an RF, its intricate nature may complicate the decision-making processes. While an RF may struggle with interpretability, it shares with some other models the need for computational resources, particularly when considering large ensembles of deep trees.

Comparatively, stochastic gradient boosting, also known as GBM, combines boosting with regularization to enhance model predictions [48]. Unlike RFs, which construct trees in parallel and focus on reducing variance among their models, GBM builds trees sequentially, concentrating on minimizing bias and improving accuracy of the model by correcting errors from previous trees.

The strength of GBM lies in its ability to combine multiple weak decision-tree-based learners, sequentially refining the model to form a more potent learner. This sequential learning allows GBM to reduce both bias and variance, making it particularly helpful in scenarios where predictive accuracy is paramount. Although each decision tree individually might not yield the best predictions in GBM, they provide improved performance collectively. The trees are constructed to minimize a specific loss function, with every new tree focusing on the previous errors of the model. A distinguishing feature of the GBM is its inclusion of a regularization component that is used to circumvent overfitting. By introducing a regularization term, GBM seeks simplicity of the model, often preferring those ones that balance complexity and performance. The GBM framework emphasizes the model by minimizing the chosen loss function, which evaluates the discrepancy between predicted and observed outcomes.

Given a training set {x_i, y_i}, for i ∈ {1, …, N}, where x_i represents the predictor variable for instance i and y_i the corresponding label, GBM employs an ensemble of decision trees to compute the prediction w_i for each instance i defined as:

Where f_j represents tree j in the ensemble, contributing scores to its leaves, and T indicates the total number of trees.

The efficacy of each tree j is evaluated using an aggregated loss function L_j, which integrates both the prediction error and a regularization component, stated as:

$$L_j=\sum _{i=1}^{N}l(y_i,w_i+f_j\left(x_i\right))+\mathrm{\Omega }\left(f_j\right),j\in \{1,\dots ,T\}$$

Where l represents the individual loss function, assessing the prediction error for each instance i, w_i is defined in Eqn. 3 and denotes the aggregate prediction up to tree j for the sample x_i, while Ω introduces a regularization term to penalize the complexity of the model. Common examples of l include the squared error for regression tasks and the logistic loss for classification tasks. They aim to tune the model predictions to closely align with the observed labels. The function f_j(x_i) minimizes L, leading to constructing each tree, iteratively refining the performance of the model.

The choice between GBM and RF hinges on specific requirements. The trade-off involves considerations of interpretability, computational efficiency, and predictive accuracy. RF may be preferred for its simplicity and robustness, especially when interpretability and computational resources are restricted. Conversely, GBM might be favored in applications where achieving the highest predictive accuracy is paramount and the data are complex.

2.3 Fuzzy methods

Within the continuum of modeling techniques that address data complexities, FRBCS emerge as a versatile alternative to ensemble models like GBM and RF, as well as a broad array of traditional and advanced methodologies [49, 50, 51]. By leveraging the concept of fuzziness, FRBCS excel in managing uncertainties and ambiguities that traditional algorithms often fail to address, making it particularly valuable in scenarios where conventional (crisp) logic proves inadequate [39]. The widespread application of fuzzy logic across disciplines such as artificial intelligence, computer science, and mathematics has driven advancements, showing its flexibility and broad applicability.

FRBCS provide good solutions to both classification and regression problems by integrating human-like reasoning with computational precision. At the core of FRBCS is the fuzzification interface, which converts precise inputs into degrees of membership using various functions. This conversion accommodates the inherent vagueness present in data, making FRBCS powerful tools for handling complex scenarios where traditional crisp logic fall short.

In the fuzzification process, we begin with a universe of discourse, χ say, representing the complete set of all possible values that a variable can assume. An element v within this universe represents a specific value of the variable. A fuzzy set A is defined as a subset of this universe, where each element v is assigned a membership degree by the function μ_A(v). This degree, ranging from 0 to 1, indicates the extent to which the element v belongs to the fuzzy set A, thereby transforming precise input values into a range of possibilities that can be processed by the fuzzy logic system.

Various types of membership functions can be used in the context of the present study. The triangular membership function is defined as:

$${\mu }_A(v)=\{\begin{array}{ll}0,& v\le a\text{ or }v\ge b;\\ \frac{v-a}{c-a},& a\le v<c;\\ \frac{b-v}{b-c},& c<v\le b;\end{array}$$

Where a is the lower limit, b the upper limit, and c the modal value with a ≤ c ≤ b.

The trapezoidal membership function is established as:

$${\mu }_A(v)=\{\begin{array}{ll}0,& v\le a\text{ or }v\ge d;\\ \frac{v-a}{c-a},& a\le v<c;\\ 1,& c\le v\le c^{\mathrm{\prime }};\\ \frac{d-v}{d-c^{\mathrm{\prime }}},& c^{\mathrm{\prime }}<v\le d;\end{array}$$

Where a is the lower limit, d the upper limit, and [c, c′] the interval modal value with a ≤ c ≤ c′ ≤ d.

The Gaussian membership function is expressed as:

$${\mu }_A\left(v\right)=\text{exp}(-\frac{{(v-c)}^2}{2{\sigma }^2})$$

Where c is the mean and σ the standard deviation of the Gaussian model.

Studies have shown that different membership functions can yield varying degrees of accuracy depending on the nature of the input data and the problem domain. The fuzzy sets obtained through fuzzification are then integrated within a fuzzy logic system, which typically includes:

• A fuzzy knowledge base consisting of a database (DB) and a rule-base (RB). The DB defines fuzzy sets with memberships expressed between 0 and 1, while the RB stores fuzzy if-then inference rules (IF-THEN rules). The premises rely on the fuzzy sets defined in the DB, while the conclusions specify the resulting actions or outcomes.

• An inference engine that processes fuzzy IF-THEN rules using methods tailored to fuzzy rule-based classification systems.

• A defuzzification process that converts fuzzy sets into actionable outputs, potentially integrated within the inference process.

In line with the basis of a fuzzy logic system, our FRBCS implementation employs several specific strategies to effectively address the complexities inherent in medical diagnosis as the following:

• Max-membership principle—This principle directly selects the class with the highest membership degree as the output, avoiding the need to convert fuzzy sets into crisp values. The principle ensures direct and actionable results, particularly in scenarios involving ambiguous or overlapping data.

• Rule weighting (FRBCS.W) strategy—In this strategy, each fuzzy rule is assigned a weight that reflects its relevance in the classification process. By prioritizing the most relevant rules, the model improves prediction accuracy, especially when dealing with complex and diverse data patterns.

• Space partitioning, the chi-square (FRBCS.CHI) method—This method applies the chi-square test to partition the input space, ensuring that the resulting fuzzy partitions are statistically significant. It enhances the alignment between fuzzy rules and the underlying data distribution, leading to more reliable and interpretable classifications.

• Genetic (GFS.GCCL) algorithms—These algorithms optimize the fuzzy rule sets, utilizing cooperative-competitive learning to evolve and refine the rules. Genetic algorithms are particularly effective in exploring large and intricate search spaces, making them ideal for optimizing fuzzy systems in complex applications.

• Pittsburgh (FH.GBML) method—This method focuses on evolving entire populations of fuzzy IF-THEN rule sets simultaneously, optimizing the entire rule set collectively. This holistic method is effective in managing complex classification tasks, ensuring cohesive and accurate outcomes.

In our FRBCS models, defuzzification is seamlessly integrated into the inference process, with the max-membership principle serving as the primary method. This guarantees that outputs are actionable by directly selecting the class with the highest membership degree, without converting fuzzy sets into crisp values. These FRBCS strategies are implemented using a comprehensive software package, as referenced in [52], which supports the deployment of these advanced fuzzy logic methods for complex real-world applications.

The performance of the FRBCS models in medical diagnosis is closely tied to the quality of data and complexity of decision boundaries. The interaction of symptoms demands flexible rule sets that adapt to evolving data patterns. By incorporating machine learning and big data, our methodology aims to enhance the adaptability of fuzzy models, ensuring precision and relevance in dynamic environments – essential aspects for successful real-world applications.

2.4 Performance metrics and validation methods

Evaluating the predictive performance and reliability of a trained model is essential. This evaluation encompasses a variety of metrics and techniques, such as cross-validation or a simple train-test split. We provide an overview of these metrics and techniques, with further details presented in [53].

Before discussing specific metrics, it is important to define the fundamental concepts of true positive (TP), true negative (TN), false positive (FP), and false negative (FN). These concepts state the outcomes of the predictions of a binary classification model in relation to the observed values and correspond to:

• TP—Correctly predicted positive observations.

• TN—Correctly predicted negative observations.

• FP—Incorrectly predicted positive observations.

• FN—Incorrectly predicted negative observations.

The confusion matrix, a tool to visualize the performance of a binary classification model, summarizes the concepts of TP, TN, FP, and FN, as shown in Table 2.

Other metrics from Table 2 include precision, sensitivity (recall), accuracy, and the area under the curve (AUC) of the receiver operating characteristic (ROC). Precision and sensitivity focus on positive predictions and are defined as:

$$\text{precision}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}}$$

$$\text{sensitivity (or recall)}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

Accuracy provides a broader measure of performance for a model and is established as:

$$\text{accuracy (binary)}=\frac{\mathrm{T}\mathrm{N}+\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

The AUC is a metric used to evaluate the performance of binary classification models. The ROC curve plots the TP rate (sensitivity) against the FP rate (1 − specificity) at various threshold settings. The AUC is a metric that represents the probability that a randomly chosen positive instance is ranked higher by the model than a randomly chosen negative instance. An AUC of 0.5 indicates no discriminative ability (equivalent to random guessing), while an AUC of 1.0 signifies perfect discrimination. This metric is particularly helpful in medical research as it provides an aggregate measure of performance across all classification thresholds, assessing a model ability to distinguish between positive and negative cases irrespective of the decision threshold.

In addition to AUC, the area under the precision-recall curve (AUPRC) is especially informative for evaluating the performance of models on imbalanced datasets. Unlike the AUC-ROC, which considers both FPs and FNs equally, the AUPRC focuses on the performance concerning the positive class (the minority class). In imbalanced datasets, where the negative class dominates, the AUPRC offers a better measure of model ability to correctly identify the positive cases by emphasizing precision (the proportion of TPs among all predicted positives) and recall. A high AUPRC indicates that the model performs well in distinguishing the minority class from the majority class, which is crucial in medical research where the consequences of FNs can be particularly severe, such as in disease detection or risk prediction.

We also calculate bootstrapped confidence intervals (CIs) for accuracy, sensitivity, specificity, and AUC. Bootstrapping is a non-parametric statistical method that enables us to estimate the distribution of these metrics by repeatedly sampling, with replacement, from the original data. Bootstrapping allows for the calculation of CIs without relying on parametric assumptions. By generating multiple bootstrap samples and computing the desired metric for each, we obtain an estimate of the metric variability, making bootstrapping especially helpful when the data distribution is unknown or the sample size is small.

Additionally, the choice between cross-validation and a train-test split depends on the dataset and the type of modeling. Cross-validation, particularly k-fold, is preferred when seeking reliable performance, especially with small datasets, as it maximizes the use of data for both training and testing, providing a measure of model performance. However, it is computationally more intensive than a train-test split, which might be more suitable for very large datasets or when a quicker model assessment is required due to its computational efficiency. This ensures that the chosen method aligns with the specific requirements and constraints of our study, balancing accuracy, and computational resources.

2.5 Methodology and software

Next, we outline the methodology and software tools employed in our study, emphasizing computational considerations and the rationale behind our methodological choices. We utilized the R software and caret library [54] for machine learning and data modeling, focusing on mortality outcome classification. Although caret provides access to a variety of models, we tailored our methodology to address the unique challenges posed by COVID-19 data, selecting models that best fit our research questions and data characteristics.

For model training and evaluation, we employed the train function within caret, which allows the application of different modeling techniques through its versatile method parameter. Our methodology included models such as GLM with stepwise selection based on the Akaike information criterion (AIC)—glmStepAIC, glmnet, GBM (gbm), RF (rf), and SVM with a radial basis function kernel (svmRadial).

Elastic net regularization, used by the glmnet model, combines the strengths of both least absolute shrinkage and selection operator (LASSO) and ridge regressions. LASSO regression, employing the L1 norm (sum of absolute values of coefficients), is effective in variable selection by shrinking some coefficients to zero, so simplifying the model. Ridge regression, using the L2 norm (sum of squares of coefficients), helps to manage multicollinearity by distributing the coefficients among correlated predictors. The elastic net produces balancing, offering a model that generalizes well, particularly when dealing with high-dimensional datasets.

We optimized and evaluated the models based on the following criteria:

• Stopping criteria—Model training was halted using a predefined threshold for improvement in predictive accuracy on a validation set, aiming to minimize overfitting.

• Computational cost—Assessed by tracking the training time and resource usage for each model, prioritizing efficiency while maintaining performance.

• Solver for optimization—For models requiring optimization, such as glmnet and GBM, we used gradient descent and its variants, selected for their effectiveness in minimizing the cost function of the model architecture.

Given that the fuzzy methods utilized in our study do not directly provide a means to determine the relevance of individual variables due to their fuzzification rather than assigning direct weights, we employed permutation feature importance. This involves modifying the values of each predictor to assess its impact on model accuracy, allowing us to identify the variables most relevant in predicting COVID-19 mortality outcomes.

The pseudo-code for feature importance analysis has the following steps:

• Step 1—For each feature in the dataset:

- Step 1.1—Permute the values of the feature.

- Step 1.2—Evaluate model performance with the permuted data.

- Step 1.3—Calculate the change in performance.

• Step 2—Rank features based on performance change.

• Step 3—Identify key features influencing model accuracy.

These steps help to discern the relative importance of variables, thereby facilitating focused analysis and informed decision-making. For transparency and enabling replication of our findings, the detailed code, including the permutation feature importance analysis, is available upon request. The feature importance analysis is outlined in Algorithm 1 and illustrated in a flowchart in Fig. 2.

Fig. 2.Schematic representation of the permutation feature importance algorithm.

Algorithm 1: Permutation feature importance analysis.

Input: Given dataset and pre-trained fuzzy model.

Output: Ranked features based on importance.

1. Start with a pre-trained fuzzy model that has been fitted to the dataset.

2. For each feature in the dataset do:

2.1 Create a copy of the original dataset.

2.2 Permute the values of the current feature in the copied dataset.

2.3 Evaluate the model using the permuted dataset and a specified performance metric.

2.4 Record the change in performance compared to the base performance.

3. Repeat Steps 1–2 for each feature in the dataset.

4. Rank the features according to the degree of change they caused after all features have been permuted and their effect on performance recorded.

5. State as more important those features causing a larger drop in performance.

6. Determine which features are relevant based on the previous steps.

3.1 Variable selection and data preprocessing

The predictor variables were chosen for their clinical relevance, ease of access for healthcare professionals, and non-invasive nature, which are defined as:

• Age—The patient’s age in years, a basic metric for assessing risk factors associated with age-related vulnerabilities.

• Time_Adm_Sym—The duration from symptom onset to hospital admission, measured in days, reflecting the timeliness of treatment initiation.

• RDV_In_Sym—The interval from symptom onset to the start of RDV treatment, providing insights into the effectiveness of early antiviral interventions.

• Hospital_duration—The total duration of hospitalization in days. This variable is only available at discharge or death, but it is included to assess its potential impact on mortality predictions in retrospective analyses. By examining the length of stay, we can gain insights into the progression of severe cases and identify factors that may contribute to longer stays, which are often associated with higher mortality. These insights could guide future research or influence decision-making in similar patient cohorts.

• FiO2Max—The maximum fraction of inspired oxygen required during hospitalization, expressed as a percentage, serving as a non-invasive proxy for respiratory compromise severity.

These variables contained no missing data, ensuring a consistent dataset, thereby eliminating the need for imputation. To ensure comparability across variables with different scales, all continuous predictor variables were normalized using min–max scaling, transforming their values to a 0–1 range. This transformation was crucial for maintaining consistency across all models, including the fuzzy systems.

The dataset was then divided into training and testing sets using an 80/20 split. This split maintained a representative distribution of survival and death outcomes in both sets. The same training and testing sets were applied across all models to maintain consistency and comparability.

3.2 Non-fuzzy model adjustments

To optimize the performance of the non-fuzzy models, including GBM, glmnet, RF and SVM, we conducted hyperparameter tuning, which involves selecting the best values for parameters that govern the model training process. This process enhances the model ability to learn from the data. Unlike these models, logistic regression does not require hyperparameter tuning, as it directly estimates relationships between predictors and outcomes through maximum likelihood estimation, relying on the data rather than tunable settings. In SVM, sigma controls the width of the Gaussian kernel, determining the influence range of a single training example, while c represents the penalty for misclassification. For GBM, n.trees is the number of trees, interaction.depth sets the tree depth, shrinkage controls the learning rate, and n.minobsinnode specifies the minimum observations per terminal node. For RF, mtry is the number of variables randomly sampled as candidates at each split in a tree. In glmnet, alpha determines the balance between LASSO and ridge penalties in glmnet, while lambda controls the strength of the regularization applied. The most frequently used hyperparameter settings for each model are summarized in Table 3.

To ensure the robustness of our results, a 5-fold cross-validation was employed on the training set, where the data were split into five subsets: four for training and one for testing. This process was repeated five times, ensuring each subset was tested once. Such a cross-validation provides stable estimation, making a separate validation set redundant. The validation was consistently applied across all models, ensuring comparability and consistency in the results.

3.3 Logistic regression

Logistic regression was implemented using the same cross-validation framework described above, although it did not require hyperparameter tuning. This regression was directly fitted to the training dataset, and the statistically significant predictors identified are presented in Table 4, reporting the impact of each predictor on the probability of death, assuming all other variables remain constant.

The OR quantifies the change in odds of death associated with a one-unit increase in each predictor variable, holding other factors constant. The interpretation of the coefficients reported in Table 4 is as follows:

• For every additional year of age, the odds of death increase by 6.5% (OR = 1.065), holding all other variables constant.

• For every additional day increase in Time_Adm_Sym, the odds of death decrease by 18.5% (OR = 0.815), assuming all other variables remain constant.

• An increase in RDV_In_Sym by one day raises the odds of death by 11.8% (OR = 1.118), provided all other predictors are held constant.

• For every additional day of Hospital_duration, the odds of death decrease by 4.9% (OR = 0.951), holding all else constant.

• For every unit increase in FiO2Max, the odds of death increase by 7% (OR = 1.070), assuming all other variables remain constant.

After fitting the model to the training data, we applied the logistic regression model to the test set to evaluate its predictive performance. The resulting confusion matrix is shown in Table 5.

Table 6 presents the performance metrics with 95% CIs, showing strong predictive performance. The high sensitivity indicates the model ability to accurately identify survival cases, while the AUC reflects its overall discriminative ability. However, the specificity, while lower, suggests a moderate ability to identify death cases accurately.

Together, Tables 5 and 6 illustrate the effectiveness of the logistic regression model in predicting patient outcomes, with a high accuracy and robust sensitivity, making it a valuable tool for clinical decision-making.

3.4 GLMNet model and coefficient estimation

The glmnet model is guided by two crucial parameters: alpha and lambda. The alpha parameter determines the balance between LASSO and ridge regressions within the model, with values ranging from 0 (pure ridge regression) to 1 (pure LASSO regression). In our study, the alpha and lambda parameters were determined through cross-validation to optimize the balance between model complexity and predictive accuracy. Specifically, alpha was set to 0.3, which introduces a bias toward ridge regression while retaining some feature selection capabilities of LASSO regression. The value of lambda was set to 0.0136, controlling the penalty on the model coefficients to prevent overfitting (see Table 3). Table 7 shows the estimated coefficients and their corresponding 95% CIs, obtained using bootstrapping to provide good estimates.

The interpretation of the coefficients and their 95% CIs is as follows:

• The coefficient for Age (estimated at 0.0493) has a 95% CI of [0.0316, 0.0667], suggesting an association between increased age and higher mortality risk.

• The coefficient for Time_Adm_Sym (estimated at −0.0905) has a 95% CI of [−0.1734, −0.0213], indicating a protective effect of earlier hospital admission.

• The coefficient for RDV_In_Sym (estimated at 0.0320) has a 95% CI of [−0.0059, 0.0884], which includes zero, indicating that the effect of RDV initiation timing on mortality is statistically nonsignificant at 5%.

• The coefficient for Hospital_duration (estimated at −0.0279) has a 95% CI of [−0.0521, −0.0065], indicating that longer hospital stays are associated with reduced mortality risk.

• The coefficient for FiO2Max (estimated at 0.0531) has a 95% CI of [0.0468, 0.0599], indicating a positive association between higher oxygen requirements and increased mortality risk.

The performance of the glmnet model is presented in Table 8. Both the glmnet and logistic regression models demonstrate strong performance across various metrics, although each has its strengths in different areas. The performance of the other models is shown in Table 9.

To further illustrate the trade-off between sensitivity and specificity across all models, we provide the ROC curves in Fig. 3. To compare the performance of the models in the context of the imbalanced dataset, we also computed the AUPRC. The glmnet model achieved an AUPRC of 0.803, while the logistic regression model showed a higher AUPRC of 0.841. For the other models, we obtained AUPRC values of 0.678, 0.570 and 0.534 for SVM, GBM and RF, respectively. These results further emphasize the superior reliability of logistic regression model in predicting the positive class, followed by the glmnet model.

Fig. 3.ROC curves for the GBM, GLMNet, logistic regression, RF and SVM models using Portuguese COVID-19 data. GBM: gradient boosting machine; SVM: support vector machines.

Based on all the given information, each model showed substantial predictive power. However, the GLMs particularly distinguished themselves due to their interpretability. In GLMs, the coefficients directly reflect the impact of each predictor on the outcome, assuming all other variables remain constant, offering a balance of interpretability and flexibility. For example, a coefficient for Age indicates how each additional year increases or decreases the odds of mortality, which clinicians can directly use to assess risk. Therefore, even though the SVM with a radial basis function kernel displayed a high AUC, GLMs—specifically the logistic regression model—strike a balance between interpretability and flexibility, making them an appealing choice for statistical modeling where understanding the relationships between predictors and outcomes is vital. This interpretability is particularly crucial in clinical settings, where decisions must be transparent and easily understood by healthcare professionals.

3.5 Comparison of fuzzy modeling methods

Now, we implement and compared different fuzzy modeling methods, including FRBCS.W, FRBCS.CHI, FH.GBML and GFS.GCCL, using the same training dataset of Portuguese COVID-19 cases that was employed for non-fuzzy models. We evaluate the effectiveness of these fuzzy methods in predicting mortality among COVID-19 patients and we compare their performance against non-fuzzy models in terms of accuracy, sensitivity, specificity, and the interpretability of the results.

Each fuzzy modeling method was implemented using the frbs.learn function, with specific control parameters that are crucial for defining the behavior of the models. The control parameters include num.labels (number of fuzzy labels), type.mf (membership function type), type.tnorm (t-norm), type.snorm (s-norm), and type.implication.func (implication function). An example in algorithmic form detailing the procedure for one of these methods is presented in Algorithm 2 and illustrated in the flowchart of Fig. 4.

Fig. 4.Flowchart of the fuzzy model generation and evaluation process.

Algorithm 2: Fuzzy model generation and evaluation.

Input: Chosen method, training data, input data range, control parameters;

Output: Evaluation of the fuzzy model performance and identification of relevant variables.

1. Initialize the fuzzy model using the chosen method, training data, input data range, and control parameters.

2. Use the code: control = list (num.labels = 7, type.mf = “GAUSSIAN”, type.tnorm = “MIN”, type.snorm = “MAX”, type.implication.func = “ZADEH”).

3. Train the fuzzy model on the training data.

4. Predict outcomes for the test data utilizing the trained fuzzy model.

5. Display a confusion matrix to assess the model performance.

6. Summarize the results of the FRBCS model, including relevant variables.

7. Use the summary to determine which variables most impact the model accuracy.

The control parameters used in the fuzzy models, including the example shown in Algorithm 2, play a crucial role in defining the model behavior. Specifically, num.labels = 7 divides the data into seven fuzzy categories, providing a finer granularity compared to simpler models. For instance, the fuzzy sets might represent levels such as “very low”, “low”, “slightly below average”, “average”, “slightly above average”, “high”, and “very high”. Choice of the membership function (type.mf = “GAUSSIAN”) ensures smooth transitions between these levels, while the t-norm and s-norm (type.tnorm = “MIN”, type.snorm = “MAX”) define the logic operations for intersection and union within the fuzzy system. Lastly, type.implication.func = “ZADEH” specifies how implications are handled within the model, following Zadeh fuzzy logic principles. These control parameters are essential in tailoring the fuzzy models to reflect the complexity of the input data, especially when dealing with subtle differences in the dataset.

It is important to note that the control parameters may vary across different fuzzy methods to best suit each method characteristics and the specific application within the context of COVID-19 data. When predicting outcomes for new and test data using the generated fuzzy logic, the model applies a set of fuzzy rules to each input case. For example, a rule might be: IF Age is “high” AND FiO2Max is “high”, THEN the outcome is 80% “death” and 20% “alive”. Here, the fuzzy output reflects the degrees of membership to each category, representing the uncertainty in the classification. The set of fuzzy rules generated by the model is then used to classify the test data, where each input is evaluated according to these rules, resulting in a fuzzy output that indicates the degrees of membership to categories such as “death” and “alive”.

To illustrate how the model classifies a new case, consider the following example. A new patient has the characteristics: Age classified as “very high” FiO2Max as “high”, Time_Adm_Sym as “moderate”, RDV_In_Sym as “early”, and Hospital_duration as “short”. The fuzzy model might apply a combination of rules, such as: IF Age is “very high” AND FiO2Max is “high” AND Time_Adm_Sym is “moderate” AND RDV_In_Sym is “early” AND Hospital_duration is “short”, THEN the outcome is 85% “death” and 15% “alive”.

In practice, the model may not have an exact rule that perfectly matches the new case. Instead, it evaluates the case based on the available rules that most closely align with the patient characteristics. The model aggregates the results of these rules, producing a fuzzy output that represents the degrees of membership to each category, such as “death” and “alive”. For this example, the model might assign an 85% membership to the category “death” and 15% to the category “alive”, reflecting the combined uncertainty and assessment of the patient risk based on the five variables.

If we need to decide, such as determining a crisp classification, the fuzzy output can be defuzzified. For instance, using the centroid method, we calculate a score based on the weighted average of the membership degrees, which then lead to a decision threshold (if the result is above 50%, classify as “death”). This allows the fuzzy model to handle uncertainty and provide predictions, which are then translated into actionable classifications.

After applying the fuzzy rules and obtaining a fuzzy output, the output can be defuzzified to a crisp value, such as “death” or “alive”. In our example, defuzzification is performed using the centroid method, where the weighted average of the membership degrees is utilized to determine the definitive category. Defuzzification is necessary in our study to enable a direct comparison with non-fuzzy models in terms of accuracy, sensitivity, and specificity.

The accuracy, sensitivity, and specificity for each fuzzy method were obtained as follows:

• FRBCS.W—Accuracy = 0.818, sensitivity = 0.829, and specificity = 0.650.

• FRBCS.CHI—Accuracy = 0.857, sensitivity = 0.828, and specificity = 1.00.

• FH.GBML—Accuracy = 0.837, sensitivity = 0.829, and specificity = 0.692.

• GFS.GCCL—Accuracy = 0.745, sensitivity = 0.829, and specificity = 0.525.

From the comparison of the fuzzy methods used, it is clear that FRBCS.CHI had the highest accuracy (0.857), specificity (1.0), and AUC (0.929). Both FRBCS.W and FH.GBML had similar accuracy (0.818 and 0.837, respectively) and sensitivity (0.829 for both), while GFS.GCCL had the lowest accuracy (0.745) and specificity (0.525). Based on these results, FRBCS.CHI shows the best overall performance, with high accuracy and specificity.

Regarding the results of the relevance of variables in the FRBCS.CHI model, FiO2Max emerges as the most relevant variable, accounting for approximately 44.42% of the predictive power of the model. This means that FiO2Max contributes the most to the model ability to make accurate predictions. Following FiO2Max, Age contributes 19.47%, RDV_In_Sym a 13.32%, Time_Adm_Sym a 13.20%, and Hospital_duration a 9.58%. These percentages reflect the relative contribution of each variable to the overall predictive accuracy of the model, as determined by the permutation feature importance analysis, which measures how much the model prediction performance decreases when each variable values are randomly permuted. The feature importance of predictor variables based on permutation method with Portuguese COVID-19 data is plotted in Fig. 5.

Fig. 5.Bar plot of the feature importance of predictor variables based on permutation method with Portuguese COVID-19 data. FiO2Max: maximum fraction of inspired oxygen required during hospitalization; Hospital_duration: length of stay in days; OR: odds ratio; RDV_In_Sym: time from symptom onset to RDV start in days; Time_Adm_Sym: time from symptom onset to hospital admission in days.

When comparing the fuzzy results and the logistic regression model, the consistency of the selected predictors in both models stands out. In the logistic regression, all the five variables are statistically significant at 5%. While fuzzy models provide insights into which variables are important, they do not inherently indicate the direction of the relationship (as whether an increase in a variable increases or decreases the probability of an outcome). This contrasts with GLMs, where coefficients directly indicate how changes in predictors influence outcomes, both in direction and magnitude. For example, a positive estimated coefficient for Age in a GLM indicates that as age increases, the estimated probability of death also increases.

3.6 Discussion on methodological choices

We next explain our choice of analytical methods and their alignment with the goal of precisely predicting COVID-19 mortality.

In our dataset, the relationships between predictor variables and clinical outcomes may not be strictly linear or binary. While some relationships could be approximately linear or binary under certain conditions, others might involve more complex, non-linear, or gradual transitions. Non-linear relationships can be effectively modeled using various machine learning techniques, such as GBM, GLMs, RF and SVM, with non-linear terms. However, conventional models may struggle with data that lack clear boundaries. For example, a slight change in FiO2Max or Age may not immediately alter the estimated survival probability but could gradually shift the associated risk. Fuzzy logic models are well-suited for such cases, as they allow for the modeling of such uncertainties by assigning degrees of membership to different categories, such as “high chance of survival” or “moderate risk of death”. The fuzzy models are particularly valuable in clinical decision-making, where it is more informative to express that a patient has, for example, a 70% chance of recovery and a 30% risk of deterioration, rather than categorizing the patient as simply “survival” or “non-survival”.

The fuzzy models used in this study—FRBCS.W, FRBCS.CHI, FH.GBML and GFS.GCCL—are specifically designed to handle these inherent uncertainties in COVID-19 clinical data. Such models generate fuzzy rules based on the input data, where each rule encapsulates a relationship between predictor variables and the outcome. Additionally, each of these fuzzy models applies different strategies for rule generation, optimization, and handling of uncertainties, which distinguishes them from one another. For instance, FRBCS.W emphasizes rule weighting, FRBCS.CHI focuses on generating rules using the chi-square method, FH.GBML applies genetic-based optimization techniques to enhance rule quality, and GFS.GCCL uses genetic algorithms with cooperative coevolutionary learning to optimize the fuzzy rules. These distinct models allow each of them to handle uncertainties and capture complex relationships in COVID-19 data differently, providing a comprehensive analysis of the factors influencing patient outcomes.