10 min read

Explorando dados ambientais em análises multivariadas

Neste roteiro, pretendo apresentar e explorar algumas ferramentas de visualização e exploração de dados multivariados, em especial variáveis ambientais coletadas em diferentes locais. Estes são dados muito utilizados em Ecologia, e portanto merecem nossa atenção, mas isso não quer dizer que as técnicas aqui apresentadas não sejam úteis a outros tipos de dados multivariados, muito pelo contrário!

Para esse roteiro, eu utilizei os livros de Borcard, Gillet, and Legendre (2011) (capítulo 2), Manly and Navarro Alberto (2017) (Capítulo 3), Everitt and Hothorn (2011) (Capítulo 2). Vale a pena a leitura destes capítulos para aprofundamento!

Explorando variáveis ambientais ao longo de um rio

Vamos utilizar os dados disponibilizados em Borcard, Gillet, and Legendre (2011) (baixe a planilha de dados ambientais nesse link e a de dados das coordenadas geográficas dos pontos nesse link), com dados ambientais coletados em diferentes pontos de um rio.

env <- read.table("DoubsEnv.csv", header=T, sep=",")
spa <- read.table("DoubsSpa.csv", header=T, sep=",", row.names = 1)
Table 1: Visualização das primeiras linhas dos dados ambientais.
X das alt pen deb pH dur pho nit amm oxy dbo
1 0.3 934 48.0 0.84 7.9 45 0.01 0.20 0.00 12.2 2.7
2 2.2 932 3.0 1.00 8.0 40 0.02 0.20 0.10 10.3 1.9
3 10.2 914 3.7 1.80 8.3 52 0.05 0.22 0.05 10.5 3.5
4 18.5 854 3.2 2.53 8.0 72 0.10 0.21 0.00 11.0 1.3
5 21.5 849 2.3 2.64 8.1 84 0.38 0.52 0.20 8.0 6.2
6 32.4 846 3.2 2.86 7.9 60 0.20 0.15 0.00 10.2 5.3

Vamos observar a planilha de dados ambientais:

str(env)
## 'data.frame':    30 obs. of  12 variables:
##  $ X  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ das: num  0.3 2.2 10.2 18.5 21.5 32.4 36.8 49.1 70.5 99 ...
##  $ alt: int  934 932 914 854 849 846 841 792 752 617 ...
##  $ pen: num  48 3 3.7 3.2 2.3 3.2 6.6 2.5 1.2 9.9 ...
##  $ deb: num  0.84 1 1.8 2.53 2.64 2.86 4 1.3 4.8 10 ...
##  $ pH : num  7.9 8 8.3 8 8.1 7.9 8.1 8.1 8 7.7 ...
##  $ dur: int  45 40 52 72 84 60 88 94 90 82 ...
##  $ pho: num  0.01 0.02 0.05 0.1 0.38 0.2 0.07 0.2 0.3 0.06 ...
##  $ nit: num  0.2 0.2 0.22 0.21 0.52 0.15 0.15 0.41 0.82 0.75 ...
##  $ amm: num  0 0.1 0.05 0 0.2 0 0 0.12 0.12 0.01 ...
##  $ oxy: num  12.2 10.3 10.5 11 8 10.2 11.1 7 7.2 10 ...
##  $ dbo: num  2.7 1.9 3.5 1.3 6.2 5.3 2.2 8.1 5.2 4.3 ...

Também vale usar as funções:

library(skimr)
skim(env)
# e 
summary(env)

Mapas das variáveis ambientais

Se você tem as coordenadas geográficas das variáveis, pode ser muito útil plotar as variáveis ambientais no mapa.

No nosso exemplo, utilizei o script disponível em Borcard, Gillet, and Legendre (2011) para plotar algumas variáveis ao longo do rio1. Neste mapa estamos utilizando uma técnica chamada “bubble maps,” no qual o tamanho dos círculos é proporcional ao valor da variável analisada.

