The methods developed in the literature for correcting verification bias are usually based on the following assumptions that limit their usefulness in applications. First, most currently available methods assume that the diagnostic test's scale is ordinal (Begg and Greenes, 1983; Gray et al., 1984; Baker, 1995; Toledano and Gatsonis 1996; Zhou, 1996, 1998a,b,c; Rodenberg and Zhou, 2000). Yet many important emerging tests are measured on a continuous scale. Second, with the exception of few articles which only deal with dichotomous diagnostic tests (e.g. Baker, 1995; Zhou, 1993, 1994; Kosinski and Barnhart, 2003), currently available methods assume that the decision to send a patient to verification is conditionally independent of the true disease status of the patient given the test results and possibly other observed covariates (see Alonzo and Pepe, 2005, for continuous scales), or equivalently, that the missing disease status is missing at random (MAR) (Rubin,1976). However, usually the doctor's decision to send a patient to verification will be based on his/her detailed information on the patient's health, which can hardly ever be accurately summarized by the test result and other measured covariates. The aforementioned methods can yield biased estimates in this case because the assumption of MAR data no longer holds since non-response and true disease status are dependent even after adjusting for measured variables. In such a case, the missing process on the disease status is known as non-ignorable.

As an example, consider a study run by the Nuclear Imaging Group at Cedars Sinai Medical Center discussed by Rotnitzky et al. (2006). A diagnostic test for coronary artery disease, electron beam computed tomography (EBCT) (Braun et al., 1996), was performed on 5130 males, of which only 379 of these had disease status checked with the more expensive Dual Isotope Myocardial Perfusion Single Photon Emission Computed Tomography (SPECT) test. Of the verified subjects only 28 were found to be diseased. The EBCT distribution of non-verified subjects is markedly skewed to the right (when compared to the verified subjects) indicating that subjects with lower values of the marker are less likely than others to be verified. It is indeed very likely that doctors based their decision to send a patient to have a SPECT test not just on the values of the EBCT marker, but on other clinical variables not available for data analysis.

Suppose that a subject i were to be classified as diseased if T_i = I (Y_i > c) is equal to 1 and as non-diseased if T_i = 0, for a given constant c. The sensitivity, Se(c) = Pr(Y > c|X = 1), and specificity, Sp(c) = Pr(Y < c|X = 0), would then be the probability of diagnosing a diseased and a healthy person correctly, respectively. For notational convenience, in what follows the c will on occasions be dropped. The theoretical ROC curve is the plot of sensitivity against 1-specificity for all possible values of c. Consequently estimation of the ROC curve follows immediately from estimation of the θ(c) = (Se(c), Sp(c)) pairs. Supposing that all patients are verified (i.e. R_i = 1), the empirical estimator is consistent for θ(c). If not all subjects undergo verification, i.e. if R_i is not 1 for all i, but the process of selection to verification is completely random, i.e. if R_i and X_i are independent, then calculated from verified subjects only remains consistent for θ(c). However, if R_i and X_i are correlated, then this naive estimator will lead to inconsistent estimates of θ.

The key issue that complicates the analysis is that both Se(c) and Sp(c) are unidentified from the observed data O_i, i = 1, ..., n. That is, there exist many curves Se(c) and Sp(c), c R, compatible with the distribution of O_i, so that even in the ideal situation in which the distribution of the observed data were perfectly known, there would remain uncertainty about the true ROC curve. To identify the ROC curve and associated summary measures one needs to make untestable assumptions, i.e. assumptions that cannot be rejected by any statistical test no matter how large n is. Rotnitzky et al. (2006) argued that the area under the ROC curve is identified under the following untestable assumption, where q (Y, V) is an arbitrary specified (i.e. known) function and h (Y, V) is an arbitrary but unknown function. Indeed, this assumption suffices also to identify θ(c) for any fixed c.

