2022-zelenko-paper[92250].pdf

submitted to Reviews in Cardiovascular Medicine

Frequentist and Bayesian analysis methods for case series data and application to early outpatient COVID-19 treatment case series of high risk patients Eleftherios Gkioulekas, Ph.D. Professor, School of Mathematical and Statistical Sciences, University of Texas Rio Grande Valley, Edinburg, TX, United States∗

Peter A McCullough, MD, MPH Chief Medical Advisor, Truth for Health Foundation, Tucson AZ, United States†

Vladimir Zelenko, MD Affiliate Physician, Columbia University Irving Medical Center, New York City, NY, United States‡ When confronted with a public health emergency, significant innovative treatment protocols can sometimes be discovered by medical doctors at the front lines based on repurposed medications. We propose a very simple hybrid statistical framework for analyzing the case series of patients treated with such new protocols, that enables a comparison with our prior knowledge of expected outcomes, in the absence of treatment. The goal of the proposed methodology is not to provide a precise measurement of treatment efficacy, but to establish the existence of treatment efficacy, in order to facilitate the binary decision of whether the treatment protocol should be adopted on an emergency basis. The methodology consists of a frequentist component that compares a treatment group against the probability of an adverse outcome in the absence of treatment, and calculates an efficacy threshold that has to be exceeded by this probability, in order to control the corresponding p-value, and reject the null hypothesis. The efficacy threshold is further adjusted with a Bayesian technique, in order to also control the false positive rate. A selection bias threshold is then calculated from the efficacy threshold to control for random selection bias.

Exceeding the efficacy threshold establishes efficacy by the preponderance of evidence, and exceeding the more demanding selection bias threshold establishes efficacy by the clear and convincing evidentiary standard. The combined techniques are applied to case series of high-risk COVID-19 outpatients, that were treated using the early Zelenko protocol and the more enhanced McCullough protocol. The resulting efficacy thresholds are then compared against our prior knowledge of mortality and hospitalization rates of untreated high-risk COVID-19 patients, as reported in the research literature. Keywords: COVID-19; SARS-CoV-2; ambulatory treatment; early treatment, mortality; hospitalization; epidemiology; biostatistics; drug repurposing.

1.

INTRODUCTION

In medical research, the efficacy of new drugs or treatment protocols is established by controlled studies in which a treatment group is compared against a control group. A case series is one half of a controlled study consisting only of the treatment group. At the beginning of the COVID-19 pandemic, practicing medical doctors were confronted with having no treatment to offer to their patients that can prevent or minimize hospitalization and/or death. In response, some doctors were compelled to innovate and discover, on their own, treatment protocols using repurposed off-label medications. Most notable examples, amongst several others, include Didier Raoult [1] in the IHU Méditerranée Infection hospital in Marseilles France, Vladimir Zelenko [2] in upstate New York, George Fareed and Brian Tyson [3] in California, Shankara Chetty [4] in South Africa, Jackie Stone [5] in Zimbabwe, and Paul

∗ Electronic address: drlf@hushmail.com † Electronic address: peteramccullough@gmail.com ‡ Electronic address: zz613@hotmail.com

Marik’s group [6, 7], which was in the beginning based at the Eastern Virginia Medical School. Their efforts to treat patients generated case series of successfully treated patients that constitute Real-World Evidence [8]. The goal of this paper is to present a statistical framework for analyzing systematic case series data of early treatment protocols with binary endpoints (e.g. hospitalization or death), and comparing them against our prior knowledge of the likelihood of adverse outcomes in the absence of treatment. Under certain conditions, which we shall elaborate on below, the proposed methodology can be used to establish the existence of treatment efficacy, but falls short of precisely measuring the corresponding odds ratio. Nevertheless, this can be enough evidence to justify a positive recommendation to adopt these treatment protocols, on an emergency basis, while more detailed clinical research is in progress. Following the recommendation of the American Statistical Association statement on statistical significance and p-values [9], the proposed approach combines use of the p-value, which enables one to reject the null hypothesis, with a Bayesian factor analysis framework [10–14] for controlling the false positive rate [15]. Empirically, we have found that the frequentist p-value framework has done a pretty good job on its own, at least for the

2 analysis of the case series data considered in this paper. However, complementing it with Bayesian factor analysis is a reasonable precaution as it can help raise the red flag when dealing with small sample sizes and/or weak signals. We apply the proposed framework to the processing of available case series data [2, 16–20] that support proposed early outpatient treatment protocols for COVID-19 patients, such as the original Zelenko triple-drug protocol [2] and the more advanced McCullough protocol [21–23]. The original Zelenko protocol was first announced on March 23, 2020 [24]. The proposed approach was to risk-stratify patients into two groups (low-risk vs high-risk), provide supportive care to the low-risk group, and treat the high-risk group with a tripledrug protocol (hydroxychloroquine, azithromycin, zinc sulfate). Zelenko’s risk stratification criteria included in the highrisk group the following: all patients older than 60 were classified as high-risk; all patients of any age with at least one comorbidity or BMI > 30 kg/m2 ; all patients of any age that presented with shortness of breath. Results were reported in an April 28, 2020 letter [16] and a June 14, 2020 letter [17], and the lab-confirmed subset of the April data was published in a formal case-control study [2]. Zelenko’s letters have been attached to our supplementary material document [25].

The rationale for the triple-drug therapy was based on the following mechanisms of action: Hydroxychloroquine prevents the virus from binding with the cells, and also acts as a zinc ionophore that brings the zinc ions inside the cells, which in turn inhibit the RDRP (RNA Dependent RNA Polymerase) enzyme used by the virus to replicate [26, 27]. Azithromycin’s role is to guard against a secondary infection, but we have since learned that it also has its own anti-viral properties [28– 30], and a signal of the efficacy of adding azithromycin on top of hydroxychloroquine can be clearly discerned in a study of nursing home patients in Andorra, Spain [31]. It is interesting that chloroquine was shown in vitro to have antiviral properties against the previous SARS-CoV-1 virus [32], and that there is an anecdotal report from 1918 [33] about the successful use of quinine dihydrochloride injections as an early treatment of the Spanish flu. In hindsight, it is now known that influenza viruses also use the RDRP protein to replicate [34], which can be inhibited with intracellular zinc ions [26, 27]. Consequently, there is a mechanism of action that can explain why we should anticipate the combination of zinc with a zinc ionophore (i.e. hydroxychloroquine, or quercetin [35], or EGCG [36]) to inhibit the replication of the influenza viruses.

