Ranking algorithms order items according to the belief that they possess some desired feature. When the presence/absence of the desired feature is known, ranking algorithms are often evaluated by “testing” items according to the relative rank or testing order. Ideally, all early tests reveal the desired feature. The statistics and machine learning communities use a number of performance curve variants to evaluate ranking algorithms, including the hit enrichment curve and the enrichment factor or lift curve. Popular software such as the R package caret (Kuhn 2008), SAS Enterprise Miner (SAS Institute Inc. 2020), and JMP (SAS Institute Inc. 2020) can be used to construct these curves. Some curves are used extensively in the evaluation of virtual screens of chemical compounds for drug discovery (Geppert et al. 2010). They are also used in a number of other applications such as the evaluation of marketing campaigns (Rosset et al. 2001).

In the context of virtual screening, the desired feature is often a biological activity. Typically the desired activity is binding to a protein target, so we will refer to the chemical compounds as ligands. Ligands are scored, where the scores are provided by one of the many ranking algorithms available, such as molecular docking algorithms, pharmacophore models, or quantitative structure-activity relationship (QSAR) models.

Empereur-Mot et al. (2016) recently provided a web application for constructing performance curves to evaluate virtual screening methods. One of the performance curves they provide is the hit enrichment curve. The software provides many nice interactive features for exploring performance metrics at points along a curve. It also provides utilities for combining the scores from multiple methods into a consensus score that may improve ranking performance.

Of tremendous importance in a world where virtual screening data is tightly guarded, Empereur-Mot et al. (2016) make some of their data freely available. We use one of the case studies as a demonstration dataset. The target was the protein regulating gene peroxisome proliferator-activated receptor gamma (PPARg), which has been linked to several diseases such as obesity, diabetes, atherosclerosis, and cancer (NCBI 2021). While this is a small dataset for a retropsective screening evaluation by present-day standards, the rarity of actives (ˆ0265 is the observed fraction of active ligands) and the small testing fractions are very consistent with larger screening studies (Zhu et al. 2013). The primary goal in such a study is to discover scoring methods that are able to correctly identify active ligands, and to do so very early in the testing phase (Truchon and Bayly 2007; Jain and Nicholls 2008; Geppert et al. 2010).

The hit enrichment curve is commonly used to summarize effectiveness of a screening campaign. It plots the proportion of active ligands identified (i.e., the recall) as a function of the fraction of ligands tested, where testing order is based on the score produced by a virtual screening method. Larger recall values are preferred and, going a step further, these are more relevant when they occur at small testing fractions. In other words, one hopes to demonstrate improvement in early enrichment.

While the hit enrichment curve is designed to show results for all testing fractions, it is common to focus on fractions below 0.1 and even below 0.001 (Zhu et al. 2013). The size of our demonstration dataset limits us to consider 0.001 as the smallest testing fraction (resulting in just three ligands tested), but present-day screening campaigns can easily involve millions of ligands and hence testing fractions below 0.001 are reasonably considered. For a given testing fraction, one may want to compare observed recall values for competing scoring methods.

For the PPARg study, Empereur-Mot et al. (2016) use three popular docking methods to score ligands: Surflex-dock, ICM, and Vina. We scale all docking scores so that a larger value is consistent with active ligands; as a result, ICM and Vina scores have been negated. Empereur-Mot et al. (2016) also evaluate several methods for constructing consensus scores, and we limit investigations to two of their best performers: maximum of the z-scores from Surflex-dock and ICM, and the minimum of the ranks from Surflex-dock and ICM.

Figure 1 shows estimated hit enrichment curves for the five scoring methods over the full range of testing fractions, along with the ideal hit enrichment curve that would be expected if the 85 active ligands were identified in the first 85 tests, and the hit enrichment curve consistent with random identification of active ligands. This figure differs slightly from that produced by the web application of Empereur-Mot et al. (2016) due to a difference in how ties are handled; we use inverse distribution functions (strategy discussed below) to avoid arbitrary indications of better performance due simply to random ordering.

Figure 1: Hit enrichment curves for the PPARg application, comparing five scoring methods over the full range of testing fractions. Ideal and random hit enrichment curves are also shown.

Comparing consensus methods to individual docking methods, is the observed improvement in recall at testing fraction 0.1 significant? With recall being the proportion of active ligands identified at this testing fraction of 321 tests, it may be tempting to apply a procedure based on comparing independent binomial proportions. There are, however, two complications. The first is that determination of testing order requires scores from all ligands, and hence this introduces correlation between the 321 tests that are applied for a single scoring method. The second complication is that competing scoring methods are very likely positively correlated and so the uncertainty associated with differences between hit enrichment curves may possibly be reduced, thus improving the power to uncover differences. This paper develops appropriate techniques for comparing hit enrichment curves that result from competing scoring methods. Our inferential procedures trivially extend to enrichment factors, as well. The need for proper inferential methods for these early enrichment metrics has been emphasized by many recent papers (Nicholls 2014; Robinson et al. 2020).

Section 2 introduces hit enrichment curves constructed from ranking algorithms.Section 3 presents four approaches for hypothesis testing and confidence intervals to compare two hit enrichment curves, along with simulation studies to compare effectiveness of the approaches; EmProc is newly proposed here, while three other approaches are applied in new ways. Section 4.1 presents confidence band procedures and simulation results for an entire hit enrichment curve from a single algorithm, while Section 4.2 considers bands for the difference between two hit enrichment curves. Section 5 revisits the PPARg application. Section 6 includes a discussion of general findings from this study, and also makes connections to the broader task of assessment of virtual screening campaigns.

