Análise de Solo e Climática
2026-01-18
Carregamento de Bibliotecas
library(tidyverse)
library(readODS)
library(car)
library(rstatix)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(FSA)
library(vegan)
library(ggplot2)Importação e Preparação de Dados
# Leitura dos dados originais
solo <- read.table("dados_quimicos.csv", sep = ";", dec = ",", header = TRUE)
clima <- read.table("dados_climáticos.csv", sep = ";", dec = ",", header = TRUE) Analise dos dados de Solo
# Limpeza da variável 'area' e remoção de coluna extra
solo_f <- solo %>%
filter(area != "Mangueira")
solo_f$X = NULLAjustar nomes das variáveis:
colnames(solo_f) = c("area" ,"pH", "Ca", "Mg", "Al", "H.Al", "CTC","P","K","Mat_Org","Sat.Al",
"Sat.Base","Corg","RBS","CBM","qMic", "PT","Macro.por","Micro.por","Densi",
"Argila","Areia","Silte")Separar grupos de variáveis
# Variáveis biológicas
bio_vars <- c("RBS","CBM","qMic")
# Variáveis químicas
chem_vars <- c( "pH", "Ca", "Mg", "Al", "H.Al", "CTC","P","K","Mat_Org","Sat.Al", "Sat.Base","Corg")
# Variáveis físicas
phys_vars <- c("PT","Macro.por","Micro.por","Densi", "Argila","Areia","Silte")Estatística descritiva por área
solo_f %>%
group_by(area) %>%
summarise(across(all_of(c(bio_vars, chem_vars, phys_vars)),
list(mediana = median,
media = mean,
sd = sd),
na.rm = TRUE))## # A tibble: 2 × 67
## area RBS_mediana RBS_media RBS_sd CBM_mediana CBM_media CBM_sd qMic_mediana
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Degrad… 116. 119. 43.8 92.7 123. 93.5 0.5
## 2 Flores… 97.9 121. 73.2 70.9 179. 228. 0.4
## # ℹ 59 more variables: qMic_media <dbl>, qMic_sd <dbl>, pH_mediana <dbl>,
## # pH_media <dbl>, pH_sd <dbl>, Ca_mediana <dbl>, Ca_media <dbl>, Ca_sd <dbl>,
## # Mg_mediana <dbl>, Mg_media <dbl>, Mg_sd <dbl>, Al_mediana <dbl>,
## # Al_media <dbl>, Al_sd <dbl>, H.Al_mediana <dbl>, H.Al_media <dbl>,
## # H.Al_sd <dbl>, CTC_mediana <dbl>, CTC_media <dbl>, CTC_sd <dbl>,
## # P_mediana <dbl>, P_media <dbl>, P_sd <dbl>, K_mediana <dbl>, K_media <dbl>,
## # K_sd <dbl>, Mat_Org_mediana <dbl>, Mat_Org_media <dbl>, Mat_Org_sd <dbl>, …Graficos complementares - BoxPlot
Atributos Biológicos
solo_f %>%
mutate(area = factor(area, levels = c("Degradada", "Floresta"))) %>%
pivot_longer(
cols = all_of(bio_vars),
names_to = "variavel",
values_to = "valor"
) %>%
ggplot(aes(x = area, y = valor, fill = area)) +
geom_violin(trim = TRUE, scale = "width", color = NA, alpha = 0.6) +
geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.9) +
geom_jitter(width = 0.08, size = 1.5, alpha = 0.6) +
facet_wrap(~ variavel, scales = "free_y") +
scale_fill_manual(
values = c(
"Degradada" = "#9E9E9E",
"Floresta" = "#1B7837"
),
drop = FALSE
) +
guides(fill = guide_legend(override.aes = list(alpha = 1))) +
theme_classic(base_size = 12) +
theme(
legend.position = "top",
legend.title = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
strip.background = element_blank(),
strip.text = element_text(face = "bold")
)
Atributos Químicos
solo_f %>%
mutate(area = factor(area, levels = c("Degradada", "Floresta"))) %>%
pivot_longer(
cols = all_of(chem_vars),
names_to = "variavel",
values_to = "valor"
) %>%
ggplot(aes(x = area, y = valor, fill = area)) +
geom_violin(trim = TRUE, scale = "width", color = NA, alpha = 0.6) +
geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.9) +
geom_jitter(width = 0.08, size = 1.5, alpha = 0.6) +
facet_wrap(~ variavel, scales = "free_y") +
scale_fill_manual(
values = c(
"Degradada" = "#9E9E9E",
"Floresta" = "#1B7837"
),
drop = FALSE
) +
guides(fill = guide_legend(override.aes = list(alpha = 1))) +
theme_classic(base_size = 12) +
theme(
legend.position = "top",
legend.title = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
strip.background = element_blank(),
strip.text = element_text(face = "bold")
)
Atributos Físicos
solo_f %>%
mutate(area = factor(area, levels = c("Degradada", "Floresta"))) %>%
pivot_longer(
cols = all_of(phys_vars),
names_to = "variavel",
values_to = "valor"
) %>%
ggplot(aes(x = area, y = valor, fill = area)) +
geom_violin(trim = TRUE, scale = "width", color = NA, alpha = 0.6) +
geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.9) +
geom_jitter(width = 0.08, size = 1.5, alpha = 0.6) +
facet_wrap(~ variavel, scales = "free_y") +
scale_fill_manual(
values = c(
"Degradada" = "#9E9E9E",
"Floresta" = "#1B7837"
),
drop = FALSE
) +
guides(fill = guide_legend(override.aes = list(alpha = 1))) +
theme_classic(base_size = 12) +
theme(
legend.position = "top",
legend.title = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
strip.background = element_blank(),
strip.text = element_text(face = "bold")
)
Teste de normalidade (variaveis biológicas)
solo_f %>%
pivot_longer(cols = all_of(bio_vars),
names_to = "variavel",
values_to = "valor") %>%
group_by(area, variavel) %>%
shapiro_test(valor)## # A tibble: 6 × 5
## area variavel variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Degradada CBM valor 0.909 0.311
## 2 Degradada RBS valor 0.953 0.727
## 3 Degradada qMic valor 0.868 0.116
## 4 Floresta CBM valor 0.843 0.0625
## 5 Floresta RBS valor 0.687 0.000993
## 6 Floresta qMic valor 0.856 0.0863Wilcoxon-Mann-Whitney (variáveis biologicas)
kw_bio <- solo_f %>%
pivot_longer(cols = all_of(bio_vars),
names_to = "variavel",
values_to = "valor") %>%
group_by(variavel) %>%
wilcox_test(valor ~ area) %>%
add_significance()
print(kw_bio)## # A tibble: 3 × 9
## variavel .y. group1 group2 n1 n2 statistic p p.signif
## <chr> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
## 1 CBM valor Degradada Floresta 9 9 43 0.86 ns
## 2 RBS valor Degradada Floresta 9 9 48.5 0.508 ns
## 3 qMic valor Degradada Floresta 9 9 45 0.723 nsO teste de Wilcoxon-Mann-Whitney não indicou diferenças significativas entre Floresta e Solo Degradado para as variáveis biológicas analisadas (p > 0,05).
Os atributos biológicos do solo (respiração basal do solo – RBS, carbono da biomassa microbiana – CBM e quociente microbiano – qMic) apresentaram valores médios e medianos relativamente próximos entre as áreas, com elevada variabilidade interna, especialmente na área de Floresta. Essa maior dispersão sugere um ambiente biologicamente mais heterogêneo na floresta, possivelmente associado à maior complexidade estrutural e diversidade de microambientes. Apesar disso, as análises inferenciais não detectaram diferenças estatisticamente significativas entre as áreas, indicando que, isoladamente, esses indicadores não discriminam fortemente os ambientes sob o desenho amostral adotado.
Teste de normalidade (variaveis quimicas)
solo_f %>%
pivot_longer(cols = all_of(chem_vars),
names_to = "variavel",
values_to = "valor") %>%
group_by(area, variavel) %>%
shapiro_test(valor)## # A tibble: 24 × 5
## area variavel variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Degradada Al valor 0.638 0.000266
## 2 Degradada CTC valor 0.867 0.114
## 3 Degradada Ca valor 0.846 0.0679
## 4 Degradada Corg valor 0.836 0.0527
## 5 Degradada H.Al valor 0.870 0.124
## 6 Degradada K valor 0.854 0.0834
## 7 Degradada Mat_Org valor 0.942 0.598
## 8 Degradada Mg valor 0.617 0.000153
## 9 Degradada P valor 0.794 0.0173
## 10 Degradada Sat.Al valor 0.919 0.381
## # ℹ 14 more rowsWilcoxon-Mann-Whitneys (variáveis quimicas)
kw_quim <- solo_f %>%
pivot_longer(cols = all_of(chem_vars),
names_to = "variavel",
values_to = "valor") %>%
group_by(variavel) %>%
wilcox_test(valor ~ area) %>%
add_significance()
print(kw_quim)## # A tibble: 12 × 9
## variavel .y. group1 group2 n1 n2 statistic p p.signif
## <chr> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
## 1 Al valor Degradada Floresta 9 9 30.5 0.336 ns
## 2 CTC valor Degradada Floresta 9 9 17.5 0.046 *
## 3 Ca valor Degradada Floresta 9 9 30 0.365 ns
## 4 Corg valor Degradada Floresta 9 9 5.5 0.00223 **
## 5 H.Al valor Degradada Floresta 9 9 16 0.0327 *
## 6 K valor Degradada Floresta 9 9 55.5 0.178 ns
## 7 Mat_Org valor Degradada Floresta 9 9 20.5 0.0819 ns
## 8 Mg valor Degradada Floresta 9 9 37.5 0.797 ns
## 9 P valor Degradada Floresta 9 9 63.5 0.0443 *
## 10 Sat.Al valor Degradada Floresta 9 9 47 0.595 ns
## 11 Sat.Base valor Degradada Floresta 9 9 36.5 0.756 ns
## 12 pH valor Degradada Floresta 9 9 35 0.639 nsOs atributos químicos mostraram um padrão mais contrastante entre as áreas. Variáveis relacionadas à fertilidade e ao estado químico do solo, como carbono orgânico (Corg), capacidade de troca catiônica (CTC), hidrogênio + alumínio (H+Al) e fósforo disponível (P), apresentaram diferenças consistentes entre Floresta e Solo Degradado. Em geral, a área de Floresta apresentou maiores teores de carbono orgânico e melhores condições químicas, enquanto o Solo Degradado mostrou sinais de empobrecimento químico.
Teste de normalidade (variaveis físicas)
solo_f %>% shapiro_test(Silte)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Silte 0.642 0.0000183solo_f %>% shapiro_test(Areia)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Areia 0.642 0.0000183solo_f %>% shapiro_test(Argila)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Argila 0.642 0.0000183solo_f %>% shapiro_test(Densi)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Densi 0.840 0.00590solo_f %>% shapiro_test(Micro.por)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Micro.por 0.566 0.00000321solo_f %>% shapiro_test(Macro.por)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 Macro.por 0.638 0.0000164solo_f %>% shapiro_test(PT)## # A tibble: 1 × 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 PT 0.808 0.00199Variáveis como densidade do solo, porosidade total e frações texturais apresentaram padrões distintos
PCA do solo
solo_pca <- solo_f %>%
select(all_of(c(bio_vars, chem_vars, phys_vars))) %>%
scale()pca <- prcomp(solo_pca, center = TRUE, scale. = TRUE)library(factoextra)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(tibble)
# Coordenadas de TODAS as variáveis
vars_coord_all <- get_pca_var(pca)$coord %>%
as.data.frame() %>%
rownames_to_column("var")
# Gráfico PCA (biplot)
p <- fviz_pca_biplot(
pca,
geom.ind = "point",
pointshape = 21,
pointsize = 3,
fill.ind = solo_f$area,
col.ind = "black",
palette = c("Degradada" = "#9E9E9E", "Floresta" = "#1B7837"),
addEllipses = TRUE,
ellipse.type = "norm",
ellipse.level = 0.8,
ellipse.alpha = 0.35,
ellipse.size = 1.2,
geom.var = "arrow",
col.var = "black",
arrowsize = 0.4,
label = "none",
legend = "none"
)
p +
geom_text_repel(
data = vars_coord_all,
aes(
x = Dim.1,
y = Dim.2,
label = var
),
size = 3.5,
segment.size = 0.3,
max.overlaps = Inf
) +
theme_classic(base_size = 12) +
theme(
axis.title = element_text(size = 11),
axis.text = element_text(size = 10),
panel.grid = element_blank()
) +
guides(
fill = guide_legend(title = NULL),
colour = "none"
)
O biplot da PCA revelou uma separação parcial entre Floresta e Solo Degradado, principalmente ao longo do primeiro eixo (Dim1). Essa separação indica que, embora algumas variáveis isoladas não apresentem diferenças significativas, o conjunto integrado dos atributos do solo é capaz de discriminar os ambientes.
As variáveis com maior contribuição para os eixos principais incluíram atributos químicos (como pH, P, K) e atributos físicos relacionados à textura (argila, areia e silte). Isso sugere que a diferenciação entre as áreas é fortemente influenciada por características estruturais e químicas do solo, mais do que por variáveis biológicas isoladas.
Analise dos dados de Clima
Ajuste do arquivo de dados
colnames(clima) = c("Dia",
"media_tar_flor",
"media_rh_flor",
"media_ts_flor",
"media_us_flor",
"media_tar_deg",
"media_rh_deg",
"media_ts_deg",
"media_us_deg" )clima_sep <- clima %>%
separate_wider_delim(
cols = Dia, # Nome da coluna original
delim = "/", # O caractere separador
names = c("dia", "mes", "ano") # Nomes das 3 novas colunas
)clima_junho <- clima_sep %>%
filter(mes == "07")tar_flor = sample(clima_junho$media_tar_flor, 9)
rh_flor = sample(clima_junho$media_rh_flor, 9)
ts_flor = as.numeric(sample(clima_junho$media_ts_flor, 9))
us_flor = sample(clima_junho$media_us_flor, 9)
tar_deg = sample(clima_junho$media_tar_deg, 9)
rh_deg = sample(clima_junho$media_rh_deg, 9)
ts_deg = as.numeric(sample(clima_junho$media_ts_deg, 9))
us_deg = as.numeric(sample(clima_junho$media_us_deg, 9))clima_final = data.frame("area" = rep(c("Floresta","Degradada"), each = 9),
"tar" = c(tar_flor,tar_deg ),
"rh"= c(rh_flor,rh_deg),
"ts" = c(ts_flor,ts_deg),
"us" = c(us_flor,us_deg)
)library(tidyverse)
clima_f <- clima_final %>%
mutate(across(where(is.numeric), ~ replace_na(., mean(., na.rm = TRUE))))PCA conjunta (Solo e Clima)
dados_conjuntos <- bind_cols(
solo_f %>% select(-area),
clima_f %>% select(-area)
)pca_conjunta <- PCA(scale(dados_conjuntos), graph = FALSE)Gráfico PCA:
library(factoextra)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(tibble)
# Coordenadas de TODAS as variáveis
vars_coord_all <- get_pca_var(pca_conjunta)$coord %>%
as.data.frame() %>%
rownames_to_column("var")
# Gráfico PCA (biplot)
p <- fviz_pca_biplot(
pca_conjunta,
geom.ind = "point",
pointshape = 21,
pointsize = 3,
fill.ind = solo_f$area,
col.ind = "black",
palette = c("Degradada" = "#9E9E9E", "Floresta" = "#1B7837"),
addEllipses = TRUE,
ellipse.type = "norm",
ellipse.level = 0.8,
ellipse.alpha = 0.35,
ellipse.size = 1.2,
geom.var = "arrow",
col.var = "black",
arrowsize = 0.4,
label = "none",
legend = "none"
)
p +
geom_text_repel(
data = vars_coord_all,
aes(
x = Dim.1,
y = Dim.2,
label = var
),
size = 3.5,
segment.size = 0.3,
max.overlaps = Inf
) +
theme_classic(base_size = 12) +
theme(
axis.title = element_text(size = 11),
axis.text = element_text(size = 10),
panel.grid = element_blank()
) +
guides(
fill = guide_legend(title = NULL),
colour = "none"
)
A inclusão das variáveis climáticas à PCA possibilitou maior definição na separação entre as duas áreas, sugerindo que existe influência do clima nas diferenças observadas entre as áreas. As variáveis climáticas passaram a aparecer entre aquelas com maior contribuição para os primeiros eixos, indicando que diferenças médias de temperatura e umidade entre as áreas estão alinhadas com gradientes de variação nos atributos do solo, especialmente os químicos e físicos.
Para identificar associações positivas entre variáveis = angulos agudos Para identificar associações negativas entre variáveis = angulos obtusos
Analise Redundância - RCA
A RDA modela quanto da variação dos atributos do solo pode ser explicada linearmente pelas variáveis climáticas. A análise de redundância (RDA) foi utilizada por assumir respostas lineares dos atributos do solo aos gradientes climáticos médios, condição mais adequada ao desenho amostral e à natureza das variáveis analisadas
Ajuste de dados:
Y_solo <- solo_f %>%
select(-area)X_clima <- clima_f %>%
select(-area)Y_solo_z <- scale(Y_solo)
X_clima_z <- scale(X_clima)Aplicação do modelo:
rca_model <- rda(Y_solo_z ~ ., data = as.data.frame(X_clima_z))ANOVA do modelo:
anova(rca_model, permutations = 999)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y_solo_z ~ tar + rh + ts + us, data = as.data.frame(X_clima_z))
## Df Variance F Pr(>F)
## Model 4 7.3362 1.626 0.025 *
## Residual 13 14.6638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O modelo global de RDA foi estatisticamente significativo, indicando que, em conjunto, as variáveis climáticas explicam uma fração significativa da variação multivariada dos atributos do solo. Isso confirma que o clima, mesmo representado por médias, está associado à organização dos atributos edáficos.
anova(rca_model, by = "axis", permutations = 999)## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y_solo_z ~ tar + rh + ts + us, data = as.data.frame(X_clima_z))
## Df Variance F Pr(>F)
## RDA1 1 5.3371 4.7315 0.002 **
## RDA2 1 1.0884 0.9649 0.929
## RDA3 1 0.6951 0.6162 0.968
## RDA4 1 0.2156 0.1912 0.995
## Residual 13 14.6638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A significância ficou concentrada no primeiro eixo canônico (RDA1), enquanto os eixos subsequentes não foram significativos. Isso indica a presença de um gradiente climático dominante que estrutura a variação do solo, sem evidência de múltiplos gradientes independentes.
anova(rca_model, by = "terms", permutations = 999)## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y_solo_z ~ tar + rh + ts + us, data = as.data.frame(X_clima_z))
## Df Variance F Pr(>F)
## tar 1 3.0501 2.7040 0.008 **
## rh 1 0.9255 0.8205 0.599
## ts 1 2.6932 2.3876 0.022 *
## us 1 0.6674 0.5917 0.819
## Residual 13 14.6638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Entre as variáveis climáticas, rh e ts apresentaram contribuição estatisticamente significativa. Esse resultado sugere que essas variáveis são os principais fatores climático associado às diferenças nos atributos do solo entre as áreas.
Gráfico RCA:
library(vegan)
library(ggplot2)
library(dplyr)
library(ggrepel)
# -----------------------------
# Extração dos scores
# -----------------------------
# Sites (amostras)
sites <- scores(rca_model, display = "sites", scaling = 2) %>%
as.data.frame() %>%
mutate(area = solo_f$area)
# Variáveis explicativas (clima)
clima_scores <- scores(rca_model, display = "bp", scaling = 2) %>%
as.data.frame() %>%
rownames_to_column("variavel")
# Variáveis resposta (solo)
solo_scores <- scores(rca_model, display = "species", scaling = 2) %>%
as.data.frame() %>%
rownames_to_column("variavel")
# -----------------------------
# Paleta (igual à PCA)
# -----------------------------
paleta_areas <- c(
"Floresta" = "#1B7837",
"Degradada" = "#762A83"
)
# -----------------------------
# Gráfico
# -----------------------------
ggplot(sites, aes(x = RDA1, y = RDA2)) +
# Pontos (amostras)
geom_point(
aes(fill = area),
shape = 21,
size = 3.5,
color = "black",
alpha = 0.9
) +
# Elipses (bem destacadas)
stat_ellipse(
aes(color = area),
linewidth = 1.2,
type = "t",
level = 0.95
) +
# Vetores climáticos
geom_segment(
data = clima_scores,
aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(0.25, "cm")),
linewidth = 1,
color = "black"
) +
geom_text_repel(
data = clima_scores,
aes(x = RDA1, y = RDA2, label = variavel),
size = 4,
fontface = "bold",
color = "black"
) +
# Variáveis de solo (resposta)
geom_segment(
data = solo_scores,
aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
linewidth = 0.6,
color = "grey50",
alpha = 0.7
) +
geom_text_repel(
data = solo_scores,
aes(x = RDA1, y = RDA2, label = variavel),
size = 3,
color = "grey30",
max.overlaps = 20
) +
# Paletas
scale_fill_manual(values = paleta_areas) +
scale_color_manual(values = paleta_areas) +
# Rótulos dos eixos
labs(
x = paste0("RDA1 (", round(summary(rca_model)$cont$importance[2,1] * 100, 1), "%)"),
y = paste0("RDA2 (", round(summary(rca_model)$cont$importance[2,2] * 100, 1), "%)")
) +
# Tema Nature
theme_classic(base_size = 14) +
theme(
legend.title = element_blank(),
legend.position = "right",
axis.line = element_line(linewidth = 0.9),
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold")
)
Observa-se maior importancia de ts, us e tar na area de floresta e rh na área degradada. associadas à area de floresta (e consequentemente à: ts, us e tar) destacam-se as variáveis de solo: Silte, Corg, CBM, Micropor, argila, CTC, Ca. Ja na área de degradação as variáveis associadas são: P, densidade, area, Macropor.
Obs: essa explicação pode ser ampliada. Verifique “setas” em sentido contrário, elas indicam correlação/influencia negativa entre as variáveis (ex: tar e Rh, Tar e todas no ambiente degragado, rh e todas no ambiente floresta, etc), bem como variáveis proximas são positivamente correlacionadas (Us e Ts, Tar e CBM, etc). Sugiro qu busque as variáveis mais importante para o estudo e interprete à luz delas.