Confidence Intervals and Sample Size to Compare the Predictive Values of Two Diagnostic Tests

Abstract

A binary diagnostic test is a medical test that is applied to an individual in order to determine the presence or the absence of a certain disease and whose result can be positive or negative. A positive result indicates the presence of the disease, and a negative result indicates the absence. Positive and negative predictive values represent the accuracy of a binary diagnostic test when it is applied to a cohort of individuals, and they are measures of the clinical accuracy of the binary diagnostic test. In this manuscript, we study the comparison of the positive (negative) predictive values of two binary diagnostic tests subject to a paired design through confidence intervals. We have studied confidence intervals for the difference and for the ratio of the two positive (negative) predictive values. Simulation experiments have been carried out to study the asymptotic behavior of the confidence intervals, giving some general rules for application. We also study a method to calculate the sample size to compare the parameters using confidence intervals. We have written a program in R to solve the problems studied in this manuscript. The results have been applied to the diagnosis of colorectal cancer.

1. Introduction

A diagnostic test is medical test that is applied to an individual in order to determine the presence of a certain disease. Binary diagnostic tests are a very common type of diagnostic test in clinical practice. A binary diagnostic test (BDT) is a diagnostic test whose possible result is positive or negative. A positive result indicates the presence of the disease, and a negative result indicates the absence. Mammography for the diagnosis of breast cancer is an example of BDT. Accuracy of a BDT is measured in terms of two fundamental parameters: sensitivity and specificity. Sensitivity (Se) is the probability of the result of the BDT being positive when the individual has the disease, and specificity (Sp) is the probability of the result of the BDT being negative when the individual does not have the disease. Therefore, Se and Sp are probabilities of getting the disease diagnosis right, and they represent the intrinsic accuracy of the BDT, since these parameters depend on the physical, chemical, or biological properties upon which the BDT is developed. Other parameters that are used to assess and compare two BDTs are the positive and negative predictive values. Positive predictive value ($\tau$) is the probability of an individual having the disease when the result of the BDT is positive, and the negative predictive value ($\upsilon$) is the probability of an individual not having the disease when the result of the BDT is negative. Predictive values represent the accuracy of the diagnostic test when it is applied to a cohort of individuals, and they are measures of the clinical accuracy of the BDT. Predictive values depend on Se, Sp, and on the disease prevalence (p), and are easily calculated applying Bayes theorem, i.e.,

$$\tau =\frac{p\times Se}{p\times Se+\left(1-p\right)\times \left(1-Sp\right)}and\upsilon =\frac{\left(1-p\right)\times Sp}{p\times \left(1-Se\right)+\left(1-p\right)\times Sp}.$$

(1)

The accuracy of a BDT is assessed in relation to a gold standard. A gold standard (GS) is a medical test that determines without error whether or not an individual has the disease. Biopsy for the diagnosis of breast cancer is an example of GS.

On the other hand, the comparison of parameters of two BDTs is an important topic in the study of statistical methods for diagnosis in Medicine. The most frequent sample design to compare the parameters of two BDTs is the paired design. The paired design consists of applying the two BDTs and the GS to all individuals of a random sample sized n. The comparison of the predictive values of two BDTs subject to a paired design has been the subject of different studies. Bennett [1,2], Leisenring et al. [3], Wang et al. [4], Kosinski [5], Tsou [6], and Takahashi and Yamamoto [7] have studied hypothesis tests to compare the two positive predictive values and the negative predictive values independently. Roldán-Nofuentes et al. [8] studied a global hypothesis test to simultaneously compare the positive and negative predictive values of two BDTs. However, the comparison of predictive values using confidence intervals has been little studied. If the hypothesis test is significant to an $\alpha$ error, the confidence interval (CI) allows determining how much one predictive value is greater than the other. Moskowitz and Pepe [9] have proposed a Wald-type CI for the ratio of the two positive (negative) predictive values. A Wald-type CI for the difference of the two positive (negative) predictive values is easily obtained by inverting the contrast statistic of the hypothesis test studied by Wang et al. [4].

The objective of this manuscript is study CIs to compare the positive (negative) predictive values of two BDTs subject to a paired design. For this, we have studied CIs for the difference and for the ratio of the two positive (negative) predictive values. If a CI for the difference (ratio) does not contain the zero (one) value, then we reject the equality between the two positive (negative) predictive values, and we estimate how much bigger one predictive value is than another one. The problem of calculating the sample size to compare the two positive (negative) predictive values through a CI is also studied.

This manuscript is structured in the following way. In Section 2, the existing CIs are presented and other new CIs are proposed, both for the ratio and for the difference between the positive (negative) predictive values. In Section 3, simulation experiments are carried out to study the coverage probabilities and the average lengths of the CIs. In Section 4, a method to calculate the sample size to compare the parameters through CIs is proposed. In Section 5, we present the program called “cicpvbdt”, which is a program written in R that solves the problems studied in this manuscript. In Section 6, the results were applied to an example on the diagnosis of colorectal cancer, and in Section 7, the results are discussed.

2. Confidence Intervals

Let us consider two BDTs that are assessed in relation to the same GS. Let $T_i$ be the variable that models the result of the ith BDT, with $i=1,2$. $T_i=0$ indicates that the test result is negative, and $T_i=1$ indicates that the test result is positive. Let D be the random variable that models the result of the GS, so that $D=1$ when the individual is diseased and $D=0$ when the individual is non-diseased. Let $S{e}_i$ and $S{p}_i$ be the sensitivity and specificity of the ith BDT.

Table 1. Observed frequencies subject to a paired design.

	${\mathit{T}}_1=1$		${\mathit{T}}_1=0$
	${\mathit{T}}_2=1$	${\mathit{T}}_2=0$	${\mathit{T}}_2=1$	${\mathit{T}}_2=0$	Total
$\mathrm{Diseased}(D=1$)	$s_{11}$	$s_{10}$	$s_{01}$	$s_{00}$	$s$
$\mathrm{Non}\text{-}\mathrm{diseased}(D=0$)	$r_{11}$	$r_{10}$	$r_{01}$	$r_{00}$	$r$
Total	$s_{11}+r_{11}$	$s_{10}+r_{10}$	$s_{01}+r_{01}$	$s_{00}+r_{00}$	$n$

shows the observed frequencies obtained subject to a paired design. The frequencies $s_{jk}$ and $r_{jk}$ are the product of a multinomial distribution whose probabilities are $p_{jk}=P\left(D=1,T_1=j,T_2=k\right)$ and $q_{jk}=P\left(D=0,T_1=j,T_2=k\right)$, with $j,k=0,1$. Applying the conditional dependence model of Vacek [10], probabilities $p_{jk}$ and $q_{jk}$ are written as

$$p_{jk}=p\left[S{e}_1^j{\left(1-S{e}_1\right)}^{1-j}S{e}_2^k{\left(1-S{e}_2\right)}^{1-k}+{\delta }_{jk}{\epsilon }_1\right]$$

(2)

and

$$q_{jk}=q\left[S{p}_1^{1-j}{\left(1-S{p}_1\right)}^{j}S{p}_2^{1-k}{\left(1-S{p}_2\right)}^k+{\delta }_{jk}{\epsilon }_0\right],$$

(3)

where ${\delta }_{jk}=1$ if $j=k$ and ${\delta }_{jk}=-1$ if $j\ne k$, with $j,k=0,1$, ${\epsilon }_1$ is the dependence factor between the two BDTs when $D=1$ and ${\epsilon }_0$ is the dependence factor between the two BDTs when $D=0$. It is verified that $0\le {\epsilon }_1\le Min\left\{S{e}_1\left(1-S{e}_2\right),S{e}_2\left(1-S{e}_1\right)\right\}$ and $0\le {\epsilon }_0\le Min\left\{S{p}_1\left(1-S{p}_2\right),S{p}_2\left(1-S{p}_1\right)\right\}$. If ${\epsilon }_1={\epsilon }_0=0$, then the two BDTs are conditionally independent on the disease. This assumption is not realistic, so in practice, it is verified that ${\epsilon }_1>0$ and/or ${\epsilon }_0>0$. Let $\pi ={\left(p_{11},p_{10},p_{01},p_{00},q_{11},q_{10},q_{01},q_{00}\right)}^T$ be the vector of probabilities of the multinomial distribution, $p={\displaystyle \sum _{i,j=0}^1{p}_{ij}}$ and $q=1-p={\displaystyle \sum _{i,j=0}^1{q}_{ij}}$. The maximum likelihood estimators of $p_{jk}$ and $q_{jk}$ are:

$${\widehat{p}}_{jk}=\frac{s_{jk}}{n}and{\widehat{q}}_{jk}=\frac{r_{jk}}{n}.$$

From Equation (1), the sensitivity and specificity of each BDT are written, in terms of the predictive values and of p, as

$$S{e}_i=\frac{{\tau }_i\left({\upsilon }_i-q\right)}{p{Y}_i}andS{p}_i=\frac{{\upsilon }_i\left({\tau }_i-p\right)}{q{Y}_i}$$

(4)

where $q=1-p$ and $Y_i={\tau }_i+{\upsilon }_i-1$. Then

$$0\le {\epsilon }_1\le Min\left\{\frac{{\tau }_1\left({\upsilon }_1-q\right)}{p{Y}_1}\left(1-\frac{{\tau }_2\left({\upsilon }_2-q\right)}{p{Y}_2}\right),\frac{{\tau }_2\left({\upsilon }_2-q\right)}{p{Y}_2}\left(1-\frac{{\tau }_1\left({\upsilon }_1-q\right)}{p{Y}_1}\right)\right\}$$

(5)

and

$$0\le {\epsilon }_0\le Min\left\{\frac{{\upsilon }_1\left({\tau }_1-p\right)}{q{Y}_1}\left(1-\frac{{\upsilon }_2\left({\tau }_2-p\right)}{q{Y}_2}\right),\frac{{\upsilon }_2\left({\tau }_2-p\right)}{q{Y}_2}\left(1-\frac{{\upsilon }_1\left({\tau }_1-p\right)}{q{Y}_1}\right)\right\}.$$

(6)

In terms of predictive values, Equations (2) and (3) are written as

$$p_{jk}=p\left[\frac{{\tau }_1^j{\left({\upsilon }_1-q\right)}^j{\left({\tau }_1+p\right)}^{1-j}{\left(1-{\upsilon }_1\right)}^{1-j}}{p^j{p}^{1-j}{Y}_1^j{Y}_1^{1-j}}\times \frac{{\tau }_2^k{\left({\upsilon }_2-q\right)}^k{\left({\tau }_2+p\right)}^{1-k}{\left(1-{\upsilon }_2\right)}^{1-k}}{p^k{p}^{1-k}{Y}_2^k{Y}_1^{1-k}}+{\delta }_{jk}{\epsilon }_1\right]$$

