Como @amoeba mencionado nos comentários, o PCA examinará apenas um conjunto de dados e mostrará os principais padrões (lineares) de variação nessas variáveis, as correlações ou covariâncias entre essas variáveis e os relacionamentos entre as amostras (as linhas ) no seu conjunto de dados.

O que normalmente se faz com um conjunto de dados de espécies e um conjunto de possíveis variáveis explicativas é ajustar uma ordenação restrita. No PCA, os principais componentes, os eixos no biplot do PCA, são derivados como combinações lineares ideais de todas as variáveis. Se você executou isso em um conjunto de dados de química do solo com variáveis pH, , TotalCarbon, talvez descubra que o primeiro componente foi $\mathrm{Ca^{2+}}$

0.5 \times p H + 1.4 \times C a^{2 +} + 0.1 \times T o t a l C a r b o n

$0.5 \times \mathrm{pH} + 1.4 \times \mathrm{Ca^{2+}} + 0.1 \times \mathrm{TotalCarbon}$

e o segundo componente

2.7 \times p H + 0.3 \times C a^{2 +} - 5.6 \times T o t a l C a r b o n

$2.7 \times \mathrm{pH} + 0.3 \times \mathrm{Ca^{2+}} - 5.6 \times \mathrm{TotalCarbon}$

Esses componentes são livremente selecionáveis a partir das variáveis medidas, e as que são escolhidas são aquelas que explicam sequencialmente a maior quantidade de variação no conjunto de dados e que cada combinação linear é ortogonal (não correlacionada com) as outras.

Em uma ordenação restrita, temos dois conjuntos de dados, mas não temos liberdade para selecionar as combinações lineares do primeiro conjunto de dados (os dados químicos do solo acima) que desejamos. Em vez disso, temos que selecionar combinações lineares das variáveis no segundo conjunto de dados que melhor explicam a variação no primeiro. Além disso, no caso do PCA, o único conjunto de dados é a matriz de resposta e não há preditores (você pode pensar na resposta como se ela própria estivesse se prevendo). No caso restrito, temos um conjunto de dados de resposta que desejamos explicar com um conjunto de variáveis explicativas.

Embora você não tenha explicado quais variáveis são a resposta, normalmente se deseja explicar a variação na abundância ou composição dessas espécies (ou seja, as respostas) usando as variáveis explicativas ambientais.

A versão restrita do PCA é uma coisa chamada Análise de Redundância (RDA) nos círculos ecológicos. Isso pressupõe um modelo de resposta linear subjacente para as espécies, que não é apropriado ou apenas apropriado se você tiver gradientes curtos ao longo dos quais as espécies respondem.

Uma alternativa ao PCA é uma coisa chamada análise de correspondência (CA). Isso é irrestrito, mas possui um modelo de resposta unimodal subjacente, que é um pouco mais realista em termos de como as espécies respondem ao longo de gradientes mais longos. Observe também que a CA modela abundâncias ou composição relativa , o PCA modela as abundâncias brutas.

Existe uma versão restrita do CA, conhecida como análise de correspondência restrita ou canônica (CCA) - que não deve ser confundida com um modelo estatístico mais formal conhecido como análise de correlação canônica.

Tanto na RDA quanto na CCA, o objetivo é modelar a variação na abundância ou composição de espécies como uma série de combinações lineares das variáveis explicativas.

Pela descrição, parece que sua esposa quer explicar a variação na composição (ou abundância) das espécies de milípedes em termos das outras variáveis medidas.

Algumas palavras de aviso; RDA e CCA são apenas regressões multivariadas; O CCA é apenas uma regressão multivariada ponderada. Tudo o que você aprendeu sobre regressão se aplica, e existem outras dicas:

