Improving estimation of the volume under the ROC surface when data are missing not at random

TL;DR

This paper introduces a new method for accurately estimating the volume under the ROC surface in diagnostic tests with missing data that are not at random, using a mean score equation approach and instrumental variables.

Contribution

It proposes a novel verification bias correction method for VUS estimation under nonignorable missingness, incorporating instrumental variables and a parametric verification model.

Findings

The new estimators are consistent and asymptotically normal.

Simulation studies show improved finite sample performance.

Application to ovarian cancer data demonstrates practical utility.

Abstract

In this paper, we propose a mean score equation-based approach to estimate the the volume under the receiving operating characteristic (ROC) surface (VUS) of a diagnostic test, under nonignorable (NI) verification bias. The proposed approach involves a parametric regression model for the verification process, which accommodates for possible NI missingness in the disease status of sample subjects, and may use instrumental variables, which help avoid possible identifiability problems. In order to solve the mean score equation derived by the chosen verification model, we preliminarily need to estimate the parameters of a model for the disease process, but its specification is required only for verified subjects under study. Then, by using the estimated verification and disease probabilities, we obtain four verification bias-corrected VUS estimators, which are alternative to those recently proposed by To Duc et al. (2019), based on a full likelihood approach. Consistency and asymptotic normality of the new estimators are established. Simulation experiments are conducted to evaluate their finite sample performances, and an application to a dataset from a research on epithelial ovarian cancer is presented.

Figures

Tables

Table 1. Table 1: Monte Carlo means for the estimators of the parameters in the verification model and for the corresponding four bias–corrected VUS estimators. The true values are: for scenario I, β 0 = 2 , β 1 = 0.5 , β 2 = − 1.2 , λ 1 = − 2 , λ 2 = − 1 formulae-sequence subscript 𝛽 0 2 formulae-sequence subscript 𝛽 1 0.5 formulae-sequence subscript 𝛽 2 1.2 formulae-sequence subscript 𝜆 1 2 subscript 𝜆 2 1 \beta_{0}=2,\beta_{1}=0.5,\beta_{2}=-1.2,\lambda_{1}=-2,\lambda_{2}=-1 and VUS = 0.791 absent 0.791 =0.791 ; for scenario II, β 0 = 1 , β 1 = 1 , λ 1 = − 2 , λ 2 = − 1 formulae-sequence subscript 𝛽 0 1 formulae-sequence subscript 𝛽 1 1 formulae-sequence subscript 𝜆 1 2 subscript 𝜆 2 1 \beta_{0}=1,\beta_{1}=1,\lambda_{1}=-2,\lambda_{2}=-1 and VUS = 0.843 absent 0.843 =0.843 .

		$n=150$		$n=250$		$n=500$
		MSEq	logLike	MSEq	logLike	MSEq	logLike
Scenario I	${\hat{\beta }}_0$	2.081	0.750	1.992	1.189	2.274	1.806
	${\hat{\beta }}_1$	0.595	0.782	0.550	0.658	0.491	0.546
	${\hat{\beta }}_2$	$-$1.281	$-$1.271	$-$1.255	$-$1.236	$-$1.230	$-$1.217
	${\hat{\lambda }}_1$	$-$1.870	$-$0.127	$-$1.919	$-$0.943	$-$2.192	$-$1.763
	${\hat{\lambda }}_2$	$-$0.132	$-$0.122	$-$0.341	$-$0.435	$-$0.902	$-$0.867
	FI	0.775	0.761	0.776	0.771	0.784	0.783
	MSI	0.772	0.757	0.774	0.769	0.783	0.783
	IPW	0.778	0.765	0.776	0.770	0.783	0.782
	PDR	0.773	0.757	0.773	0.765	0.781	0.782
Scenario II	${\hat{\beta }}_0$	3.154	5.071	2.119	5.217	1.553	4.011
	${\hat{\beta }}_1$	1.108	1.870	1.078	1.452	1.039	1.154
	${\hat{\lambda }}_1$	$-$4.131	$-$5.173	$-$3.115	$-$5.721	$-$2.549	$-$4.895
	${\hat{\lambda }}_2$	$-$0.580	$-$4.365	$-$0.899	$-$4.757	$-$1.274	$-$3.913
	FI	0.844	0.799	0.845	0.822	0.843	0.835
	MSI	0.841	0.795	0.843	0.820	0.841	0.833
	IPW	0.841	0.822	0.842	0.829	0.841	0.836
	PDR	0.839	0.801	0.842	0.825	0.841	0.835

	Bias(%)	MCSD	ASD	CP.Asy(%)
FI	$-$0.6	0.070	0.075	95.8
MSI	$-$0.6	0.070	0.075	95.8
IPW	0.1	0.074	0.069	93.2
PDR	$-$0.6	0.070	0.075	95.8

Table 3. Table 2: Relative biases (%), Monte Carlo standard deviations (MCSD), means of estimated standard deviation based on ( 3.2 ) (ASD), means of bootstrap standard deviation based on 250 bootstrap replications (BSD), empirical coverage probabilities (based on asymptotic theory and ( 3.2 ) (CP.Asy) or asymptotic theory and bootstrap standard deviations (CP.Bst)) of 95% confidence intervals for the VUS (%). The working models are correctly specified.

		Scenario II: VUS = 0.843						Scenario III: VUS = 0.457
		Bias	MCSD	ASD	BSD	CP.Asy	CP.Bst	Bias	MCSD	ASD	BSD	CP.Asy	CP.Bst
$n=150$	FI	0.2	0.054	0.055	0.057	88.9	92.2	$-$1.4	0.070	0.144	0.075	96.8	94.5
	MSI	$-$0.3	0.055	0.072	0.059	90.0	92.8	$-$1.4	0.070	0.144	0.075	96.9	94.5
	IPW	$-$0.3	0.059	0.064	0.061	88.9	91.3	$-$0.6	0.083	0.083	0.087	93.8	94.3
	PDR	$-$0.4	0.061	0.088	0.073	89.1	92.8	$-$1.4	0.070	0.144	0.075	96.9	94.5
$n=250$	FI	0.3	0.041	0.038	–	91.0	–	$-$0.1	0.052	0.063	–	96.1	–
	MSI	$-$0.0	0.042	0.040	–	92.3	–	$-$0.1	0.052	0.062	–	96.1	–
	IPW	$-$0.0	0.044	0.045	–	91.6	–	$-$0.2	0.063	0.064	–	94.5	–
	PDR	$-$0.1	0.044	0.049	–	91.2	–	$-$0.1	0.052	0.063	–	96.2	–
$n=500$	FI	$-$0.0	0.028	0.026	–	92.5	–	$-$0.8	0.037	0.040	–	95.6	–
	MSI	$-$0.2	0.028	0.028	–	93.5	–	$-$0.8	0.037	0.039	–	95.6	–
	IPW	$-$0.2	0.030	0.029	–	93.0	–	$-$0.5	0.043	0.045	–	95.3	–
	PDR	$-$0.2	0.030	0.028	–	92.3	–	$-$0.8	0.036	0.039	–	95.7	–
$n=1000$	FI	0.1	0.019	0.019	–	95.5	–	$-$0.2	0.025	0.029	–	95.6	–
	MSI	0.0	0.019	0.020	–	95.8	–	$-$0.2	0.025	0.029	–	95.5	–
	IPW	0.0	0.020	0.021	–	95.4	–	0.0	0.031	0.032	–	95.6	–
	PDR	0.0	0.020	0.020	–	94.8	–	$-$0.2	0.025	0.029	–	95.5	–