(7)

and

$$q_{jk}=q\left[\frac{{\left(1-{\tau }_1\right)}^j{\left({\upsilon }_1-q\right)}^j{\left({\tau }_1-p\right)}^{1-j}{\upsilon }_1^{1-j}}{q^j{q}^{1-j}{Y}_1^j{Y}_1^{1-j}}\times \frac{{\left(1-{\tau }_2\right)}^k{\left({\upsilon }_2-q\right)}^k{\left({\tau }_2-p\right)}^{1-k}{\upsilon }_2^{1-k}}{q^k{q}^{1-k}{Y}_2^k{Y}_2^{1-k}}+{\delta }_{jk}{\epsilon }_0\right].$$

(8)

The estimator of sensitivities and specificities are

$$\widehat{S}{e}_1=\frac{s_{11}+s_{10}}{s},\widehat{S}{e}_2=\frac{s_{11}+s_{01}}{s},\widehat{S}{p}_1=\frac{r_{01}+r_{00}}{r}\mathrm{and}\widehat{S}{p}_2=\frac{r_{10}+r_{00}}{r},$$

(9)

and applying the delta method, the estimators of the variances-covariances of $\widehat{S}{e}_i$ and $\widehat{S}{p}_i$ are

$$\begin{array}{c}\widehat{V}ar\left(\widehat{S}{e}_i\right)=\frac{\widehat{S}{e}_i\left(1-\widehat{S}{e}_i\right)}{s},Var\left(\widehat{S}{p}_i\right)=\frac{\widehat{S}{p}_i\left(1-\widehat{S}{p}_i\right)}{r},\\ \widehat{C}ov\left(\widehat{S}{e}_1,\widehat{S}{e}_2\right)=\frac{{\widehat{\epsilon }}_1}{s}and\widehat{C}ov\left(\widehat{S}{p}_1,\widehat{S}{p}_2\right)=\frac{{\widehat{\epsilon }}_0}{r},\end{array}$$

where ${\widehat{\epsilon }}_1=\frac{n{\widehat{p}}_{11}}{s}-\widehat{S}{e}_1\widehat{S}{e}_2=\frac{s_{11}{s}_{00}-s_{10}{s}_{01}}{s^2}$ and ${\widehat{\epsilon }}_0=\frac{n{\widehat{q}}_{00}}{r}-\widehat{S}{p}_1\widehat{S}{p}_2=\frac{r_{11}{r}_{00}-r_{10}{r}_{01}}{r^2}$, and the estimator of the disease prevalence is

$$\widehat{p}=\frac{s}{n}.$$

(10)

Let $Q_i=pS{e}_i+q\left(1-S{p}_i\right)$ be the probability that the result of the ith BDT is positive and let ${\overline{Q}}_i=1-Q_i=p\left(1-S{e}_i\right)+qS{p}_i$ be the probability that the result is negative. Its estimators are:

$${\widehat{Q}}_1=\frac{s_{10}+s_{11}+r_{10}+r_{11}}{n}\mathrm{and}{\widehat{Q}}_2=\frac{s_{01}+s_{11}+r_{01}+r_{11}}{n},$$

and

$${\widehat{\overline{Q}}}_1=\frac{s_{01}+s_{00}+r_{01}+r_{00}}{n}\mathrm{and}{\widehat{\overline{Q}}}_2=\frac{s_{10}+s_{00}+r_{10}+r_{00}}{n},$$

respectively. With respect to the predictive values, their estimators are:

$${\widehat{\tau }}_1=\frac{s_{11}+s_{10}}{s_{11}+s_{10}+r_{11}+r_{10}},{\widehat{\tau }}_2=\frac{s_{11}+s_{01}}{s_{11}+s_{01}+r_{11}+r_{01}},{\widehat{\upsilon }}_1=\frac{r_{01}+r_{00}}{s_{01}+s_{00}+r_{01}+r_{00}}$$

and

$${\widehat{\upsilon }}_2=\frac{r_{10}+r_{00}}{s_{10}+s_{00}+r_{10}+r_{00}}.$$

Applying the delta method, the estimators of the variances–covariances of ${\widehat{\tau }}_i$ and ${\widehat{\upsilon }}_i$ are [8]:

$$\begin{gathered} \widehat{V}ar\left({\widehat{\tau }}_1\right)=\frac{\left(s_{10}+s_{11}\right)\left(r_{10}+r_{11}\right)}{{\left(s_{10}+s_{11}+r_{10}+r_{11}\right)}^3},\widehat{V}ar\left({\widehat{\upsilon }}_1\right)=\frac{\left(s_{00}+s_{01}\right)\left(r_{00}+r_{01}\right)}{{\left(s_{00}+s_{01}+r_{00}+r_{01}\right)}^3} \\ \widehat{V}ar\left({\widehat{\tau }}_2\right)=\frac{\left(s_{01}+s_{11}\right)\left(r_{01}+r_{11}\right)}{{\left(s_{01}+s_{11}+r_{01}+r_{11}\right)}^3},\widehat{V}ar\left({\widehat{\upsilon }}_2\right)=\frac{\left(s_{00}+s_{10}\right)\left(r_{00}+r_{10}\right)}{{\left(s_{00}+s_{10}+r_{00}+r_{10}\right)}^3}, \\ \widehat{C}ov\left({\widehat{\tau }}_1,{\widehat{\tau }}_2\right)=\frac{s_{11}{r}_{10}{r}_{01}+r_{11}\left[s_{01}\left(s_{10}+s_{11}\right)+s_{11}\left(s_{11}+s_{10}+r_{11}+r_{10}+r_{01}\right)\right]}{{\left(s_{01}+s_{11}+r_{01}+r_{11}\right)}^2{\left(s_{10}+s_{11}+r_{10}+r_{11}\right)}^2} \end{gathered}$$

and

$$\widehat{C}ov\left({\widehat{\upsilon }}_1,{\widehat{\upsilon }}_2\right)=\frac{s_{10}\left(s_{00}+s_{01}\right)r_{00}+s_{00}\left[r_{00}^2+r_{01}{r}_{10}+r_{00}\left(s_{00}+s_{01}+r_{10}+r_{01}\right)\right]}{{\left(s_{00}+s_{01}+r_{00}+r_{01}\right)}^2{\left(s_{00}+s_{10}+r_{00}+r_{10}\right)}^2}.$$

When two parameters are compared in Statistics, the interest is to study the difference or the ratio between them. Then, we compare the positive (negative) predictive values of two BDTs through CIs for the difference, i.e., ${\delta }_{\tau }={\tau }_1-{\tau }_2$ and ${\delta }_{\upsilon }={\upsilon }_1-{\upsilon }_2$, and for the ratio, i.e., ${\rho }_{\tau }={\tau }_1/{\tau }_2$ and ${\rho }_{\upsilon }={\upsilon }_1/{\upsilon }_2$.

2.1. CIs for the Difference

Three CIs for each difference ${\delta }_{\tau }$ and ${\delta }_{\upsilon }$ are studied: Wald CI, bias-corrected bootstrap CI, and Monte Carlo Bayesian CI.

2.1.1. Wald CI

Wang et al. [4] have studied the comparison of the PVs of two BDTs through the weighted least square method. The test statistics for $H_0:{\tau }_1={\tau }_2$ and $H_0:{\upsilon }_1={\upsilon }_2$ are

$$z_{\tau }=\frac{{\widehat{\delta }}_{\tau }}{\sqrt{\widehat{V}ar\left({\widehat{\delta }}_{\tau }\right)}}and{z}_{\upsilon }=\frac{{\widehat{\delta }}_{\upsilon }}{\sqrt{\widehat{V}ar\left({\widehat{\delta }}_{\upsilon }\right)}},$$

respectively. Both test statistics follow a standard normal distribution where

$$\widehat{V}ar\left({\widehat{\delta }}_{\tau }\right)=\widehat{V}ar\left({\widehat{\tau }}_1\right)+\widehat{V}ar\left({\widehat{\tau }}_2\right)-2\widehat{C}ov\left({\widehat{\tau }}_1,{\widehat{\tau }}_2\right)$$

and

$$\widehat{V}ar\left({\widehat{\delta }}_{\upsilon }\right)=\widehat{V}ar\left({\widehat{\upsilon }}_1\right)+\widehat{V}ar\left({\widehat{\upsilon }}_2\right)-2\widehat{C}ov\left({\widehat{\upsilon }}_1,{\widehat{\upsilon }}_2\right)$$

are the estimators of the variances of ${\widehat{\delta }}_{\tau }$ and ${\widehat{\delta }}_{\upsilon }$, respectively. Inverting the two test statistics, the Wald CIs for ${\delta }_{\tau }$ and for ${\delta }_{\upsilon }$ are

$${\delta }_{\tau }\in {\widehat{\delta }}_{\tau }\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left({\widehat{\delta }}_{\tau }\right)}and{\delta }_{\upsilon }\in {\widehat{\delta }}_{\upsilon }\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left({\widehat{\delta }}_{\upsilon }\right)},$$

respectively, where $z_{1-\alpha /2}$ is the $100\left(1-\alpha /2\right)\mathrm{th}$ percentile of the standard normal distribution.

2.1.2. Bias-Corrected Bootstrap CI