par(mfrow=c(2,2))
plot(spa, asp=1, main="Altitude", pch=21, col="white", bg="red",
    cex=5*env$alt/max(env$alt), xlab="coord x", ylab="coord y")
lines(spa, col="light blue")
plot(spa, asp=1, main="Discharge", pch=21, col="white", bg="blue",
    cex=5*env$deb/max(env$deb), xlab="coord x", ylab="coord y")
lines(spa, col="light blue")
plot(spa, asp=1, main="Oxygen", pch=21, col="white", bg="green3",
    cex=5*env$oxy/max(env$oxy), xlab="coord x", ylab="coord y")
lines(spa, col="light blue")
plot(spa, asp=1, main="Nitrate", pch=21, col="white", bg="brown",
    cex=5*env$nit/max(env$nit), xlab="coord x", ylab="coord y")
lines(spa, col="light blue")

Abaixo vai o código para os mesmos plots usando o pacote ggplot2 para quem preferir:

library(ggplot2)
library(cowplot) # melhoras no acabamento do gráfico
library(patchwork) # compõe vários gráficos em um painel

ggplot(spa, aes(x=x, y=y), col="red") +
  geom_point(size=10*env$alt/max(env$alt), col="red") +
  ggtitle("Altitude") +
  geom_path()  +
ggplot(spa, aes(x=x, y=y)) +
  geom_point(size=10*env$deb/max(env$deb), col="blue") +
  ggtitle("Discharge") +
  geom_path() +
ggplot(spa, aes(x=x, y=y)) +
  geom_point(size=10*env$oxy/max(env$oxy), col="green") +
  ggtitle("Oxygen") +
  geom_path() +
ggplot(spa, aes(x=x, y=y)) +
  geom_point(size=10*env$nit/max(env$nit), col="brown") +
  ggtitle("Nitrate") +
  geom_path()  +
plot_layout(ncol=2)

Qual destas variávies mostra um gradiente montante-jusante? Como você poderia explicar o padrão das outras variáveis? Faça os gráficos para as outras variáveis.

Agora vamos examinar a variação de algumas variáveis ao longo do rio:

par(mfrow=c(2,2))
plot(env$das, env$alt, type="l", xlab="Distância (km)", 
    ylab="Altitude (m)", col="red", main="Altitude")
plot(env$das, env$deb, type="l", xlab="Distância (km)", 
    ylab="Discharge (m3/s)", col="blue", main="Discharge")
plot(env$das, env$oxy, type="l", xlab="Distância (km)", 
    ylab="Oxygen (mg/L)", col="green3", main="Oxygen")
plot(env$das, env$nit, type="l", xlab="Distância (km)", 
    ylab="Nitrate (mg/L)", col="brown", main="Nitrate")

Correlações entre as variáveis

Para explorar graficamente a relação bivariada entre as variáveis ambientais2, podemos usar a função pairs(), que vai desenhar todas as combinações de variáveis duas a duas.

pairs(env[,3:7],
    main="Relações entre variáveis")

Entretanto, esta função pode ser melhorada para nos dar mais informação em um só gráfico. Para fazer o plot abaixo, baixe esse código(disponibilizado por Borcard et al. 2011, de Francois Gillet) e rode a função panelutils ou use a função source para carrega-la no seu workspace.

source("panelutils.R")
pairs(env[,3:7], diag.panel = panel.hist, upper.panel = panel.cor,
      lower.panel = panel.smooth,   main="Relações entre variáveis")

Esse gráfico nos dá os histogramas de cada variável na diagonal, os gráficos de dispersão nos painéis debaixo e os valores das correlações com seu níveis de significância (astericos abaixo do valor da correlação). O histogramas são úteis para se avaliar a distribuição dos dados, e os gráficos de dispersão para ver a forma das relações entre as variáveis (relações lineares, monotônicas, etc.)

Existe também o pacote corrplot que faz gráficos de painéis bem legais. Vale a pena olhar a vinheta do pacote (está no help)! Abaixo vai um exemplo, no qual ordenamos as variáveis fazendo um agrupamento hierárquico (você escolhe o número de grupos e também o métodos de agrupamento) para ver quais variáveis são mais relacionadas entre si.

