Treatment effect of remdesivir on the mortality of hospitalised COVID-19 patients in Switzerland across different patient groups: a tree-based model analysis

explored whether the treatment response to remdesivir differed by patient characteristics. METHODS: We analysed data collected from a hospital surveillance study conducted in 21 referral hospitals in Switzerland between 2020 and 2022. We applied model-based recursive partitioning to group patients by the association between treatment levels and mortality. We included either treatment (levels: none, remdesivir within 7 days of symptom onset, remdesivir after 7 days, or another treatment), age and sex, or treatment only as regression variables. Candidate partitioning variables included a range of risk factors and comorbidities (and age and sex unless included in regression). We repeated the analyses using local centring to correct the results for the propensity to receive treatment. RESULTS: Overall (n = 21,790 patients), remdesivir within 7 days was associated with increased mortality (adjusted hazard ratios 1.28–1.54 versus no treatment). The CURB-65 score caused the most instability in the regression parameters of the model. When adjusted for age and sex, patients receiving remdesivir within 7 days of onset had higher mortality than those not treated in all identified eight patient groups. When age and sex were included as partitioning variables instead, the number of groups increased to 19–20; in five to six of those branches, mortality was lower among patients who received early remdesivir. Factors determining the groups where remde-sivir was potentially beneficial included the presence of oncological comorbidities, male sex, and high age. CONCLUSIONS: Some subgroups of patients, such as individuals with oncological comorbidities or elderly males, may benefit from remdesivir.


Introduction
Remdesivir is an antiviral medication originally intended to treat hepatitis C and respiratory syncytial virus infections [1]. Although unsuccessful for its original purpose, remdesivir proved effective against the Ebola virus disease and Marburg virus. Positive results on remdesivir's in vitro effectiveness against coronaviruses made it a promising treatment candidate for coronavirus disease 2019   [2,3]. Between May and July 2020, several regulatory agencies issued an emergency or conditional authorisation to treat patients hospitalised for COVID-19 with remdesivir [4]. Remdesivir also became the first fully licensed drug against COVID-19 by the United States Food and Drug Administration based on preliminary results from a randomised controlled trial (RCT) [5,6]. However, the initial results were controversial: another RCT found no overall benefit [7], and there were concerns about the manufacturer's role in the decision-making process [8].
In the context of the early RCT results, the World Health Organization (WHO) expert groups conducted a set of mortality trials ("Solidarity") on four repurposed drugs. Remdesivir, hydroxychloroquine, lopinavir, and interferon beta-1a were included to determine whether any of them could reduce in-hospital mortality in patients hospitalised with COVID-19 [9]. Remdesivir was the only drug not discontinued after the interim report was published. The results showed a small benefit only among patients who were not ventilated during drug administration. The negative findings from the Solidarity trial were a major contributor to the WHO recommendation in November 2020 against using remdesivir to treat COVID-19 [10]. However, in April 2022, a weak conditional recommendation was issued to administer remdesivir for patients with non-severe COVID-19 but at the highest risk of hospitalisation [11].
These RCTs had small sample sizes and thus could not assess the potential modifying effect of characteristics such as age, sex, and comorbidities on the treatment outcome. Although the trials found no differences between remdesivir and the comparators overall, some results hinted that the effect of remdesivir may be heterogeneous across the population. Therefore, we analysed our large COVID-19 hospital database in Switzerland [12] to assess mortality among patients who received remdesivir and identify subpopulations that may benefit from remdesivir treatment.