The bias-corrected bootstrap CI is calculated from B random samples with replacement generated from the sample of n individuals. In each of the B samples, we calculate ${\widehat{\tau }}_{1b}$, ${\widehat{\tau }}_{2b}$, ${\widehat{\upsilon }}_{1b}$, ${\widehat{\upsilon }}_{2b}$, ${\widehat{\delta }}_{\tau b}={\widehat{\tau }}_{1b}-{\widehat{\tau }}_{2b}$ and ${\widehat{\delta }}_{\upsilon b}={\widehat{\upsilon }}_{1b}-{\widehat{\upsilon }}_{2b}$, with $b=1,\dots ,B$. Then, the average differences are calculated as ${\widehat{\overline{\delta }}}_{\tau B}=\frac{1}{B}{\displaystyle \sum _{b=1}^B{\widehat{\delta }}_{\tau b}}$ and ${\widehat{\overline{\delta }}}_{\upsilon B}=\frac{1}{B}{\displaystyle \sum _{b=1}^B{\widehat{\delta }}_{\tau b}}$. Assuming that the bootstrap statistics ${\widehat{\overline{\delta }}}_{\tau B}$ and ${\widehat{\overline{\delta }}}_{\upsilon B}$ can be transformed to a normal distribution, the bias-corrected bootstrap CIs [11] are calculated in the following way. Let $A_{\tau }=\#\left({\widehat{\delta }}_{\tau b}<{\widehat{\delta }}_{\tau }\right)$ be the number of bootstrap estimators ${\widehat{\delta }}_{\tau b}$ that are lower than the maximum likelihood estimator (MLE) ${\widehat{\delta }}_{\tau }$, and let $A_{\upsilon }=\#\left({\widehat{\delta }}_{\upsilon b}<{\widehat{\delta }}_{\upsilon }\right)$ be the number of bootstrap estimators ${\widehat{\delta }}_{\upsilon b}$ that are lower than the MLE ${\widehat{\delta }}_{\upsilon }$. Let ${\widehat{z}}_{\tau }={\mathsf{\Phi }}^{-1}\left(A_{\tau }/B\right)$ and ${\widehat{z}}_{\upsilon }={\mathsf{\Phi }}^{-1}\left(A_{\upsilon }/B\right)$, where ${\mathsf{\Phi }}^{-1}(\cdot )$ is the inverse function of the standard normal cumulative distribution function. Let ${\alpha }_{1\tau }=\mathsf{\Phi }\left(2{\widehat{z}}_{\tau }-z_{1-\alpha /2}\right)$, ${\alpha }_{2\tau }=\mathsf{\Phi }\left(2{\widehat{z}}_{\tau }+z_{1-\alpha /2}\right)$, ${\alpha }_{1\upsilon }=\mathsf{\Phi }\left(2{\widehat{z}}_{\upsilon }-z_{1-\alpha /2}\right)$, and ${\alpha }_{2\upsilon }=\mathsf{\Phi }\left(2{\widehat{z}}_{\upsilon }+z_{1-\alpha /2}\right)$; then, the bias-corrected bootstrap CI for ${\delta }_{\tau }$ is

$${\delta }_{\tau }\in \left({\widehat{\delta }}_{\tau B}^{\left({\alpha }_{1\tau }\right)},{\widehat{\delta }}_{\tau B}^{\left({\alpha }_{2\tau }\right)}\right)$$

and the bias-corrected bootstrap CI for ${\delta }_{\upsilon }$ is

$${\delta }_{\upsilon }\in \left({\widehat{\delta }}_{\upsilon B}^{\left({\alpha }_{1\upsilon }\right)},{\widehat{\delta }}_{\upsilon B}^{\left({\alpha }_{2\upsilon }\right)}\right)$$

where ${\widehat{\delta }}_{\tau B}^{\left(\gamma \right)}$ is the $\gamma \mathrm{th}$ quantile of the distribution of the B bootstrap estimations of ${\delta }_{\tau }$, and ${\widehat{\delta }}_{\upsilon B}^{\left(\gamma \right)}$ is the $\gamma \mathrm{th}$ quantile of the distribution of the B bootstrap estimations of ${\delta }_{\upsilon }$.

2.1.3. Monte Carlo Bayesian CI

The number of diseased individuals ($s$) is the product of a binomial distribution, i.e., $s\to B\left(n,p\right)$. Conditioning on $D=1$, it is verified that

$$s_{11}+s_{10}\to B\left(s,S{e}_1\right)and{s}_{11}+s_{01}\to B\left(s,S{e}_2\right).$$

The number of non-diseased individuals ($r$) is the product of a binomial distribution, i.e., $r\to B\left(n,q\right)$. Conditioning on $D=0$, it is verified that

$$r_{01}+r_{00}\to B\left(r,S{p}_1\right)and{r}_{10}+r_{00}\to B\left(r,S{p}_2\right).$$

On the other hand, the estimators $\widehat{S}{e}_i$, $\widehat{S}{p}_i$, and $\widehat{p}$ (Equations (9) and (10)) are estimators of binomial proportions. Therefore, for these estimators we propose conjugate beta prior distributions, i.e.,

$$\widehat{S}{e}_i\to Beta\left({\alpha }_{S{e}_i},{\beta }_{S{e}_i}\right),\widehat{S}{p}_i\to Beta\left({\alpha }_{S{p}_i},{\beta }_{S{p}_i}\right)and\widehat{p}\to Beta\left({\alpha }_p,{\beta }_p\right).$$

(11)

Let $n=\left(s_{11},s_{10},s_{01},s,r_{11},r_{10},r_{01},n-s\right)$ be the vector of observed frequencies, with $s_{00}=s-s_{11}-s_{10}-s_{01}$, $r=n-s$, and $r_{00}=n-s-r_{11}-r_{10}-r_{01}$. Then, the posteriori distributions for the estimators of $S{e}_i$, $S{p}_i$, and p are:

$$\begin{array}{c}\widehat{S}{e}_1|n\to Beta\left(s_{11}+s_{10}+{\alpha }_{S{e}_1},s-s_{11}-s_{10}+{\beta }_{S{e}_1}\right),\\ \widehat{S}{p}_1|n\to Beta\left(r_{01}+r_{00}+{\alpha }_{S{p}_1},n-s-r_{01}-r_{00}+{\beta }_{S{p}_1}\right),\\ \widehat{S}{e}_2|n\to Beta\left(s_{11}+s_{01}+{\alpha }_{S{e}_2},s-s_{11}-s_{01}+{\beta }_{S{e}_2}\right),\\ \widehat{S}{p}_2|n\to Beta\left(r_{10}+r_{00}+{\alpha }_{S{p}_2},n-s-r_{10}-r_{00}+{\beta }_{S{p}_2}\right),\\ \widehat{p}|n\to Beta\left(s+{\alpha }_p,n-s+{\beta }_p\right).\end{array}$$

(12)

The posteriori distribution for the positive (negative) predictive value of each BDT, and for ${\delta }_{\tau }$ and ${\delta }_{\upsilon }$, can be approximated applying the Monte Carlo method. The Monte Carlo method is a computational method that consists of generating M values from the posteriori distributions (12). In the mth iteration, the values generated for $S{e}_j^{\left(m\right)}$, $S{p}_j^{\left(m\right)}$, and $p^{\left(m\right)}$ are plugged in the equations

$${\tau }_i^{\left(m\right)}=\frac{p^{\left(m\right)}\times S{e}_i^{\left(m\right)}}{p^{\left(m\right)}\times S{e}_i^{\left(m\right)}+\left(1-p^{\left(m\right)}\right)\times \left(1-S{p}_i^{\left(m\right)}\right)}$$

and

$${\upsilon }_i^{\left(m\right)}=\frac{\left(1-p^{\left(m\right)}\right)\times S{p}_i^{\left(m\right)}}{p^{\left(m\right)}\times \left(1-S{e}_i^{\left(m\right)}\right)+\left(1-p^{\left(m\right)}\right)\times S{p}_i^{\left(m\right)}},$$

with $i=1,2$, and then, ${\delta }_{\tau }^{\left(m\right)}={\tau }_1^{\left(m\right)}-{\tau }_2^{\left(m\right)}$ and ${\delta }_{\upsilon }^{\left(m\right)}={\upsilon }_1^{\left(m\right)}-{\upsilon }_2^{\left(m\right)}$ are calculated. As estimators of ${\delta }_{\tau }$ and ${\delta }_{\upsilon }$, we calculate the average of the M estimations of the differences, i.e., ${\widehat{\overline{\delta }}}_{\tau Bay}=\frac{1}{M}{\displaystyle \sum _{m=1}^M{\delta }_{\tau }^{\left(m\right)}}$ and ${\widehat{\overline{\delta }}}_{\upsilon Bay}=\frac{1}{M}{\displaystyle \sum _{m=1}^M{\delta }_{\upsilon }^{\left(m\right)}}$. Finally, based on the M values of ${\delta }_{\tau }^{\left(m\right)}$ and of ${\delta }_{\upsilon }^{\left(m\right)}$, we propose CIs based on quantiles. Therefore, the $100\times \left(1-\alpha \right)\%$ CI for ${\delta }_{\tau }$ is

$${\delta }_{\tau }\in \left(q_{\tau Bay}^{\left(\alpha /2\right)},q_{\tau Bay}^{\left(1-\alpha /2\right)}\right)$$

and the $100\times \left(1-\alpha \right)\%$ CI for ${\delta }_{\upsilon }$ is

$${\delta }_{\upsilon }\in \left(q_{\upsilon Bay}^{\left(\alpha /2\right)},q_{\upsilon Bay}^{\left(1-\alpha /2\right)}\right)$$

where $q_{\tau Bay}^{\left(\gamma \right)}$ $\left(q_{\upsilon Bay}^{\left(\gamma \right)}\right)$ is the $\gamma \mathrm{th}$ quantile of the distribution of the M values ${\delta }_{\tau }^{\left(m\right)}$ $\left({\delta }_{\upsilon }^{\left(m\right)}\right)$.

2.2. CIs for the Ratio

Five CIs for each ratio ${\rho }_{\tau }={\tau }_1/{\tau }_2$ and ${\rho }_{\upsilon }={\upsilon }_1/{\upsilon }_2$ are studied: Wald CI, logarithmic CI, Fieller CI, bias-corrected bootstrap CI, and Monte Carlo Bayesian CI.

2.2.1. Wald CI

Moskowitz and Pepe [9] have studied a Wald-type confidence interval for the ratio of the two positive (negative) predictive values. The $100\times \left(1-\alpha \right)\%$ Wald CI for ${\delta }_{\tau }$ is

$${\rho }_{\tau }\in {\widehat{\rho }}_{\tau }\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left({\widehat{\rho }}_{\tau }\right)}$$

and the $100\times \left(1-\alpha \right)\%$ Wald CI for ${\delta }_{\upsilon }$ is

$${\rho }_{\upsilon }\in {\widehat{\rho }}_{\upsilon }\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left({\widehat{\rho }}_{\upsilon }\right)},$$

where $\widehat{V}ar\left({\widehat{\rho }}_{\tau }\right)$ and $\widehat{V}ar\left({\widehat{\rho }}_{\upsilon }\right)$, obtained applying the delta method, are

$$\widehat{V}ar\left({\widehat{\rho }}_{\tau }\right)=\frac{{\widehat{\tau }}_2^2\widehat{V}ar\left({\widehat{\tau }}_1\right)+{\widehat{\tau }}_1^2\widehat{V}ar\left({\widehat{\tau }}_2\right)-2{\widehat{\tau }}_1{\widehat{\tau }}_2\widehat{C}ov\left({\widehat{\tau }}_1,{\widehat{\tau }}_2\right)}{{\widehat{\tau }}_2^4}$$