Let S denote the score from a ranking algorithm, where larger values of S suggest stronger belief that a ligand is active. S is reasonably regarded as a random variable. Activity of a ligand may also be regarded as a random variable: indicator function. That is, X = 1 when a ligand belongs to the active class (+) and X = 0 when a ligand belongs to the inactive class (. Given that a ligand is active, S has cumulative distribution function ), and given that a ligand is not active, S has cumulative distribution function ). Combining the scores from both classes results in the mixture distribution

Once ligands have been ranked according to their score, S, a threshold on the score, t, will prioritize a top fraction of the data set for testing. Let this top fraction be r = P(S > t), which is the x-axis of the population hit enrichment curve. The hit enrichment curve (also known as the enrichment curve, accumulation curve, or percent captured response curve) is often used when an entire curve is used to evaluate a virtual screening campaign. Population hit enrichment curves plot P(S > t|+) on the y-axis, where P(S > t|+) is known as recall at threshold t.

So far, we have described the population level distributions of random variables X and S, that is, the distribution of ligand activity for the population of drug candidates from which a data set is sampled and the distribution of scores that a ranking algorithm would assign to them. We consider a data set under examination to be a random sample of activity and score pairs from these population level distributions. Let be the scores that were sampled from the + class mixture component, be the scores that were sampled from the class mixture component,

The empirical hit enrichment curve plots the cumulative fraction of actives on the y-axis, identified as a function of the top r fraction of ranked ligands. This means all compounds with scores beyond the percentile 100(1 ) are “tested” and the cumulative fraction of actives determined. Another way of determining this percentile would be to choose a threshold ˆsuch that the fraction of items with instead of t to denote that this threshold defines a fraction of the sample data and not the population.

Specifically, we define, ˆ) to be the empirical cumulative distribution functions (cdfs) for all and + scores:

Given a testing fraction r, we ideally choose a threshold ˆto be the score percentile selected for testing such that ). To accommodate the possible existence of ties in the observed data on scores, we define ˆ. Estimated recall at testing fraction r is the fraction of the active compounds that are correctly predicted to be active (i.e. ) and is thus obtained as )). Thus, the empirical hit enrichment curve plots the pairs

Jiang and Zhao (2015) have shown that if 1), then the empirical hit enrichment curve is an unbiased estimator of the population hit enrichment curve.

We wish to determine whether one ranking algorithm has significantly better performance than another at a given testing fraction denote the true population-level recall for a ranking algorithm at testing fraction r. Let be a random sample from the ligand activity distribution and the score distributions of ranking algorithm 1 and 2. We assume that the triplets (are independent across are likely correlated. The amount of correlation will depend on the extent to which the scores are measures of the same characteristics of the ligands. For example, competing docking scoring functions are often parameterized in similar ways (e.g., Glide SP and Glide XP; see Friesner et al. (2004)), and competing QSAR models often utilize highly correlated sets of descriptors. Using the random sample, we estimate the difference in performance between two algorithms at a given and perform a hypothesis test to determine if the difference is significant.

First considering a single algorithm, let ) represent the number of active ligands that are examined at testing fraction r. We estimate recall at testing fraction , noting that both the numerator and denominator are random variables. The activity rate is reasonably estimated as ˆan alternative expression for estimated recall is ˆ). Because ˆis estimated using the entire dataset through the empirical cdf ˆis not binomially distributed. To properly account for this, Jiang and Zhao (2015) take an empirical process approach to derive asymptotic normality of ˆ. Their result is that

as , with corresponding asymptotic variance expression

where

is the simple binomial variance, and Λ) is a threshold-specific activity rate. This result assumes that 0, and that the conditional densities are positive and continuously differentiable in a neighborhood of ; henceforth called Conditions 1 and 2.

When it comes to comparing estimated recall across two competing algorithms, there are two sources of correlation that should be addressed. One source is correlation induced by needing to estimate , and that is addressed by using the result of Jiang and Zhao (2015). The other source of correlation arises because competing algorithm scores are derived using some common data, and this source of correlation has not been previously addressed in the literature. By accounting for both types of correlation, we expect to

improve the power to detect real differences in algorithmic performance.

The following subsections present four methods of testing for significant differences in recall for two competing algorithms. First, we extend the approach of Jiang and Zhao (2015) to use an empirical process approach that accounts for the correlation between algorithms, in addition to the correlation induced within a particular algorithm by estimation of this method is called EmProc. Second, we present details on how a McNemar procedure for correlated proportions may be applied to the application; as far as we know, this has not been previously done. Third, we apply the Jiang and Zhao (2015) result for hypothesis testing; this method is called IndJZ and is only optimized to address correlation within each algorithm but not correlation between algorithms. And fourth, we treat as correlated binomial random variables, thus ignoring correlations induced by estimating but accounting for correlations between algorithms; this method is called CorrBinom. The four methods are compared using a simulation study in Section 3.5.

All four methods are based on asymptotic normality of a test statistic of the form

We reject 05 level test.

Pointwise confidence intervals are obtained as

The methods differ in their approach to estimating

3.1 EmProc: Adjust for correlation between algorithms and correlation within each algorithm

Taking an empirical process approach, the functional delta method employed by Jiang and Zhao (2015) was extended to derive the following asymptotic normality result concerning

Theorem 1: Given that Conditions 1 and 2 are satisfied for both

as . Furthermore, the asymptotic variance expression is