library(corrplot)

M <- cor(env[,2:11]) # primeiro é preciso fazer as correlações
corrplot(M, method="circle", insig = "p-value", order="hclust", addrect=4)

A principal desvantage do painel de plots bivariados é que aspectos dos dados que aparecem apenas quando 3 ou mais variáveis são consideradas em conjunto não ficam aparentes.

Estes gráficos também só são adequados para comparar as variáveis e não os locais (objetos ou unidades amostrais, as linhas da matriz).

Visualizando a comparação entre locais

As faces de Chernoff

Existem algumas alternativas para visualizar a relação entre os locais de acordo com as variáveis ambientais.

Uma destas são as faces de Chernoff, com a função faces (ou faces1) disponível no pacoteTeachingDemos`. Nela podemos comparar os locais de coleta de dados olhando para as faces de cada local - faces mais mais parecidas significam locais mais parecidos. A ideia é usar as variáveis como definidoras das características das faces, logo uma das variáveis vai controlar o tamanho dos olhos, outra o formato da cabeça, o tamnho da boca, assim por diante. É uma forma bem divertida de explorar os dados!

São 15 parâmetros de implementação das faces na função faces2, usamos o argumento which se quisermos especificar qual parâmetro para qual variável:

  1. Altura da face
  2. Largura da face
  3. Formato da face
  4. Altura da boca
  5. Largura da boca
  6. Curvatura da boca
  7. Altura dos olhos
  8. Largura dos olhos
  9. Altura dos cabelos
  10. Largura dos cabelos
  11. Estilo dos cabelos
  12. Altura do nariz
  13. Largura do nariz
  14. Largura das orelhas
  15. Altura das orelhas
library(TeachingDemos)
canmat <- data.matrix(env[,2:7])
faces2(canmat, labels = row.names(env), which= c(3,5,6,8,11,13))

#outro pacote e função para desenhar as faces:
library(aplpack)
fa <- faces(env[,2:7],plot=FALSE)
plot(fa,face.type=2)

Aqui vemos novamente que pontos de amostragem no rio parecem formar um gradiente de condições ambientais, ou seja, os pontos mais próximos possuem as faces mais parecidas.

As estrelas

As estrelas são uma forma similiar de observar a relação entre os locais de amostragm baseando-se nas variáveis ambientais. A função stars está no pacote base do R.

stars(env[,3:8], labels = row.names(env), draw.segments = TRUE, nrow=5, ncol=7,
      key.loc = c(10,2))

Eu particularmente gosto muito desta visualização, pois nos permite observar as variáveis mais importantes (com os maiores valores) em cada local e os locais mais parecidos ao mesmo tempo. Por exemplo, o ponto 1 está em alta altitude e declividade (pen), e seguindo o curso do rio (indo em direção ao local 30), temos que a descarga (deb) aumenta junto com a concentração de cálcio (dur).

Estes plots são outra forma de observar o gradiente ambiental das variáveis, como feito nos mapas inicialmente.

Perfil das Variáveis

Essa é uma forma de representar os locais (objetos) todos de uma vez, e é especialmente útil quando temos variáveis em uma mesma unidade, por exemplo dados de diferentes medidas morfológicas de animais, o que não é o caso dos nossos dados. Porque os dados ambientais são dimensionalmente heterogêneos (diferentes unidades e escalas), precisamos que eles sejam padronizados para z-score (média zero e unidade de variância - centraliza e escalona). Assim, removemos a dimensionalidade dos dados, permitindo as comparações.

env.p <- scale(env[,2:12], center = T, scale = T)
env.p <- data.frame(env$X, env.p)
# para re-organizar os dados vamos usar o pacote tidyverse
library(tidyverse)
env.p %>% 
  gather("variaveis", "valores", 2:12) %>% # re-organizando os dados
  filter(env.X %in% seq(1,30,3)) %>% # pegando alguns locais apenas para visualização
ggplot(aes(x=variaveis, y=valores, group=env.X, col=env.X)) + #plot com ggplot2
  geom_line() +
  scale_color_continuous("locais", breaks = c(1,10,20,28), labels=c(1,10,20,30))

Além de compararmos os locais em relação às variáveis ambientais, também conseguimos ver a variablidade dos dados. Por exemplo, vemos que a variável pen tem pouca variabilidade, mas possui um valor extremo. O pH em comparação às outras variáveis também tem pouca variabilidade. Certas análises multivariadas tem como premissa a variância homogênea entre variáveis, e esse plot de perfil pode ajudar a diagnosticar isso.

Boxplots bivariado para comparar variáveis

Everitt & Hothorn 2011 apresentam uma maneira bem interessante de observar as relações entre 2 variáveis, o boxplot bivariado. Estes boxplots são também extremamente úteis na identificação de valores extremos (outliers).

Os boxplots bivariados são baseados nos cálculos “robustos” das medidas de locação, escala e correlação e consiste em um par de elipses concêntricas, uma das quais inclui 50% dos dados e a outra que delinea o limite dos dados e deixa de fora possíveis outliers. Também são desenhadas linhas de regressão de x em y e de y em x, com sua interseção mostrando o parâmetro de locação (média das duas variáveis). O ângulo entre as linhas de regressão serão tão pequenos quanto maior for a correlação entre as variáveis.

Vamos usar a função bvbox do pacote MVA que acompanha o livro de Everitt & Hothorn 2011, para observar as relações entre altitute e as outras variáveis:

library(MVA)
par(mfrow=c(3,3))
bvbox(env[,c(3,4)], xlab = colnames(env)[3], ylab = colnames(env)[4], pch="")
text(env[,3], env[,4], env$X)
bvbox(env[,c(3,5)], xlab = colnames(env)[3], ylab = colnames(env)[5], pch="")
text(env[,3], env[,5], env$X)
bvbox(env[,c(3,6)], xlab = colnames(env)[3], ylab = colnames(env)[6], pch="")
text(env[,3], env[,6], env$X)
bvbox(env[,c(3,7)], xlab = colnames(env)[3], ylab = colnames(env)[7], pch="")
text(env[,3], env[,7], env$X)
bvbox(env[,c(3,8)], xlab = colnames(env)[3], ylab = colnames(env)[8], pch="")
text(env[,3], env[,8], env$X)
bvbox(env[,c(3,9)], xlab = colnames(env)[3], ylab = colnames(env)[9], pch="")
text(env[,3], env[,9], env$X)
bvbox(env[,c(3,10)], xlab = colnames(env)[3], ylab = colnames(env)[10], pch="")
text(env[,3], env[,10], env$X)
bvbox(env[,c(3,11)], xlab = colnames(env)[3], ylab = colnames(env)[11], pch="")
text(env[,3], env[,11], env$X)
bvbox(env[,c(3,12)], xlab = colnames(env)[3], ylab = colnames(env)[12], pch="")
text(env[,3], env[,12], env$X)
par(mfrow=c(1,1))

Nestes plots, já podemos observar que existem pares de variáveis com outliers bivariados em alguns locais (identificados com os números nos plots). Podemos ver que, quando as variáveis não são correlacionadas, o formato do boxplot bivariado é praticamente um círculo, como acontece entre altitude e pH.

Referências consultadas

Borcard, D, F Gillet, and P Legendre. 2011. Numerical Ecology with R.
Everitt, B, and T Hothorn. 2011. An Introduction to Applied Multivariate Analysis with R.
Manly, Bryan F. J., and Jorge A. Navarro Alberto. 2017. Multivariate Statistical Methods: A Primer. Fourth edition. Boca Raton: CRC Press.

  1. Mantive os nomes das variáveis em inglês para não fazer traduções erradas, pois isso não importa pra gente neste momento.↩︎

  2. Para facilitar a visualização, à vezes uso apenas um subconjunto de variáveis.↩︎