and

$$\widehat{V}ar\left({\widehat{\rho }}_{\upsilon }\right)=\frac{{\widehat{\upsilon }}_2^2\widehat{V}ar\left({\widehat{\upsilon }}_1\right)+{\widehat{\upsilon }}_1^2\widehat{V}ar\left({\widehat{\upsilon }}_2\right)-2{\widehat{\upsilon }}_1{\widehat{\upsilon }}_2\widehat{C}ov\left({\widehat{\upsilon }}_1,{\widehat{\upsilon }}_2\right)}{{\widehat{\upsilon }}_2^4}.$$

These CIs are for ${\rho }_{\tau }={\tau }_1/{\tau }_2$ and ${\rho }_{\upsilon }={\upsilon }_1/{\upsilon }_2$. If we want to calculate the CI for the ratio ${\tau }_2/{\tau }_1=1/{\rho }_{\tau }$ and for the ratio ${\upsilon }_2/{\upsilon }_1=1/{\rho }_{\upsilon }$, then we have to divide the CI for ${\rho }_{\tau }$ by ${\widehat{\rho }}_{\tau }^2$ and the CI for ${\rho }_{\upsilon }$ by ${\widehat{\rho }}_{\upsilon }^2$. For example, if $\left(L_{\tau },U_{\tau }\right)$ is the Wald CI for ${\tau }_1/{\tau }_2$, then $\left(L_{\tau }/{\widehat{\rho }}_{\tau }^2,U_{\tau }/{\widehat{\rho }}_{\tau }^2\right)$ is the Wald CI for ${\tau }_2/{\tau }_1$.

2.2.2. Logarithmic CI

Assuming the asymptotic normality of the Napierian logarithm of ${\widehat{\rho }}_{\tau }$ and of ${\widehat{\rho }}_{\upsilon }$, i.e., $ln\left({\widehat{\rho }}_{\tau }\right)\to N\left(ln\left({\delta }_{\tau }\right),Var\left[ln\left({\delta }_{\tau }\right)\right]\right)$ and $ln\left({\widehat{\rho }}_{\upsilon }\right)\to N\left(ln\left({\delta }_{\upsilon }\right),Var\left[ln\left({\delta }_{\upsilon }\right)\right]\right)$ when n is large, an asymptotic CI for $ln\left({\rho }_{\tau }\right)$ is

$$ln\left({\widehat{\rho }}_{\tau }\right)\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\tau }\right)\right]}$$

and an asymptotic CI for $ln\left({\rho }_{\upsilon }\right)$ is

$$ln\left({\widehat{\rho }}_{\upsilon }\right)\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\upsilon }\right)\right]}.$$

Taking exponential in each of the previous expressions, the logarithmic CI for ${\delta }_{\tau }$ is

$${\rho }_{\tau }\in {\widehat{\rho }}_{\tau }\times exp\left\{\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\tau }\right)\right]}\right\}$$

and the logarithmic CI for ${\delta }_{\upsilon }$ is

$${\rho }_{\upsilon }\in {\widehat{\rho }}_{\upsilon }\times exp\left\{\pm z_{1-\alpha /2}\sqrt{\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\upsilon }\right)\right]}\right\},$$

where $\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\tau }\right)\right]$ and $\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\upsilon }\right)\right]$, obtained applying the delta method, are:

$$\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\tau }\right)\right]=\frac{\widehat{V}ar\left[{\widehat{\tau }}_1\right]}{{\widehat{\tau }}_1^2}+\frac{\widehat{V}ar\left[{\widehat{\tau }}_2\right]}{{\widehat{\tau }}_2^2}-\frac{2\widehat{C}ov\left[{\widehat{\tau }}_1,{\widehat{\tau }}_2\right]}{{\widehat{\tau }}_1{\widehat{\tau }}_2}$$

and

$$\widehat{V}ar\left[ln\left({\widehat{\rho }}_{\upsilon }\right)\right]=\frac{\widehat{V}ar\left[{\widehat{\upsilon }}_1\right]}{{\widehat{\upsilon }}_1^2}+\frac{\widehat{V}ar\left[{\widehat{\upsilon }}_2\right]}{{\widehat{\upsilon }}_2^2}-\frac{2\widehat{C}ov\left[{\widehat{\upsilon }}_1,{\widehat{\upsilon }}_2\right]}{{\widehat{\upsilon }}_1{\widehat{\upsilon }}_2}.$$

If we want to calculate the logarithmic CI for the ratio ${\tau }_2/{\tau }_1$, then the CI is obtained by calculating the inverse of each boundary of CI for ${\rho }_{\tau }={\tau }_1/{\tau }_2$. In a similar way, the CI for ${\upsilon }_2/{\upsilon }_1$ is calculated.

2.2.3. Fieller CI

The method of Fieller [12] is a classic method used to estimate the ratio of two parameters. In order to apply this method, it is necessary to assume that ${\widehat{\tau }}^T\to N\left(\tau ,{\sum }_{\tau }\right)$ and that ${\widehat{\upsilon }}^T\to N\left(\upsilon ,{\sum }_{\upsilon }\right)$; i.e., it is necessary to assume that the estimators of the positive (negative) predictive values are distributed according to a normal bivariate distribution, and where $\widehat{\tau }=\left({\widehat{\tau }}_1,{\widehat{\tau }}_2\right)$, $\widehat{\upsilon }=\left({\widehat{\upsilon }}_1,{\widehat{\upsilon }}_2\right)$,

$${\sum }_{\tau }=\left(\begin{array}{cc}Var\left({\tau }_1\right)& Cov\left({\tau }_1,{\tau }_2\right)\\ Cov\left({\tau }_1,{\tau }_2\right)& Var\left({\tau }_2\right)\end{array}\right)=\left(\begin{array}{cc}{\sigma }_{\tau 11}& {\sigma }_{\tau 12}\\ {\sigma }_{\tau 21}& {\sigma }_{\tau 22}\end{array}\right)$$

and

$${\sum }_{\upsilon }=\left(\begin{array}{cc}Var\left({\upsilon }_1\right)& Cov\left({\upsilon }_1,{\upsilon }_2\right)\\ Cov\left({\upsilon }_1,{\upsilon }_2\right)& Var\left({\upsilon }_2\right)\end{array}\right)=\left(\begin{array}{cc}{\sigma }_{\upsilon 11}& {\sigma }_{\upsilon 12}\\ {\sigma }_{\upsilon 21}& {\sigma }_{\upsilon 22}\end{array}\right).$$

Applying the method of Fieller, it is verified that ${\widehat{\tau }}_1-{\rho }_{\tau }{\widehat{\tau }}_2\to N\left(0,{\sigma }_{\tau 11}+{\rho }_{\tau }^2{\sigma }_{\tau 22}-2{\rho }_{\tau }{\sigma }_{\tau 12}\right)$ and that ${\widehat{\upsilon }}_1-{\rho }_{\upsilon }{\widehat{\upsilon }}_2\to N\left(0,{\sigma }_{\upsilon 11}+{\rho }_{\upsilon }^2{\sigma }_{\upsilon 22}-2{\rho }_{\upsilon }{\sigma }_{\upsilon 12}\right)$ when n is large. The Fieller CI for ${\rho }_{\tau }$ is obtained by solving the inequality

$$\frac{{\left({\widehat{\tau }}_1-{\rho }_{\tau }{\widehat{\tau }}_2\right)}^2}{{\widehat{\sigma }}_{\tau 11}+{\rho }_{\tau }^2{\widehat{\sigma }}_{\tau 22}-2{\rho }_{\tau }{\widehat{\sigma }}_{\tau 12}}<z_{1-\alpha /2}^2,$$

and the Fieller CI for ${\rho }_{\upsilon }$ is obtained by solving the inequality

$$\frac{{\left({\widehat{\upsilon }}_1-{\rho }_{\upsilon }{\widehat{\upsilon }}_2\right)}^2}{{\widehat{\sigma }}_{\upsilon 11}+{\rho }_{\upsilon }^2{\widehat{\sigma }}_{\upsilon 22}-2{\rho }_{\upsilon }{\widehat{\sigma }}_{\upsilon 12}}<z_{1-\alpha /2}^2.$$

Finally, the Fieller CI for ${\rho }_{\tau }$ is

$${\rho }_{\tau }\in \frac{{\widehat{\beta }}_{\tau 12}\pm \sqrt{{\widehat{\beta }}_{\tau 12}^2-{\widehat{\beta }}_{\tau 11}{\widehat{\beta }}_{\tau 22}}}{{\widehat{\beta }}_{\tau 22}}$$

where ${\widehat{\beta }}_{\tau ij}={\widehat{\tau }}_i{\widehat{\tau }}_j-{\widehat{\sigma }}_{\tau ij}{z}_{1-\alpha /2}^2$ with $i,j=1,2$ and verifying that ${\widehat{\beta }}_{\tau 12}={\widehat{\beta }}_{\tau 21}$. This CI is valid when ${\widehat{\beta }}_{\tau 12}^2>{\widehat{\beta }}_{\tau 11}{\widehat{\beta }}_{\tau 22}$ and ${\widehat{\beta }}_{\tau 22}\ne 0$. Similarly, the Fieller CI for ${\rho }_{\upsilon }$ is

$${\rho }_{\upsilon }\in \frac{{\widehat{\beta }}_{\upsilon 12}\pm \sqrt{{\widehat{\beta }}_{\upsilon 12}^2-{\widehat{\beta }}_{\upsilon 11}{\widehat{\beta }}_{\upsilon 22}}}{{\widehat{\beta }}_{\upsilon 22}}$$

where ${\widehat{\beta }}_{\upsilon ij}={\widehat{\upsilon }}_i{\widehat{\upsilon }}_j-{\widehat{\sigma }}_{\upsilon ij}{z}_{1-\alpha /2}^2$ with $i,j=1,2$, and ${\widehat{\beta }}_{\upsilon 12}={\widehat{\beta }}_{\upsilon 21}$. This CI is valid when ${\widehat{\beta }}_{\upsilon 12}^2>{\widehat{\beta }}_{\upsilon 11}{\widehat{\beta }}_{\upsilon 22}$ and ${\widehat{\beta }}_{\upsilon 22}\ne 0$. The Fieller CI for ${\tau }_2/{\tau }_1$ $\left({\upsilon }_2/{\upsilon }_1\right)$ is calculated by inverting the limits of the CI for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$.

2.2.4. Bias-Corrected Bootstrap CI