where: ) is as given in equation (1) and applied for each algorithm;

is the covariance between binomial counts; ) and determines the threshold for algorithm +) is the conditional probability that both algorithms result in testing an active ligand because it is highly ranked by both algorithms; ) is the unconditional probability that a ligand is highly ranked by both algorithms; and Λ2. Details and derivations are in the Appendix.

To estimate ), we replace the population parameters with consistent estimators. The frequency distribution for the tested/not-tested status according to both ranking algorithms for the active ligands is shown in Table 1, where counts the number of active ligands tested by algorithm j for j = 1, 2, and counts the number of active ligands tested by both algorithms. As previously discussed, we estimate tional estimates are obtained as and is obtained using Nadaraya-Watson kernel regression with the “rule-of-thumb” bandwidth selector (Li and Racine 2007).

Table 1: Frequency distribution for the tested/not-tested status of active ligands according to both ranking algorithms, at a fixed testing fraction r.

Under the null hypothesis , we could alternatively use a pooled estimator of ) to replace both in variance expres- sion (4). We consider variance estimates using both the unpooled and pooled approaches.

3.2 McNemar’s test for difference in recall

When estimating recall, the same set of active ligands simultaneously serves as the set of “trials” for both ranking algorithms, with the decision being whether or not each algorithm selects the active ligand for testing. The consequence is that the data for both algorithms may be viewed as fully paired. The standard test used for paired proportions is McNemar’s test (Agresti 2007; Fagerland et al. 2013). Table 1 shows how the number of active ligands tested by either ranking algorithm can be written as a 2 2 contingency table. The estimated recall values present themselves as the marginal probabilities of testing an active ligand for each ranking algorithm. Consequently, the asymptotic McNemar test is based

on test statistic

and it assumes that discordant counts () are large. In a

simulation study comparing the asymptotic McNemar test to a number of other tests for paired nominal data, the asymptotic McNemar test was found to be the most powerful across simulation scenarios, though slightly liberal in terms of type I error (Fagerland et al. 2013).

While the asymptotic McNemar test enforces the null condition to replace with zero in the variance expression, an alternative approach is needed to obtain pointwise confidence intervals. Pointwise Wald confidence intervals use the following standard error expression in equation (3):

Unfortunately, several studies (Newcombe 1998; Fagerland et al. 2013; Rodriguez de Gil et al. 2013) demonstrate inadequate coverage properties of the Wald interval. The BonettPrice (Bonett and Price 2012) adjusted interval is a simple modification of the Wald interval, and it has been shown to have good coverage properties (Fagerland et al. 2013; Rodriguez de Gil et al. 2013). We refer to the Bonett-Price adjustment as a “plus” adjustment because it adds one unit to each of the discordant counts shown in Table 1, then applies the Wald formula. More precisely, discordant count +1 and discordant count + 1, thus adding one to each of the marginal counts

and two to the overall total. As a result, the Bonett-Price plus interval is

noting that both the center point and the standard error have been adjusted.

3.3 IndJZ: Adjust for correlation within each algorithm but not correlation between algorithms

If we assume that are independent, then ) is obtained as in equation (1) for j = 1, 2. For testing equality of recall for the algorithms, either a pooled or unpooled estimator of the variance could be used, as previously discussed.

When the competing algorithms have scores that are highly positively correlated, it is expected that the IndJZ approach will lead to standard errors that are unnecessarily large, resulting in an underpowered test.

3.4 CorrBinom: Adjust for correlation between algorithms but not correlation within each algorithm

In this approach, we treat as if it follows a simple binomial distribution, even though it does not. As a result, the relevant variance expression is ) is defined in equation (2) and defined in equation (5).

For testing equality of recall for the algorithms, either a pooled or unpooled estimator of the variance could be used, as previously discussed.

3.5 Simulation Results

3.5.1 Simulation of Benchmarking Data Sets

Geppert et al. (2010) and Xia et al. (2015) have recently reviewed the standard data sets used to benchmark virtual screening tools. The goal in designing a benchmark data set is to mimic real world chemical collections – this means that the score and activity distributions of the compounds in the benchmark should resemble these populations.

Prospective virtual screens are typically conducted on large databases such as ZINC (Sterling and Irwin 2015). These databases can be considered random samples from “drug like” chemical space. The performance curves estimated by retrospective virtual screens will likely mis-estimate performance in a prospective virtual screen if the benchmark data set is not also a random sample. To this date, directory of useful decoys (DUD) is the most widely used collection of benchmarking data sets used in the evaluation of retrospective virtual screens, however each of the biases mentioned have been observed in these data sets (Irwin 2008; Good and Oprea 2008; Xia et al. 2015). The directory of useful decoys, enhanced (DUD-E) data sets (Mysinger et al. 2012) were developed to address some biases in these data sets, though the data sets still lack the experimental testing of decoys (i.e., there is still false negative bias) and there is an unrealistic frequency of actives included in each of the data sets. The MUV benchmarking data sets (Rohrer and Baumann 2009) have also been developed with the intention of minimizing these biases. A clear advantage over DUD-E is that the decoys in MUV have been tested experimentally. The authors of MUV collected 18 primary high-throughput screen assays from PCBioAssay (NCBI Resource Coordinators. 2016). Actives were further confirmed with low throughput assays to minimize the number of false positives, and additional checks for false negatives were performed. We modeled our simulations on the MUV benchmarking data sets, because we believe these data sets to be the most representative of the population of drug candidates.

