Estatística multivariada no R

Este exemplo realiza uma série de análises de comunidades no R, como exemplos introdutórios.

Pacotes necessários

# Instale os pacotes, se ainda não tiver:
install.packages(c("readxl", "vegan", "cluster", "dendextend", "pheatmap"))

Carregando os pacotes.

library(readxl)      # para ler arquivos Excel
library(vegan)       # para cálculo da dissimilaridade (Bray-Curtis)
library(cluster)     # funções de clustering
library(dendextend)  # para melhorar a visualização do dendrograma
library(pheatmap)    # plotar as matrizes

Carregar os dados

# Substitua o caminho abaixo pelo caminho do seu arquivo
dados <- read_excel("Data/dado_comunidade_gradiente.xlsx")

# Separar variável categórica (Habitat) e dados de abundância
habitat <- as.factor(dados$Habitat)
sites <- dados$Sites

comunidade <- as.data.frame(dados[, -(1:2)])  # remove a coluna 'Habitat' e'site'

Cálculo da matriz de dissimilaridade (Bray-Curtis)

Use a função vegdist do pacote vegan para calcular as distâncias.

dist_brays <- vegdist(comunidade, method = "bray")

Agrupamento hierárquico (UPGMA)

Rodar a análise.

upgma <- hclust(dist_brays, method = "average")

Você pode visualizar a matriz de distância para ver se está tudo ok.

# Converter em matriz simétrica
mat_dist <- as.matrix(dist_brays)
mat_dist[lower.tri(mat_dist)] <- NA # Eliminar a parte superior
rownames(mat_dist) <- colnames(mat_dist) <- paste(habitat, sites, sep = "_")
heatmap(mat_dist,
        Rowv = NA,
        Colv = NA,
        main = "Matriz de dissimilaridade (Bray-Curtis)",
        scale = "none",
        margins = c(8, 2))

Visualização do dendrograma colorido por tipo de habitat.

Para o plot há várias opções, aqui usaremos uma simples usando as funções base do R.

# Converter para objeto dendrograma
dend <- as.dendrogram(upgma)

# Atribuir cores aos rótulos com base no habitat
labels_colors(dend) <- as.numeric(as.factor(habitat))[order.dendrogram(dend)]

# Plotar
plot(dend, main = "UPGMA - Bray-Curtis", ylab = "Dissimilaridade")
legend("topright", legend = levels(as.factor(habitat)),
       fill = 1:3, title = "Habitat", border = NA)

Cálculo da matriz cofenética e sua correlação

O coeficiente de roelação cofenética mede a correlação entre as distâncias originais (ex. Bray-Curtis) e as distâncias representadas no dendrograma.Valores próximos de 1 indicam boa fidelidade da estrutura hierárquica ao padrão original de dissimilaridade. Valores > 0.75 geralmente são considerados aceitáveis.

coph <- cophenetic(upgma)
cor(dist_brays, coph)  # correlação cophenética
## [1] 0.9694109

PCoA

A PCoA (Principal Coordinates Analysis) é uma técnica de ordenação que visa representar as relações de dissimilaridade entre objetos (neste caso, comunidades amostradas em diferentes locais) em um espaço de menor dimensão, normalmente 2D ou 3D.

pcoa_result <- cmdscale(dist_brays, k = 2, eig = TRUE)  # k = número de eixos

Coordenadas

As colunas (PCoA1, PCoA2, etc.) são os eixos ordinais — construídos para preservar ao máximo as distâncias ecológicas (ex: Bray-Curtis) entre os pontos.

As coordenadas indicam a posição relativa de cada local nesse novo espaço (esses valores também são chamados de scores), onde: - Locais mais próximos no gráfico têm composições de espécies mais semelhantes; - Locais mais distantes apresentam maior dissimilaridade.

Ao contrário de PCA, a PCoA pode trabalhar com qualquer matriz de dissimilaridade, mesmo quando ela não é euclidiana (como Bray-Curtis ou Jaccard).

coord <- as.data.frame(pcoa_result$points)
colnames(coord) <- c("PCoA1", "PCoA2")
coord$Habitat <- habitat

Percentual de variação explicada (opcional, para eixos)

A função cmdscale(..., eig = TRUE) retorna os autovalores (eigenvalues) correspondentes aos eixos. Esses valores refletem quanto da estrutura original das distâncias é preservada por cada eixo. - O primeiro eixo (PCoA1) é o que melhor preserva as distâncias originais. - O segundo eixo (PCoA2) é ortogonal ao primeiro e explica o máximo restante da variação. - A soma da variância explicada pelos dois eixos indica o grau de fidelidade da visualização 2D em relação ao espaço original.

eig_values <- pcoa_result$eig
var_exp <- round(100 * eig_values[1:2] / sum(eig_values[eig_values > 0]), 1)