Table 4. Table 3: Relative biases (%), Monte Carlo standard deviations (MCSD), means of estimated standard deviation based on ( 3.2 ) (ASD), empirical coverage probabilities (based on asymptotic theory and ( 3.2 ) (CP.Asy) of 95% confidence intervals for the VUS (%).The working verification model (scenarios IV and V) and the working conditional disease model (scenario VI) are misspecified.

		Scenario IV: VUS = 0.843				Scenario V: VUS = 0.74				Scenario VI: VUS = 0.728
		Bias	MCSD	ASD	CP	Bias	MCSD	ASD	CP	Bias	MCSD	ASD	CP
$n=150$	FI	$-$0.8	0.059	0.055	89.9	$-$3.6	0.090	0.063	89.6	$-$5.2	0.081	0.068	90.7
	MSI	$-$1.1	0.059	0.060	91.4	$-$3.7	0.090	0.062	89.6	$-$5.4	0.082	0.073	91.8
	IPW	$-$0.1	0.064	0.063	88.9	$-$2.2	0.070	0.066	95.4	$-$5.5	0.083	0.071	91.1
	PDR	$-$1.1	0.067	0.066	90.3	$-$3.4	0.091	0.079	90.9	$-$3.7	0.097	0.074	91.2
$n=250$	FI	$-$0.7	0.044	0.041	92.8	$-$3.6	0.081	0.046	89.7	$-$5.2	0.067	0.049	87.8
	MSI	$-$1.0	0.044	0.041	93.8	$-$3.7	0.081	0.045	89.2	$-$5.4	0.068	0.055	90.3
	IPW	$-$0.1	0.050	0.050	91.7	$-$2.2	0.059	0.055	94.8	$-$5.5	0.065	0.060	90.0
	PDR	$-$1.1	0.051	0.051	92.3	$-$3.4	0.083	0.054	90.4	$-$3.7	0.068	0.056	90.7
$n=500$	FI	$-$0.9	0.031	0.028	92.7	$-$1.8	0.050	0.032	94.4	$-$5.0	0.047	0.036	81.9
	MSI	$-$1.1	0.031	0.030	94.9	$-$1.8	0.050	0.031	93.6	$-$5.2	0.047	0.040	87.2
	IPW	$-$0.1	0.034	0.032	91.7	$-$1.3	0.039	0.045	97.1	$-$5.4	0.049	0.040	85.0
	PDR	$-$1.1	0.036	0.033	92.9	$-$1.7	0.053	0.036	94.1	$-$3.4	0.048	0.035	87.5
$n=1000$	FI	$-$0.8	0.020	0.020	94.1	$-$1.1	0.023	0.022	96.2	$-$5.2	0.031	0.025	69.1
	MSI	$-$0.9	0.020	0.021	95.1	$-$1.1	0.023	0.021	94.9	$-$5.3	0.031	0.029	78.0
	IPW	0.2	0.022	0.023	93.9	$-$1.0	0.024	0.035	97.5	$-$5.1	0.034	0.029	74.8
	PDR	$-$0.8	0.023	0.023	95.4	$-$0.9	0.028	0.024	95.1	$-$3.2	0.032	0.025	82.0

Table 5. Table 4: Application to IGF-II biomarker: estimated coefficients (Est.), standard errors (SE) and corresponding p 𝑝 p -values for (conditional) disease and verification models; VUS estimates, along with related standard errors and 95% confidence intervals.

	Disease model
	Benign disease			Early stage
	Est.	SE	$p$-value	Est.	SE	$p$-value
Intercept	$-$1.832	0.786	0.020	0.528	0.218	0.015
$T$	$-$0.922	0.204	0.001	$-$0.880	0.205	0.001
$A_1$	$-$1.289	0.444	0.004	$-$0.625	0.168	0.001
$A_2$	$-$17.131	4.325	0.001	$-$0.026	0.141	0.852
	Verification model
	Logit link			Probit link
	Est.	SE	$p$-value	Est.	SE	$p$-value
Intercept	6.186	1.455	0.001	2.966	0.621	0.001
$T$	$-$1.338	0.376	0.001	$-$0.659	0.162	0.001
$D_1$	$-$5.733	1.501	0.001	$-$2.649	0.639	0.001
$D_2$	4.207	0.961	0.001	1.899	0.631	0.003
	bias–corrected VUS
	Logit link			Probit link
	Est.	SE	95% CI	Est.	SE	95% CI
FI	0.317	0.025	$(0.268,0.366)$	0.317	0.025	$(0.268,0.366)$
MSI	0.343	0.028	$(0.287,0.398)$	0.343	0.028	$(0.287,0.398)$
IPW	0.243	0.087	$(0.073,0.413)$	0.299	0.056	$(0.188,0.410)$
PDR	0.342	0.027	$(0.289,0.395)$	0.341	0.030	$(0.283,0.399)$

Equations

$$\mu =\frac{\mathbb{E}\left({\mathrm{I}}_{i\ell r}{D}_{1i}{D}_{2\ell }{D}_{3r}\right)}{\mathbb{E}\left(D_{1i}{D}_{2\ell }{D}_{3r}\right)},$$

$${\hat{\mu }}_{\mathrm{NP}}=\frac{\sum _{i=1}^n\sum _{\ell =1,\ell \ne i}^n\sum _{\stackrel{r=1}{r\ne \ell ,r\ne i}}^n{\mathrm{I}}_{i\ell r}{D}_{1i}{D}_{2\ell }{D}_{3r}}{\sum _{i=1}^n\sum _{\ell =1,\ell \ne i}^n\sum _{\stackrel{r=1}{r\ne \ell ,r\ne i}}^n{D}_{1i}{D}_{2\ell }{D}_{3r}}.$$

$$\pi =\mathrm{Pr}(V=1|T,\mathbf{A},D_1,D_2)=\frac{\exp \left({\beta }_0+{\beta }_{1}T+{\mathbf{\beta }}_2^{\top }\mathbf{A}+{\lambda }_1{D}_1+{\lambda }_2{D}_2\right)}{1+\exp \left({\beta }_0+{\beta }_{1}T+{\mathbf{\beta }}_2^{\top }\mathbf{A}+{\lambda }_1{D}_1+{\lambda }_2{D}_2\right)},$$

$${\rho }_k=\frac{\exp \left({\eta }_0^{(k)}+{\eta }_1^{(k)}T+{\mathbf{A}}^{\top }{\mathbf{\eta }}_2^{(k)}\right)}{1+\exp \left({\eta }_0^{(1)}+{\eta }_1^{(1)}T+{\mathbf{A}}^{\top }{\mathbf{\eta }}_2^{(1)}\right)+\exp \left({\eta }_0^{(2)}+{\eta }_1^{(2)}T+{\mathbf{A}}^{\top }{\mathbf{\eta }}_2^{(2)}\right)},$$

$${\hat{\mu }}_{\ast }=\frac{\sum _{i=1}^n\sum _{\ell =1,\ell \ne i}^n\sum _{\stackrel{r=1}{r\ne \ell ,r\ne i}}^n{\mathrm{I}}_{i\ell r}{\hat{D}}_{1i,\ast }{\hat{D}}_{2\ell ,\ast }{\hat{D}}_{3r,\ast }}{\sum _{i=1}^n\sum _{\ell =1,\ell \ne i}^n\sum _{\stackrel{r=1}{r\ne \ell ,r\ne i}}^n{\hat{D}}_{1i,\ast }{\hat{D}}_{2\ell ,\ast }{\hat{D}}_{3r,\ast }}.$$

$$\bar{S}(\mathbf{\varphi })$$

$$\bar{S}(\mathbf{\eta })$$

$$\pi$$

$$\bar{S}(\mathbf{\gamma })$$

$$h(\mathbf{\gamma };V_i,T_i,{\mathbf{A}}_{1i},D_{1i},D_{2i})=Z_i\left\{V_i-\pi (T_i,{\mathbf{A}}_{1i},D_{1i},D_{2i};\mathbf{\gamma })\right\}$$

$${\mathbb{E}}_0\left\{h(\mathbf{\gamma };V_i,T_i,{\mathbf{A}}_{1i},D_{1i},D_{2i})|T_i,{\mathbf{A}}_i\right\}$$

$${\rho }_{1(0)}({\lambda }_1,{\lambda }_2)$$

$${\rho }_{2(0)}({\lambda }_1,{\lambda }_2)$$

$${\rho }_{3(0)}({\lambda }_1,{\lambda }_2)$$

$${\rho }_{k(1)}$$

$${\rho }_{1(0)}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)$$