Basing our simulations on MUV, we simulated data sets with and skew = 499 (where skew is ratio of inactives to actives, (1 ). In general, this is representative of a typical virtual screening data set with large sample sizes and extreme class imbalance.

3.5.2 Study Design

A study was conducted to investigate the power of EmProc, McNemar, IndJZ, and CorrBinom. For each of EmProc, IndJZ, and CorrBinom, a total of four modifications were considered according to whether pooling was applied or whether the Bonett-Price plus adjustment was applied. Although the Bonett-Price plus adjustment was proposed for the Wald confidence interval to improve coverage probability, we wondered if it might also be used for hypothesis testing. McNemar was not modified. In fact, McNemar is identical to CorrBinom with pooling and without plus adjustment.

A study was also conducted to investigate coverage probabilities of the confidence intervals associated with EmProc, IndJZ, CorrBinom, and the Bonett-Price plus-adjusted Wald interval (the latter being referred to as the McNemar inverval, for brevity). There is no justification for pooling when considering confidence intervals, so no pooling is applied. We continue to investigate the impact of the Bonett-Price plus adjustment. The McNemar interval is identical to the CorrBinom interval with plus adjustment.

Four data-generating mechanisms were considered: binormal or bibeta; and correlation of 0.9 or 0.1. In the binormal model, both algorithms have as the standard normal distribution function, while Algorithm 1 has as the normal with mean 02 and variance one, and Algorithm 2 has as the normal with mean 02 and variance one. Algorithm 1 has a relatively large separation (Cohen’s D = 0.8) between the score distributions of the + and classes; Algorithm 2 has a slightly diminished performance (Cohen’s D = 0.6). To capture the fact that both algorithms are scoring the same compounds, we simulated the positive scores from a bivariate normal with marginals as described above and correlation parameter 9. The negative scores were simulated from a separate bivariate normal with the described marginals and the same correlation parameter. In the bibeta model, both algorithms have as the beta distribution with = 5, while Algorithm 1 has as the beta(5,2) and Algorithm 2 has as the beta(4,2). Sampling was done using the bivariate beta distributions to incorporate correlations similar to the binormal model. Sampling was conducted using the copula R package (Hofert et al. 2020).

There is much greater separation between in the bibeta model, so the true hit enrichment curves are higher than in the binormal model, and we expect greater correlation of scores across early testing fractions. But for both the binormal and bibeta models, parameters were chosen to result in very similar hit enrichment curves for the two competing algorithms, thus creating a challenging task for hypothesis testing of differences between hit enrichment curves.

Studies were conducted using 10,000 Monte Carlo replicates. We estimated the type I error rate for the hypothesis test methods assuming that both ranking algorithms had either the score distributions of Algorithm 1 or Algorithm 2.

3.5.3 Results

Hypothesis Testing

Pooling has no impact on either the power or protection from type I errors for EmProc, IndJZ, or CorrBinom (results not shown for brevity), so we limit further discussion to the versions of these tests constructed without pooling. The impact of the Bonett-Price plus adjustment on power is mixed (results not shown). The plus adjustment had no noticeable impact on power under the bibeta model. But under the binormal model, the plus adjustment caused a noticable decrease in power for a medium range of tests performed (approximately 30 to 500). For this reason we limit further discussion to the versions of these tests constructed without plus adjustment. As a reminder, McNemar is based on pooling and no plus adjustment.

Figure 2 shows estimated power curves for correlation of 0.1 (A) and 0.9 (B) under the bibeta model. We focus on results for number of tests (nr) ranging from two tests to testing ten percent of the total sample size; in practice, screening campaigns will interrogate only a tiny fraction of the available virtual screening library (Zhu et al. 2013). The CorrBinom and McNemar approaches are noticeably suboptimal. In the presence of weak correlation, EmProc and IndJZ have comparable performance. EmProc dominates in the presence of strong correlation between scores. Indeed, EmProc is the only approach that is designed to address both the correlation between scores of competing algorithms, and the correlation that is induced within a particular algorithm as a result of having to estimate thresholds. The binormal model (results not shown) yielded similar findings. Type I error rates (results not shown) are well controlled to their nominal values of 0.05.

For hypothesis testing, we recommend EmProc without pooling and without plus adjustment, because EmProc has the greatest power compared to IndJZ, CorrBinom, and McNemar, while maintaining control of type I error rates. If one chooses to use either IndJZ or CorrBinom, the unpooled and non-plus-adjusted versions should be used.

Figure 2 shows estimated coverage probabilities (C) and average widths (D) for confi-dence intervals EmProc, IndJZ, and CorrBinom, based on the binormal model with correlation of 0.9. The McNemar interval is equivalent to the CorrBinom interval with plus adjustment, so while the figure does not explicitly include the label of McNemar, it is included.

The most obvious finding is that the Bonett-Price plus adjustment dramatically improves coverage probabilities when the number of tests is small; this is because the plus adjustment results in wider intervals when the number of tests is small. As the number of tests increase, the plus and no-plus versions converge, and they approach nominal coverage. By not accounting for correlation across competing algorithms, IndJZ standard errors are unnecessarily large, resulting in wide intervals that provide conservative coverage. The plus-adjusted versions of EmProc and CorrBinom (and hence also McNemar) provide conservative coverage when the number of tests is small, but coverage approaches the nominal level as number of tests increase.

For confidence intervals, we recommend the plus-adjusted version of EmProc, because it is best able to balance achieving nominal coverage rates while minimizing the width of confidence intervals. And if other procedures are used, the plus adjustment should be used as well.

Figure 2: Comparison of EmProc, CorrBinom, IndJZ, and McNemar in terms of hypothesis testing and pointwise confidence intervals to compare hit enrichment curves for competing algorithms. (A-B): Estimated power of the hypothesis test to detect differences between two competing algorithms, where each algorithm follows a bibeta model and scores are correlated with 9 (B). (C-D): Estimated coverage probability (C) and average width (D) of pointwise confidence intervals for the difference in hit enrichment curves for two competing algorithms, where each algorithm follows a binormal model and scores are correlated with 9. Simulations were conducted with 10,000 Monte Carlo replicates. Shading show the Monte Carlo estimate 1.96 times the Monte Carlo standard error.

4.1 Bands for a Single Algorithm

4.1.1 Methods

An ideal scenario would be to accompany hit enrichment curves, such as those shown in Figure 1, with confidence regions. Non-overlapping regions would provide an alternative justification for claiming significant differences between competing algorithms. Let the vector of recall values from a single algorithm at the vector ordered testing fractions . We seek a 100(1) percent confidence region for . While the pointwise confidence interval approach of equation (3) could be modified using a Bonferroni adjustment, such corrections are known to be conservative when k is large, leading to unnecessarily wide intervals.

In their technical report, Jiang and Zhao (2014) suggested an alternate confidence band estimation procedure, and gave brief comments on simulation results, but some details were omitted. We complete these details to state the following result. Under the previously mentioned Conditions 1 and 2, as is the vector of recall estimators as previously defined, and . Moreover,

V

Derivation details are omitted because they are similar to the steps in the appendix; see

Ash (2020) for further details. To provide a working distribution for , an estimator of V

is obtained by replacing population parameters with consistent estimators. This working distribution is the basis of our approximate confidence regions.

Our most straight-forward approach is to use a Wald 100(1) percent confidence ellipsoid, defined aspercentile

of the chi-squared distribution with k degrees of freedom. But the Wald confidence ellipsoid does not produce regions that are of the rectanguloid formWe have chosen to use confidence regions with a rectangular structure (i.e., a confidence band) and not ellipsoids because this allows confidence regions of high dimensions to be easily visualized.

Clearly, Bonferroni regions are rectanguloid, with . We mention two additional rectanguloid regions, following the naming conventions in Montiel Olea and Plagborg-Møller (2019): the -projection and sup-t bands. -projection bands are obtained by identifying the smallest rectanguloid that contains the Wald ellipsoid, and results in . Upon further inspection, it becomes clear that -projection bands are always at least as wide as Bonferroni bands, so they are not considered further. On the other hand, sup-t bands are the smallest rectanguloid that maintains the simultaneous coverage probability of 1 , and are expected to have good performance. Their critical value q

must be obtained using Monte Carlo sampling. Briefly,

Monte Carlo sampling is used to estimate )100 percentile for the distribution

