{'fraction_non_alphanumeric': 0.0350457984866587, 'fraction_numerical': 0.013504217805293074, 'mean_word_length': 4.267353165522502, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymptotic distribution of principal component scores for pervasive, high-dimensional eigenvectors

Plots of scores from principal component analysis are a popular approach to visualize and explore high-dimensional genetic data. However, the inconsistency of the high-dimensional eigenvectors has discredited classical principal component analysis and helped motivate sparse principal component analysis where the eigenvectors are regularized. Still, classical principal component analysis is extensively and successfully used for data visualization, and our aim is to give an explanation of this paradoxical situation. We show that the visual information given by the relative positions of the scores will be consistent, if the related signal can be considered to be pervasive. Firstly, we argue that pervasive signals lead to eigenvalues scaling linearly with the dimension, and we discuss genetic applications where such pervasive signals are reasonable. Secondly, we prove within the high-dimension low sample size regime, that when eigenvalues scale linearly with the dimension, the sample component scores will appear as scaled and rotated versions of the population scores. In consequence, the relative positions and visual information conveyed by the score plots will be consistent.


Introduction
Principal component analysis (PCA) is the workhorse of variable reduction in applied data analysis. It is used to construct a small number of informative scores from the original data, and these scores are then used further in visualization or in conventional classification, clustering or regression methods. This is highly useful in the context of modern highdimensional data analysis, where the number of measured variables p exceeds the sample size n. Genomics is an application area are where the first step in exploring data is often to visually investigate the first few principal component (PC) scores.
The asymptotic behavior of high-dimensional PCA has attracted a substantial amount of attention the last few years. It has been shown, by Paul (2007) and Johnstone and Lu (2009), that the population eigenvalues and -vectors in PCA are not consistently estimated by the sample eigenvalues and -vectors under the finite γ regime, where p/n = γ as p, n → ∞. The inconsistency of the eigenvectors is quantified in terms of the inner product between the sample and the population eigenvectors, which then does not converge to 1. In view of this, Johnstone and Lu (2009) suggest that one could either conduct an initial dimension reduction, from the original number of variables to a value less than n, before applying PCA, or introduce a sparse penalty on the eigenvectors, giving rise to the sparse PCA methodology (Witten et al., 2009;Zou et al., 2006). However, in an applied setting, the behavior of the principal component scores is also of interest, in addition to the eigenvectors and eigenvalues.
Until now, only few papers have focused on the asymptotic behavior of the principal component scores. An exception is Lee et al. (2010), who note that: "Inconsistency of the sample eigenvectors does not necessarily imply poor performance of PCA". The success of applied PC scores in genomics suggests that they can be considered suitable for an analysis, in spite of the inconsistency of the eigenvectors. In this paper, we explore this somewhat paradoxical situation further, and try to bridge the gap between the theoretical problems of PCA and the practical usefulness of the method.
We first review the main structure of high-dimensional PCA and earlier asymptotic results. Next, we introduce the concept of pervasive effects and demonstrate that this leads to population eigenvalues scaling linearly with the data dimension. Under this assumption about the eigenvalues, we derive the asymptotic limiting distribution of the ratio between the estimated and the true principal component scores. The implications of our findings are explored theoretically and by simulations. We show that all sample PC scores are subject to a common scaling and a small noise term, such that the relative positions of the sample and population scores are essentially unchanged by the scaling.
2 Principal component analysis 2.1 Methods and notation PCA reduces the data dimension by constructing orthogonal linear combinations of variables, which explain their variability. The first component is the normalized linear combination of variables with the highest variance, while the second component will be the linear combination, orthogonal to the first, with the highest variance, and so on. The mathematical basis of PCA is the eigendecomposition of the data covariance matrix. Let X = [x 1 , . . . , x n ] be a p × n data matrix, where x i = [x i1 , . . . , x ip ] T are independent and identically distributed with E x i = 0 and var x i = Σ. The eigendecomposition of the covariance matrix is given by where Λ is the diagonal matrix of the eigenvalues λ 1 ≥ · · · ≥ λ p and V = [v 1 , . . . , v p ] is the matrix of eigenvectors. The weights of the orthogonal linear combinations are given by the eigenvectors, usually referred to as loadings. We denote the vector of the resulting population component scores by The eigenvalues express the variance of the component scores so that the vector of standardized population component scores is given by where the jth vector of scores is z T j = [z j1 , . . . , z jn ] and Z = [z 1 , . . . , z p ] T . An applied data analysis is based on the sample covariance matrix denoted byΣ = 1 n XX T , with the eigendecompositionΣ =VDV T , Here D = diag(d 1 , . . . , d p ) contains the sample eigenvalues andV = [v 1 , . . . ,v p ] the corresponding sample eigenvectors. Following the earlier notation, we construct the sample component scores asŝ T j =v T j X, and the sample standardized scores asẑ We further assume that the population eigenvalues follow the spiked eigenvalue model introduced by Johnstone (2001), where the first m population eigenvalues are substantially larger than the remaining non-spiked eigenvalues.

Brief summary of earlier results
The question of consistency is central in statistics, as the sample estimates should converge to the population parameters when the sample size increases. Anderson (1963) showed that the sample eigenvectors and -values,v and d, will consistently estimate the population eigenvectors and -values, v and λ, when p is fixed and n → ∞.
However, this is not true in the high-dimensional setting, where p > n. Starting with Paul (2007) and Johnstone and Lu (2009), it has been shown that the sample eigenvalues and -vectors are not asymptotically consistent when p, n → ∞ at a constant ratio p/n = γ > 0 and the population eigenvalues are fixed. Paul (2007) showed that the inner product between the sample and the population eigenvector converges, when λ j > 1 + √ γ, to A different asymptotic setting starts from the geometrical structure of the data in a high-dimensional space (Ahn et al., 2007;Hall et al., 2005). Jung and Marron (2009) introduced the high dimension, low sample size (HDLSS) setting, where n is fixed and the spiked eigenvalues grow with the dimension p, according to λ i = σ 2 i p α , i = 1, . . . , m. In this asymptotic setting, the consistency of PCA depends on α, as p → ∞. Eigenvectors are estimated consistently when α > 1, while the estimates are strongly inconsistent when α < 1. In the boundary case α = 1, a situation explored by Jung et al. (2012), the sample eigenvectors are neither consistent nor strongly inconsistent, but reach a limiting distribution depending on n. In the case m = 1, where there is a single spiked eigenvalue, we have the following: The main focus of the above-mentioned papers has been on the eigenvector inconsistency, and few results are concerned with principal component scores. On exception is Lee et al. (2010), which established the asymptotic limit of the inner product between the sample and the population scores. They extended the result to prediction and found a theoretical asymptotic shrinkage factor for predicted scores. This can be applied as a bias adjustment, which turns out to be useful in the context of genetic population stratification problems. Also Aoshima (2009, 2012) explore the consistency of scores as p → ∞ and n → ∞. Leek (2011) showed that the estimated right singular vectors, which corresponds to the PC scores, in a low-dimensional conditional factor model converge (fixed n and p → ∞) to a set of vectors which span the same column space as the true factors.
Further, Shen et al. (2013Shen et al. ( , 2012 investigated the ratio between the individual sample and the population scoresẑ i /z i , i = 1, . . . , n, instead of the inner product between the score vectors. Following the regime of Jung and Marron (2009), they showed that for α > 1 the ratio converges to a random variable independent of i. This implies that a two-dimensional plot of the sample scores is asymptotically only a scaled version of the population score plot. The visual information contained in the samples scores will therefore remain the same as in the population scores. In this paper, we investigate the same problem as Shen et al. (2012), but in the situation where α = 1. We first motivate why this is an interesting assumption to make, and then prove the asymptotic behavior of the sample scores under this assumption.

Data structure and eigenvalues
The initial results concerning the consistency of sample eigenvalues and -vectors were derived on the basis of random matrix theory (Bai and Silverman, 2010). This requires the ratio p/n to remain constant as p, n → ∞ and the population eigenvalues to be fixed. In contrast, the HDLSS regime of Jung and Marron (2009) considers situations where the population eigenvalues depend asymptotically on the dimension p, according to λ i ∼ p α , α > 0.
Which of these two settings is the more appropriate remains an open question. As Lee et al. (2010) conclude: "It may be argued that for real data where p/n is "large," we should follow the paradigm of [the HDLSS regime]. However, for any real study, it is unclear how to test whether p increases at a faster rate than λ i or vice versa, making the application of [the HDLSS regime] difficult in practice." Establishing a natural connection between the eigenvalue model and real data problems is not an easy task. We believe it is appropriate to assume the eigenvalues to scale linearly with the dimension p, corresponding to the case of Jung et al. (2012) where α = 1. We argue in favor of this by introducing an assumption on the eigenvector coefficients. Our aim is to translate the assumption about the eigenvalues into an assumption regarding the latent structure and the data-generating mechanism, as this is generally easier to relate to.
First, let the observations x i be generated by a Gaussian latent variable model with an additive, isotropic error and a single factor: where the scalar z i ∼ N (0, 1) and the noise vector ǫ i ∼ N (0, σ 2 I). Tipping and Bishop (1999) established the connection between the Gaussian latent variable model under isotropic noise and PCA. The population covariance matrix of x i is then given by As v is an eigenvector of vv T when normalized, it is also an eigenvector of Σ, such that the corresponding eigenvalue of vv T is given by the normalizing constant λ 1 (vv T ) = p j=1 v 2 j . Thus the largest eigenvalue of Σ is given by The relationship between the largest population eigenvalue of Σ and the dimension can therefore be determined by the structure of v. For instance will the following three structures set p j=1 v 2 j and the eigenvalue to scale linearly with p, as p → ∞: i) the values of the v j are fixed, and the number of non-zero effects scales with p ii) the number of non-zero effects is fixed, and some values of the v j scale linearly in p iii) a combination of i) and ii), where the combined rate is linear It can be difficult to find realistic examples which would fulfill the settings ii) or iii). However the first situation can be interpreted in terms of pervasiveness (Fan et al., 2011): In the field of high-dimensional approximated factor models, Fan et al. (2011) refer to Definition 1 as the pervasiveness assumption. Here the number of non-zero entries in v is a non-vanishing proportion of the dimension p, as p increases. This stands in contrast to a sparse signal, where the number of non-zero effects is fixed, such that the proportion converges to zero.
If we assume the vector v to be pervasive with fixed values, meaning that the latent factor z i will have a pervasive effect on the observed variable x i , we have the following: and v is assumed pervasive with fixed values, the largest population eigenvalue λ 1 of the population covariance matrix Σ fulfills the bound The remaining population eigenvalues are given by λ i = σ 2 for i = 2, . . . , p.
The result follows from the existence of two constants 0 < c 1 ≤ c 2 < ∞, depending on r p and the minimum and maximum of the non-zero square loadings v 2 j respectively, such that the following bound is satisfied If the observations are given by m components where the v k , k = 1, . . . , m are orthogonal and pervasive with fixed values and p j=1 v 2 1j ≥ · · · ≥ p j=1 v 2 mj , the covariance matrix of x i will have m eigenvalues, which scale linearly with the dimension λ i ∼ p, i = 1, . . . , m.
It is also possible to interpret the pervasiveness in terms of the covariance matrix We can group the non-zero v i together into blocks, where the dimension of the blocks depends on the proportion r. We illustrate the situation in Example 1, where all effects are equal and the population covariance matrix consists of separate clusters, where all variables within each cluster are equally correlated, while the different clusters are independent. Each cluster corresponds to an eigenvector proportional to v j = [0, . . . , 0, 1, . . . , 1, 0, . . . , 0] with ones for the variables within the cluster and zeros for others. If the cluster sizes are not fixed, but a proportion of the total number of variables, the eigenvectors will be pervasive.
Each of the three largest eigenvalues represents one cluster and the importance of the cluster is determined by the proportion r j , how many variables that are represented, and the degree of correlation within the cluster ρ j . More strongly correlated variables will exhibit a clearer signal, while the remaining eigenvalues represent the noise. The pervasive eigenvectors can therefore be interpreted as variable clusters, where the cluster size is a percentage of the total number of variables.