Data sources
The COVID-19 Hospital Surveillance System (CH-SUR) database is coordinated by the Federal Office of Public Health (FOPH) in collaboration with the University of Geneva [12]. It contains data on hospitalised patients diagnosed with COVID-19. To date, 21 university and cantonal hospitals across Switzerland participate in the system, covering most tertiary hospitals in Switzerland. Data collection was approved by the Ethics Committee of the Canton of Geneva, Switzerland (CCER, 2020-00827) and all local ethics committees. Study planning, conduct, and reporting were aligned with the Declaration of Helsinki (2013 revision). An up-to-date cohort description, protocol, codebook, questionnaire set, and further detail is available on the website https://www.unige.ch/medecine/hospital-covid/.
Two stays of the same patient were considered two episodes if the time between the discharge and new admission dates is more than 30 days or if the second hospitalisation was in a different hospital without data being shared. Various information was collected for each episode in a standardised manner, including demographics, admission details, complementary clinical details (comorbidities, complications, admission to an intermediate care or intensive care unit, and treatments), and follow-up status (death, discharge, or transfer).
We included adults (aged ≥18 years) with a laboratoryconfirmed COVID-19 diagnosis who were hospitalised for more than 24 hours and recorded in CH-SUR between 1 st February 2020 and 7 th February 2022. The follow-up period was from the day of hospitalisation until the final discharge date of the first episode, death, or database closure.
For subjects with more than one episode, only the first episode was included. We included patients who received either no treatment, remdesivir as the only (antiviral) treatment, or any treatment but not remdesivir. Patients who received remdesivir were split into those starting treatment within one week of the onset of symptoms (or hospitalisation, if hospitalised prior to symptom onset or data on symptoms were missing); and those starting later [13]. Patients who received remdesivir with other treatments were excluded since it would be challenging to differentiate between the effect of remdesivir and other drugs. Records were excluded if missing a patient's outcome, outcome date, age, or sex or if records had inconsistencies in key dates. Treatment other than antivirals, interferon, and chloroquine (e.g. corticosteroids, antibodies, and immunomodulatory treatment) were not included in the analysis.

Tree model
We split the patients into groups with similar treatment responses using model-based recursive partitioning (from this point forward, tree model) and the model4you R package (see "Technical description of the model-based recursive partitioning" and supplementary figure S1 in the appendix) [14][15][16]. In this method, a multivariable Cox proportional hazards regression analysis is performed on the entire dataset first. It adjusts for selected covariables, including the main variables of interest (treatment). Next, the patients are iteratively partitioned into groups by assessing the instabilities in the model parameters according to a pre-defined list of patient characteristics. A parameter is considered unstable if the partial derivative of its contribution to the partial likelihood correlates with at least one patient characteristic. The stability is assessed by splitting the dataset into two parts based on the values of each variable on the list using the Bonferroni-adjusted permutation test of statistics of a quadratic form [17]. If any parameter with significant instability (p <0.01) is found, we split the dataset into two parts according to the parameter with the lowest p-value. The same process is repeated for both groups until we find no parameter with instability.
In real-world settings, the decision to treat a patient with remdesivir or another drug usually depends on the patient's characteristics, risk factors, comorbidities, and clinical signs and symptoms. To reduce the risk of confounding by indication, we considered a method called local centring [18][19][20]. The treatment variable T 0 (which usually would be either 1 or 0, depending on whether the patient received the treatment or not) is replaced by T = T 0 -P(X), where P(x) is the propensity of receiving treatment based on the covariables X. In other words, the higher the probability of receiving remdesivir, the higher the benefit needs to be to contribute to the regression with the same effect. We estimated the propensities of all three treatment options by random forests using the partykit package in R (see "Technical description of the model-based recursive partitioning" in the appendix) [21].