à medida que você aumenta o número de variáveis explicativas, as restrições tornam-se cada vez menores e você não está realmente extraindo componentes / eixos que explicam a composição da espécie de maneira ideal e
com o CCA, à medida que aumenta o número de fatores explicativos, você corre o risco de induzir um artefato de curva na configuração de pontos no gráfico do CCA.
a teoria subjacente à RDA e à ACC é menos desenvolvida do que os métodos estatísticos mais formais. Só podemos escolher razoavelmente quais variáveis explicativas continuar usando a seleção por etapas (o que não é ideal por todos os motivos pelos quais não gostamos como um método de seleção em regressão) e precisamos usar testes de permutação para fazer isso.

então meu conselho é o mesmo que com a regressão; pense com antecedência quais são suas hipóteses e inclua variáveis que refletem essas hipóteses. Não basta jogar todas as variáveis explicativas na mistura.

Exemplo

Ordenação sem restrições

PCA

Vou mostrar um exemplo comparando PCA, CA e CCA usando o pacote vegan para R, que ajudo a manter e que foi projetado para atender a esses tipos de métodos de ordenação:

library("vegan")                        # load the package
data(varespec)                          # load example data

## PCA
pcfit <- rda(varespec)
## could add `scale = TRUE` if variables in different units
pcfit

> pcfit
Call: rda(X = varespec)

              Inertia Rank
Total            1826     
Unconstrained    1826   23
Inertia is variance 

Eigenvalues for unconstrained axes:
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
983.0 464.3 132.3  73.9  48.4  37.0  25.7  19.7 
(Showed only 8 of all 23 unconstrained eigenvalues)

vegan não padroniza a inércia, ao contrário de Canoco, então a variação total é 1826 e os valores próprios estão nessas mesmas unidades e somam 1826

> cumsum(eigenvals(pcfit))
      PC1       PC2       PC3       PC4       PC5       PC6       PC7       PC8 
 982.9788 1447.2829 1579.5334 1653.4670 1701.8853 1738.8947 1764.6209 1784.3265 
      PC9      PC10      PC11      PC12      PC13      PC14      PC15      PC16 
1796.6007 1807.0361 1816.3869 1819.1853 1821.5128 1822.9045 1824.1103 1824.9250 
     PC17      PC18      PC19      PC20      PC21      PC22      PC23 
1825.2563 1825.4429 1825.5495 1825.6131 1825.6383 1825.6548 1825.6594

Também vemos que o primeiro valor próprio é cerca da metade da variação e, com os dois primeiros eixos, explicamos ~ 80% da variação total

> head(cumsum(eigenvals(pcfit)) / pcfit$tot.chi)
      PC1       PC2       PC3       PC4       PC5       PC6 
0.5384240 0.7927453 0.8651851 0.9056821 0.9322031 0.9524749

Um biplot pode ser obtido a partir das pontuações das amostras e espécies nos dois primeiros componentes principais

> plot(pcfit)

insira a descrição da imagem aqui

Há duas questões aqui

A ordenação é essencialmente dominada por três espécies - essas espécies estão mais distantes da origem -, pois são os taxa mais abundantes no conjunto de dados
Existe um forte arco de curva na ordenação, sugestivo de um gradiente único longo ou dominante que foi dividido nos dois principais componentes principais para manter as propriedades métricas da ordenação.

CA

Uma AC pode ajudar com ambos os pontos, pois lida melhor com o gradiente longo devido ao modelo de resposta unimodal e modela a composição relativa das espécies, não a abundância bruta.

O código vegan / R para fazer isso é semelhante ao código PCA usado acima

cafit <- cca(varespec)
cafit

> cafit <- cca(varespec)
> cafit
Call: cca(X = varespec)

              Inertia Rank
Total           2.083     
Unconstrained   2.083   23
Inertia is mean squared contingency coefficient 

Eigenvalues for unconstrained axes:
   CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
0.5249 0.3568 0.2344 0.1955 0.1776 0.1216 0.1155 0.0889 
(Showed only 8 of all 23 unconstrained eigenvalues)

Aqui explicamos cerca de 40% da variação entre sites em sua composição relativa

> head(cumsum(eigenvals(cafit)) / cafit$tot.chi)
      CA1       CA2       CA3       CA4       CA5       CA6 