Application of Bayes rule shows that (2.2) holds for some h (Y, V) if and only if where c (Y, V) = E {exp (q (Y, V) X) |R = 1, Y, V} (Rotnitzky et al., 2006). The function q (Y, V) measures the degree of the residual association between R and X within levels of Y and V. For example, the choice q (Y, V) = −1 indicates that the residual association is constant within levels of (Y, V) and such that among subjects with identical values of (Y, V), the odds of being sent to verification is exp (1) = 2.71 times larger for diseased subjects than for non-diseased subjects. The choice q (Y, V) = 0 for all V and Y corresponds to the assumption that the missing disease status is missing at random. In general, when q (Y, V) ≠ 0, the selection process is said to be non-ignorable.

We refer to the model that assumes just (2.2) with specified q (Y, V) and h (Y, V) unknown (or equivalently equation (2.3) with specified q (Y, V)) as model . Arguing as in Scharfstein et al. (1999), it can be shown that imposing a given choice of q does not restrict the observed data distribution. That is, model is a non-parametric model for the distribution of O. Since all values of the selection bias parameter determine the same model for the observed data distribution, the selection bias parameter is untestable, because any q will perfectly fit the observed data. Following Robins et al. (2000) we therefore recommend conducting inference about θ(c) by repeating the estimation under different plausible choices for the selection bias function q as a form of sensitivity analysis to different degrees of selection bias. This raises the question of how to choose the selection bias functions in practice. We suggest that one chooses a collection of simple selection bias functions indexed by one or two parameters β that are to be varied in a sensitivity analysis, with values of the parameters equal to zero corresponding to the function q = 0. For example, the choice q (Y, V) = β is tantamount to assuming that the odds of verification of diseased vs non-diseased subjects is constant across all values of the marker and the auxiliary V and equal to exp (–β). This is the choice used in our simulation study in Section 4. One could consider the choice q (Y, V) = β₀ + β₁Y if the investigator believed that the odds of verification of diseased vs non-diseased subjects is modified by the values of Y, this modification being monotone in Y.

The identity (2.3) implies that Consequently, writing Se(c) as E {E (X|Y, V, R) T} / E {E (X|Y, V, R)} and Sp(c) as E {[1 – E (X|Y, V, R)] (1 – T)} / E {[1 – E (X|Y, V, R)]} we obtain that under model and Likewise, letting π (Y, V) denote the quantity (1 + exp {h (Y, V) + q (Y, V) X})⁻¹, the identity (2.2) implies that P (R = 1|X, Y, V) = π (Y, V). Consequently, under modelg , Se(c) and Sp(c) can be expressed also as

Model implicitly assumes that all subjects have a positive probability of being selected for verification regardless of the values of their test results and covariates. See Rotnitzky et al. (2006) for a discussion of the possibility of relaxing this condition. For technical reasons related to the finiteness of the variance of the estimators that we will propose later, we further need to assume that

Unfortunately, though sufficient for identification, model is insufficient for estimation of θ(c) due to the curse of dimensionality. Specifically, expressions (2.4), (2.5) and (2.6) imply that if one were to estimate θ(c) under model , one would need to non-parametrically estimate either the function h (Y, V) or the conditional expectation E (X|Y, V, R = 1) using smoothing techniques, a practically unfeasible task when the dimension of V is large. Nevertheless, the expressions (2.4), (2.5) and (2.6) suggest two alternative dimension reducing strategies for estimating θ(c).