Analyses
We performed four analyses using the tree model described above. First, we conducted two analyses where treatment, age, and sex were included in the initial Cox model as in-dependent variables. The levels of treatment were: none; remdesivir within a week of symptom onset; remdesivir later; treatment with drugs other than remdesivir. Analysis I was without, and analysis II was with local centring.
Second, we conducted two analyses where treatment was the only independent variable in the Cox model, also without and with local centring (analyses III and IV, respectively). Remdesivir was included as a time-depending covariable by splitting the person-time into two records: from enrolment to treatment initiation and from treatment initiation to the end of follow-up. Other treatment was included as constant, ignoring the starting time of treatment.
We selected a set of metabolic markers and comorbidities for the partitioning that may be associated with the prognosis of COVID-19. The following variables, all measured at baseline, were included (table 1): body mass index (BMI) [22]; high blood urea nitrogen level (>19 mg/dl) [23]; high respiratory rate (>30/min) [24]; low blood pressure (diastolic <60 mmHg or systolic <90 mmHg) [23]; having chronic obstructive pulmonary disease (COPD) [23]; having a chronic cardiovascular disease [23]; having a chronic renal disease [25]; and having an oncological condition (cancer, tumour) [23]. We also included a composite variable based on a dichotomised CURB-65 score (≥2 vs <2). The CURB-65 score was defined as meeting at least two of the following conditions: confusion (abbreviated mental test score <9); high urea nitrogen level; high respiratory rate; low blood pressure; and age ≥65 years, (all cut-offs defined above) [26]. Age and sex were included as partitioning variables in analyses III and IV. Age and BMI were considered continuous, and all others as dichotomous variables. The determination of propensity scores was per-formed using all the variables mentioned above. Missing data were imputed using the R package mice and all included explanatory variables.
We present the results of the Cox models as adjusted hazard ratios (aHR) with corresponding 95% confidence intervals (CI). The tree model results are presented graphically, showing the process from the full dataset into the final partitioning. For each leaf node (i.e., a subgroup defined according to the partitioning factors), the primary outcome is the regression coefficient of remdesivir use within seven days of symptom onset on mortality. The secondary outcomes were the regression coefficients for other covariables (late remdesivir use, use of another treatment, and age and sex if not included as partitioning variables). All analyses were conducted in R (version 4.2.1).

Sensitivity analyses
We conducted four sensitivity analyses. The first three were restricted to different time periods to account for the differences in variants, clinical practice, and other factors. In sensitivity analysis 1, we included the first wave (1 st February 2020 until 31 st October 2020). Sensitivity analysis 2 included the first two waves (until 28 th February 2021). Sensitivity analysis 3 included the first to fourth waves (until 20 th October 2021). Finally, we conducted sensitivity analysis 4, where the CURB-65 score (a composite of other variables included in the tree model) was excluded from the partitioning.

Role of the funding source
The funders had no role in the study design, the collection, analysis, or interpretation of data, the writing of the man- Dichotomous >19 mg/dl blood urea nitrogen: yes or no T T T T High respiratory rate [24] Dichotomous >30/minute: yes or no T T T T Low blood pressure [23] Dichotomous Diastolic <60 mmHg or systolic <90 mmHg: yes or no T T T T High CURB65 score [26] Dichotomous At least two (vs 1 or none) of the following indicators: age ≥65 years, high urea nitrogen, high respiratory rate, low blood pressure, abbreviated mental test score <9.

T T T T
Respiratory disease [23] Dichotomous The patient has a chronic obstructive pulmonary disease: yes or no T T T T Cardiovascular disease [23] Dichotomous The patient has a chronic cardiovascular disease: yes or no T T T T Chronic renal disease [25] Dichotomous The patient has a chronic renal disease: yes or no T T T T Oncological condition [23] Dichotomous The patient has been diagnosed with cancer, a tumour, or other oncological pathological finding: yes or no T T T T C, included in the initial Cox model as a covariable; T, included in the tree model as a partitioning variable.
* Local centring refers to an adjustment of the variable indicator, where it is defined as T = T 0 -P(x) where T 0 is the crude indicator (1 or 0), and P(x) is the probability that T 0 = 1 given the covariables x.
uscript, or in the decision to submit the paper for publication.

Patient characteristics
The initial dataset included 30,653 records of 29,639 patients. Eighty-four records representing subsequent episodes with at least a 30-day gap were excluded. The records of patients returning to the hospital within 30 days were merged. As mentioned, records were excluded if there was a missing outcome or sex or had missing or inconsistent dates related to hospital admission or treatment initiation (n = 4176, 14.1%). They were also excluded for discharge on the day of entry (n = 143, 0.5%), receiving both remdesivir with other tratment (n = 198, 0.9%), and age <18 years (n = 921, 3.2%; supplementary figure S2). All 2411 patients (8.1%) from one hospital were excluded because of the disproportionately high risk of a missing outcome.
Ultimately, 21,790 patients were included (table 2): 1347 (6.2%) received remdesivir within the first seven days after symptom onset; 946 (4.3%) started remdesivir more than seven days after symptom onset; and 2517 (11.6%) received another treatment (chloroquine, ribavirin, interferon, and/or lopinavir combined with ritonavir). Patients who received remdesivir were younger than those who did not receive treatment, particularly those who started treat-  ceiving treatment other than remdesivir was not associated with mortality.
The partitioning did not differ between the models without or with local centring (figure 2; supplementary tables S1-S2). CURB-65 score caused the highest instability. Patients with high CURB-65 scores were next split by res-piratory rate and those with low CURB-65 scores by oncological conditions. Ultimately, eight leaf nodes were identified. Almost half of the observations were in one leaf node (low CURB-65 score, no oncological conditions or renal disease, normal blood urea nitrogen). Early remdesivir use was not associated with lower mortality in any of