0.2519837 0.4232578 0.5357951 0.6296236 0.7148866 0.7732393

A parcela conjunta das pontuações de espécies e locais é agora menos dominada por algumas espécies

> plot(cafit)

insira a descrição da imagem aqui

Qual PCA ou CA você escolhe deve ser determinado pelas perguntas que deseja fazer dos dados. Geralmente, com dados de espécies, estamos mais interessados na diferença no conjunto de espécies, de modo que a CA é uma escolha popular. Se tivermos um conjunto de dados de variáveis ambientais, como a química da água ou do solo, não esperamos que eles respondam de maneira unimodal ao longo de gradientes, para que a CA seja inadequada e o PCA (de uma matriz de correlação, usando scale = TRUEa rda()chamada) seja mais apropriado.

Ordenação restrita; CCA

Agora, se tivermos o segundo conjunto de dados que desejamos usar para explicar padrões no primeiro conjunto de dados de espécies, devemos usar uma ordenação restrita. Geralmente, a escolha aqui é CCA, mas o RDA é uma alternativa, assim como o RDA após a transformação dos dados para permitir que ele lide melhor com os dados das espécies.

data(varechem)                          # load explanatory example data

Reutilizamos a cca()função, mas fornecemos dois quadros de dados ( Xpara espécies e Ypara variáveis explicativas / preditivas) ou uma fórmula de modelo listando a forma do modelo que desejamos ajustar.

Para incluir todas as variáveis que poderíamos usar varechem ~ ., data = varechemcomo fórmula para incluir todas as variáveis - mas como eu disse acima, essa não é uma boa ideia em geral

ccafit <- cca(varespec ~ ., data = varechem)

> ccafit
Call: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn +
Zn + Mo + Baresoil + Humdepth + pH, data = varechem)

              Inertia Proportion Rank
Total          2.0832     1.0000     
Constrained    1.4415     0.6920   14
Unconstrained  0.6417     0.3080    9
Inertia is mean squared contingency coefficient 

Eigenvalues for constrained axes:
  CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10  CCA11 
0.4389 0.2918 0.1628 0.1421 0.1180 0.0890 0.0703 0.0584 0.0311 0.0133 0.0084 
 CCA12  CCA13  CCA14 
0.0065 0.0062 0.0047 

Eigenvalues for unconstrained axes:
    CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8     CA9 
0.19776 0.14193 0.10117 0.07079 0.05330 0.03330 0.01887 0.01510 0.00949

O trio da ordenação acima é produzido usando o plot()método

> plot(ccafit)

insira a descrição da imagem aqui

Obviamente, agora a tarefa é descobrir qual dessas variáveis é realmente importante. Observe também que explicamos cerca de 2/3 da variação de espécies usando apenas 13 variáveis. Um dos problemas de usar todas as variáveis nessa ordenação é que criamos uma configuração em arco nas pontuações de amostras e espécies, que é puramente um artefato do uso de muitas variáveis correlacionadas.

Se você quiser saber mais sobre isso, consulte a documentação vegana ou um bom livro sobre análise de dados ecológicos multivariados.

Relação com regressão

É mais simples ilustrar o link com o RDA, mas o CCA é o mesmo, exceto que tudo envolve somas marginais de tabela e linha bidirecional como pesos.

No fundo, o RDA é equivalente à aplicação do PCA a uma matriz de valores ajustados a partir de uma regressão linear múltipla ajustada aos valores de cada espécie (resposta) (abundância), por exemplo, com preditores dados pela matriz de variáveis explicativas.

Em R, podemos fazer isso como

## centre the responses
spp <- scale(data.matrix(varespec), center = TRUE, scale = FALSE)
## ...and the predictors
env <- as.data.frame(scale(varechem, center = TRUE, scale = FALSE))

## fit a linear model to each column (species) in spp.
## Suppress intercept as we've centred everything
fit <- lm(spp ~ . - 1, data = env)

