Neste exemplo, realizamos uma Análise de Componentes Principais (PCA) com base em variáveis ambientais medidas em diferentes unidades amostrais. O objetivo é explorar padrões de correlação entre variáveis e reduzir a dimensionalidade dos dados. Em seguida, usamos os componentes principais obtidos como variáveis preditoras em uma regressão linear múltipla para explicar a diversidade de espécies em cada local.
library(readxl)
library(vegan)
library(dplyr)
dados_bioclim <- readxl::read_excel("Data/dados_bioclimaticos_pca.xlsx")
# Selecionar apenas as variáveis ambientais (excluindo o nome dos sites)
env_vars <- dados_bioclim[, -1]
A função prcomp() realiza a PCA com centramento e escalonamento das variáveis. Isso é importante porque variáveis ambientais podem estar em escalas muito distintas (por exemplo, °C vs mm), o que poderia enviesar a análise se não fossem padronizadas.
# Escalonar os dados (centro e reduz)
pca_bioclim <- prcomp(env_vars, center = TRUE, scale. = TRUE)
summary(pca_bioclim)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.6519 1.1684 1.0910 0.9505 0.69220 0.51092 0.26873
## Proportion of Variance 0.3898 0.1950 0.1700 0.1290 0.06845 0.03729 0.01032
## Cumulative Proportion 0.3898 0.5849 0.7549 0.8839 0.95239 0.98968 1.00000
O scree plot mostra a variância explicada por cada componente principal. A “joelhada” no gráfico sugere o número ideal de componentes a serem retidos. Em geral, mantêm-se os PCs que explicam conjuntamente mais de ~70% da variância.
plot(pca_bioclim, type = "l", main = "Scree plot da PCA")
Os loadings (também chamados de coeficientes, pesos ou cargas) indicam quanto cada variável original contribui para cada componente principal. Eles são fundamentais para interpretar ecologicamente os eixos da PCA.
pca_bioclim$rotation
## PC1 PC2 PC3 PC4
## BIO1_TempMediaAnual 0.57668544 -0.16550915 -0.03487808 -0.005585306
## BIO4_VarTemp -0.16475256 -0.63399152 -0.27726475 0.344236212
## BIO5_TempMax 0.55404413 -0.09213906 0.11537372 -0.039660260
## BIO6_TempMin 0.53326216 -0.19045986 0.07243651 0.119605112
## BIO12_PrecipAnual -0.20078566 -0.58573073 0.39224732 0.239396992
## BIO13_PrecipMax 0.05644168 0.42678423 0.22098970 0.868446922
## BIO14_PrecipMin 0.07386830 0.02591606 -0.83705151 0.232553089
## PC5 PC6 PC7
## BIO1_TempMediaAnual -0.09903333 -0.135898030 -0.78135592
## BIO4_VarTemp 0.61019459 0.007612314 -0.05605040
## BIO5_TempMax 0.18294461 -0.615405029 0.50741441
## BIO6_TempMin -0.06946241 0.748199849 0.30850614
## BIO12_PrecipAnual -0.61177228 -0.164969976 0.06289148
## BIO13_PrecipMax 0.06416385 -0.060080903 -0.06250074
## BIO14_PrecipMin -0.44852970 -0.110102708 0.16073012
O biplot mostra simultaneamente os sites (pontos) e as variáveis ambientais (vetores). A posição dos pontos reflete semelhança multivariada entre locais. As setas vermelhas indicam a direção e intensidade da correlação das variáveis com os eixos principais. Variáveis com vetores longos contribuem mais para o eixo em que estão alinhadas.
# Calcular variância explicada (%)
var_exp <- round(100 * (pca_bioclim$sdev^2 / sum(pca_bioclim$sdev^2)), 1)
# Biplot melhorado
biplot(pca_bioclim,
scale = 0, # mantém escala original dos pontos
cex = 0.8, # tamanho dos rótulos
col = c("gray30", "red"), # cor dos pontos e dos vetores
xlab = paste0("PC1 (", var_exp[1], "%)"),
ylab = paste0("PC2 (", var_exp[2], "%)"))
Aqui extraímos os scores dos dois primeiros componentes, que contêm a maior parte da variação ambiental sintetizada. Esses scores são usados como variáveis preditoras, pois são não correlacionados (ortogonais) e representam combinações lineares das variáveis originais.
# Obter os scores dos dois primeiros PCs
scores_pca <- as.data.frame(pca_bioclim$x[, 1:2])
colnames(scores_pca) <- c("PC1", "PC2")
Aqui carregamos os dados de comunidade para calcular a diversidade.
# 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'
A diversidade é calculada usando a função diversity()
do
pacote vegan
. O Default dessa função é calcular a
diversidade de Shannon.
scores_pca$Diversidade <- diversity(comunidade) # Div de Shannon
Agora podemos fazer nossa regressão multipla e interpretar os resultados.
modelo <- lm(Diversidade ~ PC1 + PC2, data = scores_pca)
summary(modelo)
##
## Call:
## lm(formula = Diversidade ~ PC1 + PC2, data = scores_pca)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16278 -0.07039 0.01250 0.04525 0.20192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.361848 0.016030 209.726 <2e-16 ***
## PC1 0.001413 0.009870 0.143 0.887
## PC2 -0.007562 0.013954 -0.542 0.592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0878 on 27 degrees of freedom
## Multiple R-squared: 0.0115, Adjusted R-squared: -0.06172
## F-statistic: 0.1571 on 2 and 27 DF, p-value: 0.8554
plot(Diversidade ~ PC1 + PC2, data = scores_pca)
A PCA permite explorar os principais gradientes ambientais presentes nos dados e sintetizá-los em componentes ortogonais. Usar esses componentes como preditores em um modelo de regressão reduz problemas de multicolinearidade e permite testar se combinações de fatores ambientais estão associadas a padrões de diversidade. Essa estratégia é particularmente útil quando se deseja reduzir a complexidade dos dados antes de modelagens mais formais.