The bias-corrected bootstrap CI for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$ is obtained in a similar way to that of ${\delta }_{\tau }$ $\left({\delta }_{\upsilon }\right)$. In each sample with a replacement obtained, we calculate ${\widehat{\tau }}_{1b}$, ${\widehat{\tau }}_{2b}$, ${\widehat{\upsilon }}_{1b}$, ${\widehat{\upsilon }}_{2b}$, ${\widehat{\rho }}_{\tau b}={\widehat{\tau }}_{1b}/{\widehat{\tau }}_{2b}$, and ${\widehat{\rho }}_{\upsilon b}={\widehat{\upsilon }}_{1b}/{\widehat{\upsilon }}_{2b}$, with $b=1,\dots ,B$. Then, based on the B ratios, we estimate the average ratios as ${\widehat{\overline{\rho }}}_{\tau B}=\frac{1}{B}{\displaystyle \sum _{b=1}^B{\widehat{\rho }}_{\tau b}}$ and ${\widehat{\overline{\rho }}}_{\upsilon B}=\frac{1}{B}{\displaystyle \sum _{b=1}^B{\widehat{\rho }}_{\upsilon b}}$. Assuming that these statistics can be transformed to a normal distribution, the bias-corrected bootstrap CI [11] for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$ is calculated in a similar way as the bias-corrected bootstrap CI for ${\delta }_{\tau }$ $\left({\delta }_{\upsilon }\right)$, considering that $A_{\tau }=\#\left({\widehat{\rho }}_{\tau b}<{\widehat{\rho }}_{\tau }\right)$ and that $A_{\upsilon }=\#\left({\widehat{\rho }}_{\upsilon b}<{\widehat{\rho }}_{\upsilon }\right)$. Finally, the bias-corrected bootstrap CI for ${\rho }_{\tau }$ is

$${\rho }_{\tau }\in \left({\widehat{\rho }}_{\tau B}^{\left({\alpha }_1\right)},{\widehat{\rho }}_{\tau B}^{\left({\alpha }_2\right)}\right),$$

where ${\widehat{\rho }}_{\tau B}^{\left(\gamma \right)}$ is the $\gamma \mathrm{th}$ quantile of the distribution of the B bootstrap estimations of ${\rho }_{\tau }$. Similarly, the bias-corrected bootstrap CI for ${\rho }_{\upsilon }$ is

$${\rho }_{\upsilon }\in \left({\widehat{\rho }}_{\upsilon B}^{\left({\alpha }_1\right)},{\widehat{\rho }}_{\upsilon B}^{\left({\alpha }_2\right)}\right),$$

where ${\widehat{\rho }}_{\upsilon B}^{\left(\gamma \right)}$ is the $\gamma \mathrm{th}$ quantile of the distribution of the B bootstrap estimations of ${\rho }_{\upsilon }$. The bias-corrected bootstrap CI for ${\tau }_2/{\tau }_1$ $\left({\upsilon }_2/{\upsilon }_1\right)$ is calculated by inverting the limits of the bias-corrected bootstrap CI for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$.

2.2.5. Monte Carlo Bayesian CI

The Monte Carlo Bayesian CI for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$ is obtained in a similar way as the Monte Carlo Bayesian CI for ${\delta }_{\tau }$ $\left({\delta }_{\upsilon }\right)$. Considering the same distributions (10) and (11), in the mth iteration, we calculate the ratios ${\rho }_{\tau }^{\left(m\right)}={\tau }_1^{\left(m\right)}/{\tau }_2^{\left(m\right)}$ and ${\rho }_{\upsilon }^{\left(m\right)}={\upsilon }_1^{\left(m\right)}/{\upsilon }_2^{\left(m\right)}$. As estimators, we calculate ${\widehat{\overline{\rho }}}_{\tau Bay}=\frac{1}{M}{\displaystyle \sum _{m=1}^M{\rho }_{\tau }^{\left(m\right)}}$ and ${\widehat{\overline{\rho }}}_{\upsilon Bay}=\frac{1}{M}{\displaystyle \sum _{m=1}^M{\rho }_{\upsilon }^{\left(m\right)}}$. Finally, we calculate the CIs based on quantiles, i.e.,

$${\rho }_{\tau }\in \left(q_{\tau Bay}^{\left(\alpha /2\right)},q_{\tau Bay}^{\left(1-\alpha /2\right)}\right)and{\rho }_{\upsilon }\in \left(q_{\upsilon Bay}^{\left(\alpha /2\right)},q_{\upsilon Bay}^{\left(1-\alpha /2\right)}\right)$$

where $q_{\tau Bay}^{\left(\gamma \right)}$ $\left(q_{\upsilon Bay}^{\left(\gamma \right)}\right)$ is the $\gamma \mathrm{th}$ quantile of the distribution of the M values ${\rho }_{\tau }^{\left(m\right)}$ $\left({\rho }_{\upsilon }^{\left(m\right)}\right)$. The Monte Carlo Bayesian CI for ${\tau }_2/{\tau }_1$ $\left({\upsilon }_2/{\upsilon }_1\right)$ is calculated by inverting the limits of the Monte Carlo Bayesian CI for ${\rho }_{\tau }$ $\left({\rho }_{\upsilon }\right)$.

3. Simulation Experiments

The CIs studied in Section 2 are approximate, and therefore, it is necessary to study their asymptotic behaviors. For this, Monte Carlo simulation experiments have been carried out to study the coverage probabilities and the average lengths of the CIs studied, considering a confidence level of 95%. These experiments have consisted of generating $N=10,000$ random samples with multinomial distribution sized $n=\left\{50,100,200,500,1000\right\}$, and whose probabilities have been calculated from Equations (7) and (8). The experiments have been designed from the predictive values of both BDTs. As the value of disease prevalence, we have taken $p=\left\{10\%,25\%,50\%,75\%\right\}$, and as values of predictive values, we have taken the values ${\tau }_i,{\upsilon }_i=\left\{0.70,0.75,\dots ,0.90,0.95\right\}$, which are realistic values in clinical practice. Next, using these values, we have calculated the maximum values of the dependence factors ${\epsilon }_1$ and ${\epsilon }_0$ (Equations (5) and (6)). As values of ${\epsilon }_1$ and ${\epsilon }_0$, we have taken intermediate and high values, i.e., 50% of the maximum value of ${\epsilon }_i$ and 90% of the maximum value of ${\epsilon }_i$, respectively. Finally, we have calculated the probabilities of the multinomial distributions using Equations (7) and (8). In each scenario, we have calculated all the CIs for each of the N random samples.

For the bias-corrected bootstrap CIs, for each one of the N random samples, $B=2000$ samples with replacement have been generated, and from these B samples, the bias-corrected bootstrap CIs have been calculated.

For the Monte Carlo Bayesian CI, we have considered a $Beta\left(1,1\right)$ distribution as a priori distribution for the estimators of sensitivities, specificities and prevalence. The $Beta\left(1,1\right)$ distribution is a non-informative distribution, which is flat for every possible value of the sensitivities, specificities, and prevalence. Therefore, the impact of the $Beta\left(1,1\right)$ distribution on the posteriori distributions is minimal. Moreover, for each one of the N random samples, $M=10,000$ random samples have been generated, and the Monte Carlo Bayesian CIs have been calculated from them.

In each of the N samples generated, all the CIs have been calculated. Furthermore, it has been checked whether each CI contains the value of the parameter (difference or ratio, depending on the type of CI). The coverage probability has been calculated by dividing the number of samples in which the CI contains the parameter by the total number of samples. For each CI, its length (upper limit minus the lower limit) has also been calculated, and finally, the average length of each CI has been calculated.

3.1. CIs for the Differences and Ratios of Positive Predictive Values

Table 2. Asymptotic behaviors of the CIs for the difference of the two positive predictive values.

$\begin{array}{c}{\mathit{\tau }}_1=0.75 {\mathit{\tau }}_2=0.75 {\mathit{\upsilon }}_1=0.95 {\mathit{\upsilon }}_2=0.95 {\mathit{\delta }}_{\mathit{\tau }}=0\\ {\mathit{\epsilon }}_1=0.124 {\mathit{\epsilon }}_0=0.010 \mathit{p}=0.10\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	1	0.749	1	0.743	1	0.787
100	1	0.538	1	0.527	1	0.683
200	0.988	0.409	1	0.397	1	0.454
500	0.953	0.261	0.990	0.254	0.998	0.363
1000	0.944	0.185	0.957	0.186	0.993	0.259
$\begin{array}{c}{\tau }_1=0.90 {\tau }_2=0.85 {\upsilon }_1=0.95 {\upsilon }_2=0.90 {\delta }_{\tau }=0.05\\ {\epsilon }_1=0.021 {\epsilon }_0=0.044 \mathit{p}=0.50\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	0.982	0.251	1	0.242	0.999	0.326
100	0.951	0.174	0.966	0.182	0.993	0.223
200	0.954	0.122	0.952	0.126	0.991	0.156
500	0.941	0.077	0.933	0.077	0.981	0.099
1000	0.952	0.055	0.948	0.055	0.988	0.070
$\begin{array}{c}{\tau }_1=0.85 {\tau }_2=0.75 {\upsilon }_1=0.95 {\upsilon }_2=0.85 {\delta }_{\tau }=0.10\\ {\epsilon }_1=0.037 {\epsilon }_0=0.024 \mathit{p}=0.25\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	0.998	0.513	1	0.499	1	0.602
100	0.981	0.354	1	0.343	0.999	0.445
200	0.941	0.250	0.987	0.243	0.991	0.318
500	0.956	0.158	0.959	0.159	0.989	0.204
1000	0.953	0.112	0.954	0.113	0.989	0.145

CP: coverage probability. AL: average length. Wald: Wald CI. BCB: bias-corrected bootstrap CI. MCB: Monte Carlo Bayesian CI.

shows some of the results obtained for the three CIs for the difference ${\delta }_{\tau }$ for four different scenarios and for intermediate values of ${\epsilon }_1$ and ${\epsilon }_0$. When the sample size is small $\left(n=50\right)$ or moderate $\left(n=100\right)$, the CIs for ${\delta }_{\tau }$ have a coverage probability close to 1. For the difference ${\delta }_{\tau }$, in very general terms, the Wald CI is the interval that has a coverage probability with better fluctuations around 95%, especially when n is moderate or large $\left(n\ge 200\right)$. The bias-corrected bootstrap CI has a very similar behavior to the Wald CI, especially when the sample size is large. In general terms, the Monte Carlo Bayesian CI has a coverage probability greater than that of the other two intervals, even when the coverage probability of the other two intervals fluctuates around 95%.

Regarding the CIs for the ratio ${\rho }_{\tau }$,

Table 3. Asymptotic behaviors of the CIs for the ratio of the two positive predictive values.