Other RNA viruses, including the respiratory syncytial virus (RSV) [37] and the highly pathogenic Marburg and Ebola viruses [38, 39], are also using the RDRP protein to replicate. It remains a compelling hypothesis, deserving further investigation, that the zinc/zinc ionophore concept could play an important role as part of a broader multi-drug treatment or prophylaxis protocol against these serious infectious diseases. Zelenko’s protocol was soon extended into a sequenced multi-drug approach, known as the McCullough protocol [21– 23], which is based on the insight that COVID-19 is a triphasic illness that manifests in three phases: (1) an initial viral replication phase, in which the virus infects cells and

uses them to replicate and make new viral particles, during which patients present with flu-like symptoms; (2) an inflammatory hyper-dysregulated immune-modulatory florid pneumonia, that presents with a cytokine storm, coughing, and shortness of breath, triggered by the toxicity of the spike protein [40], when it is released, as viral particles are destroyed by the immune system, triggering release of interleukin-6 and a wave of cytokines; (3) a thromboembolic phase, during which microscopic blood clots develop in the lungs and the vascular system, causing oxygen desaturation, and very damaging complications that can include embolic stroke, deep vein thrombosis, pulmonary embolism, myocardial injury, heart attacks, and damage to other organs. The rationale of the original Zelenko protocol was that early intervention to stop the initial viral replication phase could prevent the disease from progressing to the second and third phase, and, in doing so, prevent hospitalizations or death. The McCullough protocol [21–23] extends the Zelenko protocol by using multiple drugs in combination sequentially to mitigate each of the three phases of the illness, depending on how they present for each individual patient.

McCullough’s therapeutic recommendations for handling the cytokine injury phase and the thrombosis phase of the COVID-19 illness are, for the most part, standard on-label treatments for treating hyper-inflammation and preventing blood clots. The most noteworthy innovations to the antiviral part of the protocol are the addition of ivermectin, which has 20 mechanisms of action against COVID-19 [41], as an antiviral medication [42–47] to be used as an alternative or in conjunction with hydroxychloroquine, the addition of a nutraceutical bundle [48– 50] combined with a zinc ionophore (quercetin [35] or EGCG [36]) for both low-risk and high-risk patients, and lowering the age threshold for high-risk patients to 50 years. The MATH+ protocol [6, 7], developed for hospitalized patients by Marik’s group, follows the same principles of a sequenced multi-drug treatment. A similar treatment protocol, based on similar insights, was independently discovered and published on May 2020 by Chetty [4] in South Africa. McCullough’s protocol [21–23] was adopted by some treatment centers throughout the United States and overseas, but has not been endorsed by the United States public health agencies, ostensibly due to lack of support of the entire sequenced treatment algorithm by a randomized controlled trial (hereafter RCT).

In spite of the urgent need for safe and effective early outpatient treatment protocols for COVID-19, there has been no attempt to conduct any such trials of any comprehensive multi-drug outpatient treatment protocols throughout the pandemic. Instead, the prevailing approach has been to try to build treatment protocols, one drug at a time, after validating each drug with an RCT. Because COVID-19 is a multifaceted tri-phasic illness, and there is no a priori reason to expect that a single drug alone will work for all 3 phases of the disease, this orthodox approach is not an optimal research strategy. The first priority should be to validate the efficacy of treatment protocols that use multiple drugs in combination, since this is what is actually going to be used in practice to treat patients. On that front, there are published observational data [18–20] on the efficacy of the McCullough protocol, in

3 addition to Zelenko’s reported outcomes [2, 16, 17, 24], and there is an abundance of unpublished real world data from several treatment centers, from around the world, that have yet to be analyzed [51]. The statistical framework proposed in this paper is the missing link for conducting a formal analysis of the available observational data on early COVID-19 treatment protocols. The broader context in which the proposed statistical methodology is situated is as follows. Shortly before COVID19 was declared a pandemic by the World Health Organization, an article [52] was published on February 23, 2020 in the New England Journal of Medicine arguing that “the replacement of randomized trials with non-randomized observational status is a false solution to the serious problem of ensuring that patients receive treatments that are both safe and effective”. The opposing viewpoint was published earlier in 2017 by Frieden [53], highlighting the limitations of RCTs and the need to leverage and overcome the limitations of all available sources of evidence, including real world evidence [8], in order to make lifesaving public health decisions.

In particular, Frieden [53] stressed that the very high cost of RCTs and the long timelines needed for planning, recruiting patients, conducting the study, and publishing it, are limitations that “affect the use of randomized controlled trials for urgent health issues, such as infectious disease outbreaks for which public health decisions must be made quickly on the basis of limited and imperfect data.” Deaton and Cartwright [54] presented the conceptual framework that underlies RCTs and highlighted several limitations. Among them, they have stressed that randomization requires very large samples on both arms of the trial, otherwise, an RCT should not be presumed to be methodologically superior to a corresponding observational study. For example, the randomized controlled trial study of hydroxychloroquine by Dubee et al [55], was administratively stopped after recruiting 250 patients, with 124 in the treatment group and 123 in the control group. Although a two-fold mortality rate reduction was observed by day 28, the study failed to reach statistical significance, due to the small sample size. Even if the mortality rate reduction signal had been strong enough to achieve statistical significance, the small sample size would have still prevented sufficient randomization.

Consequently, although the study has gone through the motions of an RCT, it is not methodologicaly superior to a retrospective observational study. There are several other RCT studies of hydroxychloroquine with similar shortcomings [56]. Furthermore, although a properly conducted RCT has internal validity, in that the inferences are applicable to the specific group of patients that participated in the trial, the external validity of the RCT outcomes needs to be justified conceptually on the basis of prior knowledge, which is either observational, or based on a deeper understanding of the underlying mechanisms of action. Because COVID-19 mortality risk in the absence of early treatment can span three orders of magnitude (from 0.01% to more than 10%) [57–63], depending on age and comorbidities, trials using low-risk patient cohorts are not informative about expected outcomes on the high-risk patient cohorts and vice versa. Likewise, the timing of treatment will