$${\rho }_{2(0)}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)$$

$${\rho }_{3(0)}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)$$

$$\bar{S}(\mathbf{\gamma };\mathbf{\eta })$$

$$\mathbb{E}\left[V_i{\rho }_{k(1)i}(\mathbf{\eta })+(1-V_i){\rho }_{k(0)i}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)|T_i,{\mathbf{A}}_i\right]=\mathrm{Pr}(D_k=1|T_i,{\mathbf{A}}_i)$$

$$D_{ki,\mathrm{FI}}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)=V_i{\rho }_{k(1)i}(\mathbf{\eta })+(1-V_i){\rho }_{k(0)i}(\mathbf{\eta },{\lambda }_1,{\lambda }_2),$$

$$D_{ki,\mathrm{MSI}}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)=V_i{D}_{ki}+(1-V_i){\rho }_{k(0)i}(\mathbf{\eta },{\lambda }_1,{\lambda }_2),$$

$$D_{ki,\mathrm{IPW}}(\mathbf{\gamma })=\frac{V_i{D}_{ki}}{{\pi }_i(\mathbf{\gamma })}\mathrm{and}{D}_{ki,\mathrm{PDR}}(\mathbf{\eta },\mathbf{\gamma })=\frac{V_i{D}_{ki}}{{\pi }_i(\mathbf{\gamma })}-{\rho }_{k(0)i}(\mathbf{\eta },{\lambda }_1,{\lambda }_2)\frac{V_i-{\pi }_i(\mathbf{\gamma })}{{\pi }_i(\mathbf{\gamma })}.$$