## Collect fitted values for each species and do a PCA of that
## matrix
pclmfit <- prcomp(fitted(fit))

Os valores próprios para essas duas abordagens são iguais:

> (eig1 <- unclass(unname(eigenvals(pclmfit)[1:14])))
 [1] 820.1042107 399.2847431 102.5616781  47.6316940  26.8382218  24.0480875
 [7]  19.0643756  10.1669954   4.4287860   2.2720357   1.5353257   0.9255277
[13]   0.7155102   0.3118612
> (eig2 <- unclass(unname(eigenvals(rdafit, constrained = TRUE))))
 [1] 820.1042107 399.2847431 102.5616781  47.6316940  26.8382218  24.0480875
 [7]  19.0643756  10.1669954   4.4287860   2.2720357   1.5353257   0.9255277
[13]   0.7155102   0.3118612
> all.equal(eig1, eig2)
[1] TRUE

Por alguma razão, não consigo obter as pontuações do eixo (cargas) correspondentes, mas, invariavelmente, elas são dimensionadas (ou não), por isso preciso examinar exatamente como elas estão sendo feitas aqui.

Nós não fazemos o RDA via rda()como mostrei no lm()etc, mas usamos uma decomposição QR para a parte do modelo linear e, em seguida, SVD para a parte do PCA. Mas os passos essenciais são os mesmos.

Gavin Simpson
fonte

+1 e procurando continuação! Vários comentários: (1) no seu exemplo PC1 não é ortogonal ao PC2; você pode, por exemplo, mudar de para para corrigi-lo. (2) Provavelmente faria sentido editar o título do OP para refletir o conteúdo da sua resposta; a edição atual do título é minha, mas eu não tinha ideia do que estava falando do OP. (3) A "resposta" é geralmente univariada ou multivariada? Parece o último, mas como é multivariada, por exemplo, a abundância de milípedes? Abundância de várias espécies? (4) Como esses métodos diferem da regressão? Você pode incluir alguns indicadores matemáticos?

+ 1.3

$+1.3$

- 5.6

$-5.6$

Ameba

Obrigado pela sugestão e acompanhamento - não me ocorreu tornar ortogonais os exemplos de combinações lineares, mas eu os atualizei. Re 2), fiz uma presunção, mas dado que há c. 12.000 espécies de milípedes, suspeito que a resposta aqui seja a observação da abundância de espécies em cada um dos locais de amostragem. Nesse sentido, o RDA ou o CCA modelariam uma matriz de resposta multivariada de dimensão . Vou tentar lidar com 4 depois, depois de colocar as crianças na cama.

m

$m$

n

$n$

n \times m

$n \times m$

Gavin Simpson

@amoeba Desculpe pelo atraso, mas adicionei uma seção à minha resposta para tentar mostrar o link com a regressão e como o RDA pode ser visto como um PCA de valores ajustados de uma série de regressões lineares, uma por variável de resposta.

Gavin Simpson

@amoeba Fazemos SVD de (valores ajustados) e não dos coeficientes , pelo menos é o que fornece: . Portanto, a RDA é freqüentemente chamada de regressão de classificação reduzida.

X β

$\mathbf{X}\beta$

β

$\beta$ fitted()

X β

$\mathbf{X}\beta$

Gavin Simpson

A origem da RDA deve-se a Rao (1964), que é um artigo estatístico, portanto deve ser apropriado.

Gavin Simpson

Quais critérios usar para separar variáveis em variáveis explicativas e respostas para métodos de ordenação em ecologia?

Respostas:

Exemplo

Ordenação sem restrições

PCA

CA

Ordenação restrita; CCA

Relação com regressão

Quais critérios usar para separar variáveis ​​em variáveis ​​explicativas e respostas para métodos de ordenação em ecologia?

Respostas:

Exemplo

Ordenação sem restrições

PCA

CA

Ordenação restrita; CCA

Relação com regressão

Quais critérios usar para separar variáveis em variáveis explicativas e respostas para métodos de ordenação em ecologia?