confound the results of an RCT. In general, better results are expected when treatment is initiated earlier, rather than later. As was noted by Risch [64], when interpreting evidence from RCTs, and more broadly from any study, we should bear in mind that results of efficacy or toxicity of a treatment regimen on hospitalized patients cannot be extrapolated to outpatients and vice versa. Likewise, Risch [64] noted that evidence of efficacy or lack of efficacy of a single drug do not necessarily extrapolate to using several drugs in combination. This latter point is further amplified when there is an algorithmic overlay governing, which drugs should be used and when, based on the individual patient’s medical history and ongoing response to treatment. Consequently, RCTs that compare a single drug monotherapy against supportive care are not always informative about whether the drug should be included in a multi-drug protocol. In addition to all that, we are also confronted with an ethical concern. If the available observational evidence are sufficiently convincing, then there is a crossover point where it is no longer ethical to justify randomly refusing treatment to a large number of patients, in order to have a sufficiently large control group.

The corresponding mathematical challenge is being able to quantify the quality of our observational evidence in order to determine whether or not we are already situated beyond this ethical crossover point. Just as the quality of evidence provided by randomized controlled trials is fluid, with respect to successful randomization and external validity, the same is true about the quality of real world evidence [8] that will inevitably become available from the initial scientific mobilization response to an emerging pandemic. The first element of pandemic response is to assess and monitor the situation by prospectively collecting data, needed to construct predictive models of the probability of hospitalization and death, in the absence of treatments that have yet to be discovered, as a function of the patient’s medical profile/history. These models do not necessarily need to be sophisticated at the early stages of pandemic response. It could be sufficient to be able to predict good lower bounds for the hospitalization or mortality probabilities, as opposed to more precise estimates. This early data can be used to identify the predictive factors for hospitalization or death and risk stratify the patients into low risk and high risk categories.

They can also be used as a historical control group that establishes our prior knowledge of expected outcomes, in the absence of treatment, that has yet to be discovered. In parallel with gathering and analyzing data, which is the primary duty and responsibility of our public health and academic institutions, medical doctors have an ethical responsibility to use the emerging scientific understanding of the new disease and it’s mechanisms of actions to try to save the lives of as many patients as possible. Under article 37 of the 2013 Helsinki declaration [65], it is ethically appropriate for physicians to “use an unproven intervention, if in the physician’s judgment it offers hope of saving life, re-establishing health or alleviating suffering”, provided, there is informed consent from the patient, and “where proven interventions do not exist, or other known interventions have been ineffective”. When this effort leads to the discovery of a treatment proto-

4 col, with an empirical signal of benefit and acceptable safety, and using the treatment protocol results in a case series of treated patients, then the confluence of the following conditions makes it possible to statistically establish the existence of treatment efficacy: First, the proposed treatment protocols should use repurposed drugs [66] with a known excellent safety record. When testing new drugs, we have no prior knowledge of the risks involved and a rigorous controlled study is required to determine the balance of risks and benefits. Both Zelenko’s triple drug therapy [2] as well as the expanded multi-drug McCullough protocol [21–23] for the early outpatient treatment of COVID-19 are based exclusively on safe repurposed medications. Second, we need data that give us prior knowledge of the probability risk of the relevant binary endpoints (i.e. hospitalization and/or death) in the absence of treatment, as a function of the relevant stratification parameters. Third, and most importantly, the case series corresponding to treated patients should exhibit a very strong signal of benefit, relative to our prior experience with untreated patients, prior to the discovery of the respective treatment protocol. Under these conditions, the idea that is proposed in this paper works as follows.

Our input is the number N of highrisk patients treated, the number of patients a with an adverse outcome (i.e. hospitalization or death) and selection criteria for extracting the high-risk cohort under consideration, from which we can deduce, based on prior knowledge, that the unknown probability x of an adverse outcome without treatment is bounded by p1 < x < p2 . We also choose the desired level of p-value upper bound p0 , which is typically p0 = 0.05 (95% confidence), although we shall also consider p0 = 0.01 and p0 = 0.001. The output is an efficacy threshold x 0 (N, a, p0 ) that gives us the following rigorous mathematical statement: if x > x 0 (N, a, p0 ), then we have more than 1 − p0 confidence that the treatment is effective relative to the standard of care. This statement has to be paired with the subjective assessment of our prior knowledge, based on which we need to show that x 0 (N, a, p0 ) < p1 . The upper bound p2 is used by the Bayesian factor technique as part of finalizing the calculation of the efficacy threshold x 0 (N, a, p0 ). When there is a large gap between p1 and x 0 (N, a, p0 ), and furthermore, when the treatment relies on repurposed drugs with known excellent safety record, then we have clear and convincing evidence that the treatment is effective.

On the other hand, when new medications, as opposed to repurposed drugs, are introduced into a preexisting treatment protocol, they should be rigorously studied both for safety and efficacy with prospective RCTs. From the sample size N and the finalized efficacy threshold x 0 , we also calculate a selection bias threshold x 1 (N, x 0, p0 ), higher than x 0 , that quantifies how large the gap between p1 and x 0 needs to be, in order to mitigate with 1− p0 confidence, any possible random selection bias in the case-series sample (N, a). As a result, we can assert the existence of treatment efficacy using two distinct standards of evidence. If we can establish p1 > x 0 , then the preponderance of evidence is in favor of the existence of treatment efficacy, and this can justify its provisional adoption on an emergency basis, in order to gather more evidence. If we can establish that x 1 < p1 , then the

evidence becomes clear and convincing, and if these results are replicated by multiple treatment centers, then it becomes ethically questionable to deny patients access to the treatment protocol, for the purpose of conducting an RCT, or simply due to therapeutic nihilism by public health authorities. The paper is organized as follows. On Section 2 we present the technique for calculating the efficacy threshold and the selection bias threshold. We also explain the relationship of the proposed technique with the exact Fisher test and with the binomial proportion confidence interval problem. On Section 3, we present a Bayesian technique for adjusting the efficacy thresholds in order to also control the corresponding false positive rate. In Section 4, we illustrate an application of both techniques to the Zelenko case series [2, 16, 17] as well as the Procter [18, 19] and Raoult [20] case series. Discussion and conclusions are given in Section 5. With the exception of Section 3, which is mainly relevant to a more careful analysis by biostatisticians, we have strived to make Section 2 and Section 4 of the paper relevant and accessible to both clinicians and biostatisticians, by minimizing the mathematical details. Material that is relevant only to biostatisticians is relegated to the appendices.

The computer code and the corresponding calculations are included in the supplementary data document [25].

2.

FREQUENTIST METHODS FOR CASE SERIES ANALYSIS