Realistic examples from genomics
We present two situations in genomics, an area with several types of high-dimensional data, where the biological processes suggest the pervasiveness assumption to be reasonable. One example is genetic markers such as SNPs, single-basepair polymorphic genetic loci, i.e. having at least two alleles with an associate allelic frequency in a population. The neutral theory of molecular evolution states that allele frequencies at most loci (SNPs) change due to two stochastic processes; mutation and random drift.
If the main variation in the data sample stems from differences between ethnic populations, random allelic drift is the main driver behind changes in the genetic markers. This will give many and randomly distributed differences and when new markers are included, we expect a certain proportion to be informative with respect to the ethnicity. This corresponds to our notion of pervasive effects, and if the effects are fixed, the corresponding eigenvalue will scale linearly with total number of included variables. The longer two populations have been separated, the larger degree of SNPs expressing differences we expect, as observed by Yamaguchi-Kabata et al. (2008) when comparing Europeans and Japanese to subgroups within the Japanese population.
Another example is microarray expression data, which quantify the amount of a gene product called mRNA, whose expression is necessary for making proteins. Cancer is a relevant diseases in this respect, as it can be considered to have a systemic effect on gene expression (Perou et al., 2000). We therefore expect to observe many more or less randomly distributed differentially expressed genes between groups of cancer patients and healthy individuals. The meta-analysis of Kondrakhin et al. (2008) showed that 5% of 24726 genes are differentially expressed between cases of breast cancer and controls. This situation also corresponds to our notion of pervasive effects.