$\begin{array}{c}{\mathit{\tau }}_1=0.75 {\mathit{\tau }}_2=0.75 {\mathit{\upsilon }}_1=0.95 {\mathit{\upsilon }}_2=0.95 {\mathit{\rho }}_{\mathit{\tau }}=1\\ {\mathit{\epsilon }}_1=0.124 {\mathit{\epsilon }}_0=0.010 \mathit{p}=0.10\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	1	1.326	1	1.450	1	2.046	1	1.348	1	2.183
100	0.999	0.966	1	1.018	1	1.311	1	0.973	1	1.534
200	0.989	0.630	0.992	0.643	0.994	0.652	1	0.648	1	0.978
500	0.962	0.359	0.966	0.361	0.953	0.369	0.988	0.364	0.998	0.535
1000	0.956	0.250	0.952	0.251	0.944	0.253	0.956	0.256	0.993	0.363
$\begin{array}{c}{\tau }_1=0.90 {\tau }_2=0.85 {\upsilon }_1=0.95 {\upsilon }_2=0.90 {\rho }_{\tau }=1.06 \\ {\epsilon }_1=0.021 {\epsilon }_0=0.044 \mathit{p}=0.50\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	0.986	0.326	0.994	0.328	0.993	0.334	1	0.347	1	0.448
100	0.949	0.216	0.952	0.216	0.950	0.218	0.957	0.232	0.994	0.288
200	0.952	0.151	0.955	0.151	0.954	0.152	0.953	0.157	0.990	0.196
500	0.941	0.095	0.941	0.095	0.940	0.095	0.935	0.095	0.982	0.122
1000	0.950	0.067	0.951	0.067	0.949	0.067	0.946	0.067	0.987	0.086
$\begin{array}{c}{\tau }_1=0.85 {\tau }_2=0.75 {\upsilon }_1=0.95 {\upsilon }_2=0.85 {\rho }_{\tau }=1.13\\ {\epsilon }_1=0.037 {\epsilon }_0=0.024 \mathit{p}=0.25\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	0.997	0.950	1	0.958	0.998	0.961	1	0.992	1	1.407
100	0.972	0.591	0.983	0.596	0.979	0.636	1	0.689	0.999	0.841
200	0.941	0.390	0.943	0.392	0.940	0.396	0.988	0.398	0.989	0.528
500	0.950	0.241	0.954	0.241	0.957	0.242	0.960	0.244	0.989	0.314
1000	0.951	0.169	0.953	0.170	0.950	0.171	0.953	0.171	0.988	0.218

CP: coverage probability. AL: average length. Wald: Wald CI. Log: logarithmic CI. Fieller: Fieller CI. BCB: bias-corrected bootstrap CI. MCB: Monte Carlo Bayesian CI.

shows the results obtained for the same scenarios as in Table 2. When the sample size is small $\left(n=50\right)$, the five CIs for ${\rho }_{\tau }$ have a coverage probability close to 1. In general terms, there is not an important difference between the coverage probabilities and the average lengths of the Wald, logarithm, and Fieller CIs, especially when $n\ge 100$. When the sample size is small, the logarithmic CI and the Fieller CI have an average length slightly greater than the Wald CI. The bias-corrected bootstrap CI has a very similar behavior to the Wald, logarithmic, and Fieller CIs, especially when the sample size is large. In general terms, the Monte Carlo Bayesian CI has a coverage probability greater than that of the other four intervals.

3.2. CIs for the Differences and Ratios of Negative Predictive Values

Table 4. Asymptotic behaviors of the CIs for the difference of the two negative predictive values.

$\begin{array}{c}{\mathit{\tau }}_1=0.75 {\mathit{\tau }}_2=0.75 {\mathit{\upsilon }}_1=0.95 {\mathit{\upsilon }}_2=0.95 {\mathit{\delta }}_{\mathit{\upsilon }}=0\\ {\mathit{\epsilon }}_1=0.124 {\mathit{\epsilon }}_0=0.010 \mathit{p}=0.10\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	1	0.127	1	0.119	1	0.154
100	0.999	0.072	1	0.069	0.999	0.095
200	0.989	0.046	1	0.042	0.999	0.063
500	0.949	0.028	0.968	0.028	0.994	0.040
1000	0.946	0.020	0.943	0.020	0.993	0.028
$\begin{array}{c}{\tau }_1=0.90 {\tau }_2=0.85 {\upsilon }_1=0.95 {\upsilon }_2=0.90 {\delta }_{\upsilon }=0.05\\ {\epsilon }_1=0.021 {\epsilon }_0=0.044 \mathit{p}=0.50\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	0.999	0.264	1	0.276	1	0.344
100	0.966	0.169	0.952	0.170	0.995	0.222
200	0.950	0.115	0.932	0.119	0.989	0.147
500	0.949	0.073	0.948	0.073	0.983	0.090
1000	0.952	0.052	0.953	0.052	0.984	0.063
$\begin{array}{c}{\tau }_1=0.85 {\tau }_2=0.75 {\upsilon }_1=0.95 {\upsilon }_2=0.85 {\delta }_{\upsilon }=0.10\\ {\epsilon }_1=0.037 {\epsilon }_0=0.024 \mathit{p}=0.25\end{array}$
	Wald		BCB		MCB
n	CP	AL	CP	AL	CP	AL
50	0.936	0.207	0.720	0.182	0.948	0.218
100	0.938	0.142	0.874	0.133	0.960	0.151
200	0.948	0.099	0.937	0.096	0.967	0.107
500	0.957	0.062	0.961	0.062	0.975	0.068
1000	0.946	0.044	0.947	0.044	0.964	0.048

CP: coverage probability. AL: average length. Wald: Wald CI. BCB: bias-corrected bootstrap CI. MCB: Monte Carlo Bayesian CI.

shows the results for the three CIs for the difference ${\delta }_{\upsilon }$ for the same scenarios as in Table 2 and Table 3. In general terms, when the sample size is small, the three CIs have a coverage probability close to 1, although in some situations, the bias-corrected bootstrap CI may have a coverage probability well below 95%. In general terms, the Monte Carlo Bayesian CI has a coverage probability that almost always fluctuates above 95%. For the difference ${\delta }_{\upsilon }$, the Wald CI is the interval that has a coverage probability with better fluctuations around 95%, especially when the sample size is moderate or large.

Regarding the CIs for the ratio ${\rho }_{\upsilon }$,

Table 5. Asymptotic behaviors of the CIs for the ratio of the two negative predictive values.

$\begin{array}{c}{\mathit{\tau }}_1=0.75 {\mathit{\tau }}_2=0.75 {\mathit{\upsilon }}_1=0.95 {\mathit{\upsilon }}_2=0.95 {\mathit{\rho }}_{\mathit{\upsilon }}=1\\ {\mathit{\epsilon }}_1=0.124 {\mathit{\epsilon }}_0=0.010 \mathit{p}=0.10\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	1	0.144	1	0.144	1	0.145	1	0.128	1.000	0.173
100	0.999	0.076	1	0.079	0.999	0.080	1	0.074	0.999	0.103
200	0.988	0.046	0.991	0.047	0.992	0.048	1	0.045	0.999	0.068
500	0.950	0.030	0.950	0.030	0.949	0.030	0.971	0.029	0.994	0.042
1000	0.948	0.021	0.947	0.021	0.946	0.021	0.946	0.021	0.993	0.030
$\begin{array}{c}{\tau }_1=0.90 {\tau }_2=0.85 {\upsilon }_1=0.95 {\upsilon }_2=0.90 {\rho }_{\upsilon }=1.06\\ {\epsilon }_1=0.021 {\epsilon }_0=0.044 \mathit{p}=0.50\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	1	0.324	1	0.326	0.999	0.332	1	0.349	1	0.462
100	0.954	0.201	0.964	0.201	0.962	0.202	0.964	0.219	0.995	0.274
200	0.945	0.134	0.946	0.134	0.947	0.135	0.921	0.138	0.989	0.175
500	0.946	0.085	0.945	0.085	0.946	0.085	0.947	0.085	0.982	0.105
1000	0.954	0.060	0.952	0.060	0.952	0.060	0.950	0.060	0.983	0.074
$\begin{array}{c}{\tau }_1=0.85 {\tau }_2=0.75 {\upsilon }_1=0.95 {\upsilon }_2=0.85 {\rho }_{\upsilon }=1.12\\ {\epsilon }_1=0.037 {\epsilon }_0=0.024 \mathit{p}=0.25\end{array}$
	Wald		Log		Fieller		BCB		MCB
n	CP	AL	CP	AL	CP	AL	CP	AL	CP	AL
50	0.936	0.271	0.934	0.272	0.933	0.275	0.724	0.283	0.933	0.291
100	0.939	0.184	0.936	0.184	0.935	0.185	0.849	0.191	0.961	0.199
200	0.945	0.129	0.947	0.129	0.945	0.129	0.923	0.125	0.963	0.140
500	0.957	0.081	0.958	0.081	0.958	0.081	0.959	0.081	0.974	0.088
1000	0.949	0.057	0.950	0.057	0.950	0.057	0.948	0.057	0.964	0.062

CP: coverage probability. AL: average length. Wald: Wald CI. Log: logarithmic CI. Fieller: Fieller CI. BCB: bias-corrected bootstrap CI. MCB: Monte Carlo Bayesian CI.

shows the results obtained for the same scenarios as in Table 4. When the sample size is small $\left(n=50\right)$, the five CIs for ${\rho }_{\tau }$ fail or have a coverage probability close to 1. In general terms, the conclusions for the bias-corrected bootstrap CI and for the Monte Carlo Bayesian CI are very similar to those obtained for the corresponding intervals for the difference ${\delta }_{\upsilon }$. With respect to the other intervals, there is not an important difference between the coverage probabilities and the average lengths of the Wald, logarithm, and Fieller CIs, especially when $n\ge 100$. When the sample size is small, the logarithmic CI and the Fieller CI have an average length slightly greater than that of the Wald CI.

Similar conclusions are obtained when ${\epsilon }_1$ and ${\epsilon }_0$ take high values. Therefore, the dependency factors ${\epsilon }_1$ and ${\epsilon }_0$ do not have an important effect on the behavior of the CIs for the difference (ratio) of the two negative predictive values.

As a conclusion, the following general rules of application can be given depending on the sample size, since the sample size is the only parameter controlled by the researcher: (a) apply the Wald CI for the difference of the positive (negative) predictive values whatever the sample size; (b) apply the Wald CI for the ratio of the two positive (negative) predictive values when the sample size is small, and apply the Wald CI, the logarithmic CI, the Fieller CI, or the bias-corrected bootstrap CI when the sample size is moderate or high.