$$G_{i\ell r,\text{FI}}(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$G_{i\ell r,\text{MSI}}(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$G_{i\ell r,\text{IPW}}(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$G_{i\ell r,\text{PDR}}(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$G_{\ast }(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$e_{\ast }(\mu ,\mathbf{\eta },\mathbf{\gamma })$$

$$\sqrt{n}\left({\tilde{\mu }}_{\ast }-{\mu }_0\right)\stackrel{d}{⟶}\mathcal{N}(0,{\sigma }_{\ast }^2),$$

$${\Lambda }_{i,\ast }$$

$$Q_{i,\ast }$$

$$0$$

$$\sqrt{n}{G}_{\ast }({\mu }_0,{\mathbf{\eta }}_0,{\mathbf{\gamma }}_0)\stackrel{d}{⟶}\frac{1}{\sqrt{n}}\sum _{i=1}^n{\Lambda }_{i,\ast }.$$

$$\frac{\partial e_{*}(\mu,\boldsymbol{\eta}_{0},\boldsymbol{\gamma}_{0})}{\partial\mu}\bigg{|}_{\mu=\mu_{0}}=-\Pr(D_{1}=1)\Pr(D_{2}=1)\Pr(D_{3}=1),$$

Full text

Key words: Diagnostic test; Instrumental variable; Mean score equation; Nonignorable missing data mechanism; ROC analysis; Verification bias.

1 Introduction and background

The evaluation of the accuracy of diagnostic tests (or biomarkers) is an increasingly relevant issue in modern medicine. Usually, this issue is solved by means of a receiver operating characteristic (ROC) analysis. In particular, the volume under the ROC surface (VUS) is often used for measuring the overall accuracy of a diagnostic test when the possible disease status belongs to one of three ordered categories. Under correct ordering, values of VUS vary from 1/6, suggesting that the test is no better than chance alone, to 1, that implies a perfect test, i.e. a test that perfectly discriminates among the three categories (Scurfield,, 1996).

In medical studies, the VUS of a new test is typically estimated through a sample of measurements obtained by some suitable sample of patients, for which the true disease status is assessed by means of a gold standard (GS) test. However, in many cases, due to the expensiveness and/or invasiveness of the GS test, only a subset of patients undergoes disease verification. In such situations, statistical inference based only on verified subjects is typically biased. This bias is known as verification bias.

In order to correct for verification bias, the researchers must formulate some assumptions about the selection mechanism for the disease verification. When the decision to send a subject to verification may be directly based on the presumed subject’s disease status, or, more generally, the selection mechanism may depend on some unobserved covariates related to disease, the missing data mechanism is called nonignorable (NI).

Generally speaking, NI data missingness is a challenge for inference, and, specifically, the issue of correcting for NI verification bias in ROC surface analysis is still scarcely considered in the statistical literature. To the best of our knowledge, only Zhang and Alonzo, (2018) and To Duc et al., (2019) gave a contribution in this direction. In particular, To Duc et al., (2019) extended the model discussed in Liu and Zhou, (2010) to match the case of three-category disease status, and proposed a likelihood-based approach to derive four bias–corrected estimators of the VUS, namely, full imputation (FI), mean score imputation (MSI), inverse probability weighting (IPW) and pseudo doubly robust (PDR) estimators.

Suppose we need to evaluate the predictive ability of a new continuous diagnostic test in a context where the disease status of a patient can be described by three ordered categories, “non–diseased”, “intermediate” and “diseased”, say. Consider a sample of $n$ subjects and let $T$, ${\mathcal{D}}$ and $\mathbf{A}$ denote the test result, the disease status and a vector of covariates for each subject, respectively. Hereafter, we assume, without loss of generality, that test values are positively associated with the severity degree of the disease. The disease status ${\mathcal{D}}$ can be modeled as a trinomial random vector ${\mathcal{D}}=(D_{1},D_{2},D_{3})^{\top}$, such that $D_{k}$ is a Bernoulli random variable having mean $\theta_{k}=\Pr(D_{k}=1)$ where $\theta_{1}+\theta_{2}+\theta_{3}=1$. Hence, $\theta_{k}$ represents the probability that a generic subject, classified according to its disease status, belong to the class $k$. The accuracy of $T$ is measured by its VUS, which is defined as

[TABLE]

where the indices $i$, $\ell$, $r$ refer to three different subjects, $\mathrm{I}_{i\ell r}=\mathrm{I}(T_{i}<T_{\ell}<T_{r})+\frac{1}{2}\mathrm{I}(T_{i}<T_{\ell}=T_{r})+\frac{1}{2}\mathrm{I}(T_{i}=T_{\ell}<T_{r})+\frac{1}{6}\mathrm{I}(T_{i}=T_{\ell}=T_{r})$ and $\mathrm{I}(\cdot)$ is the indicator function (Nakas and Yiannoutsos,, 2004). When the disease status $\mathcal{D}$ is available for all subjects, a natural nonparametric estimator of $\mu$ is given by

[TABLE]

Suppose now that the disease status ${\mathcal{D}}$ is missing for a subset of patients in the study. Let $V_{i}$ be the verification status for the $i$-th subject in study, that is $V_{i}=1$ if $\mathcal{D}_{i}$ is observed and $V_{i}=0$ otherwise. To account for a possible nonignorable missing data mechanism for the disease status, To Duc et al., (2019) considered the verification model

[TABLE]

where the parameters $\lambda_{1}$ and $\lambda_{2}$ describe the type of missingness mechanism. In particular, the missing data mechanism is ignorable (MAR) if $\lambda_{1}=\lambda_{2}=0$; NI, otherwise. Following Liu and Zhou, (2010), to estimate the parameter $\boldsymbol{\phi}=(\beta_{0},\beta_{1},\boldsymbol{\beta}^{\top}_{2},\lambda_{1},\lambda_{2})^{\top}$ in the verification model (1.3), To Duc et al., (2019) proposed to use a likelihood-based approach. Specifically, the authors assumed that the disease process follows a multinomial logistic model for the whole sample, i.e.,

[TABLE]

where $\rho_{k}=\Pr(D_{k}=1|T,\mathbf{A}),k=1,2$. Then, they jointly considered the verification model (1.3) and the disease model (1.4) to construct the (observed) log-likelihood function, $\ell(\boldsymbol{\phi},\boldsymbol{\eta})$ say, where $\boldsymbol{\eta}=(\eta^{(1)}_{0},\eta^{(1)}_{1},\boldsymbol{\eta}^{(1)\top}_{2},\eta^{(2)}_{0},\eta^{(2)}_{1},\boldsymbol{\eta}^{(2)\top}_{2})^{\top}$ is the vector of parameters of the disease model. The authors proved that the joint model based on (1.3) and (1.4) is identifiable. Then, the maximum likelihood estimator (MLE) of $(\boldsymbol{\phi}^{\top},\boldsymbol{\eta}^{\top})^{\top}$ can be obtained by maximizing $\ell(\boldsymbol{\phi},\boldsymbol{\eta})$ or solving the corresponding score equation. Based on imputation and re-weighting methods, four estimators for each disease status can be obtained: FI, with $\widehat{D}_{ki,\mathrm{FI}}=\hat{\rho}_{ki}=\widehat{\Pr}(D_{k}=1|T_{i},\mathbf{A}_{i})$; MSI, with $\widehat{D}_{ki,\mathrm{MSI}}=V_{i}D_{ki}+(1-V_{i})\hat{\rho}_{k(0)i}$, where $\hat{\rho}_{k(0)i}=\widehat{\Pr}(D_{ki}=1|T_{i},\mathbf{A}_{i},V_{i}=0)$; IPW, with $\widehat{D}_{ki,\mathrm{IPW}}=V_{i}D_{ki}/\hat{\pi}_{i}$ and PDR, with $\widehat{D}_{ki,\mathrm{PDR}}=V_{i}D_{ki}/\hat{\pi}_{i}-\hat{\rho}_{k(0)i}(V_{i}-\hat{\pi}_{i})/\hat{\pi}_{i}$, with $\hat{\pi}_{i}=\widehat{\Pr}(V_{i}=1|T_{i},\mathbf{A}_{i},D_{1i},D_{2i})$. Therefore, bias-corrected VUS estimators are derived by using quantities $\widehat{D}_{ki,*}$ in (1.2), i.e.,

[TABLE]

Here, the symbol $*$ stands for type of estimator, i.e., FI, MSI, IPW and PDR.

Despite its usefulness, the above reviewed approach presents some limitations. Firstly, it is scarcely flexible, since identifiability is proven only for the specific parametric assumptions about the verification and the disease processes. Secondly, because of NI missigness, it does not allow a statistical check of appropriateness of such assumptions (see also Riddles et al.,, 2016, Morikawa et al.,, 2017 and Yu et al.,, 2018). Finally, as pointed out by the authors themselves, MLE of ignorable/nonignorable parameters $\lambda_{1}$ and $\lambda_{2}$ are typically highly biased, in cases of moderate or even high sample sizes.

In the present paper, we try to overcome such limitations. More precisely, here we consider the mean score function as a tool to estimate the parameters in the verification model, and avoid identifiability problems by resorting to the introduction of instrumental variables (Wang et al.,, 2014; Riddles et al.,, 2016). As we show in what follows, in order to solve the mean score equation derived by the chosen verification model, we preliminarily need to estimate the parameters of a model for the disease process, but its specification is required only for verified subjects. Thus, the new proposed approach essentially avoids the specification of a parametric regression model for the disease status of the whole sample, and presents the advantage that it allows a statistical evaluation, based on the available sample data, of the adequacy of the required conditional (to $V=1$) disease model. Furthermore, as confirmed by our simulation experiments, the reduction in complexity of the adopted model translates into a greater stability of estimates of the parameters.

Based on the new approach, we provide four new VUS estimators (FI, MSI, IPW and PDR), which accommodate for nonignorable missingness of the disease status. We prove consistency and asymptotic normality of the proposed estimators, and also provide estimators of their variances.

The paper is organised as follows. In Section 2, we present the proposed approach, and give the new bias–corrected VUS estimators, for which, in Section 3, we discuss the asymptotic behaviour. In Section 4, we present the results of the simulation study, and in Section 5 we give an illustrative application. Section 6 contains some final remarks. Some technical details can be found in Appendix.

2 The proposal

It is worth noting that MLE of $(\boldsymbol{\phi}^{\top},\boldsymbol{\eta}^{\top})^{\top}$, on which estimators (1.5) ultimately rely, could be also obtained as solution of the so-called mean score equation (Louis,, 1982; Kim and Shao,, 2013; Riddles et al.,, 2016). Let $S(\boldsymbol{\phi})$ and $S(\boldsymbol{\eta})$ denote the score functions derived from the verification model (1.3) and from the disease model (1.4), respectively. For the joint model, the mean score function is $\bar{S}(\boldsymbol{\phi},\boldsymbol{\eta})=(\bar{S}(\boldsymbol{\phi}),\bar{S}(\boldsymbol{\eta}))^{\top}$ (Kim and Shao,, 2013), with

[TABLE]

where $h(\boldsymbol{\phi};\cdot)$ and $u(\boldsymbol{\eta};\cdot)$ denote the standard contributions to $S(\boldsymbol{\phi})$ and $S(\boldsymbol{\eta})$, respectively. Terms $\mathbb{E}_{0}(\cdot|T_{i},\mathbf{A}_{i})=\mathbb{E}(\cdot|T_{i},\mathbf{A}_{i},V_{i}=0)$, represent the contribution to the scores, arising from unverified subjects.

Hereafter, we consider an alternative estimation procedure which solves the mean score equation derived from the verification model only, and uses instrumental variables. To this aim, we assume that the vector of covariates $\mathbf{A}$ can be decomposed as $\mathbf{A}=\left(\mathbf{A}^{\top}_{1},\mathbf{A}^{\top}_{2}\right)^{\top}$ such that the verification status $V$ does not depend on $\mathbf{A}_{2}$ given $T$, $\mathbf{A}_{1}$, $D_{1}$ and $D_{2}$. Variable $\mathbf{A}_{2}$ is known as instrumental variable and helps to make the verification model identifiable (see Wang et al.,, 2014 and Riddles et al.,, 2016) without further assumptions on the disease model. By using the instrumental variable $\mathbf{A}_{2}$, we can reduce the verification model (1.3) to

[TABLE]

where $\boldsymbol{\gamma}=\left(\beta_{0},\beta_{1},\boldsymbol{\beta}^{\top}_{2,1},\lambda_{1},\lambda_{2}\right)^{\top}$. Then, MLE estimate of $\boldsymbol{\gamma}$ can be obtained by solving the mean score equation $\bar{S}(\boldsymbol{\gamma})=\boldsymbol{0}$ where

[TABLE]

with

[TABLE]

for $Z_{i}=(1,T_{i},\mathbf{A}_{1i}^{\top},D_{1i},D_{2i})^{\top}$ and

[TABLE]

Let $\rho_{k(v)}$ denote $\Pr(D_{k}=1|V=v,T,\mathbf{A})$ for $k=1,2,3$ and $v=0,1$. By an application of Bayes’ rule (see Appendix A), it can be shown that

[TABLE]

Hence, the mean score function (2.3) depends on the disease probabilities for verified units only. If (just to keep the parallel with the proposal in To Duc et al., (2019)) we use a multinomial logistic regression model for the conditional disease probabilities $\rho_{k(1)}=\Pr(D_{k}=1|V=1,T,\mathbf{A})$, i.e.,

[TABLE]

$k=1,2$, the functions $\rho_{k(0)}(\lambda_{1},\lambda_{2})$ depend also on $\boldsymbol{\eta}$, and the expressions in (2.5) become

[TABLE]

By using (2.4) and (2.7), the mean score function (2.3) can be rewritten as

[TABLE]

Then, once an estimate $\widehat{\boldsymbol{\eta}}$ is available, we can get an empirical form of the mean score function (2.3), say $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})$. The empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$ can be solved to obtain an estimate $\widehat{\boldsymbol{\gamma}}$ of $\boldsymbol{\gamma}$.

Observe that $\widehat{\boldsymbol{\gamma}}$ would be the MLE estimate only if $\widehat{\boldsymbol{\eta}}$ was solution of equation (2.1). It should be noted, however, that in our simplified approach, $\widehat{\boldsymbol{\eta}}$ is the solution of a score equation that depends on the parametric model specification for the conditional (to $V=1$) disease process. The verification mechanism may be modelled by means of a generalised linear regression model for binary data with a fixed link function (e.g., logistic, probit, log-log or complementary log-log). Whichever choice is taken, identifiability of the model holds thanks to the use of instrumental variables. In what follows, for the purposes of presentation, we adopt a standard logistic model for the verification process, and a multinomial logistic model for the conditional (to verification) disease process, a solution that also allow to keep the parallel with the approach in To Duc et al., (2019).

Let $U(\boldsymbol{\eta})=\sum\limits_{i=1}^{n}V_{i}u_{i}(\boldsymbol{\eta})$, with $u_{i}=u(\boldsymbol{\eta};T_{i},\mathbf{A}_{i},D_{1i},D_{2i})$, denote the estimating function for $\boldsymbol{\eta}$ based on the verified units in the sample. Since our empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$ is a specific case of a more general one discussed by Riddles et al., (2016), it follows that, under certain regularity conditions given in Appendix B, $\widehat{\boldsymbol{\gamma}}$ and $\widehat{\boldsymbol{\eta}}$ are consistent and asymptotically normal estimators (see Theorem 1 in Riddles et al.,, 2016).

In order to obtain four bias-corrected VUS estimators (FI, MSI, IPW and PDR), we follow the arguments in To Duc et al., (2019), starting by the fact that

[TABLE]

and considering the quantities

[TABLE]

[TABLE]

[TABLE]

in which we stressed the fact that the probabilities of selection for verification $\pi_{i}$ depend on $\boldsymbol{\gamma}$. Plugging in the above expressions the estimates $\widehat{\boldsymbol{\eta}}$ and $\widehat{\boldsymbol{\gamma}}$, yields the estimates $\widehat{D}_{ki,\mathrm{FI}}$, $\widehat{D}_{ki,\mathrm{MSI}}$, $\widehat{D}_{ki,\mathrm{IPW}}$ and $\widehat{D}_{ki,\mathrm{PDR}}$. Then, the bias-corrected VUS estimators are obtained by replacing each disease status in (1.2) with the new estimates $\widehat{D}_{ki,*}$, where the symbol $*$ indicates, again, the kind of estimator, i.e., FI, MSI, IPW and PDR.

Note that, the parameters $\boldsymbol{\gamma}$ of the verification model (2.2) and $\boldsymbol{\eta}$ of the conditional (to $V=1$) disease model (2.6) are estimated separately (see the estimation process in Appendix C). However, the accuracy of the estimates $\widehat{\pi}_{i}$ in the verification model (2.2) may suffer from possible fitting’s problems in the disease model (2.6), because the empirical mean score equation uses $\widehat{\boldsymbol{\eta}}$. On the other hand, the estimates $\widehat{\rho}_{k(0)i}$ also may be affected by possible poor accuracy in the estimation of ${\lambda}_{1}$, ${\lambda}_{2}$ or ${\boldsymbol{\eta}}$. Ultimately, the fact that the proposed method is a parametric method has as a consequence that it could produce VUS estimates even strongly influenced by incorrect specifications of the involved models.

3 Asymptotic behaviour

As for the estimators discussed in To Duc et al., (2019), the new proposed VUS estimators can be found as solutions of appropriate estimating equations, solved along with the score equations $U(\boldsymbol{\eta})=\boldsymbol{0}$ and the empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$. The estimating functions for FI, MSI, IPW and PDR estimators have generic term (corresponding to a generic triplet of sample units), respectively,

[TABLE]

for $i\neq\ell$, $\ell\neq r$ and $r\neq i$. In the following, we will use the general notation $G_{i\ell r,*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})$, where $*$ stands for FI, MSI, IPW and PDR, and denote by $\widetilde{\mu}_{*}$ the corresponding new VUS estimator. Moreover, let $\boldsymbol{\eta}_{0}$, $\boldsymbol{\gamma}_{0}$ and $\mu_{0}$ be the true parameter values. We assume that

1. (C1)

the U–process, $\sqrt{n}\left\{G_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})-e_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})\right\}$, is stochastically equicontinuous, where