Asymptotic results
In the following, we present two results regarding the consistency of PC scores of highdimensional data. The asymptotic framework follows the high-dimension low sample size regime for the case where α = 1 as considered by Jung et al. (2012), and we state the same general conditions for the distribution of the component scores and for the structure of the population eigenvalues.
Firstly, the assumption of independent and normally distributed z ij can be relaxed, as x i has zero mean and covariance matrix Σ, to the following distributional condition: Condition 1. The standardized principal component scores z i have finite fourth moments and are uncorrelated but possibly dependent fulfilling the ρ-mixing condition.
The ρ-mixing condition is satisfied if the maximal correlation coefficient approach zero, and F L K is the σ-field of events generated by the variables z i , K ≤ i ≤ L. Secondly, the structure of the non-spiked eigenvalues λ m+1 , . . . , λ p in the spiked covariance model can be generalized by the following condition: Condition 2. For the eigenvalues λ m+1 , . . . , λ p , it must hold that Condition 2 insures that the non-spiked eigenvalues do not decrease too fast and that the mean converges to τ 2 . The constant non-spiked eigenvalues λ m+1 = · · · = λ p = τ 2 in the spiked covariance model is the simplest situation which fulfils condition 2.
Finally, we assume the spiked eigenvalues to scale linearly with the dimension: Assumption 1 (Linearity). For the m spiked components, the eigenvalues depend on the dimension p according to λ 1 = σ 2 1 p, λ 2 = σ 2 2 p · · · λ m = σ 2 m p, where σ 2 1 ≥ · · · ≥ σ 2 m > 0 represent the signal strength. Theorem 1 determines the asymptotic limiting distributions of the ratio between the sample and the population principal component scores under the Conditions 1 and 2, and Assumption 1. The results depend on the stochastic behavior of the eigenvalues and the eigenvectors of an m × m matrix W =Z T 1:mZ 1:m , whereZ 1:m = [σ 1 z 1 , . . . , σ m z m ]. We denote the jth eigenvalue of W by φ j (W) and the jth eigenvector by v j (W).
Theorem 1. Under Conditions 1 and 2, and Assumption 1 for m ≥ 1, the ratio between the sample and the population principal component scores converges in distribution to the following limit, as p → ∞: where the ratio R j is distributed as and ε ij is distributed as Remark 1. If the standardized component scores are assumed to be iid normally distributed z ij ∼ N (0, 1), i = 1, . . . , n, j = 1, . . . , p, W will be an m × m Wishart distributed matrix W ∼ W m diag(σ 2 1 , . . . , σ 2 m ), n − 1 , and φ j (W) and v j (W) will be asymptotically (as n → ∞) independent and normally distributed (Jolliffe, 2002).
Proof of Theorem 1. Theorem 1 follows from Lemmas 1 and 2, which are given by the results of Jung et al. (2012).
The normalized sample principal component scores can be decomposed by the expres- k v k z ik . By using the expression in (2), the ratio between the sample PC scores and the population PC scores can be decomposed and we can insert the m spiked population eigenvalues to give: for i = 1, . . . , n and j = 1, . . . , m. When p → ∞, the scaled sample eigenvalue converges according to Lemma 1 to p −1 d i d → φ i (W)/n + τ 2 /n, while the term consisting of the non-spiked eigenvalues can be rewritten as the vector By Lemma 2, a version of the law of large numbers, we have for a fixed m that 1 p p k=m+1 λ k z ik z lk → 0 for l = i and 1 p p k=m+1 λ k z 2 ik → τ 2 . Therefore this vector converges to the unit vector e i multiplied by τ 2 at position i and zero everywhere else: Then, according to Jung et al. (2012), the results in Lemma 1, the ratio between the sample and population scores converges tô The expression σ 2 z T kZ 1:m v j (W) in the first term corresponds to the kth row of W, and due to the eigen-equation Wv j (W) = φ j (W)v j (W), this term can be rewritten as while the unit vector in the second term gives