4.1.2 Simulation Results

A study compared coverage probabilities achieved by confidence bands constructed using sup-t and Bonferroni approaches. Results are shown for both the standard (formulas shown in Section 4.1.1) and plus-adjusted versions of sup-t and Bonferroni bands. The familiar trick of “add two successes and add two failures” (Agresti and Coull 1998) before estimating proportions is what we refer to as the plus adjustment for bands corresponding to a single hit enrichment curve, not the Bonett-Price plus adjustment for comparing two algorithms that was used earlier. We computed bands for using a grid of number of compounds tested between two and 15,000. We considered 25 points on the grid, defined as: 28; 105, 300, 1500, 15000.

Figure 3 shows estimated coverage probabilities (A) and average widths (B) of confi-dence bands created based on five distributional cases. The distributional cases represent varying degrees of separation between + and classes, and are chosen to mimic real-world scenarios. Case 1 is a binormal model, with equal unit variance and means zero and 1.4. Case 2 is another binormal model, with equal unit variance and means zero and 0.5; Case 2 offers much less separation than Case 1, so Case 2 results in lower values of recall. Case 3 is a bibeta model, with Beta(2,5) and Beta(5,2) distributions. Case 4 is another, more separated, bibeta model, with Beta(1,20) and Beta(20,1) distributions. Case 5 is made up of distributions of limited extent, namely uniform on (0,0.75) and uniform on (0.25,1).

The plus-adjusted Bonferroni bands have the highest coverage, but they are also the widest. The plus-adjusted sup-t bands are not as wide, yet have excellent coverage. As such, for confidence bands applied to a single hit enrichment curve, we recommend the plus-adjusted sup-t bands.

4.2 Bands for the Difference Between Competing Algorithms

4.2.1 Methods

While the pointwise confidence intervals of Section 3 offer effective comparisons of competing algorithms at a few selected testing fractions, it may be more desirable to perform comparisons across a large range of testing fractions. This may be accomplished by converting the pointwise confidence intervals into confidence bands, in much the same way that confidence bands were obtained in Section 4.1.

Letting denote the vector of recall values from Algorithm 2), the method is

Figure 3: Comparison of sup-t and Bonferroni confidence bands. (A-B): Estimated coverage probability (A) and average width (B) of confidence bands for hit enrichment curves for a single algorithm, where the algorithm is generated from five different cases. (C-D): Estimated coverage probability (C) and average width (D) of confidence bands for the difference between two hit enrichment curves generated under four scenarios. Simulations were conducted with 10,000 Monte Carlo replicates. Error bars show the Monte Carlo estimate 1.96 times the Monte Carlo standard error.

based on the asymptotic result