Original article
Swiss Med Wkly. 2023;153:40095 the nodes. Late administration of remdesivir was beneficial in two cases: in patients with oncological conditions and a high CURB-65 score but a normal respiratory rate; and in patients with a low CURB-65 score and a high urea nitrogen level but no oncological or renal conditions.

Analyses III and IV adjusted for treatment in the initial model
The initial Cox model was adjusted only for treatment. Early treatment with remdesivir remained significantly associated with increased mortality in the overall population both without (analysis III; figure 3A; aHR = 1.29, 95% CI 1.13-1.48) and with (analysis IV; figure 3B; aHR = 1.28, 95% CI 1.11-1.47) local centring, although the effect size was smaller than in analyses I and II. Neither late remdesivir treatment nor treatment with other antivirals was associated with mortality.
CURB-65 score again had the highest instability ( figure 4). Age had the next highest instability, split at 78 for those with high CURB-65 scores and 73 years for those with low scores (figure 5; supplementary tables S3-S4). The trees with and without local centring were almost identical. The number of leaf nodes was 20 in analysis III and 19 in analysis IV. One leaf node was clearly the largest, covering about a fourth of the observations. It included patients with a low CURB-65 score, age ≤59 years, no oncological conditions, and a normal blood urea nitrogen level. Early remdesivir use was beneficial in five nodes in both analyses. However, the sample sizes were relatively small, ranging from 192 to 663, and regression coefficients ranged between -0.26 and -0.91 without local centring and between -0.25 and -0.80 with local centring. The conditions determining the largest group were the presence of oncological comorbidity, a low CURB-65 score, and an age below 59 years. In addition, without local centring, patients with a high CURB-65 score and age ≤64 years had a lower risk of death when treated early with remdesivir.

Sensitivity analyses
In sensitivity analysis 1 (restricted to the first wave of the pandemic), early remdesivir use was associated with reduced mortality in the models without age or sex, although the association was significant only without local centring (aHR 0.66, 95% CI 0.44-0.97; supplementary table S5). In the model with age and sex, mortality was reduced moderately but not significantly. Early remdesivir use was particularly beneficial among those with a high CURB-65 score; patients with a low CURB-65 score and oncological conditions had higher mortality when remdesivir was administered early. When age and sex were used as splitting variables, a benefit of early remdesivir was seen, particularly among those over 78 years old.
When the second wave was included (sensitivity analysis 2), the protective effect of remdesivir disappeared: mortality was higher among those receiving remdesivir early (supplementary table S6). The same was seen in sensitivity analysis 3 when the period was extended to the end of the fourth wave. The trees were similar to those in the main analysis although not identical: the main components of the partitioning, and the branches where remdesivir was beneficial, were the same as in the main analysis (supplementary table S7).

Original article
Swiss Med Wkly. 2023;153:40095 In sensitivity analysis 4 (without CURB-65 score), blood urea nitrogen became the key partitioning variable. In the analyses with age and sex in the Cox model, all treatment levels were beneficial among patients who had high urea nitrogen and low blood pressure but a normal respiratory rate; however, the effect was minor among those treated early with remdesivir compared to other treatment options.
In the analyses without age and sex in the Cox model, there were benefits of remdesivir use in two groups: older patients with a high urea nitrogen level and high respiratory rate; and older women with renal disease but normal urea nitrogen level and normal respiratory rate (supplementary  table S8).