The first option is to assume that the unknown function h (Y, V) follows a parametric model where h (Y, V; γ) is smooth in γ and γ⁰ is an unknown p_h × 1 finite dimensional parameter. Assumption (2.8) and assumption (2.2) with specified q (Y, V) define a semiparametric model for the observed data which for ease of reference we refer to as . Note that this model imposes just a, possibly non-linear, logistic regression model on the selection probabilities with offset q (Y, V) X. Expressions (2.6) imply that one can obtain a consistent and asymptotically normal estimator estimator of θ(c), throughout referred to as inverse probability weighted (IPW) estimator, under model B (q) by replacing X_i in (2.1) with X_i, where π (Y, V; γ) = (1 + exp {h (Y, V; γ) + q (Y, V) X})⁻¹ and is a consistent estimator of γ⁰. Unfortunately, even though under model , R follows a logistic regression on Y, X and V, we cannot find from a standard logistic regression fit with offset because the offset q (Y, V) X is not observed when R = 0. Instead, following Rotnitzky et al. (2006) we can compute as the solution to where and u (Y, V) is an arbitrary, user specified, column vector function of the same dimension as γ. The choice of u impacts the efficiency with which we estimate γ. Using the theoretical results in Rotnitzky and Robins (1997), it can be shown that if u (Y, V) is equal to then the solution of (2.9) has asymptotic variance equal to the semiparametric variance bound for estimators of γ under model . Note, in particular, that for the choice q (Y, V) = 0, B (γ; u) reduces to which, at γ⁰, is equal to the score, i.e. the derivative of the log-likelihood for γ, in the logistic regression model logit P (R = 1|Y, V) = –h (Y, V; γ).

A second option is to assume that E (X|Y, V, R = 1) = Pr (X = 1|R = 1, Y, V) follows a parametric model where m (Y, V; μ) is a known function, smooth in μ, and μ⁰ is an unknown p_m × 1 parameter vector. Assumption (2.11) is a model for the disease probabilities among the verified subjects. The model defined by (2.3) with specified q (Y, V) and (2.11) is a semiparametric model for the law of the observed data which we will denote with . The Semi-Parametric Maximum Likelihood (SPML) estimator of θ(c) under model is obtained by replacing each X_i in (2.1) with , where with Pr (X = 1|R = 1, Y, V; μ) = exp {m (Y, V; μ)} / [1 + exp {m (Y, V; μ)}], and solves the score equations with .

Neither option is entirely satisfactory because option 1 results in inconsistent estimation if model is incorrect and option 2 results in inconsistent estimation if model is misspecified. There is yet a third option which provides an estimator consistent for θ(c) under either model or model but not necessarily both. Specifically, define where Let be the estimator of θ obtained by replacing each X_i in (2.1) with . That is,

The estimator is said to be double-robust because, as stated in Theorem 1 below, it is consistent and asymptotically normal for θ so long as model , condition (2.7) and one of the models (2.11) or (2.8) holds, but not necessarily both. The key to the consistency of lies in the fact that if model holds then, if model (2.11) additionally holds, or if model (2.8) additionally holds,

In contrast to the first two estimators IPW and SPML, the double-robust estimator gives the analyst two chances, instead of only one, to get nearly correct inference. Of course, there can be an efficiency price to using the DR estimator. If the disease regression model (2.11) is correct both the DR estimator and the SPML estimator will be consistent but the DR estimator will generally have larger (asymptotic) variance. Clearly, the efficiency gains over are obtained at the cost of potential severe bias.

An interesting remark is that, invoking the theoretical results of Rotnitzky and Robins, 1997, it can be shown that if we use u_opt (Y, V) to compute then the that uses this is locally semiparametric efficient under model at the local model . This means that, if both models and are correct, has asymptotic variance that attains the semiparametric variance bound for estimators bof θ(c) under model . This is unfeasible because to compute u_opt (Y, V) we need to know γ⁰ and, when qB(Y, V) ≠ 0, we also need the conditional probabilities P (X = 1|R = 1, Y, V). However, a feasible double-robust estimator that retains the local efficiency properties of the unfeasible can be obtained with the following two-stage algorithm. At the first stage we compute and a preliminary estimator that solves (2.9) using any u (Y, V). We then calculate the function (Y, V) defined like u_opt (Y, V) in (2.10) but with replacing γ⁰ and with the bexpectations in the numerator and denominator computed under P (X = 1|R = 1, Y, V; ). At the second stage, we compute solving (2.9) with u replaced by and finally compute using this . This remark raises the interesting point of how much efficiency is gained by the two-stage locally efficient compared to the single stage computed using and the preliminary . Our simulation study in Section 4 explores this issue using the natural, easy to compute, choice of u (Y, V) = h (Y, V; γ) /γ for the preliminary estimator . For this choice, the theoretical results of Rotnitzky and Robins, 1997, imply that the single stage and the two-stage estimators are equally asymptotically efficient if q (Y, V) = 0 but the two-stage estimator has strictly smaller asymptotic variance if q (Y, V) ≠ 0. For ease of reference, we call the single stage estimator that uses with u (Y, V) = h (Y, V; γ) /γ the DR1 estimator and we call the two-stage estimator using this as the preliminary estimator of γ, the DR2 estimator.