Vizualização

O resultado mostra a posição relativa das unidades amostrais no espaço multivariado, preservando ao máximo as distâncias originais. É particularmente útil para visualizar padrões de substituição de espécies ao longo de gradientes.

# Cores por tipo de habitat
cores <- c("Seca" = "orange", "Intermediaria" = "skyblue", "Umida" = "darkgreen")

# Plot
plot(coord$PCoA1, coord$PCoA2,
     col = cores[coord$Habitat],
     pch = 19, cex = 1.5,
     xlab = paste0("PCoA1 (", var_exp[1], "%)"),
     ylab = paste0("PCoA2 (", var_exp[2], "%)"),
     main = "Análise de Coordenadas Principais (PCoA)")

legend("bottomleft", legend = names(cores), col = cores, pch = 19, title = "Habitat")

Análise de Correspondência (CA)

A Análise de Correspondência (CA) é uma técnica de ordenação multivariada exploratória indicada para analisar tabelas de contagem ou abundância de espécies, especialmente quando os dados são composicionais e esparsos (com muitos zeros). Ela projeta simultaneamente unidades amostrais e espécies em um mesmo espaço ordinado, utilizando a distância qui-quadrado como métrica, o que permite identificar padrões de coocorrência e associação entre linhas e colunas. Diferentemente da PCA, a CA não requer normalização dos dados e é particularmente útil quando se trabalha com dados não gaussianos. Suas principais vantagens incluem a facilidade de interpretação gráfica em contextos ecológicos e a robustez diante de dados de abundância, mas pode ser sensível à dominância de espécies muito frequentes e não permite a inclusão direta de variáveis ambientais, sendo, nesse caso, preferível a CCA. Em geral, a CA é recomendada quando se busca explorar padrões de composição entre amostras e espécies, sem a necessidade de hipóteses rígidas sobre distribuição ou linearidade dos dados.

ca <- cca(comunidade)
# Obter scores de amostras (sites) e espécies
sites_coord <- scores(ca, display = "sites")
species_coord <- scores(ca, display = "species")

# Plot vazio com limites ajustados
plot(ca, type = "n", scaling = 2,
     main = "Análise de Correspondência (CA)",
     xlab = paste0("CA1 (", round(ca$CA$eig[1] / sum(ca$CA$eig) * 100, 1), "%)"),
     ylab = paste0("CA2 (", round(ca$CA$eig[2] / sum(ca$CA$eig) * 100, 1), "%)"))

# Adicionar os locais (amostras)
points(sites_coord, pch = 21, bg = cores[habitat], cex = 1.5)

# Adicionar nomes das espécies com cor neutra
text(species_coord, labels = rownames(species_coord), col = "gray30", cex = 0.6)

# Adicionar legenda
legend("topright", legend = unique(habitat),
       fill = cores[unique(habitat)], border = "black",
       title = "Habitat", cex = 0.9)

PERMANOVA

Para testar se grupos de amostras (e.g., tipos de habitat: “Seca”, “Intermediária”, “Úmida”) diferem estatisticamente em termos de composição de espécies, é necessário usar um teste como uma PERMANOVA. PERMANOVA é um teste não-paramétrico baseado em permutações que compara as distâncias intra e entre grupos usando uma matriz de dissimilaridade (como Bray-Curtis). Implementado em R via a função adonis2() do pacote vegan.

resultado <- adonis2(dist_brays ~ habitat, permutations = 999)
print(resultado)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist_brays ~ habitat, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## Model     2   3.4816 0.70527 32.304  0.001 ***
## Residual 27   1.4550 0.29473                  
## Total    29   4.9366 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A PERMANOVA assume homogeneidade da dispersão dos grupos (distâncias intra-grupo semelhantes). Se a dispersão for significativamente diferente, a significância da PERMANOVA pode refletir isso — não necessariamente diferença real nas médias de composição. Isso pode ser testado com:

mod <- betadisper(dist_brays, habitat)
anova(mod)  # teste de homogeneidade
## Analysis of Variance Table
## 
## Response: Distances
##           Df    Sum Sq    Mean Sq F value Pr(>F)
## Groups     2 0.0018925 0.00094626  1.0145  0.376
## Residuals 27 0.0251835 0.00093272

CCA

A CCA é uma técnica de ordenação constrangida, ou seja, ela relaciona explicitamente a variação na composição de espécies a um conjunto de variáveis ambientais. Diferente da CA e da PCoA, que são exploratórias e baseadas apenas nos dados de comunidade, a CCA permite testar e visualizar como as variáveis ambientais estruturam as comunidades.

Aqui vamos carregar um conjunto de variáveis ambientais simuladas e utilizá-las para explicar os padrões de composição de espécies.