[TABLE]

and

[TABLE] 2. (C2)

$e_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})$ is differentiable in $(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})$, and \dfrac{\partial e_{*}(\mu,\boldsymbol{\eta}_{0},\boldsymbol{\gamma}_{0})}{\partial\mu}\Bigg{|}_{\mu=\mu_{0}}\neq 0; 3. (C3)

$G_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})$ and $\frac{\partial G_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial(\boldsymbol{\eta},\boldsymbol{\gamma})^{\top}}$ converge uniformly (in probability) to $e_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})$ and $\frac{\partial e_{*}(\mu,\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial(\boldsymbol{\eta},\boldsymbol{\gamma})^{\top}}$, respectively.

Let $\mathbf{O}_{i}=(T_{i},\mathbf{A}^{\top}_{i},D_{1i},D_{2i},D_{3i},V_{i})^{\top}$, $\mathcal{I}_{\boldsymbol{\eta}}(\boldsymbol{\eta})=\mathbb{E}\left\{-\frac{\partial u_{i}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}^{\top}}\right\}$, $\mathcal{I}_{s,\boldsymbol{\gamma}}(\boldsymbol{\eta},\boldsymbol{\gamma})=\mathbb{E}\left\{-\frac{\partial s_{i}(\boldsymbol{\gamma};\boldsymbol{\eta})}{\partial\boldsymbol{\gamma}^{\top}}\right\}$ and $\mathcal{I}_{s,\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\gamma})=\mathbb{E}\left\{-\frac{\partial s_{i}(\boldsymbol{\gamma};\boldsymbol{\eta})}{\partial\boldsymbol{\eta}^{\top}}\right\}$.

Theorem 3.1**.**

Suppose that conditions (C1)–(C3), and (D1)–(D9) in Appendix B, hold. Then, under the verification model (2.2) and the (conditional) disease model (2.6), $\widetilde{\mu}_{*}$ is consistent and asymptotically normal, i.e.,

[TABLE]

where  $\sigma_{*}^{2}=\dfrac{\mathbb{V}\mathrm{ar}\left\{\Lambda_{i,*}+Q_{i,*}\right\}}{\left[\Pr(D_{1}=1)\Pr(D_{2}=1)\Pr(D_{3}=1)\right]^{2}}$  and

[TABLE]

Proof.

(sketch). Under the verification model (2.2) and the disease model (2.6), we can show that quantities $G_{i\ell r,*}(\mu_{0},\boldsymbol{\eta}_{0},\boldsymbol{\gamma}_{0})$ have 0 expectation. Then, the consistency of $\widetilde{\mu}_{*}$ can be proven by arguments similar to those used in Theorem 1 of To Duc et al., (2019). As for the asymptotic normality, we follow the steps in the proof of Theorem 2 in To Duc et al., (2019). In particular, using also equation (B.1), since $e_{*}(\mu_{0},\boldsymbol{\eta}_{0},\boldsymbol{\gamma}_{0})=0$, we can write

[TABLE]

An application of a standard result about the limit distribution of U-statistics (van der Vaart,, 2000), yields to

[TABLE]

This, together with the fact that

[TABLE]

and equation (B.1), show that (3.1) is equivalent to

[TABLE]

Then, we have that

[TABLE]

and the asymptotic normality of $\widetilde{\mu}_{*}$ follows by the Central Limit Theorem. ∎

The asymptotic variances $\sigma_{*}^{2}$ can be consistently estimated by

[TABLE]

where

[TABLE]

and

[TABLE]

Clearly, alternative approaches to estimate the variance of the proposed VUS estimators are possible. For instance, one may resort to the classical bootstrap method.

4 Simulation study