In this section, we present the technique for comparing a treatment group case series of high-risk patients against the expected probability x of an adverse outcome without treatment, based on prior knowledge. Since our prior knowledge bounds the probability x inside an interval p1 < x < p2 but the precise value of x is unknown, we calculate the minimum value (efficacy threshold) that this probability has to exceed in order to be able to reject the null hypothesis, that the treatment has no efficacy. The proposed technique is equivalent to an exact Fisher test where we take the limit of an infinitely large control group with probability of an adverse outcome set equal to x. We also explain the relationship of the proposed approach with the binomial confidence interval problem, and provide evidence that the corresponding coverage probability is conservative. The assumption of conservative coverage is used, in turn, to derive a selection bias threshold that x should exceed in order reject the possibility of a false positive result due to random selection bias.

2.1.

Comparing treatment group against expected adverse event rate without treatment.

Suppose that we have a treatment group of high-risk patients in which N patients have received treatment, and a patients have had an adverse outcome. Let us also assume that all N patients in the case series satisfy precise selection criteria, used to classify them as high-risk patients, from which we can infer, from our prior knowledge, that in the absence of

5 treatment, the probability x of an adverse outcome is bounded in the interval p1 < x < p2 . To establish the existence of treatment efficacy, we assume the null hypothesis, that the treatment has no effect and that consequently, the probability of an adverse outcome in the treatment group is also equal to x. Under this null hypothesis the probability of observing a patients with an adverse outcome out of a total of N patients is given by ! N a pr(N, a|x) = x (1 − x) N −a, (1) a which corresponds to a binomial distribution. The first factor gives the number of combinations for choosing the a patients that have an adverse outcome out of all N patients. The second factor x a is the probability that the chosen a patients have an adverse outcome, under the assumption of the null hypothesis. The third factor (1 − x) N −a is likewise the probability that the remaining N − a patients will not have an adverse outcome. Consequently, the product of the three factors is the probability of seeing the event (N, a) under the null hypothesis. The corresponding p-value is calculated by adding to the probability of the event (N, a), the probability of all other events with smaller or equal probability, and it reads p(N, a, x) =

N X

pr(N, n|x)H (pr(N, a|x) − pr(N, n|x)),

(2)

n=0

