Computational methods to simultaneously compare the predictive values of two diagnostic tests with missing data: EM-SEM algorithms and multiple imputation

1 Computational methods to simultaneously compare the predictive values of two diagnostic tests with missing data: EM-SEM algorithms and multiple imputation J.A. Roldán-Nofuentes Biostatistics, School of Medicine, University of Granada, 18016, Spain Email: [email protected] This is an Accepted Manuscript of an article published by Taylor & Francis in the Journal of Statistical Computation and Simulation on 2021, available at https://www.tandfonline.com/doi/abs/10.1080/00949655.2021.1926461 Abstract. Predictive values are measures of the clinical accuracy of a binary diagnostic test, and depend on the sensitivity and the specificity of the diagnostic test and on the disease prevalence among the population being studied. This article studies hypothesis tests to simultaneously compare the predictive values of two binary diagnostic tests in the presence of missing data. The hypothesis tests were solved applying two computational methods: the expectation maximization and the supplemented expectation maximization algorithms, and multiple imputation. Simulation experiments were carried out to study the sizes and the powers of the hypothesis tests, giving some general rules of application. Two R programmes were written to apply each method, and they are available as supplementary material for the manuscript. The results were applied to the diagnosis of Alzheimer’s disease. Key words: EM and SEM algorithms, Missing data, Multiple imputation, Partial verification, Predictive values. Mathematics Subject Classification: 62P10, 6207. 2 1. Introduction A diagnostic test is a medical test that is applied to a patient to determine the presence or absence of a certain disease. When the result of a diagnostic test is positive or negative, the diagnostic test is called a binary diagnostic test (BDT). The mammography for breast cancer is an example of a BDT. The clinical effectiveness of a BDT is measured in terms of two parameters: the positive predictive value and the negative predictive value. Positive predictive value (  ) is the probability of a patient having the disease when the result of the BDT is positive, and the negative predictive value (  ) is the probability of the patient not having the disease when the result of the BDT is negative. The predictive values (PVs) depend on the sensitivity (Se) and on the specificity (Sp) of the BDT and on the disease prevalence (p) among the population studied, i.e.     and 11 p Se q Sp τυ p Se q Sp p Se q Sp           , (1) where 1qp . While Se and Sp quantify how well the BDT reflexes the true disease status, the PVs quantify the clinical value of the BDT, since the patient is more interested in knowing the probability of having or not having the disease given a result of the diagnostic test. The parameters of a BDT are estimated in relation to a gold standard (GS), which is a medical test which determines without any errors whether or not the patient has the disease. A biopsy for breast cancer is an example of a GS. In clinical practice, the most common sample design to compare the PVs of two BDTs is paired design [1, 2]. This type of design consists of applying the two BDTs to all of the individuals in a sample sized n whose disease status is known through the application of a GS. The comparison of the PVs of two BDTs subject to paired design has been the subject of different studies in statistics literature. Leisenring et al [3], Wang et al [4], Kosinski [5] and Tsou [6] have studied asymptotic methods to compare the two positive PVs and the negative 3 PVs independently, i.e. solving the two hypothesis tests 0 1 2 :H   and 0 1 2 :H   each one of them to an α error. The Kosinski method has a better asymptotic performance (in terms of type I error and power) than the methods of Leisenring et al and of Wang et al. The method of Tsou leads to the same results as the Kosinski method. Roldán Nofuentes et al [7] studied a global hypothesis test to simultaneously compare the PVs of two BDTs, i.e. solving the global hypothesis test   0 1 2 1 2 : and H      vs   1 1 2 1 2 : and/or H      , and proposed a method based on chi-squared distribution and multiple comparisons. These authors have demonstrated that the comparison of the positive PVs (negative PVs) of two BDTs subject to a paired design must be carried out simultaneously, solving the global hypothesis test   0 1 2 1 2 : and H      . They have also demonstrated that the comparison of the PVs is made independently, i.e. solving the tests 0 1 2 :H   and 0 1 2 :H   each one of them to an  error, the results may be mistaken. In Appendix A these methods are summarized. When assessing or comparing parameters of BDTs it is not common for the GS to be applied to all of the individuals in the sample, leading to the problem known as partial disease verification [8, 9]. Therefore, if the GS consists of a costly test or one which means some risk for the individual, then it will not be applied to all of the individuals in the sample, and consequently the true disease status is unknown for a subgroup of individuals. When comparing parameters of two BDTs in the presence of partial disease verification, it is common to assume that the verification process is missing at random (MAR). This assumes that the process to verify the disease status of an individual through the application of the GS only conditionally depends on the results of the two BDTs and it does not depend on the disease status of the individual. Subject to the MAR assumption, there are many different studies in statistics literature which compare parameters of two or more BDTs in the presence of partial verification. Zhou [9] studied a hypothesis test to compare the sensitivities (specificities) of two BDTs applying the method of maximum likelihood. Roldán Nofuentes 4 and Luna [10] studied the individual comparison of the PVs of two BDTs applying the maximum likelihood (ML) method. Marín-Jiménez and Roldán-Nofuentes [11] extended the study of Roldán-Nofuentes and Luna [10] to the case of more than two BDTs comparing the PVs simultaneously applying the ML method. Harel and Zhou [12] compared the sensitivities (specificities) of two BDTs through confidence intervals applying multiple imputation (MI). Roldán-Nofuentes and Luna [13] compared the sensitivities and the specificities independently, as well as the PVs, of two BDTs applying the expectation maximization (EM) algorithm and the supplemented expectation maximization (SEM) algorithm. In this article, hypothesis tests are studied to simultaneously compare the PVs of two BDTs when in the presence of partial disease verification the missing data mechanism is MAR, applying two computational methods: the EM and SEM algorithms, and MI. The EM algorithm is a classic method for estimating parameters in the presence of missing data. The advantage of the EM algorithm over the ML method [11] is that the EM algorithm can be applied when some observed frequency is zero, while the ML method [11] cannot be applied in this situation. Regarding the MI, this method offers better results than the ML method when comparing the sensitivities (specificities) of two BDTs in the presence of missing data. Therefore, it is convenient to evaluate this method to solve the problem presented here. MI cannot be applied when some observed frequency is zero. Therefore, with both types of computational methods we seek to solve the global hypothesis global test to simultaneously compare the positive PVs and the negative PVs of the two BDTs, i.e. 0 1 2 1 2 1 1 2 1 2 and and/or H:τ τ υ υ H:τ τ υ υ .   (2) In Section 2, we solve the global test applying the EM and SEM algorithms, and in Section 3 the same test is solved applying MI. In Section 4, simulation experiments are carried out to study the type I error and the power of the previous tests with each one of the two methods (EM-SEM algorithms and MI), and some general rules of application are given. In Section 5, 5 two programmes written in R are presented to solve the problem posed by applying each computational method. In Section 6, the results were applied to a real example of the diagnosis of Alzheimer’s disease, and in Section 7 the results obtained are discussed. 2. EM and SEM algorithms Let us consider two BDTs that are applied to all of the individuals in a random sample sized n, and let us also consider a GS that is only applied to a subgroup of the sample. This situation leads to the observed frequencies in Table 1a, where the variable h T models the result of the hth BDT ( 1 h T when the result is positive and 0 h T when it is negative), the variable V models the verification process ( 1V when the disease status of an individual is verified with the GS and 0V when it is not), and the variable D models the result of the GS ( 1D when the individual verified has the disease and 0D when this individual does not). In Table 1a, ij a is the number of individuals with the disease among whom 1 Ti and 2 Tj , ij b is the number of individuals without the disease among whom 1 Ti and 2 Tj , and ij c is the number of individuals with an unknown disease status among whom 1 Ti and 2 Tj , with , 0,1ij . Let 1 ,0ij ij aa   , 1 ,0 ij ij bb   , 1 ,0 ij ij cc   , ij ij ij ij n a b c   and 1 ,0ij ij nn   . For the hth BDT, let   11 hh Se P T D   ,   00 hh Sp P T D   ,   11 hh P D T     and   00 hh P D T     , with 1,2h . Let   1p P D be the disease prevalence and   10q p P D    . From the expressions (1) each Se and Sp is written, in terms of the PVs and of p, as     and h h h h hh hh qp Se Sp pY qY       , (3) where 1 h h h Y     . 6 ==== INSERT TABLE 1 HERE ==== In the presence of partial disease verification, the verification probabilities are designed as   12 1 , , ijk P V T i T j D k      , i.e. ijk  is the probability of verifying the GS the disease status of an individual for whom 1 Ti , 2 Tj and Dk , with , , 0,1i j k  . Assuming that the missing data mechanism is MAR, i.e. that the probability of verifying the disease status of an individual only conditionally depends on the results of the two BDTs and does not depend on the result of the GS, it is verified that   12 1, ijk ij P V T i T j      . Supposing the MAR assumption, the data from Table 1a are the product of a multinomial distribution sized n and the probabilities are written in terms of the PVs as                             12 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 2 1 12 11 11 1 1 1 1 2 2 2 2 11 11 1, 1, , 11 , 1, 0, , 11 ij i i i j j j ij ij ij i i i i j j j j ij i i i j j j i ij i i i i P V D T i T j q p q p pp p Y Y p p Y Y P V D T i T j q p q p qq q Y Y                                                                    0 11 22 12 , 1 0, , , j ij j j j j ij ij ij ij ij q q Y Y P V T i T j                    (4) where     1 2 1 2 1 1212 1qq p YY           and       1 2 1 2 0 0212 1 1 1qq q YY             . with 1 ij   if ij and 1 ij   if ij . Parameters 1  and 0  are the covariances [14] between the two BDTs when 1D and when 0D respectively, verifying that     1 1 1 2 2 12 1 1 max , qq pY pY          and     0 1 1 2 2 12 1 1 max 1 , 1 pp qY qY                         . 7 If 10 1   , then the two BDTs are conditionally independent on the disease, a situation which is not realistic in practice, and therefore it must be verified that 11   and/or 01   . An EM algorithm is then proposed to estimate the PVs of the two BDTs. 2.1. EM algorithm In the 34 table of observed frequencies, the missing information is the true disease status of the individuals who are not verified with the GS, i.e. the missing data is the value of the variable D for the individuals among whom 0V . This information is reconstructed in the E step of the algorithm and in the M step the values of the maximum likelihood estimators are imputed. Let us assume that that among the ij c individuals who are not verified   0V , ij d have the disease and ij ij cd do not have it, with , 0,1ij . Then the table of observed frequencies can be expressed in the form of a 24 table with frequencies ij ij ad for 1D and ij ij ij b c d for 0D (Table 1b). Let   1 1 2 2 1 0 , , , , , , T p       θ be the vector of parameters. From the table of complete data, the log-likelihood function based on n individuals is           11 , 0 , 0 log log , ij ij ij ij ij ij ij i j i j l a d b c d         θ (5) where   12 , , 1 ij P T i T i D     and   12 , , 0 ij φ P T i T i D    . The expressions of these probabilities in terms of the PVs are shown in Appendix B of the supplementary material. The components of the vector θ are going to be estimated applying the EM algorithm. Therefore, let us suppose that   k ij d is the value of ij d in the kth iteration of the EM algorithm, and     1 ,0 kk ij ij dd   . The values of the MLEs in the kth iteration are calculated through the following equations 8                                                     11 1 1 1 1 0 0 0 00 10 11 00 01 11 2 1 1 2 0 0 0 00 01 11 00 10 11 11 111 1 1 1 1 00 0 11 ˆ ˆ, , 11 ˆ ˆ, , ˆ ˆ, , ˆ k k k k j j j j j jj k k k k i i i i i ii kk k kk kk i i j j ij k k a d b c d n n n n a d b c d n n n n a d a d ad pna d a d b c d b                                               11 11 11 11 1 1 1 1 1 1 00 . k kk i i i j j j ij cd b c d b c d                  (6) The estimators in the   1k th iteration of the algorithm are calculated applying equations (6) substituting the superindex k with 1k , where         1ˆ, , 0,1 ˆˆ k kij ij ij kk ij ij d c i j     , when   ˆk ij  and   ˆk ij  are the estimators of the probabilities ij  and ij  in the kth iteration of the algorithm, and which are calculated substituting in the expressions of ij  and ij  (see Appendix B of supplementary material) the parameters with their respective estimators obtained in the kth iteration. As initial value   0 ij d one can take any value between 0 and ij u . The EM algorithm stops when the difference between the values of the log-likelihood functions of two consecutive iterations is lower than a value  , for example 10 10    or 12 10    . If the EM algorithm has converged in K iterations, we denote through   1 1 2 2 1 0 ˆˆ ˆ ˆ ˆ ˆ ˆ ˆ , , , , , , T p       θ the final estimators obtained. The estimators of the PVs obtained by applying the EM algorithm converge to the maximum likelihood estimators [10, 11] (proof can be seen in Appendix C of the supplementary material). The variances-covariances of ˆ θ are then estimated applying the SEM algorithm [15]. 9 2.2. SEM algorithm The estimation of the matrix of asymptotic variances-covariances of ˆ θ can be obtained applying the SEM algorithm [15], which is a computational method which estimates the variances-covariances matrix of a vector of estimators from the calculations performed in the application of the EM algorithm. Let ˆ θ be the variances-covariances matrix of ˆ θ , Dempster et al [16] demonstrated that   1 1 ˆoc I I DM      θ , (7) where I is the identity matrix and 1 mis oc DM I I  , when oc I is the Fisher information matrix of the complete data and mis I is the Fisher information matrix of the missing data. The SEM algorithm consists of three phases: 1) assessment of the matrix 1 oc I , 2) assessment of the matrix DM, and 3) assessment of the matrix ˆ θ . The main phase is to calculate the elements of the DM matrix. The following three phases are then analysed. The first phase consists of assessing the 1 oc I . This matrix is the inverse of the Fisher information matrix of the complete data, i.e.   2 oc ij l I    θ , where   lθ is the function (5) and each i  is one of the parameters of θ . This matrix is calculated from the last table after the application of the EM algorithm, substituting the parameters with their corresponding estimations obtained in the last iteration of the EM algorithm. If the EM algorithm has converged in K iterations, then the frequencies of the last 24 table are   K ij ij ad for 1D and   K ij ij ij b c d for 0D . The second part of the SEM algorithm consists of calculating the DM matrix. The elements of this matrix, denoted as ij r , , 1,...,7ij , are obtained applying the following algorithm: 16 and the combined p-value is   2, 2l p value P F F   with     2 31 2 2 1 1 M l M r     . c) Combined likelihood-ratio tests A third method to solve the global test is combining likelihood-ratio tests [25]. Let             ( ) ( ) ( ) ( ) 11 10 01 00 11 10 01 00 , , , , , , , T m mmmmmmmm x x x x y y y yz be the vector of imputed frequencies in the mth complete dataset. Let   m ij  and   m ij  be the probabilities corresponding to each cell of the imputed 24 table, whose expressions are similar to those given in Appendix A of supplementary material adding the superindex (m) to all of the parameters. Let                     11 10 01 00 11 10 01 00 , , , , , , , T m m m m m m m m m         ψ be the corresponding vector of probabilities. The complete-data log-likelihood function is                 11 , 0 , 0 ; log log m m m m m m m ij ij ij ij i j i j L x y             ψz . Maximizing this function it holds that     ˆmm ij ij xn   and     ˆmm ij ij yn   , which are the nonrestricted estimators of   m ψ (i.e.   ˆm ψ ), with , 0,1ij . If in the mth set of imputed data the null hypothesis           0 1 2 1 2 : and m m m m H      is true, then it is easy to show that the loglikelihood function is                                         * * * 0 0 11 11 10 01 10 00 00 * * * 11 11 10 01 10 00 00 , log log log log log log , m m m m m m m m m m m m m m m m m L x x x x y y y y                                       ψz where   *m ij  and   *m ij  are the probabilities subject to the null hypothesis. Maximizing this function it holds that                     ** 10 01 10 01 , if , if ˆˆ and 2 , if 2 , if , mm ii ii mm ij ij m m m m x n i j y n i j x x n i j y y n i j                     17 which are the estimators of   m ψ subject to the null hypothesis (i.e.   0 ˆm ψ ), with , 0,1ij . Performing algebraic operations, the likelihood-ratio test statistic for the global test           0 1 2 1 2 : and m m m m H      is                                                   3 0 0 10 01 10 01 10 01 10 01 10 01 10 01 10 01 10 01 ˆˆ 2 ; ; 2 2 2 2 2 log log log log . m m m m m m m m m m m m m m m m m m m m m m m F L L x x y y x x y y x x x x y y y y                                                 ψ z ψ z Let   33 1 1Mm m FF M  be the average of likelihood-ratio statistics, and let   1 1ˆ Mm m M  ψψ and   00 1 1ˆ Mm m M  ψψ be the vectors whose components are the measures of the estimators in the M sets of imputed data subject to the non-restricted model and subject to the null hypothesis, respectively. Let             * 3 0 0 1 2;; Mm m m m m F L L M    ψ z ψ z , which is the measure of the likelihood-ratio statistics each one of which is assessed in ψ and 0 ψ , and let     * 3 3 3 1 21 M r F F M    . Finally, the likelihood-ratio test statistic for the global hypothesis test is [25]   * 3 3 3 21 F Fr  , which is distributed according to a F-distribution with 2 and             2 3 2 1 3 2 4 2 6 1 , if 2 1 4 1 31 1 , if 2 1 4 2 M MM Mr l M r M                   degrees of freedom. 18 3.3. Individual tests As with the EM-SEM algorithms, the global hypothesis test can also be solved from the individual tests applying MI with methods that compare the positive PVs and the negative PVs independently or with a method of multiple comparisons. The methods that are going to be considered to individually compare the PVs are those of Leisenring et al [3], Wang et al [4] and Kosinski [5]. For the method by Wang et al and for the method by Kosinski, the combination of results is achieved applying the rules of Rubin [19]. For the method by Leisenring et al, the combination of results is achieved calculating the average statistic. The test statistics of the method by Wang et al and of the method by Kosinski are of the type   ˆˆ ˆ Var  , and therefore the combination of results is achieved applying the rules of Rubin. Nevertheless, in the case of the method by Leisenring et al, the test statistic to compare the equality of the two positive (negative) PVs is not of the type   ˆˆ ˆ Var  , and therefore the rules of Rubin cannot be applied. For the method of Wang et al and for the method of Kosinski the results obtained are then combined applying the rules of Rubin. Firstly, we calculate the overall estimate of the difference between the two positive PVs, i.e.   1 1ˆ Mm m M    , where       12 ˆ ˆ ˆ m m m     . The variance of  is     1 ˆˆ1 Var Var B M    , where       1 1ˆ ˆˆ Mm m Var Var M    is the within imputation variance (the average of the complete data variance estimates) and     2 1 1ˆ 1 Mm m BM     is the between imputation variance (the variance of the complete data point estimates). Finally, the test statistic for the test 0 1 2 :H   is   ˆ Var   , whose 19 distribution is (Rubin, 1987) a t-distribution with     ˆ 11 1 Var M vM MB          degrees of freedom. The comparison of the negative PVs is made in a similar way substituting  with  . For the method by Leisenring et al, in the mth complete dataset the test statistic for 0 1 2 :H   is   m z  (its expression can be seen in Appendix A), whose distribution is a normal standard one when the sample size is large. Then for the Central Limit Theorem, the average of all of the test statistics   1 1Mm m zz M    has a normal standard distribution when M is large. The process for the test 0 1 2 :H   is similar to the previous one. 4. Simulation experiments Monte Carlo simulation experiments were carried out to study the type I errors and the powers of the hypothesis tests studied in Sections 2 and 3, as well as the relative biases of the estimators of the PVs obtained with both methods. These experiments consisted of generating 10000N random samples of multinomial distributions sized   50,100,200,500,1000,2000n , and whose probabilities were calculated from equations (4). As PVs we considered the values   0.70,0.75,0.80,0.85,0.90,0.95 , which are values that appear quite frequently in clinical practice, and as disease prevalence we took the values   25%,50%,75%p . Once the PVs and p are set, the Se and the Sp of each BDT were calculated from equations (3). As values of the covariances   i  we considered intermediate and high values. Finally, the probabilities of the multinomial distributions were calculated applying equations (4). Therefore, the probabilities of the multinomial distributions were calculated from the PVs, and there was no previous setting of the values of sensitivities and specificities of the BDTs. The simulation experiments were designed in such a way that in all 20 of the samples generated it is possible to apply the EM-SEM algorithms and MI. Therefore, if in a sample a frequency ij a (or ij b ) is zero then it is not possible to apply MI (it is not possible to apply logistic regression to impute the missing data), then this sample was discarded and another one was generated instead until completing the N samples. Regarding the EM-SEM algorithms, we established as a stop criterion 12 10    and 6 10    respectively, and as initial values of the EM algorithm the values   02 ij ij dc are used. Regarding MI, for each one of the N random samples 20M complete data sets were generated and 100 cycles were performed. In the first phase simulations were made considering 20M and 50M and performing 100 and 200 cycles in each case, obtaining very similar results; therefore, to reduce computation time we finally considered 20M and 100 cycles. For all of the study, we set as the nominal error 5%   , and considered that a method overwhelms the nominal error or exceeds it too much when its type I error is higher than 7%. The simulation experiments were carried out with R [26] and for the MI we used the “mice” library [27]. Therefore, in the simulation experiments we studied and compared the type I errors and the powers of sixteen different methods to solve the global hypothesis test (2). Of the sixteen methods, four are based on the EM-SEM algorithms and the other twelve on MI. The methods based on the EM-SEM algorithms are: (a) a global hypothesis global test based on the chisquared distribution with 5%   ; (b) an individual comparison of the positive PVs and the negative PVs with 5%   , Bonferroni and Holm. The twelve methods based on MI are: (a) a global hypothesis test applying the Wald method, the combination of p-values and the combined likelihood ratio tests, all of them with 5%   ; (b) an individual comparison of the positive PVs and the negative PVs applying the methods of Leisenring et al, Wang et al and Kosinski, each of them with 5%   , Bonferroni and Holm. 21 4.1. EM-SEM algorithms Table 2 shows some of the results obtained for the type I errors of the different methods applying the EM-SEM algorithms. In this Table, EM-SEM-Global consists of solving the global hypothesis test, EM-SEM-Individual consists of comparing the PVs solving the individual hypothesis tests each of them to an error 5%   , and EM-SEM-Bonferroni consists of comparing the PVs solving the individual hypothesis tests along with the Bonferroni method to an error 5%   . The results obtained applying the Holm method to an error 5%   are not shown as they are practically the same as those obtained with the Bonferroni method. In general terms, the type I errors of all of the methods increase when the verification probabilities increase, and decrease when the covariances i α increase. In general terms, all of the methods are very conservative when the sample size is small   50n or moderate   100 200n . When the sample size is large   500n , depending on the verification probabilities and on the covariances i  , the global test has a type I error that fluctuates around the nominal error. On some occasions, above all when the verification probabilities are low or the covariances are high, the type I error of the global test may slightly exceed the nominal error without actually overwhelming it. Method EM-SEM- Individual may overwhelm the nominal error, above all when the sample size is large. Method EM-SEM-Bonferroni has a type I error whose behaviour is very similar to that of the global test. ==== INSERT TABLE 2 HERE ==== Regarding the powers of these methods, Table 3 shows some results. The power of these methods increases when there is an increase in in the verification probabilities, whereas the 22 covariances i  do not have a clear effect upon the power. All of the methods have a very small power when the sample size is small   50n or moderate   100 200n , and it is necessary to have a large sample size (depending on the verification probabilities) so that the power is higher (over 80%). Although there are no clear rules, in general terms the global test is normally more powerful than Method 2. Regarding the method EM-SEM-Individual, there are also no clear rules, sometimes the global test is more powerful and on other occasions the method EM-SEM-Individual is more powerful (it is a method that easily overwhelms the nominal error), depending on the values that the PVs take. ==== INSERT TABLE 3 HERE ==== From the results of the simulation experiments applying the EM-SEM algorithms, the method to compare the PVs of the two BDTs with the best asymptotic behaviour is the global test, since its type I error does not exceed the nominal error too much and, in general terms, it has more power than the method EM-SEM-Bonferroni (this is a method whose type I error also does not exceed the nominal error too much). Method EM-SEM-Individual may exceed the nominal error too much and, therefore, lead to false significances. The same previous results are obtained if the global test is solved by applying the ML method [10, 11], because the estimators obtained through the EM-SEM algorithms converge to the ML estimators. 4.2. Multiple Imputation Table 4 shows the results obtained for the type I errors through MI for the same scenarios as Table 2. In this table, Leisenring-Individual, Wang-Individual and Kosinski-Individual, refers to the individual comparison of the PVs applying the method of Leisenring et al with 5%   , 23 the method of Wang et al with 5%   and the method of Kosinski with 5%   , respectively. Similarly, Leisenring-Bonferroni, Wang-Bonferroni and Kosinski-Bonferroni, refers to the individual comparison of the PVs applying the method of Leisenring et al, Wang et al and Kosisnki, respectively, along with the Bonferroni method and 5%   . The results obtained applying the Holm method are not shown as they are practically identical those obtained with Bonferroni. As with the EM-SEM algorithms, the type I errors of all the methods based on MI increase when the verification probabilities increase, and decrease when the covariances i  increase. Regarding the global tests, the type I error of the test based on the combination of p-values is very similar to the type I error of the combined likelihood-ratio tests, both of which fluctuate around the nominal error when the sample size is large. The global test based on the Wald test is very conservative (even when the sample size is large), and its type I error is smaller than that of the other two methods. The methods based on the individual comparisons to an error 5%   (Leisenring- Individual, Wang-Individual and Kosinski-Individual) have type I errors that may exceed the nominal error too much, above all when the sample size is large. Therefore, these methods may lead to an excess of false significances. Regarding the methods based on the individual tests along with Bonferroni, the Leisenring-Bonferroni method has a type I error which may exceed the nominal error too much when the sample size is large. Kosinski-Boferroni method has a type I error with better fluctuations around the nominal error (with exceeding it too much) than the Wang-Bonferroni method, above all when the sample size is large. In general terms, there is no important difference between the type I error of the Kosinski-Bonferroni method and the type I error of the global test based on the combination of p-values (or combined likelihood-ratio tests), above all when the sample size is large. 24 ==== INSERT TABLE 4 HERE ==== Regarding the powers, Table 5 shows the results for the same scenarios given in Table 3. The power of all of the methods increases when the verification probabilities increase, whereas the covariances i  do not have a clear effect upon the power. In general terms, the global test based on the combination of p-values is more powerful than the combined likelihood-ratio tests, and the difference is greater when the verification probabilities are low than when they are high. Moreover, both methods are more powerful than the Wald test (as this test is very conservative in relation to the other two). Comparing the global test based on the combination of p-values in relation to the Leisenring-Individual, Wang-Individual and Kosinski-Individual methods, there are no clear rules about their behaviour. Sometimes the global test is more powerful and on other occasions these methods are more powerful (methods which may clearly overwhelm the nominal error), depending on the verification probabilities and on the values that the PVs take. Regarding the Leisenring-Bonferroni method, in general terms the global test based on the combination of p-values is less powerful, due to the fact that the type I error of the Leisenring-Bonferroni method (which may clearly overwhelm the nominal error) is greater than that of the global test (which does not overwhelm the nominal error). Regarding the Wang-Bonferroni method and the Kosinski-Bonferroni method, there is no important difference between their powers, and these powers are a little lower than those of the global test based on the combination of p-values, above all when the sample size is large. Having analysed the results of the simulation experiments applying MI, the method to compare the PVs of two BDTs with the best asymptotic behaviour is the global test based on the combination of p-values, since its type I error does not overwhelm the nominal error and 25 its power is somewhat higher than that of the other methods which do not overwhelm the nominal error. ==== INSERT TABLE 5 HERE ==== 4.3. Relative biases Table 6 shows some results for the relative biases of the estimators of the PVs applying the EM algorithm and applying MI. The relative biases decrease when the verification probabilities increase whereas the covariances i  have practically no effect on the estimators obtained applying both methods. In general terms, the difference between the relative biases obtained with both methods is very small, and therefore both methods lead to estimations which on average are very similar to each other. ==== INSERT TABLE 6 HERE ==== 4.4. EM-SEM algorithms when some frequency is zero The EM-SEM algorithms can be applied when some frequency ij a or ij b is equal to zero. Simulation experiments have been carried out to study the asymptotic behavior of these algorithms in this situation. These experiments have been designed in a similar way to the previous case, but the samples in which some frequency ij a or ij b is equal to zero have not been eliminated. Table 7 shows the type I errors and the powers for the same scenarios given in Tables 2 and 3. In general terms, the conclusions are the same as those given in Section 4.1. ==== INSERT TABLE 7 HERE ==== 32 20M complete datasets to save computation time. The simulation experiments demonstrated that the global test based on the combination of p-values (which is the test based on MI with the best asymptotic behaviour) has a type I error which is very similar to the global test based on the EM-SEM algorithms. Regarding power, that of the global test based on the combination of p-values is a little higher than that of the global test based on the EMSEM algorithms when the sample is small or moderate, and they are very similar when the sample is large. Therefore, the number of complete datasets was sufficiently large, and did not have any negative effect on the size and the power of the global test. Disclosure statement No potential conflict of interest was reported by the author. References 1. Pepe, M.S. (2003). The Statistical Evaluation of Medical Tests for Classification and Prediction. New York: Oxford University Press. 2. Zhou, X.H., Obuchowski, N.A. and McClish, D.K. (2011). Statistical Methods in Diagnostic Medicine (Second Edition). New Jersey: John Wiley & Sons. 3. Leisenring, W., Alonzo, T. and Pepe, M.S. (2000). Comparisons of predictive values of binary medical diagnostic tests for paired designs. Biometrics, 56, 345-351. 4. Wang, W., Davis, C.S. and Soong, S.J. (2006). Comparison of predictive values of two diagnostic tests from the same sample of subjects using weighted least squares. Statistics in Medicine, 25, 2215-2229. 5. Kosinski, A.S. (2013). A weighted generalized score statistic for comparison of predictive values of diagnostic tests. Statistics in Medicine, 32, 964-977. 33 6. Tsou, T.S. (2018). A new likelihood approach to inference about predictive values of diagnostic tests in paired designs. Statistical Methods in Medical Research, 27, 541-548. 7. Roldán-Nofuentes, J.A., Luna del Castillo, J.D. and Montero-Alonso, M.A. (2012). Global hypothesis test to simultaneously compare the predictive values of two binary diagnostic tests. Computational Statistics and Data Analysis, 56, 1161-1173. 8. Begg, C.B. and Greenes, R.A. (1983). Assessment of diagnostic tests when disease verification is subject to selection bias. Biometrics, 39, 207-215. 9. Zhou, X.H. (1998). Comparing accuracies of two screening tests in a two-phase study for dementia. The Journal of the Royal Statistical Society, Series C Applied Statistics, 47, 135- 147. 10. Roldán Nofuentes, J.A. and Luna del Castillo, J.D. (2008). The effect of verification bias on the comparison of predictive values of two binary diagnostic tests. Journal of Statistical Planning and Inference, 138, 950-963. 11. Marín-Jiménez, A.E. and Roldán-Nofuentes, J.A. (2014). Global hypothesis test to compare the likelihood ratios of multiple binary diagnostic tests with ignorable missing data. SORT - Statistics and Operations Research Transactions, 38, 305-324. 12. Harel, O. and Zhou, X.H. (2007). Multiple imputation for the comparison of two screening tests in two-phase Alzheimer studies. Statistics in Medicine, 26, 2370-2388. 13. Roldán Nofuentes, J.A. and Luna del Castillo, J.D. (2008). EM algorithm for comparing two binary diagnostic tests when not all the patients are verified. Journal of Statistical Computation and Simulation, 78, 19-35. 14. Berry G., Smith, C., Macaskill, P. and Irwig L. (2002). Analytic methods for comparing two dichotomous screening or diagnostic tests applied to two populations of differing disease prevalence when individuals negative on both tests are unverified. Statistics in Medicine, 21, 853-862. 34 15. Meng, X. and Rubin, D.B. (1991). Using EM to obtain asymptotic variance-covariance matrices: the SEM algorithm. Journal of the American Statistical Association, 86, 899-909. 16. Dempster, A., Laird, N. and Rubin, D.B. (1977). Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society Series B, 39, 1-38. 17. Bonferroni, C.E. (1936). Teoria statistica delle classi e calcolo delle probabilità, Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze, 8, 3-62. 18. Holm, S. (1979). A simple sequential rejective multiple testing procedure. Scandinavian Journal of Statistics, 6, 65-70. 19. Rubin, D.B. (1987). Multiple Imputation for Nonresponse in Surveys. New York: Wiley. 20. Schafer, J.L. (1997). Analysis of Incomplete Multivariate Data. London: Chapman and Hall. 21. Little, R.J.A. and Rubin, D.B. (2002). Statistical analysis with missing data (Second Edition). New Jersey: Wiley. 22. White, I.R., Royston, P. and Wood, A.M. (2011). Multiple imputation using chained equations: issues and guidance for practice. Statistics in Medicine, 30, 377-399. 23. Li, K.H, Raghunathan, T.E. and Rubin, D.B. (1991). Large-sample significance levels from multiply imputed data using moment-based statistics and an F reference distribution. Journal of the American Statistical Association, 86, 1065-1073. 24. Li, K.H, Meng, X.L., Raghunathan, T.E. and Rubin, D.B.(1991). Significance levels from repeated p-values with multiply-imputed data. Statistica Sinica, 1, 65-92. 25. Meng, X.L., Raghunathan, T.E. and Rubin, D.B. (1992). Performing likelihood ratio tests with multiply-imputed data sets. Biometrika, 79, 103-111. 35 26. R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/. 27. van Buuren, S. and Groothuis-Oudshoorn, K. (2011). Mice: multivariate imputation by chained equations in R. Journal of Statistical Software, 45, 1-67. 28. Hall, K.S., Ogunniyi, A.O., Hendrie, H.C., Osuntokun, B.O., Hui, S.L., Musick, B., Rodenberg, C.S., Unverzagt, F.W., Guerje, O., and Baiyewu, O. (1996). A cross-cultural community based study of dementias: methods and performance of survey instrument. International Journal of Methods in Psychiatric Research, 6, 129-142. 36 Supplementary material of the manuscript: Computational methods to simultaneously compare the predictive values of two diagnostic tests with missing data: EM-SEM algorithms and multiple imputation Appendix A Let us consider that two BDTs are applied to all the individuals in a random sample sized n, whose disease status (present or absent) is known through the application of a GS. Let ij x   ij y be the numbers of diseased (non-diseased) individuals among whom Test 1 leads to a result i and Test 2 leads to result j, with , 0,1ij (0 indicates a negative result and 1 a positive result). We now summarize the methods of Leisenring et al [3], Wang et al [4] and Kosinski [5]. The Tsou method [6] is not considered since it is equivalent to the Kosinski method. Method of Leisenring et al Leisenring et al [3] studied the comparison of the positive and negative PVs of two binary tests through marginal regression models, and they were able to estimate these models separately or jointly using GEE models. Leisenring et al deduced score statistics to compare the positive and negative PVs of two binary tests in paired designs. Using the notation from the previous Section, the score statistic for the test 0 1 2 :H   is                   11 1 01 1 10 1 2 2 2 2 2 2 2 2 2 2 2 2 11 1 1 01 1 1 10 1 1 11 1 1 01 1 1 10 1 1 1 2 1 1 1 2 1 1 1 1 2 1 x Z x Z x Z zx D Z x D Z x D Z y D Z y D Z y D Z                   and the score statistic to compare the test 0 1 2 :H   is                   00 2 10 2 01 2 2 2 2 2 2 2 2 2 2 2 2 2 00 2 2 10 2 2 01 2 2 00 2 2 10 2 2 01 2 2 1 2 1 1 1 2 1 1 1 1 2 1 y Z y Z y Z zy D Z y D Z y D Z x D Z x D Z x D Z                   . 37 Score statistics have a normal distribution when the null hypothesis is true, and where 11 01 11 01 1 11 01 10 11 10 01 22 x x y y Zx x x y y y         , 11 01 10 1 11 01 10 11 10 01 2 22 x x x Dx x x y y y       . 00 10 00 10 2 00 01 10 00 01 10 22 x x y y Zx x x y y y         and 00 01 10 2 00 01 10 00 01 10 2 22 y y y Dx x x y y y       . Method of Wang et al Wang et al [4] studied the comparison of the PVs of two binary tests through a weighted least square method and compared their method to that of Leisenring et al, before recommending the comparison of the PVs using the weighted least square method based on the difference between the two positive (negative) PVs. The test statistics for 0 1 2 :H   and 0 1 2 :H   are       12 1 2 1 2 ˆˆ ˆ ˆˆ ˆ ˆ ˆ ˆ 2, zVar Var Cov         and       12 1 2 1 2 ˆˆ ˆ ˆˆ ˆ ˆ ˆ ˆ 2, zVar Var Cov         respectively, where 11 10 1 11 10 11 10 ˆxx x x y y      , 11 01 2 11 01 11 01 ˆxx x x y y      , 01 00 1 01 00 01 00 ˆyy x x y y      and 10 00 2 10 00 10 00 ˆyy x x y y      . Both test statistics follow a standard normal distribution, and the variances are estimated by applying the delta method (the expressions are shown in the method of Roldán-Nofuentes et al [7] which will now be summarized). 38 Method of Kosinski Kosinski [5] proposed a weighted generalized score statistic to solve the hypothesis test of comparison of the PVs. The weighted generalized score statistic for the test 0 1 2 :H   is     12 10 11 01 11 ˆˆ 11 ˆˆ 12 p p p z Cn n n n             , and the weighted generalized score statistic for the test 0 1 2 :H   is     12 00 01 00 10 ˆˆ 11 ˆˆ 12 p p p z Cn n n n             , which has a standard normal distribution when the null hypothesis is true, where 11 10 01 11 10 01 2 ˆ2 p x x x n n n    and 00 01 10 00 01 10 2 ˆ2 p y y y n n n    are the pooled positive PV and pooled negative PV respectively, and   22 11 11 11 10 01 ˆˆ 1 2 pp p xy Cn n n    ,   2 2 00 00 00 01 10 ˆˆ 1 2 pp p xy Cn n n    and ij ij ij n x y . Method of Roldán-Nofuentes et al Roldán-Nofuentes et al [7] studied the simultaneous comparison of the PVs of two BDTs subject to a paired design. The simultaneous comparison of the PVs of two binary tests consists of solving the hypothesis test   0 1 2 1 2 : and H      vs   1 1 2 1 2 : and/or H      , Applying the delta method, the estimated variances-covariances of the estimators of the PVs are: 39        10 11 10 11 13 10 11 10 11 ˆˆx x y y Var n x x y y      ,        00 01 00 01 13 00 01 00 01 ˆˆx x y y Var n x x y y             01 11 01 11 23 01 11 01 11 ˆˆx x y y Var n x x y y      ,        00 10 00 10 23 00 10 00 10 ˆˆx x y y Var n x x y y      ,           01 10 11 11 01 10 11 11 01 10 11 10 11 12 22 01 11 01 11 10 11 10 11 ˆˆˆ ,x x y x y y y y x x x y y Cov x x y y x s y y                 ,             00 10 11 10 10 10 10 11 00 10 10 00 10 11 12 22 00 10 00 10 10 11 10 11 ˆˆ ˆ,x x x y x y x x y r x y r y Cov x x y y x x y y                ,             00 01 11 01 01 01 01 11 00 01 01 00 01 11 21 22 00 01 00 01 01 11 01 11 ˆˆ ˆ,x x x y x y x x y y x y y y Cov x x y y x x y y                ,           2 00 00 01 10 00 00 01 10 00 01 10 00 01 12 22 00 01 00 01 00 10 00 10 ˆˆˆ ,x y y y y y x x x x x y y Cov x x y y x x y y                 ,     1 1 2 2 ˆˆ ˆˆ ˆˆ , , 0Cov Cov      . The test statistic for   0 1 2 1 2 : and H      is   1 2ˆ ˆˆ T T T Q η A A A Aη , where   1 1 1 2 ˆ ˆ ˆ ˆˆ , , , T     η , ˆ  is the estimated variance-covariance matrix of ˆ η and A is the design matrix, i.e. 1 0 1 0 0 1 0 1      A . The test statistic 2 Q is distributed asymptotically according to a central chi-square distribution with two degrees of freedom if 0 H is true. When all of these methods are used applying multiple imputation, all the equations are valid for the mth complete dataset, adding superindex   m to all of the terms of the equations. 40 Appendix B Probabilities   12 , , 1 ij P T i T i D     and   12 , , 0 ij φ P T i T i D    are written in terms of the PVs as:                                  1 1 2 1 2 2 1 1 2 1 2 11 10 1 2 1 2 2 2 1 1 1 1 01 12 1 1 2 2 1 1 2 2 2 2 1 00 12 1 1 2 0 2 2 0 1 2 1 2 11 10 12 , , , 1 1 1 1 , 11 11 , q pY q qq pYY pYY q pY q pYY q p p p Y pYY q qY q qq qYY                                                                                            12 2 2 1 0 1 1 01 12 2 0 1 2 1 2 1 2 00 12 2 1 1 1 1 1 2 1 2 12 , 11 , 11 1 2 1 . qYY q qY q qYY q q q YY qYY p p p Y qYY                                             Appendix C The maximum likelihood estimator of PVs in the presence of partial verification of the disease are [10, 11] 111 10 10 11 1 1 1 ˆjj jjj na n n a b     and 100 10 00 01 0 0 1 ˆjj jjj nb n n a b     for test 1, and 111 20 01 11 1 1 1 ˆii iii na n n a b     and 100 20 00 10 0 0 1 ˆii iii nb n n a b     for test 2. From the table of complete data (Table 1b), the log-likelihood function based on n individuals is 41           11 , 0 , 0 log log . ij ij ij ij ij ij ij i j i j l a d b c d         θ From this function it is obtained that ij ij ij ad ˆn    and ij ij ij ij b c d ˆn    . In order to demonstrate that the EM algorithm converges to the ML estimators, we are going to follow the same steps as Little and Rubin [7]. With the EM algorithm, the estimator of 1  is calculated as             11 11 1 1 1 1 1 1 00 10 11 10 11 11 ˆ 11 ˆˆˆ k kk j j j j j kk jj jj a d a c n n n n                . Therefore it is necessary to show that     11 11 0 k jj j ad     converges to 111 011 jj jjj na ab   . Taking     111 1 1 1 ˆ ˆ ˆ kk jj j j j ad n        ,     11 1 1 1 1 1 ˆ ˆ ˆ kk j j j j j j b c d n        and     11 1 1 1 1 11 ˆ ˆˆ kk j j j j j jj d d d c        , then it is obtained that 11 1 11 jj j jj ac dab  . Then     1 1 1 1 11 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 kj j j j j j j j j j j j j j j j j j j j a c a c n a a d a a a b n c a b                                   . Therefore,   1 1 ˆk   converges to 1 ˆ  . The convergence of the rest of the estimators of PVs is demonstrated in a similar way.