Once some general rules of application have been established, what is better: to use a CI for the difference or a CI for the ratio? Simulation experiments have shown that the Wald CIs for the difference and the Wald CIs for the ratio have a very similar coverage probability. Furthermore, the Wald CI for the difference has a coverage probability very similar to that of the Fieller CI when the sample size is large. The Wald CIs for the difference are obtained by inverting the Wald test statistics of the tests $H_0:{\tau }_1-{\tau }_2=0$ and $H_0:{\upsilon }_1-{\upsilon }_2=0$, and the Wald CIs for the ratio are obtained by inverting the Wald test statistics of the tests $H_0:{\tau }_1/{\tau }_2=1$ and $H_0:{\upsilon }_1/{\upsilon }_2=1$. Wang et al. [4] have shown through simulation experiments that the hypothesis tests $H_0:{\tau }_1-{\tau }_2=0$ and $H_0:{\upsilon }_1-{\upsilon }_2=0$ have better asymptotic behavior than the tests $H_0:{\tau }_1/{\tau }_2=1$ and $H_0:{\upsilon }_1/{\upsilon }_2=1$. Furthermore, Wang et al. recommend using the difference-based approach as it is more straightforward and more understandable for researchers. Therefore, we recommend using a CI for the difference instead of a CI for the ratio.

4. Sample Size

The calculation of the sample size to compare parameters has great interest in Statistics. Next, we propose a procedure to determine the sample size necessary to estimate the difference between the two positive (negative) predictive values with a precision ${\varphi }_{\tau }$ $\left({\varphi }_{\upsilon }\right)$ and a confidence $100\left(1-\alpha \right)\%$. This procedure is based on the Wald CI for the difference ${\delta }_{\tau }$ $\left({\delta }_{\upsilon }\right)$, since in general terms, this CI is the interval with the best asymptotic behavior. This procedure requires having a pilot sample (or another study) in order to estimate the predictive values and their differences. If the pilot sample is not small and the Wald CI for the difference ${\delta }_{\tau }$ $\left({\delta }_{\upsilon }\right)$ contains the value 0, then the null hypothesis of equality of the predictive values is not rejected, and it does not make sense to calculate the sample size. However, if the sample is small, it may be necessary to calculate the sample size, since the Wald CI will be very wide and may contain the value 0 even if the predictive values are different. Let us considerer that ${\tau }_1\ge {\tau }_2$ $\left({\upsilon }_1\ge {\upsilon }_2\right)$ and therefore ${\delta }_{\tau }\ge 0$ $\left({\delta }_{\upsilon }\ge 0\right)$, and let ${\varphi }_{\tau }$ and ${\varphi }_{\upsilon }$ be the precisions set by the researcher ($\varphi$ must be a small value if the researcher wants high precision). Based on the asymptotic normality of ${\widehat{\delta }}_{\tau }={\widehat{\tau }}_1-{\widehat{\tau }}_1$ and of ${\widehat{\delta }}_{\upsilon }={\widehat{\upsilon }}_1-{\widehat{\upsilon }}_2$, it is verified that

$${\widehat{\delta }}_{\tau }\in {\delta }_{\tau }\pm z_{1-\alpha /2}\sqrt{Var\left({\widehat{\delta }}_{\tau }\right)}\mathrm{and}{\widehat{\delta }}_{\upsilon }\in {\delta }_{\upsilon }\pm z_{1-\alpha /2}\sqrt{Var\left({\widehat{\delta }}_{\upsilon }\right)},$$

i.e., the probability of obtaining an estimator ${\widehat{\delta }}_{\tau }$ $\left({\widehat{\delta }}_{\upsilon }\right)$ is in this interval with a probability $100\left(1-\alpha \right)\%$.

For positive predictive values, the method is as follows. Setting a precision ${\varphi }_{\tau }$, the sample size is calculated from the equation

$${\varphi }_{\tau }=z_{1-\alpha /2}\sqrt{Var\left({\widehat{\delta }}_{\tau }\right)},$$

(13)

where the variance is

$$Var\left({\widehat{\delta }}_{\tau }\right)=\frac{pq{Q}_2{\tau }_1{\overline{\tau }}_1+pq{Q}_1{\tau }_2{\overline{\tau }}_2-2\left(p{q}^2{\tau }_1{\tau }_2{\epsilon }_0+p^{2}q{\overline{\tau }}_1{\overline{\tau }}_2{\epsilon }_1+{\tau }_1{\tau }_2{\overline{\tau }}_1{\overline{\tau }}_2{Q}_1{Q}_2\right)}{npq{Q}_1{Q}_2}.$$

(14)

The proof can be seen in the Appendix A. This variance depends on the positive predictive values $\left({\tau }_i\right)$, on the disease prevalence $\left(p\right)$, on the probability of a positive result of each test $\left(Q_i\right)$, on the dependency factors $\left({\epsilon }_i\right)$, and on the sample size n. Substituting in Equation (13) the parameters with their estimators and clearing n, the sample size to estimate the difference ${\delta }_{\tau }$ with precision ${\varphi }_{\tau }$ and a confidence $100\left(1-\alpha \right)\%$ is

$$n_{\tau }=n=\frac{z_{1-\alpha /2}^2}{{\varphi }_{\tau }^2}\times \frac{\widehat{p}\widehat{q}{\widehat{\tau }}_1{\widehat{\overline{\tau }}}_1{\widehat{Q}}_2+\widehat{p}\widehat{q}{\widehat{\tau }}_2{\widehat{\overline{\tau }}}_2{\widehat{Q}}_1-2\left(\widehat{p}{\widehat{q}}^2{\widehat{\tau }}_1{\widehat{\tau }}_2{\widehat{\epsilon }}_0+{\widehat{p}}^2\widehat{q}{\widehat{\overline{\tau }}}_1{\widehat{\overline{\tau }}}_2{\widehat{\epsilon }}_1+{\widehat{\tau }}_1{\widehat{\tau }}_2{\widehat{\overline{\tau }}}_1{\widehat{\overline{\tau }}}_2{\widehat{Q}}_1{\widehat{Q}}_2\right)}{\widehat{p}\widehat{q}{\widehat{Q}}_1{\widehat{Q}}_2}.$$

(15)

Once the equation for the sample size is obtained, the method to calculate the sample size consists of the following steps:

(1)

Take pilot samples sized ${n^{\prime }}_{\tau }$, and from this sample, calculate ${\widehat{\tau }}_i$, ${\widehat{\upsilon }}_i$, ${\widehat{\epsilon }}_i$, $\widehat{p}$, ${\widehat{Q}}_i$ and the Wald CI for the difference ${\delta }_{\tau }$. If the Wald CI has a precision ${\varphi }_{\tau }$, then with the pilot sample, the precision has been reached, and the process ends. In this situation, the difference ${\delta }_{\tau }$ has been estimated with a precision ${\varphi }_{\tau }$ to a confidence $100\left(1-\alpha \right)\%$. Otherwise, go to the next step.

(2)

From the estimates obtained with the pilot sample, calculate the sample size $n_{\tau }$ applying Equation (15).

(3)

Take the sample sized $n_{\tau }$ ($n_{\tau }-{n^{\prime }}_{\tau }$ individuals are added to the initial pilot sample), and from this sample, calculate all the estimators and the Wald CI for the difference ${\delta }_{\tau }$. If the Wald CI has a precision ${\varphi }_{\tau }$, then the process ends (the precision has been reached with the new sample). If the Wald CI does not have the precision ${\varphi }_{\tau }$, then this sample is considered as a pilot sample and go to step 1.

This method to calculate the sample size n is an iterative method, which depends on the initial pilot sample and therefore does not guarantee that the difference between the positive predictive values will be estimated with the precision ${\varphi }_{\tau }$.

Sample size to estimate the difference ${\delta }_{\upsilon }$ is calculated in a similar way. In this case,

$$Var\left({\widehat{\delta }}_{\upsilon }\right)=\frac{pq{\upsilon }_1{\overline{\upsilon }}_1{\overline{Q}}_2+pq{\upsilon }_2{\overline{\upsilon }}_2{\overline{Q}}_1-2\left(p{q}^2{\overline{\upsilon }}_1{\overline{\upsilon }}_2{\epsilon }_0+p^{2}q{\upsilon }_1{\upsilon }_2{\epsilon }_1+{\upsilon }_1{\upsilon }_2{\overline{\upsilon }}_1{\overline{\upsilon }}_2{\overline{Q}}_1{\overline{Q}}_2\right)}{npq{\overline{Q}}_1{\overline{Q}}_2}$$

(16)

and the sample size $n_{\upsilon }$ to estimate the difference ${\delta }_{\upsilon }$ with precision ${\varphi }_{\upsilon }$ and a confidence $100\left(1-\alpha \right)\%$ is

$$n_{\upsilon }=\frac{z_{1-\alpha /2}^2}{{\varphi }_{\upsilon }^2}\times \frac{\widehat{p}\widehat{q}{\widehat{\upsilon }}_1{\widehat{\overline{\upsilon }}}_1{\widehat{\overline{Q}}}_2+\widehat{p}\widehat{q}{\widehat{\upsilon }}_2{\widehat{\overline{\upsilon }}}_2{\widehat{\overline{Q}}}_1-2\left(\widehat{p}{\widehat{q}}^2{\widehat{\overline{\upsilon }}}_1{\widehat{\overline{\upsilon }}}_2{\widehat{\epsilon }}_0+{\widehat{p}}^2\widehat{q}{\widehat{\upsilon }}_1{\widehat{\upsilon }}_2{\widehat{\epsilon }}_1+{\widehat{\upsilon }}_1{\widehat{\upsilon }}_2{\widehat{\overline{\upsilon }}}_1{\widehat{\overline{\upsilon }}}_2{\widehat{\overline{Q}}}_1{\widehat{\overline{Q}}}_2\right)}{\widehat{p}\widehat{q}{\widehat{\overline{Q}}}_1{\widehat{\overline{Q}}}_2}.$$

(17)

If the researcher wants to estimate ${\delta }_{\tau }$ with precision ${\varphi }_{\tau }$ and also wants to estimate ${\delta }_{\upsilon }$ with precision ${\varphi }_{\upsilon }$, at the same level of confidence, then the final sample size is $n=Max\left(n_{\tau },n_{\upsilon }\right)$. Using the largest sample size, can guarantee that the CI for the difference of the two positive predictive values and that the CI for the difference of the two negative predictive values verify the set precision for each of them.