where H is the modified Heaviside function given by ( 1, if x ≥ 0 H (x) = . 0, if x < 0

(3)

The Heavyside function factor in Eq. (2) selects the events (N, n) that are less probable than the observed event (N, a) for inclusion in the probability sum, as per the formal definition of the p-value. In order to reject the null hypothesis, we need to construct a convincing argument that establishes that p(N, a, x) < p0 , with p0 = 0.05 in order to achieve 95% confidence. Such an argument, in effect, is a hypothesis test that compares the treatment group (N, a) outcome against a fixed probability x for an adverse outcome in the absence of treatment. Our proposal for doing so is conceptually very simple. First we calculate an efficacy threshold x 0 (N, a, p0 ) such that x 0 (N, a, p0 ) < x ≤ 1 =⇒ p(N, a, x) < p0 .

(4)

In doing so, we are seeking the smallest possible value of x 0 that satisfies Eq. (4). If our prior belief about x is that it satisfies p1 < x < p2 , then it follows that if we show that our prior belief about the lower bound p1 of the probability x of an adverse outcome without treatment exceeds the efficacy threshold x 0 , then we have a statistically significant signal of benefit in favor of the proposed treatment protocol. This is, in turn, sufficient to recommend to other physicians to consider using the treatment protocol, on an emergency basis, in order to save as many patients as possible, as soon as possible.

We stress again that implicit in this reasoning is the assumption that all observed adverse events in the treatment group case series have been caused by the disease and not by the treatment. For this reason, this methodology has to be limited only to the evaluation of treatment protocols using repurposed medications [66] with previously known acceptable safety. Furthermore, in order to have a prior belief constraining the probability x of an adverse outcome for high-risk patients, in the absence of treatment, it is necessary for public health agencies and academic institutions to prospectively collect data on the predictive factors for hospitalizations or death, as soon as possible at the beginning of an emerging new disease. These data can then be used both to define the selection criteria for identifying patients as high-risk and to constrain the corresponding probability x within the interval p1 < x < p2 .

2.2.

Comments on the proposed hypothesis testing technique

We now make the following comments about the proposed hypothesis testing technique. First, we note that the p-value p(N, a, x), corresponding to a comparison of a case series (N, a) of a treatment group against the probability x of an adverse outcome without treatment, as given by Eq. (2), can be also obtained by running an exact Fisher test with an artificial control group (M, b) of M patients with b adverse outcomes, with x = b/M, in the limit where the size of this artificial control group goes to infinity. In appendix A, we give a mathematical proof of this claim and also explain the mathematically precise for formulation of the statement. This convergence property is in fact, a consequence of a known relationship [67–69] between the hypergeometric distribution, used in the calculation of the exact Fisher test p-value, and the binomial distribution used in the calculation of p(N, a, x). Second, the proposed hypothesis technique is also mathematically related with the well-researched binomial confidence interval problem [70]. Given the case series (N, a) of a treatment group of N patients, with a patients having an adverse outcome, the challenge of the binomial confidence interval problem is to identify a probability interval (q1, q2 ), such that we can assert, with 1 − p0 confidence, that the probability of an adverse event with treatment is inside the interval (q1, q2 ).

If the null hypothesis is satisfied, then the probability of an adverse outcome with treatment is equal to the probability of an adverse event without treatment, and it follows that the intervals (p1, p2 ) and (q1, q2 ) have to intersect. The contrapositive of this deduction is that if the intervals (p1, p2 ) and (q1, q2 ) do not intersect, then the null hypothesis is false. This argument shows that the upper endpoint q2 is the efficacy threshold x 0 (N, a, p0 ) that has to be exceeded by all probabilities in the interval (p1, p2 ) in order to reject the null hypothesis and claim a signal of benefit. More specifically, the method proposed in the preceding section for calculating the efficacy threshold x 0 (N, a, p0 ) is equivalent to calculating the upper endpoint of the Sterne interval [71] for the corresponding binomial confidence interval problem.

6 It is worth noting that although several alternative techniques have been proposed for solving the binomial proportion confidence interval problem, none of them has coverage consistent with the desired statistical confidence and most of them do not have conservative coverage [72]. This means that, given a case series (N, a) for a treatment group, the obtained 95% confidence interval (q1, q2 ) for the probability of an adverse event, with treatment, could be wider or narrower than it should be, depending on the sample size N and the unknown true value of that probability. Furthermore, it has already been proven that no solution to the binomial confidence interval problem exists with perfect coverage [73]. For our purposes, a solution technique with conservative coverage that always overestimates the efficacy threshold x 0 (N, a, p0 ) is acceptable, and to be preferred over techniques that will sometimes overestimate and sometimes underestimate the efficacy threshold. The coverage of a specific solution technique is quantified via the coverage probability c(N, p0 |x), which is defined as the conditional probability of observing a case series (N, a), given a fixed sample size N, for which our solution technique will yield a confidence interval that includes the true probability of an adverse event, under the condition that this true probability is equal to x.

Here, 1 − p0 is the desired level of confidence. For very large sample sizes, c(N, p0 |x) is calculated using computer simulations, however for smaller samples, it can be calculated analytically [74] from the equation, c(N, p0 |x) =

N X

I (N, n, x, p0 ) pr(N, n|x)

(5)

n=0

where I (N, n, x, p0 ) is an indicator function, such that I (N, n, x, p0 ) = 1 if only if x is in the confidence interval obtained by the proposed solution technique for a given case series (N, n) corresponding to 1 − p0 confidence. Otherwise we set I (N, n, x, p0 ) = 0. Conservative coverage requires that our solution technique satisfy c(N, p0 |x) ≥ 1 − p0 for all values 0 ≤ x ≤ 1 of the probability x. The Clopper-Pearson interval [75] is a very well known solution technique to the binomial confidence interval problem that is known to have conservative coverage. However, an efficacy threshold, defined as the upper limit of the ClopperPearson interval, is not equal to what we would have obtained from an exact Fisher test in the limit of an infinite control group, unlike with the Sterne interval [71]. In Fig. 1, we show the coverage probability for the Sterne interval for sample sizes N = 20 and N = 100 and note that it also has conservative coverage, which is very desirable in the context of hypothesis testing. In Fig.

2, we compare the coverage probability of the Clopper-Pearson interval against the coverage probability of the Sterne interval and note that although they are both conservative, the Sterne interval has less conservative coverage probability than the Clopper-Pearson, over the same sample size.This suggests that the Sterne interval is the better choice, both from the standpoint of coverage probability and also due to its relationship with the exact Fisher test. Our third comment concerns the numerical calculation of the efficacy thresholds x 0 (N, a, p0 ) from the function

FIG. 1: Coverage probability for the Sterne interval [71] with sample sizes N = 20 and N = 100. The black curve corresponds to N = 20 and the blue curve, which is situated below the black curve, corresponds to N = 100. The coverage probabilities were calculated using 0.01 increments.

p(N, a, x). To illustrate this calculation with an example, on Fig. 3, we plot the p-value p(N, a, x) against the expected mortality rate x without early outpatient treatment of COVID19, based on Procter’s combined case series [19] of 869 high risk patients that received an early treatment protocol with 2 reported deaths. The figure has vertical lines marking the crossover to 95%, 99% and 99.9% confidence. The corresponding efficacy thresholds are located at the points where the zigzag graph of the function p(N, a, x) intersects with the vertical lines. The discontinuous behavior of p(N, a, x) may seem paradoxical, since we would have expected it to be monotonically decreasing with respect to x, but we have found that it is caused by some of the right-tail contributions to the p-value sum given by Eq. (2). If we make an unwarranted approximation, replacing the right-tail sum with the left-tail sum, we obtain the smooth curve shown in Fig. 3. The intersection points of the smooth curve with the vertical lines, give the upper endpoint of the Clopper-Pearson interval [75]. Similar graphs for all case series considered in the study have been included in our supplementary data document [25].

We have found empirically, at least for the case series being studied here, that both the correct zigzag curve and the approximate smooth curve give almost the same values for all relevant efficacy thresholds.

2.3.

Selection bias mitigation and selection bias thresholds

The idea of hypothesis testing, in which we compare a case series of treated patients against the historical population level (or possibly a more limited) control group, is vulnerable to the

7

FIG. 2: Comparison of the coverage probability for the ClopperPearson interval [75] versus the Sterne interval [71] with sample size N = 100. The black curve shows the coverage probability for the Clopper Pearson interval, and the blue curve, which is situated below the black curve, shows the coverage probability for the Sterne interval. The coverage probabilities were calculated using 0.01 increments on the horizontal axis.

criticism of possible selection bias in the treatment case series in favor of establishing treatment efficacy. Some of the selection bias could be systemic, but even in the absence of any systemic bias, some selection bias will inevitably occur randomly, as a consequence of using a small sample of patients for the treatment case series, randomly chosen out of the general population. We propose the following idea for mitigating random selection bias, and then we discuss the problem of selection bias more broadly. Suppose that we have a case series (N, a), of N treated patients with a adverse outcomes, and suppose that we have calculated the efficacy threshold x 0 for this case series. We can choose x 0 to be either set equal to x 0 (N, a, p0 ), or we can choose to have it further increased, if necessary, using the Bayesian technique of Section 3.

In either case, if we have a prior belief that the probability of an adverse event, without treatment, in the high-risk part of the general population, under the same high-risk patient selection criteria used to form the case series, is equal to x, then, when selecting a random sample of N high-risk patients out of the general population, we can have 1 − p0 confidence that the true rate x 0 of adverse events without treatment for that particular sample will range according to a confidence interval m1 (N, x, p0 )/N ≤ x 0 ≤ m2 (N, x, p0 )/N. Here, m1 (N, x, p0 ) is the minimum number of adverse events and m2 (N, x, p0 ) is the maximum number of adverse events that we expect to see in any one particular sample of N high-risk patients, in the absence of treatment, with confidence 1 − p0 . The possible criticism of our approach is that, perhaps, for the specific sample of patients in our case series, the true adverse event rate x 0, without treatment, could happen to be below the ef-

FIG. 3: Relationship between p-value and expected mortality rate for high risk patients without early treatment, based on the case series data from Procter’s dataset of 869 high-risk patients [19]. The zigzag curve follows p(N, a, x) given by Eq. (2), whereas the smooth curve approximates the right tail terms in the p-value sum by replacing them with the left-tail terms on the horizontal axis.

ficacy threshold x 0 , in spite of the corresponding population level adverse event rate x exceeding the efficacy threshold x 0 . This raises the question of how big does the gap between x and x 0 need to be, to ensure that the entire confidence interval of x 0 lies above the efficacy threshold x 0 ? Since the lower endpoint of this confidence interval is m1 (N, x, p0 )/N, the answer to this question defines a new higher threshold x 1 (N, x 0, p0 ), which we shall call the selection bias threshold such that x 1 (N, x 0, p0 ) < x ≤ 1 =⇒ x 0 < m1 (N, x, p0 )/N

(6)

In Appendix B we prove that this selection bias threshold can be calculated by choosing the smallest possible value of x 1 that satisfies the implication x 1 (N, x 0, p0 ) < x ≤ 1 =⇒ p(N, dx 0 Ne, x) < p0

(7)

Here, the notation dx 0 Ne represents rounding the number x 0 N upwards towards the nearest integer. We note that the calculation of the confidence interval for x 0, as given by Appendix B, uses the assumption that the Sterne interval solution [71] of the binomial confidence interval problem has conservative coverage. Given the selection bias threshold x 1 , and the prior knowledge that in a similarly high-risk cohort at the population level, the probability x of an adverse outcome in the absence of treatment, ranges between p1 ≤ x ≤ p2 , establishing that x 1 < p1 with some gap between x 1 and p1 can be used to rule out random selection bias, with 1 − p0 confidence, as the sole cause of a signal of efficacy. Generally speaking, it is more likely than not that a strong efficacy signal cannot be caused in a particular observation, solely as a result of random selection bias, as long as the center

8 x of the confidence interval for x 0 exceeds the efficacy threshold x 0 . Even if part of the x 0 confidence interval is below x 0 , more than half of the interval will be above x 0 . As a result, the efficacy threshold x 0 and the selection bias threshold x 1 quantify two levels of evidence. Showing x 0 < p1 establishes the existence of treatment efficacy by the preponderance of evidence. Meeting this evidentiary standard should be sufficient for communicating the proposed treatment protocol to other physicians with a caveat that it is still investigational, and that more data is needed before making a definitive claim. Showing that x 1 < p1 establishes the existence of treatment efficacy by the clear and convincing evidentiary standard, and if the gap between x 1 and p1 exceeds one order of magnitude, then we should be able to claim the existence of treatment efficacy beyond reasonable doubt. Our view is that exceeding the selection bias threshold x 1 , for a treatment protocol with acceptable safety, is an objective milestone, beyond which therapeutic nihilism and even the denial of treatment for research purposes becomes unethical.

With regards to the broader problem of systemic bias, there are multiple possibilities to consider: There may be some population-level geographic bias in the patients that live in the geographic area served by a particular treatment center; there may be reporting bias, in that we hear about case series because of the good outcomes, without these outcomes being representative of the actual outcomes at the national or international level; there may be bias in the patient demographics (ratio of low vs high risk patients), and with respect to the timing of treatment (early vs. late treatment). The latter concern can be addressed by stratifying the case series with respect to risk and/or timing of treatment. Geographic bias can be addressed by investigating case series across multiple geographic locations. Outcome reporting bias can be minimized, if we have consecutive case series from the same treatment center, where the initially reported outcomes are replicated by subsequent results. Further mitigation of bias is possible by establishing a large gap between the selection bias threshold and the lower bound for adverse outcomes in the historical control statistics. Last, but not least, statistical quantitative evidence can be corroborated and amplified with more qualitative evidence based on the Bradford Hill criteria [76].

Particularly relevant, in that regard, are: (1) plausibility, i.e. the existence of a known biological mechanism of action that explains why the treatment protocol is expected to work; (2) consistency, i.e. observing the same effect in different treatment centers in different locations; (3), biological gradient, i.e. observing improved outcomes with increased medication dosage or length of treatment, additional medications, or by initiating treatment earlier, rather than later. (4) temporality, i.e. immediate improvement in symptoms, following the administration of the treatment protocol. The statistical evidence alone, speak in support of only the strength of association criterion, but that is only one of the several criteria proposed by Bradford Hill [76]. If we can establish that these additional criteria are satisfied, then that constitutes additional evidence on top of the statistical evidence, that treatment efficacy exists, and that the signal of benefit

cannot be explained away, as a result solely caused by selection bias.

3.

BAYESIAN FACTOR ANALYSIS OF EFFICACY THRESHOLDS

The methodology that we proposed in Section 2 is also vulnerable to the criticism that rejecting the null hypothesis, solely on the basis that the p-value satisfies p < 0.05, is not sufficient for asserting that treatment efficacy is statistically significant. This is indeed the position of the recent American Statistical Association statement on statistical significance and p-values [9]. The problem is that p-values only measure how incompatible the data are with the null hypothesis. However, this measure does not always do a good job at controlling the probability of a false positive result [77]. To estimate the latter probability we would have to formulate the appropriate alternate hypothesis and consider how much the data is compatible or incompatible with that alternate hypothesis. This has prompted recommendations to lower the p-value threshold down to 0.01 or 0.001 [77, 78]. However, this is only a stopgap measure that does not fundamentally address the problem. In this section, we supplement the p-value based analysis of Section 2, with a proposal for a Bayesian factor analysis [10– 14]. The Bayesian factor compares the alternate hypothesis (treatment efficacy), against the null hypothesis and can be used to calculate the probability of a false positive result [15]. We do not mean to suggest that the Bayesian factor should replace the p-value in hypothesis testing.

Our view is that we need to use both. That is, use the p-value to reject the null hypothesis, and then use the Bayesian factor to assess the strength of the evidence in favor of the alternate hypothesis. This viewpoint is similar to earlier proposals for conditional frequentist testing [10]. In the following, we will briefly review the Bayesian factor framework and then outline our specific proposal for validating and adjusting, as needed, the efficacy threshold x 0 (N, a, p0 ). We note that the calculation of the selection bias threshold x 1 is independent of the technique used to calculate the efficacy threshold x 0 . In terms of procedure, one could initially calculate the efficacy threshold x 0 using only the technique of section 2 and then use that to calculate the corresponding selection bias threshold x 1 . Alternatively, a more detailed analysis would involve first calculating the efficacy threshold using the technique of section 2, then adjusting it using the technique presented in this section, and then using the adjusted efficacy threshold x 0 to calculate the corresponding selection bias threshold x 1 .

3.1.

Bayesian factor and the false positive rate

Let ,A, B, be two arbitrary events in some probability space. From the definition of conditional probability, we obtain the

9 Bayes rule, noting that p( A ∩ B) p(B ∩ A) = p( A) p( A) p( A|B)p(B) . = p( A)

p(B| A) =

(8) (9)

Let D represent our data, H0 represent the null hypothesis, and H1 represent the alternate hypothesis. In the Bayesian statistics framework, we assign probabilities p(H0 ), p(H1 ) to the hypotheses H0, H1 representing our prior belief about how likely each hypothesis is, and then calculate the updated probabilities p(H0 |D) and p(H1 |D) on the condition of observing the data D. In this way, Bayesian statistics is distinct from frequentist statistics where probabilities are not assigned to the hypotheses themselves. From the Bayes rule we have, p(D|H1 )p(H1 ) , p(D) p(D|H0 )p(H0 ) p(H0 |D) = , p(D)

p(H1 |D) =

(10) (11)

(12)

The Bayes factor B(D|H1, H0 ) is defined to read B(D|H1, H0 ) =

p(D|H1 ) , p(D|H0 )

(13)

and it is the numerical factor that amplifies our prior belief about the odds ratio b(H1, H0 ) = p(H1 )/p(H0 ) after seeing the data D. Here, p(D|H1 ) is the probability of seeing the data D if H1 is true and p(D|H0 ) is likewise the probability of seeing the data D if H0 is true. To interpret the meaning of the Bayesian factor, the following argument is used to calculate the posterior probabilities p(H1 |D) and p(H0 |D) in terms of B(D|H1, H0 ) and b(H1, H0 ) = p(H1 )/p(H0 ). We assume that H0, H1 satisfy p(H0 ) + p(H1 ) = 1 and p(H0 |D) + p(H1 |D) = 1. Combining the second equation with Eq. (10) and Eq. (11) gives the Bayes theorem p(D) = p(D|H0 )p(H0 ) + p(D|H1 )p(H1 ),

p(H0 |D) < p0 ⇐⇒ b(H1, H0 ) >

(14)

and it follows that the probability of a false positive result is given by p(D|H0 )p(H0 ) (15) p(D) p(D|H0 )p(H0 ) = (16) p(D|H0 )p(H0 ) + p(D|H1 )p(H1 ) p(D|H0 )p(H0 ) = p(D|H0 )p(H0 )[1 + B(D|H1, H0 )b(H1, H0 )] (17) 1 = . (18) 1 + B(D|H1, H0 )b(H1, H0 )

1 − p0 , p0 B(D|H1, H0 )

(19)

which relates an upper bound p0 on the probability p(H0 |D) with a corresponding lower bound bmin (p0, B) on the prior likelihood ratio b(H1, H0 ), which is given by bmin (p0, B) =

and dividing the two equations gives p(H1 |D) p(D|H1 ) p(H1 ) = . p(H0 |D) p(D|H0 ) p(H0 )

We see that the false positive probability approximately scales as the inverse of the Bayes factor B(D|H1, H0 ). On the other hand, the dependence of p(H0 |D)) on the prior likelihood ratio b(H1, H0 ), which measures our subjective belief about the odds ratio between H1 and H0 , before seeing the data D, is uncomfortable. There are three ways to cope with that: First, one can simply join the frequentist camp, consider probabilities based on beliefs as meaningless, and forget about the whole thing. Second, one can use an uninformed prior, meaning that we assume that both hypotheses H0 and H1 are equally probable, not having any prior knowledge that favors one over the other, and choose p(H0 ) = p(H1 ) = 1/2, which corresponds to b(H1, H0 ) = 1. An interesting third way is to use the reverse Bayesian analysis technique proposed by Colquhoun [15], which is based on the equivalence