Discussion
Hospitalised COVID-19 patients in Switzerland treated with remdesivir had higher mortality than those who received no antiviral treatment, but the association was not uniform between sub-populations. Overall, the CURB-65 score -a composite variable that combines older age (≥65 years) and the presence of selected clinical and laboratory findings -caused the most instability in the regression parameters. However, it was challenging to find patterns in the factors potentially modifying the treatment effect of remdesivir. The factors frequently appearing on the pathways to effective remdesivir included oncological comorbidities, male sex, and old age.
The CURB-65 score has been validated to predict mortality in community-acquired pneumonia [26], but its performance for COVID-19 outcomes is limited [27]. Nevertheless, CURB-65 was the most important variable in the partitioning process, possibly because it is a composite variable that combines several predictors. When CURB-65 was omitted from the analysis, urea nitrogen level and respiratory rate became the key components causing instability. Younger (≤59 years) individuals with oncological co-

Original article
Swiss Med Wkly. 2023;153:40095 morbidities were the largest group among those with low CURB-65 scores (≤1) to see benefits with remdesivir; their regression coefficients corresponded with at least 22% lower odds of death when they received early remdesivir than without treatment. Cancer patients usually have a weakened immune system because of cancer therapy and are thus more susceptible to SARS-CoV-2 infection. Studies from the early stages of the pandemic also suggested higher mortality among cancer patients than non-cancer patients [28,29]. Nevertheless, the evidence related to antiviral treatment for COVID-19 in cancer patients is limited [30].
Despite the promising results in the early stages of the pandemic, evidence of the benefits of remdesivir is scarce. In April 2022, the WHO Living Guideline on Therapeutics and COVID-19 was updated with a weak conditional rec-ommendation supporting remdesivir for patients with nonsevere COVID-19 but at a high risk of being hospitalised [11]. An update is also expected for patients with severe COVID-19. The recommendation is based on up-to-date evidence from five trials with 2,710 patients. Although the differences in mortality and recovery speed are minimal, remdesivir clearly reduces the risk of hospitalisation; it prevented 73 hospital admissions among 1,000 patients, which is of the same magnitude as other treatments including nirmatrelvir, molnupiravir, sotrovimab, and casirivimab-imdevimab [11].
Another living systematic review has identified five RCTs assessing remdesivir treatment, two overlapping with the WHO review [31]. The conclusions were similar: no difference in mortality was observed, but there may be benefits in non-mortality outcomes. Nevertheless, to our knowl- with local centring. CURB-65, the composite score for severe pneumonia depending on confusion, urea nitrogen, respiratory rate, blood pressure, and age; Resp. rate, respiratory rate per minute; BUN, blood urea nitrogen level (mg/dL); age, age in years.

Original article
Swiss Med Wkly. 2023;153:40095 edge, no RCT has shown any harmful outcomes associated with remdesivir. Remdesivir has also not been associated with any major adverse events [32]. Instead, a recent nonrandomised study found a significant reduction in mortality in patients receiving remdesivir and dexamethasone versus dexamethasone alone [33]. Despite the lack of randomisation, the study arms were similar, at least regarding known patient characteristics.
Previous studies have suggested that remdesivir treatment is beneficial when given at an early stage [11,13]. Our findings did not support this claim. The previous studies also had limitations that mitigate the strength of this finding. To our knowledge, the only RCT that assessed the role of treatment timing was conducted in China in the early days of the epidemic with less than 300 patients. The remdesivir arm had lower mortality than the placebo arm when treatment/placebo was started within 10 days of symptom onset, but if started later, remdesivir resulted in higher mortality than placebo [7]. However, these associations were not significant. A retrospective cohort study from the United Arab Emirates found lower mortality among patients receiving remdesivir within the first seven days than those treated later [34]. In another prospective study from Italy, receiving remdesivir within five days versus later also reduced mortality [35]. However, the reason for delayed treatment was the delayed admission to the hospital in most cases, so it is questionable to what extent these differences can be attributable to remdesivir. A large cohort study from Hong Kong compared the timing of remdesivir initiation relative to dexamethasone and found lower mortality if remdesivir was initiated before or at the same time as dexamethasone [36]. The overall mortality in our cohort was also comparable to that in other studies where remdesivir was much more common [37].

