Impacto da restauração de florestas ciliares na recuperação do solo
Autor¹: Antonio Barros; Autor²: João Netto; Autor³: Leonardo Takaku
2022-07-15
Contextualização
Pesquisou-se o impacto da restauração de florestas ciliares na recuperação do solo, avaliando seus atributos químicos. Para isto, ele conduziu uma amostragem de solo em oito áreas de restauro na região de córrego do feijão sendo, três áreas em restauros novos, com cinco anos de idade e seis áreas restauradas a 20 anos. Amostras foram conduzidas também em outros dois tipos de uso do solo usados como áreas de referência, sendo quatro áreas de pastagem e quatro áreas de remanescentes florestais. Em todas as áreas, ele marcou parcelas 10 × 10 m e dentro de cada parcela ele retirou seis subamostras de 0 a 20cm de profundidade (figura ao lado). Ele então homogeneizou as subamostras em uma bandeja e retirou 300g de solo, de modo que cada amostra representasse uma parcela que era sua unidade de observação. Estime a SB (soma de bases), CTC (Capacidade de Troca Catiônica) e o V% (Saturação de bases) , estime as correlações entre os parâmetros químicos do solo e conduza uma Análise de Componentes Principais.
Importando e Tratando o Conjunto de Dados - Qualidade do Solo
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows# Diretório de trabalho
setwd("C:/Users/Ultron/Documents/Restauro")
# Importar dados de qualidade do solo
rSolo <- read.csv('Data.csv', sep = ';')
# Renomear os fatores
rSolo$Tipo = factor(x = rSolo$Tipo,
levels = c('RFl', 'ResAnt', 'Pas', 'ResNovo'),
labels = c('Remanescente Florestal', 'Restauro Florestal Antigo',
'Pastagem', 'Restauro Florestal Novo'))# Lendo o conjunto de dados
rSolo## Area Tipo pH MO N P K Ca Mg Al
## 1 C1 Remanescente Florestal 3.8 68 5600 22 2.5 7 8 135
## 2 C1 Remanescente Florestal 3.8 59 5320 24 2.6 9 11 149
## 3 C1 Remanescente Florestal 3.6 60 2800 16 1.6 9 4 121
## 4 C2 Remanescente Florestal 3.4 77 1330 10 2.0 2 1 109
## 5 C3 Remanescente Florestal 5.1 68 1792 17 1.6 79 24 22
## 6 C3 Remanescente Florestal 4.4 59 3780 14 1.7 29 12 58
## 7 C3 Remanescente Florestal 5.0 56 700 15 1.9 57 22 34
## 8 C3 Remanescente Florestal 4.7 70 3080 11 1.6 35 16 58
## 9 C4 Remanescente Florestal 6.1 56 3150 17 2.0 86 36 15
## 10 C4 Remanescente Florestal 5.8 55 700 13 1.6 72 19 16
## 11 C4 Remanescente Florestal 5.0 45 4900 13 1.7 38 17 25
## 12 C4 Remanescente Florestal 4.9 54 4200 15 1.6 53 20 31
## 13 R6 Restauro Florestal Antigo 3.7 32 2695 6 2.1 4 3 64
## 14 R6 Restauro Florestal Antigo 3.8 51 2660 8 3.1 7 6 71
## 15 R1 Restauro Florestal Antigo 4.6 47 2499 11 1.4 15 10 58
## 16 R1 Restauro Florestal Antigo 4.5 49 1932 13 1.5 20 12 64
## 17 R1 Restauro Florestal Antigo 4.9 45 2660 11 1.6 23 13 52
## 18 R1 Restauro Florestal Antigo 4.1 45 2660 8 2.0 7 5 71
## 19 R2 Restauro Florestal Antigo 4.9 53 1400 10 1.5 26 20 42
## 20 R2 Restauro Florestal Antigo 4.8 59 1064 8 1.6 25 18 58
## 21 R2 Restauro Florestal Antigo 4.9 52 1015 6 1.1 18 16 52
## 22 R3 Restauro Florestal Antigo 4.6 55 2625 3 1.4 14 10 64
## 23 R3 Restauro Florestal Antigo 4.6 54 2800 7 1.6 18 12 64
## 24 R3 Restauro Florestal Antigo 4.5 52 2170 7 1.2 16 15 64
## 25 R3 Restauro Florestal Antigo 4.4 51 2275 8 1.5 9 9 58
## 26 R4 Restauro Florestal Antigo 4.8 58 2450 12 2.0 28 24 42
## 27 R4 Restauro Florestal Antigo 4.7 49 3150 12 2.0 27 19 47
## 28 R4 Restauro Florestal Antigo 4.8 50 3115 10 1.8 25 23 42
## 29 R5 Restauro Florestal Antigo 5.0 53 2485 9 1.6 27 15 31
## 30 R5 Restauro Florestal Antigo 5.2 50 1953 15 1.7 34 16 31
## 31 R5 Restauro Florestal Antigo 4.9 49 1064 9 2.1 27 15 34
## 32 R5 Restauro Florestal Antigo 5.0 50 1820 25 2.1 32 20 31
## 33 A1 Pastagem 4.3 30 1085 7 1.6 6 5 22
## 34 A1 Pastagem 4.4 27 1120 7 1.0 8 7 31
## 35 A1 Pastagem 4.2 30 1190 12 1.4 6 5 31
## 36 A1 Pastagem 4.2 20 1400 3 2.1 4 4 34
## 37 A2 Pastagem 4.1 23 1050 6 1.2 4 4 28
## 38 A2 Pastagem 4.3 19 1015 10 1.4 4 4 28
## 39 A2 Pastagem 4.1 19 1365 9 1.3 4 4 28
## 40 A2 Pastagem 3.9 25 1820 7 1.2 3 3 34
## 41 A3 Pastagem 4.3 32 1855 9 1.2 16 5 34
## 42 A3 Pastagem 4.4 32 1715 8 1.1 17 6 31
## 43 A4 Pastagem 3.7 23 2800 3 1.3 1 1 31
## 44 A4 Pastagem 3.4 33 1176 7 1.4 1 1 52
## 45 R8 Restauro Florestal Novo 3.6 52 2100 12 1.4 2 3 71
## 46 R8 Restauro Florestal Novo 3.5 46 980 9 1.6 1 1 71
## 47 R9 Restauro Florestal Novo 4.6 40 1715 11 1.0 15 7 34
## 48 R9 Restauro Florestal Novo 4.6 39 1540 7 1.6 16 10 31
## 49 R9 Restauro Florestal Novo 4.4 34 1505 3 1.3 9 7 31
## 50 R9 Restauro Florestal Novo 4.5 36 1190 3 1.2 8 7 28
## 51 R7 Restauro Florestal Novo 3.8 27 2590 7 1.0 5 3 42
## 52 R7 Restauro Florestal Novo 3.7 30 1085 3 1.1 4 2 52
## 53 R7 Restauro Florestal Novo 4.1 33 1680 3 1.3 6 4 31
## 54 R7 Restauro Florestal Novo 3.8 36 1225 6 1.3 6 2 79Capacidade de Troca de Cátions (CTC), Soma de Bases (SB) e Saturação por Bases (V%)
Em razão da superfície eletricamente carregada que apresentam as argilas coloidais, as substâncias húmicas e os sesquióxidos de ferro e alumínio (principais componentes da fração mineral dos solos sob condições tropicais), os íons e moléculas polarizadas são atraídos ligando-se a estes componentes de forma reversível. As argilas minerais, as substâncias húmicas e os óxidos de ferro e alumínio possuem determinada superfície de troca e são os principais coloides responsáveis pela capacidade de troca de cátions (CTC) dos solos sob condições tropicais. Em razão do maior número de cargas negativas do que positivas desses coloides, a adsorção é principalmente de cátions. No entanto, há alguns sítios nestes coloides com cargas positivas que podem atrair ânions (principalmente nos óxidos de ferro e alumínio).
A CTC pode ser expressa como “CTC Efetiva” quando considerar todos os cátions permutáveis do solo, expressa sem considerar o íon \(H^{^+}\): \[CTC_{efetiva} = Ca^{2^+} + Mg^{2^+} + K^{^+} + Al^{3^+}\] A soma de bases trocáveis (SB) de um solo, argila ou húmus representa a soma dos teores de cátions permutáveis, exceto \(H^{^+}\): \[SB = Ca^{2^+} + Mg^{2^+} + K^{^+}\]
Denomina-se saturação por bases (V%) a soma das bases trocáveis expressa em porcentagem de capacidade de troca de cátions: \[V_{\%} =\frac {100 \times (Ca^{2^+} + Mg^{2^+} + K^{^+})}{Ca^{2^+} + Mg^{2^+} + K^{^+} + Al^{3^+}} \]
# Criar variáveis adicionais
Solo <- rSolo %>%
mutate(
SB = K + Ca + Mg, # soma de bases
CTC = SB + Al, # capacidade de troca catiônica
V = (SB * 100) / CTC # saturação por bases em porcentagem
)
# Consolidar as variáveis em uma matriz
Solo <- cbind(rSolo, Solo[, c("SB", "CTC", "V")])
# Visualizando os dados
head(Solo)## Area Tipo pH MO N P K Ca Mg Al SB CTC
## 1 C1 Remanescente Florestal 3.8 68 5600 22 2.5 7 8 135 17.5 152.5
## 2 C1 Remanescente Florestal 3.8 59 5320 24 2.6 9 11 149 22.6 171.6
## 3 C1 Remanescente Florestal 3.6 60 2800 16 1.6 9 4 121 14.6 135.6
## 4 C2 Remanescente Florestal 3.4 77 1330 10 2.0 2 1 109 5.0 114.0
## 5 C3 Remanescente Florestal 5.1 68 1792 17 1.6 79 24 22 104.6 126.6
## 6 C3 Remanescente Florestal 4.4 59 3780 14 1.7 29 12 58 42.7 100.7
## V
## 1 11.475410
## 2 13.170163
## 3 10.766962
## 4 4.385965
## 5 82.622433
## 6 42.403178Análise Preliminar da Medida de centralidade: Média Aritmética
As medidas de centralidade são uma forma de resumir um conjunto de dados numéricos, fornecendo informações sobre valores típicos ou “centrais” do conjunto. As principais medidas de centralidade incluem a média, mediana e moda. A média é calculada somando todos os valores do conjunto e dividindo pelo número de observações. Ela é sensível a valores extremos e pode não ser representativa se o conjunto tiver valores muito discrepantes.
O cálculo da média aritmética é dada por: \[\bar{\mu} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}\]
# Calculando a Média Aritmétca
df_resumo <- Solo %>%
group_by(Tipo) %>%
summarise(pH = mean(pH),
MO = mean(MO),
N = mean(N),
P = mean(P),
K = mean(K),
Ca = mean(Ca),
Mg = mean(Mg),
Al = mean(Al),
SB = mean(SB),
CTC = mean(CTC),
V = mean(V))
# Confeccionando uma tabela para apresentação dos resultados
kable(df_resumo,
align = c("l", rep("c", 18)),
caption = "Medidas de centralidade: Média",
booktabs = TRUE,
digits = 2) %>%
kable_styling(full_width = FALSE)| Tipo | pH | MO | N | P | K | Ca | Mg | Al | SB | CTC | V |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Remanescente Florestal | 4.63 | 60.58 | 3112.67 | 15.58 | 1.87 | 39.67 | 15.83 | 64.42 | 57.37 | 121.78 | 49.78 |
| Restauro Florestal Antigo | 4.64 | 50.20 | 2224.60 | 9.90 | 1.75 | 20.10 | 14.05 | 52.00 | 35.90 | 87.89 | 40.70 |
| Pastagem | 4.11 | 26.08 | 1465.92 | 7.33 | 1.35 | 6.17 | 4.08 | 32.00 | 11.60 | 43.60 | 26.10 |
| Restauro Florestal Novo | 4.06 | 37.30 | 1561.00 | 6.40 | 1.28 | 7.20 | 4.60 | 47.00 | 13.08 | 60.08 | 23.99 |
A saturação por bases é um excelente indicativo das condições gerais de fertilidade do solo, sendo utilizada até como complemento na nomenclatura dos solos. Os solos podem ser divididos de acordo com a saturação por bases: solos eutróficos (férteis) = \(V_{\%}\) ≥ 50%; solos distróficos (pouco férteis) = V% < 50%. Alguns solos distróficos podem ser muito pobres em \(Ca^{2^+}\), \(Mg^{2^+}\) e \(K^{^+}\) e apresentar teor de alumínio trocável muito elevado, chegando a apresentar saturação em alumínio (m%) superior a 50% e nesse caso são classificados como solos álicos (muito pobres): Al trocável ≥ 3 \(mmolc.dm^{-3}\) e \(m_{\%}\) ≥ 50%. Neste caso, nenhum solo da região apresentou características eutróficas. Entretanto, é comumente aceito na bibliografia que fragmentos remanescentes (\(V_\%\) = 49.80) devem ter qualidade de solo superior a restauros recentes e de média idade, um padrão observado nos remanescentes florestais da região. Este declínio na qualidade pode estar associado às características fitofisionômicas da região de Cerrado (rico em Alumínio) e pela proximidade com unidades urbanas e agrícolas de alta pressão antrópica. É importante mencionar que, de acordo com a Ecologia de Paisagem, o tamanho do fragmento é um fator determinante para sua qualidade eutrófica.
Pressupostos do Modelo da Análise de Componentes Principais (PCA)
Os principais pressupostos para modelagem de componentes principais estão:
* Multinormalidade: os dados devem apresentar distribuição normal e, para identificar a normalidade das variáveis, podemos utilizar o Teste de Shapiro-Wilk, Teste Kolmogorov-Smirnov, entre outros. Neste estudo, optou-se pelo Teste de Shapiro-Wilk pela baixa quantidade de amostras.
* Colinearidade: a proximidade entre os pontos indica variáveis que estão altamente correlacionadas, tornando o modelo menos confiável.
Teste de Normalidade: Shapiro-Wilk
O teste de Shapiro-Wilk testa a hipótese nula de que uma amostra \(x_{1}, x_{2}, ..., x_{n}\) veio de uma população normalmente distribuída. A estatística do teste é \[W = \frac{(\sum_{i = 1}^{n}a_ix_{(i)})^2}{\sum_{i=1}^{n}(x_i-\overline{x})^2}\] \(x_{(i)}\) é a \(i\)-ésima estatística de ordem, ou seja, \(i\)-ésimo menor número da amostra; \(\overline{x} = \frac{1}{n}(x_{1} + x_{2}, ... + x_{n})\).
Os coeficiente de \(a_i\) são dados por: \[(a_{1}, a_{2}, ..., a_{n}) = \frac{m^TV^{-1}}{C},\] onde \(C\) é uma norma de vetor: \[C = ||V^{1}-m|| = (m^TV^{-1}V^{-1}m)^{\frac{1}{2}}\] e o vetor de \(m\), \[m = (m_1, m_2, ..., m_n)\] é feito dos valores esperados das estatísticas de ordem de variáveis aleatórias independentes e distribuídas de forma idêntica, amostradas a partir da distribuição normal padrão. Finalmente, \(V\) é a matriz de covariância dessas estatísticas de ordem normal.
O teste de hipótese é orientado por:
\(H_0\) para todo \(p-value\) > 0.05, evidencia-se
distribuição normal; e
\(H_1\) para todo \(p-value\) < 0.05, evidencia-se a não
distribuição normal.
shapiro.test(Solo$pH)##
## Shapiro-Wilk normality test
##
## data: Solo$pH
## W = 0.96469, p-value = 0.1118Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de pH são provenientes de distribuições normais com \(p-value > 0.05\).
shapiro.test(Solo$MO)##
## Shapiro-Wilk normality test
##
## data: Solo$MO
## W = 0.9631, p-value = 0.09488Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Matéria Orgânica são provenientes de distribuições normais com \(p-value > 0.05\).
shapiro.test(Solo$N)##
## Shapiro-Wilk normality test
##
## data: Solo$N
## W = 0.88473, p-value = 8.793e-05Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Nitrogênio (N) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$P)##
## Shapiro-Wilk normality test
##
## data: Solo$P
## W = 0.91957, p-value = 0.001439Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Fósforo (P) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$K)##
## Shapiro-Wilk normality test
##
## data: Solo$K
## W = 0.91391, p-value = 0.0008859Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Potássio (K) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$Ca)##
## Shapiro-Wilk normality test
##
## data: Solo$Ca
## W = 0.78616, p-value = 1.935e-07Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Cálcio (Ca) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$Mg)##
## Shapiro-Wilk normality test
##
## data: Solo$Mg
## W = 0.9139, p-value = 0.0008849Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Magnésio (Mg) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$SB)##
## Shapiro-Wilk normality test
##
## data: Solo$SB
## W = 0.84702, p-value = 6.571e-06Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Soma de Bases (SB) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$CTC)##
## Shapiro-Wilk normality test
##
## data: Solo$CTC
## W = 0.95148, p-value = 0.02888Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Capacidade de Troca de Cátions (CTC) são provenientes de distribuições não-normais com \(p-value < 0.05\).
shapiro.test(Solo$V)##
## Shapiro-Wilk normality test
##
## data: Solo$V
## W = 0.95605, p-value = 0.0459Ao nível de significância \(p-value < 0.05\), concluímos que as amostras de Saturação de Bases (V) são provenientes de distribuições não-normais com \(p-value < 0.05\).
BoxCox: Normalizando Dados
Aplicamos a transformação de BoxCox nas variáveis que obtiveram \(p-value < 0.05\). As transformações de dados designadas para atingir uma finalidade especificada, por exemplo, estabilidade de variância, aditividade de efeitos e simetria da densidade. Se alguém for bem sucedido em encontrar uma transformação adequada, o método ordinário de análise estará disponível. Entre as muitas transformações paramétricas, a família em [a1] é comumente utilizada.
\[ X ^ {( \lambda ) } = \left \{ \begin{array}{l} { { \frac{X ^ \lambda - 1 } \lambda } \ \textrm{ for } \lambda \neq0, } \\ { { \mathop{\rm log} } X \ \textrm{ for } \lambda = 0. } \end{array} \right . \] A fórmula \(\frac{(x^λ−1)}{λ}\) é escolhido de modo que \(x(λ)\) é contínuo como \(λ\) tende a zero e monótono aumentando em relação a \(x\) para qualquer \(λ\). O parâmetro power \(λ\) é estimado por uma técnica gráfica ou pelo método de máxima verossimilhança. Infelizmente, um formulário fechado para o estimador \(\hat{λ}\) raramente podem ser encontrados. Assim, o gráfico da máxima verossimilhança contra \(λ\) é útil. O valor de \(\hat{λ}\) obtido dessa maneira é tratado como se fosse um valor verdadeiro e, em seguida, ajusta-se o modelo aos dados transformados. Tal abordagem pode ser facilmente realizada, e uma teoria assintótica associada a outros parâmetros é útil.
Este tratamento tem, no entanto, algumas dificuldades porque \(\hat{λ}\) tem uma variabilidade e depende dos próprios dados fornecidos. Sabe-se que a estimativa de \(λ\) pela máxima verossimilhança e pelos testes de razão de verossimilhança conexos podem ser fortemente influenciados por valores atípicos (Outlier). Além disso, em certas situações, a teoria limitante usual baseada no conhecimento \(λ\) não se sustenta no caso desconhecido. Por conseguinte, foram propostos vários procedimentos de estimativa sólidos.
# Pacotes de dados necessários
library(car)
library(fpp)
# Preparando a função que recepcionará os dados
transform <- function(x){
x <- BoxCox(x,
lambda = BoxCox.lambda(
x, method = c('loglik'), lower = -3, upper = 3))
}
# Normalizando os dados com base na função criada e adaptada
Normalizado <- cbind(Solo[, 1:4], # Código para substituição dos valores nas colunas
lapply(Solo[, 5:13], transform)) # Código de transformação
# Visualizando os dados
head(Normalizado, 15)## Area Tipo pH MO N P K Ca
## 1 C1 Remanescente Florestal 3.8 68 5.781275 5.571932 0.7197571 2.1481404
## 2 C1 Remanescente Florestal 3.8 59 5.759580 5.832579 0.7431979 2.4573094
## 3 C1 Remanescente Florestal 3.6 60 5.478483 4.682902 0.4141682 2.4573094
## 4 C2 Remanescente Florestal 3.4 77 5.129037 3.539203 0.5763270 0.7177346
## 5 C3 Remanescente Florestal 5.1 68 5.272123 4.844601 0.4141682 5.4797060
## 6 C3 Remanescente Florestal 4.4 59 5.612160 4.338619 0.4602111 4.0036033
## 7 C3 Remanescente Florestal 5.0 56 4.806141 4.514493 0.5407938 4.9826104
## 8 C3 Remanescente Florestal 4.7 70 5.521373 3.756175 0.4141682 4.2694359
## 9 C4 Remanescente Florestal 6.1 56 5.531426 4.844601 0.5763270 5.6116873
## 10 C4 Remanescente Florestal 5.8 55 4.806141 4.154376 0.4141682 5.3367469
## 11 C4 Remanescente Florestal 5.0 45 5.724564 4.154376 0.4602111 4.3872689
## 12 C4 Remanescente Florestal 4.9 54 5.658147 4.514493 0.4141682 4.8739935
## 13 R6 Restauro Florestal Antigo 3.7 32 5.461168 2.492009 0.6092085 1.4869835
## 14 R6 Restauro Florestal Antigo 3.8 51 5.455231 3.058657 0.8423199 2.1481404
## 15 R1 Restauro Florestal Antigo 4.6 47 5.426767 3.756175 0.3071734 3.1101942
## Mg Al SB CTC V
## 1 2.886887 2.568127 3.313854 11.727871 4.440983
## 2 3.510455 2.590446 3.658746 12.270594 4.867080
## 3 1.719055 2.542576 3.074807 11.206407 4.252623
## 4 0.000000 2.517407 1.746189 10.468927 2.100097
## 5 5.315193 2.014615 5.920370 10.909909 13.976531
## 6 3.691453 2.347394 4.556025 9.964705 9.775996
## 7 5.092358 2.176062 5.516539 10.501526 12.851493
## 8 4.324656 2.347394 4.862730 10.344121 10.411789
## 9 6.433854 1.854050 6.193554 11.314830 14.545365
## 10 4.729818 1.882416 5.727551 10.269447 14.207832
## 11 4.465209 2.064231 4.974706 9.156176 12.753919
## 12 4.854854 2.143544 5.391250 10.155610 12.874101
## 13 1.301297 2.376085 2.471082 8.746306 4.689612
## 14 2.372566 2.405433 3.203303 9.398302 6.035323
## 15 3.317541 2.347394 3.872681 9.278595 8.240752Colinearidade e Dimensionalidade dos dados
Dimensionalidade e problemas relacionados
A dimensionalidade dos dados é definida pelo número de variáveis independentes, preditivas ou explicativas no conjunto de dados. Dito isto, muitas vezes é observado que a dimensionalidade afeta o desempenho computacional mais do que o número de observações. Não apenas aumenta a complexidade do conjunto de dados, mas também limita a capacidade de explorar e modelar relacionamentos entre variáveis.
Causas da dimensionalidade
Dois fatores podem contribuir para a alta dimensão. Primeiro, ter variáveis categóricas com vários níveis e, segundo, variáveis redundantes em um conjunto de dados. Quando incluímos variáveis categóricas na modelagem, usamos variáveis fictícias, e quanto mais o número de níveis em uma variável categórica, mais variáveis fictícias criamos, o que aumenta a dimensionalidade. Isso pode ser resolvido por uma abordagem conhecida como “Collapsing the levels”, que será discutida em um documento separado. Este documento se concentrará em lidar com variáveis redundantes e multicolinearidade.
A multicolinearidade (colinearidade) pode ser definida como um fenômeno no qual duas ou mais variáveis em uma regressão multivariada são altamente correlacionadas, ou seja, uma variável pode ser prevista linearmente a partir das outras com um grau substancial de precisão. As variáveis colineares contêm as mesmas informações sobre a variável dependente. Se medidas nominalmente “diferentes” realmente quantificam o mesmo fenômeno, então elas são variáveis redundantes.
A presença de variáveis redundantes e multicolinearidade entre as variáveis pode afetar as análises de variáveis de várias maneiras, algumas das quais são compartilhadas abaixo:
- Desestabilizar as estimativas dos parâmetros
- Interpretações confusas do modelo
- Aumento do risco de overfitting
- Aumento do tempo de computação
Pequenas mudanças nos dados de entrada (como remover ou adicionar uma única variável) levam a uma grande mudança no modelo, resultando até mesmo em mudanças de sinal dos parâmetros. Além disso, a presença de multicolinearidade aumenta a variância ou erro padrão das estimativas dos coeficientes, tornando-os sensíveis a pequenas mudanças, resultando em dificuldade de interpretação. Quando falamos sobre o melhor modelo de regressão para um conjunto de dados, isso significa que o modelo considera variáveis preditoras em que cada preditor se correlaciona com a variável alvo/dependente, mas se correlaciona quase minimamente entre si. Tal modelo é chamado de modelo de baixo ruído e é estaticamente robusto e pode ser usado para pontuação Agora, se a correlação entre as variáveis em si for alta e atingirmos alto poder preditivo no conjunto de dados de treinamento, é possível que não alcancemos o mesmo no conjunto de dados de teste, a menos que exista a mesma relação colinear entre as variáveis no conjunto de dados de teste, caso em que a multicolinearidade não afetará a capacidade preditiva de seu modelo.
Detectando Multicolinearidade
Depois de entender como a multicolinearidade afeta nossa análise e capacidade preditiva do modelo, é muito importante aprender a detectar a presença de multicolinearidade em nossos dados. Os seguintes métodos podem ser usados para detectar multicolinearidade: * A análise exibe sinal de multicolinearidade - como estimativas de coeficientes variando de modelo para modelo * O valor do quadrado R é grande, mas nenhum dos pesos beta é estatisticamente significativo, ou seja, o teste F para o modelo geral é significativo, mas os testes t para estimativas de coeficientes individuais não são. * Correlação entre pares de variáveis são grandes. * Fator de inflação de variância.
Correlação de Pearson
Utilizaremos neste estudo a Correlação de Pearson para determinar se há colinearidade entre os dados, aceitando um coeficiente de correlação > .75 como critério de exclusão. Em estatística descritiva, o coeficiente de correlação de Pearson, também chamado de “coeficiente de correlação produto-momento” ou simplesmente de ρ de Pearson mede o grau da correlação (e a direção dessa correlação - se positiva ou negativa) entre duas variáveis de escala métrica (intervalar ou de rácio/razão). Calcula-se o coeficiente de correlação de Pearson segundo a seguinte fórmula:
\[ R_{x, y} = \frac{(\sum^{n}_{i=1} x_{i} - \overline{x}) – (\sum^{n}_{i=1} y_{i} - \overline{y}) }{\sqrt{\sum_{i=1}^{n}(x_{i}-\overline{x})^2 \times {\sqrt{\sum_{i=1}^{n}(y_{i}-\overline{y})^2}}}}\] \(x_{1}, x_{2}, ..., x_{n}\) e \(y_{1}, y_{2}, ..., y_{n}\) são os valores medidos de ambas as variáveis, e \(\overline{x} = \frac{1}{n}\times \sum_{i = 1}^{n}x_{i}\) e \(\overline{y} = \frac{1}{n}\times \sum_{i = 1}^{n}y_{i}\) são as médias aritméticas de ambas as variáveis.
# Armazenando os valores de correlação entre variáveis em um objeto Pearson
Pearson <- cor(Normalizado[3:13], method = c('pearson'))
# Visualizando os dados
head(round(Pearson,2))## pH MO N P K Ca Mg Al SB CTC V
## pH 1.00 0.31 -0.01 0.31 0.11 0.90 0.89 -0.60 0.91 0.29 0.95
## MO 0.31 1.00 0.35 0.57 0.52 0.53 0.54 0.38 0.56 0.90 0.25
## N -0.01 0.35 1.00 0.35 0.36 0.18 0.23 0.34 0.21 0.44 0.01
## P 0.31 0.57 0.35 1.00 0.48 0.49 0.49 0.08 0.52 0.66 0.33
## K 0.11 0.52 0.36 0.48 1.00 0.23 0.32 0.24 0.29 0.57 0.09
## Ca 0.90 0.53 0.18 0.49 0.23 1.00 0.94 -0.36 0.99 0.55 0.92# Carregando o pacote para análise gráfica da correlação de Pearson
library(corrplot)
# Visualizando a matriz de correlação de Pearson por análise gráfica
corrplot(Pearson, method="circle")
# Criando um novo objeto com as dimensões reduzidas
Solo_Dim <- Normalizado[, c(1, 2, 4, 5, 6, 7, 12, 13)]
# Visualizando o objeto
head(Solo_Dim)## Area Tipo MO N P K CTC
## 1 C1 Remanescente Florestal 68 5.781275 5.571932 0.7197571 11.727871
## 2 C1 Remanescente Florestal 59 5.759580 5.832579 0.7431979 12.270594
## 3 C1 Remanescente Florestal 60 5.478483 4.682902 0.4141682 11.206407
## 4 C2 Remanescente Florestal 77 5.129037 3.539203 0.5763270 10.468927
## 5 C3 Remanescente Florestal 68 5.272123 4.844601 0.4141682 10.909909
## 6 C3 Remanescente Florestal 59 5.612160 4.338619 0.4602111 9.964705
## V
## 1 4.440983
## 2 4.867080
## 3 4.252623
## 4 2.100097
## 5 13.976531
## 6 9.775996Análise de Componentes Principais
Muitas vezes, não é útil ou informativo apenas olhar para todas as variáveis em um conjunto de dados para correlações ou covariâncias. Uma abordagem preferível é derivar novas variáveis das variáveis originais que preservam a maior parte das informações fornecidas por suas variâncias. A análise de componentes principais é um método estatístico amplamente utilizado e popular para reduzir dados com muitas dimensões (variáveis), projetando os dados com menos dimensões usando combinações lineares das variáveis, conhecidas como componentes principais. As novas variáveis projetadas (componentes principais) não são correlacionadas entre si e são ordenadas de forma que os primeiros componentes retenham a maior parte da variação presente nas variáveis originais. Assim, o PCA também é útil em situações em que as variáveis independentes estão correlacionadas entre si e pode ser empregado na análise exploratória de dados ou na construção de modelos preditivos. A análise de componentes principais também pode revelar características importantes dos dados, como outliers e desvios de uma distribuição multinormal.
Definindo Componentes Principais
O primeiro passo na definição dos componentes principais de \(p\) variáveis originais é encontrar uma função linear \(a′_1y\) , onde \(a_1\) é um vetor de \(p\) constantes, para os vetores de observação que possuem variância máxima. Esta função linear é definida como: \[a′_1y = a_{11} x_1 + a_{12} x_2 + ⋯ + a_{1p} x_p = \sum_{j=1}^{p} a_{1j}x_j\] A análise de componentes principais continua a encontrar uma função linear \(a'_{2}y\) que não está correlacionada com \(a'_{1}y\) com variância maximizada e assim por diante até \(k\) componentes principais.
Derivação de Componentes Principais
Os componentes principais de um conjunto de dados são obtidos da matriz de covariância de amostra \(S\) ou da matriz de correlação \(R\). Embora os componentes principais obtidos de \(S\) sejam o método original de análise de componentes principais, os componentes de \(R\) podem ser mais interpretáveis se as variáveis originais tiverem unidades diferentes ou amplas variações (Rencher 2002, pp. 393). Por enquanto, \(S\) será referido como \(\sum\) (denota uma matriz de covariância conhecida) que será usada na derivação. O objetivo da derivação é encontrar \(a′_{k}y\) que maximize a variância de \(a′_{k}y \sum a_{k}\). Para isso, vamos considerar o primeiro vetor \(a′_{1}y\) que maximiza \(Var(a′_{1}y) = a′_{1}y \sum a_{1}\). Para fazer essa maximização, precisaremos de uma restrição para controlar valores desnecessariamente grandes de \(a_{1}\). A restrição neste exemplo é o vetor de comprimento unitário \(a′_{1}a_{1} = 1\). Essa restrição é empregada com um multiplicador de Lagrange λ para que a função seja maximizada em uma restrição de igualdade de \(g(x) = 0\). Assim, a função lagrangiana é definida como: \[a′_{1} \sum a_{1} – λ(a′_{1}ya_{1} – 1)\]
Breve Aparte: Multiplicadores de Lagrange
O método do multiplicador de Lagrange é usado para encontrar um máximo ou mínimo de uma função multivariada com alguma restrição nos valores de entrada. Como estamos interessados em maximização, o problema pode ser resumido como ‘maximize \(f(x)\) sujeito a \(g(x) = c\)’. Neste exemplo, \(g(x) = a′_{1}ya_{1} = 1\) e \(f(x) = a′_{1}y \sum a_{1}\). O multiplicador de Lagrange, definido como \(λ\), permite a combinação de f(x) e g(x) em uma nova função \(L(x, λ)\), definida como: \[L(a′_{1}y, λ) = f(a′_{1}y) − λ(g(a′_{1}y) − c)\]
O sinal de \(λ\) pode ser positivo ou negativo. A nova função é então resolvida para um ponto estacionário, neste caso 0 , usando derivadas parciais: \[\frac{∂L(a′_{1}y, λ)}{∂a′_{1}y} = 0\]
Voltando à análise de componentes principais, diferenciamos \(L(a_{1}) = a′_{1} \sum a_{1} − λ(a′_{1}ya_{1} − 1)\) em relação a \(a_{1}\): \[\frac {∂L}{∂a_{1}} = 2 \sum a_{1} − 2λa_{1} = 0\] \[\sum a_{1} – λa_{1} = 0\]
Expressando o que foi dito acima com uma matriz identidade, I: \[(\sum − λI)a_{1} = 0\]
O que mostra que \(λ\) é um autovetor da matriz de covariância \(\sum\) e \(a_{1}\) é o autovetor correspondente. Como afirmado anteriormente, estamos interessados em encontrar \(a′_{1}y\) com variância máxima. Portanto \(λ\) deve ser o maior possível que segue \(a_{1}\) é o autovetor correspondente ao maior autovalor de \(\sum\). Os componentes principais restantes são encontrados de maneira semelhante e correspondem ao \(k\)-ésimo componente principal. Assim, a segunda componente principal é \(a′_{2}y\) e é equivalente ao autovetor do segundo maior autovalor de \(\sum\), e assim por diante.
# Carregando pacotes
library(FactoMineR)
library(factoextra)
# Gerando a matriz de autovalores e autovetores
Solo.PCA <- PCA(Solo_Dim[, 3:8], graph = F)
# Visualizando a matriz de autovalores e autovetores
summary(Solo.PCA)##
## Call:
## PCA(X = Solo_Dim[, 3:8], graph = F)
##
##
## Eigenvalues
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## Variance 3.216 1.041 0.678 0.541 0.439 0.085
## % of var. 53.595 17.348 11.307 9.009 7.318 1.423
## Cumulative % of var. 53.595 70.943 82.250 91.259 98.577 100.000
##
## Individuals (the 10 first)
## Dist Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr
## 1 | 4.447 | 3.843 8.506 0.747 | -2.110 7.921 0.225 | 0.007 0.000
## 2 | 4.510 | 3.869 8.620 0.736 | -1.980 6.976 0.193 | 0.082 0.018
## 3 | 2.745 | 1.959 2.210 0.509 | -1.339 3.190 0.238 | -0.518 0.733
## 4 | 3.477 | 1.581 1.439 0.207 | -1.512 4.069 0.189 | -2.496 17.003
## 5 | 3.094 | 2.505 3.613 0.656 | 1.711 5.210 0.306 | -0.171 0.080
## 6 | 2.159 | 1.963 2.219 0.826 | -0.218 0.085 0.010 | 0.771 1.621
## 7 | 3.162 | 1.381 1.098 0.191 | 2.119 7.986 0.449 | -1.775 8.597
## 8 | 2.407 | 2.067 2.460 0.738 | 0.121 0.026 0.003 | 0.261 0.186
## 9 | 3.331 | 2.957 5.035 0.788 | 1.109 2.189 0.111 | 0.826 1.863
## 10 | 3.111 | 0.935 0.504 0.090 | 2.577 11.815 0.686 | -1.424 5.536
## cos2
## 1 0.000 |
## 2 0.000 |
## 3 0.036 |
## 4 0.515 |
## 5 0.003 |
## 6 0.127 |
## 7 0.315 |
## 8 0.012 |
## 9 0.061 |
## 10 0.210 |
##
## Variables
## Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr cos2
## MO | 0.876 23.881 0.768 | 0.041 0.158 0.002 | -0.251 9.310 0.063 |
## N | 0.563 9.861 0.317 | -0.472 21.445 0.223 | 0.669 66.030 0.448 |
## P | 0.797 19.733 0.635 | 0.173 2.864 0.030 | 0.038 0.212 0.001 |
## K | 0.724 16.279 0.523 | -0.230 5.063 0.053 | -0.176 4.568 0.031 |
## CTC | 0.925 26.628 0.856 | -0.035 0.121 0.001 | -0.182 4.876 0.033 |
## V | 0.341 3.617 0.116 | 0.856 70.349 0.732 | 0.319 15.004 0.102 |# Guardando os autovalores em um objeto
Solo_Eig <- get_eigenvalue(Solo.PCA)
# Visualizando os autovalores
Solo_Eig## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.21568300 53.594717 53.59472
## Dim.2 1.04089109 17.348185 70.94290
## Dim.3 0.67843349 11.307225 82.25013
## Dim.4 0.54051879 9.008647 91.25877
## Dim.5 0.43908709 7.318118 98.57689
## Dim.6 0.08538654 1.423109 100.00000Os dois primeiros componentes principais respondem por \(70.94\%\) da variância total. Um gráfico de plotagem dos autovalores pode ser traçado para visualizar a proporção de variância explicada por cada autovalor subsequente.
#fviz_eig(Solo.PCA, addlabels = T, ylim = c(0, 90), main = 'Gráfico de Autovalores')
fviz_eig(Solo.PCA,
geom = 'line',
width=0.8,
addlabels=T,
linecolor = 'black',
ylim = c(0, 60),
main = 'Plotagem de Autovalores')
Os elementos dos autovetores de S são os ‘coeficientes’ ou ‘cargas’ dos componentes principais.
# Criando um objeto para extração dos autovetores para variáveis
Solo.Var <- get_pca_var(Solo.PCA)
# Criando um objeto para extração dos autovetores para amostras
Solo.Ind <- get_pca_ind(Solo.PCA)
# Extraindo os autovetores: 'coord'
Solo.Var$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## MO 0.8763261 0.04061008 -0.25132436 -0.31194128 0.18774509
## N 0.5631184 -0.47245948 0.66930534 -0.07909167 0.07193295
## P 0.7965893 0.17265936 0.03791857 0.08650356 -0.57028984
## K 0.7235232 -0.22956742 -0.17603949 0.60092487 0.17790287
## CTC 0.9253556 -0.03546154 -0.18187842 -0.23946493 0.05802222
## V 0.3410424 0.85571840 0.31904928 0.10497107 0.19600513Os dois primeiros componentes principais são assim:
\[z_{1} = a′_{1}y = .876y_{1} + .563y_{2} + .797y_{3} + .724y_{4} + .925y_{5} + .341y_{6}\]
\[z_{2} = a′_{2}y = .041y_{1} - .472y_{2} + .173y_{3} - .230y_{4} - .035y_{5} + .856y_{6}\]
Círculo de correlação
A correlação entre uma variável e um componente principal (CP) é usada como as coordenadas da variável no PC. A representação das variáveis difere do gráfico das observações: as observações são representadas por suas projeções, mas as variáveis são representadas por suas correlações (Abdi e Williams 2010).
# Gerando gráfico de autovetores para as dimensões 1 e 2
fviz_pca_var(Solo.PCA, col.var = 'blue')
\[ \]
Nesse caso, é possível identificar que a MO, CTC, K e P apresentam
contribuições similares para o Eixo 1, que pode ser identificado pela
distância angular dos vetores dessa variável com esse eixo. Além disso,
é possível identificar uma forte correlação entre as variáveis de CTC,
MO, P e K demonstrada pela formação de ângulos agudos entre estes
vetores. Por sua vez, V% apresenta maior correlação com o Eixo 2. Além
disso, é possível inferir que as variáveis V% e N apresentam correlação
próxima da nulidade ou muito fraca, que graficamente é representada pela
formação de ângulo próximo a 90º.
Ainda é possível observar que o solo das unidades de Remanescentes Florestais apresenta forte correlação com MO, P, K e CTC que pode ser identificada pela orientação e comprimento desses vetores e pela tabela de resumo dos dados, considerando também a posição desses pontos em relação ao eixo das abcissas e ordenadas. Logo, é possível inferir que o solo dos remanescentes florestais apresenta maior fertilidade. Conforme Genu et al. (2013) e Begon et al. (2007), áreas com maiores índices (frações) de argila possibilitam maior retenção de cátions, tais como (Ca2+), magnésio (Mg2+), potássio (K+), sódio (Na+), alumínio (Al3+) e hidrogênio (H+) e, consequentemente, aumentando a fertilidade do solo (CTC). Brady-Weil (2013) descreve que a camada superficial do solo é responsável pela adsorção de 50 a 90% dos cátions que viabilizam índices mais elevados de Ca, Mg e K, que estão diretamente associadas com o CTC. Ou seja, trata-se de uma área homogênea que passou por diferentes estágios sucessionais, e, hoje, abriga uma fitofisionomia clímax. Essa correlação entre variáveis pode ser identificada pela matriz de correlação de Pearson, que define a colinearidade.
Biplot entre autovetores e unidades de restauro amostradas (Correspondência Canônica)
fviz_pca_biplot(Solo.PCA, habillage = as.factor(Normalizado$Tipo), label = 'var',
addEllipses=TRUE, ellipse.level=0.95, pointsize = 2, title = 'PCA de Restauro Florestal')
O Restauro Florestal com 20 anos (Restauro Antigo), apresenta um comportamento de solo muito semelhante às unidades de controle, porém com menores concentrações de MO, P, K e CTC. Portanto, é possível identificar uma menor correlação linear que graficamente é representada pela maior proximidade e dispersão dos dados em relação ao centro biplot. Essa unidade de estudo apresentou maior correlação com o N, que é uma variável independente de V%. Portanto, essa interação indica que o solo está passando por um processo de fixação de Nitrogênio e, consequentemente, apresenta menor fertilidade que a unidade de controle.
Conclusão
Apenas depois de 20 anos de restauro florestal é possível identificar qualidade de solo semelhante às vegetações remanescentes.