In order to compare the proposed mean score equation-based approach with the fully likelihood-based method in (To Duc et al.,, 2019) and to evaluate the performance of the new bias-corrected VUS estimators, we carry out several simulation experiments. The number of replications is 1000 in each experiment, and the considered sample sizes are $150,250,500$ and $1000$.

4.1 Simulation setup

With regard to data generation, we consider the following six different scenarios:

1. I.

For each unit, the test result $T$ and a covariate $A$ are generated by the model

[TABLE]

The disease status $\mathcal{D}=(D_{1},D_{2},D_{3})^{\top}$ is generated by the following multinomial logistic model:

[TABLE]

The verification status $V$ is simulated by the following logistic model

[TABLE]

Under this setting, the verification rate is roughly 0.57 and the true VUS is 0.791. This setting is the same as Scenario I in To Duc et al., (2019). 2. II.

For each unit, the disease status $\mathcal{D}=(D_{1},D_{2},D_{3})^{\top}$ is generated from a multinomial distribution with $\theta_{1}=0.7$, $\theta_{2}=0.2$ and $\theta_{3}=0.1$ (recall that $\theta_{k}=\Pr(D_{k}=1),$ $k=1,2,3$). The test results $T$ and a covariate $A$ are generated as $T|D_{k}=1\sim\mathcal{N}(k-1,0.5^{2})$ and $A|D_{k}=1\sim\mathcal{N}(0.5(k-1),0.5^{2})$, while the verification status $V$ is simulated through the model

[TABLE]

The verification rate is roughly 0.44 and the true VUS is 0.843. 3. III.

The disease status $\mathcal{D}$ is generated as in scenario I. The test results $T$ and a covariate $A_{1}$ are generated as $T|D_{k}=1\sim\mathcal{N}(0.4(k-1),0.5^{2})$ and $A_{1}|D_{k}=1\sim\mathcal{N}(0.5(k-1),0.5^{2})$, for $k=1,2,3$. A second covariate $A_{2}$ is simulated as $A_{2}|D_{1}=1\sim\mathcal{U}(-2,-1)$, $A_{2}|D_{2}=1\sim\mathcal{U}(-1,1)$ and $A_{2}|D_{3}=1\sim\mathcal{U}(1,2)$. The true VUS is 0.457. The verification status $V$ is simulated by a probit model

[TABLE]

where $\Phi(\cdot)$ denotes cumulative distribution function of a standard normal. In this case, the verification rate is roughly 0.47. 4. IV.

The disease status $\mathcal{D}$, the test result $T$ and a covariate $A$ are generated in the same way as in scenario II. The verification status $V$ is simulated through the model

[TABLE]

The verification rate is roughly 0.42, and the true VUS still is 0.843. 5. V.

We generate the test result $T\sim\mathcal{U}(-3,3)$, two covariates $A_{1}\sim\mathcal{N}(0,1)$ and $A_{2}\sim\mathrm{Ber}(0.6)$. The disease status $\mathcal{D}=(D_{1},D_{2},D_{3})^{\top}$ is generated by using the multinomial logistic model

[TABLE]

The verification status $V$ is simulated by the following model:

[TABLE]

In this scenario, the verification rate is roughly $0.56$. The true VUS value is $0.74$. 6. VI.

We generate the test result $T$, and two covariates $A_{1},A_{2}$ as in scenario IV. The disease status $\mathcal{D}=(D_{1},D_{2},D_{3})^{\top}$ is generated by using the multinomial logistic model

[TABLE]

The verification status $V$ is simulated by the following model:

[TABLE]

The verification rate is roughly $0.46$, and the true VUS value is $0.728$.

Scenarios I and II aim to compare the discussed mean score equation-based approach with the fully likelihood-based one, in terms of accuracy of parameter estimators in the verification model, and VUS estimators. As pointed out above, scenario I is identical to Scenario I in To Duc et al., (2019) and does not involve instrumental variables. On the contrary, in scenario II the covariate $A$ plays the role of instrumental variable. Moreover, scenarios II to VI aim to assess the finite-sample behavior of new proposed bias-corrected VUS estimators. In particular, as explained below, scenarios IV, V and VI serve to evaluate the effects on the estimators of some types of misspecifications in the working models.

By applying the Bayes’ rule, for all considered scenarios, the true conditional (to $V=1$) disease process takes the following form:

[TABLE]

where $\mathbf{A}$ is simply the covariate $A$ in the scenarios I, II, IV and V, $\mathbf{A}=(A_{1},A_{2})^{\top}$ in scenarios III and VI and $\pi_{10}$, $\pi_{01}$ and $\pi_{00}$ denote $\Pr(V=1|T,\mathbf{A},D_{1}=1,D_{2}=0)$, $\Pr(V=1|T,\mathbf{A},D_{1}=0,D_{2}=1)$ and $\Pr(V=1|T,\mathbf{A},D_{1}=0,D_{2}=0)$, respectively. Clearly, in practice this is a multinomial logistic model with unknown offset terms (intercepts).

In making simulations, when considering the new mean score equation-based approach, we fit the the conditional disease model as a multinomial logistic one, i.e., we specify the working conditional disease model as follows:

[TABLE]

where $\mathbf{A}=A$ for scenarios I, II and IV, $\mathbf{A}=(A_{1},A_{2})^{\top}$ for scenarios III and V and $\mathbf{A}=(A_{1}^{2},A_{2})^{\top}$ for scenario VI. As for the verification model, in the estimation procedure we use the correct logistic model in scenarios I, II, III and VI, while in scenarios IV we use a logistic models involving only $T$, and in scenario V we use a probit model involving only $T$ and the covariate $A_{1}.$ Therefore, the working models are correctly specified in scenarios I, II and III, while they present various types of misspecification in scenarios IV, V and VI. More precisely: in scenario IV, we omit the covariate $A$ in the verification model and use it as instrumental variable; in scenario V, we omit the covariate $A_{2}$ in the verification model, use it as instrumental variable and incorrectly define the model (probit instead of logit); in scenario VI, we use $A_{1}^{2}$ instead of $A_{1}$ and omit the interaction term in the conditional disease model. Finally, as for the fully likelihood-based approach, in scenarios I and II, the working (unconditional) disease model and the working verification model are both correctly specified.

4.2 Computational aspects

In order to obtain the new VUS estimates, in our simulation we need to solve, at each replication, the score equation $U(\boldsymbol{\eta})=\boldsymbol{0}$, and then the empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$. The first equation refers to the conditional (to $V=1$) disease model and, since we work with a multinomial logistic model in all scenarios, we resort to the routine multinom() (of R package nnet) to solve it. Given $\widehat{\boldsymbol{\eta}}$, the empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$ could be solved in $\boldsymbol{\gamma}$ by using some numeric algorithm, such the ones deployed in the R routine nleqslv(). However, numerical algorithms for solving equations are often poorly stable and reliable, or, better, are less reliable than optimization algorithms. Therefore, in order to obtain $\widehat{\boldsymbol{\gamma}}$, we minimize the squared Euclidean norm of $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})$, i.e., $\left\|\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})\right\|^{2}$. The optimization is performed by using the L-BFGS-B algorithm supplied in the R routine optim().

4.3 Results

Simulation results are given in tables 1-3. Table 1 refers to the comparison between the mean score equation-based approach (MSEq) and the fully likelihood-based approach (logLike). For three different sample sizes, i.e., 150, 250 and 500, the table provides Monte Carlo means for the estimators of the parameters in the verification model, together with Monte Carlo means for the four (FI, MSI, IPW and PDR) VUS estimators, for scenarios I and II. Recall that, in such scenarios, the considered working models are correctly specified. The results clearly show that the new approach achieves better performance, providing less biased estimators, already at the smallest sample size.

Results in Table 2 allow us to evaluate the behavior of the proposed VUS estimators, and of the corresponding variance estimators based on (3.2), when the considered working models are correctly specified. For both scenarios II and III and four different sample sizes (150, 250, 500 and 1000), the table contains Monte Carlo relative biases, Monte Carlo standard deviations and (means of) estimated standard deviations based on (3.2), for the new VUS estimators. The table also provides the empirical coverages of the 95% confidence intervals for the VUS, obtained through the normal approximation approach applied to each estimator and the use of the corresponding variance estimator (3.2).