# Carregar dados ambientais
dados_env <- read_excel("Data/dados_ambientais_para_cca.xlsx")

# Filtrar só variáveis numéricas para a CCA
env_vars <- dados_env[, c("Umidade", "Sombra", "pH_solo")]

Com os dados carregados, executamos a CCA relacionando a matriz de abundância à matriz ambiental. O objeto resultante contém os eixos canônicos que maximizam a variação da composição explicada pelas variáveis ambientais.

cca_result <- cca(comunidade ~ ., data = env_vars)
summary(cca_result)
## 
## Call:
## cca(formula = comunidade ~ Umidade + Sombra + pH_solo, data = env_vars) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total         1.07181    1.00000
## Constrained   0.07415    0.06919
## Unconstrained 0.99766    0.93081
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CCA1    CCA2    CCA3    CA1    CA2     CA3     CA4
## Eigenvalue            0.03745 0.02395 0.01276 0.3654 0.3146 0.04070 0.02939
## Proportion Explained  0.03494 0.02234 0.01191 0.3409 0.2935 0.03797 0.02742
## Cumulative Proportion 0.03494 0.05728 0.06919 0.4101 0.7036 0.74154 0.76897
##                           CA5     CA6     CA7     CA8     CA9    CA10    CA11
## Eigenvalue            0.02485 0.02467 0.02265 0.02223 0.01758 0.01684 0.01535
## Proportion Explained  0.02319 0.02301 0.02113 0.02074 0.01640 0.01571 0.01432
## Cumulative Proportion 0.79215 0.81516 0.83629 0.85703 0.87344 0.88915 0.90347
##                          CA12    CA13    CA14     CA15     CA16     CA17
## Eigenvalue            0.01450 0.01275 0.01180 0.010363 0.009960 0.008487
## Proportion Explained  0.01353 0.01189 0.01101 0.009669 0.009292 0.007918
## Cumulative Proportion 0.91700 0.92889 0.93990 0.949573 0.958865 0.966783
##                           CA18     CA19     CA20     CA21     CA22     CA23
## Eigenvalue            0.006978 0.006397 0.004744 0.004468 0.004346 0.003375
## Proportion Explained  0.006511 0.005968 0.004426 0.004168 0.004055 0.003149
## Cumulative Proportion 0.973294 0.979262 0.983688 0.987856 0.991911 0.995060
##                           CA24     CA25      CA26
## Eigenvalue            0.002615 0.001648 0.0010310
## Proportion Explained  0.002440 0.001538 0.0009619
## Cumulative Proportion 0.997500 0.999038 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                          CCA1    CCA2    CCA3
## Eigenvalue            0.03745 0.02395 0.01276
## Proportion Explained  0.50498 0.32294 0.17208
## Cumulative Proportion 0.50498 0.82792 1.00000

Agora ajustamos os vetores ambientais com a função envfit(), que testa a significância estatística da correlação entre cada variável e os eixos da ordenação. O r² representa o grau de correlação e o p-valor é obtido por permutação. É essencial observar quais variáveis são significativamente associadas à ordenação (geralmente p < 0.05), pois isso indica que essas variáveis estruturam a composição das espécies.

fit_env <- envfit(cca_result, env_vars, permutations = 999)
fit_env
## 
## ***VECTORS
## 
##             CCA1     CCA2     r2 Pr(>r)  
## Umidade -0.84607  0.53307 0.1787  0.067 .
## Sombra  -0.89641  0.44323 0.1565  0.093 .
## pH_solo  0.62679 -0.77919 0.1943  0.052 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

Agora vamos ver o gráfico:

# Coordenadas escaladas
ord_scores <- scores(cca_result, display = "sites", scaling = 2)

# Plot base
plot(cca_result, type = "n", scaling = 2,
     main = "CCA - Composição de espécies vs variáveis ambientais")
points(ord_scores, col = cores[habitat], pch = 21, bg = cores[habitat], cex = 1.5)

# Vetores ambientais (sem scaling)
plot(fit_env, p.max = 0.5, col = "red", add = TRUE)

# Elipses
ordiellipse(cca_result, groups = habitat, kind = "ehull", draw = "polygon",
            border = cores,
            col = adjustcolor(cores, alpha.f = 0.3),
            label = TRUE, scaling = 2)

legend("topright", legend = names(cores),
       fill = cores, border = "black", title = "Habitat", cex = 0.9)

Esse gráfico permite observar, por exemplo, se determinados habitats estão associados a certos valores de umidade, sombra ou pH do solo, e se há sobreposição ou separação entre os habitats em termos de composição. A direção dos vetores mostra para onde os valores das variáveis aumentam, e sua projeção sobre os grupos ajuda a entender quais variáveis ambientais influenciam mais fortemente cada tipo de comunidade.