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

9

Eu tenho diferentes variáveis ​​que interagem dentro de uma população. Basicamente, tenho feito um inventário de milípedes e medido alguns outros valores do terreno, como:

  • As espécies e a quantidade de espécimes coletados
  • Os diferentes ambientes onde os animais estão
  • o pH
  • A porcentagem de material orgânico
  • a quantidade de P, K, Mg, Ca, Mn, Fe, Zn, Cu
  • Relação Ca + Mg / K

Basicamente, eu gostaria de usar o PCA para determinar quais variáveis ​​direcionam a variabilidade das amostras e tornam a floresta (ambientes) diferente; quais variáveis ​​devo usar para "variáveis" e quais para "indivíduos"?

Leonardo
fonte
Eu acho que você pode estar confuso sobre o PCA. Todas as variáveis ​​podem (é claro) ser apenas "variáveis". Você provavelmente fez várias medições em locais diferentes (ou em momentos diferentes); esses locais (ou horários) são seus "indivíduos", ou melhor, "amostras".
Ameba
Além disso, não posso deixar de perguntar: seu perfil diz que você é um fundador de startups; é uma startup trabalhando com milípedes? Uau!
Ameba
na verdade, @amoeba é minha esposa, trabalhando bem nisso, sou bom em cálculo, mas não me desenvolvo muito bem em estatística. E ela queria que eu perguntasse.
Leonardo
11
Até agora, isso é realmente uma questão estatística? Embora pareça confuso quanto à terminologia estatística, na medida em que é difícil decodificar, se isso for resolvido, o conselho é usar o julgamento científico.
Nick Cox
2
Isso pode muito bem ser uma questão estatística, apenas em um contexto diferente e com termos diferentes dos que são mais comuns em estatística. Eu acredito que você está perguntando sobre os métodos de ordenação da ecologia. Este site pode ser útil para você. Relativamente poucos membros ativos aqui são especialistas nisso, mas o @GavinSimpson pode ajudá-lo se conseguirmos chamar a atenção dele.
gung - Restabelece Monica

Respostas:

18

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 foiCa2+

0.5×pH+1.4×Ca2++0.1×TotalCarbon

e o segundo componente

2.7×pH+0.3×Ca2+5.6×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

  1. 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
  2. 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
4
+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? - 5,6+1.35.6
Ameba
2
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. n n × mmnn×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
11
@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 βXββfitted()Xβ
Gavin Simpson
11
A origem da RDA deve-se a Rao (1964), que é um artigo estatístico, portanto deve ser apropriado.
Gavin Simpson