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 rio¹. 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 ambientais², 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

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)

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

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))

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.