Como as crianças conseguem reunir seus pais em uma projeção PCA de um conjunto de dados do GWAS?

Pegue 20 pontos aleatórios em um espaço de 10.000 dimensões com cada ID de coordenada de . Divida-os em 10 pares ("casais") e adicione a média de cada par ("um filho") ao conjunto de dados. Em seguida, faça o PCA nos 30 pontos resultantes e plote PC1 vs PC2.$\mathcal N(0,1)$

Uma coisa notável acontece: cada "família" forma um trio de pontos que estão todos juntos. É claro que toda criança está mais próxima de cada um de seus pais no espaço original de 10.000 dimensões, de modo que se poderia esperar que ele estivesse perto dos pais também no espaço da PCA. No entanto, no espaço da PCA, cada par de pais também está próximo, embora no espaço original sejam apenas pontos aleatórios!

Como as crianças conseguem reunir os pais na projeção do PCA?

$\quad\quad\quad\quad$

Pode-se preocupar que isso seja de alguma forma influenciado pelo fato de os filhos terem uma norma mais baixa que os pais. Isso não parece ser a questão: se eu produzir as crianças como onde e são pontos dos pais, então eles vão ter, em média, a mesma norma como os pais. Mas ainda observo qualitativamente o mesmo fenômeno no espaço da APC:$(x+y)/\sqrt{2}$$x$$y$

$\quad\quad\quad\quad$

Essa pergunta está usando um conjunto de dados de brinquedos, mas é motivada pelo que observei em um conjunto de dados do mundo real de um estudo de associação genômica (GWAS) em que as dimensões são polimorfismos de nucleotídeo único (SNP). Esse conjunto de dados continha trios mãe-pai-filho.

Código

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)

def generate_families(n = 10, p = 10000, divide_by = 2):
    X1 = np.random.randn(n,p)    # mothers
    X2 = np.random.randn(n,p)    # fathers
    X3 = (X1+X2)/divide_by       # children
    X = []
    for i in range(X1.shape[0]):
        X.extend((X1[i], X2[i], X3[i]))
    X = np.array(X)

    X = X - np.mean(X, axis=0)
    U,s,V = np.linalg.svd(X, full_matrices=False)
    X = U @ np.diag(s)
    return X

n = 10
plt.figure(figsize=(4,4))
X = generate_families(n, divide_by = 2)
for i in range(n):
    plt.scatter(X[i*3:(i+1)*3,0], X[i*3:(i+1)*3,1])
plt.tight_layout()
plt.savefig('families1.png')

plt.figure(figsize=(4,4))
X = generate_families(n, divide_by = np.sqrt(2))
for i in range(n):
    plt.scatter(X[i*3:(i+1)*3,0], X[i*3:(i+1)*3,1])
plt.tight_layout()
plt.savefig('families2.png')

Durante a discussão com @ttnphns nos comentários acima, percebi que o mesmo fenômeno pode ser observado com muito menos de 10 famílias. Três famílias ( n=3no meu trecho de código) aparecem aproximadamente nos cantos de um triângulo equilátero. De fato, basta considerar apenas duas famílias ( n=2): elas acabam separadas ao longo de PC1, com cada família projetada aproximadamente em um ponto.

O caso de duas famílias pode ser visualizado diretamente. Os quatro pontos originais no espaço de 10.000 dimensões são quase ortogonais e residem em um subespaço 4-dimensional. Então eles formam um 4-simplex. Após a centralização, eles formarão um tetraedro regular, que é uma forma em 3D. Aqui está como ele se parece:

Antes que os filhos sejam adicionados, o PC1 pode apontar para qualquer lugar; não há direção preferida. No entanto, após duas crianças serem posicionadas no centro de duas arestas opostas, o PC1 passará por elas! Esse arranjo de seis pontos foi descrito por @ttnphns como um "haltere":

essa nuvem, como um haltere, tenderá a puxar os PCs principais para que eles perfurem as regiões pesadas

Observe que as arestas opostas de um tetraedro regular são ortogonais entre si e também ortogonais à linha que liga seus centros. Isso significa que cada família será projetada em um único ponto no PC1.

$\sqrt{2}$