Ordenação
Paulo Paiva
14/05/2020
Análises Multivariadas de Ordenação
As análises multivariadas recebem este nome pois envolvem duas ou mais variáveis dependentees. Uma das análises mais comuns, especialmente em estudos ecológicos e de morfometria, é a análise de componentes principais.
Análise de Componentes Principais
Para a análise de componentes principais vamos usar outro dataset, iris, que envolve medidas de largura e comprimento de sépalas e pétalas de flores de três espécies do gênero Iris.
data(iris)
?iris
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaExistem várias funções no R para análise de componentes principais (ACP, ou PCA em inglês), duas delas, prcomp ou princomp, ambos do pacote stats são carregadas automaticamente pelo R. Existem outros pacotes que também fazem a mesma análise, mas isso é assunto para outro momento.
iris.pca <- prcomp(iris[,1:4])O comando prcomp foi aplicado apenas as colunas 1 a 4 e em todas as suas linhas - lembre da notação entre os colechetes após o nome do data frame [,1:4]. Neste caso usamos apenas o comando prcomp com o seu default - para saber qual é o default do comando, basta ver o help (?prcomp). Neste caso o parâmetro scale=FALSE signfica que, para realizar a análise foi utilizada uma matriz de covariância. Caso parâmetro fosse scale=TRUE, seria efetuada uma análise com uma matriz de correlação. Na primeira, as variáveis que tem uma variação maior tem um peso maior na análise, enquanto na segunda, as variáveis tem o mesmo peso, todas são padronizadas para ter média igual a 0 e variância igual a 1.
Veja com calma em um livro sobre análise multivariada qual o tipo de matriz para o seu caso, mas lembre-se nunca use matriz de covariância (scale=FALSE) quando as variáveis foram medidas em escalas diferentes. Aqui no caso dos dados de iris as medidas são todas expressas na mesma escala (centímetros).
Vamos ver agora uma síntes dos resultados com os comandos:
summary(iris.pca)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 2.0563 0.49262 0.2797 0.15439
## Proportion of Variance 0.9246 0.05307 0.0171 0.00521
## Cumulative Proportion 0.9246 0.97769 0.9948 1.00000Na janela temos a contribuição de cada componente (PC1 a PC4), a proporção da variância total explicada por cada um e, abaixo, a proporção cumulativa (esta deve somar um total de 1.000 ou 100% da variação, pois é a variância total de todos os componentes).
A contribuição pode ser vista graficamente com o comando plot
plot(iris.pca)
Já o comando print fornece outros resultados:
print(iris.pca)## Standard deviations (1, .., p=4):
## [1] 2.0562689 0.4926162 0.2796596 0.1543862
##
## Rotation (n x k) = (4 x 4):
## PC1 PC2 PC3 PC4
## Sepal.Length 0.36138659 -0.65658877 0.58202985 0.3154872
## Sepal.Width -0.08452251 -0.73016143 -0.59791083 -0.3197231
## Petal.Length 0.85667061 0.17337266 -0.07623608 -0.4798390
## Petal.Width 0.35828920 0.07548102 -0.54583143 0.7536574Neste são fornecidos novamente os desvios padrões (a contribuiçõ de cada componente) e as coordenadas de cada uma das variáveis originais (FL, RW, CL, CW e BD) nos eixos PC1 a PC4.
A visualização destes resultados para os dois primeiros componentes principais (PC1 e PC2) é fornecida através do comando:
biplot(iris.pca)
Onde as variáveis são representadas por eixos (cujas extremidades são as coordenadas que vimos em print) e as observações como pontos dispersos na figura.
As posições de cada observação na figura também são armazenadas na memória dentro do objeto iris.pca. Quando chamamos o •summary* ou o print apenas parte dos resultados são mostrados. Para saber quais os resultados que fazem parte desta análise, temos que ver a estrutura do objeto iris.pca:
str(iris.pca)## List of 5
## $ sdev : num [1:4] 2.056 0.493 0.28 0.154
## $ rotation: num [1:4, 1:4] 0.3614 -0.0845 0.8567 0.3583 -0.6566 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
## $ center : Named num [1:4] 5.84 3.06 3.76 1.2
## ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## $ scale : logi FALSE
## $ x : num [1:150, 1:4] -2.68 -2.71 -2.89 -2.75 -2.73 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
## - attr(*, "class")= chr "prcomp"Note que existem vários objetos dentro de iris.pca, todos eles aparecem após o sinal de cifrão ($sdev, $rotation, $center, $scale, $x). Vamos ver um destes objetos e o que significa:
head(iris.pca$x)## PC1 PC2 PC3 PC4
## [1,] -2.684126 -0.3193972 0.02791483 0.002262437
## [2,] -2.714142 0.1770012 0.21046427 0.099026550
## [3,] -2.888991 0.1449494 -0.01790026 0.019968390
## [4,] -2.745343 0.3182990 -0.03155937 -0.075575817
## [5,] -2.728717 -0.3267545 -0.09007924 -0.061258593
## [6,] -2.280860 -0.7413304 -0.16867766 -0.024200858Este objeto é um data frame que tem justamente as coordenadas de cada individuo para cada componente principal, ou seja, no nosso biplot a posição dos pontos são as coordenadas de cada ponto (observação) em PC1 e PC2.
Como o banco de dados que usamos se refere a três espécies diferentes, espécies estas que nós não utilizamos na análise, nós podemos, de forma análoga que fizemos no plot da regressão anterior, colorir as três espécies para procurar ver se estas apesentam padrões diversos. Primeiro vamos separar os escores dos dois primeiros componentes que estão no objeto iris.pca$x. Como os dois primerios componentes são as duas primeiras colunas, vamos criar um novo objeto da seguinte forma:
escores<-as.data.frame(iris.pca$x[,1:2])
head(escores)## PC1 PC2
## 1 -2.684126 -0.3193972
## 2 -2.714142 0.1770012
## 3 -2.888991 0.1449494
## 4 -2.745343 0.3182990
## 5 -2.728717 -0.3267545
## 6 -2.280860 -0.7413304Agora vamos plotar os dois colorindo as diferentes espécies com cores diferentes
plot(escores$PC2~escores$PC1, pch=(19), col=c("red","blue", "green") [as.numeric(iris$Species)])
col=c(“red”“,”blue“,”green“)[as.numeric(crabs$sp)])
Ou de uma forma mais bonitinha usando o ggplot2
escores$Species<-iris$Species
head(escores)## PC1 PC2 Species
## 1 -2.684126 -0.3193972 setosa
## 2 -2.714142 0.1770012 setosa
## 3 -2.888991 0.1449494 setosa
## 4 -2.745343 0.3182990 setosa
## 5 -2.728717 -0.3267545 setosa
## 6 -2.280860 -0.7413304 setosaNote que pegamos a variável(=coluna) Species que estava no data frame iris e “colamos” ela na matriz escores . Com isso, agora podemos plotar:
ggplot(escores, aes(PC1,PC2, colour=Species))+geom_point()Note que plotamos apenas os dois primeiros componentes, mas não necessariamente estes dois componentesos únicos a serem interpretados. Tudo depende da signficância dos componentes. Existe uma série de abordagens e técnicas para decidir quais componente são interpretáveis. Duas delas, são apresentada abaixo através de um função criada pelos autores do livro Boccard et al. 2010. Numerical Ecology with R:
evplot <- function(ev)
{
# Broken stick model (MacArthur 1957)
n <- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
# Plot eigenvalues and % of variation for each axis
op <- par(mfrow=c(2,1))
barplot(ev, main="Eigenvalues", col="bisque", las=2)
abline(h=mean(ev), col="red")
legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")
barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE,
main="% variation", col=c("bisque",2), las=2)
legend("topright", c("% eigenvalue", "Broken stick model"),
pch=15, col=c("bisque",2), bty="n")
par(op)
}Você pode salvar esta função no seu computador colando ela num script e salvando depois com o nome: evplot.R. Toda vez que quiser usar, basta chamar a função:
source("evplot.R")Aplicando a função aos nosso resultado acima para os valores do desvio padrão dos componentes principais (PCs), temos:
evplot(iris.pca$sdev)Na primeira figura temos o plot, novamente, dos componentes mostrando o valor médio dos desvios sugerindo que apenas o primeiro componente é interpretável, por ser maior que a média. O segundo plot é o do modelo Broken-Stick, em que a porcentagem de variação dos componentes é comparada com a probabilidade dos mesmos valores seerem frutos do acaso, representada aqui pelas barras vermelhas. Note que, novamente, o segundo componente é menor que o acaso, portanto apenas o primeiro é interpretável.
Análise Discriminante
Na Análise Discriminante as diversas variáveis são ordenadas a partir de um critério de classificação. Por exemplo, no caso acima em que efetuamos a ACP, as medidas das pétalas e sépalas das flores foram ordenadas baseadas apenas nas suas medidas. Por esta razão, estes tipos de análises multivariadas são também denominados de análises de ordenação. Só depois da análise é que nós plotamos as espécies com cores diferentes. Entretanto, uma outra estratégia seria avaliar se é possível discriminar (=distinguir) as diferentes espécies utilizando-se as mesmas medidas morfométricas usadas na ACP. Para isso temos que dizer antes de realizar a análise qual é a variável que vai ser usada na distinção, no caso, a variável Species . Vamos usar então o dataset iris novamente e a função lda do pacote MASS.
library(MASS)Primeiro vamos separar as variáveis numéricas, as medidas das pétalas e sépalas, fazendo uma nova matriz delas que chamaremos de medidas
medidas<-as.matrix(iris[,1:4])
head(medidas)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 5.1 3.5 1.4 0.2
## [2,] 4.9 3.0 1.4 0.2
## [3,] 4.7 3.2 1.3 0.2
## [4,] 4.6 3.1 1.5 0.2
## [5,] 5.0 3.6 1.4 0.2
## [6,] 5.4 3.9 1.7 0.4Agora vamos fazer a análise discrimante colocando do lado esquerdo a variável que vai distinguir (no caso a coluna de número 5 que é Species) e do lado direito as variáveis numéricas, que salvei como a matriz medidas
meu.discr<-lda(iris[,5]~medidas)E agora plotando as duas primeiras dimensões, ou seja as duas primeiras funções discriminates, semelhante ao que fizemos na ACP:
plot (meu.discr)
Ou se preferir mais bonitinho usando o ggplot
attach(iris)Mas para isso eu preciso, do valor de cada observação para as funções discriminantes geradas:
attach(iris)
preds<- predict(meu.discr)Plotando estas funções com as observações:
library(ggplot2)
funcoes<-as.data.frame(preds$x)
ggplot(funcoes, aes(LD1,LD2, colour=iris$Species))+geom_point()
E vendo os valores das funções e das médias de cada uma das espécies para cada uma das variáveis e as coordenadas nos dois primeiros eixos das medidas, mostrando qual a medida que mais influencia na separação:
meu.discr## Call:
## lda(iris[, 5] ~ medidas)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## medidasSepal.Length medidasSepal.Width medidasPetal.Length
## setosa 5.006 3.428 1.462
## versicolor 5.936 2.770 4.260
## virginica 6.588 2.974 5.552
## medidasPetal.Width
## setosa 0.246
## versicolor 1.326
## virginica 2.026
##
## Coefficients of linear discriminants:
## LD1 LD2
## medidasSepal.Length 0.8293776 0.02410215
## medidasSepal.Width 1.5344731 2.16452123
## medidasPetal.Length -2.2012117 -0.93192121
## medidasPetal.Width -2.8104603 2.83918785
##
## Proportion of trace:
## LD1 LD2
## 0.9912 0.0088Note que as variáveis Petal.Length e Petal.Width são negativas em LD1 com valores muito baixos (<-2) enquanto Sepal.Width é positiva. Este resultado indica que as pétalas são mais longas e largas no lado esquerdo da figura, caracterizando as espécies Iris virginica e Iris versicolor. Já as sépalas são mais largas no lado direito, caracterizando a espécie Iris setosa. Você pode conferir este padrão olhando as médias destas medidas para cada uma das espécies.
A figura pode ser alterada, mas não abordaremos este assunto agora.
Mas um dos resultados mais importantes da análise discriminante é a chamada Matriz de Classificação, que é uma matriz onde a classificação original das observações (ou seja, a qual espécie aquele indivíduo pertence) é comparada com o resultado da análise discriminante, ou seja, considerando a combinação de medidas, a qual espécie este indivíduo seria mais provavelmente alocado. Note pela tabela abaixo que quase sempre há uma coincidênica, mas que no caso de I. versicolor dois indivíduos foram previstos (predicted) como pertencentes a I.virginica. Por outro lado um exemplar de I.virginica foi alocado em I. versicolor. Esta alocação se deve a probabilidade posterior (aquela fruto da análise discriminante), indicar outra espécie do que a que ele realmente pertence.
Um resultados obtido quando se aplica a função discriminante é a classe a qual ele tem maior probabilidade de pertencer em função das medidas tomadas. Abaixo podemos comparar esta classificação que é obtida pela probabilidade posterior com a classe original (= espécie) a qual aquele indivíduo pertencia. Na comparação abaixo fica claro que na grande maioria dos casos os dados são concordantes, embora nos indivíduos de número 71, 84 e 134 as classificações divergem.
comparacao<-as.data.frame(cbind(as.character(iris[,5]), as.character(preds$class)))Você pode ver a probabilidade posterior de cada um destes casos “suspeitos” abaixo, notando que, por exemplo, o primeiro caso (71) foi alocado como I.virginica 74,7% de probabilidade embora fosse originalmente I.versicolor ao qual foi alocado com apenas 25,32% de probabilidade
suspeitos<-rbind(preds$posterior[71,], preds$posterior[84,],preds$posterior[134,])
rownames(suspeitos)<-c("71","84", "134")
suspeitos## setosa versicolor virginica
## 71 7.408118e-28 0.2532282 0.7467718
## 84 4.241952e-32 0.1433919 0.8566081
## 134 1.283891e-28 0.7293881 0.2706119tabelaDA <- table(iris[,5], preds$class)
tabelaDA##
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 1 49Este tipo de resultado normalmente é apresentado na forma de porcentagem onde o número dos indivíduos classificados corretamente em cada uma das espécies, ou seja, aqueles cuja classificação posterior coincide com o grupo (espécie) a qual o indivíduo originalmente pertence. Assim temos uma porcentagem de classificação correta de cada uma das espécies e uma porcentagem total de acertos onde o total dos indivíduos classificadoss como pertencentes a própria espécie (ou seja a digonal da matriz) é dividido pelo número total de indivíduos.
diag(prop.table(tabelaDA, 1))## setosa versicolor virginica
## 1.00 0.96 0.98sum(diag(prop.table(tabelaDA)))## [1] 0.98Neste caso tivemos uma porcentagem de acertos de 98%.
NMDS - Escalonamento não-métrico Multidimensional
Muitas vezes, na ordenação, o objetivo é apenas ordenar as observações baseadas apenas nas suas similaridades representadas pela proximidade na figura. Ou seja, observações mais similares são mais próximas entre si. Estas similaridades são ordenadas da maior para a menor e a proximidade não é diretamente proporcional à distância, por isso se trata de um método não métrico.
O escalonamento multidimensional não-métrico (NMDS), às vezes denominado apenas como MDS, é uma técnica muito utilizada para este tipo de plotagem, especialmente em estudos de impacto ambiental
Um pacote muito útil para diversas analises multivariadas, incluindo várias não abordadas aqui, é o vegan
library("vegan")Aqui, vamos usar os mesmos dados das medidas das sépalas e pétalas de iris com a função metaMDS:
iris.mds<-metaMDS(iris[,1:4], "euclidean")## Run 0 stress 0.02525059
## Run 1 stress 0.03844478
## Run 2 stress 0.03083074
## Run 3 stress 0.03569637
## Run 4 stress 0.04275574
## Run 5 stress 0.0466533
## Run 6 stress 0.03777019
## Run 7 stress 0.03211894
## Run 8 stress 0.04068994
## Run 9 stress 0.02528294
## ... Procrustes: rmse 0.001652091 max resid 0.01652819
## Run 10 stress 0.04194455
## Run 11 stress 0.02683015
## Run 12 stress 0.03205697
## Run 13 stress 0.03863817
## Run 14 stress 0.04773489
## Run 15 stress 0.02528829
## ... Procrustes: rmse 0.00150647 max resid 0.01593428
## Run 16 stress 0.04298034
## Run 17 stress 0.02682913
## Run 18 stress 0.03224204
## Run 19 stress 0.03846268
## Run 20 stress 0.03957805
## *** No convergence -- monoMDS stopping criteria:
## 19: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrminO parâmetro “euclidean” fornecido indica que, as distâncias entre pontos foram calculadas através do método de distância euclideana, já o default da função é a distância de Bray-Curtis, usando o parâmetro “bray”.
A qualidade da ordenação no MDS é estimada através do cálculo do stress:
iris.mds$stress## [1] 0.02525059Valores menores do que 0.1 são considerados adequados, como é o nosso caso aqui (0.0253).
O plot desta análise produz apenas uma nuvem de pontos sem especificar qual indivíduo pertence a qual espécie (até por que esta informação não foi incluída na análise). Para associar os pontos (indivíduos) às espécies, vamos salvar as coordenadas dos pontos, nas duas primeiras dimensões em uma nova matriz e, como efetuado na PCA, colorir os indivíduos das três espécies indicando que os 50 primeiros dados (setosa) são vermelhos, os 50 seguintes (versicolor) são azuis e os 50 últimos (virginica), verdes.
As coordenadas dos dois eixos (NMDS1 e NMDS2), estão dentro do objeto iris.mds produzido quando rodamos a análise e de lá que vamos tirar as duas cordenadas:
iris.mds$pointscoord.pontos<-iris.mds$points[1:150,]Posso também colocar uma legenda com esta informação, neste caso escolhi o canto inferior (“bottomright”), para não sobrepor com os dados.
plot(coord.pontos, col=c("red","blue", "green") [as.numeric(iris$Species)], main="NMDS")
legend("bottomright", c("setosa","versicolor","virginica"), pch=1, col=c("red", "blue", "green"))
Você pode ainda ver como as diferentes medidas diferem entre as espécies. Por exemplo vamos plotar as bolinhas das observações com tamanhos diferentes que sejam proporcionais ao valor do comprimento petal (coluna 3 do nosso data frame), isto é dado pelo comando cex=(iris[,3]). onde cex é o tamanho e iris[,3] se refere a terceira coluna de iris. Você poderia também escrever iris$Petal.Length que é a mesma coisa.
coord.pontos<-iris.mds$points[1:150,]plot(coord.pontos, xaxt='n', yaxt='n',col=c("red","blue", "green") [as.numeric(iris$Species)], main="NMDS", cex=(iris[,3]))
legend("bottomright", c("setosa","versicolor","virginica"), pch=1, col=c("red", "blue", "green"))
Você também pode colocar a legenda em qualquer lugar do plot. Para isso, dentro da legenda, você deve incluir o comando locator(1) e clicar com o mouse no local onde você quer colocar a figura:
plot(coord.pontos, xaxt='n', yaxt='n',
col=c("red","blue", "green") [as.numeric(iris$Species)], main="NMDS", cex=(iris[,3]))
legend(locator(1), c("setosa","versicolor","virginica"), pch=1, col=c("red", "blue", "green"))Análise de Variância Multivariada Permutacional (PERMANOVA)
Tanto o PCA quanto o MDS são análises que representam espacialmente a combinação linear (métrica no caso do PCA e não-métrica no caso do MDS), portanto elas não são consideradas como testes estatísticos clássicos por não haver testes de hipótese associados. No caso do MDS, visto acima. as várias medidas tomadas são de grupos de espécies e o resultado se baseia apenas na variação individual, os grupos foram incluídos depois sobre o plot da figura. Será que esta diferença apontada entre os grupos é estatisticamente signficativa? É muito comum a MDS (ou NMDS) ser acompanhada de um teste, este é efetuado através da chamada PERMANOVA, que é um tipo de análise de variância multivariada que usa um procedimento de embaralhamento (permutação) para aplicar um teste de signficância. Ela testa se as similaridades entre grupos é maior do que a similaridade entre indivíduos do mesmo grupo.
No R o comando que faz esta análise é o adonis2 (o comando adonis também pode ser utilizado) do pacote vegan
adonis2(iris[,1:4]~iris$Species, method="euclidean")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = iris[, 1:4] ~ iris$Species, method = "euclidean")
## Df SumOfSqs R2 F Pr(>F)
## iris$Species 2 592.07 0.86894 487.33 0.001 ***
## Residual 147 89.30 0.13106
## Total 149 681.37 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note que tanto o comando, colocando as variáveis dependentes antes do til (** ~ **) como a variável categórica (após o til) é igual ao comando da ANOVA da qual ela é análoga. No caso inclui qual o método (Bray-Curtis), mas não era necessário pois a distância de Bray-Curtis já é o default do comoando adonis. O resultado mostra que a diferença entre espécies que aparece na primeira linha difere signficativamente (p<0.001) com uma grande porcentagem de classificação (R2=0.88). Esta análise foi efetuada a partir de 999 permutações mas este número pode ser alterado acrescentando, entre vírgulas outro número de permutações (permutations=9999, por exemplo).
Todas as análises que vimos acima envolvem apenas uma matriz (ou data.frame) com várias variáveis. No caso da AD e da PERMANOVA, também foram incluídos grupos a serem testados. Entretanto, existem análises de ordenação, denominadas de Análises Diretas ou Análises Constrangidas, que envolvem duas matrizes, uma de variáveis dependentes contínuas e uma de variáveis independentes (contínuas ou não) que condicionam as primeiras. No entanto, neste momento, não vamos explicar este tipo de análise que pode ser encontrada na literatura (CCA, RDA, db-RDA, CAP, etc.) ou mesmo no help do pacote vegan.
Na próxima apostila vamos ver sobre análises multivariadas de classificação.