The method for calculating the sample size depends on the values of the estimators obtained from the pilot sample. As the values of the estimators depend on each sample (and therefore vary from one sample to another), it is necessary to study how the values of the estimators affect the calculation of the sample size. Therefore, we have carried out simulation experiments to study the effect of the values of the estimators on the calculation of the sample size. These simulation experiments consisted of the following steps:

(1)

Calculate the sample size $n_{\tau }$ $\left(n_{\upsilon }\right)$, Equations (14) and (16), from the values of the parameters in the scenarios considered (Table 2 and Table 4). Therefore, these equations have been applied using the values of the parameters instead of the values of the estimators.

(2)

Generate $N=10,000$ multinomial random samples sized $n_{\tau }$ $\left(n_{\upsilon }\right)$ and whose probabilities have been calculated from Equations (7) and (8), using the parameters of the scenarios considered, and as values ${\epsilon }_i$, intermediate (50%) and high (90%) values have been considered. From each one of the N random samples, all the estimators (${\widehat{\tau }}_i$, ${\widehat{\upsilon }}_i$, ${\widehat{\epsilon }}_i$, $\widehat{p}$ and ${\widehat{Q}}_i$) have been calculated, and then, the sample size ${n^{\prime }}_{\tau }$ $\left({n^{\prime }}_{\upsilon }\right)$ has been calculated applying Equation (14) (Equation (16)).

(3)

In each scenario considered, the average sample size and relative bias have been calculated.

Table 6. Sample size for estimated the difference between the positive (negative) predictive values.

Positive Predictive Values
	$\begin{array}{c}{\tau }_1=0.90{\tau }_2=0.85{\upsilon }_1=0.95{\upsilon }_2=0.90\\ {\delta }_{\tau }=0.05{\epsilon }_1=0.021{\epsilon }_0=0.044p=0.50\end{array}$		$\begin{array}{c}{\tau }_1=0.85{\tau }_2=0.75{\upsilon }_1=0.95{\upsilon }_2=0.85\\ {\delta }_{\tau }=0.10{\epsilon }_1=0.037{\epsilon }_0=0.024p=0.25\end{array}$
	${\varphi }_{\tau }=0.025$	${\varphi }_{\tau }=0.05$	${\varphi }_{\tau }=0.025$	${\varphi }_{\tau }=0.05$
Sample size	1203	301	5048	1262
Average sample size	1204	302	5054	1267
Relative bias (%)	0.17	0.33	0.12	0.40
Negative predictive values
	$\begin{array}{c}{\tau }_1=0.90{\tau }_2=0.85{\upsilon }_1=0.95{\upsilon }_2=0.90\\ {\delta }_{\upsilon }=0.05{\epsilon }_1=0.021{\epsilon }_0=0.044p=0.50\end{array}$		$\begin{array}{c}{\tau }_1=0.85{\tau }_2=0.75{\upsilon }_1=0.95{\upsilon }_2=0.85\\ {\delta }_{\upsilon }=0.10{\epsilon }_1=0.037{\epsilon }_0=0.024p=0.25\end{array}$
	${\varphi }_{\upsilon }=0.025$	${\varphi }_{\upsilon }=0.05$	${\varphi }_{\upsilon }=0.025$	${\varphi }_{\upsilon }=0.05$
Sample size	1079	270	782	196
Average sample size	1080	272	783	198
Relative bias (%)	0.09	0.74	0.13	1.02

shows the results obtained for different precision values (2.5% and 5%, which are values that can be considered as high precision values) and $1-\alpha =0.95$. The relative biases are very small, the equations of the sample sizes provide robust values, and therefore, the pilot sample has little effect on the calculation of the sample sizes.

5. Program Cipvbdt

We have written a program in R [13] to solve the problems raised in this manuscript. The program is called “cicpvbdt” (confidence intervals to compare the predictive values of binary diagnostic tests), and it allows calculating all CIs and sample sizes. The program runs with the command “$\mathrm{cicpvbdt}\left(s_{11},s_{10},s_{01},s_{00},r_{11},r_{10},r_{01},r_{00},{\delta }_{\tau },{\delta }_{\upsilon }\right)$”. By default, the level of confidence is 95%. The program does not calculate the sample sizes when ${\delta }_{\tau }=0$ and ${\delta }_{\upsilon }=0$, and only calculates the sample sizes when ${\delta }_{\tau }>0$ and/or ${\delta }_{\upsilon }>0$. In this last situation, the program checks if the set precision is reached. The program checks that all values are valid (e.g., that there are no negative observed frequencies, etc.). The program also checks that all the parameters and their variances–covariances can be estimated. For the bias-corrected bootstrap CIs, 2000 samples with replacement are generated, and for the Monte Carlo Bayesian CIs, 10,000 random samples are generated. The results obtained are saved in a file called “results_cicpvbt.txt” in the folder from where the program is run. The program is available as Supplemental Material of this manuscript.

6. Example

The results obtained have been applied to a study on the diagnosis of colorectal cancer, using two diagnostic tests: Fecal Immunochemical Testing (FIT) and Fecal Occult Blood Testing (FOBT). The GS for the diagnosis of colorectal cancer is the biopsy.

Table 7. Observed frequencies and CIs.

Observed Frequencies
	FIT: positive		FIT: negative
Biopsy	FOBT: positive	FOBT: negative	FOBT: positive	FOBT: negative	Total
Cancer	68	18	1	13	100
Normal	4	1	2	61	68
Total	72	19	3	74	168
Results
	Positive predictive value		Negative predictive value
FIT	$0.945\pm 0.024$		$0.818\pm 0.044$
FOBT	$0.920\pm 0.031$		$0.667\pm 0.049$
$p$	${\epsilon }_1$	${\epsilon }_0$	$Q_1$	$Q_2$
0.595	0.087	0.052	0.542	0.446
$\mathrm{CIs}\mathrm{for}{\delta }_{\tau }={\tau }_1-{\tau }_2$
Wald		BCB		MCB
$\left(-0.016,0.066\right)$		$\left(-0.013,0.073\right)$		$\left(-0.045,0.105\right)$
$\mathrm{CIs}\mathrm{for}{\rho }_{\tau }={\tau }_1/{\tau }_2$
Wald	Log	Fieller	BCB	MCB
$\left(0.981,1.073\right)$	$\left(0.982,1.074\right)$	$\left(0.983,1.076\right)$	$\left(0.985,1.084\right)$	$\left(0.952,1.124\right)$
$\mathrm{CIs}\mathrm{for}{\delta }_{\upsilon }={\upsilon }_1-{\upsilon }_2$
Wald		BCB		MCB
$\left(0.081,0.222\right)$		$\left(0.089,0.231\right)$		$\left(0.049,0.248\right)$
$\mathrm{CIs}\mathrm{for}{\rho }_{\upsilon }={\upsilon }_1/{\upsilon }_2$
Wald	Log	Fieller	BCB	MCB
$\left(1.101,1.353\right)$	$\left(1.108,1.350\right)$	$\left(1.112,1.368\right)$	$\left(1.121,1.382\right)$	$\left(1.069,1.420\right)$

Wald: Wald CI. Log: logarithmic CI. Fieller: Fieller CI. BCB: bias-corrected bootstrap CI. MCB: Monte Carlo Bayesian CI.

shows the observed frequencies when the two BDTs and the GS have been applied to a sample of 168 adult men suspected of having colorectal cancer. Using the program “cicpvbdt” with the command “$\mathrm{cicpvbdt}\left(68,18,1,13,4,1,2,61,0,0\right)$”; all the results shown in Table 7 are obtained.

The estimated positive predictive values of FIT and of FOBT are 94.5% and 92.0%, and the estimated negative predictive values are 81.8% and 66.7%, respectively. Using the recommendations given in Section 3, the 95% Wald CI for the difference between the two positive predictive values contains the value zero, and therefore (with $\alpha =5\%$), the equality of the two positive predictive values is not rejected.

Regarding the negative predictive values, the 95% Wald CI does not contain the value zero, and therefore, we reject the equality of both negative predictive values. Therefore, negative predictive value of FIT is significantly greater than the negative predictive value of FOBT. The negative predictive value of FIT is (with a confidence of 95%) a value between 8.1% and 22.2% greater than the negative predictive value of FOBT. The same conclusions are obtained using the other CIs.

To illustrate the method of calculating the sample size, we are going to suppose that the clinician is interested in calculating the sample size necessary to estimate the difference between the two negative predictive values with a precision ${\varphi }_{\upsilon }=0.05$ and $1-\alpha =0.95$. The 95% Wald CI for ${\delta }_{\upsilon }={\upsilon }_1-{\upsilon }_2$ is $\left(0.081,0.222\right)$, and the precision is 0.0705 $\left(=\left(0.222-0.081\right)/2\right)$. Since ${\varphi }_{\upsilon }=0.05<0.0705$, with the sample of 168 individuals, the precision has not been reached, and therefore, the sample size must be calculated. Using the sample of 168 patients as a pilot sample and executing the command $“\mathrm{cicpvbdt}\left(68,18,1,13,4,1,2,61,0,0.05\right)”$, it is obtained that $n_{\upsilon }=338$. A sample of 338 patients is necessary to estimate the difference between the two negative predictive values with a precision ${\varphi }_{\upsilon }=0.05$ and a confidence of 95%. To the sample of 168 patients, another 170 new patients must be added. The two BDTs and the biopsy should be applied to new patients. Finally, it is necessary to recalculate the CIs with the sample of 338 patients and check if the set precision is verified.

7. Discussion

Comparison of the predictive values of two medical tests is a topic of interest in biostatistics. There are several articles that have studied hypothesis tests to solve these problems; however, the comparison of predictive values through confidence intervals has been little studied. In this manuscript, we have studied confidence intervals for the difference and for the ratio of the positive (negative) predictive values of two diagnostic tests under a paired design. We have carried out simulation experiments to study the asymptotic behaviors of the CIs, and we have given general rules of application. These rules are based on the sample size, since this is the only parameter that is set by the researcher, and also on the practical interpretation of the CIs. As a general conclusion, we recommend using the Wald interval for the difference of the two positive (negative) predictive values.

We have also proposed a method, based on the Wald CI for the difference, to calculate the sample size to estimate the difference between the two positive (negative) predictive values with a determined precision and confidence. This method starts from an initial pilot sample, and then the sample size is calculated from the estimators obtained with the initial sample. This method depends on the estimators of the pilot sample, so we have carried out simulation experiments to study the effect of the pilot sample on the sample size. The results obtained in these experiments have shown that the pilot sample does not have any important effect on the calculation of the sample size, and that therefore, the method has practical validity.