Alonzo and Pepe (AP), 2005, derived a locally-efficient double-robust estimator that, except for a minor difference in the computation of , is computed identically to our DR2 estimator in the case q = 0. Note that when q = 0, the function u_opt (Y, V) reduces to . The distinction in the computation of between our procedure and that of AP is that our solves (2.9) with u (Y, V) replaced by while AP's solves (2.9) but with u (Y, V) replaced by the function u (Y, V; γ) = π (Y, V; γ) × h (Y, V; γ) / γ, depending on the unknown γ. It is easy to show that both estimators of γ are asymptotically equivalent and efficient. Consequently, when q = 0, our DR2 estimator is asymptotically equivalent to AP's estimator of θ.

The following Theorem whose proof is sketched in the Appendix establishes the asymptotic properties of the estimator computed using any fixed function u (Y, V) . In particular, it establishes the double-robustness property of . It also provides a consistent variance estimator when model , condition (2.7) and either model (2.11) or model (2.8) holds. Further details can be found in Fluss (2006). To state the theorem define θ₁ = Se (c) and θ₂ = Sp (c), κ = E (X), where γ* and μ* are the probability limits of and . Note that under standard regularity conditions for L-estimators (see for example, van der Vaart (1998)), γ* = γ⁰ when model (2.8) holds and μ* = μ⁰ when model (2.11) holds.

Theorem 1 Suppose that model and (2.7) hold. Under standard regularity conditions for L-estimators, if model (2.8) and (2.7) hold or if model (2.11) holds, , where Ω = Cov(M). Ω can be consistently estimated with where

In our simulations we computed under the following four scenarios for the working models:

It can be shown that DR2 and DR1 are algebraically identical when they are computed under the working models considered in (iv) . They are also identical under the setting (iii) if the estimator assumes β = 0 (see proof in Fluss (2006)). The performances of the estimators of θ(c) were assessed by means of Monte Carlo bias and standard error (MCSE). In addition, we examined the agreement between the Monte Carlo mean of the estimated standard errors (computed with the formulae given in Theorem 1) and the MC standard error of the DR estimators. As formulae for estimates of the standard error of the DR-PAV estimators are not presently available, for these estimators we report results on the bootstrap estimator of standard error. Due to time considerations the bootstrap was only computed for DR2 and n = 200. The symbol ‘−’ stands for unavailable data.

