Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Zhang, Cassie Yue
Master of Science - MSc
University of British Columbia
This thesis addresses a significant problem in numerous scientific fields - the challenge of determining whether two sets of data are statistically identical, a procedure known as a two-sample homogeneity test. Various methods, including well-known parametric tests like t-tests and analysis of variance (ANOVA), have been employed for two-sample homogeneity tests. However, these tests have certain limitations. They rely heavily on specific assumptions about the data's distribution. If these assumptions are violated, the results can be misleading. To overcome these limitations, non-parametric methods that do not require strict assumptions about the data's nature have been explored. However, these too have their shortcomings. For instance, techniques like the Wilcoxon rank-sum test can only be applied effectively in univariate cases, as multivariate data lacks a natural ranking. This complicates the analysis and limits the utility of such methods for multivariate data. In response to identified challenges, this work introduces an innovative approach to rank multivariate data through the concept of data depth. It encompasses a thorough comparison of eleven widely-used depth functions and the development of two novel two-sample homogeneity tests, significantly enhancing the accuracy of homogeneity assessments in multivariate datasets. Further, inspired by the quality index [25] used in quality assurance contexts to measure a population's overall 'outlyingness' relative to another, the proposed test statistics are designed to yield a 𝒳²(1) asymptotic null distribution. Simulation studies underscore the robust performance of these tests, and the discussion also extends into their potential application in multivariate multisample contexts, demonstrating the versatility and broad applicability of depth-based methodologies in statistical analysis. The increasing number of samples correspondingly increases the number of pairwise quality indexes, an aspect considered and explored as part of the study. Future works include application on high-dimensional data and change point analysis. With its focus on boosting the power of depth-based tests, this thesis marks a significant step forward in the field of data analysis.
Ranking Functions of DatabyCassie Yue ZhangB.Sc., Thompson Rivers University, 2022M.Sc., The University of British Columbia, 2023A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)November 2023© Cassie Yue Zhang 2023The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Ranking Functions of Datasubmitted by Cassie Yue Zhang in partial fulfilment of the requirementsof the degree of Master of ScienceDr. Xiaoping Shi, Department of Computer Science, Mathematics, Physics and StatisticsSupervisorDr. John Braun, Department of Computer Science, Mathematics, Physics and StatisticsSupervisory Committee MemberDr. Lengyi Han, Department of Computer Science, Mathematics, Physics and StatisticsSupervisory Committee MemberDr. Julian Cheng, School of EngineeringUniversity ExamineriiAbstractThis thesis addresses a significant problem in numerous scientific fields - thechallenge of determining whether two sets of data are statistically identical,a procedure known as a two-sample homogeneity test.Various methods, including well-known parametric tests like t-tests andanalysis of variance (ANOVA), have been employed for two-sample homo-geneity tests. However, these tests have certain limitations. They rely heav-ily on specific assumptions about the data’s distribution. If these assump-tions are violated, the results can be misleading.To overcome these limitations, non-parametric methods that do not re-quire strict assumptions about the data’s nature have been explored. How-ever, these too have their shortcomings. For instance, techniques like theWilcoxon rank-sum test can only be applied effectively in univariate cases,as multivariate data lacks a natural ranking. This complicates the analysisand limits the utility of such methods for multivariate data.In response to identified challenges, this work introduces an innovativeapproach to rank multivariate data through the concept of data depth. Itencompasses a thorough comparison of eleven widely-used depth functionsand the development of two novel two-sample homogeneity tests, signifi-cantly enhancing the accuracy of homogeneity assessments in multivariatedatasets. Further, inspired by the quality index [25] used in quality as-surance contexts to measure a population’s overall ‘outlyingness’ relative toanother, the proposed test statistics are designed to yield a χ2(1) asymptoticnull distribution. Simulation studies underscore the robust performance ofthese tests, and the discussion also extends into their potential application inmultivariate multisample contexts, demonstrating the versatility and broadapplicability of depth-based methodologies in statistical analysis.The increasing number of samples correspondingly increases the numberof pairwise quality indexes, an aspect considered and explored as part ofthe study. Future works include application on high-dimensional data andchange point analysis. With its focus on boosting the power of depth-basedtests, this thesis marks a significant step forward in the field of data analysis.iiLay SummaryThis research enhances our ability to compare two sets of data and test ifthey’re statistically identical. Traditional techniques, such as t-tests andANOVA, may be inaccurate when certain data assumptions aren’t met.While non-parametric methods avoid this issue, but are mostly limited tounivariate data. This work introduces ’data depth’, a concept allowing mul-tivariate data to be ranked. This thesis presents a new method based ondata depth, proving its effectiveness through simulation studies and real-lifeexamples. This research examines existing methods, introduces data depth,and explores this new technique, making it easier for more accurate datacomparisons across various disciplines.iiiPrefaceThis thesis encapsulates rigorous research and discussions undertaken lastsummer with Dr. Xiaoping Shi and Dr. Cindy Fu. A substantial segment ofthis work, especially the content in Chapter 6, has been formally published.I am proud to have co-authored the following paper, which presents keyaspects of the research found within this thesis:• Shi, X., Zhang, Y., & Fu, Y. (2023). Two-Sample Tests Based on DataDepth. Entropy, 25(2), 238. [1]In my role as a co-author, I contributed extensively to the manuscript,including the drafting of the code for numerical simulations. The powercomparison plots presented in the aforementioned chapter were developedby Dr. Xiaoping Shi.The findings displayed here are the result of many insightful group dis-cussions. Particularly, our brainstorming over the last summer was crucialin formulating the new test statistics that are the foundation of this work. Iundertook the task of running numerical simulations and gathering multiplereal data sets to test and validate our theories.While writing this thesis, I had the opportunity to present our work ata poster session during the “Statistics Conference in Genomics, Pharmaceu-tical Science, and Health Data Science,” held at the University of Victo-ria. Additionally, I gave an oral presentation at the conference on “RecentResearch Progress in Advanced Statistical Methods” at Squamish. Thesepresentations led to valuable feedback and the generation of new ideas. Ialso want to highlight that the work conducted for this thesis did not involveany human or animal subjects, and therefore did not require ethics approvalfrom the UBC Research Ethics Board.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Parametric test of the multivariate homogeneity tests . . 43.1 Introduction to MANOVA . . . . . . . . . . . . . . . . . . . 43.2 Limitations and challenges . . . . . . . . . . . . . . . . . . . 63.3 Robust MANOVA . . . . . . . . . . . . . . . . . . . . . . . . 64 Non-parametric univariate homogeneity tests . . . . . . . . 94.1 Introduction to the Wilcoxon Rank-Sum Tests . . . . . . . . 94.2 Extensions of Wilcoxon rank-sum test . . . . . . . . . . . . . 104.3 Real data analysis . . . . . . . . . . . . . . . . . . . . . . . . 115 Non-parametric multivariate homogeneity tests . . . . . . 145.1 Depth functions and one-sample centrality test . . . . . . . . 145.2 Selection of depth functions [65] . . . . . . . . . . . . . . . . 415.2.1 Choosing statistical depth measures . . . . . . . . . . 415.2.2 Computational feasibility . . . . . . . . . . . . . . . . 48vTable of Contents5.3 Depth-based two-samples homogeneity test . . . . . . . . . . 505.3.1 Depth-based Kruskal-Wallis-type statistics H [46] . . 505.3.2 Depth-based quality index Q [25] . . . . . . . . . . . 556 A novel depth-based pairwise quality index . . . . . . . . . 586.1 Two-sample pairwise quality indexes Wm,n(w1, w2) and Mm,n 586.2 Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . 606.2.1 Two-sample cases with the Mahalanobis, spatial andprojection depth [25] . . . . . . . . . . . . . . . . . . 606.2.2 Two sample power comparison with other depths . . 666.2.3 Multi-sample comparison . . . . . . . . . . . . . . . . 716.2.4 High-dimension cases . . . . . . . . . . . . . . . . . . 736.3 Real data analysis . . . . . . . . . . . . . . . . . . . . . . . . 756.3.1 Iris data . . . . . . . . . . . . . . . . . . . . . . . . . 756.3.2 Skull data . . . . . . . . . . . . . . . . . . . . . . . . 766.3.3 leukemia data . . . . . . . . . . . . . . . . . . . . . . 787 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82AppendicesA R code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90A.1 Parametric test of the multivariate homogeneity tests . . . . 90A.1.1 Table 3.2 Classical and robust MANOVA comparison 90A.2 Non-parametric univariate homogeneity tests . . . . . . . . . 91A.2.1 Sec 4.3 Wilcoxon rank-sum test real data analysis . . 91A.3 Non-parametric multivariate homogeneity tests . . . . . . . . 93A.3.1 Fig 5.3 Half space depth demo . . . . . . . . . . . . . 93A.3.2 Fig 5.6 Projection depth demo . . . . . . . . . . . . . 94A.3.3 Fig 5.8 Simplicial depth demo . . . . . . . . . . . . . 95A.3.4 Contour plots for the eleven depth functions with theexample 6-points dataset . . . . . . . . . . . . . . . . 97A.3.5 Contour plots for the eleven depth functions with 100standard normal bivariate random variable . . . . . . 98A.3.6 Fig 5.16 Simulated examples of different type of datasymmetry . . . . . . . . . . . . . . . . . . . . . . . . 100viTable of ContentsA.3.7 Fig 5.17 Comparison of Mahalanobis depth contoursand density . . . . . . . . . . . . . . . . . . . . . . . . 100A.4 A novel depth-based pairwise quality index . . . . . . . . . . 101A.4.1 Fig 6.1 Type I errors for the two-sample case . . . . . 101A.4.2 Fig 6.5, 6.6, 6.7, 6.2, 6.3, 6.4 Two-sample power com-parisons . . . . . . . . . . . . . . . . . . . . . . . . . 104A.4.3 Fig 6.8, 6.9 Three-sample power comparisons . . . . . 104A.4.4 Power Quantification . . . . . . . . . . . . . . . . . . 106A.4.5 Fig 6.10 Simulation of high-dimension cases . . . . . 107A.4.6 Sec 6.3.3 Data analysis . . . . . . . . . . . . . . . . . 109B An asymptotic normality theorem for the quality index [53] 113B.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . 113B.2 Op and op notation . . . . . . . . . . . . . . . . . . . . . . . 115B.3 Auxiliary lemmas . . . . . . . . . . . . . . . . . . . . . . . . 116B.4 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . 117B.5 Verification of assumptions with Mahalanobis depth (MHD) 118viiList of Tables3.1 Data arrangement for the multivariate case . . . . . . . . . . 53.2 Summary of MANOVA Results . . . . . . . . . . . . . . . . . 74.1 Results of the Standard Wilcoxon rank-sum Test . . . . . . . 114.2 Results of the Wilcoxon rank-sum Test Stratified by Center . 125.1 Data arrangement for t d-variate distributions . . . . . . . . . 525.2 Depth ranks based on pooled sample . . . . . . . . . . . . . . 525.3 Depth ranks based on sample k . . . . . . . . . . . . . . . . . 546.1 Estimated p-values for maximum statistic M30,30,30 and DbRstatistic [46] for each depth function under 10,000 repetitions. 78viiiList of Figures5.1 Presents a contour plot of Mahalanobis depth derived fromthe 100 observations generated from a standard bivariate nor-mal distribution. . . . . . . . . . . . . . . . . . . . . . . . . . 165.2 Left graph: A contour plot using the Mahalanobis depthof the five points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Rightgraph: A contour plot (using the same depth function as theleft graph) of 100 random variates with standard bivariatenormal distribution. . . . . . . . . . . . . . . . . . . . . . . . 195.3 The halfspace depth is the smallest fraction of data points inthe half-space defined by each line. Left: One case when theline passes through the point of interest (red dot) vertically;Middle: One case when the line passes through the point ofinterest horizontally; Right: One case when the line passesthrough the point of interest diagonally. . . . . . . . . . . . . 215.4 Left graph: A contour plot using the half-space depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 225.5 Left graph: A contour plot using the L2 depth of the fivepoints (−0.5, 5),(1, 3.5),(4, 6),(2, 3),(1, 1),(2, 4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255.6 Three-dimensional Representation of Points and their pro-jections. The figure presents four points (A, B, C, D) in athree-dimensional space, with the point of interest (D) high-lighted in red. projections of the points onto the X-axis aredisplayed as blue lines, demonstrating one possible orienta-tion for calculating the projection depth. . . . . . . . . . . . 26ixList of Figures5.7 Left graph: A contour plot using the projection depth of thesix points (−0.5, 5),(1, 3.5),(4, 6),(2, 3),(1, 1),(2, 4) while theleft is of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275.8 Geometrical representation of all possible convex hull of size 3formed by the six points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1),(2, 4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285.9 Left graph: A contour plot using the simplicial depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305.10 Left graph: A contour plot using the zonoid depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.11 Left graph: A contour plot using the spatial depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.12 Demonstration of manual computation of Lens depth of thepoint (1, 3.5) w.r.t to the six points (−0.5, 5), (1, 3.5), (4, 6),(2, 3), (1, 1), (2, 4). . . . . . . . . . . . . . . . . . . . . . . . . 375.13 Demonstration of manual computation of β-skeleton depth ofthe point x w.r.t two other points X1 and X2. . . . . . . . . . 385.14 Left graph: A contour plot using the β-skeleton depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph:A contour plot (using the same depth function as the leftgraph) of 100 random variates with standard bivariate normaldistribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.15 A visual representation of onion depth of 100 random vari-ates with standard bivariate normal distribution. Each colorrepresents a different layer of the “onion”. . . . . . . . . . . . 425.16 Simulated example of different types of data symmetry. Left:A two-drug example to illustrate central symmetry; Middle:A particles diffusion example to illustrate angular symmetry;Right: A weight and hight example to illustrate ellipticalsymmetry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43xList of Figures5.17 Comparison of Mahalanobis depth contours (colours of points)and density (solid line) for a standard bivariate normal dis-tribution for a sample of n = 500. . . . . . . . . . . . . . . . 445.18 Comparison of Mahalanobis depths calculated using momentand MCD methods for a 2-dimensional dataset (n=50). Left:The dataset includes a single outlier at (5,5), marked in red;Middle: The contour lines represent depth regions derivedfrom Mahalanobis depth; Right: The contour lines representdepth regions derived from Mahalanobis depth with MCDestimate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 456.1 Comparison of empirical upper 95% quantiles of three statis-tics Wm,n(0.5, 0.5) (1st row), Wm,n( nm+n ,mm+n) (2nd row),and Mm,n (3rd row) for m = 100, 200, . . . , 1000 and n = m(1st column) or n = m/2 (2nd column). . . . . . . . . . . . . 626.2 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N(0, I2×2 + 0.5I˜2×2) for m = 100, 200, . . . , 1000.Different rows represent distinct depth functions applied, andtwo columns differentiate cases n = m and n = m/2. . . . . 636.3 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.35, 0.35)⊤, I2×2). . . . . . . . . . . . . . . 646.4 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.3, 0.3)⊤, I2×2 + 0.4I˜2×2). . . . . . . . . . . . 656.5 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N(0, I2×2 + 0.5I˜2×2) for m = 100, 200, . . . , 600.Different rows represent distinct depth functions applied, andtwo columns differentiate cases n = m and n = m/2. . . . . 686.6 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.35, 0.35)⊤, I2×2). . . . . . . . . . . . . . . 696.7 Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.3, 0.3)⊤, I2×2 + 0.4I˜2×2). . . . . . . . . . . . 706.8 Comparisons of power of two statisticsMm,n and DbR statis-tic in [46] under alternative hypothesis F (1) = F (2) = N(0, I2×2)and F (3) = N((0, 0)⊤, I2×2 + 0.5I˜2×2). . . . . . . . . . . . . . 726.9 Comparisons of power of two statisticsMm,n and DbR statis-tic in [46] under alternative hypothesis F (1) = N(0, I2×2),F (2) = N((0.3, 0.3)⊤, I2×2), and F (3) = N((0, 0)⊤, I2×2 +0.5I˜2×2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73xiList of Figures6.10 Comparison of power of M statistics using extended Maha-lanobias depth with MBD ranking test with 5000 replicationswhere p = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . 766.11 Scale curves for three species of iris under Mahalanobis depth. 776.12 Scale curves of male Egyptian skulls from the last three dif-ferent time periods under Mahalanobis depth. . . . . . . . . . 786.13 Scale curves for ALL and AML samples under the extendedMahalanobis depth, showing that the ALL group has a largervariability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79xiiAcknowledgementsI wish to express my sincere gratitude to Dr. Xiaoping Shi and Dr. Cindy Fufor the consistent guidance, encouragement, and expertise provided through-out this research journey.xiiiChapter 1IntroductionMultivariate data analysis is a pivotal technique employed across disciplines,such as psychology, biology, medicine, and the broader social sciences [2–4]. When researchers accumulate diverse types of data, advanced statisticalmethods become imperative for accurate interpretation. One of importantmethods is the two-sample homogeneity test, which discerns if two datasets are from the same distribution. For illustrative purposes, within therealm of medicine, these tests might ascertain the efficacy of a medication bycontrasting data prior to and subsequent to its administration. Similarly, insociological studies, these tests might facilitate the comparison of behavioralpatterns across distinct demographic cohorts.Historically, methods such as the t-test and multivariate analysis of vari-ance (MANOVA) have been the preferred tools for such analytical endeavors[10].This thesis is structured to provide a comprehensive examination of thesetests. It will juxtapose conventional methodologies with emerging ones toascertain their respective efficacies under varying circumstances. Addition-ally, this thesis will introduce and elaborate on the concept of ”depth func-tions,” a novel approach poised to enhance multivariate data analysis. Byits conclusion, this study aims to:1. Furnish an expansive overview of both traditional and novel two-sample homogeneity tests.2. Delve into the utility and applicability of depth functions in multivari-ate data analysis.3. Propose a depth-based methodology tailored to optimize two-samplehomogeneity tests in multivariate data scenarios.1Chapter 2Literature reviewThe evolution of MANOVA can be traced back to Ronald A. Fisher’s in-troduction of ANOVA in 1936 [11]. Hotelling then expanded on Fisher’swork, presenting the two-sample multivariate T 2-test [12]. Another notableadvancement came with the introduction of Wilks’ lambda, a criterion de-veloped for testing the equality of population mean vectors using MANOVA[13]. This tool found its application across disciplines. In behavioral andsocial sciences, for instance, it gauges the effects of treatments on humanbehavior [2]. Meanwhile, ecologists use it to discern the impacts of envi-ronmental variables on multiple outcomes [3], and in medical research, itmeasures the efficacy of treatments across multiple health outcomes [4].One of the key features of tools like MANOVA or the t-test is their re-liance on specific assumptions. The standard t-test, for example, presumesnormality, independence, and homogeneity of variances across groups [14].When these assumptions aren’t met, the Welch’s t-test serves as a morerobust alternative, especially when faced with unequal variances and sam-ple sizes [15]. However, when these assumptions are greatly violated, it canreduce the reliability of the results [5]. Furthermore, the F-test is used totest for equality of variances between groups. Its robustness, however, isalso dependent on the assumption of normality [6]. In situations where thesample sizes are small and the data distributions are unknown, the exactparametric tests come into focus. The Fisher’s Exact Test, for instance,is often mentioned for its use in analyzing contingency tables when samplesizes are too small for the chi-squared test to be reliable [7]. Another criticalavenue of parametric testing in the context of two-sample comparisons in-volves the use of Mixed Models, which do not strictly require the assumptionof independence among observations and can manage unbalanced data [8].Lastly, the Hotelling’s T-square test extends the t-test to the multivariatecase, comparing the means of two groups on multiple dependent variablessimultaneously [9]Understanding the limitations of parametric tests, researchers turned tonon-parametric methods for more flexibility. These methods, such as theMann-Whitney U test [20], Wilcoxon signed-rank test [21], and Kruskal-2Chapter 2. Literature reviewWallis test [22, 23], don’t rely heavily on the stringent assumptions requiredby their parametric counterparts [24]. The versatility of non-parametrictests is also seen in tools like the Crame´r test [26] and the Energy Distancetest [27]. The recent work connecting the multivariate Wasserstein test tothe Energy Distance test exemplifies the ongoing evolution and integrationof these methods [28].Despite their flexibility, non-parametric tests come with trade-offs. They’retypically less powerful than parametric tests when the assumptions for thelatter hold true [29].The power of a statistical test is defined as the probability that the testcorrectly rejects the null hypothesis when the alternative hypothesis is true.Mathematically, if β is the probability of making a Type II error (failing toreject the null hypothesis when the alternative hypothesis is true), then thepower of the test is given by:Power = 1− βAdding to the repertoire of comparison tools are statistical depth func-tions. Originating with Tukey’s depth in 1975 [47], these functions gainedtraction because of their versatility. The 1990s saw the introduction of moreadvanced depth functions, some even integrating the Mahalanobis distanceconcept [25, 50]. Their role and application in two-sample tests will bedelved into in subsequent sections [25, 46].3Chapter 3Parametric test of themultivariate homogeneitytestsMultivariate two-sample homogeneity tests are statistical procedures usedto determine whether two multivariate samples are drawn from the samedistribution. Consider two multivariate random samples x1, x2, . . . , xm andy1, y2, . . . , yn drawn from distributions F and G, respectively. The objectiveof the two-sample homogeneity test is to test whether H0 : F = G. Thecorresponding alternative hypothesis suggests that some differences existbetween the two distributions. The typical alternative hypothesis is basedon location shift, scale change, or both.A wide range of research exists on this subject. For instance, the classi-cal parametric MANOVA, an extension of ANOVA, offers a tool for testingmultiple multivariate samples under the assumptions of normality and homo-geneity of covariance matrices. In this chapter, we will describe MANOVAand demonstrate its effectiveness as two-sample homogeneity test statisticsusing numerical simulations and real data analysis.3.1 Introduction to MANOVAMANOVA is a technique that measures the differences between two or moremetric-dependent variables based on a set of categorical variables acting asindependent variables. As an extension of ANOVA, MANOVA provides amore comprehensive approach to analyzing multivariate data, and it remainsthe most commonly used multivariate comparison of group means [10].Suppose there are g samples with n observations in each sample, and pdependent variables, where 1 ≤ i ≤ n and 1 ≤ k ≤ g, and 1 ≤ j ≤ p. Such amultivariate situation can typically be organized in a table, as demonstratedin Table 3.1.43.1. Introduction to MANOVATable 3.1: Data arrangement for the multivariate caseSamples1 2 ... gObservations1 Y11 Y21 ... Yg12 Y12 Y22 ... Yg2............n Y1n Y2n ... YgnLet Yki represent the i-th observation in sample k, Here,Yki =Yki1Yki2...Ykip.In order to use MANOVA, the following assumptions must be met:1. Independence: Data points are randomly and independently collectedfrom the population.2. Normality: Within each group of the independent variables (whichare categorical), dependent variables are multivariate and normallydistributed.3. Interval measurement: Each dependent variable should have an inter-val measurement scale.4. Homogeneity of covariance matrices: The data from all samples havea common covariance matrix Σ.The basic hypothesis of MANOVA is that the sample mean vectors are allequal to one another:H0 : µ1 = µ2 = · · · = µg,where µk is the mean vector of sample k. Rejecting the null hypothesissuggests that at least one sample’s mean vector is significantly differentfrom the others.Wilks’ Lambda, denoted as Λ can be expressed as:53.2. Limitations and challengesΛ = det(E)det(E+H) .where E is the Error Sum of Squares and Cross Products, and H is theHypothesis Sum of Squares and Cross Products. See the detailed descriptionof MANOVA in the textbook [17].It is commonly used for inference about the mean vectors of severalmultivariate normal populations and is a multivariate test statistic used inMANOVA to test the hypothesis. It represents the ratio of the determinantof the error sum of squares and cross-products matrix E to the determinantof the total sum of squares and cross-products matrix T. Smaller valuesof Wilks’ Lambda suggest that the group means are significantly differentfrom one another, whereas values closer to 1 imply no significant differencesamong the group means.3.2 Limitations and challengesDespite its widespread use, MANOVA has some limitations and challenges.It requires certain assumptions to be met, such as multivariate normalityand homogeneity of covariance matrices. Violation of these assumptionsmay lead to inaccurate results [16]. Also, interpreting the results of it canbe complex, mainly when dealing with multiple dependent variables andinteractions between independent variables [18]. Lastly, MANOVA may besensitive to outliers, influencing test results.3.3 Robust MANOVAWith the advent of modern statistical software and computational power,newer techniques have been developed to improve the performance and ro-bustness of MANOVA. One of these advancements is the robust MANOVAwhich has been designed to make it more robust against violations of as-sumptions, such as a robust MANOVA based on the Minimum CovarianceDeterminant (MCD) estimator [19], and another method based on the rankorders of the data [44].The Minimum Covariance Determinant (MCD) estimator introduced by[45] looks for a subset of h observations with lowest determinant of thesample covariance matrix. The size of the subset h defines the so calledtrimming proportion and depending on the desired robustness it is chosenbetween half and the full sample size. The MCD location estimate M is63.3. Robust MANOVAdefined as the mean of that subset and the MCD scatter estimate C is amultiple of its covariance matrix. The multiplication factor is selected sothat C is consistent at the multivariate normal model and unbiased at smallsamples [19].We evaluate the performance and compare the classical, MCD-basedand rank-based MANOVA in terms of their ability to detect group differ-ences and their robustness against outliers and violations of distributionalassumptions using a simulation study. We generate two random samplesx1, x2, ..., xm and y1, y2, ..., yn from two bivariate normal distributions F andG with a mean change. One sample is taken from F = N(0, I2×2) and theother sample comes from G = N((0.2, 0.2), I2×2). Next, an outlier is intro-duced by appending an extreme data point to the dataset with a [10, 10]value.We employed the Wilks.test function from the R package rrcov to con-sider three different methods of mean and variances estimation includingthe standard Wilks’ Lambda estimators, MCD estimators and rank-basedWilks’ lambda as proposed by [44].The results of the analysis are displayed in Table 3.2. The p-value is0.1702 for the classical MANOVA, and less than 2.2e-16 for the MCD-basedMANOVA, indicating a highly significant difference between the groups.For the rank-based MANOVA, the p-value is 4.489e-05, also demonstratinga significant difference between the groups.Table 3.2: Summary of MANOVA ResultsMethod Wilks’ Lambda Chi Square Value DF p-valueClassical 0.96451 3.5417 2.0000 0.1702MCD-based 0 Inf 1.9001 < 2.2× 10−16Rank-based 0.81521 20.022 2.000 4.489× 10−5The degree of freedom of MCD-based method shown in Table 3.2 is notan integer. This is because the distribution of the robust Wilks’ Lambdastatistic based on MCD differs from the classical one, it is necessary to finda good approximation for this distribution. The author of the methodologybehind R package rrcov constructed an approximate distribution based ona Monte Carlo study and examine its accuracy [19]. This implies that thenon-integer degree of freedom might be due to the usage of an approximatedistribution, derived from a Monte Carlo study, to better suit the robustWilks’ Lambda statistic based on MCD.In summary, the robust MANOVA methods, especially the MCD-based73.3. Robust MANOVAapproach, appear to be more effective in detecting group differences andshow greater robustness against outliers and violations of distributional as-sumptions compared to the classical MANOVA method.MANOVA has been a valuable statistical tool for analyzing multivari-ate data across various fields. Despite its limitations, ongoing research andadvancements in statistical methodology have addressed some of these chal-lenges, further enhancing the utility and applicability of MANOVA in con-temporary research.8Chapter 4Non-parametric univariatehomogeneity testsIn this chapter, we study a widely-used nonparametric method, theWilcoxonrank-sum test. We first examine its histories and methodologies followed bya real-world data analysis. Finally, we conduct a review of recent advance-ments in this method.4.1 Introduction to the Wilcoxon Rank-SumTestsThe Wilcoxon rank-sum test, also known as the Mann-Whitney U test, is anon-parametric statistical test that was initially proposed by [20] to deter-mine whether two independent samples originate from the same distribution,specifically focusing on the tendency of one group to have larger values thanthe other. As a non-parametric test, the Wilcoxon rank-sum test has abroad range of applications, particularly when dealing with non-normallydistributed data caused by small sample sizes, ordinal scales, and outliers.• Small Sample Sizes: Researchers across various domains, from neu-roscience to ecology, sometimes have to grapple with small samplesizes [30, 31]. This could potentially affect the power and validity ofparametric tests. As a non-parametric method, the Wilcoxon rank-sum test is robust and reliable in such circumstances, thus making ita preferred choice.• Ordinal Data: Sociology and psychology studies often use ordinalscales, such as Likert scales in surveys or questionnaires [32]. TheWilcoxon rank-sum test is particularly suitable for ordinal data asit compares median ranks rather than actual data values, enablingmeaningful analyses of ordinal data.• Robustness to Outliers: Outliers can significantly influence theresults of parametric tests. In areas such as environmental science94.2. Extensions of Wilcoxon rank-sum testor economics, extreme observations can occur due to various unpre-dictable factors [33]. As the Wilcoxon rank-sum test relies on ranksrather than actual data values, it is more robust to outliers.The Wilcoxon rank-sum test is performed by assigning ranks to datavalues from both groups, ordered from smallest to largest. See textbook [34]for a detailed description of Wilcoxon rank-sum test.In summary, the Wilcoxon rank-sum test provides a robust nonpara-metric procedure for the comparison of two independent groups, provingespecially useful when the assumptions of parametric tests are not met.4.2 Extensions of Wilcoxon rank-sum testSince the proposal of the original Wilcoxon rank-sum test, various extensionshave been developed for more complex scenarios. One notable extension isthe stratified Wilcoxon rank-sum test. Introduced by [35], this test allows forstratification factors in the data, providing a means to control for confound-ing variables. In practice, stratification is widely discussed as an essentialdesign feature in randomized clinical trials (RCTs) that improves the com-parability of groups and increases the statistical power of the study [36–38].RCTs often involve different strata, such as geographic locations (e.g., dif-ferent hospitals or centers in a multi-center trial), patient characteristics(e.g., age, sex, or disease severity), or other factors believed to influence theprimary outcome. Stratification by these variables helps ensure that thetreatment groups are balanced concerning these factors.Another variant of the original Wilcoxon rank-sum test is the clusteredWilcoxon rank-sum test that takes into account intra-cluster dependence[39–41]. Unlike the stratified test which is used when there are potentialconfounding variables that need to be controlled for, the cluster test wouldbe appropriate if observations are grouped and possibly correlated. Similarto the stratified test, these extensions are particularly useful in practicalscenarios in medical studies, where data are often clustered or stratified.Moreover, [41] proposed an extension that is based on within-cluster resam-pling, remaining valid when the cluster sizes are informative. This idea wasthen extended to accommodate situations where the number of members ina group within a cluster is informative [42].Regarding implementation, the base R software provides the wilcox.testfunction for conducting the standard Wilcoxon rank-sum test. For the strat-ified Wilcoxon rank-sum test, R package coin is used as a general tool for104.3. Real data analysisnon-parametric statistical testing, while R package sanon offers specific so-lutions for the analysis of stratified clinical trials with two treatments andpotential covariate adjustment. The clusrank package enables performingthe Wilcoxon rank-sum test on clustered data.4.3 Real data analysisWe now illustrate the usage of the Wilcoxon rank-sum test using the realdataset resp from R package sanon. resp is derived from a randomized clini-cal trial for a respiratory disorder. Patients are randomly assigned to eitherthe group receiving the treatment under investigation (A) or to a groupreceiving a placebo (P). The trial involved two clinical centers where thetrial was conducted and 111 observed patients at four post-baseline visits.Post-baseline visits are follow-up appointments after the initial measurement(baseline) has been taken. A rating of symptom control on an ordinal scale(ranging from 0 = terrible to 4 = excellent) is provided for each patient ateach visit.In our analysis, we aim to compare the rating of symptom control be-tween the test treatment group and the placebo group at each visit. TheWilcoxon rank-sum test was appropriate for this analysis because the ratingof symptom control is an ordinal variable, and this non-parametric test issuitable for such types of data. In Table 4.1, the Z-Value and correspondingp-value for each visit are presented.Table 4.1: Results of the Standard Wilcoxon rank-sum TestVisit Z-Value P-Value Significant at 0.05 levelVisit1 1823.5 0.08147 NoVisit2 2192 6.487e-05 YesVisit3 2041.5 0.002127 YesVisit4 1900.5 0.02709 YesAt Visit 1, no significant difference between treatment groups was found(Z=1823.5, p=0.08147). However, significant differences were observed atVisit 2 (Z=2192, p=6.487e-05), Visit 3 (Z=2041.5, p=0.002127), and Visit4 (Z=1900.5, p=0.02709).These results indicate that the treatment effects differ significantly acrossthe latter three visits, suggesting a notable effect of the test treatment onpatient-reported symptom control over time, emphasizing the value of theWilcoxon rank-sum test in identifying treatment differences.114.3. Real data analysisOur next question is, do the centers have a large differential effect onthe treatment outcomes? To answer this question, we conduct a stratifiedWilcoxon rank-sum test stratified by center. By performing the test withineach center separately for each visit, we aim to account for potential differ-ences in patient populations or practices between centers. Our results weredisplayed in Table 4.2.Table 4.2: Results of the Wilcoxon rank-sum Test Stratified by CenterVisit Z-Value P-Value Significant at 0.05 level1 1.7922 0.0731 No2 4.0304 5.569e-05 Yes3 3.0674 0.00216 Yes4 2.2614 0.02373 YesAt the first visit, the results (Z = 1.7922, p = 0.0731) did not providesignificant evidence of a difference between the two groups. However, atsubsequent visits, the tests yielded significant results indicating a distinctdifference in symptom control ratings between the treatment and placebogroups, with Z and p-values as follows: Visit 2 (Z = 4.0304, p = 5.569e-05),Visit 3 (Z = 3.0674, p = 0.00216), and Visit 4 (Z = 2.2614, p = 0.02373).The results from the Wilcoxon rank-sum tests on the rating of symp-tom control by treatment, with and without stratification by visit, indicatesubtle differences. For Visit 1, stratification resulted in a slightly lower p-value (0.08094) compared to the non-stratified approach (0.08147). Thistrend continues for Visits 2, 3, and 4, where stratification yields marginallysmaller p-values (6.403e-05, 0.002106, and 0.02688 respectively) than theirnon-stratified counterparts (6.487e-05, 0.002127, and 0.02709 respectively).These results suggest that considering the ’visit’ as a stratifying variablemight allow for a more sensitive detection of treatment effects. However,the relatively minor differences in p-values between the stratified and non-stratified approaches indicate that the impact of visit stratification on theassessment of treatment differences is minimal in this specific dataset.Consequently, further considerations might be given to the incorporationof other potential stratifying factors such as the treatment center, age, orsex, which could have a more significant influence on the observed effects.In conclusion, stratification in clinical trials is an important tool to con-trol for confounding variables, reduce bias, and improve the comparabilityof treatment groups. In practice, both analyses can provide valuable in-sights. The non-stratified analysis gives an overall picture, while the strat-124.3. Real data analysisified analysis helps us understand the treatment effect within specific sub-groups, and provides valuable insights when there are underlying differencesin the population that might affect the response variable. This can be es-pecially important in clinical trials where there are often many factors thatcan influence the outcomes. It is worth noting that, while stratification hasclear benefits, it is essential to carefully consider which variables to stratifyby. Over-stratification can lead to small numbers of patients within strata,potentially limiting the effectiveness of randomization and reducing the sta-tistical power of the study [43].13Chapter 5Non-parametric multivariatehomogeneity testsIn the two-sample homogeneity test, nonparametric methods provide robustand flexible tools. The Wilcoxon rank-sum test is a good example when ap-plied to univariate cases. However, when it comes to multivariate data, thedirect generalization of the Wilcoxon rank-sum test is not straightforward.Traditional non-parametric methods, such as the Wilcoxon Rank-Sum test,do not easily extend to the multivariate case due to the complexity of rank-ing.Depth function is one method to solve this problem because it offers away to extend the concept of “ranking” from univariate data to multivariatedata by defining a measure of centrality. Each point in the data space isassigned a “depth” value, indicating how centrally it is located within thedata cloud. The closer the point is to the “center” of the data, the higher itsdepth value, thereby providing an intuitive, geometric way of understandingmultivariate data. This method is particularly beneficial in detecting outliersand understanding the shape of the data.In this chapter, we will first introduce a wide range of depth functionswhich identify the “center” of a multivariate sample. Then, we critically eval-uate their respective strengths and weaknesses. Later, we present two two-sample homogeneity testing methods based on the concept of data depth.Specifically, the Kruskal-Wallis-type statistics H [46] and a quality index Q[25].5.1 Depth functions and one-sample centralitytestThere are some advantages of using statistical depth:1. Depth is free of strong distributional assumptions. Unlike the MANOVAtest, the depth-based test does not require the normality assumption.Moreover, there is a robust version of depth obtained with minimum145.1. Depth functions and one-sample centrality testcovariance determinant (MCD) estimators [51, 52]. The work in [53]listed four weak conditions for the depth function.2. There are many different versions of depth that are easy to access.For example, the R package ddalpha includes Mahalanobis depth [25],spatial depth [54–56], projection depth [56, 57], etc.3. Depth can provide a ranking of multivariate data. The authors of[22, 23] developed the so-called Kruskal–Wallis test based on the rankof univariate samples. Unlike univariate data, multivariate data lack anatural ranking. Recently, the authors of [46] generalized the Kruskal–Wallis test to a multivariate multisample homogeneity test and pro-posed a depth-based rank (DbR) test implemented in the R packageddalpha.Let F (x) be a distribution in d-dimensional space and x ∈ Rd. LetD(x;F ) denote a depth function which is a transformation from Rd into[0, 1] and refers to the depth of the point x with respect to the probabilitydistribution F . The depth function serves as a measure of the centrality or“outlyingness” of a point in relation to a reference distribution or dataset.Here, the F is such a reference distribution. The distribution F plays acrucial role in determining the depth of a point because the depth of a pointx is effectively a measure of its centrality with respect to F . A larger depthvalue indicates that the sample point is closer to the centrality of a givendistribution or data cloud. The work in [56] summarized four postulates ofdepth. Before discussing the four postulates of depth function, let us firstintroduce several key concepts related to the field of data depth.• Level Sets: In Figure 5.1, level sets are regions that contain points ofa similar depth. Since we use a color gradient to represent depth, thelevel sets are contours of the same color. The level sets (= centralregions) Dα(X), α ∈]0, 1], are bounded. That is, each of these regionsencircles points of a specific depth and does not stretch out indefinitely,but rather is contained within a finite area.• Median: Any point x∗ that possesses maximum depth is called a me-dian. It’s the most “central” point in the dataset or distribution, ac-cording to the depth function. For example, in Figure 5.1, the medianwould be a point within the 0.9 contours.Central Symmetry: Given a dataset X exhibiting central symmetrywith respect to a certain point z∗ ∈ Rd, it logically follows that any155.1. Depth functions and one-sample centrality test 0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.9 −2 −1 0 1 2−2−10123Depth values00.20.40.60.81Sample size: 100Figure 5.1: Presents a contour plot of Mahalanobis depth derived from the100 observations generated from a standard bivariate normal distribution.depth function will result in z∗ being identified as the median. Theproperty of central symmetry of X about z∗ implies that the distribu-tion of X − z∗ is identical to the distribution of z∗ −X.Depth functions serve as an effective instrument to detect the cen-trality of points in a distribution. Considering a centrally symmetricdistribution, z∗ is inherently the most central point, as it is the pointabout which the remaining data of the distribution are symmetricallyarranged. As a consequence, it is logical to expect that any depthfunction with an acceptable degree of accuracy would identify z∗ ashaving the highest depth, and thereby, qualify as the median.This characteristic underlines the utility of depth functions. Especiallyin the context of centrally symmetric distributions, depth functionsprovide an effective method to locate the center of symmetry.Now we discuss the four postulates of depth: D is a depth function ifit is invariant against Rd-transformations in some class T , null at infinity,monotone decreasing on rays from its maximum, and upper continuous. Forx ∈ Rd and a distribution F :• T -invariance: D(T (x);T (F )) = D(x;F ) for all T ∈ T .T is specified as the class of affine transformations of Rd. An affine165.1. Depth functions and one-sample centrality testtransformation on Rd is defined as:T (x) = Ax+ bwhere A is a non-singular d× d matrix and b is a vector in Rd.Generally, the statistical depth of a sample remains unaffected by alter-ations in the measurement scale. For instance, converting a sample’sweight unit from kilograms to grams would not impact the statisticaldepth.• Null at infinity: lim||x||→∞D(x;F ) = 0.As ||x|| approaches infinity, the statistical depth of point x should be0, regardless of the center point’s influence on this property.• Monotone on rays: If a point x∗ has maximal depth, that isD(x∗;F ) =maxy∈RdD(x;F ). Then for any r in the unit sphere of Rd, the functionγ 7→ D(x∗+ γr;F ) does not increse with γ > 0. Monotonicity on raysalso means that they are star-shaped about x∗, hence nested. TheMonotonicity ensures that the level sets are star-shaped around themedian x∗.Generally, if a set is star-shaped with respect to a particular point, itmeans that any point within the set can be connected to that pointby a line segment that remains entirely within the set. Applied to thelevel sets in the context of data depth, this means that for each levelset, any point within the set can be connected to the median by a linethat does not cross the boundary of the level set.For example, in Figure 5.1, as we draw any line outward starting fromthe median, we would notice that the points we cross become progres-sively decreasing in depth. We can also see that each level set containsthe next one as we move further from the median, meaning they arenested.• Upper semicontinuous: The upper level setsDα(X) = {z ∈ Rd|D(z;X) ≥ α}are closed for all α ∈ [0, 1].A level set corresponding to a depth value α is the set of all pointsin the data space that have a depth greater than or equal to α. Theproperty of being “Upper semicontinuous” means that these level sets175.1. Depth functions and one-sample centrality testare closed for all α ∈ [0, 1]. A set S is defined as “closed” if it containsall its limit points. A point x is a limit point of S if every openneighborhood around x contains a point in S other than x itself.Numerous concepts of depth have been introduced in existing research,although not all comply with the postulates stated earlier. The followingcontent of this chapter provides a review of eleven widely recognized depthnotions.Mahalanobis depth [50]Mahalanobis depth, based on Mahalanobis distance, offers a measure of cen-trality for multivariate data. Consider a random sample X = {x1, x2, ..., xn}drawn from a distribution F in Rd. The Mahalanobis depth of any point xin Rd with respect to this sample is given by:DMah(x;F ) = (1 + ||x− µX ||2ΣX )−1, (5.1)where µX is the sample mean vector and ΣX is the sample covariance ma-trix. Given z = x − µX , the term ||z||2ΣX = zTΣ−1X z denotes the squaredMahalanobis distance of z in Rd.Example 1. Let us consider six points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4)and compute the Mahalanobis depth for point (1, 3.5).First, we calculate the mean of the dataset:µ = (µx, µy) = ((−0.5+1+4+2+1+2)/6, (5+3.5+6+3+1+4)/6) = (1.5, 3.75).The covariance matrix Σ =( 2.24 0.880.88 2.98). Calculate the Mahalanobis depth forpoint (1, 3.5) as follows:DMah = (1 +√(x− µ)′Σ−1(x− µ))−1 = 0.87We can verify this result using the depth.Mahalanobis function in the Rpackage ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.Mahalanobis(pt, data))The function depth.contours produces a contour plot Fig 5.2 (Left) forthese six points. Fig 5.2 (Right) displays a contour graph of 100 randomvariates with standard bivariate normal distribution for better visualization.185.1. Depth functions and one-sample centrality test 0.2 0.2 0.3 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 2 3 4123456 Depth values00.20.40.60.81Sample size: 6 0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.9 −2 −1 0 1 2−2−10123Depth values00.20.40.60.81Sample size: 100Figure 5.2: Left graph: A contour plot using the Mahalanobis depth ofthe five points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contourplot (using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.Halfspace depth (or location depth or Tukey depth)[47]The halfspace depth [47] of a data point x in Rd with distribution F can bedescribed asDH(x;F ) = inf{P (H) : H is a closed halfspace, x ∈ H}, (5.2)where P is a probability measure, usually based on how many observationsare in the half-space. For a given point x, the half-space depth is defined asthe smallest fraction of the data that a half-space containing x can encom-pass.For univariate data, the halfspace depth of a specific point x is definedas:DH(x;F ) = min{F (x), 1− F (x)}. (5.3)In this scenario, the half-space depth of a point x is simply the minimumof the cumulative distribution function value at x and its complement. Thisis because in a univariate setting, a half-space corresponds to an intervalthat stretches from x to either the leftmost or the rightmost end of the realline. In a two-dimensional space, a half-space is everything on one side ofa line (including the line itself). If we move to a three-dimensional space, ahalf-space is everything on one side of a plane (including the plane itself). Inhigher dimensions, the concept is the same but harder to visualize; it wouldbe everything on one side of a hyperplane.195.1. Depth functions and one-sample centrality testExample 2. To understand how halfspace depth works, we consider a simpleexample in two-dimensional space R2. Consider six points: (−0.5, 5), (1, 3.5),(4, 6), (2, 3), (1, 1), (2, 4). We aim to determine the half-space depth for (1, 3.5).For a given point x in R2, the halfspace depth function calculates the mini-mum number of data points contained in any closed halfspace that includespoint x. A closed halfspace is a region in R2 that is formed by drawing astraight line through x, and then considering one side of the line (includingthe line itself) as the region of interest. The halfspace depth of x is theminimum number of data points found in any of these halfspaces, dividedby the total number of data points.Computing the halfspace depth of x requires analyzing all potential linesthrough x and counts data points within each corresponding closed halfs-pace. Given the infinite possibilities of lines, the R function depth.halfspacein the ddalpha package offers two approaches for depth computation. Theapproximate method calculates the Tukey depth, closely following [73], byprojecting data onto lines uniformly dispersed over the (d− 1)-dimensionalsphere. This determines the depth as the minimum univariate Tukey depthfrom these projections. For more precise calculations, the exact method,as detailed by [49], computes the Tukey depth by considering the minimaldepth from the sum in two orthogonal complementary affine spaces. Thismethod optimizes depth across various data point combinations, based onthe chosen algorithmic dimensionality. Both strategies adeptly handle thechallenge of infinite lines, focusing on randomized projections or orthogonalspaces. A graph illustration of half-sapce depth is given in 5.3.205.1. Depth functions and one-sample centrality test−1 0 1 2 3 4 5 601234567xyVertical Line−1 0 1 2 3 4 5 601234567xyHorizontal Line−1 0 1 2 3 4 5 601234567xyBisector LineFigure 5.3: The halfspace depth is the smallest fraction of data points in thehalf-space defined by each line. Left: One case when the line passes throughthe point of interest (red dot) vertically; Middle: One case when the linepasses through the point of interest horizontally; Right: One case when theline passes through the point of interest diagonally.In this example, we first consider the vertical line passing through (1, 3.5).The points on the left half-space are (0.5, 5), (1, 1) and (1, 3.5), and thepoints on the right half-space are (1, 1), (1, 3.5), (2, 3), (2, 4), (4, 6). Thefraction of points on the left half space is 3/6 and the one on the right halfspace is 5/6. The smaller one is 3/6 = 0.5.Next, let’s consider the horizontal line passing through (1, 3.5). Similarto the previous case, the upper half space contains 4 points while the lowercontains 3 points. The minimum fraction under this scenario is 3/6 = 0.5.Let’s also consider the line that is the bisector of the first quadrant,passing through (1, 3.5). The minimum fraction is the fraction of points inthe upper left half space which is 2/6 = 0.33.The three possible lines are just examples. We have to consider allpossible lines to conduct the half-space depth which is the minimum of theminimum fractions. The result is 0.33 in this case. We can verify this resultusing function depth.halfspace in R package ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.halfspace(pt, data))The function depth.contours produces a contour plot Fig 5.4 (Left) for these215.1. Depth functions and one-sample centrality testsix points. Fig 5.4 (Right) displays a contour graph of 100 random variateswith standard bivariate normal distribution for better visualization. 0.05 0.05 0.1 0.15 0.2 0.25 0 1 2 3 4123456 Depth values00.050.10.150.20.250.30.35Sample size: 6 0.05 0.1 0.15 0.15 0.2 0.25 −2 −1 0 1 2−2−10123Depth values00.10.20.30.40.5Sample size: 100Figure 5.4: Left graph: A contour plot using the half-space depth of the fivepoints (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contour plot(using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.Note that the point with the largest half-space depth would be consideredthe most central point of the data. Intuitively, points towards the “middle”of the data distribution will have a higher half-space depth because no matterhow we draw a line (or a plane, or a hyperplane) through them, we willalways encompass a lot of data points. On the contrary, a lower valuesuggests that x is an outlier.Lp Depth [56]The Lp depth of a point x with respect to a distribution F is defined as:DLp(x;F ) = (1 + EF ||x−X||p)−1, (5.4)where X is a random variable with distribution F and 1 ≤ p ≤ ∞. Thedepth DLp is based on the expected “outlyingness” of a point, as measuredby the Lp distance represented by ||x−X||p. The Lp distance, also knownas Lp norm, Minkowski distance, or generalised distance, is a metric in anormed vector space, Rn, which generalizes the Euclidean distance (p = 2)and Manhattan distance (p = 1). Given two points, x = (x1, x2, ..., xn) andy = (y1, y2, ..., yn) in Rn, the Lp distance between x and y is defined as:dLp(x, y) =(n∑i=1|xi − yi|p) 1p, (5.5)225.1. Depth functions and one-sample centrality testfor 1 ≤ p < ∞. For p = ∞, dL∞(x, y) is defined as the maximum absolutedifference between the components of x and y:dL∞(x, y) = maxi|xi − yi|. (5.6)Example 3. Consider points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4) inR2. We will compute the Lp depth of the point (1, 3.5) for p = 1 and p = 2.We first consider p = 1. The Lp distance between the point of interest(1, 3.5) and each point in the data is calculated as|1− 0.5|+ |3.5− 5| = 2, |1− 1|+ |3.5− 3.5| = 0,|1− 4|+ |3.5− 6| = 5.5, |1− 2|+ |3.5− 3| = 1.5,|1− 1|+ |3.5− 1| = 2.5, |1− 2|+ |3.5− 4| = 1.5.The average of these distances is (2 + 0 + 5.5 + 1.5 + 2.5 + 1.5)/6 = 2.17.The L1 depth is therefore (1 + 2.17)−1 = 0.32. Now, consider p = 2. TheL2 distances are given by:√(1− 0.5)2 + (3.5− 5)2 = 1.80,√(1− 1)2 + (3.5− 3.5)2 = 0,√(1− 4)2 + (3.5− 6)2 = 3.81,√(1− 2)2 + (3.5− 3)2 = 1.12,√(1− 1)2 + (3.5− 1)2 = 2.5,√(1− 2)2 + (3.5− 4)2 = 1.12.The average of these distances is (1.80+0+3.81+1.12+2.5+1.12)/6 = 2.63.The L2 depth is therefore (1 + 1.72)−1 = 0.28. We can calculate L2 depthusing function depth.L2 in R package ddalpha: We can confirm it usingfunction depth.L2 in R package ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.L2(pt, data))It is worth noting that the result of the function depth.L2 is 0.47, which isdifferent from the manually calculated result 0.28. This is because depth.L2calculates depth based on the Mahalanobis distances rather than raw Eu-clidean distances. See more detailed computation in the depth.L2 sourcecode. The following code is how we can use R to help us manually calculateL2 depth with Mahalanobis distance. We can see the result is 4.66, whichis very close to the result of the function depth.L2.235.1. Depth functions and one-sample centrality test# Given pointsdata <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)# Point of interestpt <- c(1, 3.5)# Calculate Mahalanobis distancesdistances <- sqrt(mahalanobis(data, pt, cov(data)))# Calculate average distanceavg_distance <- mean(distances)# Calculate L2 depthL2_depth <- 1 / (1 + avg_distance)print(L2_depth)[1] 0.4661882The function depth.contours produces a contour plot Fig 5.5 (Left) forthese six points. Fig 5.5 (Right) displays a contour graph of 100 randomvariates with standard bivariate normal distribution for better visualization.Both the Mahalanobis depth and the Lp depth with Euclidean distanceare distance-based measures of data depth. They quantify the notion of“centrality” or “outlyingness” of a data point relative to a distribution ordataset using some form of distance metric. While the Mahalanobis depthis calculated based on a single point’s distance from the mean or center ofthe distribution, the Lp depth is calculated based on the expected distanceof a point from all other points in the distribution.Projection depth [56]Let MAD denotes the median absolute deviation. The projection depth ofa point x with respect to a distribution F is defined as:DProj(x;F ) = (1 + supp∈Sd−1| < p, x > −med(< p,X >)|MAD(< p,X >) )−1,245.1. Depth functions and one-sample centrality test 0.28 0.3 0.32 0.32 0.34 0.36 0.36 0.36 0.38 0.4 0.42 0.44 0.46 0 1 2 3 4123456 Depth values0.250.30.350.40.450.5Sample size: 6 0.2 0.2 0.22 0.24 0.26 0.26 0.28 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 −2 −1 0 1 2−2−10123Depth values0.150.20.250.30.350.40.45Sample size: 100Figure 5.5: Left graph: A contour plot using the L2 depth of the five points(−0.5, 5),(1, 3.5),(4, 6),(2, 3),(1, 1),(2, 4). Right graph: A contour plot (usingthe same depth function as the left graph) of 100 random variates withstandard bivariate normal distribution.where X is a random variable with distribution F . Sd−1 is the unit sphere ind dimensions indicating all possible directions in which we could project thepoint x. < p, x > and < p,X > are dot products between the vector p andthe vectors x (or X). < p, x > and < p,X > represent the projection of thepoint x (orX) onto the direction defined by p. med(V ) denotes the median ofa univariate random variable V . MAD(< p,X >) = med(V −med(V )) is themedian absolute deviation of the projections of the data. The quantity | <p, x > −med(< p,X >)| represents the absolute deviation of the projectionof the point x from the median of the projections of the data.By its definition, the projection depth seeks to find the “deepest” locationof the point x relative to the dataset by considering all possible projectionlines p in the d-dimensional space. p ranges over the unit sphere Sd−1 in thed-dimensional space, which essentially means we are considering all possibledirections in which we can project the point x and the data. The supremumover all p in Sd−1 ensures that we are taking the “worst-case” scenario intoaccount. In other words, we are considering the direction in which the pointx is the most “outlying” relative to the data, which will give us the smallestdepth value. By taking the worst case, we ensure that our depth measure issensitive to outliers. This is part of what makes projection depth a robuststatistic. Even if a point is close to the median in most projections butvery far from the median in one particular projection, it will still have a lowprojection depth. The projection depth is not straightforward to compute255.1. Depth functions and one-sample centrality testby hand due to the need to consider all possible projection lines in the space.However, I will visually illustrate a simplified example, projecting onto onlyone axis.Example 4. Let A, B, C, and D be four points in a 3-dimensional spacewhere A = (1, 2, 3), B = (2, 3, 4), C = (3, 2, 1), D = (4, 4, 4). We wish tocalculate the projection depth of point D. Here, we calculate one projectiononto the x-axis. See a visualization in Figure 5.6Figure 5.6: Three-dimensional Representation of Points and their projec-tions. The figure presents four points (A, B, C, D) in a three-dimensionalspace, with the point of interest (D) highlighted in red. projections of thepoints onto the X-axis are displayed as blue lines, demonstrating one possi-ble orientation for calculating the projection depth.Example 5. To better observe the depth contour, we considering anotherexample that determine the projection depth for (1, 3.5) with respect to sixpoints: (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4).A <- c(-0.5, 5)B <- c(1, 3.5)C <- c(4, 6)D <- c(2, 3)E <- c(1, 1)F <- c(2, 4)data <- rbind(A, B, C, D, E, F)265.1. Depth functions and one-sample centrality testP <- c(1, 3.5)print(depth.projection(P, data))The function depth.contours produces a contour plot Fig 5.7 (Left) for thesesix points. Fig 5.7 (Right) displays a contour graph of 100 random variateswith standard bivariate normal distribution for better visualization. 0.15 0.2 0.2 0.2 0.25 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0 1 2 3 4123456 Depth values0.10.20.30.40.50.60.70.8Sample size: 6 0.2 0.2 0.3 0.4 0.5 −2 −1 0 1 2−2−10123Depth values00.20.40.60.81Sample size: 100Figure 5.7: Left graph: A contour plot using the projection depth of thesix points (−0.5, 5),(1, 3.5),(4, 6),(2, 3),(1, 1),(2, 4) while the left is of 100random variates with standard bivariate normal distribution.In practice, the projection depth calculation would involve consideringall possible lines and is usually done with the help of computational softwaredue to the complexity involved.Simplicial depth [58]Simplicial depth assesses centrality using the idea of a convex hull withina multivariate distribution. The term ”conv” denotes the convex hull, themost compact convex set encompassing all given points.Simplicial depth measures a point’s centrality by examining the likeli-hood that it falls inside the convex hull of a randomly selected sample withd+ 1 points, where d represents the number of dimensions. In this context,a simplex requires d+ 1 points. For instance, in a two-dimensional setting,we analyze the probability that a point is enclosed within a triangle definedby three random points. Given X as a random variable with a distributionF , the simplicial depth of a point x is outlined as:275.1. Depth functions and one-sample centrality testDsim(x;F ) = P [x ∈ conv({X1, ..., Xd+1})],Let’s take a look at the univariate equivalent of simplicial depth, as presentedby [76], namely,Dsim(x;F ) = p(x ∈ X1, X2), (5.7)Here x is in R1, X1 and X2 are two independent observation from a univari-ate c.d.f. F and X1, X2 represents the closed line segment connecting X1and X2. When F is continuous,D(x) = 2F (x)[1− F (x)]. (5.8)It follows that any point which maximizes D(x) is a population median. Themaximum value of D(·) is 12 in this case, and D(x) decreases monotonicallyto 0 and x becomes increasingly distant from the median [76].Example 6. In a 2-dimensional case, the convex hull of three points is thetriangle formed by those points. In a 3-dimensional case, the convex hull offour points is the tetrahedron formed by those points. Let’s consider points(−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4). and we aim to determine thesimplicial depth for (2, 1).ABCDEF246−1 0 1 2 3 4xyFigure 5.8: Geometrical representation of all possible convex hull of size 3formed by the six points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4).285.1. Depth functions and one-sample centrality testLooking at these combinations, we can see that the point (1, 3.5) lies withinthe triangles formed by points in 10 combinations:(A,B,C), (A,B,D), (A,B,E), (A,B, F ), (B,C,D),(B,C,E), (B,C, F ), (B,D,E), (B,D,F ), (B,E, F ).For a 2-dimensional setting, the simplicial depth is defined as the number oftriangles formed by any three points from the set that contain the point ofinterest. Figure 5.8 shows all 20 possible triangles, made from combinationsof the six points. Therefore, the simplicial depth of the point (1,3.5) is 10/20= 0.5. However, calculating simplicial depth is not computationally easy,especially when the number of dimensions and data size grow, as it involvescombinatorial calculations of all subsets of size d + 1. We can calculatethe simplicial depth of the point (1,3.5) using function depth.simplicial in Rpackage ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.simplicial(pt, data))The function depth.contours produces a contour plot Fig 5.9 (Left) for thesesix points. Fig 5.9 (Right) displays a contour graph of 100 random variateswith standard bivariate normal distribution for better visualization.Simplicial volume depth (or Oja depth) [59]For any points v1, ..., vd+1 ∈ Rd, the convex hull conv({v1, ..., vd+1}) hasdimensional volumevold(conv{v1, ..., vd+1}) = 1d! |det((1, vT1 )T , ..., (1, vTd+1)T )|. (5.9)The simplicial volume depth or Oja depth is defined byDOja(x;F ) = (1 + E[vold(conv{x,X1, ..., Xd})])−1, (5.10)where X1, ..., Xd are independent copies of X.The simplicial depth quantifies the centrality of a point by counting howmany simplices (formed by the data points) include the point. On the other295.1. Depth functions and one-sample centrality test 0.05 0.05 0.1 0.1 0.15 0.15 0.2 0.25 0.3 0.35 0 1 2 3 4123456 Depth values00.10.20.30.40.5Sample size: 6 0.02 0.04 0.06 0.08 0.1 0.12 0.16 0.2 −2 −1 0 1 2−2−10123Depth values00.050.10.150.20.250.3Sample size: 100Figure 5.9: Left graph: A contour plot using the simplicial depth of the fivepoints (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contour plot(using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.hand, the Oja depth measures centrality by determining the average volumeof the simplices that are created using the data points in combination withthe specified point.Example 7. Let’s consider points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1),(2, 4). and we aim to determine the Oja depth for (1, 3.5).The volume in two dimensions is simply the area. According to theformula, the volume of the convex hull formed by the point (1,3.5) and anypair of other points is the area of the triangle formed by these points. Thearea A of a polygon with n vertices can be calculated using the ShoelaceFormula as follows:A = 12n∑i=1xi(yi+1 − yi−1) (5.11)Here, (xi, yi) are the coordinates of the vertices, n is the number ofvertices, and the vertices are ordered in a counterclockwise direction. y0 andyn are interpreted as yn and y1, respectively. In the case of a triangle withvertices (x1, y1), (x2, y2), and (x3, y3), the area can be calculated directly:Area = 12 |x1(y2 − y3) + x2(y3 − y1) + x3(y1 − y2)| (5.12)305.1. Depth functions and one-sample centrality testConsider Figure 5.8. The area of the triangles formed by each pair of pointstogether with (1,3.5) (Point B) are:Area(ABC) = 4.13,Area(ABD) = 0.38,Area(ABE) = 1.88,Area(ABF ) = 1.13,Area(BCD) = 2.00,Area(BCE) = 3.75,Area(BCF ) = 0.5,Area(BDE) = 1.25,Area(BDF ) = 0.5,Area(BEF ) = 1.25,The expected area is the average of these ten triangles is 1.25. Therefore,the Oja depth of the point (1,3.5) is:DOja((1, 3.5);F ) = (1 + 1.12)−1 ≈ 0.47. (5.13)In more intuitive terms, the Oja depth of a point is high if, on average, ittends to form large-volume simplices with randomly chosen data points. Inother words, if we think of the given point as being in the center of the data,it would be included in a lot of large simplices. However, if the point is nearthe edge of the data, the simplices it forms would generally have smallervolumes, resulting in a lower Oja depth. In this sense, both the Oja depthand the simplicial depth attempt to measure the same underlying concept,but they do so differently.Zonoid depth [60]Before delving into the Zonoid depth, we must first clarify some fundamentalconcepts in measure theory.Sigma-Algebra, F : A sigma-algebra on a set X is defined as a collectionof subsets of X that meets the following criteria:1. The entire set X belongs to F .2. For any subset A in F , its complement Ac is also included in F .3. For a countable series A1, A2, ... in F , their union also lies within F .Borel Sigma-Algebra, B(R): The Borel sigma-algebra on R, denotedas B(R), is the smallest sigma-algebra that contains all open intervals inR. In more formal terms, it’s the sigma-algebra created by the set of openintervals:B(R) = σ ({(a, b) : a, b ∈ R and a < b}) ,where σ(·) refers to the sigma-algebra generated by a given collection of sets.315.1. Depth functions and one-sample centrality testMeasurability: A function g : Rd → [0, 1/α] is termed measurable if, forany Borel set in [0, 1/α], its inverse is recognized as a Borel set in Rd.For 0 < α ≤ 1,DαZon(X) = {E[Xg(X)] : g : Rd → [0, 1/α] measurable and E[g(X)] = 1}(5.14)is the zonoid α-region of X. For α = 0 set D0Zon(X) = Rd. The zonoidα-region of X, denoted by DαZon(X), relies on the concept of an α-region,which in general statistical terms, is a region that contains α fraction of theprobability mass or probability density of the data. While α-region is likea “shape” that contains a certain proportion (α) of the data, the zonoidα-region is a specific way to define this shape using a particular weightingfunction (g). The function g(X) is designed in a way that allocates weightsto different points in the dataset. The points inside the α-region will getmore weight (close to 1/α) and those outside this region will get less weight(close to 0). The zonoid depth is defined asDZon(x;F ) = sup{α : x ∈ DαZon(X)}. (5.15)For any given point x, we find its zonoid depth by identifying the largest αfor which the zonoid α-region contains x. For an empirical distribution X,the zonoid depth is then calculated asDZon(x;F ) = sup{α : αλi ≤ 1/n, x =n∑i=1λixi,n∑i=1λi = 1, λi ≥ 0∀i}, (5.16)where λi is the weight. DZon(x;F ) find the largest α for which x can still berepresented as a weighted average of the points in the empirical distributionunder the given constraints.Example 8. Consider the six points: (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4),we calculate the zonoid depth of the point (1, 3.5).Let λ1, ..., λ6 ≥ 0, and λ1+ ...+λ6 = 1, we find a convex combination of thepoint (1, 3.5) asλ1(−0.5, 5) + λ2(1, 3.5) + λ3(4, 6) + λ4(2, 3) + λ5(1, 1) + λ6(2, 4).Then, we find a combination of coefficients that gives the point of interest(1, 3.5) where all the coefficients are less than or equal to 1/n, where n is325.1. Depth functions and one-sample centrality testthe number of points. Therefore, the aim is to solve the following system ofequations:−0.5λ1 + λ2 + 4λ3 + 2λ4 + λ5 + 2λ6 = 1,5λ1 + 3.5λ2 + 6λ3 + 3λ4 + λ5 + 4λ6 = 3.5,λ1, ..., λ6 ≥ 0, and λ1, ..., λ6 ≤ 1/6.In this context, α is used to control the maximum weight that each λi canhave. When n = 6, each weight λi must be less than or equal to α/6. Oncewe find a solution to this system, the maximum value of α such that thissystem has a solution will be the zonoid depth of the point (1, 3.5). We cancalculate the zonoid depth of the point (1,3.5) using function depth.zonoidin R package ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.zonoid(pt, data))The function depth.contours produces a contour plot Fig 5.10 (Left) for thesesix points. Fig 5.10 (Right) displays a contour graph of 100 random variateswith standard bivariate normal distribution for better visualization. 0.1 0.2 0.3 0.3 0.3 0.3 0.4 0.5 0.6 0.7 0.8 0 1 2 3 4123456 Depth values00.20.40.60.81Sample size: 6 0.1 0.2 0.3 0.4 0.5 0.6 −2 −1 0 1 2−2−10123Depth values00.20.40.60.81Sample size: 100Figure 5.10: Left graph: A contour plot using the zonoid depth of the fivepoints (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contour plot(using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.335.1. Depth functions and one-sample centrality testSpatial depth [61]The spatial depth is defined asDSpa(x;F ) = 1− ||E[ x−X||x−X|| ]||, (5.17)where 0/0=0. Consider a point x and another point X randomly selectedfrom the data set. The vector from x to X is given by X − x. If x is thecenter of the distribution, then on average, the vectors from x to all otherpoints in the distribution will point in all directions equally, cancelling eachother out when summed. This results in the expectation or average vectorE[ X−x||X−x|| ] being close to the zero vector, so its norm ||E[ X−x||X−x|| ]|| is close tozero, and the spatial depth 1− ||E[ X−x||X−x|| ]|| is close to one, indicating thatx is a central point.Example 9. For the six points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4),we compute the spatial depth for (1, 3.5).We calculate the spatial depth for the point (1, 3.5) as follows:E[ x−X||x−X|| ] =16[(1.5√2.25,−1.5√2.25) + (0, 0) + ( −3√10,2.5√10)+ (−1√1,−0.5√1) + ( 0√6.25,2.5√6.25) + ( −1√0.25,0.5√0.25)]= (0.31, 0.06)⇒ DSpa(1,3.5) = 1− ||E[x−X||x−X|| ]|| = 0.68We can verify this result using function depth.spatial in R package ddalpha:data <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.spatial(pt, data))It is worth noting that there might be a slight difference between the spatialdepth value between the manual computation and R function depth.spatial.The function depth.spatial takes into account clusters. It computes spatial345.1. Depth functions and one-sample centrality testdepth for every data point with respect to each specified cluster. Conse-quently, the output is a matrix where each row is aligned with a distinct datapoint and each column corresponds to a cluster. In contrast, our manualcomputation calculates the spatial depth for a single data point in relationto all other data points and does not incorporate the concept of clusters.In addition, depth.spatial interprets spatial depth using the Mahalanobisdistance, which involves the covariance matrix of the data. By contrast,the manual computation calculates spatial depth using Euclidean distance.Hence, depth.spatial is more suitable for multivariate data, while the ap-proach used in the manual calculation is more straightforward and easier tounderstand.The function depth.contours produces a contour plot Fig 5.11 (Left) forthese six points. Fig 5.11 (Right) displays a contour graph of 100 randomvariates with standard bivariate normal distribution for better visualization. 0.1 0.1 0.2 0.2 0.3 0.3 0.4 0.5 0.6 0.7 0 1 2 3 4123456 Depth values00.20.40.60.81Sample size: 6 0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 −2 −1 0 1 2−2−10123Depth values00.20.40.60.81Sample size: 100Figure 5.11: Left graph: A contour plot using the spatial depth of the fivepoints (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contour plot(using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.Lens depth [62]The Lens depth is a measure that uses a dissimilarity function to evaluatehow “central” a point is in relation to other points in a dataset. A dissim-ilarity function is a symmetric function that calculates how different twopoints are from each other. It returns a zero value if and only if the twopoints are the same.355.1. Depth functions and one-sample centrality testIn the context of Lens depth, the dissimilarity function is used to identifywhich one of three given points is the most “central”. Specifically, the point xis considered the most central of the three points x, y, z if it is less dissimilarto each of the other two points than they are to each other. This relationshipcan be represented as:v(x, y) < v(y, z) and v(x, z) < v(y, z), (5.18)where v(., .) is the dissimilarity function. In simpler terms, if point x is closerto both points y and z than they are to each other, then x is considered themost central.Using this concept, the Lens depth of a point is defined as the probabilitythat a given point is the most central among three randomly chosen points.Mathematically, this is expressed as:DLens(x;F ) = P [x ◁ (X1, X2)], (5.19)where X1 and X2 are independent random samples from the data. x ◁(X1, X2) stands for the event v(x, y) < v(y, z) and v(x, z) < v(y, z). If theEuclidean distance is used as the dissimilarity function, the Lens depth of apoint is the probability that the given point is closer to each of two randomlyselected points than they are to each other:DEuclLens(x;F ) = P [max{||x−X1||, ||x−X2||} < ||X1 −X2||]. (5.20)Example 10. Calculate the Lens depth for the point (1, 3.5) with respect to(−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4).To compute the Lens depth, we evaluate all potential pairs of points (ex-clude (1, 3.5)) within our dataset: (−0.5, 5), (4, 6); (−0.5, 5), (2, 3); (−0.5, 5),(1, 1); (−0.5, 5), (2, 4); (4, 6), (2, 3); (4, 6), (1, 1); (4, 6), (2, 4); (2, 3), (1, 1);(2, 3), (2, 4); (1, 1), (2, 4).The Lens depth is then defined as the proportion of these pairs for whichthe point of interest (1,3.5) is more similar to each point in the pair (i.e.,less distant) than the two points in the pair are to each other. We utilize theEuclidean distance as our measure of dissimilarity, denoted as Dist(x, y) fortwo given points x and y. We take the first pair as an example to illustratethe calculation process. As shown in Figure 5.12, point (1,3.5) is labelled inred and the selected points (-0.5,5) and (4,6) are labelled in green. Now wecompute the distance between each pairs of the points:365.1. Depth functions and one-sample centrality testDist((1, 3.5), (−0.5, 5)) = 2.12,Dist((1, 3.5), (4, 6)) = 3.91,Dist((−0.5, 5), (4.6)) = 4.61.In this scenario, the largest distance from (1, 3.5) to either of the twopoints in the pair, maxDist((1, 3.5), (−0.5, 5)),Dist((1, 3.5), (4, 6)) = 3.91, isindeed less than the distance between the two points themselves, Dist((−0.5, 5), (4, 6)) =4.61, labelled in red. This indicates that (1, 3.5) is indeed more central inrelation to the pair (−0.5, 5), (4, 6).Figure 5.12: Demonstration of manual computation of Lens depth of thepoint (1, 3.5) w.r.t to the six points (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1),(2, 4).This procedure of distance comparison can be replicated for every possi-ble pair in our dataset. Consequently, the Lens depth of (1, 3.5) is calculatedas the proportion of such pairs wherein (1, 3.5) emerges as the most centralpoint, according to the principles laid out by the Lens depth calculation.This process, hence, provides a measure of the centrality of a point withina multivariate dataset.375.1. Depth functions and one-sample centrality testβ-skeleton depth [63]The β-skeleton depth is a measure of the centrality of a data point, dictatedby a parameter β ≥ 1. It’s calculated in relation to two other points, X1and X2. The β-skeleton depth is defined asDSkelβ(x;F ) = P [||x−X1β2|| < β/2||X1 −X2||] (5.21)and ||x−X2β1|| < β/2||X1 −X2||],where X1β2 = β2X1 + (1− β2 )X2. It is a point that lies between X1 and X2,calculated as a weighted average of X1 and X2. The β-skeleton depth of xis the probability that x is closer to both X1β2 and X2β1 than X1 and X2are to each other (with their distance scaled by β/2).See Figure 5.13 for a visual demonstration. Consider x, X1, and X2(Labeled in red). X1, and X2 are randomly chosen points from the dataset.Both X2β1 and X1β2 are labeled in green. The point X1β2 is calculated as aweighted average of X1 and X2, i.e., X1β2 = β2X1 + (1− β2 )X2. Here, X1β2lies on the line segment between X1 and X2, closer to X1 if β > 1 and closerto X2 if β < 1. X2β1 is calculated in a similar way as X1β2 but with weightsreversed.Figure 5.13: Demonstration of manual computation of β-skeleton depth ofthe point x w.r.t two other points X1 and X2.385.1. Depth functions and one-sample centrality testBoth the distance between x and X1β2 and the distance between x andX2β1 are labelled in yellow. Neither of them is longer than the distancebetween X1 and X2 scaled by β/2 labelled in blue.Now imagine all other possible choices of X1 and X2 (Not shown in Fig5.13) in a dataset, the β-skeleton depth is the probability that x is closer toboth X1β2 and X2β1 than X1 and X2 are to each other scaled by β/2.Example 11. Calculate the β-skeleton depth for the point (1, 3.5) with re-spect to the dataset: (−0.5, 5), (1, 3.5), (4, 6), (2, 3), (1, 1), (2, 4), whereβ = 2 (chosen arbitrarily for this example).To calculate the β-skeleton depth of a point, we again need to examine allpossible pairs of points (exclude (1, 3.5)) in our dataset: (−0.5, 5), (4, 6);(−0.5, 5), (2, 3); (−0.5, 5), (1, 1); (−0.5, 5), (2, 4); (4, 6), (2, 3); (4, 6), (1, 1);(4, 6), (2, 4); (2, 3), (1, 1); (2, 3), (2, 4); (1, 1), (2, 4).The β-skeleton depth is the proportion of these pairs for which the point(1, 3.5) is more similar (i.e., less distant) to each point in the pair than theyare to each other (scaled by β/2). We can use the Euclidean distance as ourmeasure of dissimilarity denoted as Dist(x, y) for two given points x and y.For example, for the first pair, we first compute the points X1β2 and X2β1,which lie between the points (−0.5, 5) and (1, 3.5), and are calculated as aweighted average of the two points. Let β = 2, compute X1β2 and X2β1 as:X1β2 =β2 (−0.5, 5) + (1−β2 )(−0.5, 5),X2β1 =β2 (1, 3.5) + (1−β2 )(4, 6).Then we compute:Dist((1, 3.5), X1β2) = 2.12,Dist((1, 3.5), X2β1) = 3.91.Finally, we compare these with β/2 × Dist((−0.5, 5), (4, 6)) = 4.61. Sincethe distance from (1, 3.5) to X1β2 and X2β1 is less than β/2 times the dis-tance between (−0.5, 5) and (4, 6), then (1, 3.5) is more central for the pair(−0.5, 5), (1, 3.5).We need to repeat the process for all possible pairs of points in thedataset and calculate the proportion of times the inequality is satisfied tofind the overall β-skeleton depth of the point (1,3.5). We can calculate itusing function depth.halfspace in R package ddalpha:395.1. Depth functions and one-sample centrality testdata <- matrix(c(-0.5, 5, 1, 3.5, 4, 6, 2, 3, 1, 1, 2, 4),ncol = 2, byrow = TRUE)pt <- c(1, 3.5)print(depth.betaSkeleton(pt, data))The function depth.contours produces a contour plot Fig 5.14 (Left) for thesesix points. Fig 5.14 (Right) displays a contour graph of 100 random variateswith standard bivariate normal distribution for better visualization. 0.1 0.1 0.2 0.2 0.2 0.3 0.4 0.4 0.5 0.6 0.7 0 1 2 3 4123456 Depth values00.20.40.60.8Sample size: 6 0.05 0.05 0.05 0.1 0.1 0.15 0.2 0.2 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −2 −1 0 1 2−2−10123Depth values00.10.20.30.40.50.6Sample size: 100Figure 5.14: Left graph: A contour plot using the β-skeleton depth of thefive points (-0.5,5),(1,3.5),(4,6),(2,3),(1,1),(2,4). Right graph: A contourplot (using the same depth function as the left graph) of 100 random variateswith standard bivariate normal distribution.Convex hull peeling depth (or onion depth) [64]The convex hull peeling depth, also known as the onion depth, is a depthmeasure designed specifically for an empirical distribution of a random vari-able X, represented as a finite set SX .In order to understand this depth measure, imagine peeling an onionlayer by layer. At the start, the first “layer” encompasses the entirety of thedata set SX as a convex hull, denoted as C1(X) = conv(SX). A convex hullof a dataset is the smallest convex polygon that encloses all the points inthe dataset.In each subsequent step, we peel away the outermost layer, leaving theinterior of the previous convex hull. Then we form a new convex hull around405.2. Selection of depth functions [65]the remaining points. This new, slightly smaller convex hull is the next layer,denoted as Cj+1(X) for the j-th layer. It can be expressed asCj+1(X) = conv((SX) ∩ intCj(X)), j = 1, 2, ... (5.22)We continue peeling away layers and forming new convex hulls until we reachthe point of interest x. The onion depth of point x is then defined as thenumber of layers we peeled off before reaching x. In other words, it’s thesmallest index j at which the point x is contained within the convex hullCj . The onion depth is then formally represented as:Donion(x;F ) =∑j≥1IintCj(X)(x), (5.23)with IS denoting the indicator function of a set S.Example 12. The figure 5.15 (Right) visually represents the concept ofonion depth in a 2D dataset. Each layer of the “onion” is a convex hullconsisting of points that form the boundary of the current dataset. Layersare “peeled off” sequentially starting from the outermost layer, and a distinctcolor is used to represent each layer. Each peeled layer corresponds to anincrease in onion deptThe function depth.contours produces a contour plot Fig 5.15 of 100random variates with standard bivariate normal distribution.5.2 Selection of depth functions [65]5.2.1 Choosing statistical depth measuresThis section discusses the selection of a statistical depth measure based onthe properties of the data and its generation process. For instance, certaincharacteristics of the data generation process might guide such selection,such as the existence of symmetries or the propensity of the process togenerate outliers. In addition, some data points could be missing or thedistances might only be scaled in an ordinary manner.SymmetryAlthough the actual data might not display perfect symmetry due to noise,measurement errors, and other uncontrolled factors, datasets with symmet-ric characteristics widely exist in various fields and assuming a certain typeof symmetry can often lead to models that are simpler. In these cases, itmakes sense to choose a depth measure that respects this symmetry.415.2. Selection of depth functions [65] 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 −2 −1 0 1 2−2−10123Depth values00.020.040.060.080.10.12Sample size: 100Figure 5.15: A visual representation of onion depth of 100 random variateswith standard bivariate normal distribution. Each color represents a differ-ent layer of the “onion”.• Central symmetry indicates that the data points are symmetric abouta central point. In other words, for every point p in the data, the pointthat is diametrically opposite to p (through the center of symmetry) isalso in the data. Imagine a study measuring the effects of two drugs,where one drug increases a certain metric, and another decreases itby the same amount, leading to a central symmetric pattern in Figure5.16 (Left). In the literature, this type of symmetry is particularlynotable in scenarios involving radial inflow or outflow, such as in thefields of fluid mechanics and heat transfer [66]. In the domain of im-age processing, central symmetry is a common assumption underlyingmany algorithms designed for operations like blurring [67], sharpening[68], and edge detection[69].• Angular symmetry (or rotational symmetry) is a similar concept but inthe circular or spherical sense. Imagine a study of the diffusion of par-ticles in a liquid. The particles are expected to move in all directionsequally, which might result in angular symmetry in the data shown inFigure 5.16 (Middle). This form of symmetry is a common occurrencein astronomical and physical data, where celestial objects and atomicstructures frequently demonstrate rotational symmetry [70]. In thesocial sciences and psychology, angular symmetry could manifest in425.2. Selection of depth functions [65]0 50 100 150 200−3−2−10123Central SymmetryIndexDrug effect−2 −1 0 1 2−3−2−101234Angular Symmetryxy150 160 170 180 1905060708090Elliptical SymmetryHeightWeightFigure 5.16: Simulated example of different types of data symmetry. Left:A two-drug example to illustrate central symmetry; Middle: A particlesdiffusion example to illustrate angular symmetry; Right: A weight and hightexample to illustrate elliptical symmetry.circular models of human behavior or attitudes. A notable exampleis the Circumplex Model of Emotions, which posits that emotionalexperiences are distributed in a circular manner, with diametricallyopposed emotions lying across a central neutral point [71].• Elliptical symmetry means that the data points lie on an ellipse (in2 dimensions) or an ellipsoid (Spherical symmetry in 3 or more di-mensions). An example can be a study on the height and weight ofindividuals. The data might be elliptically distributed shown in Figure5.16 (Right). A multitude of natural and social phenomena modeledusing multivariate normal distributions exhibit elliptical symmetry.For instance, in the financial sector, the returns of multiple assets arefrequently modeled as multivariate normal, thereby showcasing ellip-tical contours [72]. Similarly, in medical statistics, various biometricmeasurements on individuals (like height, weight, blood pressure, etc.)often take on an elliptically symmetric distribution.Affine invariance refers to a property or measure remaining unchangedunder affine transformations. When the data follows an elliptical distribu-tion, it is appropriate to use any affine-invariant depth measure. This isbecause such measures will all yield central regions that are ellipsoidal inshape, which matches the shape of the data distribution, and “recognize”the inherent shape of the data. This is particularly useful in multivariate435.2. Selection of depth functions [65]analysis, where understanding the shape and orientation of the data canhelp. In Figure 5.17, we can see that the points in the center of the datadistribution are coloured with red colour (deepest depth value), This illus-trates the concept that, for an elliptical distribution, the central regions(the points of highest depth) coincide with the density level sets (the areasof highest density).XY−3 −2 −1 0 1 2 3−3−2−10123Figure 5.17: Comparison of Mahalanobis depth contours (colours of points)and density (solid line) for a standard bivariate normal distribution for asample of n = 500.Many multivariate depths are affine invariant, which is often consid-ered as a standard requirement. Examples are the Mahalanobis, halfspace,projection, simplicial, zonoid, and onion depths. However, among all affine-invariant depth measures, the Mahalanobis depth is often preferred in prac-tical applications for a couple of reasons. Firstly, the Mahalanobis depthdirectly accounts for the covariance structure of the data, which makes itparticularly well-suited to the task of measuring depth in multivariate dis-tributions. Secondly, the Mahalanobis depth is computationally efficient tocalculate, especially when compared to other depth measures, which can bemore computationally demanding.445.2. Selection of depth functions [65]Possible OutliersOutliers can heavily influence statistical analyses, leading to skewed results.To reduce their impact, robust depths are often employed to minimize theinfluence of outliers.A robust depth measure can be used for this purpose. Combinatorialdepth measures like halfspace or simplicial depth are recommended.If the data distribution is assumed to be elliptical, with the exceptionof outliers, the Mahalanobis depth can be used. If the outlier exists, a ro-bust covariance estimator, such as the Minimum Covariance Determinant(MCD), is used instead of the traditional covariance matrix. Fig 5.18 (Mid-dle) illustrates the lack of robustness of Mahanobis depth, and Fig 5.18(Right) the improvement with MCD.−1 0 1 2 3 4 5−2024−1 0 1 2 3 4 5−2024Figure 5.18: Comparison of Mahalanobis depths calculated using momentand MCD methods for a 2-dimensional dataset (n=50). Left: The datasetincludes a single outlier at (5,5), marked in red; Middle: The contour linesrepresent depth regions derived from Mahalanobis depth; Right: The con-tour lines represent depth regions derived from Mahalanobis depth withMCD estimate.In Figure 5.18, contours become stretched and distorted due to the out-lier’s influence. The contour lines are considerably less influenced by outliersfollowing optimization using the MCD estimate. The Mahalanobis depthwith the MCD is computationally cheaper than the combinatorial depthmeasures.The spatial depth is also mentioned as a robust measure. However, it isnot affine invariant. To address this, the data can be transformed using arobust scatter transformation, a type of affine transformation. The scattertransformation is applied before the spatial depth calculation, to make the455.2. Selection of depth functions [65]results more reliable.When considering the robustness of a depth measure, the asymptoticbreakdown point of its median is a key factor. The breakdown point is ameasure of the estimator’s robustness against outliers.The breakdown point is a key concept in robust statistics, and it providesa formal way to quantify and compare the robustness of different statisticalmethods. The concept of the breakdown point was introduced in the contextof robust statistics to measure the robustness of an estimator when data maycontain outliers or be contaminated in other ways.Let X = x1, x2, ..., xn be a dataset, and suppose we have an estimator T ,which is a function from our data to [0, 1]. We define the ϵ-contaminationneighborhood of X, denoted as N(X, ϵ), as the set of all datasets Y thatcan be obtained by changing an ϵ proportion of points in X. Here, ϵ is avalue between 0 and 1.Then, the breakdown point of T atX, denoted by ϵ(T,X), is the smallestproportion of data points we can change to make the estimator go to infinity.Example 13. For example, the mean has a breakdown point of 0, becausechanging a single data point can make the mean arbitrarily large. On theother hand, the median has a breakdown point of 0.5, because we need tochange at least half of the data points to make the median arbitrarily large.Thus, the median is considered a more robust measure of central tendencythan the mean.Different depth measures can have different breakdown points. For in-stance, the breakdown point of the Tukey median is at least 1/(d+1), whilethe Oja median and median of simplicial volume regions, onion regions,zonoid regions, and weighted mean regions all have a breakdown point of0, meaning they are not robust against even a single outlier. On the otherhand, the spatial median has a breakdown point of 1/2, making it morerobust.Furthermore, the direction of outliers can also be important. Some depthmeasures, like the zonoid depth, have a projection property, which meansthat they can identify the directions in which a point is an outlier. In suchcases, a point is marked as an outlier if it has a low depth value when mini-mized over the relevant directions. However, these measures are not robust,and their use may be suitable only when the interest is in the direction ofan outlier rather than its depth.465.2. Selection of depth functions [65]Missing dataIn general, when calculating data depth, denoted as D(x;F ), we’re evaluat-ing the measure at a particular point x as well as at n points that representthe support of X. Typically, we would expect these data to be given as a full(n+1)×dmatrix of real numbers. However, it’s important to note that someconcepts of depth can still be applied even when the data is incomplete.If the coordinate values of x are missing, we could limit the depth cal-culation to a subspace of dimensions that correspond to non-missing values.There are certain depth measures that possess what’s known as the projec-tion property:D(x;F ) = infP∈Sd−1D(< p, x >;Fp), x ∈ Rd, (5.24)Here, Sd−1 represents a unit sphere in a d-dimensional space. X followsdistribution F and Fp represents the distribution of the random variable <p,X >, meaning the projection of X on a ray from the origin in the directionof p. When a depth measure has the projection property, it provides anupper limit for the depth value that would otherwise be unknown due tomissing data.On the other hand, when we encounter missing coordinate values in thepoints that make up the support of X, a method to estimate a lower limiton the zonoid depth, DZon(x;F ), can be applied. Assume that values of thej∗th coordinate, xij∗ , are missing at observation units i ∈ J ⊂ 1, 2, ..., n. Wecan then replace each missing value with the average of non-missing values,x¯j∗ = avei/∈Jxij∗ .This imputation (fill in) approach is based on the followingrestriction in the definition of DZon:yj∗ =n∑i=1λixij∗ , that is, yj∗ − x¯j∗ =n∑i=1λi(xij∗ − x¯j∗).The effect of this imputation is to lower at least one of the weights, λi, be-cause their associated data points no longer contribute to the sum. However,this does not decrease the zonoid depth, DZon(x;F ) = infi(nλi)−1. This isbecause infi(nλi)−1 selects the smallest of these inverted values, which cor-responds to the inverse of the largest weight after scaling by n, and makinga weight smaller can’t make the largest weight larger. Hence, despite themissing data and their imputation, we can still get a valid estimate of thezonoid depth.475.2. Selection of depth functions [65]5.2.2 Computational feasibilityIn theory, data depth can be accurately computed for any given point. How-ever, the computational effort required to do so escalates with the numberof points (denoted by n) and even more sharply with the dimensionality ofthe data (denoted by d). This issue is particularly severe for combinatorialdepths.Due to these computational constraints, statisticians sometimes resort toapproximate techniques. When dealing with a large n, a possible approachis to thin out the data by randomly selecting a subset of points. Thisreduces the complexity of the task without sacrificing too much information.On the other hand, handling high-dimensional data (large d) is much morechallenging.This section explores strategies to manage these issues, along with exactand approximate methods for calculating point depth and depth regions.Exact CalculationsIn general, Mahalanobis, L2, and spatial depth calculations are typicallyefficient and rapid, with computation times that are largely independent ofmoderate values for the number of data points (n) and dimensions (d). Ma-halanobis, Lp, and spatial depth exhibit the best possible time complexity.In the context of data depth measures, these depths have the optimal timecomplexity of O(n), where n is the number of data points and the computa-tion time grows linearly with n. This is considered optimal because, at thevery least, each data point must be considered once in these calculations.In a two-dimensional space (d = 2), all types of depth can be efficientlycalculated. However, as the dimensionality increases to (d ≥ 3), computa-tional efficiency drops for the combinatorial invariant depth measures, suchas simplicial and halfspace depth. In these scenarios, the computation time texpands with the size of the dataset as t = nconstd , with a further exponentialgrowth as the number of dimensions increases. Similar trends are observedfor the simplicial volume depth. Among these depth measures, simplicialdepth is the slowest, while halfspace depth is comparatively quicker. Thisperformance trend corresponds to their respective complexities, which areof order nd+1, nd and nd−1log(n).• nd+1: The complexity grows polynomially with the increase in thenumber of data points (n) and their dimensionality (d). It’s repre-sentative of high computational complexity, meaning the computationtime will increase significantly as n or d increases.485.2. Selection of depth functions [65]• nd: The complexity grows proportionally to the number of data pointsto the power of their dimensions. The computation time will increase,but not as rapidly as in the nd+1 scenario.• nd−1log(n): A lower order of complexity compared to the previous two.The log(n) term indicates that as the number of data points increases,the computation time will increase, but at a decreasing rate.The zonoid depth algorithm can handle larger datasets without a sig-nificant increase in computation time and is relatively efficient in handlinghigh-dimensional data as well. Comparatively, the onion depth consistentlyrequires more computation time than the zonoid depth, but exhibits a sim-ilar growth pattern.Approximate CalculationsWhile exact procedures for calculating data depth are the preferred “goldstandard”, practical constraints such as large dimensions (d), large samplesize (n), or frequent evaluations can make these infeasible.In such scenarios, approximate calculations can be implemented. Forinstance, simplicial and simplicial volume depths can be approximated byconsidering a fixed number or a constant portion of all simplices. For theonion depth, the convex hull is often approximately calculated due to highdimensions, leading to potential ambiguities.If a depth satisfies the projection property, it can be approximated usingfinite directions to yield an upper bound of the depth. However, not alldepths can be approximated this way due to nonconvex regions or the needfor drawing simplices. The random Tukey depth [73], which uses randomlychosen directions on the unit sphere, is another approximation method, butthe number of necessary directions may increase exponentially with d [74].The precision of such approximations is also important, with absoluteand relative errors being significant considerations. For example, the randomTukey depth can often result in substantial relative errors, particularly inhigher dimensions.In some cases, exact depths may be more computationally efficient thanrepeatedly approximating values. However, in scenarios such as depth-basedclassification or nonparametric testing based on data depth, where the preci-sion of individual depth values is less important, approximations may suffice.Several other approximate algorithms for Tukey depth exist, offeringpotentially faster solutions for large d [75]. The appropriateness of such495.3. Depth-based two-samples homogeneity testapproximations should be considered in relation to the overall procedureand the acceptable level of error or “noisiness”.Large and High-Dimensioned DataThe application of depth statistics to large or high-dimensional data setsexist in many research filed, especially in the gene expression data [86]. Dif-ferent depth notions may be suitable for different scenarios, and sometimespre-processing might be required.The computational feasibility of a particular depth notion depends notonly on the size (n) and dimension (d) of the data, but also their relativesize. When n significantly exceeds d, many data points will lie on the borderof the data cloud’s convex hull, rendering depth measures like halfspace,simplicial, zonoid, and onion depths ineffective. For scenarios where n < dor the sample covariance matrix is poorly conditioned, adjustments mayneed to be made to depth measures like Mahalanobis depth.When both n and d are large with n > d, Mahalanobis depth can typi-cally be calculated with relative ease. For moderate d, zonoid, spatial, lens,and Lp depth are computationally feasible. However, if these depth proper-ties are not well-suited to the problem, n and/or d may need to be reduced.Reducing n can be achieved by “thinning” or selecting a representative sub-set of the sample, while reducing d involves preselecting relevant features orattributes based on additional information.5.3 Depth-based two-samples homogeneity test5.3.1 Depth-based Kruskal-Wallis-type statistics H [46]A study in 1999 [76] presents the use of the depth-depth (DD) plot, a two-dimensional graphic tool developed for visually comparing two d-dimensionaldistributions utilizing their respective samples. Their work demonstratedthat discrepancies in distributions, such as differences in location, scale,skewness, and kurtosis, could be linked to unique patterns identifiable withinDD-plots.Subsequently, [77] and [78] presented findings which indicated that spe-cific depth functions can serve to effectively characterize distributions.These studies form a foundation for the formulation of formal depth-based test statistics for the comparison of two or more multivariate distri-butions.505.3. Depth-based two-samples homogeneity testFor a univariate sample, the process of data ranking is straightforward.If we have a random sample X1, ..., Xn from a continuous random variableX, we can arrange these data points in ascending order, such as X(1) < ... <X(n). The rank of Xi is defined by the number of data points that are lessthan or equal to Xi, denoted as:R(Xi) = #{Xj : Xj ≤ Xi}Here, #A denotes the cardinality or the number of elements in the set A.This represents the linear ordering and ranking of univariate data.In contrast, for a multivariate sample X1, ...Xn where d ≥ 2, a naturallinear ranking does not exist. However, we can employ the concept of datadepth to formulate a depth ranking. Given a depth function D and a dis-tribution function F , it’s possible to calculate the depths of all the samplepoints X1, ...,Xn. These depths are denoted as D(X1, F ), ..., D(Xn, F ), andthey can be arranged according to increasing depth values, thereby creat-ing a depth ranking of the sample points. The rank of Xi in multivariatesamples is defined as:R(Xi) = #{j;D(Xj , F ) ≤ D(Xi, F )}. (5.25)A larger rank corresponds to a more central position with respect to thedata cloud.In terms of data ranking, it’s critical to note that different distributionfunctions will yield distinct rankings. The empirical depth, D(·, Fn) is ofspecific interest when generating depth-based rankings.Consider a scenario where X1, ...,Xn is a random sample extracted froman unidentified distribution Fn. The ranks of Xi based on the empiricaldepth D(·, Fn) and theoretical depth D(·, F ) will be referred to as Ri andR∗i , respectively.It is reasonable to expect that Ri serves as an empirical approximationof R∗i . However, this approximation will not be accurate without certainadditional regularity conditions.Suppose there are t samples with nt observations in each sample, andd variables, where 1 ≤ i ≤ nt and 1 ≤ j ≤ t, and 1 ≤ k ≤ d. Such amultivariate situation can typically be organized in a table, as demonstratedin Table 5.1.515.3. Depth-based two-samples homogeneity testTable 5.1: Data arrangement for t d-variate distributionsTreatments1 2 · · · tX11 X12 · · · X1tX21 X22 · · · X2t.........Xn11 Xn22 · · · XnttIn the context of t-sample tests, it’s assumed under H0 that there existsa common distribution, denoted by F . For each j-th sample denoted byXj = {X1j , ...,Xnjj}, we calculate an empirical distribution, representedby Fnj . The aggregate empirical distribution, Fn, is then computed as aweighted sum of F jn asFn = n−1t∑j=1njFnj ,where n = ∑nj . Furthermore, the grand sample set, X , is defined as theunion of all the individual sample sets Xj for all treatments asX = ∪tj=1Xj .Finally, let Dij = D(Xij ;Fn) be the depth of the Xij for i = 1, 2, ..., nj ,and j = 1, 2, ..., t with respect to Fn. Write D∗ij = D(Xij ;F ).Let Rij be the depth rank of Xij with respect to the depth functionD(·, Fn). These ranks are shown in Table 5.2.Table 5.2: Depth ranks based on pooled sampleTreatments1 2 · · · tR11 R12 · · · R1tR21 R22 · · · R2t.........Rn11 Rn22 · · · RnttUnder the null hypothesis H0 : F1 = ... = Ft = F , the probability of aparticular observation Xij (indexed by i in the jth sample) obtaining arank r isPrH0(Rij = r) ={n−1 if r = 1, ..., n0 otherwise,525.3. Depth-based two-samples homogeneity testand for two different observations (i, j) ̸= (i′, j′), the joint probability of(i, j) ̸= (i′, j′)) having certain ranks (r, r′) isPrH0(Rij , Ri′j′ = r′) =1n(n−1) if r ̸= r′0 otherwise,Now letR·j =nj∑i=1Rij and R¯·j =R·jnj ,(5.26)so it can be shown thatEH0(R·j) =nj(n+ 1)2 ,VarH0(R·j) =nj(n− nj)(n+ 1)12 ,CovH0(R·j , R·j′) =−njnj′(n+ 1)12 , j ̸= j′.Now considerKj =[R¯·j − EH0(R¯·j)]2VarH0(R¯·j)(5.27)= 12nj(n− nj)(n+ 1)(R¯·j −n+ 12 )2.We take inspiration from the univariate Kruskal-Wallis test and define ourtest statistic as the weighted mean of the Kj values.K =t∑j=1(1− njn)Kj (5.28)= 12n(n+ 1)t∑j=1R2·jnj− 3(n+ 1).Under the null hypothesis, it’s important to mention that the test statisticK is distribution-free. However, a significant limitation of the test basedon K is that it has low power when trying to detect alternative hypothesesinvolving shifts in location [79–81].To develop a test with high power against alternatives involving lo-cation shifts, we should explore ranking procedures that display bettersensitivity to location shifts. As a substitute, let’s examine the depths535.3. Depth-based two-samples homogeneity testDij(k) = D(Xij ;Fnk) and ranksRij(k) for i = 1, 2, ..., nj and j, k = 1, 2, ..., t.Table 5.3 displays the ranks Rij(k), for any k = 1, ..., t.Table 5.3: Depth ranks based on sample kTreatments1 2 · · · tR11(k) R12(k) · · · R1t(k)R21(k) R22(k) · · · R2t(k).........Rn11(k) Rn22(k) · · · Rntt(k)The ranks Rij provide a basic measure of the relative position of theobservation Xij within the entire pooled sample X . This ranking method-ology is commonly employed in statistical tests, such as the Kruskal-Wallistest, to investigate whether several independent samples are from the samedistribution.On the other hand, Rij(k) ranks the observation Xij with respect to thek-th sample only, rather than the pooled sample. This different perspectiveallows for more sensitivity to differences between individual samples, whichcan be crucial in identifying location shifts.Now, it is reasonable to form a test statistic, denoted as H(k), based onthe ranks Rij(k), in a way that is similar to the creation of K as illustratedin formulas 5.27 and 5.28 above. Following that, we might want to pool thetest statistics H(1), ...,H(t) to form a comprehensive test statistic for thenull hypothesis that all samples come from the same distribution. Based onthese considerations, we propose to defineHj(k) =12nj(n− nj)(n+ 1)(R¯·j(k)−n+ 12 )2 (5.29)which is similar to formula 5.28. ThenH(k) =t∑j=1(1− njn)Hj(k) (5.30)= 12n(n+ 1)t∑j=1nj(R¯·j(k)− n+ 12 )2is the proposed text statistics. However, its limitation is that it depends on545.3. Depth-based two-samples homogeneity testthe choice of k. To solve this problem, defineH = 1tt∑k=1H(k) (5.31)= 12n(n+ 1)tt∑k=1t∑j=1nj(R¯·j(k)− n+ 12 )2= 12n(n+ 1)tt∑k=1t∑j=1R2·j(k)nj− 3(n+ 1).A large value of the test statistic H indicates the rejection of the null hy-pothesis.In conclusion, the depth-based Kruskal-Wallis-type statisticsH addressesthe limitation of the test based on K, which suffers from low power whendetecting alternative hypotheses involving location shifts. The concept ofdata ranking was refined and then extended to the consideration of t sam-ples with nt observations for each sample, and d variables. A new statistic,H(k), is formed based on the ranks Rij(k), which provides a measure of therelative position of the observation Xij within the k-th sample. The teststatistic H is then defined as the average of H(k) over all k. This compre-hensive statistic intends to account for the sensitivity to differences betweenindividual samples, which can be crucial in identifying location shifts, andmitigates the dependency on the choice of k.5.3.2 Depth-based quality index Q [25]Based on statistical depth, [25] approached the two-sample problem fromthe angle of quality assurance and proposed the quality index Q. The Qstatistic is designed for the quality control problem, where F represents the“good” population, and G represents a future “check” population. In thecontext of quality control, F is naturally used as a reference distribution,and the Q statistic measures the overall “outlyingness” of G relative to thereference distribution F .Let X ∼ F and Y ∼ G be two independent random variables in Rd.Let D(y,H) be a depth function of a given distribution H in Rd evaluatedat point y. Here is a measure [25] to quantifies the relative outlyingness ofpoint y with respect to distribution F :R(y;F ) = PF (X : D(X;F ) ≤ D(y;F )),555.3. Depth-based two-samples homogeneity test. In the context of two d-dimensional distributions F and G, the Q-statisticis defined as:Q(F,G) =∫R(y;F )dG(y)= P (D(X;F ) ≤ D(Y ;F ) | X ∼ F, Y ∼ G).R(y;F ) represents the proportion of the F population that has a depthlower than point y. Intuitively, R(y;F ) measures the relative outlyingnessof y with respect to F . Therefore, Q(F,G) measures the probability that arandomly drawn point from F will have a depth value that is smaller thanthe depth value of a randomly drawn point from G, when both points areevaluated with respect to the distribution F . The Q(F,G) index can beemployed to detect treatment effects or quality deterioration. When F andG are unknown, the Q-statistic can be estimated by:Q(Fm, Gn) =1nn∑j=1R(yj ;Fm), (5.32)where Fm and Gn represent the empirical distributions of F and G, re-spectively. R(yj ;Fm) is defined as the sample proportion of x1, x2, . . . , xmhaving D(Xi;Fm) ≤ D(yj ;Fm). When considering the empirical distribu-tion, a larger value of Q(Fm, Gn) indicates that Fm has a higher depth valuethan Gn on average, and vice versa.Under the null hypothesis F = G, [25] proved in one dimension d = 1that:((1/m+ 1/n)/12)−1/2(Q(Fm, Gn)− 1/2) d−→ N(0, 1). (5.33)In higher dimensions, [25] conjectured that the same limiting distribu-tion holds for general depth functions and in the general multivariate setting.Regularity conditions [53] were provides for the limiting distribution to holdand verifies the conditions for some commonly used depth functions. Addi-tionally, [53] proves an asymptotic normality theorem and auxiliary lemma,contributing to a better understanding of the Q-statistic and its applicationin various settings.In a more practical situation, F and G are unknown. Let x1, ..., xm bea sample from a reference distribution Fm and y1, ..., yn be a sample fromGn. Then, an estimate of Q isQ(Fm, Gn) =1nn∑i=1R(Fm; yi), (5.34)565.3. Depth-based two-samples homogeneity testwhere R(Fm; yi) is the proportion of xj having D(xj ;Fm) ≤ D(yi;Fm). Itis showed [25] that[ 112(1m+ 1n)]−1 [Q(Fm, Gn)− 12]2 d−→ χ2(1). (5.35)In conclusion, the depth-based quality index Q, as introduced by [25], isa measure to quantify the outlyingness of one population relative to another,specifically designed for quality control applications. In real-world scenar-ios, empirical distributions are used to estimate this statistic. A limitingdistribution [25] is also conjectured for the Q-statistic, which was further ex-plored and confirmed by [53] under certain conditions. Thus, the Q statisticis a useful tool for assessing relative outlyingness between two distributions,widely applicable in fields such as quality control and statistical analysis.57Chapter 6A novel depth-based pairwisequality indexIn this chapter, we investigate the use of the quality index Q beyond thecontext of quality control but to the general two-sample homogeneity prob-lem. In the general two-sample problem, either of the two samples can beregarded as a reference distribution; we need to consider pairwise qualityindexes in order to capture the disparity between the two populations effec-tively. Hence, we propose two new depth-based test statistics constructedfrom pairwise quality indexes for the general two-sample homogeneity test.Interestingly, due to the asymptotical symmetry of the pairwise quality in-dexes, we show that the two proposed test statistics share the same asymp-totic null distribution. Simulations studies have demonstrated the superiorperformance of the proposed tests, including the generalization in multivari-ate multisample cases. The test procedure has been illustrated through tworeal data examples.6.1 Two-sample pairwise quality indexesWm,n(w1, w2) and Mm,nBased on the two pairwise quality indexes Q(Fm, Gn) and Q(Gn, Fm), weconstruct the following two new test statistics, the weighted average statis-ticsWm,n(w1, w2) and the maximum statisticsMm,n to capture the disparitybetween the two populations efficiently:Wm,n(w1, w2) =[ 112(1m+ 1n)]−1w1(Q(Fm, Gn)− 12)2 + w2(Q(Gn, Fm)− 12)2w1 + w2, (6.1)586.1. Two-sample pairwise quality indexes Wm,n(w1, w2) and Mm,nwhere w1, w2 ≥ 0 and w1 + w2 = 1, andMm,n =[ 112(1m+ 1n)]−1max{(Q(Fm, Gn)− 12)2, (Q(Gn, Fm)− 12)2}. (6.2)Note that prior information like the sample size is needed to calculatethe weighted average statistic Wm,n(w1, w2). When w1 = 1, w2 = 0 orw1 = 0, w2 = 1, we treat G or F as a reference distribution. One may setthe weights according to the sample sizes. When m = n, one may use equalweights (i.e. w1 = w2 = 0.5); otherwise set w1 = n/(m + n) = 1 − w2. Inthe general two-sample problem, both proposed test statistics can achievea higher power than the original Q statistic under various alternatives hy-pothesis, as demonstrated by our simulation studies. Recall that the powerof a statistical test is defined as the probability that the test correctly rejectsthe null hypothesis when the alternative hypothesis is true. Mathematically,if β is the probability of making a Type II error (failing to reject the nullhypothesis when the alternative hypothesis is true), then the power of thetest is given by:Power = 1− β.We use two methods to evaluate the discrepancies in the power of differ-ent statistical tests: the Kullback-Leibler (KL) divergence and the averagedifference measure. The KL divergence provides a nuanced assessment ofthe disparity between the expected outcomes of the tests, while the averagedifference gives a direct, mean-based comparison of the test powers.We recommend using the maximum statistic Mm,n when no prior infor-mation like sample size on the weights is available.For a general depth function, [53] studied sufficient conditions to guaran-tee the asymptotic null distributions in (5.33) and (5.34). Under the sameconditions, we establish the asymptotically symmetrical properties of thepairwise quality indexes, which is a key result that leads to the asymptoticnull distributions of the proposed test statistics. Under Conditions A1-A4(See Appendix B), we show that Q(Fm, Gn)−1/2 ≈ 1/2−Q(Gn, Fm). Basedon the asymptotic symmetry property of the pairwise quality indexes, weshow that both proposed test statistics are asymptotic pivotal and have avery simple χ2(1) asymptotic null distribution, as stated in the followingtheorem.Theorem 1. Given two random samples x1, x2, . . . , xm and y1, y2, . . . , yndrawn from distributions F and G, respectively. Let Fm and Gn be thecorresponding empirical distribution function. Consider the two statistics596.2. Simulation studiesMm,n and Wm,n defined in (6.1) and (6.2), respectively. If min(m,n)→∞and m/n tends to a positive constant, under Conditions A1-A4 and whennull hypothesis is true, we have Wm,n(w1, w2) d−→ χ2(1) and Mm,n d−→ χ2(1).Proof. Define I(x, y,H) = I{D(x;H) ≤ D(y;H)}. Under ConditionsA1-A4, [53] showed thatQ(Fm, Gn)−Q(F,G) =(Q(Fm, Gn)−Q(F,Gn)) + (Q(F,Gn)−Q(F,G))=∫∫I(x, y, F )dG(y)d(Fm(x)− F (x))+∫∫I(x, y, F )dF (x)d(Gn(y)−G(y)) + op(m−1/2).Let ξm,n =∫∫I(x, y, F )dG(y)d(Fm(x)−F (x))+∫∫I(x, y, F )dF (x)d(Gn(y)−G(y)), and we write Q(Fm, Gn) − Q(F,G) = ξm,n + op(m−1/2). Note thatF and G have the same distribution under the null hypothesis, then wehave Q(F,G) = Q(G,F ) = 12 and I(y, x,G) = 1 − I(x, y, F ) by definitionof I(x, y,H). Therefore, the pairwise quality indexes possess the followingasymptotic symmetry propertyQ(Gn, Fm)− 1/2 = −ξm,n + op(n−1/2)= 1/2−Q(Fm, Gn) + op(n−1/2) + op(m−1/2). (6.3)By [25, 53], ξ˜m,n := [( 1n +1m)112 ]−1/2(ξm,n)d−→ N(0, 1), which meansξ˜2m,nd−→ χ2(1). Sincem/n tends to a positive constant, [( 1n+ 1m) 112 ]−1[Q(Gn, Fm)−1/2]2 = ξ˜2m,n + op(1) and [( 1n +1m)112 ]−1[Q(Fm, Gn) − 1/2]2 = ξ˜2m,n + op(1),which imply that Wm,n(1/2, 1/2) d−→ χ2(1) and Mm,n d−→ χ2(1). For otherchoices of weights, we can similarly obtainWm,n(w1, w2) d−→ χ2(1) since m/nconverges to a positive constant.6.2 Simulation studies6.2.1 Two-sample cases with the Mahalanobis, spatial andprojection depth [25]In this section, simulation studies are conducted to examine the finite sampleperformance of the proposed tests in the general two-sample problem. Weevaluated the type one error and power of the three proposed statistics,comparing them to the statistics presented in [25].606.2. Simulation studiesThe Kullback–Leibler (KL) divergence can be employed to measure howone power curve differs from another - essentially providing a numericalvalue that reflects the divergence in their respective probabilities. This canbe a powerful tool for comparing the efficacy of different statistical tests orexperimental designs. Let x be a random variable, p(x) and q(x) are twocontinuous probability density functions, then the Kullback-Leibler diver-gence of q from p is defined as:DKL(p ∥ q) = Ep log[p(x)q(x)].For an in-depth description of the KL divergence and its applications, referto the textbook [83]. Alternatively, we can calculate the average of theabsolute differences at corresponding points to compare two power curves.If Pi and Qi represent the power at point i for each curve, then the averagedifference is given by:DAvg =1nn∑i=1|Pi −Qi|.This method provides a straightforward measure of the overall disparitybetween the curves.To maintain alignment with the experiment in [25], we utilized the samethree depth functions: Mahalanobis, spatial, and projection depth. Wegenerate two random samples x1, x2, . . . , xm and y1, y2, . . . , yn from distri-butions F and G, respectively.We first consider the Type I error of three statistics Wm,n(12 ,12),Wm,n( nm+n ,mm+n), and Mm,n. Let F = G = N(0, I2×2), where N(0, I2×2)represents the bivariate normal distribution with mean vector 0 and two-by-two identity covariance matrix. We set m = 100, 200, . . . , 1000 and n = mor m/2. By Theorem 1, all these statistics have the same χ2(1) asymp-totic null distribution. We consider the upper 95% quantile of χ2(1), 3.84,and compare it with the empirical quantiles of three statistics for differentconfigurations of m and n with 100,000 repetitions. Figure 6.1 shows theconvergence rate of empirical quantiles to the theoretical one. When m = n,Wm,n(0.5, 0.5) and Wm,n( nm+n ,mm+n) are the same and have a fast rate ofconvergence; otherwise both rates are slower as shown in the first two rows.The convergence rate forMm,n is the slowest. Other simulations (not shownhere for the sake of brevity) show that the quantile of Mm,n approaches thenominal value when m exceeds 5000. Overall, all three depth functions lead616.2. Simulation studies200 400 600 800 10004.04.55.05.5m (n=m)Critical value for equal weightsMahalanobisSpatialProjection200 400 600 800 10004.04.55.05.5m (n=m)Critical value for unequal weightsMahalanobisSpatialProjection200 400 600 800 10004.04.55.05.5m (n=m)Critical value for maximum MahalanobisSpatialProjection200 400 600 800 10004.04.55.05.5m (n=m/2)Critical value for equal weightsMahalanobisSpatialProjection200 400 600 800 10004.04.55.05.5m (n=m/2)Critical value for unequal weightsMahalanobisSpatialProjection200 400 600 800 10004.04.55.05.5m (n=m/2)Critical value for maximum MahalanobisSpatialProjectionFigure 6.1: Comparison of empirical upper 95% quantiles of three statisticsWm,n(0.5, 0.5) (1st row), Wm,n( nm+n ,mm+n) (2nd row), and Mm,n (3rd row)for m = 100, 200, . . . , 1000 and n = m (1st column) or n = m/2 (2ndcolumn).to comparable results, except that Mm,n has a faster rate of convergencewhen Mahalanobis depth is used for computation.Next, we consider the power of six test statistics: Wm,n(1, 0),Wm,n(0, 1),Wm,n(0.5, 0.5), Wm,n( nm+n ,mm+n), Mm,n, and a depth-based rank (DbR)statistic [46]. The estimated power is calculated based on the simulatedupper quantile α = 0.05 for each statistic with Mahalanobis depth, spatialdepth, and projection depth as the underlying depth function. The powerstudy is carried out for the following three alternatives:(1) Two bivariate normal distributions with a scale change. One sampleis taken from F = N(0, I2×2) , and the other sample comes from G =N(0, I2×2 + 0.5I˜2×2), where I˜2×2 = ((0, 1)⊤, (1, 0)⊤). Figure 6.2 shows the626.2. Simulation studies200 400 600 800 10000.00.40.8m (n=m)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m)Power (Projection)200 400 600 800 10000.00.40.8m (n=m/2)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m/2)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m/2)Power (Projection)Wm,n(1, 0) Wm,n(0, 1) Wm,n(12, 12) Wm,n(nm + n, mm + n) Mm,n DbRFigure 6.2: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N(0, I2×2+0.5I˜2×2) for m = 100, 200, . . . , 1000. Different rowsrepresent distinct depth functions applied, and two columns differentiatecases n = m and n = m/2.power comparison for different sample sizes and various depth functions with1,000 repetitions. Compared to the difference between balanced sample sizeand unbalanced sample size, the difference among three depth functionsis very small. The maximum statistic Mm,n and one particular weightedstatistic Wm,n(0, 1) have the largest power but Wm,n(1, 0) has the smallestpower. In this case, F is preferred to be the reference distribution. Inpractice, since it is not clear which reference distribution should be used tocapture the disparity between the two samples more effectively, the maximalstatistic is recommended. It is interesting to see that the DbR statistic iscomparable to the weighted statistic with equal weights.(2) Two bivariate normal distributions with a mean change. One sam-ple is taken from F = N(0, I2×2) and the other sample comes from G =N((0.35, 0.35)⊤, I2×2). Different from the first scenario, Figure 6.3 shows636.2. Simulation studies200 400 600 800 10000.00.40.8m (n=m)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m)Power (Projection)200 400 600 800 10000.00.40.8m (n=m/2)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m/2)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m/2)Power (Projection)Wm,n(1, 0) Wm,n(0, 1) Wm,n(12, 12) Wm,n(nm + n, mm + n) Mm,n DbRFigure 6.3: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.35, 0.35)⊤, I2×2).that the maximum statistic Mm,n outperforms all other five statistics todetect a location shift between the two samples.(3) Two bivariate normal distributions with both mean and scale changes.One sample is taken from F = N(0, I2×2), and the other sample is drawnfrom G = N((0.3, 0.3)⊤, I2×2+ 0.4I˜2×2). Figure 6.4 shows that Mm,n andWm,n(0, 1) are comparable, and both outperform all other four test statis-tics.We also consider other alternatives. For example, F = N(0, I2×2) and Gis a bivariate t distributions with mean (0, 0)⊤, the scale matrix being theidentity matrix, and the degree of freedom being 5. All six statistics havecomparable performance. The results are not presented here for the sake ofbrevity.In addition to comparing the power discrepancies through visual powercurves, we have the option to quantify these differences with two analyticalmethods: the Kullback–Leibler (KL) divergence and an average difference646.2. Simulation studies200 400 600 800 10000.00.40.8m (n=m)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m)Power (Projection)200 400 600 800 10000.00.40.8m (n=m/2)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=m/2)Power (Spatial)200 400 600 800 10000.00.40.8m (n=m/2)Power (Projection)Wm,n(1, 0) Wm,n(0, 1) Wm,n(12, 12) Wm,n(nm + n, mm + n) Mm,n DbRFigure 6.4: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.3, 0.3)⊤, I2×2 + 0.4I˜2×2).656.2. Simulation studiesapproach. These methods allow us to systematically compare the statisticalpower of the tests under various scenarios.For instance, under the alternative hypothesis F = N(0, I2×2) versusG = N(0, I2×2 + 0.5I˜2×2) with sample sizes ranging from m = 100 to 600with n = m/2, the KL divergence and average difference are calculated forboth M and W (1, 0). Moreover, we consider a second scenario under thealternative hypothesis F = N(0, I2×2) versus G = N((0.35, 0.35)⊤, I2×2)with m = n.Here, the KL divergence and average difference help quantify the dispar-ity in test power. Using the KL.plugin function from the R package entropy,we find a KL divergence of 0.003 and an average difference of 0.005 for thefirst scenario, suggesting a smaller power difference. In contrast, the secondscenario yields a higher KL value of 0.049 and an average difference of 0.022,indicating a more substantial power difference between the tests.6.2.2 Two sample power comparison with other depthsIn the previous section, we evaluated the type one error and power of thethree proposed statistics, comparing them to the statistics presented in [25].To maintain alignment with the experiment in [25], we utilized the samethree depth functions: Mahalanobis, spatial, and projection depth. In thissection, we conduct an analogous experiment to the previous one, but withthe application of different depth functions mentioned in the previous chap-ter to further assess the performance of our proposed statistics. The depthswe explore include: β-Skeleton, L2, onion, simplicial and zonoid depths.(1) Two bivariate normal distributions with a scale change. One sam-ple is taken from F = N(0, I2×2) , and the other sample comes fromG = N(0, I2×2 + 0.5I˜2×2), where I˜2×2 = ((0, 1)⊤, (1, 0)⊤). Figure 6.5 showsthe power comparison for different sample sizes and various depth func-tions with 500 repetitions. In this case, all statistics lose power with theβ-Skeleton depth. When applying the L2 and Onion(qhpeeling) depths,the maximum statistic Mm,n and a specific weighted statistic Wm,n(0, 1)demonstrate the highest power, whereas Wm,n(1, 0) has the lowest power.Conversely, when the simplicial and zonoid depths are applied, DbR demon-strates the highest power, with Mm,n and Wm,n(0, 1) showing comparableperformance. However, both instances of Mm,n outperform at least one ofthe W statistics.(2) Two bivariate normal distributions with a mean change. One sam-ple is taken from F = N(0, I2×2) and the other sample comes from G =N((0.35, 0.35)⊤, I2×2). Different from the first scenario, Figure Figure 6.6666.2. Simulation studiesdemonstrates that with the application of β-Skeleton, L2, and zonoid depths,Mm,n over-performs the other four test statistics. Conversely, when the sim-plicial and Onion(qhpeeling) depths are utilized, Wm,n( nm+n ,mm+n) has thebest performance.(3) Two bivariate normal distributions with both mean and scale changes.One sample is taken from F = N(0, I2×2), and the other sample is drawnfrom G = N((0.3, 0.3)⊤, I2×2+ 0.4I˜2×2). 6.7 demonstrates that with theapplication of β-Skeleton and L2 depths, Mm,n over-performs the other fourtest statistics. Conversely, when the simplicial and zonoid depths are uti-lized, DbR has the best performance. When the Onion(qhpeeling) depth isapplied, Mm,n and Wm,n(0, 1) showing comparable performance.676.2. Simulation studiesFigure 6.5: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N(0, I2×2 + 0.5I˜2×2) for m = 100, 200, . . . , 600. Different rowsrepresent distinct depth functions applied, and two columns differentiatecases n = m and n = m/2.686.2. Simulation studiesFigure 6.6: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.35, 0.35)⊤, I2×2).696.2. Simulation studiesFigure 6.7: Power comparison under alternative hypothesis F = N(0, I2×2)against G = N((0.3, 0.3)⊤, I2×2 + 0.4I˜2×2).706.2. Simulation studies6.2.3 Multi-sample comparisonAs demonstrated in the last section, the maximal statistic Mm,n is morepowerful than the weighted statistic in a two-sample comparison, especiallyagainst location shift alternatives. In this section, we extend the maximumstatistic Mm,n for comparing multisample multivariate distributions. Con-sider k random samples that are drawn from distributions F (j) with samplesizes nj and empirical distribution F (j)nj , for j = 1, 2, . . . , k.The generalized maximum statistic for k-sample comparison isMn1,...,nk = max1≤i,j≤k,i̸=j[112(1ni+ 1nj)]−1 [Q(F (i)ni , F(j)nj )−12]2. (6.4)Similar to the proof of Theorem 1, for x > 0 we haveP (Mn1,...,nk ≤ x)→ P{max1≤i<j≤k(ci,jZi + c˜i,jZj)2 ≤ x}, (6.5)where Z1, Z2, . . . , Zk are independent from N(0, 1), ci,j = limn−1/2i (n−1i +n−1j )−1/2, and c˜i,j = limn−1/2j (n−1i + n−1j )−1/2 with c2i,j + c˜2i,j = 1.When k = 3, we haveMn1,n2,n3 = max{[112(1n1+ 1n2)]−1[Q(F (1)n1 , F(2)n2 )−12]2,[( 112(1n1+ 1n3)]−1[Q(F (1)n1 , F(3)n3 )−12]2,[( 112(1n2+ 1n1)]−1[Q(F (2)n2 , F(1)n1 )−12]2,[( 112(1n2+ 1n3)]−1[Q(F (2)n2 , F(3)n3 )−12]2,[( 112(1n3+ 1n2)]−1[Q(F (3)n3 , F(2)n2 )−12]2,[( 112(1n3+ 1n1)]−1[Q(F (3)n3 , F(1)n1 )−12]2}.andP (Mn1,n2,n3 ≤ x)→ P{(c1,2Z1 + c˜1,2Z2)2≤ x, (c1,3Z1 + c˜1,3Z3)2 ≤ x, (c2,3Z2 + c˜2,3Z3)2 ≤ x. (6.6)The asymptotic properties of the generalized maximum statistic for k-sample comparison are left for future research. Here, we conduct simulation716.2. Simulation studies200 400 600 800 10000.00.40.8m (n=k=m)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=k=m)Power (Spatial)200 400 600 800 10000.00.40.8m (n=k=m)Power (Projection)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Spatial)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Projection)Mm,n DbRFigure 6.8: Comparisons of power of two statistics Mm,n and DbR statisticin [46] under alternative hypothesis F (1) = F (2) = N(0, I2×2) and F (3) =N((0, 0)⊤, I2×2 + 0.5I˜2×2).studies to evaluate the power performance of the maximum statisticMm,n,k,and the DbR statistic. The estimated power is based on the simulatedupper quantile α = 0.05 for each statistic with Mahalanobis depth, spatialdepth, or projection depth as the underlying depth function. We considerthe following two alternatives:(1) Three bivariate normal distributions: Two samples are taken fromF (1) = F (2) = N(0, I2×2), and the third sample comes fromF (3) = N((0, 0)⊤, I2×2 + 0.5I˜2×2). Figure 6.8 compare the power of twostatistics Mm,n,k and DbR statistic for m = 100, 200, . . . , 1000, three depthfunctions and two combinations n = k = m and n = 2k = m/2 for 1,000repetitions. We can see that the maximal statistics Mm,n,k outperforms theDbR statistic.(2) Three distinguished bivariate normal distributions: F (1) = N(0, I2×2),F (2) = N((0.3, 0.3)⊤, I2×2), and F (3) = N((0, 0)⊤, I2×2 + 0.5I˜2×2). Figure726.2. Simulation studies200 400 600 800 10000.00.40.8m (n=k=m)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=k=m)Power (Spatial)200 400 600 800 10000.00.40.8m (n=k=m)Power (Projection)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Mahalanobis)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Spatial)200 400 600 800 10000.00.40.8m (n=2k=m/2)Power (Projection)Mm,n DbRFigure 6.9: Comparisons of power of two statistics Mm,n and DbRstatistic in [46] under alternative hypothesis F (1) = N(0, I2×2), F (2) =N((0.3, 0.3)⊤, I2×2), and F (3) = N((0, 0)⊤, I2×2 + 0.5I˜2×2).6.9 shows that the maximal statistic Mm,n,k outperforms the DbR statisticfor all the cases considered.In summary, our simulation studies show that the maximal statisticMm,n,k is more powerful than the DbR statistic under various alternativesin three-sample comparisons. The choice of depth function does not signifi-cantly change the performance of the tests.6.2.4 High-dimension casesOur findings can be further extended to high-dimensional data, a commonscenario in numerous statistical applications, most notably in biologicalstudies [84–86]. Motivated by genetic sequence investigations [86], we fo-cus on a two-sample homogeneity test under conditions where the numberof dimensions (p) is large and the sample size (n) is small. This situation, of-ten referred to as the ”large p, small n” problem, presents unique challenges736.2. Simulation studiesand opportunities in statistical analysis.It is well-known that the sample covariance matrix is singular when p >n, but a valid covariance matrix must be positive-definite. This fatal flawhampers the application of sample covariance matrix in high-dimensionalmultivariate statistical analyses [88].In order to manage the issue of the sample covariance matrix, we refinedthe M-statistics using an extended version of Mahalanobis depth. Recallthat the original Mahalanobis depth is defined asDMah(x;F ) = (1 + ||x− µX ||2ΣX )−1 (6.7)where µX and ΣX denote the expectation vector and covariance matrixof X respectively. The Mahalanobis norm of z ∈ Rd is given by ||z −µX ||2ΣX = zTΣ−1X z. It measures the distance of the vector z from the meanµX , taking into account the covariance of X. The key feature of the originalMahalanobis depth is that it utilizes the empirical covariance matrix ΣX .The extended Mahalanobis depth, however, improves the estimation ofthe covariance matrix by using the oracle approximating shrinkage estima-tor, as proposed in [87]. This approach estimates the covariance matrixas:Σˆ = ρFˆ + (1− ρ)ΣX (6.8)Here, ρ ∈ (0, 1) serves as a control parameter or weight. ΣX representsthe empirical covariance matrix, while Fˆ stands for a target matrix. It’sproposed that a structured estimate with Fˆ = Tr(ΣX/p) · Ip×p is used,where Ip×p is a p-dimensional identity matrix.According to the source code, there are two steps to compute ρ. Theinitial stage is about computing two intermediary terms, denoted as term1and term2.term1 =(1− 2p)Tr(ΣX2) + Tr(Σ2X), (6.9)term2 = (n+ 1− 2p)(Tr(ΣX2)−Tr(Σ2X)p), (6.10)where p stands for the quantity of features or columns present in the dataset.With the intermediary terms, ρ is derived through a relationship betweenterm1 and term2, and its bounds fall within the 0 to 1 range. It is computedusing the following formula:746.3. Real data analysisρ = max(min(term1term2 , 1), 0). (6.11)The extended Mahalanobis depth thus refines the covariance matrix usedin the depth computation, leading to more robust depth values.We compare the power of the M statistics using the extended Maha-lanobis depth with the modified band depth (MBD) statistics [82] which isalso designed for high-dimension data. The finite-dimensional MBD for ad-dimensional point y in the sample y1, ...,yn is as follows:MBDd(y) =(n2)−1 ∑1≤i1<i2≤nd−1 (6.12)×d∑k=1I{min{yi1(k),yi2(k)}≤y(k)≤max{yi1(k),yi2(k)}}.In our simulation study, one sample is taken from N(1, 0.1|i−j|) where1 ≤ i ≤ n and 1 ≤ j ≤ p. Another sample comes from N(1, (k×0.1)|i−j|). kranges from 2 to 9. The samples N(1, 0.1|i−j|) and N(1, (k × 0.1)|i−j|) havedifferent covariance structures. The parameter k multiplies the covarianceof the second sample, increasing the covariance as k ranges from 2 to 9. TheM statistics perform well under this case where the covariance of the twosamples varies 6.10.6.3 Real data analysisIn this section, we illustrate our proposed maximum statistic for three-sample comparison in the analysis of two real-data examples.6.3.1 Iris dataThe iris dataset is a built-in dataset in R that contains measurements onsepal length and width as well as petal length and width for three speciesof iris setosa, sersicolor, and virginica. There are 50 measurements for eachspecies. The measurement unit is in centimeters. For the iris dataset, weare interested in whether the joint distributions of sepal length, sepal width,petal length, and petal width are different across three species of iris.To visualize the dispersion of distributions of three species, we comparescale curves proposed by [76] and implemented in R package DepthProc756.3. Real data analysis2 3 4 5 6 7 8 90.00.20.40.60.81.0MMBDFigure 6.10: Comparison of power of M statistics using extended Maha-lanobias depth with MBD ranking test with 5000 replications where p = 100.which allows us to graphically compare the scale of multiple distributions.Let Dα(F ) be the α-trimmed region with respect to distribution F , that isDα(F ) ={x ∈ Rd : D(x;F ) ≥ α}.Similarly, we have sample versions of D(x;Fm) and Dα(Fm). Let V (α;Fm)be the volume of convex region Dα(Fm). Define a sample scale curve bytaking the plot of 1 − α versus the volume V (α;Fm). The faster-growingscale curve is associated with a larger scale of distribution. Figure 6.11presents the scale curves for the three species of iris under Mahalanobisdepth.We can see that the scale curves of the three species are non-overlapping,which indicates the scales of those three samples are different. This mayimply that those three species may have different distributions. In fact, theestimated p-value is zero for each depth function and n1 = n2 = n3 = 50under 10,000 repetitions. The asymptotic p-value by (6.6) is also zero foreach depth function.6.3.2 Skull dataThe researchers have suggested that changes in skull size over time are ev-idence of interbreeding between a resident population and a migrant popu-lation. The R package HSAUR includes the male Egyptian skulls data with766.3. Real data analysis0.00.20.40.60.00 0.25 0.50 0.75 1.001 − alphaVolume namessetosaversicolorvirginicaFigure 6.11: Scale curves for three species of iris under Mahalanobis depth.four measurements (maximal Breadth of skull, basibregmatic Height of skull,basialveolar Length of skull, and nasal Height of skull) from five differenttime periods, 4000 B.C., 3300 B.C., 1850 B.C., 200 B.C., and 150 A.D..There are 30 measurements for each time period. We wish to determine ifthere are any differences in the skull sizes in the last three time periods.Figure 6.13 presents the scale curves of skulls from the last three differenttime periods under Mahalanobis depth. Since the values of V (α;Fm) de-pend on the units of observations, it is reasonable to have distinct ranges ofV (α;Fm) in iris and skull datasets.We can see that the scale curves of skulls from 1850 B.C. and 4000 B.C.are very close. Furthermore, we have estimated the p-values for maximumstatistic M30,30,30 and DbR statistic for each depth function under 10,000repetitions. Table 6.1 shows all estimated p-values which are greater than0.05. This confirms that there is no strong evidence to show the interbreed-ing between a resident population and a migrant population in those threetime periods. We note that the asymptotic p-value by (6.6) is 0.027, 0.015,and 0.052 for Mahalanobis, spatial, and projection depth, respectively, whichare much smaller than the corresponding estimated p-values.776.3. Real data analysis0500010000150000.00 0.25 0.50 0.75 1.001 − alphaVolume names150 A.D.1850 B.C.4000 B.C.Figure 6.12: Scale curves of male Egyptian skulls from the last three differenttime periods under Mahalanobis depth.Table 6.1: Estimated p-values for maximum statistic M30,30,30 and DbRstatistic [46] for each depth function under 10,000 repetitions.Maximum statisticM30,30,30DbR statisticMahalanobis 0.3442 0.1036spatial 0.2629 0.0677projection 0.3977 0.14576.3.3 leukemia dataThe leukemia gene expression dataset from the ReorderCluster R pack-age comprises expression levels of 100 genes across two leukemia classes:Acute lymphoblastic leukemia (ALL), represented by 47 samples, and Acutemyeloid leukemia (AML), represented by 25 samples. This dataset is a crit-ical tool for testing the homogeneity of the two leukemia populations understudy. A key element in this analysis is the scale curves, which provide astraightforward and effective mechanism for visualizing and comparing thedispersion of these high-dimensional samples. By leveraging this resource,We can effectively understand and identify the built-in differences and com-plexities between the gene expressions of the two types of leukemia.786.3. Real data analysis0.0 0.2 0.4 0.6 0.8 1.0050100150Proportion pArea A(p)Sample ASample BFigure 6.13: Scale curves for ALL and AML samples under the extendedMahalanobis depth, showing that the ALL group has a larger variability.Our findings indicate that the proposed statistic M, which is calculatedusing the empirical quantile, provides a p-value of 0.0028. This outcomealigns with the results obtained from the scale curves, successfully differen-tiating between the two groups under study.In this chapter, we have proposed two new test statistics for testing ho-mogeneity of two multivariate samples. Different from other existing depth-based tests, our proposed test statistics were inspired by the quality index inthe context of quality assurance. Constructed based on the pairwise qualityindexes, our test statistics are shown to have the same χ2(1) asymptotic nulldistribution. Simulations studies demonstrate the superior performance ofthe proposed tests. The generalization of the proposed tests into the mul-tivariate multisample situation is discussed as well, along with some simu-lation studies for three sample comparisons. When the number of samplesgrows, the number of pairwise quality indexes increases. The improvementin the performance of the proposed tests and the higher-order approximationof their asymptotic distributions will be studied in future research.79Chapter 7ConclusionThis thesis conduct an in-depth exploration of the utility and robustness ofvarious statistical tools and procedures used in multivariate analysis. Werecognized that parametric methods like MANOVA, while widely adoptedand pivotal in the field, have inherent limitations, including assumptions ofmultivariate normality and homogeneity of covariance matrices, the com-plexity of results interpretation, and the method’s sensitivity to outliers.These constraints could lead to an inaccurate result. However, the evolutionof statistical methodology has provide us with more robust solutions. Therobust MANOVA methods, particularly those utilizing the Minimum Co-variance Determinant (MCD) estimator and rank orders, showed greater re-silience against outliers and violations of distributional assumptions, therebyenhancing the accuracy of the results.Our study also led to an appreciation for nonparametric tests, like theWilcoxon rank-sum statistics, which circumvent the rigorous assumptionsoften violated in parametric methods. However, due to the challenge ofextending the concept of ranking to multivariate data, the application ofthis method is mostly restricted to univariate cases.Consequently, we introduced depth functions as a method of ranking toovercome this limitation. A review of eleven depth functions was conducted,and their applicability and computational feasibility were analyzed. For el-liptical distributions, the Mahalanobis depth emerged as a preferred measuredue to its ability to account for the covariance structure of the data and itscomputational efficiency. However, in cases of outliers, a robust covarianceestimator such as the Minimum Covariance Determinant (MCD) is recom-mended. When missing coordinate values occur, a method for estimatinga lower limit on the zonoid depth can be used. This demonstrates the effi-ciency of the zonoid depth algorithm in handling large and high-dimensionaldatasets. Furthermore, in scenarios where the direction of an outlier is morepertinent than its depth, depth measures with a projection property, likethe zonoid depth, prove beneficial. The concept of depth function help togeneralize “ranking” in univariate non-parametric methods to the multivari-ate case. Two multivariate homogeneity test statistics based on data depth80Chapter 7. Conclusionwas then introduced. The one is the Kruskal-Wallis-type statistics H, theother is the quality index Q in the context of quality assurance.Our work extended to developing two new test statistics for testing ho-mogeneity of two multivariate samples, inspired by the quality index in thecontext of quality assurance. Unlike existing depth-based tests, our pro-posed test statistics were constructed based on the pairwise quality indexesand were shown to have the same χ2(1) asymptotic null distribution. Simu-lation studies demonstrate the superior performance of the proposed tests,and their generalization to multivariate multisample situations. However, asthe number of samples increases, so does the number of pairwise quality in-dexes, which will necessitate further research to improve the performance ofthe proposed tests and attain higher-order approximations of their asymp-totic distributions.The results of our research offer essential insights and open new opportu-nities for improving the robustness and accuracy of multivariate statisticalanalysis. These findings can lead to further investigations, particularly inchange point detection. When a single change point within a time seriesdataset is known to exist, the series can be iteratively divided into two sub-datasets. By performing the two-sample homogeneity test on these subsets,the division yielding the most significant difference is likely to indicate thelocation of the change point. This iterative approach can help in locatingthe change point in time series data.Despite the inherent limitations associated with any research study, ourwork makes a significant contribution to the field and sets the groundwork forfuture research. As such, we believe that continued exploration in this fieldwill yield even more robust and efficient statistical tools, thereby advancingthe field of multivariate analysis.81Bibliography[1] Shi, X., Zhang, Y., & Fu, Y. Two-Sample Tests Based on Data Depth.Entropy, 25(2), 238.[2] Tabachnick, B. G., Fidell, L. S., & Ullman, J. B. Using multivariatestatistics. Boston, MA: pearson 2013, 6, 497–516.[3] Quinn, G. P., & Keough, M. J. Experimental design and data analysisfor biologists. Cambridge university press 2002.[4] O’Brien, R. G., & Kaiser, M. K. Manova method for analyzing re-peated measures designs: an extensive primer. Psychological bulletin1985, 97(2), 316.[5] Zimmerman, D. W. Invalidation of parametric and nonparametric sta-tistical tests by concurrent violation of two assumptions. The Journalof Experimental Education 1998, 67(1), 55–68.[6] Snedecor, G. W., & Cochran, W. G. Statistical Methods, eight edition.Iowa state University press, Ames, Iowa 1989, 1191(2).[7] Fisher, R. A. On the interpretation of χ2 from contingency tables,and the calculation of P. Journal of the royal statistical society 1922,85(1), 87–94.[8] Laird, N. M., & Ware, J. H. Random-effects models for longitudinaldata. Biometrics 1982, 38(4), 963–974.[9] Hotelling, H. The generalization of Student’s ratio. In Breakthroughsin statistics: Foundations and basic theory. New York, NY: SpringerNew York 1992, 54-65.[10] Cole, D. A., Maxwell, S. E., Arvey, R., & Salas, E. Multivariategroup comparisons of variable systems: Manova and structural equa-tion modeling. Psychological Bulletin. 1993, 114(1), 174.82Bibliography[11] Fisher, R. A. Design of experiments. British Medical Journal. 1936,1(3923), 554.[12] Harold, H. Relations between two sets of variates. Biometrika 1936,28(3/4), 321.[13] Wilks, S. S. On the independence of k sets of normally distributedstatistical variables. Econometrica 1935, 3(3), 309–326.[14] Student. The probable error of a mean. Biometrika 1908, 6(1), 1–25.[15] Welch, B. L. The generalization of ‘STUDENT’S’problem when severaldifferent population varlances are involved. Biometrika 1947, 34(1-2),28–35.[16] Stevens, J. P. Applied multivariate statistics for the social sciences.Routledge 2012.[17] Hair, J. F., Black, W. C., Babin, B. J., Anderson, R. E., & Tatham,R. L. Multivariate data analysis 6th Edition 2006.[18] Meyers, L. S., Gamst, G., & Guarino, A. J. Applied multivariate re-search: Design and interpretation. Sage publications 2016.[19] Todorov, V., & Filzmoser, P. Robust Statistic for the One-wayManova. Computational Statistics & Data Analysis 2010, 54(1), 37–48.[20] Wilcoxon, F. Individual comparisons by ranking methods. SpringerNew York 1992, 196–202.[21] Woolson, R. F. Wilcoxon signed-rank test. Wiley encyclopedia of clin-ical trials 2007, 1–3.[22] Kruskal, W. A nonparametric test for the several sample problem. TheAnnals of Mathematical Statistics 1952, 23, 525–540.[23] Kruskal, W. and Wallis, A. Use of ranks in one-criterion varianceanalysis. Journal of the American Statistical Association 1952, 47,584–621.[24] Grawe, P. H. Nonparametric Statistics for the Behavioral Sciences.Numeracy: Advancing Education in Quantitative Literacy 2016 9(1).83Bibliography[25] Liu, R. Y. and Singh, K. A quality index based on data depth andmultivariate rank tests. Journal of the American Statistical Associa-tion 1993, 88(421), 252-260.[26] Baringhaus, L.; Franz, C. Fast nonparametric classification based ondata depth. J. Multivar. Anal. 2004, 88, 190–206.[27] Sze´kely, G.J.; Rizzo, M.L. Testing for equal distributions in high di-mension. InterStat 2004, 5, 1–6.[28] Ramdas, A.; Trillos, N.G.; Cuturi, M. On Wasserstein two-sampletesting and related families of nonparametric tests. Entropy 2017, 19,47.[29] Lumley, T., et al. The importance of the normality assumption in largepublic health data sets. Annual review of public health 2002, 23(1),151–169.[30] Button, K. S., Ioannidis, J. P., Mokrysz, C., Nosek, B. A., Flint, J.,Robinson, E. S., & Munafo`, M. R. Power failure: why small sample sizeundermines the reliability of neuroscience.Nature reviews neuroscience2013, 14(5), 365–376.[31] Hernandez, P. A., Graham, C. H., Master, L. L., & Albert, D. L.The effect of sample size and species characteristics on performanceof different species distribution modeling methods. Ecography 2006,29(5), 773–785.[32] Winship, C., & Mare, R. D. Regression models with ordinal variables.American sociological review 1984, 512–525.[33] Balke, N. S., & Fomby, T. B. Large shocks, small shocks, and economicfluctuations: outliers in macroeconomic time series. Journal of AppliedEconometrics 1994, 9(2), 181–200.[34] Conover, W. J. Practical nonparametric statistics. john wiley & sons1999, 350.[35] Rosner, B., Glynn, R. J., & Lee, M. L. T. The Wilcoxon signed ranktest for paired comparisons of clustered data. Biometrics 2006, 62(1),185–192.[36] Britto, C., Pollard, A. J., Voysey, M., & Blohmke, C. J. An appraisalof the clinical features of pediatric enteric fever: systematic review84Bibliographyand meta-analysis of the age-stratified disease occurrence. Clinical In-fectious Diseases 2017, 64(11), 1604–1611.[37] Costantino, A., Rinaldi, V., Moffa, A., Luccarelli, V., Bressi, F., Cas-sano, M., & Baptista, P. Hypoglossal nerve stimulation long-termclinical outcomes: a systematic review and meta-analysis. Sleep andBreathing 2020, 24, 399–411.[38] Jouneau, S., Crestani, B., Thibault, R., Lederlin, M., Vernhet, L.,Yang, M., & Cottin, V. Post hoc analysis of clinical outcomes inplacebo-and pirfenidone-treated patients with IPF stratified by BMIand weight loss. Respiration 2020, 101(2), 142–154.[39] Rosner, B., & Grove, D. Use of the Mann–Whitney U-test for clustereddata. Statistics in medicine 1999, 18(11), 1387–1400.[40] Rosner, B., Glynn, R. J., & Ting Lee, M. L. Incorporation of cluster-ing effects for the Wilcoxon rank-sum test: a large-sample approach.Biometrics 2003, 59(4), 1089–1098.[41] Datta, S., & Satten, G. A. Rank-sum tests for clustered data. Journalof the American Statistical Association 2005, 100(471), 908–915.[42] Dutta, S., & Datta, S. A rank-sum test for clustered data when thenumber of subjects in a group within a cluster is informative. Biomet-rics 2016, 72(2), 432–440.[43] Pocock, S. J. The combination of randomized and historical controlsin clinical trials. Journal of chronic diseases 1976, 29(3), 175–188.[44] Nath, R., & Pavur, R. A new statistic in the one-way multivariateanalysis of variance. Computational Statistics & Data Analysis 1985,2(4), 297–315.[45] Rousseeuw, P.J. Multivariate estimation with high-breakdown point.Mathematical Statistics and Applications 1985, 283-–297.[46] Chenouri, S. and Small, C. G. A nonparametric multivariate multisam-ple test based on data depth. Electronic Journal of Statistics 2012, 6,760–782.[47] Tukey, J. W. Mathematics and the picturing of data. Canadian Math-ematical Congress 1974, 2, 523–531.85Bibliography[73] Cuesta-Albertos, J.A. & Nieto-Reyes, A. The random Tukey depth.Computational Statistics and Data Analysis 2008, 52, 4979–4988.[49] Dyckerhoff, R. & Mozharovskyi, P. Exact computation of the halfspacedepth. Computational Statistics and Data Analysis 2016, 98, 19–30.[50] Mahalanobis, P. C. On the generalized distances in statistics: Maha-lanobis distance. Journal Soc. Bengal 1936, 26, 541–588.[51] Rousseeuw, P.J.; Leroy, A.M. Robust Regression and Outlier Detec-tion; John Wiley & Sons: New York, NY, USA, 1987.[52] Lopuhaa, H.P.; Rousseeuw, P.J. Breakdown points of affine equivari-ant estimators of multivariate location and covariance matrices. Ann.Stat. 1991, 19, 229–248.[53] Zou, Y. and He, X. One the limiting distributions of multivariatedepth-based rank-sum statistics and related tests. The Annals ofStatistics 2006, 24(6), 2879–2896.[54] Brown,M., B. Statistical use of spatial median. J.Roy.Statist.Soc.1958, 53, 448–456.[55] Gower, C.,J. Algorithm as 78: The mediancentre. App.Statist. 1974,23, 466–470.[56] Zou, Y. and Serfling, R. General notions of statistical depth function.The Annals of Statistics 2000, 28, 461-482.[57] Liu, R. Y. Data depth and multivariate rank tests. In L1-Statisticsand Related Methods (Y. Dodge, ed.) 1992, 279-294.[58] Liu, R. Y. On a notion of data depth based on random simplices. TheAnnals of Statistics 1990, 405–414.[59] Oja, H. Descriptive statistics for multivariate distributions. Statistics& Probability Letters 1983, 1(6), 327–332.[60] Koshevoy, G., & Mosler, K. zonoid trimming for multivariate distri-butions. The Annals of Statistics 1997, 25(5), 1998–2017.[61] Serfling, R. A depth function and a scale curve based on spatial quan-tiles. Statistics in Industry and Technology: Statistical Data Analysis2002, 25–3886Bibliography[62] Liu, Z., & Modarres, R. Lens data depth and median. Journal ofNonparametric Statistics 2011, 23(4), 1063–1074.[63] Yang, M., & Modarres, R. β-Skeleton depth functions and medi-ans. Communications in Statistics-Theory and Methods 2018, 47(20),5127–5143.[64] Eddy, W. F. Graphics for the multivariate two-sample problem: Com-ment. Journal of the American Statistical Association 1981, 76(374),287–289.[65] Mosler, K. and Mozharovskyi, P. Choosing among notions of multi-variate depth statistics. Statistical Science, 2022 37, 348-368.[66] Incropera, F. P., DeWitt, D. P., Bergman, T. L., & Lavine, A. S.Fundamentals of heat and mass transfer New York: Wiley 1996, 6,408–409.[67] Biemond, J., van der Putten, F. G., & Woods, J. W. Identification andrestoration of images with symmetric noncausal blurs. IEEE Transac-tions on Circuits and Systems 1988, 35(4), 385–393.[68] Kaiser, J., & Hamming, R. Sharpening the response of a symmetricnonrecursive filter by multiple use of the same filter. IEEE Transac-tions on Acoustics, Speech, and Signal Processing 1977, 25(5), 415–422.[69] Meer, P., & Georgescu, B. Edge detection with embedded confidence.IEEE Transactions on pattern analysis and machine intelligence 2001,23(12), 1351–1365.[70] Palumbo, D. C., Wong, G. N., & Prather, B. S. Discriminating accre-tion states via rotational symmetry in simulated polarimetric imagesof m87. The Astrophysical Journal 2020, 894(2), 156.[71] Russell, J. A. A circumplex model of affect. Journal of personality andsocial psychology 1980, 39(6), 1161.[72] McNeil, A. J., Frey, R., & Embrechts, P. Quantitative risk manage-ment: concepts, techniques and tools-revised edition. Princeton uni-versity press 2015.[73] Cuesta-Albertos, J. A., & Nieto-Reyes, A. The random Tukey depth.Computational Statistics Data Analysis 2008, 52(11), 4979–4988.87Bibliography[74] Nagy, S. The halfspace depth characterization problem. In Conferenceof the International Society for Non-Parametric Statistics 2018, 379–389.[75] Zuo, Y. A new approach for the computation of halfspace depth inhigh dimensions. Communications in Statistics-Simulation and Com-putation 2019, 48(3), 900–921.[76] Liu,R. Y., Jesse, M. P. and Kesar, S. Multivariate analysis by datadepth: Descriptive statistics, graphics and inference. The Annals ofStatistics 1999, 783–858.[77] Koshevoy, G. A. Projections of lift zonoids, the Oja depth and Tukeydepth. Preprint 2001.[78] Rousseeuw, P. J., & Ruts, I. The depth function of a population dis-tribution. Metrika 1999, 49(3), 213–244.[79] Chenouri, S. Multivariate robust nonparametric inference based ondata depth [Ph. D thesis]. Canada: Dept. of Statistics & ActuarialScience, University of Waterloo 2004.[80] Chenouri, S., Small, C. G., & Farrar, T. J. Data depth-based non-parametric scale tests. Canadian Journal of Statistics 2011, 39(2),356–369.[81] Liu, R. Y., Serfling, R. J., & Souvaine, D. L. (Eds.). Data depth: ro-bust multivariate analysis, computational geometry, and applications.American Mathematical Soc. 2006, 72.[82] Lo´pez-Pintado, S., Romo, J., & Torrente, A. Robust depth-based toolsfor the analysis of gene expression data. Biostatistics 2010, 11(2),254–264.[83] MacKay, D. J. Information theory, inference and learning algorithms.Cambridge university press.[84] Bu¨hlmann, P., & Van De Geer, S. Statistics for high-dimensional data:methods, theory and applications. Springer Science & Business Media2011.[85] Donoho, D. L. High-dimensional data analysis: The curses and bless-ings of dimensionality. AMS math challenges lecture 2000, 1(2000),32.88[86] Chen, S. X., & Qin, Y. L. A two-sample test for high-dimensional datawith applications to gene-set testing. The Annals of Statistics 2010,38(2), 808–835.[87] Chen, Y., Wiesel, A., Eldar, Y. C., & Hero, A. O. Shrinkage algo-rithms for MMSE covariance estimation. IEEE transactions on signalprocessing 2010, 58(10), 5016–5029.[88] Wang, C., Du, J., & Fan, X. High-dimensional correlation matrix esti-mation for general continuous data with Bagging technique. MachineLearning 2022, 111(8), 2905–2927.[89] Rousson, V. On distribution-free tests for the multivariate two-samplelocation-scale model. Journal of multivariate analysis, 2002, 80, 43-57.[90] Seber, G. A. F. Linear Regression Analysis. Wiley, New York. 1977,MR0436482.89Appendix AR codeA.1 Parametric test of the multivariatehomogeneity testsA.1.1 Table 3.2 Classical and robust MANOVA comparisonlibrary(rrcov)set.seed(1)mu1 = c(0,0)mu2 = c(0.2,0.2)sd1 = matrix(c(1,0,0,1), nc = 2)# Generate the datasetgroup1 <- matrix(rnorm(100, mean = mu1, sd = sd1), ncol = 2)group2 <- matrix(rnorm(100, mean = mu2, sd = sd1), ncol = 2)# Combine the groups into a single datasetdata <- rbind(group1, group2)group <- factor(rep(1:2, each = 50))# Add an outlier to the datasetdata <- rbind(data, c(10, 10))group <- factor(c(group, 2))# Perform classical MANOVAclassical_manova <- Wilks.test(data, grouping = group)# Perform robust MANOVArobust_manova <- Wilks.test(data, grouping = group, method = "rank")90A.2. Non-parametric univariate homogeneity tests# Print the resultsprint(classical_manova)print(robust_manova)A.2 Non-parametric univariate homogeneitytestsA.2.1 Sec 4.3 Wilcoxon rank-sum test real data analysis# install packages if not already installedif (!require(tidyverse)) install.packages("tidyverse")if (!require(sanon)) install.packages("sanon")if (!require(coin)) install.packages("coin")# load packagelibrary(tidyverse)library(sanon)library(coin)data("resp", package = "sanon")# Reshape the data to long formatresp_long <- resp %>%pivot_longer(cols = starts_with("visit"), names_to = "visit", values_to = "global_rating")# Convert character variables to appropriate classesresp_long$visit <- as.factor(resp_long$visit)resp_long$treatment <- as.factor(resp_long$treatment)resp_long$global_rating <- as.numeric(resp_long$global_rating)## Stratified Wilcoxon rank-sum test# List to store the resultsresults_list <- list()# Unique visitsvisits <- unique(resp_long$visit)91A.2. Non-parametric univariate homogeneity tests# Run the test for each visitfor (v in visits) {# Subset the data for the specific visitdata_subset <- subset(resp_long, visit == v)# Perform the testresult <- wilcox_test(global_rating ˜ treatment | center,data = data_subset)# Store the resultresults_listpaste("visit", v) <- result}print(results_list)## Standard Wilcoxon rank-sum test# List to hold the resultsresults_unstratified <- list()# Run the test for each visitfor (i in 1:4) {# Subset the data for the current visitdata_sub <- subset(resp_long, visit == paste0("visit", i))# Run the test without stratificationresult <- wilcox.test(global_rating ˜ treatment, data =data_sub)# Store the result in the listresults_unstratifiedpaste0("visit", i) <- result}print(results_unstratified)92A.3. Non-parametric multivariate homogeneity testsA.3 Non-parametric multivariate homogeneitytestsA.3.1 Fig 5.3 Half space depth demo# Points coordinatesx <- c(-0.5,1,4,2,1,2)y <- c(5,3.5,6,3,1,4)# Half-space depth pointpt_x <- 1pt_y <- 3.5postscript(file = paste0("/.../Plot/Plot_DepthDemo/Demo_HalfSpaceDepth.eps"),width = 570/72, height = 220/72, horizontal = FALSE,paper = "special") # width and height are ininchespar(mfrow = c(1, 3))# First plot (Vertical Line)plot(x, y, pch = 19, xlim = c(-1, 6), ylim = c(0, 7), xlab = "x", ylab = "y")points(pt_x, pt_y, col = "red", pch = 19)abline(v = pt_x, col = "blue")title("Vertical Line")# Second plot (Horizontal Line)plot(x, y, pch = 19, xlim = c(-1, 6), ylim = c(0, 7), xlab = "x", ylab = "y")points(pt_x, pt_y, col = "red", pch = 19)abline(h = pt_y, col = "green")title("Horizontal Line")# Third plot (Bisector Line)plot(x, y, pch = 19, xlim = c(-1, 6), ylim = c(0, 7), xlab = "x", ylab = "y")points(pt_x, pt_y, col = "red", pch = 19)abline(a = pt_y - pt_x, b = 1, col = "orange") # Bisector linehas slope 1 and intercept determined by the point (1,3.5)93A.3. Non-parametric multivariate homogeneity teststitle("Bisector Line")# Close the postscript devicedev.off()A.3.2 Fig 5.6 Projection depth demo# Install and load scatterplot3d packageif (!require(scatterplot3d)) {install.packages("scatterplot3d")library(scatterplot3d)}# Define the pointspoints3d = matrix(c(1, 2, 3, 2, 3, 4, 3, 2, 1, 4, 4, 4), ncol= 3)postscript(file = paste0("/.../Plot/Plot_DepthDemo/Demo_ProjectionDepth.eps"),width = 570/72, height = 500/72, horizontal = FALSE,paper = "special") # width and height are ininchespar(mfrow = c(1, 1))# Create the 3D scatter plots3d <- scatterplot3d(points3d, pch = 19, color = c("red", "green", "blue", "purple"),xlim = c(0,5), ylim = c(0,5), zlim = c(0,5),xlab = "X", ylab = "Y", zlab = "Z", box =FALSE)# Add text to each point to identify themlabels = c("A", "B", "C", "D")for (i in 1:4) {text(s3d$xyz.convert(points3d[i,1], points3d[i,2], points3d[i,3]), labels[i], pos=4)}94A.3. Non-parametric multivariate homogeneity tests# Add projection lines onto x-axis for each pointfor (i in 1:4) {s3d$points3d(c(points3d[i, 1], points3d[i, 1]),c(points3d[i, 2], 0),c(points3d[i, 3], 0),type = "l", lty = "dashed", col = "black")}# Close the postscript devicedev.off()A.3.3 Fig 5.8 Simplicial depth demolibrary(ggplot2)# Create the data frame of pointspoints <- data.frame(x = c(-0.5, 1, 4, 2, 1, 2),y = c(5, 3.5, 6, 3, 1, 4),label = c("A", "B", "C", "D", "E", "F"))# Define the point of interestpoint_of_interest <- data.frame(x = 1, y = 3.5, label = "P")postscript(file = paste0("/.../Plot/Plot_DepthDemo/Demo_SimplicialDepth.eps"),width = 570/72, height = 500/72, horizontal = FALSE,paper = "special") # width and height are ininchespar(mfrow = c(1, 1))library(ggplot2)# Create the data frame of pointspoints <- data.frame(x = c(-0.5, 1, 4, 2, 1, 2),y = c(5, 3.5, 6, 3, 1, 4),label = c("A", "B", "C", "D", "E", "F"))95A.3. Non-parametric multivariate homogeneity tests# Define the point of interestpoint_of_interest <- data.frame(x = 1, y = 3.5, label = "P")# Plot the pointsp <- ggplot() +geom_point(data = points, aes(x, y), size = 3) +geom_text(data = points, aes(x, y, label = label), vjust =-1)# Add lines to x and y axesp <- p + geom_vline(xintercept = -1, color = "black", size =0.5) +geom_hline(yintercept = 0, color = "black", size = 0.5)# Add the trianglescombinations <- combn(nrow(points), 3) # all possiblecombinations of 3 pointsfor(i in seq(ncol(combinations))) {subset_points <- points[combinations[, i], ]subset_points <- rbind(subset_points, subset_points[1, ]) #close the triangleif("D" %in% subset_points$label) { # check if the point ofinterest is in the hullp <- p + geom_path(data = subset_points, aes(x, y), color ="blue")} else {p <- p + geom_path(data = subset_points, aes(x, y), color ="blue")}}# Increase the size of x and y axis numbersp <- p + theme(axis.text = element_text(size = 14), panel.background = element_rect(fill = "white"))print(p)# Close the postscript device96A.3. Non-parametric multivariate homogeneity testsdev.off()A.3.4 Contour plots for the eleven depth functions withthe example 6-points dataset# Load the necessary packageslibrary(ddalpha)library(grDevices)# Define datadata <- matrix(c(-0.5,5,1,3.5,4,6,2,3,1,1,2,4), ncol = 2,byrow = TRUE)pt <- c(1, 3.5)# List of depth functionsdepths <- c("betaSkeleton", "halfspace", "L2", "Mahalanobis","projection", "simplicial", "spatial", "zonoid", "potential", "qhpeeling")# Generate a grid of pointsx <- seq(min(data[, 1]), max(data[, 1]), length.out = 100)y <- seq(min(data[, 2]), max(data[, 2]), length.out = 100)grid <- expand.grid(x = x, y = y)# Loop over each depth functionfor (depth_fun in depths) {# Calculate depth value for the pointdepth_val <- do.call(paste0("depth.", depth_fun), list(x =pt, data = data))print(paste("Depth value for", depth_fun, "is", depth_val))# Calculate depth values for each grid pointdepth_grid <- apply(grid, 1, function(pt) do.call(paste0("depth.", depth_fun), list(x = pt, data = data)))# Open a postscript devicepostscript(file = paste0("/.../Plot/Plot_DepthContour/Con",depth_fun, ".eps"),97A.3. Non-parametric multivariate homogeneity testswidth = 430/72, height = 400/72, horizontal =FALSE, paper = "special") # width and heightare in inches# Draw colored contour plotcontour(x, y, matrix(depth_grid, nrow = length(x)), col =rainbow(10),lwd = 3,labcex = 1.5)# Add legendlegend("topright", legend = pretty(range(depth_grid)), fill= rainbow(10), title = "Depth values")# Add a note indicating the sample sizemtext(paste("Sample size: ", nrow(data)), side = 3, line =-2, adj = 0)# Close the postscript devicedev.off()}A.3.5 Contour plots for the eleven depth functions with100 standard normal bivariate random variable# Load the necessary packageslibrary(ddalpha)library(grDevices)# Define dataset.seed(123) # For reproducibilitydata <- matrix(rnorm(200), ncol = 2)# List of depth functionsdepths <- c("betaSkeleton", "halfspace", "L2", "Mahalanobis","projection", "simplicial", "spatial", "zonoid", "potential", "qhpeeling")# Generate a grid of pointsx <- seq(min(data[, 1]), max(data[, 1]), length.out = 100)y <- seq(min(data[, 2]), max(data[, 2]), length.out = 100)98A.3. Non-parametric multivariate homogeneity testsgrid <- expand.grid(x = x, y = y)# Create directory "DepthContourPlot" if it doesn’t existdir.create(file.path("...", "DepthContourPlot"), showWarnings= FALSE)# Loop over each depth functionfor (depth_fun in depths) {# Calculate depth value for a random point in the datasetpt <- data[sample(nrow(data), 1), ]depth_val <- do.call(paste0("depth.", depth_fun), list(x =pt, data = data))# Calculate depth values for each grid pointdepth_grid <- apply(grid, 1, function(pt) do.call(paste0("depth.", depth_fun), list(x = pt, data = data)))# Open a postscript devicepostscript(file = paste0(".../Con100", depth_fun, ".eps"),width = 430/72, height = 400/72, horizontal =FALSE, paper = "special") # width and heightare in inches# Draw colored contour plotcontour(x, y, matrix(depth_grid, nrow = length(x)), col =rainbow(10),lwd = 3,labcex = 1.5)# Add legendlegend("topright", legend = pretty(range(depth_grid)), fill= rainbow(10), title = "Depth values")# Add a note indicating the sample sizemtext(paste("Sample size: ", nrow(data)), side = 3, line =-2, adj = 0)# Close the postscript devicedev.off()}99A.3. Non-parametric multivariate homogeneity testsA.3.6 Fig 5.16 Simulated examples of different type of datasymmetry# Set up the layoutlayout(matrix(1:3, ncol = 3))# Central symmetryset.seed(123)drug_effect <- rnorm(100, 0, 1)data <- c(drug_effect, -drug_effect)plot(data, ylim=c(-3,3), main="Central Symmetry", ylab="Drugeffect")# Angular symmetryset.seed(123)angle <- runif(100, -pi, pi)radius <- rexp(100, rate = 1)x <- radius*cos(angle)y <- radius*sin(angle)plot(x, y, main="Angular Symmetry")# Elliptical symmetryset.seed(123)library(MASS)n <- 80mu <- c(170, 70) # mean height and weightSigma <- matrix(c(100, 80, 80, 100), nrow=2)data <- mvrnorm(n, mu, Sigma)plot(data, xlab="Height", ylab="Weight", main="EllipticalSymmetry")A.3.7 Fig 5.17 Comparison of Mahalanobis depth contoursand density# Set the seedset.seed(123)# Import required packageslibrary(MASS)100A.4. A novel depth-based pairwise quality indexlibrary(ddalpha)# Generate multivariate normally distributed datan <- 500mu <- c(0, 0) # meanSigma <- matrix(c(1, 0.8, 0.8, 1), nrow=2) # covariance matrixwith strong correlationdata <- mvrnorm(n, mu, Sigma)# Calculate Mahalanobis depthsdepths <- depth.Mahalanobis(as.matrix(data), as.matrix(data))# Calculate 2D density estimatedensity_estimate <- kde2d(data[,1], data[,2], n=50)# Plot density contourscontour(density_estimate, xlab="X", ylab="Y", lwd=1, cex.axis=1)# Add the data points, colored by depthpoints(data, pch=19, col=rainbow(10)[as.numeric(cut(depths,10))], cex=0.3, cex.axis=1.5)A.4 A novel depth-based pairwise quality indexA.4.1 Fig 6.1 Type I errors for the two-sample casem=n=c(1:20)*50Type.Depth=c(1,2,3)Q=array(NA,dim=c(length(Type.Depth),length(m),5))for(iT in 1:length(Type.Depth)){for(im in 1:length(m)){my_data=paste(".../null",Type.Depth[iT],m[im],n[im], ".rds",sep = "-")five_table=readRDS(file =my_data)Q[iT,im,] <- apply(five_table,2,function (x) quantile(x,probs= 0.95))}}101A.4. A novel depth-based pairwise quality indexpar(mfrow=c(3,2) )par(mar=c(4,4,0,0)+0.1,fig=c(0,0.5,2/3,1) )plot(m,Q[1,,3],ylim=c(min(Q[,,3]),max(Q[,,3])),type=’l’,col=’black’,lwd=2,ylab="Critical value for equal weights",xlab="m (n=m)")abline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,3],lty=2,col="blue",lwd=2)lines(m,Q[3,,3],lty=3,col="purple",lwd=2)legend("topright",c("Mahalanobis","Spatial","Projection"),lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)par(mar=c(4,4,0,0)+0.1,fig=c(0,0.5,1/3,2/3) , new=TRUE)plot(m,Q[1,,5],ylim=c(min(Q[,,5]),max(Q[,,5])),type=’l’,col=’black’,lwd=2,ylab="Critical value for unequal weights",xlab="m (n=m)")abline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,5],lty=2,col="blue",lwd=2)lines(m,Q[3,,5],lty=3,col="purple",lwd=2)legend("topright",c("Mahalanobis","Spatial","Projection"),lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)par(mar=c(4,4,0,0)+0.1,fig=c(0,0.5,0,1/3) , new=TRUE)plot(m,Q[1,,4],ylim=c(qchisq(df=1,0.95),max(Q[,,4])),type=’l’,col=’black’,lwd=2,ylab="Critical value for maximum",xlab="m (n=m)")abline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,4],lty=2,col="blue",lwd=2)lines(m,Q[3,,4],lty=3,col="purple",lwd=2)legend("topright",102A.4. A novel depth-based pairwise quality indexc("Mahalanobis","Spatial","Projection"),lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)m=c(1:20)*50n=m/2Type.Depth=c(1,2,3)Q=array(NA,dim=c(length(Type.Depth),length(m),5))for(iT in 1:length(Type.Depth)){for(im in 1:length(m)){my_data=paste(".../null",Type.Depth[iT],m[im],n[im], ".rds",sep = "-")five_table=readRDS(file =my_data)Q[iT,im,] <- apply(five_table,2,function (x) quantile(x,probs= 0.95))}}par(mar=c(4,4,0,0)+0.1,fig=c(0.5,1,2/3,1), new=TRUE)plot(m,Q[1,,3],ylim=c(min(Q[,,3]),max(Q[,,3])),type=’l’,col=’black’,lwd=2,ylab="Critical value for equal weights",xlab="m (n=m/2)")abline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,3],lty=2,col="blue",lwd=2)lines(m,Q[3,,3],lty=3,col="purple",lwd=2)legend("topright",c("Mahalanobis","Spatial","Projection"),lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)par(mar=c(4,4,0,0)+0.1,fig=c(0.5,1,1/3,2/3), new=TRUE)plot(m,Q[1,,5],ylim=c(min(Q[,,5]),max(Q[,,5])),type=’l’,col=’black’,lwd=2,ylab="Critical value for unequal weights",xlab="m (n=m/2)")103A.4. A novel depth-based pairwise quality indexabline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,5],lty=2,col="blue",lwd=2)lines(m,Q[3,,5],lty=3,col="purple",lwd=2)legend("topright",c("Mahalanobis","Spatial","Projection"),lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)par(mar=c(4,4,0,0)+0.1,fig=c(0.5,1,0,1/3), new=TRUE)plot(m,Q[1,,4],ylim=c(qchisq(df=1,0.95),max(Q[,,4])),type=’l’,col=’black’,lwd=2,ylab="Critical value for maximum",xlab="m (n=m/2)")abline(h=qchisq(df=1,0.95),col="red",lwd=2)lines(m,Q[2,,4],lty=2,col="blue",lwd=2)lines(m,Q[3,,4],lty=3,col="purple",lwd=2)legend("topright",c("Mahalanobis","Spatial","Projection"), lty=c(1,2,3),col=c("black","blue","purple"),bty=’n’,lwd=2)A.4.2 Fig 6.5, 6.6, 6.7, 6.2, 6.3, 6.4 Two-sample powercomparisonsThe numerical simulations of the two-sample power comparisons are orga-nized at Cassie’s GitHub.A.4.3 Fig 6.8, 6.9 Three-sample power comparisons# Empirical Distributionlibrary(DepthProc)set.seed(1)mu1 <- c(0,0)mu2<-c(0,0)mu3<-c(0,0)sigma1 <- matrix(c(1,0,0,1), nc = 2)sigma2 <- matrix(c(1,0,0,1), nc = 2)104A.4. A novel depth-based pairwise quality indexsigma3 <- matrix(c(1,0,0,1), nc = 2)Loop=10000m=n=k=c(1:10)*100Type.Depth=c(1,2,3)for(iT in 1:length(Type.Depth)){for(im in 1:length(m)){MV=c()for (i in 1:Loop) {Fm <- rbind(mvrnorm(m[im], mu1 ,sigma1))Gn <- rbind(mvrnorm(n[im], mu2 ,sigma2))Hk <- rbind(mvrnorm(k[im], mu3 ,sigma3))MV=rbind(MV,c(max_three(Fm,Gn,Hk,Type.Depth[iT]),H_test_three(Fm,Gn,Hk,Type.Depth[iT]) ))}# Save an object to a filemy_data=paste(".../null",Type.Depth[iT],m[im],n[im],k[im],".rds",sep = "-")saveRDS(MV, file = my_data)}}set.seed(1)mu1 <- c(0,0)mu2 <- c(0,0)sigma1 <- matrix(c(1,0,0,1), nc = 2)sigma2 <- matrix(c(1,0,0,1), nc = 2)Loop=100000m=n=c(1:10)*100Type.Depth=c(1,2,3)for(iT in 1:length(Type.Depth)){for(im in 1:length(m)){h=c()for (i in 1:Loop) {Fm <- rbind(mvrnorm(m[im], mu1 ,sigma1))Gn <- rbind(mvrnorm(n[im], mu2 ,sigma2))h=rbind(h,H_test(Fm,Gn,Type.Depth[iT]))}colnames(h)=c(’DbR’)# Save an object to a filemy_data=paste(".../null",105A.4. A novel depth-based pairwise quality indexType.Depth[iT],m[im],n[im], "DbR.rds",sep = "-")saveRDS(h, file = my_data)}}A.4.4 Power Quantificationlibrary(entropy)#------------mean/nonequal-sizes/Maha------------power_curve <- c(0.215,0.244,0.332,0.435,0.527,0.595,0.624,0.724,0.770,0.849)# Mp <- power_curve / sum(power_curve)#sum(p)power_curve <- c(0.134,0.178,0.191,0.238,0.272,0.331,0.348,0.409,0.448,0.521)# W(1,0)q <- power_curve / sum(power_curve)#sum(q)# Compute KL divergence from p to q using entropy packagekl_div_entropy <- KL.plugin(p, q)print(kl_div_entropy)# Compute the average differenceaverage_difference <- mean(abs(p - q))print(average_difference)#------------variance/nonequal-sizes/Maha------------power_curve <- c(0.233,0.321,0.432,0.553,0.685,0.722,0.785,0.860,0.887,0.918)# Mp <- power_curve / sum(power_curve)#sum(p)106A.4. A novel depth-based pairwise quality indexpower_curve <- c(0.085,0.064,0.063,0.070,0.08,0.082,0.078,0.099,0.084,0.096)# W(1,0)q <- power_curve / sum(power_curve)#sum(q)# Compute KL divergence from p to q using entropy packagekl_div_entropy <- KL.plugin(p, q)print(kl_div_entropy)# Compute the average differenceaverage_difference <- mean(abs(p - q))print(average_difference)A.4.5 Fig 6.10 Simulation of high-dimension cases# Simulation# H0m=50n=30p=100mu1=mu2=rep(1,p) # Generate multivariate normal distributionsigma1=0.1ˆ(abs(outer(1:p, 1:p, "-")))sigma2=0.1ˆ(abs(outer(1:p, 1:p, "-")))result=compare=c()for (i in 1:5000) {set.seed(NULL)Fm <- rbind(mvrnorm(m, mu1 ,sigma1))Gn <- rbind(mvrnorm(n, mu2 ,sigma2))result=rbind(result,M_HighDim(Fm,Gn))}quantile_M <- quantile(result,probs = 0.95)# Haquantile_M <- 45.2929m=50n=30p=100107A.4. A novel depth-based pairwise quality indexmu1=mu2=rep(1,p) # Generate multivariate normal distributionpower_seq = c()for(i in 2:9){sigma1=(0.1*i)ˆ(abs(outer(1:p, 1:p, "-"))) #2-9 loop 0.1ksigma2=0.1ˆ(abs(outer(1:p, 1:p, "-")))result=compare=c()iter<-5000for (i in 1:iter) {set.seed(NULL)Fm <- rbind(mvrnorm(m, mu1 ,sigma1))Gn <- rbind(mvrnorm(n, mu2 ,sigma2))compare=rbind(compare,R.test(x=Fm,y=Gn,n=10,m=10,seed=i))result=rbind(result,M_HighDim(Fm,Gn))}power <- c(sum(result>=quantile_M)/iter,sum(compare<=0.05)/iter)power_seq<- rbind(power_seq,power)# The power comparison plot# Create datadf <- data.frame(x_s1 = c(2, 3, 4, 5, 6, 7, 8, 9, 2, 3, 4, 5, 6, 7, 8, 9),type_s1 = c(rep("M", 8), rep("B", 8)),value_s1 = c(0.0544, 0.1676, 0.8612, 1, 1, 1, 1, 1,0.043, 0.038, 0.0454, 0.0446, 0.047, 0.043,0.0496, 0.0664))# Separate data by typedata_M <- df[df$type_s1 == "M",]data_B <- df[df$type_s1 == "B",]# Create a blank plotplot(1, type="n", xlab="", ylab="", xlim=c(2,9), ylim=c(0,1.1),108A.4. A novel depth-based pairwise quality indexxaxt = "n", yaxt = "n") # Here xlim and ylim set the rangesfor the axes# Add minimal axesaxis(1, at=2:9)axis(2, at=seq(0,1.1, by=0.1))# Add lineslines(data_M$x_s1, data_M$value_s1, type="b", col="red", pch=4, lwd=2)lines(data_B$x_s1, data_B$value_s1, type="b", col="darkgreen",pch=15, lwd=2)# Add labelsxlabel <- "x values"ylabel <- "Values"xlab(xlabel)ylab(ylabel)# Add legendlegend("topleft", legend=c("M", "MBD"), col=c("red", "darkgreen"), pch=c(4,15), lwd=2)}A.4.6 Sec 6.3.3 Data analysis# Scale curve---------ScaleCurve.Maha=function (x, y = NULL, alpha = seq(0, 1, 0.01), name = "X", name_y = "Y", title = "Scale Curve"){x <- na.omit(x)if (is.data.frame(x)) {x <- as.matrix(x)}if (!is.matrix(x)) {stop("x must be a matrix or data frame!")}if (!is.null(y)) {if (is.data.frame(y)) {109A.4. A novel depth-based pairwise quality indexy <- as.matrix(y)}if (!is.matrix(y)) {stop("y must be a matrix or data frame!")}}dim_x <- dim(x)[2]uxname_list <- list(u = x, X = x)m=apply(x,2,mean)covar=solve(CovEst.2010OAS(x)$S)depth_est=NULLfor(j in 1:dim(x)[1]){tmpX = x[j,] - m;dist=t(tmpX)%*%covar%*%tmpXdepth_est[j]=1.0/(1.0 + dist[1])}k <- length(alpha)vol <- 1:kalpha_border <- ecdf(depth_est)(depth_est)for (i in 1:k) {if(sum(alpha_border >= alpha[i])>1) {tmp_x <- x[alpha_border>= alpha[i], ]xx=c(1:dim_x)AUC = trapz(xx,apply(tmp_x,2,max))AUC = AUC-trapz(xx,apply(tmp_x,2,min))vol[i]=AUC}else{ vol[i]=0}}return(list(p=alpha,vol=rev(vol)))}# Empirical quantile-----------m=47n=25p=100mu1=rep(0,p) # Generate multivariate normal distributionmu2=rep(0,p)sigma1=matrix(diag(rep(1,p)),ncol=p)sigma2=matrix(diag(rep(1,p)),ncol=p)110A.4. A novel depth-based pairwise quality indexresult=c()for (i in 1:5000) {set.seed(NULL)Fm <- rbind(mvrnorm(m, mu1 ,sigma1))Gn <- rbind(mvrnorm(n, mu2 ,sigma2))result=rbind(result,M_HighDim(Fm,Gn))}# P-value----------------data(leukemia)X=as.matrix(leukemia[which(leukemia[,102]=="M"),2:101])Y=as.matrix(leukemia[-which(leukemia[,102]=="M"),2:101])r1=Q(X,Y)r2=Q(Y,X)((1/dim(X)[1]+1/dim(Y)[1])*(1/12))ˆ(-1)*max(c((r1-0.5)ˆ2,(r2-0.5)ˆ2))pvlaue <- sum(result>=48.95833)/5000# Plot for Real Data ---------------------data(leukemia)X=as.matrix(leukemia[which(leukemia[,102]=="M"),2:101])Y=as.matrix(leukemia[-which(leukemia[,102]=="M"),2:101])Re1=ScaleCurve.Maha(X)Re2=ScaleCurve.Maha(Y)leu_type <- cbind(t(rep(0,101)),t(rep(1,101)))leu_x <- cbind(t(Re1$p),t(Re2$p))leu_y <- cbind(t(Re1$vol),t(Re2$vol))leu_result <- as.data.frame(cbind(t(leu_type),t(leu_x),t(leu_y)))ggplot(leu_result, aes(x=V2, y=V3, group=as.factor(V1))) +geom_line(aes(color=as.factor(V1)),size=1)+theme(legend.position="top")+labs(x="Proportion p", y = "Area A(p)")+guides(color=guide_legend("Samples"))+111A.4. A novel depth-based pairwise quality indextheme_stata() + scale_color_stata(labels = c("AML","ALL"))112Appendix BAn asymptotic normalitytheorem for the quality index[53]Denote the Liu–Singh depth-based rank-sum statistic as Q-statistic, whichis asymptotically distribution-free under the null hypothesis. As the samplesize increases, the test becomes more robust to deviations from the conjec-tured limiting distribution, leading to more reliable results and conclusions.An attempt was made to prove the conjecture [89]. However, this attemptwas not completely successful because it did not properly address the dif-ferentiability aspect of the depth functionals. This might have led to a lessrigorous proof.This chapter is a review of [53] which addresses the existing mathematicalgap by establishing regularity conditions for the limiting distribution andverifying these conditions for various commonly used depth functions. Inaddition, it proved the limiting behavior of these statistics. This researchcontributes to the understanding and robustness of depth-based methodsand their sensitivity to outliers and noise.B.1 AssumptionsThere are four assumptions [53] that are needed for the main result. The ver-ification of the assumptions is provided in later sections for some commonlyused depth functions. Assume, without loss of generality, that m ≤ n here-after. Fm is the empirical version of F and D(·; ·) is a given depth functionwith 0 ≤ D(x;H) ≤ 1 for any point x and distribution H in Rd.(A1) P{y1 ≤ D(Y ;F ) ≤ y2} ≤ C|y2−y1| for some C and any y1, y2 ∈ [0, 1].(A2) supx∈Rd |D(x;Fm)−D(x;F )| = o(1), with probability 1, as m→∞.(A3) E(supx∈Rd |D(x;Fm)−D(x;F )|) = O(m−1/2).113B.1. Assumptions(A4) E(∑i piX(Fm)piY (Fm)) = o(m−1/2) if there exist ci such thatpiX(Fm) > 0 and piY (Fm) > 0 forpiZ(Fm) := P (D(Z;Fm) = ci|Fm), i = 1, 2, ....(A1) imposes Lipschitz continuity on the depth function D(Y ;F ), whichis a property that bounds the rate of change of a function. It puts a bound onhow quickly the depth function can change concerning the probability mea-sure P . Specifically, the probability of the depth function D(Y ;F ) fallingwithin the interval [y1, y2] is bounded by a constant C times the length ofthe interval |y2 − y1|. Here, C is a positive constant, and y1 and y2 are anytwo points in the range [0, 1]. This Lipschitz-type condition helps controlthe behavior of the depth function, ensuring it has well-behaved propertiesnecessary for proving the asymptotic behavior of Q.(A2) and (A3) describe the behavior of the difference between the em-pirical depth function D(x;Fm) and the true depth function D(x;F ) as thesample size m goes to infinity. (A2) describes the pointwise convergenceof the empirical depth function to the true depth function. It states thatthe supremum (maximum) difference between the empirical and true depthfunctions over all points x in Rd converges to zero, almost surely (withprobability 1), as the sample size m goes to infinity. This implies that theempirical depth function is a good approximation of the true depth function.(A3) is a measure of the average behavior of the convergence. It describesthe expected value of the supremum difference between the empirical andtrue depth functions. The notation O(m−1/2) tells us that the functiongrows at most as fast as m−1/2. This provides information about the av-erage rate of convergence, helping us understand how fast the convergenceoccurs on average.In the context of (A4), the probability of observing a depth of ci fora random variable Z given the empirical distribution Fm is represented bypiZ(Fm) := P (D(Z;Fm) = ci|Fm). There are some depths ci that areobserved with positive probability for both X and Y. However, the expectedvalue of the sum converges to zero at a rate faster than m−1/2. Intuitively,this means that even when there are some depths ci that are observed withpositive probability for both X and Y, the overall sum of the product of theseprobabilities converges to zero as the sample sizem increases. Consequently,the expected contribution of ci becomes smaller and smaller as the samplesize grows.Theorem 2. [53] Let X ∼ F and Y ∼ G be independent and X1, ..., Xm114B.2. Op and op notationand Y1, ..., Yn be independent samples from F and G, respectively.(σ2GF /m+ σ2FG/n)−1/2(Q(Fm, Gn)−Q(F,G)) d−→ N(0, 1), as m→∞,whereσ2FG =∫P 2(D(X;F ) ≤ D(y, F ))dG(y)−Q2(F,G),σ2GF =∫P 2(D(x;F ) ≤ D(Y, F ))dF (x)−Q2(F,G).This equation can be understood as a Central Limit Theorem (CLT)-likeresult for the Q statistic. It shows that the difference between the empiricalQ statistic and the true, population-based Q statistic, properly scaled by theinverse of the square root of the sum of the variances of the depth functiondivided by the sample sizes, converges to a standard normal distributionas the sample size m goes to infinity. σ2FG and σ2GF are the variances ofthe depth function, under the distributions F and G, respectively. σ2FGis calculated by integrating the squared probabilities of the event that thedepth of a random point X from distribution F is less than or equal to thedepth of a fixed point y from the distribution G.B.2 Op and op notationIn probability theory, the notations op and Op are used to describe thebehavior of random variables or random sequences as the sample size orsome other variable goes to infinity while taking into account the randomnature of these variables. In contrast, o and O notations are used to describethe behavior of deterministic functions as a variable goes to infinity.The little-o of probability denoted op, is used to indicate a rate of con-vergence for random variables or random sequences. If Xn is a sequence ofrandom variables and an is a sequence of constants, we say that Xn = op(an)if, for any ϵ > 0, we havelimn→∞P (|Xn/an| > ϵ) = 0.This means that Xn converges to zero “faster” than an as n goes to infinity,which is a key aspect of little-o notation in probability. On the other hand,the big-O of probability, denoted Op, is used to indicate an order of mag-nitude of convergence. We say that Xn = Op(an) if there exists a positiveconstant C such thatlim supn→∞P (|Xn/an| > C) = 0.115B.3. Auxiliary lemmasThis means that the ratio Xn/an is bounded by the constant C in probabil-ity. In other words, Xn converges to zero “at the same rate” as an up to aconstant factor, which is a key aspect of big-O notation in probability.B.3 Auxiliary lemmasLemma 1 [53]. For any distribution functions H, F1 and F2 in Rd, pointsx and y in Rd and a given depth function D(·; ·),I(x, y,H) = I{D(x;H) ≤ D(y;H)},I(x, y, F1, F2) = I(x, y, F1)− I(x, y, F2).(i)∫∫I(x, y, F )d(Gn(y)−G(y))d(Fm(x)− F (x)) = Op(1/√mn),(ii)∫∫I(x, y, Fm, F )d(Fm(x)− F (x))dG(y) = op(1/√m),(iii)∫∫I(x, y, Fm, F )dFm(x)d(Gn −G).The first double integral quantifies the joint behaviour of the samplesX and Y concerning the depth relationship under the distribution F , whiletaking into account the deviations between empirical and true distributionfunctions. It tells us that as we collect more data, the empirical distributionswill become closer to the true distributions, and the depth relationshipsbetween points in X and Y will also converge, at a rate proportional to1/√mn in probability.The second double integral is the difference between these two indicatorfunctions. It captures the change in the depth relationship between x and ywhen switching from the empirical distribution Fm to the true distributionF . It indicates that the difference in the depth relationships between thepoints in the two samples when using the empirical distribution instead ofthe true distribution decreases as the sample size increases. As the samplesize grows, the empirical distribution becomes a better approximation of thetrue distribution and the depth relationships are better preserved.Lemma 2 [53]. Assume that X ∼ F and Y ∼ G are independent. Thenunder (A4), we have∫∫I(D(x;Fm) = D(y;Fm))dF (x)dG(y) = o(m1/2).The integral represents the joint probability of x and y having the samedepth value under the empirical distribution function Fm. Under assumption116B.4. Proof of Theorem 2(A4), this integral converges to zero at a faster rate than m1/2, i.e., o(m1/2).This means that the probability of observing a pair of points (x, y) fromthe independent distributions F and G with equal depth values under theempirical distribution function Fm becomes negligible as the sample size mincreases.Lemma 3 [53]. Let X ∼ F and Y ∼ G be independent and let X1, ..., Xmand Y1, ..., Ym be independent samples from F and G, respectively. Under(A1)-(A4),Q(Fm, Gn)−Q(F,Gn) =∫∫I(x, y, F )dG(y)d(Fm(x)− F (x)) + op(m−1/2)Consequently,√m(Q(Fm, Gn)−Q(F,Gn)) d−→ N(0, σ2GF ).Lemma 3 shows that the difference between the depth-based rank statis-tic calculated using the empirical distribution function Fm and the truedistribution function F can be decomposed into two parts: the double in-tegral term and the op(m−1/2) term. the op(m−1/2) term converges to zerofaster than m−1/2, leading to the result that the product of the square rootof m and the difference between the depth-based rank statistics convergesin distribution to a normal distribution with mean 0 and variance σ2GF .B.4 Proof of Theorem 2By lemma 3, we haveQ(Fm, Gn)−Q(F,G)= (Q(Fm, Gn)−Q(F,Gn)) + (Q(F,Gn)−Q(F,G))=∫∫I(x, y, F )dG(y)d(Fm(x)− F (x))+∫∫I(x, y, F )dF (x)d(Gn(y)−G(y)) + op(m−1/2).The op(m−1/2) term in the equation represents an error term that goes tozero as m goes to infinity. This term is negligible compared to the mainterm which is why it is written in o notation (as opposed to uppercase Onotation which denotes a term that dominates all other terms).This result decomposes the difference between the test statistics into twoterms: one that captures the bias due to using the empirical distribution117B.5. Verification of assumptions with Mahalanobis depth (MHD)function Fm instead of the true distribution function F , and another thatcaptures the bias due to using the empirical distribution function Gn insteadof the true distribution function G. The additional error term, op(m−1/2),ensures the consistency of the test statistic.Now, because the samples X1, . . . , Xm and Y1, . . . , Yn are independent,we can apply the central limit theorem to the empirical distribution functionsFm and Gn. This allows us to establish the asymptotic normality of thedifference (σ2GF /m+ σ2FG/n)−1/2(Q(Fm, Gn)−Q(F,G)) as m→∞.B.5 Verification of assumptions withMahalanobis depth (MHD)The depth of a point x is defined asMHD(x;F ) = 1/(1 + (x− µ)′Σ−1(x− µ)), x ∈ Rd,where µ and Σ−1 are location and covariance measures of a given distributionF . Fm is the empirical version of F , and µm and Σ−1m are strongly consistentestimators of µ and Σ−1.Note that if convergence is almost certain then the estimator is said tobe strongly consistent. Hence,supx∈Rd |MHD(x;Fm)−MHD(x;F )|= |MHD(xm;Fm)−MHD(xm;F )|.by the continuity of MHD(x;F ) and MHD(x;Fm) in x for some xm =x(Fm, F ) ∈ Rd such that ||xm|| ≤ M ≤ ∞ for some M > 0. Here, thenotation xm = x(Fm, F ) represents a point in Rd where the absolute differ-ence between the Mahalanobis depth (MHD) with respect to the empiricaldistribution Fm and the true distribution F is maximized. Then,|MHD(xm;Fm)−MHD(xm;F )|= |µm − µ|′Σ−1m (µm + µ− 2xm) + (xm − µ)′(Σ−1m − Σ−1)(xm − µ)|(1 + ||Σ−1/2m (xm − µm)||2)(1 + ||Σ−1/2(xm − µ)||2).With the strong consistency of µm and Σm, as m → ∞. In order to showthat the Mahalanobis depth (MHD) satisfies assumption (A2), we need toprove that the supremum of the absolute difference between MHD(x;Fm)and MHD(x;F ) converges to zero almost surely as m → ∞. Strong con-sistency of an estimator means that as the sample size m goes to infinity,118B.5. Verification of assumptions with Mahalanobis depth (MHD)the estimator converges to the true parameter almost surely. In this case,µm is a strongly consistent estimator of µ, and Σm is a strongly consistentestimator of Σ. This implies that as m → ∞, µm a.s.−−→ µ and Σm a.s.−−→ Σ.Because of the strong consistency, the following terms will converge to zeroalmost surely as m→∞:supx∈Rd |MHD(x;Fm)−MHD(x;F )| a.s.−−→ 0.This result demonstrates that the MHD satisfies assumption (A2):supx∈Rd |D(x;Fm)−D(x;F )| = o(1),almost surely, as m→∞.In the context of probability and statistics, Ho¨lder’s inequality is usefulfor bounding the expected value of the product of random variables. Giventwo random variables X and Y with finite expectations, and two real num-bers p and q such that 1 ≤ p, q < ∞ and 1p + 1q = 1, Ho¨lder’s inequalitystates that:E(|XY |) ≤ p√E(|X|p) · q√E(|Y |q).Ho¨lder’s inequality and expectations of quadratic forms (page 13 of [90])yield (A3) if conditions (i) and (ii) below are met.(A4) holds trivially if (iii) holds. (iii) states that the probability massof X at any ellipsoid is 0. Due to the continuous nature of the distributionand the ellipsoidal depth contours of MHD, Mahalanobis depth satisfies(A4), and the interaction between multiple depths sharing the same valuesis negligible, so the terms piX(Fm) and piY (Fm) are very small.(i) µm and Σm are strongly consistent estimators of µ and Σ, respectively.(ii) E(µm − µ)i = O(m−1/2), E(Σ−1m −Σ−1)jk = O(m−1/2), 1 ≤ i, j, k ≤ d,where the subscripts i and jk denote the elements of a vector and amatrix, respectively.(iii) The probability mass of X at any ellipsoid is 0.In this chapter, we delved into the Liu–Singh depth-based rank-sumstatistic, also known as the Q-statistic and its asymptotic behavior. Fourprimary assumptions (A1-A4) and three auxiliary lemmas were introducedand discussed to lay the groundwork for proving the asymptotic normalitytheorem for the quality index. The theorem states that the difference be-tween the empirical Q statistic and the true, population-based Q statistic,119B.5. Verification of assumptions with Mahalanobis depth (MHD)properly scaled by the inverse of the square root of the sum of the variancesof the depth function divided by the sample sizes, converges to a standardnormal distribution as the sample size goes to infinity. This chapter deeperthe understanding of the asymptotic behavior of the depth-based rank sumstatistics and promise more reliable results and conclusions in statisticalanalysis. The mathematical investigation conducted in this chapter sets thegroundwork for future research and discussions on the asymptotic behaviorof other depth-based methods.120
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Xiaoping, Shi
Ranking functions of data
Text