1 − p0 , p0 B

(20)

with B being the value of the corresponding Bayesian factor. The meaning of Eq. (20) is that, given a desired lower bound p0 for the false positive rate and a threshold B for the Bayesian coefficient, bmin (p0, B) is the minimum prior likelihood ratio p(H1 )/p(H0 ) for our prior knowledge of the extent to which the alternate hypothesis H1 is favored over the null hypothesis H0 , for which the Bayesian threshold B can control the false positive rate and keep it below p0 . As such, given our subjective choice for bmin , one can calculate the threshold B for the Bayesian factor corresponding to the minimum tolerated false positive rate p0 . Since we wish to constrain the false positive rates to less than 0.05, in order to claim 95% statistical significance, we choose p0 = 0.05. Kass and Raftery [14] and Jeffries [79] both recommend that the threshold B > 100 be used for a decisive acceptance of the alternate hypothesis H1 over the null hypothesis H0 . Using B = 100 we find that bmin (0.05, 100) = 0.19. This means that if we associate the decisive threshold B > 100 with 95% confidence, doing so is equivalent to a prior belief that the null hypothesis is 5 times more likely than the alternate hypothesis. In turn, this prior belief can be used to deduce Bayesian factor thresholds for higher levels of confidence, consistently with our choice to associate B > 100 with 95% confidence.