In our simulations, applying the PAV algorithm generally had little effect on both the bias and standard error so here we will restrict our discussion to the simulation results for the non-PAV adjusted estimators. As predicted by theory, the simulation results show that the DR estimators greatly corrected the noticeable bias of the NAIVE estimator (especially when n=2000) in any of the first three scenarios in which either the disease or the verification processes were correctly modeled. Even when both working models were incorrectly specified (scenario (iv)) the DR method exhibited substantially smaller bias than that of the NAIVE estimator. When the true β was −1 and n = 2000, the MAR/Alonzo-Pepe estimator had bias that, even though small (roughly 7% relative bias for Se and 0.5% for Sp in scenarios (i) and (ii)), it was orders of magnitude greater than that of the DR estimator that used the true value of β. When n = 200, and under the same scenarios (i) and (ii) , the biases of the DR estimators were still smaller than those of the MAR estimators but of comparable sizes. In scenario (iii) the bias reduction of the DR estimator compared to the MAR estimator still existed but was less pronounced. Finally, in scenario (iv) in which, according to the theory none of the estimators is consistent, we see that for estimation of Se, the magnitude of the bias is still smaller for the DR compared to that of the MAR estimator but this order reverses for estimation of Sp. The reason why under scenarios (i) – (iii) the MAR estimator is biased downwards is as follows. Under our setting, for subjects with a given value of Y and V, the odds of being verified are 2.7 times greater for those diseased than for those non-diseased. Thus, within each level of Y and V the fraction of non-verified subjects will be larger in the non-diseased group than in the diseased group. However, the procedure that incorrectly assumes MAR will implicitly impute, within levels of V and Y , half of the non-verified to the diseased group and the other half to the non-diseased group. But this will spuriously inflate the left tail of the distribution of Y for disease subjects and deflate the distribution of Y for the non-diseased subjects because most of the imputation will occur for low values of Y since the chances of being verified increase with Y . The effect of this spurious distortion of the tail distribution is that both P (Y ≤ c|X = 0) and P (Y > c|X = 1) will be underestimated. Furthermore, because the diseased subjects are only 30% of the entire sample, the distortion caused on the distribution of Y by incorrectly assigning non-diseased subjects to the non-disease group will be greater on the diseased group than on the non-diseased group. Thus, we would expect (as confirmed in our simulations) the bias in the estimation of Se = P (Y > c|X = 1) to be greater than the bias in the estimation of Sp = P (Y ≤ c|X = 0) .

The MCSE was similar to the Monte Carlo mean of the estimated standard errors, especially for n=2000. For n=200, the estimator of standard error underestimated a bit. As predicted by the theory when β = −1 and both working models are correct, i.e. under scenario (i) , the Monte-Carlo variance of the DR2 estimator was less than that of the DR1 estimator (the variance reduction being, nevertheless, generally small). When both working models are correct and β = 0, the theoretical results establish that both DR1 and DR2 are asymptotically equivalent. This was confirmed in our simulations for n = 2000 since the MCSE of DR1 and DR2 were roughly the same. Yet, for estimation of Se when n = 200, DR2 still exhibited some variance reduction over that of DR1 (roughly 5% reduction).

The comparison of the MCSE of the COMPLETE and DR2 estimators raises another interesting point. The COMPLETE estimator was calculated using all n simulated data points (i.e. as though no X′s were missing). Yet, in our analysis, roughly 70% of the subjects had missing X. An estimator based just on the data points with X observed, would then be expected to have MC variance roughly equal to 1/0.3=3.33 times the MC variance of the COMPLETE estimator. Yet, the MC variance of the DR2 estimator was far smaller than that. For example, in scenario 1 and β = −1, the (MCSE(DR2)/MCSE(COMPLETE))² = 1.57. This result confirms the theoretical efficiency properties of the DR2 procedure. Specifically, the estimator that uses the DR2 estimator with both working models correctly specified, exploits all the information available in the observed data about the parameters Sp and Se. Thus, it not only exploits the information in the units with observed X′s but also in the units with missing X′s.

The naive estimate gives a significantly different ROC curve than the DR and with a smaller area under the curve. We can see there is an evident effect of the value of β on the ROC curve, but all the DR estimated curves show the marker has better diagnostic power than is shown by the naive estimator. We constructed 95% confidence regions for several points on the ROC curve. Three points corresponding to and 0.9 (using β = −1) are presented in Figure 4. Since Fluss (2006) found the logit transformation did not make much of a difference for large sample size we used the simple Bonferroni (BON) and Multivariate Normal (MVN) distribution based CR. The resulting CRs are shown in Figure 4. The MVN CRs are very similar to the Bonferroni ones. The width of the CRs in the sensitivity axes is very wide indicating on a large amount of uncertainty in the sensitivity estimation. This is due to the small number of verified disease cases.