Limitations
Our study has several limitations. First, the allocation of remdesivir was not random. Many of the variables we assessed probably influenced the decision to administer remdesivir. We attempted to mitigate the risk of the indication confounding the results through local centring; however, the results were similar with and without this correction method. Second, the variables for the partitioning were selected subjectively based on previous knowledge. Third, we excluded patients who received remdesivir together with other drugs. However, such patients were <10% of all patients treated with remdesivir, and excluding this group is unlikely to have a major impact. Fourth, we did not adjust for the centre effect; the 20 hospitals likely have different practices for administering different drugs. Fifth, we refrained from quantifying the uncertainty around the findings; although the mortality was lower in patients treated with remdesivir in some patient groups, this could be due to chance, especially in the smallest subgroups.
The results should be interpreted as a description of the CH-SUR data and not generalised to the overall population. Further, the tree model does not address associations between the partitioning variables and the outcome; it only attempts to divide the population into subgroups depending on the treatment response. The findings should thus not be misinterpreted as predictors of the effect of remdesivir. Finally, the selection of variables for the partitioning and regression model was made unsystematically.

Conclusion
In conclusion, in Swiss hospitals, COVID-19 patients treated with remdesivir had a higher risk of death than those without treatment. However, this rather controversial finding is likely a result of unmeasured confounders or artefacts related to the decisions on drug administration. While there is a growing consensus that the ability of remdesivir to reduce mortality is minimal [11,31], some individuals may still benefit from it. Our results demonstrate that the association between remdesivir use and mortality varies substantially between patient groups. In particular, it would be worth examining the role of remdesivir in treating elderly males and patients with oncological co-

Original article
Swiss Med Wkly. 2023;153:40095 morbidities. Vaccines have become widely available to control the pandemic, and the evidence on COVID-19 is constantly growing and maturing. Therefore, the treatment and management of COVID-19 are moving from an emergency response towards a precision medicine approach [38] with the help of large collaborative studies.