and for

The first two components of equation (7) are obtained using equation (6), and the latter

two components are obtained using

where θ

conditional and unconditional probabilities that both algorithms test a ligand because it is

highly ranked by both algorithms, albeit at different testing fractions (8) does not impose any restrictions between testing fractions

As described in Section 4.1, matrix V is estimated and used to construct sup-t and Bonferroni bands.

4.2.2 Simulation Results

Similar to Section 3, a study was conducted to compare coverage probabilities and average widths of confidence bands constructed using sup-t and Bonferroni approaches under four settings of two competing algorithms: binormal or bibeta, and 1 or 0.9. Bands were computed for using a grid of size 25, with number of tested compounds being: 28; 105, 300, 1500, 15000.

Figure 3 shows estimated coverage probabilities (C) and average widths (D). Results are very similar to those observed for confidence bands for a single algorithm, namely that the plus-adjusted sup-t bands provide the best balance between coverage probabilities and average width.

For confidence bands applied to the difference between two hit enrichment curves, we recommend the plus-adjusted sup-t bands for achieving nominal coverage rates and minimizing width. The covariance used in constructing these bands arise from the EmProc approach.

For testing fractions 0.001, 0.01, and 0.1, Table 2 provides details for all pairwise comparisons between scoring methods Surflex-dock (the best docking method), ICM (the worst docking method), and the maximum z-score (the best consensus method). For each of the three pairs, we provide the estimated difference between hit enrichment curves. Standard errors and resulting p-values are provided for the EmProc approach to conducting inference, and also for the remaining approaches McNemar, IndJZ, and CorrBinom.

Acknowledging the multiple-testing scenario required to compare two scoring methods, Table 2 also provides multiplicity-adjusted p-values. Given choice of a particular approach to inference (for example, EmProc), there are nine tests based on three pairs of scoring methods and three testing fractions. A simple Benjamini-Hochberg (Benjamini and Hochberg 1995) step-up procedure is used to control the false discovery rate. At testing fraction 0.1 (321 tests), the difference between Surflex-dock and the consensus method is not significant, but ICM is significantly worse than both Surflex-dock and the consensus method. These conclusions are clearly supported by all inferential approaches. At testing fraction 0.01 (32 tests), no significant differences are observed, but the EmProc procedure is seen to produce the smallest standard errors. With our relatively small dataset, testing fraction 0.001 results in only three tests, and there is too much uncertainty to make a reliable conclusion.

Ignoring the need for multiplicity adjustments, Figure 4 provides information similar to Table 2, except for many more testing fractions. The upper grid shows (unadjusted) p-values from testing for equality of hit enrichment curves, while the lower grid shows (unadjusted) pointwise 95 percent confidence intervals of differences between hit enrichment curves. EmProc is the only method that is consistently among the best performers. The Pearson correlation is largest (but still only moderate at 0.624) between scores from Surflex-dock and the consensus method, so the poor performance of IndJZ for this comparison is not surprising. On the other hand, the performance of IndJZ understandably improves when comparing Surflex-dock and ICM, with Pearson correlation of only 0.156 between scores.

The diagonal grid of Figure 4 shows estimated score densities within activity classes for each scoring method. The consensus method has the biggest separation of activity-class densities, with Kullback-Leibler divergence (Kullback and Leibler 1951) of 2.50 compared to the substantially smaller values of 0.00162 for surflex-dock and 2for ICM, and this is consistent with the consensus scoring method seeming to outperform the others.

Taking an alternative approach, Figure 5 shows plus adjusted sup-t confidence bands for the three scoring methods. These bands account for correlation between recall values at distinct testing fractions for a single curve, and offer simultaneous coverage across the curve. While they do not account for correlations between curves, they are helpful visualizations of uncertainty that go beyond simply graphing the curves alone.

A more direct approach to pairwise comparisons between curves, while adjusting for the many comparisons that occur along the curve, is shown in Figure 6. These bands for the differences between hit enrichment curves offer simultaneous coverage across the curve. They account for correlation between recall values at distinct testing fractions for a single curve, and for correlation between estimated recall from different scoring methods.

Figure 4: Comparisons of scoring methods surflex-dock, ICM, and their consensus, across 15 testing fractions using four testing procedures. The diagonal grid shows estimated score densities within activity classes for each scoring method; also shown are Kullback-Leibler divergences between estimated densities for the + and classes. The upper grid shows p-values from testing for equality of hit enrichment curves from a pair of scoring methods, with colors corresponding to different testing procedures. The lower grid shows pointwise 95 percent confidence intervals that accompany results from the upper grid. Comparisons

have not been adjusted for multiple testing as was the case in Table 2.26

Figure 5: Simultaneous 95 percent plus-adjusted sup-t confidence bands for the three scor- ing methods surflex-dock, ICM, and their consensus. Even while spreading interest across the entire range of testing fractions, significant differences are still detected between the consensus and ICM for some intermediate testing fractions.

Early enrichment is broadly recognized as a primary goal of virtual screening (Truchon and Bayly 2007). There is, however, considerable debate about how to assess achievement of this goal (Robinson et al. 2020). The receiver operating characteristic (ROC) curve, and its summary metric of area under the curve, lack sensitivity to both the early recognition goal and the rarity of active compounds that often exists in screening studies (Truchon and Bayly 2007; Saito and Rehmsmeier 2015). Partial area under the ROC curve makes a step towards addressing the early recognition goal but not the rarity of active compounds.

