Mann-Whitney type tests for microarray experiments: the R package gMWT

JSS Journal of Statistical Software June 2015, Volume 65, Issue 9. http://www.jstatsoft.org/ Mann-Whitney Type Tests for Microarray Experiments: The RPackage gMWT Daniel Fischer University of Tampere Hannu Oja University of Turku Abstract We present the Rpackage gMWT which is designed for the comparison of several treatments (or groups) for a large number of variables. The comparisons are made using certain probabilistic indices (PI). The PIs computed here tell how often pairs or triples of observations coming from different groups appear in a specific order of magnitude. Classical two and several sample rank test statistics such as the Mann-Whitney-Wilcoxon, Kruskal-Wallis, or Jonckheere-Terpstra test statistics are simple functions of these PI. Also new test statistics for directional alternatives are provided. The package gMWT can be used to calculate the variable-wise PI estimates, to illustrate their multivariate distribution and mutual dependence with joint scatterplot matrices, and to construct several classical and new rank tests based on the PIs. The aim of the paper is first to briefly explain the theory that is necessary to understand the behavior of the estimated PIs and the rank tests based on them. Second, the use of the package is described and illustrated with simulated and real data examples. It is stressed that the package provides a new flexible toolbox to analyze large gene or microRNA expression data sets, collected on microarrays or by other high-throughput technologies. The testing procedures can be used in an eQTL analysis, for example, as implemented in the package GeneticTools. Keywords: eQTL, Jonckheere-Terpstra test, Kruskal-Wallis test, Mann-Whitney test, permutation test, several samples, simultaneous testing, union-intersection test, U-statistic. 1. Introduction We consider nonparametric tests used in the analysis of gene or microRNA expression data sets with several treatments (groups). For each separate expression variable, the null hypothesis to be tested is that there is no difference between the distributions of the expression in different groups. To avoid strong (parametric) distributional assumptions, the alternatives are formulated using probabilities that pairs or triples of observations coming from different 2gMWT: Generalized Mann-Whitney Tests in R groups are in a specific order of magnitude. The interesting probabilities are called probabilistic indices (PI), see also Thas, De Neve, Clement, and Ottoy (2012). The test statistics are based on natural estimates of these PIs, that is, the corresponding two and several sample U-statistics. Classical several-sample rank test statistics such as the Kruskal-Wallis or Jonckheere-Terpstra test are special cases in this approach. Also, as the number of variables (microRNAs) is typically huge and the test statistics for different variables are dependent, we face a serious simultaneous testing problem. See Fischer, Oja, Sen, Schleutker, and Wahlfors (2014) for more details. The package gMWT (Fischer and Oja 2015) provides nonparametric tools for the comparison of several groups/ treatments when the number of variables is large, and is available from the Comprehensive RArchive Network (CRAN) at http://CRAN.R-project.org/package= gMWT. The tools are the following. (i) Computation of the PI estimates for the group comparisons. The probabilistic indices here are (a) the probability Ptt0that a random observation from group tis smaller than a random observation from group t0, and (b) the probability Ptt0t00 that observations from groups t, t0, t00 appear in this same order. The tools are also given to produce the plots of variable-wise PIs. (ii) Computation of the pvalues of some classical and some new nonparametric tests for the comparison of several groups/treatments. The tests are based on the use of the probabilistic indices Ptt0and Ptt0t00 . Classical Mann-Whitney-Wilcoxon, Kruskal-Wallis and Jonckheere-Terpstra tests are included. (iii) Tools for the simultaneous testing problem. As the package is meant for the analysis of gene expression data for example, tools to control the family-wise error rate and/or the false discovery rate are provided as plots for expected versus observed rejected null hypotheses with the Simes (improved Bonferroni) and Benjamini-Hochberg rejection lines. A list of rejected null hypotheses may be obtained as well. Some standard nonparametric methods such as the Mann-Whitney and Kruskal-Wallis tests have been implemented in the Rstats (RCore Team 2014) package. Linear rank statistics for the two and several sample location problems with ordered and unordered alternatives have been implemented also in the coin (Hothorn, Hornik, van de Wiel, and Zeileis 2008) package. Exact and permutation versions of the Jonckheere-Terpstra test are given by the package clinfun (Seshan 2014); this function is used in our package as the second option in our implementation. One contribution of our package gMWT is that these and several other nonparametric tests are collected with the same syntax under the same roof with a simultaneous testing possibility for several variables. The interfaces of the functions are tailored for large datasets with many groups and several variables, so that the application and comparisons of competing testing procedures are easier. With scatterplot matrices for the relevant PIs, it is also possible to illustrate and understand the joint variable-wise behavior of the standard tests. The structure of this paper is as follows. After a brief review of the theory in Section 2 we present some practical solutions in Section 3for the computation of the PIs and the permutational pvalues of the corresponding tests. In Section 4a general description of the package gMWT is given with a typical workflow for its use. We also discuss the calculation of Journal of Statistical Software 3 the PIs and their scatterplot matrices, and it is described how the tests are performed. Also, the tools for the multiple testing problem are described. In Section 5, the use of the package is illustrated with a simulated data set as well as with real genotype data. In the latter case, an expression quantitative trait locus (eQTL) analysis is performed with the packages gMWT and GeneticTools (Fischer 2014). 2. Statistical inference based on probabilistic indices 2.1. Null hypothesis and alternatives based on Ptt0and Ptt0t00 Consider first the univariate case and the comparison of Tgroups. Let xt1, . . . , xtNtbe a random sample from a distribution with cumulative distribution function Ft,t= 1, . . . , T, and let the samples be independent. The total sample size is then N=N1+· · · +NT. We wish to test the null hypothesis H0:F1=F2=· · · =FT. The interesting alternatives are formulated using certain probabilistic indices. As ties may often be present, we write I(x, y) = I(x<y) + 1 2I(x=y) and I(x, y, z) = I(x<y<z) + 1 2I(x=y < z) + 1 2I(x<y=z) + 1 6I(x=y=z), with I(·) being the indicator function, which is 1 if the argument (·) is true and 0 else. The interesting alternatives are then given in terms of the probabilities Ptt0=E(I(xt, xt0)) and Ptt0t00 =E(I(xt, xt0, xt00 )). Note that, as I(x, y) = I(x, y, z) + I(x, z, y) + I(z, x, y) the probabilities satisfy Ptt0=Ptt0t00 +Ptt00 t0+Pt00 tt0. Under the null hypothesis H0:F1=F2=· · · =FT, for all t, t0, t00, Ptt0=1 2and Ptt0t00 =1 6. We say that F1and F2are stochastically ordered and write F1st F2if F1(x)≥F2(x)∀x∈R. Then Ftst Ft0⇒Ptt0≥1 2 and Ftst Ft0st Ft00 ⇒Ptt0t00 ≥1 6 but the converse statements are not true. 4gMWT: Generalized Mann-Whitney Tests in R In the comparison of T= 3 treatments interesting alternatives might then be formulated, for example, as H1:P12 6=1 2or P13 6=1 2or P23 6=1 2, or H1:P12 ≥1 2or P13 ≥1 2or P23 ≥1 2with at least one strict inequality, or H1:P13 ≥1 2or P23 ≥1 2with at least one strict inequality, or H1:P123 >1 6. The tests will then be based on the estimates ˆ P12,ˆ P13,ˆ P23 and ˆ P123 and should be constructed keeping the interesting alternative in mind. 2.2. Estimation of Ptt0and Ptt0t00 The probabilities Ptt0and Ptt0t00 are naturally estimated by corresponding U-statistics ˆ Ptt0=1 NtNt0 Nt X i=1 Nt0 X i0=1 I(xti, xt0i0) and ˆ Ptt0t00 =1 NtNt0Nt00 Nt X i=1 Nt0 X i0=1 Nt00 X i00 =1 I(xti, xt0i0, xt00 i00 ). A natural statistic for the comparison between group tand other groups is ˆ Pt=1 N−NtX t06=t Nt0ˆ Ptt0. A general several-sample U-statistic theory can be used to find the (joint) limiting properties of ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 under the null hypothesis. See, e.g., Chapter 5 in Serfling (1980). 2.3. Tests based on estimates ˆ Ptt0and ˆ Ptt0t00 As seen before, we have a hierarchy nˆ Ptt0t00 o→nˆ Ptt0o→nˆ Pto and one can construct test statistics at different levels of this hierarchy. Some choices are the following. 1. Use test statistics ˆ Ptt0t00 for H0:Ft=Ft0=Ft00 vs. H1:Ptt0t00 6=1 6. Of course one-sided alternatives are possible as well. Journal of Statistical Software 5 2. Use the Mann-Whitney (MW) test statistics ˆ Ptt0for H0:Ft=Ft0vs. H1:Ptt06=1 2and the Jonckheere-Terpstra (JT) test statistics for H0:F1=· · · =FTvs. H1:Ptt0≥1 2 for all t<t0with at least one strict inequality. Note that F1st F2st · · · st FT with at least one strict inequality implies the latter H1. We have two versions of JT test statistic, namely, JT =X t<t0 NtNt0ˆ Ptt0and JT ∗=X t<t0 ˆ Ptt0. 3. For a fixed group t, use a Mann-Whitney test statistic ˆ Ptfor H0:F1=· · · =FTvs. H1:F1=... =Ft−1=Ft+1 =... =FT6=Ft. Use the Kruskal-Wallis test statistic KW =12 N(N+ 1) T X t=1 (ˆ Pt−Nt(N−Nt)/2)2 Nt for H0:F1=· · · =FTvs. H1:Ft6=Ft0for at least one pair t, t0. The alternative then implies that Ptt06=1 2for at least one pair t, t0. 4. Use a union-intersection test (UIT) to compare three groups t,t0, and t00. The test statistic is a combination of statistics ˆ Ptt00 and ˆ Pt0t00 and is meant for the alternative max(Ptt00 , Pt0t00 )>1 2. The test statistic can be found in Appendix A, see also Fischer et al. (2014) for the details. 3. Computational solutions 3.1. Fast computation of ˆ Ptt0and ˆ Ptt0t00 Consider the univariate case and write the N-vector x= (x1, x2, . . . , xN)>= (x11, . . . , x1N1, x21, . . . , x2N2, . . . , xT1, . . . , xT NT)> of observations coming from all Tgroups. The PIs are based on two N×Nmatrices Ist = Ist(x) and Ieq =Ieq(x) with the elements (Ist(x))ij =I(xi< xj) and (Ieq(x))ij =I(xi=xj), i, j = 1, . . . , N. The matrices Ist =Ist(x) and Ieq =Ieq(x) can then be decomposed as Ist =    Ist 11 Ist 12 . . . Ist 1T Ist 21 Ist 22 . . . Ist 2T . . . . . . . . . . . . Ist T1Ist T2. . . Ist T T     and Ieq =    Ieq 11 Ieq 12 . . . Ieq 1T Ieq 21 Ieq 22 . . . Ieq 2T . . . . . . . . . . . . Ieq T1Ieq T2. . . Ieq T T     , where the Ni×Njsubmatrices Ist ij and Ieq ij compare treatments iand j,i, j = 1, . . . , T. Then ˆ Ptt0=1 NtNt0 1> NtIst tt0+1 2Ieq tt01Nt0, 6gMWT: Generalized Mann-Whitney Tests in R 0 100 200 300 400 500 Simulated sample size Time in seconds 30 150 300 450 600 750 900 1050 ● ● ● ● Naive,R Submatrices,R Submatrices,C++ Naive,C++ Figure 1: Computation times for the pvalues using the permutation test version with test statistic ˆ Ptt0t00 . Three groups with equal group sizes were used, and the number of permutations in each case was 2000. Naive and submatrix approaches implemented in Rand C++ are compared. where 1kis the notation for a k-vector full of ones. For the triples we get ˆ Ptt0t00 =1 NtNt0Nt00 1> NtIst tt0Ist t0t00 +1 2Ieq tt0Ist t0t00 +1 2Ist tt0Ieq t0t00 +1 6Ieq tt0Ieq t0t00 1Nt00 . In case that no ties are present, the matrix Ieq is simply a zero matrix. Using a naive implementation, we would calculate the probabilities ˆ Ptt0t00 one by one while going through all NtNt0Nt00 triple comparisons. In our submatrix approach we thus calculate the matrices Ist and Ieq only once and then use the submatrices to find the probabilities ˆ Ptt0t00 . This leads to improved calculation times especially for permutation versions of the tests. For the computation time comparisons in Rand C++ (via Rcpp,Eddelbuettel and Fran¸cois 2011, and RcppArmadillo,Eddelbuettel and Sanderson 2014), see Figure 1. 3.2. Computation of pvalues If the null hypothesis H0:F1=· · · =FTis true, then Px ∼xfor all N×Npermutation matrices P. By Px ∼xwe mean that the distributions of Px and xare the same. Matrix P is an N×Npermutation matrix if it is obtained from an identity matrix INby permuting its rows and/or columns. The number of distinct permutation matrices is N!. Note that our test statistics are functions of the matrices Ist and Ieq and that Ist(Px) = PIst(x)P>and Ieq(Px) = PIeq(x)P>. The exact pvalue from a permutation test using a test statistic Q=Q(Ist,Ieq) is then PnQ(PIstP>,PIeqP>)≥Q(Ist,Ieq)o, Journal of Statistical Software 7 where the probability is taken over N! equally probable values of P. The pvalue can in practice be estimated by 1 M M X m=1 InQ(PmIstP> m,PmIeqP> m)≥Q(Ist,Ieq)o, where P1,...,PMis a random sample from a uniform distribution over the set of N×N permutation matrices. Naturally, the larger M, the better is the estimate of the exact pvalue. Approximate pvalues may also be based on the limiting joint normality of the estimates (U-statistics) ˆ Pt,ˆ Ptt0, and ˆ Ptt0t00 . If no ties are present, the test statistics based on the PIs are strictly distribution-free with limiting variances and covariances that are easily found. 4. The package gMWT 4.1. General features The Rpackage gMWT can be used to calculate the variable-wise probabilistic indices ˆ Pt,ˆ Ptt0, and ˆ Ptt0t00 , to illustrate their joint distributions and dependence with scatterplot matrices, and to perform various rank tests based on ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 described in Section 2.3. See Figure 2for possible workflows. A practical application for the testing procedures is an eQTL analysis for a combined analysis of microarray and genotype data. In Section 5.2 we illustrate the use of the packages GeneticTools and gMWT with the directional triple test for testing for eQTL. In the following, the input matrix Xis a data matrix with observations as rows and variables as columns. The vector gindicates the group membership; its length is then the number of rows in X. 4.2. Computation of ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 The estimated PIs, ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 , are calculated using the command estPI(X, g, type = "pair", goi, mc = 1, order = TRUE) The options type = "single","pair" or "triple" specify the PIs to be computed, that is, ˆ Pt,ˆ Ptt0or ˆ Ptt0t00 . The vector goi (”groups of interest”) specifies the groups (values of g) to be used in the comparisons. The option order specifies whether the PIs should be calculated for all possible pairs and triples (order = FALSE) or just for pairs and triples with increasing group labels. In the four group case, for example, an estPI call with the (default) options (type = "pair", order = TRUE) would calculate the estimated PIs ˆ P12,ˆ P13,ˆ P14,ˆ P23,ˆ P24,ˆ P34 and in case of using the parameters (type = "triple", order = TRUE) the estimated PIs ˆ P123,ˆ P124,ˆ P134,ˆ P234. For matrix valued X, the option mc can be used to execute the parallel calculation on mc-many cores in order to speed up the calculation (available only for Linux systems). The result of the function estPI is a list, containing a matrix probs with the PIs as rows and variables as columns. The other list items are the used parameters. 8gMWT: Generalized Mann-Whitney Tests in R Data NxG G Variables: g=1,...,G N obs. n=1,...,N Groups t=1,...T single pairs triple Calculate Probabilistic Index Function: estPI Apply test Function: gmw KW MW: 1 vs. All MW: pairs UIT JT JT* Triple Plot the test results Function: plot plotProb plotProb Figure 2: Possible calculation workflows. Journal of Statistical Software 9 4.3. Scatterplots for ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 Based on the probabilities calculated via estPI the package creates diagnostic scatterplot matrices for variable-wise PIs. The command is either plotPI(X, g, col, zoom = FALSE, highlight = NULL, hlCol = "red") if the diagnostic plots are produced directly from the data or pe <- estPI(X, g, type = "pair", goi) plot(pe, col, zoom = FALSE, highlight, hlCol) if the probabilities are first calculated by estPI. The function plotPI naturally allows also all options used in the estPI call. The type of the plots depends on the selected PIs; the pairwise scatterplots are then for all possible pairs of PIs. Additional plotting options are highlight,hlCol and zoom. Using the highlight option, indicated variables are plotted with a color specified in hlCol. If the Boolean flag zoomed is set, the plots are zoomed to the active area of the PIs. Without this flag, the bivariate plots are in [0,1] ×[0,1]. 4.4. Tests based on ˆ Pt,ˆ Ptt0and ˆ Ptt0t00 The basic call for applying the testing procedure is gmw(X, g, goi, test = "mw", type = "permutation", prob = "pair", nper = 2000, alternative = "two.sided", mc = 1, output = "min", keepPM = FALSE, mwAkw = FALSE) where only input variables Xand gare compulsory. If Xis a matrix then the chosen test is applied variable-wise to the data and the results are reported as a matrix of pvalues. Specifying the option output = "full" leads to a more detailed list output with the same length as the number of different alternatives are tested and each list item contains then again a list with as many columns as there are in X. Each entry in the encapsulated list is a test result of class ‘htest’. The vector ggives the group numbers in a natural order. The groups used in the analysis can be specified via goi. If no goi is specified, all groups are used. With the option test, the test ("uit", "triple", "mw", "kw", "jt", "jt*") can be specified. The option prob may be used only in the case of the Mann-Whitney test. For all pairwise group comparisons one uses (prob = "pair") while option (prob = "single") compares each group to the rest of the data. The option type specifies whether the permutation type test ("permutation") or the asymptotical test version ("asymptotical") is used. For some standard tests the procedures from base Rare available, and the Jonckheere-Terpstra test is implemented in the clinfun package. The option type = "external" allows the use of these test versions. A permutation type of test is available for all tests but, as the package is still under active development, the asymptotical versions are not yet available for all cases. The number of permutations is selected with the option nper. Different alternatives are available whenever they are natural and are set with alternative = "smaller","greater" and "two.sided". 16 gMWT: Generalized Mann-Whitney Tests in R Gene1 − 1 Chr 1 : 8474334 − 9574334 Chromosomal Position in MB p−value 8.5 8.72 8.94 9.16 9.38 9.6 0 0.2 0.4 0.6 0.8 1 Mono. NA ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● Figure 9: eQTL-plot. The dashed line symbolizes the location of the linked gene and each dot represents the test result for a SNP test. The dotted line refers to the chosen significance level. In case that all individuals have the same genotype no test is performed and this is marked as “Mono.”. If no genotype information was available for a SNP, this is marked as NA. The expression matrix is specified with the gex option. The option xAnnot is important for the gene annotations and is a list where each list item refers to a gene from the columns of geneEx. The names have to match and the eQTL is only performed for matching pairs. In case that no annotation xAnnot is given to the function no window is used and all combinations are considered instead. Be aware that this might lead to a very long lasting calculation. Each list item in geneAnnotations is a matrix like R> geneAnnotations $Gene1 Chr Start End 1 1 8974334 9074334 $Gene2 Chr Start End 1 12 135633062 135738062 2 12 135735062 135838062 Journal of Statistical Software 17 This takes into account that certain probes have multiple locations in the genome and our method will test for all those locations. In case that the labeling of individuals between the expression data and the genotype data differs, there is an option to give a new vector of labels, called genoSamples. Here can new labels for the individuals in the genotype data be specified, that match with the row names of geneEx. It is important to note that the order of the rows and columns of all lists and matrices does not have to match. The function takes the smallest subsets of individuals and genes and takes then those SNPs, which are in a window around that gene. The window size can be specified with the option windowSize using the unit megabases (MB). After the eQTL is performed we can visualize the results with R> plot(myEQTL, which, file, sig) The which option specifies for which genes from geneEx the plots shall be created. If no option is given, then all plots are created. Because this might lead to a vast number of pictures, there is also an option file to specify a file name such that the plots are saved in a file with this name. The output of a single plot can be seen in Figure 9with a chosen significance level sig = 0.1. This way we can see for interesting genes the behavior of the surrounding SNPs onto the gene expression. For small datasets it is also possible to check visually for interesting genes. For larger gene sets the function extractEQTL can be used to determine a set of interesting genes. We applied this function on real data e.g., in Siltanen et al. (2013). References Benjamini Y, Hochberg Y (1995). “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society B,57(1), 289–300. Clayton D (2014). snpStats:SnpMatrix and XSnpMatrix Classes and Methods.Rpackage version 1.16.0, URL http://www-gene.cimr.cam.ac.uk/clayton/. Eddelbuettel D, Fran¸cois R (2011). “Rcpp: Seamless Rand C++ Integration.” Journal of Statistical Software,40(8), 1–18. URL http://www.jstatsoft.org/v40/i08/. Eddelbuettel D, Sanderson C (2014). “RcppArmadillo: Accelerating Rwith High-Performance C++ Linear Algebra.” Computational Statistics & Data Analysis,71, 1054–1063. Fischer D (2014). GeneticTools: Collection of Genetic Data Analysis Tools.Rpackage version 0.3, URL http://CRAN.R-project.org/package=GeneticTools. Fischer D, Oja H (2015). gMWT: Generalized Mann-Whitney Type Tests.Rpackage version 1.0, URL http://CRAN.R-project.org/package=gMWT. Fischer D, Oja H, Sen PK, Schleutker J, Wahlfors T (2014). “Generalized Mann-Whitney Type Tests for Microarray Experiments.” Scandinavian Journal of Statistics,41(3), 672–692. 18 gMWT: Generalized Mann-Whitney Tests in R Hothorn T, Hornik K, van de Wiel MA, Zeileis A (2008). “Implementing a Class of Permutation Tests: The coin Package.” Journal of Statistical Software,28(8), 1–23. URL http://www.jstatsoft.org/v28/i08/. Perlman MD (1969). “One-Sided Testing Problems in Multivariate Analysis.” The Annals of Mathematical Statistics,40(2), 549–567. RCore Team (2014). R: A Language and Environment for Statistical Computing.RFoundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/. Serfling RJ (1980). Approximation Theorems of Mathematical Statistics. John Wiley & Sons. Seshan VE (2014). clinfun: Clinical Trial Design and Data Analysis Functions.Rpackage version 1.0.6, URL http://CRAN.R-project.org/package=clinfun. Siltanen S, Fischer D, Rantapero T, Laitinen V, Mpindi JP, Kallioniemi O, Wahlfors T, Schleutker J (2013). “ARLTS1 and Prostate Cancer Risk – Analysis of Expression and Regulation.” PLoS ONE,8(8), e72040. Simes RJ (1986). “An Improved Bonferroni Procedure for Multiple Tests of Significance.” Biometrika,73(3), 751–754. Thas O, De Neve J, Clement L, Ottoy JP (2012). “Probabilistic Index Models.” Journal of the Royal Statistical Society B,74(4), 623–671. Westfall PH, Young SS (1993). Resampling-Based Multiple Testing: Examples and Methods for p-Value Adjustment. John Wiley & Sons. Journal of Statistical Software 19 A. Union-intersection test for P13 and P23 We implemented two ways to calculate pvalues from the UIT, the first one is based on a permutation approach and the other one is based on asymptotical results. In both cases we first need to calculate the critical value cof the test statistic Q∗. For the test statistics S1=ˆ P13 and S2=ˆ P23,S= (S1, S2)>∼N2(µ,Σ) approximately with Σ=1ρ ρ1. Our UIT test statistic Q∗is then Q∗=I0·0 + I1·(S2−ρS1)2 1−ρ2+I2·(S1−ρS2)2 1−ρ2+I3·S>Σ−1S, with I0=I(S1≤ρS2, S2≤ρS1),I1=I(S1<0, S2> ρS1), I2=I(S1> ρS2, S2<0),I3=I(S1≥0, S2≥0). As I0+I1+I2+I3= 1, at most one of the three test statistics in the sum contributes to Q∗. An approximate pvalue is then given by the approximation (Perlman 1969) P(Q∗> c) = 1 2P(χ2 1> c) + cos−1ρ 2πP(χ2 2> c). For a permutation test version, we permute the elements of the vector of the group variable Mtimes with resulting values of test statistics Q∗ 1, Q∗ 2, . . . , Q∗ M; the approximate pvalue is then p=1 M M X m=1 I(Q∗ m≥Q∗). Affiliation: Daniel Fischer School of Health Sciences University of Tampere 33014 Tampere, Finland E-mail: Daniel.Fische[email protected] Hannu Oja Department of Mathematics and Statistics University of Turku 20014 Turku, Finland Journal of Statistical Software http://www.jstatsoft.org/ published by the American Statistical Association http://www.amstat.org/ Volume 65, Issue 9 Submitted: 2012-05-30 June 2015 Accepted: 2014-08-06