Overall, Table 2 shows that the proposed approach provides satisfactory results in terms of relative bias and coverage probability when the verification models are correctly specified. At the smallest sample size, the larger absolute relative bias is equal to $1.4\%$ (scenario III). As expected, when the sample size increases, the relative bias of the VUS estimators decreases, becoming negligible for $n\geq 250$. Moreover, especially for scenario III (in which the verification process follows a probit model), in case of smallest sample size, $n=150$, the (means of) estimated standard deviations based on (3.2) still seem to be affected by a poor accuracy in the estimation of the parameters in the working models. However, in similar cases, as shown in the table, more accurate estimates of standard deviations can be obtained via bootstrap. Bootstrap estimates also improve empirical coverage probabilities of confidence intervals, when used within the normal approximation approach (see also the results for scenario II).

We emphasize that, when working models are correctly specified, low-accuracy problems can be related to low verification rates. For example, in scenario III the verification rate is 0.47. If we change the verification process and use the following model to simulate $V_{i}$ values:

[TABLE]

the verification rate becomes roughly 0.72, and the corresponding simulation results (for $n=150$ and 1000 replications), given below,

are clearly better than those in Table 2.

[FIGURE:]

Table 3 presents the simulation results for scenarios IV-VI, in which the working models are, in various ways, misspecified. Again, for each scenario and four different sample sizes, the table gives Monte Carlo relative bias, Monte Carlo standard deviations and (means of) estimated standard deviations based on (3.2), for the new VUS estimators. The table also provides the empirical coverages of the 95% confidence intervals for the VUS, obtained through the normal approximation approach applied to each estimator and the use of the corresponding variance estimator (3.2).

In the considered scenarios, all proposed VUS estimators appear to be biased, even when the sample size increases. However, simulation results appear to stay on acceptable levels (the relative biases are not too large and the coverage probabilities are close to 95%) in scenarios IV and V, i.e., when the misspecification follows by the omission of a variable and/or a wrong choice of the link function in the working generalized linear model for the verification process. Clearly, due to the strong deviation between the working conditional disease model and the true one, the results are significantly worse in the scenario VI. Anyway, variance estimators (3.2) seem relatively more robust with respect to the considered misspecifications in the working models.

[FIGURE:]

5 Application

We illustrate an application of our method by using data from the Pre-PLCO Phase II Dataset, from the SPORE/Early Detection Network/Prostate, Lung, Colon, and Ovarian Cancer Study. The study protocol and the data are publicly available at the address 111http://edrn.nci.nih.gov/protocols/119-spore-edrn-pre-plco-ovarian-phase-ii-validation. In this study, several serum protein biomarkers are evaluated in terms of their ability to correctly classify possible epithelial ovarian cancer (EOC) cases into benign disease, early stage (I-II) and late stage (III-IV). Among such biomarkers, we are interested in insulin-like growth factor-II (IGF-II), previously studied in Mor et al., (2005) and Visintin et al., (2008) to distinguish non-diseased from cancer cases (early stage and the late stage). The data we consider refer to 156 patients with benign disease, 71 patients with early stage, and 82 with late stage. Moreover, the data present 87 unverified units; therefore, only 78% of patients receive the assessment of their true disease status.

The boxplots for the IGF-II values, stratified by the disease status and based on the verified units only, are presented in Figure 1, which show that lower values of IGF-II are associated with higher severity of disease; thus, in our analysis, we consider minus the values of the IGF-II as test values. Moreover, we consider two other serum proteins, i.e., HE4 and CA125, as covariates.

Our goal is to estimate the VUS for the IGF-II biomarker. In the analysis, we standardize the values of IGF-II, HE4 and CA125, and below, for convenience, denote them as $T$, $A_{1}$ and $A_{2}$, respectively. We also let $D_{1}$ and $D_{2}$ to be the variables corresponding to benign disease and early stage, respectively.

We fit the conditional (to $V=1$) disease model by using a multinomial logistic regression model with $T$, $A_{1}$ and $A_{2}$ as predictors and late stage as reference level. As for the verification model, we consider a logistic model and a probit model, with only $T$ as a predictor; hence, $A_{1}$ and $A_{2}$ are used as instrumental variables. The estimated parameters, estimated standard errors and corresponding $p$-values for the significance test for coefficients in the models are given in Table 4. The $p$-values refer to the Wald-type statistics, which use results in (B.4). The resulting $p$-values for ignorable/nonignorable coefficients in both models (logistic and probit), indicate that, in such case, the missing data mechanism is nonignorable.

Table 4 also provides the four bias-corrected VUS estimates, along with the standard errors and 95% confidence intervals (CI), constructed using normal approximation and (3.2). The Naïve VUS estimate, calculated from verified subjects only, is 0.460 with standard error 0.035 and 95% CI $(0.393,0.528)$. Therefore, all bias-corrected VUS estimates are significantly different to the Naïve one. Moreover, unlike the logistic model, the IPW estimate in probit model is comparable to the other VUS estimates. This may suggest that the probit link function is a good choice for the generalized regression model describing the verification process for these data.

Finally, we highlight the purely illustrative nature of this application and the fact that, in real applications, the approach proposed in this paper allows to evaluate the model chosen to describe the (conditional) disease process by diagnostic procedures and/or ad hoc statistical tests present in the litterature, as those, for instance, in Goeman and le Cessie, (2006) for the multinomial logistic regression model.

6 Conclusion

In this paper, we propose a mean score equation-based approach to estimate the volume under a ROC surface when the disease status is missing not at random. This approach can take advantage of the use of instrumental variables, which help avoid possible identifiability problems. We prove that our proposed bias–corrected VUS estimators are consistent and asymptotically normal.

Our method is essentially a parametric method, and requires the estimation of the parameters of a (parametric) model chosen to describe the conditional (to $V=1$) disease process, and then the solution of an empirical mean score equation resulting from a parametric model chosen for the verification process. The method may also be used in a semiparametric context, by resorting to a nonparametric regression approach to fit the conditional disease model, as in Morikawa et al., (2017). Clearly, this topic deserves further investigation.

Compared to the fully likelihood-based approach discussed in To Duc et al., (2019), the new method is simpler, and more flexible also thanks to the possible use of instrumental variables. Moreover, as our simulation results show, the new bias–corrected VUS estimators are generally more accurate with moderate sample sizes.

Resorting to instrumental variables may be necessary when relying on working models other than those, namely logistic and multinomial logistic, assumed by To Duc et al., (2019) to describe, respectively, the verification and disease process. In practice, however, the choice of variables that could play the role of instrumental variables is not trivial. With respect to this issue, we suggest a possible selection strategy that foresees three steps:

1. (i)

choose a working model for the conditional (to $V=1$) disease process and select (by using a standard backward stepwise regression) a vector of covariates $\mathbf{A}$, statistically significant; 2. (ii)

choose a working model for the verification process and, given the estimate $\widehat{\boldsymbol{\eta}}$ obtained in step (i), use, one at a time, the elements of $\mathbf{A}$, jointly with the test $T$, to solve the resulting empirical mean score equation $\bar{S}(\boldsymbol{\gamma};\widehat{\boldsymbol{\eta}})=\boldsymbol{0}$; 3. (iii)

take as vector of instrumental variables, the vector $\mathbf{A}_{2}$ consisting of the elements of $\mathbf{A}$ resulting statistically non-significant at the previous step (ii).

Then, the final estimate $\widehat{\boldsymbol{\gamma}}$ is obtained by solving the empirical mean score equation based on $T$ and, when $\mathbf{A}_{2}$ does not coincide with $\mathbf{A}$, the component $\mathbf{A}_{1}$ of $\mathbf{A}=\left(\mathbf{A}^{\top}_{1},\mathbf{A}^{\top}_{2}\right)^{\top}$. This is also the selection strategy that we used in the illustrative application of Section 5.

