The normalization is the whole method. Centring alone leaves each
marker weighted by its own allele frequency, so common SNPs dominate
the leading components for no genetic reason. Patterson, Price &
Reich divide instead by the drift standard deviation the binomial
model predicts, p_j = mu_j / 2 and
M_ij = (g_ij - mu_j) / sqrt(p_j (1 - p_j)), which puts every
marker on the same footing, and then take the eigenvectors of
\(X = M M' / m\) over INDIVIDUALS – an n by n problem however
many markers there are.
Value
List with eigenvalues, pcs,
variance_explained, estimate, tw_statistic,
n_eff, n_dropped, n, m.
Details
Significance of the leading eigenvalue is judged against the
Tracy-Widom law. With \(l1' = n lambda_1 / sum(lambda)\) and
\(n' = (n + 1) (sum l)^2 / ((n - 1) sum l^2 - (sum l)^2)\),
mu = (sqrt(n - 1) + sqrt(n'))^2 / n',
\(sigma = ((sqrt(n-1) + sqrt(n')) / n')
(1/sqrt(n-1) + 1/sqrt(n'))^(1/3)\), the statistic
\((l1' - mu) / sigma\) follows TW1 under the null of no
structure. \(n'\) is the EFFECTIVE number of markers and is not
m: linkage makes markers less independent than they look, and
using m overstates significance.
Monomorphic markers carry no information and would divide by zero; they are dropped and counted rather than left to produce Inf.
References
Patterson, N., Price, A. L. & Reich, D. (2006). Population structure and eigenanalysis. PLoS Genetics, 2(12), e190. doi:10.1371/journal.pgen.0020190
Examples
M <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2)
Pcasnps(M)
#> $eigenvalues
#> [1] 2.666667 0.000000
#>
#> $pcs
#> [,1] [,2]
#> [1,] 0.7071068 0.7071068
#> [2,] -0.7071068 0.7071068
#>
#> $variance_explained
#> [1] 1 0
#>
#> $estimate
#> [1] 2.666667
#>
#> $tw_statistic
#> [1] NaN
#>
#> $n_eff
#> [1] NaN
#>
#> $n_dropped
#> [1] 2
#>
#> $n
#> [1] 2
#>
#> $m
#> [1] 1
#>
#> $method
#> [1] "EIGENSTRAT genotype PCA (Patterson-Price-Reich 2006)"
#>
#> attr(,"class")
#> [1] "morie_rich_result" "list"