Data sharing statement
The anonymised data can be accessed through a multistage process described elsewhere (https://www.unige.ch/ medecine/hospital-covid/files/4015/9427/6937/ CH_SUR_Multicentric_process_final.pdf). Applicants must complete a concept sheet and send it to the study team. An Executive Committee of experts and hospital participant representatives will review the concept. Depending on the goal of the analysis, additional ethics clearance might be needed. Data will be restricted to the request and shared through a secure platform.

Supplementary Text S1. Technical description of the model-based recursive partitioning.
Model-based recursive partitioning methods, also known as tree-based methods for subgroup analyses, use splitting procedures for partitioning patients into groups with higher and lower treatment effect. The treatment effect is first estimated by an initial regression model under the assumption that the parameters (i.e. the coefficient of the covariables in the regression equation) are applicable to all patients. If presence of differing treatment effects in different subgroups is found, this parameter will be estimated for each of these subgroups separately.
The partitioning of the patients into subgroups is based on finding instabilities in the model parameters.
We consider a parameter unstable if the partial derivative of the parameter's contribution to the objective function does not fluctuate randomly around zero but correlates to at least one patient characteristic. Hence the tree-based model relies on the assumption that the differential treatment effect can be explained as a function of the patient characteristics.
Below, we describe the process of model-based recursive partitioning we used in this analysis in detail. For the practical implementation, we used the R function pmtree (package model4you version 0.9-7; https://cran.r-project.org/web/packages/model4you/index.html) with the respective Cox model (coxph, package survival version 3.2-13; https://cran.r-project.org/web/packages/survival/index.html) as the initial model and the cleaned dataset with variables listed in the main text (Table 1) as data. The full set of R scripts (without the data) is available on request from the corresponding author (janne.estill@unige.ch).

Initial model
The first step of the algorithm is to compute the general model ( , | ) for all patients of the study sample, where Y is the outcome, X=(X n ) are the included covariables, and θ = (β 0 , β 1 ,..., β n ) the set of parameters where β 0 is the intercept and β i , i=1,…,n, the coefficients of the variables X i . As an initial model we chose the semi-parametric Cox model defined by where h(t) is the hazard of the event (in our case, death) as function of time, h 0 (t) is the baseline hazard (when all covariables have the default value; h 0 can take arbitrary form); x i are the covariables and β i the corresponding parameters.
To solve the Cox model, we will find the set of parameters that minimize the objective function Ψ, in our case, the negative log-likelihood. This is equivalent with solving the score equation (equation 2), where i= (1,…,N), and N is the number of patients in the dataset. The score contributions are the partial derivatives of the log-likelihood with respect to the set of parameters θ evaluated at the N observed data points and the estimated parameters � .

Partitioning of the dataset
After solving the initial Cox model, we will start the partitioning process of the dataset. A partition B of a set D is a set of subsets of D, for which the following conditions are true: B does not include the empty set; the union of all elements of B is D; and the intersection of any two elements of B is empty. Our aim is to find the smallest partitioning B of our dataset D so that in each element of B, the parameters of the Cox model are stable in relation to any parameter from a list of potential partitioning variables.
We will construct the partitioning B in this case recursively. We begin with a partitioning B 0 = {D } consisting of the dataset D alone, a set of parameters θ which is taken from the initial Cox model, and a pre-defined list of potential partitioning variables = ( 1 , … , ). We applied Bonferroni-adjusted permutation test using test statistics of a quadratic form to test the stability of the parameters θ (see details in: Hothorn, Hornik and Zeileis 2006). If some of the variables from Z cause instability (P < 0.05), we choose the one with the strongest association (lowest P-value), denoted here as z j . The values of this variable are then divided into two disjoint sets: for dichotomous variables, this partitioning is trivial, for continuous or non-binary categorical variables the process of finding the partitioning of the variable values is presented in Section 3 of Hothorn, Hornik and Zeileis 2006. The next partitioning B 1 will then consist of two subsets: 10 = { ∈ | ( ) = 0} and 11 = { ∈ | ( ) = 1}. We will then run the Cox model defined above for both datasets of the partition, B 10 and B 11 , separately. In graphical terms, this is the start of the tree: from the root node D we will split the tree into two branches, leading to nodes B 10 and B 11 , forming the next partitioning B 1 = {B 10 , B 11 } Next, the same process will be repeated for both branches separately. If the analyses of both nodes (B 10 and B 11 ) reveal unstable parameters, they are both split into two branches, and the next partitioning B 2 will have four elements; but it is also possible that in one or both nodes no unstable parameters are found which means that the respective branch is no longer split, and the node is determined a leaf node, from which no further splitting will be made. The entire process is repeated until no unstable parameters are found in any node. The collection of all leaf nodes is the final partitioning.

Local centering
The tree model is a powerful method to detect heterogenous treatment effect but still one major problem of studies using observational data remains -adjusting for the selection bias due to the nonrandom assignment of a treatment to a patient.
In randomized experiments, the propensities P(x) = P to receive or not a treatment are constant by definition. The treatment effect can therefore be interpreted as being causal. In contrast, in the case of this study the estimation of a causal effect from observational studies may be biased by the factors that influence the assignment of treatment to each individual. To adjust for this bias and provide more reasonable model results we applied an approach called "local centering" (Robinson 1988; Athey, Tibshirani and Wager 2019; Gao and Hastie 2022; Dandl, Bender and Hothorn 2022. Briefly, the treatment indicator T is centered to T = T 0 -P(x), where T 0 is either 1 (the patient received treatment) or 0 (the patient did not receive treatment) and P(x) is the propensity of receiving treatment based on the covariables x (before the treatment effect is estimated by the tree model. This approach leads to more robustness to confounding effects because it regresses out the effect of the covariables x on the treatment indicator T. We estimated the propensities P(x) by creating a random forest with the cforest founction (package We note that for patients on remdesivir the person-time was split into two parts: the time before (I remdesivir = 0) and after (I remdesivir = 1) starting treatment. These are considered as separate records in the Cox and tree model, and in the propensity calculation. For other drugs, we only have one record per patients, assuming that he/she received treatment since enrollment.
Supplementary Early rdv, remdesivir administered within 7 days of symptom onset; late rdv, remdesivir administered later than 7 days after symptom onset.