p(H0 |D) =

3.2.

Application to hypothesis testing for case series

Now, let us consider how Bayesian factor analysis can be applied to a case series with a treatment group of N patients, where a patients have an adverse outcome. Let x 0 be the corresponding efficacy threshold, determined via the techniques of Section 2, and let x be the probability of an adverse outcome with treatment. We define a null hypothesis H0 and an

10 alternate hypothesis H1 about the value of x such that H0 : x 0 < x ≤ 1, H1 : 0 < x ≤ x 0 .

(21) (22)

We use for x 0 the upper endpoint of the binomial proportion confidence interval corresponding to the observed data (N, a). Consequently, the null hypothesis H0 has been defined to place x outside and above that interval, and the alternative hypothesis H1 considers the remaining possible values for x. Because both H0 and H1 are composite hypotheses, it is necessary to introduce prior probabilities pr(x|H0 ) and pr(x|H1 ), corresponding to H0 and H1 . It may seem tempting to just use uninformed priors, both for H0 and H1 , however, doing so would certainly not be appropriate for the null hypothesis H0 in almost all situations, since with many illnesses, we can rule out probabilities of adverse outcome beyond some upper bound p2 . We can thus use instead an uninformed prior on the interval [x 0, p2 ], given by ( 1/(p2 − x 0 ), if x ∈ [x 0, p2 ] pr(x|H0 (x 0, p2 )) = (23) 0, if x ∈ (p2, 1], and perform an appropriate sensitivity analysis on the parameter p2 . In general, increasing p2 will tend to increase the Bayes factor, since doing so will tend to increase the contrast between the null and alternate hypotheses. So we can explore how much p2 can be decreased and still maintain a decisive Bayes factor.

Likewise, for the alternate hypothesis H1 , we will use an uninformed prior on the interval [0, t] with t ≤ x 0 given by ( 1/t, if x ∈ [0, t] pr(x|H1 (x 0, t)) = (24) 0, if x ∈ (t, x 0 ]. The reason for this choice is that we have found empirically that in some cases, the Bayes factor may actually increase, if instead of an uninformed prior on [0, x 0 ] we use an uninformed prior on the shorter interval [0, t]. From an intuitive standpoint, we surmise that if the data has a very strong efficacy signal, then the contrast between the null and alternate hypotheses is increased when one eliminates the relatively unlikely values of x between t and x 0 . For this reason, we shall use the maximum value of the Bayes factor taken over all values t ∈ (0, x 0 ), on a decimal logarithmic scale, which is given by b(x 0, p2 ) = max b0 (x 0, p2, t),

(25)

b0 (x 0, p2, t) = log B(N, a|H1 (x 0, t), H0 (x 0, p2 )).

(26)

t ∈(0,x0 ]

In appendix C we prove that the function b0 (x 0, p2, t) is initially increasing and then decreasing with respect to t with a maximum in the interval [a/N, 1]. If this maximum is located in the narrower interval [a/N, x 0 ] then the optimal Bayes factor is indeed obtained when we use a choice t ∈ (0, p0 ) for the prior distribution of the alternate hypothesis H1 . If the maximum is formally located at t > x 0 , then the optimal Bayes factor is obtained at t = x 0 . The resulting metric b(x 0, p2 ) is

still dependent on the parameter p2 of the prior distribution of the null hypothesis H0 . To complete the metric definition by Eq. (25) and Eq. (26), we now show the calculation of the Bayes factor B(N, a|H1 (x 0, t), H0 (x 0, p2 )) between H1 and H0 as of function of x 0, p2, t and the data N, a. We note that the probabilities for seeing the data (N, a) under the hypotheses H1 and H0 are given by: Z 1 pr(N, a|H0 (x 0, p2 )) = dx pr(N, a|x) pr(x|H0 (x 0, p2 )) x0

(27) = =

1 p2 − x 0

Z p2 dx pr(N, a|x)

(28)

x0

1 N p2 − x 0 a

! Z p2

x a (1 − x) N −a dx,

(29)

x0

and pr(N, a|H1 (x 0, p2 )) =

Z x0

dx pr(N, a|x) pr(x|H1 (x 0, t))

0

(30) Z t

1 dx pr(N, a|x) t 0 !Z t 1 N = x a (1 − x) N −a dx, t a 0

=

(31) (32)

consequently, the corresponding Bayes factor is given by pr(N, a|H1 (x 0, p2 )) B(N, a|H1 (x 0, t), H0 (x 0, p2 )) = pr(N, a|H0 (x 0, p2 )) Z t x a (1 − x) N −a dx p2 − x 0 0 Z p2 = . t x a (1 − x) N −a dx

(33)

(34)

x0

The integrals can be calculated using exact algebra or numerically with the open source computer algebra software Maxima [80]. The exact algebra calculation takes longer to carry out, but we have confirmed that the numerical calculation using the function quad qagr is just as accurate. In order to control for the false positive rate, we propose that the efficacy thresholds x 0 (N, a, p0 ) with p0 = 0.05 should be increased, if necessary, by requiring that they also satisfy b(x 0, p2 ) ≥ 2. Since the threshold used for a decisive Bayes factor with p0 = 0.05 corresponds approximately to bmin (p0, B) = 1/5, it is reasonable to use the empirical formula ! 5(1 − p0 ) b(x 0, p2 ) ≥ log , (35) p0 to adjust the efficacy thresholds x 0 (N, a, p0 ) for an arbitrary value of demanded confidence p0 . For p0 = 0.01, this gives b(x 0, p2 ) ≥ 2.7 and for p0 = 0.001 we find b(x 0, p2 ) ≥ 3.7 as the Bayes factor thresholds corresponding to a prior likelihood ratio p(H1 )/p(H0 ) = 1/5 and, as such, they are the thresholds that we recommend imposing on the Bayes factors for the purpose of adjusting the corresponding efficacy thresholds x 0 (N, a, p0 ) for the choices p0 = 0.01 and p0 = 0.001.

11 4. APPLICATION TO THE ANALYSIS OF EARLY OUTPATIENT COVID-19 TREATMENT CASE SERIES 4.1.

Review of the Zelenko, Procter and Raoult case series

We shall now analyze, using the aforementioned methodology, the high-risk patient case series by Zelenko [2, 16, 17], Procter [18, 19], and Raoult [20]. The main reason for focusing on these case series specifically, is that they consist exclusively of high-risk patients, where early outpatient treatment is expected to make a difference. In this subsection we shall present the details of these case series, focusing on outcomes and treatment protocols used. In the Zelenko April 2020 letter [16], Zelenko reported on his outcomes based on a total of 1,450 patients that he treated for COVID-19 until April 28, 2020 in an Orthodox Jewish community in upstate New York. From this cohort, 405 patients were classified as high risk and treated with his triple-drug therapy (hydroxychloroquine, azithromycin, zinc sulfate). The reported outcomes were 6 hospitalizations and 2 deaths. From amongst the patients classified as low risk, who were given only supportive care, there were no hospitalizations or deaths. Zelenko’s criteria for risk stratification define three categories of high risk patients: (1) every patient older than 60; (2) every patient younger than 60 but with comorbidities; (3) patients younger than 60 and without comorbidities that presented with shortness of breath.

A subset of the April 28, 2020 case series was published in a case controlled study [2] that included only the treated patients with COVID-19 infection that was confirmed by a PCR test or an antibody IgG test. The remaining patients were clinically diagnosed from symptomatic presentation and via ruling out a bacterial or influenza infection. This Derwand– Scholz–Zelenko study (hereafter DSZ study) [2] included 335 patients of which, 141 patients were classified as high-risk patients and treated with the triple drug protocol with 4 hospitalizations and 1 death. Detailed demographic data is given for the high-risk patient treatment group, including a detailed breakdown in the three high-risk categories. The study also included a control group of 377 patients who were seen by other treatment centers in the same community, that were only offered supportive care and no early outpatient treatment. From this untreated group, 13 patients died and 58 patients were hospitalized. The untreated group includes both low-risk and high-risk patients, so we expect that it underestimates both the hospitalization and mortality risk for high-risk patients.

Unfortunately, demographic data was not available for the untreated group, so from a strictly methodological point of view, one cannot entirely rule out the theoretical possibility that the untreated group might have consisted of patients that are at higher risk on average than those of the high-risk treatment group. On the other hand, using a case series of untreated patients from Israel [60], with demographic data indicating a combination of low and high-risk patients, with 143 deaths reported out of 4,179 untreated patients, gives the same mortality rate as in the DSZ control group, suggesting that the DSZ control group also consists of a mixed demographic of low and high risk patients.