Appendix A Conditional probabilities of disease for unverified subjects

Applying Bayes’ rule, we have

[TABLE]

for $v=0,1$. Thus, we obtain the following ratios

[TABLE]

where

[TABLE]

with $(d_{1},d_{2})\in\left\{(1,0),(0,1),(0,0)\right\}$. These expressions are equivalent to the system of equations

[TABLE]

By using the fact that $\rho_{3(0)}=1-\rho_{1(0)}-\rho_{2(0)}$, from (A.1a), we obtain

[TABLE]

Substituting (A.2) into (A.1b), leads to the following equation

[TABLE]

After some algebra, we obtain that $\rho_{2(0)}$ is equal to

[TABLE]

and then $\rho_{1(0)}$ is

[TABLE]

In case of the logistic model (2.2), it is straightforward to show that

[TABLE]

Thus, we get the results in (2.5).

Appendix B Conditions for asymptotic normality of $\left(\widehat{\boldsymbol{\eta}}^{\top},\widehat{\boldsymbol{\gamma}}^{\top}\right)^{\top}$

Recall that $\mathcal{I}_{\boldsymbol{\eta}}(\boldsymbol{\eta})=\mathbb{E}\left\{-\frac{\partial u_{i}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}^{\top}}\right\}$, $\mathcal{I}_{s,\boldsymbol{\gamma}}(\boldsymbol{\eta},\boldsymbol{\gamma})=\mathbb{E}\left\{-\frac{\partial s_{i}(\boldsymbol{\gamma};\boldsymbol{\eta})}{\partial\boldsymbol{\gamma}^{\top}}\right\}$, $\mathcal{I}_{s,\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\gamma})=\mathbb{E}\left\{-\frac{\partial s_{i}(\boldsymbol{\gamma};\boldsymbol{\eta})}{\partial\boldsymbol{\eta}^{\top}}\right\}$ and consider the conditions:

1. (D1)

$\Pr\left(D_{k}=1|T,\mathbf{A}_{1},\mathbf{A}_{2}=\boldsymbol{a}_{2}^{(1)}\right)\neq\Pr\left(D_{k}=1|T,\mathbf{A}_{1},\mathbf{A}_{2}=\boldsymbol{a}_{2}^{(2)}\right)$ for given two values $\boldsymbol{a}_{2}^{(1)}\neq\boldsymbol{a}_{2}^{(2)}$ of the (instrumental) variable $\mathbf{A}_{2}$, and

[TABLE]

with $\mathbf{A}=\left(\mathbf{A}^{\top}_{1},\mathbf{A}^{\top}_{2}\right)^{\top}$; 2. (D2)

$\pi(T_{i},\mathbf{A}_{1i},D_{1i},D_{2i};\boldsymbol{\gamma})$ is a continuous function of $\boldsymbol{\gamma}$ with first and second derivatives continuous in an open set containing $\boldsymbol{\gamma}_{0}$ as an interior point; 3. (D3)

the variables $(T_{i},\mathbf{A}^{\top}_{i},D_{1i},D_{2i},D_{3i},V_{i})^{\top}$, $i=1,\ldots,n$, are independent and identically distributed; 4. (D4)

the parameter spaces $\Theta$ for $\boldsymbol{\eta}$ and $\Omega$ for $\boldsymbol{\gamma}$ are compact and have finite dimension; 5. (D5)

$\pi(T_{i},\mathbf{A}_{1i},D_{1i},D_{2i};\boldsymbol{\gamma})$ is bounded away from 0. 6. (D6)

the MLE , $\widehat{\boldsymbol{\gamma}}_{n}$ say, is consistent and asymptotically normal; 7. (D7)

$\mathbb{E}\left\{U(\boldsymbol{\eta})\right\}=\boldsymbol{0}$ has a unique solution $\boldsymbol{\eta}_{0}\in\Theta$; 8. (D8)

$\mathcal{I}_{\boldsymbol{\eta}}(\boldsymbol{\eta}_{0})$ exists and is invertible; 9. (D9)

there exists a neighborhood $\mathcal{N}$ of $\boldsymbol{\eta}_{0}$ such that the quantiles $\sup_{\boldsymbol{\eta}\in\mathcal{N}}\|U(\boldsymbol{\eta})\|$, $\sup_{\boldsymbol{\eta}\in\mathcal{N}}\|\frac{\partial}{\partial\boldsymbol{\eta}^{\top}}U(\boldsymbol{\eta})\|$ and $\sup_{\boldsymbol{\eta}\in\mathcal{N}}\|U(\boldsymbol{\eta})U(\boldsymbol{\eta})^{\top}\|$ have finite expectations, with $\|\mathbf{X}\|\equiv\sum_{i}\sum_{j}X_{ij}^{2}$.

By Theorem 1 in Riddles et al., (2016), under the regularity conditions (D1)–(D9), the solution $\left(\widehat{\boldsymbol{\eta}}^{\top},\widehat{\boldsymbol{\gamma}}^{\top}\right)^{\top}$ of

[TABLE]

is a consistent estimator of $\left(\boldsymbol{\eta}^{\top}_{0},\boldsymbol{\gamma}^{\top}_{0}\right)^{\top}$, and satisfies

[TABLE]

Thus, by the Central Limit Theorem, $\left(\widehat{\boldsymbol{\eta}}^{\top},\widehat{\boldsymbol{\gamma}}^{\top}\right)^{\top}$ is asymptotically normal. Moreover,

[TABLE]

where

[TABLE]

and

[TABLE]

where $\Sigma_{\boldsymbol{\eta}}=\mathcal{I}^{-1}_{\boldsymbol{\eta}}(\boldsymbol{\eta}_{0})\mathrm{Cov}(u_{i}(\boldsymbol{\eta}_{0}))\mathcal{I}^{-1}_{\boldsymbol{\eta}}(\boldsymbol{\eta}_{0})$. Consistent estimators $\widehat{\Sigma}_{\boldsymbol{\gamma}}$ of $\Sigma_{\boldsymbol{\gamma}}$ and $\widehat{\Sigma}_{\boldsymbol{\eta}}$ of $\Sigma_{\boldsymbol{\eta}}$ can be obtained by

[TABLE]

with $f_{i}(\widehat{\boldsymbol{\eta}},\widehat{\boldsymbol{\gamma}})=s_{i}(\widehat{\boldsymbol{\gamma}};\widehat{\boldsymbol{\eta}})-\left\{\sum_{i=1}^{n}\frac{\partial}{\partial\boldsymbol{\eta}^{\top}}s_{i}(\widehat{\boldsymbol{\gamma}};\widehat{\boldsymbol{\eta}})\right\}\left\{\sum_{i=1}^{n}\frac{\partial}{\partial\boldsymbol{\eta}^{\top}}V_{i}u_{i}(\widehat{\boldsymbol{\eta}})\right\}^{-1}V_{i}u_{i}(\widehat{\boldsymbol{\eta}})$.

Appendix C Estimation Procedure

Modelling Verification Process

Modelling Disease Process for Verified Subjects (eventually using instrumental veriables)

$\widehat{\boldsymbol{\eta}}$

${\rho}_{k(0)i}(\widehat{\boldsymbol{\eta}},\lambda_{1},\lambda_{2})$

Empirical Mean Score Equation

$\widehat{\beta}_{0}$, $\widehat{\beta}_{1}$, $\widehat{\boldsymbol{\beta}}_{2,1}$, $\widehat{\lambda}_{1}$, $\widehat{\lambda}_{2}$

${\rho}_{k(1)i}(\widehat{\boldsymbol{\eta}})$

${\rho}_{k(0)i}(\widehat{\boldsymbol{\eta}},\widehat{\lambda}_{1},\widehat{\lambda}_{2})$

$\widehat{\pi}_{i}$

VUS-FI Estimator

VUS-MSI Estimator

VUS-PDR Estimator

VUS-IPW Estimator