Figure 6: Simultaneous 95 percent plus-adjusted sup-t confidence bands for pairwise dif- ferences between the three scoring methods surflex-dock, ICM, and their consensus. A: consensus versus surflex-dock, B: consensus versus ICM, C: surflex-dock versus ICM. The bands are unsurprisingly wider than the pointwise EmProc intervals shown in the lower panel of Figure 4, but they are not very much wider. The ICM scoring method is signifi-cantly less effective than both consensus and suflex-dock for testing fractions between 0.02 and 0.5.

The robust initial enhancement (RIE, Sheridan et al. (2001)) and its normalized version the Boltzmann enhanced discrimination of ROC (BEDROC, Truchon and Bayly (2007)) directly address early recognition by incorporating an exponential weight that decreases for active compounds that are discovered later in testing. The versions of RIE and BEDROC that are proposed by Truchon and Bayly (2007) are computationally efficient. They are, however, global measures that are not tied to a specific testing fraction, and they depend on a parameter that must be specified by the user. A larger provides greater weight to active ligands found early, but weights are applied to all active ligands, even those found very late. Moreover, statistical inference is not straightforward.

For the PPARg study using the default value of = 20, the Empereur-Mot et al. (2016) web application reports BEDROC values as follows, where larger is better: 0.743 for the consensus; 0.687 for Surflex-dock; and 0.447 for ICM. Without uncertainty measures associated with these scores, it is difficult to conclude separation of the scoring methods. Furthermore, the results presented in Table 2 suggest conclusions that are more nuanced, with no significant differences at testing fraction 0.01 but some big differences at testing fraction 0.1. The improvement in BEDROC for the consensus method is mostly driven by the improvement near testing fraction 0.1, but this is a much larger fraction that what is typically used in screening evaluations. This illustrates how BEDROC can provide a misleading sense of improvement in early enrichment. It averages performance over all testing fractions and weights early fractions more heavily, but it is often difficult to determine from the tuning parameter alone to what extent the early fractions are weighted.

The hit enrichment curve offers multiple benefits for assessing early enrichment. First, it directly addresses the measure of interest, namely enrichment. For a given dataset, the closer the estimated hit enrichment curve is to the ideal hit enrichment curve, the more desirable is the associated algorithm. Second, the user is able to directly enforce their definition of “early” by specifying testing fractions of interest, without needing to rely on the indirect methods offered through measures such as BEDROC or partial area under the ROC curve. And third, localized assessment is possible, rather than whole-curve assessment. These are the reasons we have chosen to focus this work on hit enrichment curves.

Enrichment factors, and the associated enrichment factor curve, are also very popular for assessing early enrichment. As it turns out, all results in this paper can be trivially modified to obtain inference for the enrichment factor curve; see Ash (2020) for further details. While the hit enrichment curve plots , the enrichment factor curve plots so standard error expressions are easily modified by division by the testing fraction r. Enrichment factors (sometimes denoted ) simply focus on the enrichment factor curve at a specific testing fraction.

In this article, we provide a template for rigorously comparing competing algorithms while accounting for uncertainty, two types of correlation, and the multiple testing issue. If interest is restricted to comparing hit enrichment curves for competing algorithms at a few pre-selected testing fractions, the hypothesis testing and confidence interval procedures of Section 3 offer effective strategies. On the other hand, while it is best practice to decide testing fractions a priori, we acknowledge this is not always done or even possible, so we also provide confidence bands to compare entire hit enrichment curves. Additionally, bands allow agmented graphical presentation of the entire hit enrichment curves for competing algorithms.

The EmProc procedure is newly proposed here and is expected to perform as well or better than the other procedures considered (CorrBinom, McNemar, and IndJZ) in the presence of correlation between competing algorithms and/or correlation across different testing fractions within a single algorithm. The other procedures considered address only a single type of correlation and hence are not recommended for general-purpose use.

Finally, these inferential procedures are conveniently made available in R package

chemmodlab on CRAN and GitHub.

Appendix: Proof of Theorem 1

Agresti, A. (2007), An Introduction to Categorical Data Analysis, Wiley Series in Probability and Statistics, John Wiley & Sons, Inc., Hoboken, NJ.

Agresti, A. and Coull, B. A. (1998), ‘Approximate is better than ‘exact’ for interval estimation of binomial proportions’, The American Statistician 52(2), 119–126.

Ash, J. R. (2020), Methods Development for Quantitative Structure-Activity Relationships, North Carolina State University.

Benjamini, Y. and Hochberg, Y. (1995), ‘Controlling the false discovery rate: A practical and pow- Journal of the Royal Statistical Society. Series B (Methodological) 57, 289–300.

Bonett, D. G. and Price, R. M. (2012), ‘Adjusted Wald confidence interval for a difference of Journal of Educational and Behavioral Statistics 37(4), 479–488.

