An Eigenvalue test for spatial principal component analysis

Background The spatial Principal Component Analysis (sPCA, Jombart (Heredity 101:92-103, 2008) is designed to investigate non-random spatial distributions of genetic variation. Unfortunately, the associated tests used for assessing the existence of spatial patterns (global and local test; (Heredity 101:92-103, 2008) lack statistical power and may fail to reveal existing spatial patterns. Here, we present a non-parametric test for the significance of specific patterns recovered by sPCA. Results We compared the performance of this new test to the original global and local tests using datasets simulated under classical population genetic models. Results show that our test outperforms the original global and local tests, exhibiting improved statistical power while retaining similar, and reliable type I errors. Moreover, by allowing to test various sets of axes, it can be used to guide the selection of retained sPCA components. Conclusions As such, our test represents a valuable complement to the original analysis, and should prove useful for the investigation of spatial genetic patterns. Electronic supplementary material The online version of this article (10.1186/s12859-017-1988-y) contains supplementary material, which is available to authorized users.


INTRODUCTION
The principal component analysis (PCA ;Pearson 1901;Hotelling 1933) is one of the most common multivariate approaches in population genetics (Jombart et al 2009). Although PCA is not explicitly accounting for spatial information, it has often been used for investigating spatial genetic patterns (Novembre and Stephens 2008). As a complement to PCA, the spatial principal component analysis (sPCA; Jombart et al. 2008) has been introduced to explicitly include spatial information in the analysis of genetic variation, and gain more power for investigating spatial genetic structures.
sPCA finds synthetic variables, the principal components (PCs), which maximise both the genetic variance and the spatial autocorrelation as measured by Moran's I (Moran 1950).
As such, PCs can reveal two types of patterns: 'global' structures, which correspond to positive autocorrelation typically observed in the presence of patches or clines, and 'local' structures, which correspond to negative autocorrelation, whereby neighboring individuals are more genetically distinct than expected at random (for a more detailed explanation on the meaning of global and local structures see Jombart et al.. 2008). The global and local tests have been developed for detecting the presence of global and local patterns, respectively . Unfortunately, while these tests have robust type I error, they also typically lack power, and can therefore fail to identify existing spatial genetic patterns . Moreover, they can only be used to diagnose the presence or absence of spatial patterns, and are unable to test the significance of specific structures revealed by sPCA axes.
In this paper, we introduce an alternative statistical test which addresses these issues.
This approach relies on computing the cumulative sum of a defined set of sPCA

Test statistic
As in most multivariate analyses of genetic markers, our approach analyses a table of centred allele frequencies (i.e. set to a mean frequency of zero), in which rows represent individuals or populations, and columns correspond to alleles of various loci Jombart et al 2009;Jombart et al 2010). We note X the resulting matrix, and n the number of individuals analysed. In addition, the sPCA introduces spatial data in the form of a n by n matrix of spatial weights L, in which the i th row contains weights reflecting the spatial proximity of all individuals to individual i. The PCs of sPCA are then found by the eigen-analysis of the symmetric matrix ): We note λ the corresponding non-zero eigenvalues. We differentiate the r positive eigenvalues λ + , corresponding to global structures, and the 's' negative eigenvalues λ -, corresponding to local structures, so that λ = {λ + ,λ -}. Without loss of generality, we assume both sets of eigenvalues are ordered by decreasing absolute value, so that λ1 + > λ2 + > … > λr + and |λ1 -| > |λ2 -| > … > |λs -|. Simply put, each eigenvalue quantifies the magnitude of the spatial genetic patterns in the corresponding PC: larger absolute values indicate stronger global (respectively local) structures. We note V + = {v1 + , …, vr + } and V -= {v1 -, …, vs -} the sets of corresponding PCs.The most natural choice of test statistic to assess whether a given PC contains significant structure would seem to be the corresponding eigenvalue. This would, however, not account for the dependence on previous PCs: vj + (respectively vj -) can only be significant if all previous PCs {v1 + , …, vj-1 + } are also significant. To account for this, we define the test statistic for vj + as:

Permutation procedure
fi + and fibecome larger in the presence of strong global or local structures in the first i th global / local PCs. Therefore, they can be used as test statistics against the null hypotheses of absence of global or local structures in these PCs. The expected distribution of fi + and fiin the absence of spatial structure is not known analytically.
Fortunately, it can be approximated using a Monte-Carlo procedure, in which at each permutation individual genotypes are shuffled to be assigned to a different pair of coordinates than in the observed original dataset and fi + and fiare computed. Note that the original values of the test statistic are also included in these distributions, as the initial spatial configuration is by definition a possible random outcome. The p-values are then computed as the relative frequencies of permuted statistics equal to or greater than the initial value of fi + or fi -.
To guide the selection of global and local PCs to retain, the simulated values of each eigenvalue (from most positive to most negative), which make up the fi + and fistatistics, are also recorded during the permutation procedure. In this way, if global or local structures are detected to be significant, an observed p-value for each observed eigenvalue can be estimated by comparison with its simulated eigenvalue distribution.
Note that the number of eigenvalues produced by an sPCA does not change between the observed and permutated datasets, so each observed eigenvalue can be compared with the distribution of the corresponding simulated one. This testing procedure can be used with increasing numbers of retained axes. Because each test is conditional on the previous tests, incremental Bonferroni correction is used to avoid the inflation of type I error, so that the significance level for the i th PC will be α / i, where α is the target type I error. Hence, the correction implies that if the most positive (or negative) eigenvalue is significant in regards with the chosen p-value threshold, the second eigenvalue is tested for a p-value threshold that is the half of the previous and so on. The entire testing procedure is implemented in the function spca_randtest in the package adegenet (Jombart 2008; Jombart and Ahmed 2011) for R (R Core Team 2017). A flow chart of the test procedure is shown in Figure 1.

Simulation study
To assess the performance of our test, we simulated genetic data under three migration models: island (IS) and stepping stone (SS), using the software and homogeneous rates (0.005, 0.01, and 0.1). We performed 100 independent runs for each of the three migration rates, for a total of 300 simulated dataset per migration model.
To quantify type I error rates for the spca_randtest, global and local tests, we extracted 100 random coordinates from 10 square 2D grids, using the function spsample from the spdep package (Bivand et al. 2013). In order to evaluate the rate of false negatives for global patterns, we manually generated 10 sets of 100 pairs of coordinates simulating gradients and/or patches from 2D grids. An example of simulated global patterns is presented in Figure 2. To test for the rate of false negatives for local patterns, we perform

Statistical power of the spca_randtest
We compared the performances of the spca_randtest with the global and local tests in three settings: in the absence of spatial structure, and in the presence of global, and local structures. The results obtained in the absence of spatial structure show that all tests have reliable type I errors (Table 1 and 2). The spca_randtest exhibited consistently better performances for detecting existing structures in the data than both global and local tests (Table 1 and 2). Although our simulated local spatial patterns turned out more difficult to detect than global patterns, the spca_randtest is twice to five times more effective than the local test (Table 1 and 2). Generally, the underlying migration model, the migration rate and the number of loci affect the ability of all tests to detect non-random spatial patterns.
Both spca_randtest and global and local tests have in fact a lower sensitivity in presence of island migratory schemes, while results for stepping stone and isolation by distance models are more satisfying (Table 1 and 2). Increasing migration rates lead to a higher rates of false negatives for all tests, which can be overcome using more loci (Table 1 and 2).
Significant eigenvalues are assessed using a hierarchical Bonferroni correction which accounts for non-independence of eigenvalues and multiple testing ( Figure 2). Strong patterns (e.g. IBD) tend to produce a higher number of significant components than weak patterns (e.g. island models with high migration rates), which are otherwise captured by fewer to no components.

Application to real data
We have run the sPCA to compare the new spca_randtest and previous tests to a real

DISCUSSION
We introduced a new statistical test associated to the sPCA to evaluate the statistical significance of global and local spatial patterns. Using simulated data, we show that this new approach outperforms previously implemented tests, having greater statistical power (lower type II errors) whilst retaining consistent type I errors. Our simulations also suggest that demographic settings and migratory models can substantially impact the ability to detect spatial patterns. Indeed, high migration rates, non-hierarchical migration models, such as island model, and low amount of loci can hamper or worsen the performance of the test, preventing the detection of actual spatial patterns. In lack of previous information on the demographic history and/or the movement ecology of the population under study, it is certainly useful to exploit all the available genetic information. In this regards, our simulations show how an increased number of loci does improve the ability of the test to provide meaningful results.
The impact of specific factors such as the effective population size or the number of individuals sampled per population remain to be investigated. A more extensive simulation study, possibly comparing different non-model based methods such as sPCA, would clarify the extent of the spatial information that can be obtained with such methods without comparing explicit evolutionary hypotheses. In fact, the sPCA and the associated spca_randtest cannot distinguish between explicit migration models. However, the possibility to detect which eigenvalues contain the spatial information provides the user with further information to interpret the biological meaning of the spatial structure, by focusing on few meaningful dimensions.
Our data application seems to confirm that the spca_randtest is more effective than global