This simplifies expression (3) to
By splitting the sum, we get the result and we have used the simple notations

Implications for the visualization of scores
In the application of PCA, the first few sample scores are used for visualization and in conventional classification and regression methods. Besides the in itself valuable ability to visualize high-dimensional data in a two-dimensional (or 3D) fashion, the score plot can be useful for comparing observations, detecting subgroups, and for identifying outliers and bad data quality. PCA is often viewed as the canonical first step in an applied high-dimensional analysis. However, when it is known that eigenvectors are inconsistently estimated, will a plot of the sample scores give valid information about the population scores? We use Theorem 1 to answer this question in two steps. Firstly, we show by simulation and by providing a supporting theoretical argument, that ε ij is considerably smaller than R j . Therefore we refer to ε ij as noise. Secondly, we highlight the fact that the R j are independent of i. As the R j express ratios, the relative positions of the sample scores will be more or less the same as for the population scores. To illustrate these two points, we take a detailed look at the situation with two components, m = 2, in Example 3.

Example 2.
If there is only one component, m = 1, W is a scalar such that the estimated eigenvalue is given by φ 1 (W) = σ 2 1 z T 1 z 1 and the eigenvector is constant v 1 (W) = 1. Therefore, as ǫ i1 is zero, the limiting distribution of the ratio between the normalized sample and population scores is according to Theorem 1 given as .
The distributions depend on the two eigenvectors v j (W) = v j1 v j2 and eigenvalues φ j (W), First, we illustrate within Example 3 that ε ij can be considered as noise compared to R j . This is done by simulating from the distributions in (4) and (5) when the normalized scores are assumed to be standard normally distributed z ij ∼ N (0, 1). Tables 1 and 2 display the simulated mean and the standard deviation from the distributions of R 1 , ε i1 and R 2 , ε i2 respectively for different signal strength ratios θ = σ 2 2 /σ 2 1 and sample sizes n. The tables are shown graphically in Figure 1. We observe from Tables 1 and 2 that the expectation of ε ij is zero, whereas the R j are expected to be one. When taking the variability into account, we see that for moderate sample sizes (from 80 and upwards), a noise level of two SD is around 15-20 % of R 1 , still quite small. For larger sample sizes, the error drops to 2-8 % of R 1 . The noise also decreases when the separation between the signals increases. Only when σ 2 1 and σ 2 2 are close to each other and n is small, could ε ij be comparable to R j . The distribution of the noise is also illustrated graphically in Figure 3 by the 90 % probability contour for each observation.
To better understand the structure of the noise, we explore in Result 2 the distribution of ε ij for large n, a different asymptotic setting then considered earlier.
Three parameters have an effect on the noise distribution; the sample size n, the signal strength σ 2 , and the number of components m. From Result 2, we see the role of these parameters for the case of m = 2. As the sample size increases, the asymptotic variance of ε ij will decrease, and if we use the standard deviation as a measure of magnitude of the noise, it will decrease as n −1/2 . This is observed in Figure 1b), which displays graphically the simulated standard deviation of ε ij (as circles) for increasing sample size n together with the asymptotic values (dashed lines). The fit of the simulated standard deviation to the asymptotic scaling of n −1/2 necessarily becomes better as n increases.
The impact of the two signal strengths σ 2 1 and σ 2 2 is, as shown by Result 2, in terms of the ratio θ = σ 2 2 /σ 2 1 . The standard deviation of ε ij scales with the relative signal strength as θ/(1 − θ) for ε i1 . This is seen in Figure 1a), which displays the simulated standard deviation (as circles) for increasing relative signal strength also together with the asymptotic standard deviation (dashed lines). If σ 2 1 and σ 2 2 are close, the variability of the noise increases sharply, which can be interpreted as an overlap or interaction between the signals, making them difficult to distinguish. As observed in Tables 1 and 2, when the strength of the first signal increases relative to the second, the noise decreases.
Because there are only two components, v 21 (W) and v 12 (W) will have the same absolute value, but with opposite signs; hence they have a perfect negative correlation and this is seen in Figure 3. The noise for the second component is larger due to the scaling of σ 1 , which is reflected by greater extent of the contours in the vertical direction. Also the effect of the score ratios, z i1 /z i2 , on the noise is seen as the contours are wider closer to the x-axis.
The second key observation is that R 1 and R 2 are independent of i and thereby common to all observations. As they are ratios, they express a common scaling for the scores, which can be seen graphically as a shift, outwards or inwards. The consequence is that the relative positions of sample scores, and thereby most of the visual information, will be consistent with the population scores. We observe this in Figure 2, which displays the sample and population first and second PC scores for one simulated sample. A radial shift, which is not exact due to the noise, is evident when comparing the sample and population scores. Also the fact that R 2 will be generally larger than R 1 can be observed in Figure 2. The score plot will therefore be an appropriate tool to explore the population features, even though eigenvectors and absolute score values are not correctly estimated, asymptotically.

Conclusion
The use of high-dimensional PCA suffers from the somewhat paradoxical situation that theoretically, the eigenvectors and -values are not correctly estimated, but the method is highly successful in practice. The results in this paper attempt to bridge this gap by showing that the relative positions, and thereby the visual content, in a PC score plot are more or less the same for the true and the estimated scores. The assumption is that the leading eigenvalues scale linearly with the dimension. This assumption is fulfilled if the variability is caused by an latent factor with pervasive effects on the variables. This situation is reasonable in genetic markers from different ethnic populations and in microarray expression data from cancer cases and controls.
Future work should consider the implication of these results, when using principal component scores in further analyzes. The same asymptotic framework can be considered for regression, clustering and classification. Especially the effect of the limiting distribution on regression coefficients and misclassification rates would be of interest.