Empereur-Mot, C., Zagury, J.-F. and Montes, M. (2016), ‘Screening Explorer–An Interactive Journal of Chemical Information and Modeling 56(12), 2281–2286. (Web application at http://stats.drugdesign.fr).

Fagerland, M. W., Lydersen, S. and Laake, P. (2013), ‘The McNemar test for binary matched- pairs data: mid-p and asymptotic are better than exact conditional’, BMC Medical Research Methodology 13(1), 91.

Friesner, R. A., Banks, J. L., Murphy, R. B., Halgren, T. A., Klicic, J. J., Mainz, D. T., Repasky, M. P., Knoll, E. H., Shelley, M., Perry, J. K., Shaw, D. E., Francis, P. and Shenkin, P. S. (2004), ‘Glide: A new approach for rapid, accurate docking and scoring. 1. method and assessment of docking accuracy’, Journal of Medicinal Chemistry 47(7), 1739–1749.

Geppert, H., Vogt, M. and Bajorath, J. (2010), ‘Current trends in ligand-based virtual screen- ing: Molecular representations, data mining methods, new application areas, and performance Journal of Chemical Information and Modeling 50(2), 205–216.

Good, A. C. and Oprea, T. I. (2008), ‘Optimization of camd techniques 3. virtual screening en- richment studies: a help or hindrance in tool selection?’, Journal of Computer-Aided Molecular Design 22(3-4), 169–178.

Hofert, M., Kojadinovic, I., Maechler, M. and Yan, J. (2020), copula: Multivariate Dependence with Copulas. R package version 1.0-1. URL: https://CRAN.R-project.org/package=copula

Irwin, J. J. (2008), ‘Community benchmarks for virtual screening’, Journal of Computer-Aided Molecular Design 22(3-4), 193–199.

Jain, A. N. and Nicholls, A. (2008), ‘Recommendations for evaluation of computational methods’, Journal of Computer-Aided Molecular Design 22(3-4), 133–139.

Jiang, W. and Zhao, Y. (2014), Some technical details on confidence intervals for lift measures in data mining, Technical report.

Jiang, W. and Zhao, Y. (2015), ‘On asymptotic distributions and confidence intervals for lift Journal of the American Statistical Association 110(512), 1717– 1725.

Kuhn, M. (2008), ‘Building predictive models in R using the caret package’, Journal of Statistical Software 28(5), 1–26.

Kullback, S. and Leibler, R. A. (1951), ‘On information and sufficiency’, The Annals of Mathematical Statistics 22(1), 79–86.

Li, Q. and Racine, J. S. (2007), Nonparametric econometrics: theory and practice, Princeton University Press.

Montiel Olea, J. L. and Plagborg-Møller, M. (2019), ‘Simultaneous confidence bands: Theory, implementation, and an application to svars’, Journal of Applied Econometrics 34(1), 1–17.

Mysinger, M. M., Carchia, M., Irwin, J. J. and Shoichet, B. K. (2012), ‘Directory of useful decoys, enhanced (dud-e): Better ligands and decoys for better benchmarking’, Journal of Medicinal Chemistry 55(14), 6582–6594.

NCBI (2021). https://www.ncbi.nlm.nih.gov/gene/5468.

NCBI Resource Coordinators. (2016), ‘Database resources of the national center for biotechnology information’, Nucleic Acids Research 44(D1), D7–D19.

Newcombe, R. G. (1998), ‘Improved confidence intervals for the difference between binomial proportions based on paired data’, Statistics in Medicine 17, 2635–2650.

Nicholls, A. (2014), ‘Confidence limits, error bars and method comparison in molecular modeling. part 1: the calculation of confidence intervals’, Journal of computer-aided molecular design 28(9), 887–918.

Robinson, M. C., Glen, R. C. et al. (2020), ‘Validating the validation: reanalyzing a large-scale comparison of deep learning and machine learning models for bioactivity prediction’, Journal of computer-aided molecular design pp. 1–14.

Rodriguez de Gil, P., Pham, J. R. T., Nguyen, D., Kromrey, J. D. and Kim, E. S. (2013), SAS Macros CORR-P and TANGO: Interval Estimation for the Difference Between Correlated Proportions in Dependent Samples, in ‘Proceedings of the SouthEast SAS Users Group 2013’.

Rohrer, S. G. and Baumann, K. (2009), ‘Maximum unbiased validation (muv) data sets for virtual screening based on pubchem bioactivity data’, Journal of Chemical Information and Modeling 49(2), 169–184.

Rosset, S., Neumann, E., Eick, U., Vatnik, N. and Idan, I. (2001), Evaluation of Prediction Models for Marketing Campaigns, ACM Press, New York, New York, pp. 456–461.

Saito, T. and Rehmsmeier, M. (2015), ‘The precision-recall plot is more informative than the roc plot when evaluating binary classifiers on imbalanced datasets’, PloS one 10(3), e0118432.

SAS Institute Inc. (2020), ‘SAS Enterprise Miner 15.1’.

Sheridan, R. P., Singh, S. B., Fluder, E. M. and Kearsley, S. K. (2001), ‘Protocols for bridging the peptide to nonpeptide gap in topological similarity searches’, Journal of Chemical Information and Computer Sciences 41(5), 1395–1406.

Sterling, T. and Irwin, J. J. (2015), ‘Zinc 15 – ligand discovery for everyone’, Journal of Chemical Information and Modeling 55(11), 2324–2337.

Truchon, J.-F. and Bayly, C. I. (2007), ‘Evaluating virtual screening methods: Good and bad metrics for the “early recognition” problem’, Journal of Chemical Information and Modeling 47(2), 488–508.

Xia, J., Tilahun, E. L., Reid, T.-E., Zhang, L. and Wang, X. S. (2015), ‘Benchmarking methods and data sets for ligand enrichment assessment in virtual screening’, Methods 71, 146–157.

Zhu, T., Cao, S., Su, P.-C., Patel, R., Shah, D., Chokshi, H. B., Szukala, R., Johnson, M. E. and Hevener, K. E. (2013), ‘Hit identification and optimization in virtual screening: Practical recommendations based on a critical literature analysis’, Journal of Medicinal Chemistry 56(17), 6560–6572.