script_aplicacao
DielsonMarcioFonsecaMaues
2026-02-26
============================================================================= ##### CAPITULO DE LIVRO ##### ##### ESTATÍSTICA MULTIVARIADA APLICADA AO CONTEXTO AGROPECUÁRIO #### =============================================================================
ANÁLISE DE VARIÂNCIA (ANOVA)
# 1) Inserindo os dados no R:
# Criando os vetores
Bloco = factor(rep(1:5, each = 4))
Dose_N = factor(rep(c(0, 50, 100, 150), times = 5))
Produtividade = c(
4200, 5100, 5800, 6000,
4300, 5200, 5900, 6100,
4100, 5000, 5700, 5900,
4250, 5150, 5850, 6050,
4350, 5300, 6000, 6200
)
# Criando o dataframe
dados = data.frame(Bloco, Dose_N, Produtividade)
# Visualizando os dados
head(dados)## Bloco Dose_N Produtividade
## 1 1 0 4200
## 2 1 50 5100
## 3 1 100 5800
## 4 1 150 6000
## 5 2 0 4300
## 6 2 50 5200# 2) Ajustando o modelo de ANOVA:
modelo = aov(Produtividade ~ Dose_N + Bloco, data = dados)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## Dose_N 3 10045375 3348458 26788 < 2e-16 ***
## Bloco 4 185500 46375 371 1.85e-12 ***
## Residuals 12 1500 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 3) Verificação dos pressupostos:
# Normalidade dos resíduos
shapiro.test(residuals(modelo))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.69434, p-value = 3.328e-05# Homogeneidade das variâncias
bartlett.test(Produtividade ~ Dose_N, data = dados)##
## Bartlett test of homogeneity of variances
##
## data: Produtividade by Dose_N
## Bartlett's K-squared = 0.1169, df = 3, p-value = 0.9897# 4) Transformação de dados - logarítmica:
dados$Prod_log = log(dados$Produtividade)
modelo_log = aov(Prod_log ~ Dose_N + Bloco, data = dados)
summary(modelo_log)## Df Sum Sq Mean Sq F value Pr(>F)
## Dose_N 3 0.3890 0.12968 18875.8 < 2e-16 ***
## Bloco 4 0.0067 0.00168 244.1 2.22e-11 ***
## Residuals 12 0.0001 0.00001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 5) Verificando os pressupostos novamente:
shapiro.test(residuals(modelo_log))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_log)
## W = 0.92377, p-value = 0.1171bartlett.test(Prod_log ~ Dose_N, data = dados)##
## Bartlett test of homogeneity of variances
##
## data: Prod_log by Dose_N
## Bartlett's K-squared = 0.21732, df = 3, p-value = 0.9747str(dados)## 'data.frame': 20 obs. of 4 variables:
## $ Bloco : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ Dose_N : Factor w/ 4 levels "0","50","100",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ Produtividade: num 4200 5100 5800 6000 4300 5200 5900 6100 4100 5000 ...
## $ Prod_log : num 8.34 8.54 8.67 8.7 8.37 ...# 6) Teste de comparação de médias (Tukey):
TukeyHSD(modelo_log, "Dose_N")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Prod_log ~ Dose_N + Bloco, data = dados)
##
## $Dose_N
## diff lwr upr p adj
## 50-0 0.19445185 0.18953028 0.19937343 0
## 100-0 0.32193924 0.31701767 0.32686082 0
## 150-0 0.35556536 0.35064379 0.36048693 0
## 100-50 0.12748739 0.12256582 0.13240897 0
## 150-50 0.16111351 0.15619194 0.16603508 0
## 150-100 0.03362612 0.02870454 0.03854769 0# 7) Gráfico Ilustrativo:
boxplot(dados$Prod_log ~ Dose_N,
xlab = "Doses de Nitrogênio (kg ha-1)",
ylab = "Produtividade (kg ha-1)",
col = "lightgreen")
boxplot(Produtividade ~ Dose_N,
xlab = "Doses de Nitrogênio (kg ha-1)",
ylab = "Produtividade (kg ha-1)",
col = "green")
=============================================================================
ANÁLISE DE CORRELAÇÃO (CA)
# 1) Construindo Dataset:
set.seed(456)
n <- 45 # número de talhões ou amostras
# Solo
pH <- rnorm(n, 5.5, 0.4)
MO <- rnorm(n, 25, 5) # g/dm³
P <- rnorm(n, 12, 4) # mg/dm³
K <- rnorm(n, 0.25, 0.08) # cmolc/dm³
Ca <- rnorm(n, 3.5, 0.8)
Mg <- rnorm(n, 1.2, 0.3)
# Planta (nutrição foliar)
N_foliar <- rnorm(n, 28, 3) # g/kg
K_foliar <- rnorm(n, 18, 2)
Ca_foliar <- rnorm(n, 12, 2)
# Produtividade (influenciada por solo e nutrição)
Produtividade <- 20 +
(MO*0.3) +
(P*0.2) +
(Ca*0.5) +
(N_foliar*0.4) +
rnorm(n, 0, 3)
# Qualidade (nota de bebida)
Qualidade <- 70 +
(MO*0.2) +
(K_foliar*0.5) -
(pH*1.5) +
rnorm(n, 0, 2)
# Dataframe final
dados_cafe <- data.frame(
pH, MO, P, K, Ca, Mg,
N_foliar, K_foliar, Ca_foliar,
Produtividade, Qualidade
)
head(dados_cafe)## pH MO P K Ca Mg N_foliar
## 1 4.962591 25.60426 10.107055 0.16745827 2.389449 1.3800834 20.85503
## 2 5.748710 25.59075 9.323176 0.08423626 3.010651 1.3594260 28.45511
## 3 5.820350 28.85027 17.434107 0.21591463 3.754796 0.7531238 22.28668
## 4 4.944443 19.12299 17.651672 0.18784195 3.311675 1.6069946 32.96279
## 5 5.214257 27.04519 9.068905 0.23545502 4.360042 1.4009892 29.49625
## 6 5.370376 21.67525 11.677792 0.18101541 3.419364 1.6973276 27.66870
## K_foliar Ca_foliar Produtividade Qualidade
## 1 18.74925 13.62733 43.43299 75.85353
## 2 19.90963 13.83555 44.87583 78.54782
## 3 18.27132 10.66261 46.04109 76.63333
## 4 19.48889 13.64474 49.03802 74.52419
## 5 21.60345 13.48275 46.73890 78.08979
## 6 17.90111 10.32615 37.30197 76.13576# 2) Matriz de correlação:
correlacao <- cor(dados_cafe, method = "pearson")
round(correlacao, 2)## pH MO P K Ca Mg N_foliar K_foliar
## pH 1.00 -0.04 0.02 0.03 0.08 0.01 -0.18 0.10
## MO -0.04 1.00 -0.32 0.06 -0.04 -0.07 -0.12 0.07
## P 0.02 -0.32 1.00 0.06 0.10 0.09 -0.14 -0.11
## K 0.03 0.06 0.06 1.00 0.25 0.12 -0.13 -0.05
## Ca 0.08 -0.04 0.10 0.25 1.00 0.08 0.04 0.07
## Mg 0.01 -0.07 0.09 0.12 0.08 1.00 0.11 0.29
## N_foliar -0.18 -0.12 -0.14 -0.13 0.04 0.11 1.00 0.17
## K_foliar 0.10 0.07 -0.11 -0.05 0.07 0.29 0.17 1.00
## Ca_foliar -0.02 0.09 0.03 0.11 -0.14 0.40 -0.02 0.21
## Produtividade -0.06 0.23 0.10 -0.07 0.13 -0.19 0.18 0.21
## Qualidade 0.00 0.52 -0.26 0.06 0.01 0.02 0.11 0.13
## Ca_foliar Produtividade Qualidade
## pH -0.02 -0.06 0.00
## MO 0.09 0.23 0.52
## P 0.03 0.10 -0.26
## K 0.11 -0.07 0.06
## Ca -0.14 0.13 0.01
## Mg 0.40 -0.19 0.02
## N_foliar -0.02 0.18 0.11
## K_foliar 0.21 0.21 0.13
## Ca_foliar 1.00 -0.11 -0.14
## Produtividade -0.11 1.00 0.05
## Qualidade -0.14 0.05 1.00# 3) Teste de significância
install.packages("Hmisc")## Instalando pacote em 'C:/Users/diels/AppData/Local/R/win-library/4.6'
## (como 'lib' não foi especificado)## package 'Hmisc' successfully unpacked and MD5 sums checked
##
## Os pacotes binários baixados estão em
## C:\Users\diels\AppData\Local\Temp\Rtmp8YBV6a\downloaded_packageslibrary(Hmisc)## Registered S3 method overwritten by 'htmlwidgets':
## method from
## print.htmlwidget tools:rstudio## Registered S3 method overwritten by 'data.table':
## method from
## print.data.table##
## Anexando pacote: 'Hmisc'## Os seguintes objetos são mascarados por 'package:base':
##
## format.pval, unitsresultado <- rcorr(as.matrix(dados_cafe))
resultado$r # correlação## pH MO P K
## pH 1.000000000 -0.03677500 0.02409038 0.02847610
## MO -0.036774996 1.00000000 -0.31950598 0.05897776
## P 0.024090385 -0.31950598 1.00000000 0.06057842
## K 0.028476102 0.05897776 0.06057842 1.00000000
## Ca 0.079873939 -0.03803321 0.10486283 0.24718274
## Mg 0.007687873 -0.07463606 0.08961305 0.11629821
## N_foliar -0.176984774 -0.11527467 -0.14078644 -0.13427244
## K_foliar 0.103364138 0.07365793 -0.11411009 -0.05450790
## Ca_foliar -0.019435162 0.09006604 0.02547473 0.11253380
## Produtividade -0.064270343 0.23244938 0.09951363 -0.06878564
## Qualidade 0.003151234 0.51960199 -0.25673045 0.06245014
## Ca Mg N_foliar K_foliar
## pH 0.07987394 0.007687873 -0.17698477 0.10336414
## MO -0.03803321 -0.074636061 -0.11527467 0.07365793
## P 0.10486283 0.089613049 -0.14078644 -0.11411009
## K 0.24718274 0.116298209 -0.13427244 -0.05450790
## Ca 1.00000000 0.082393364 0.03929497 0.06596230
## Mg 0.08239336 1.000000000 0.10513728 0.29292114
## N_foliar 0.03929497 0.105137275 1.00000000 0.16735915
## K_foliar 0.06596230 0.292921138 0.16735915 1.00000000
## Ca_foliar -0.14045951 0.401505672 -0.02351090 0.21259886
## Produtividade 0.13290536 -0.186186969 0.17846976 0.20725596
## Qualidade 0.01032238 0.015589698 0.11301564 0.13286794
## Ca_foliar Produtividade Qualidade
## pH -0.01943516 -0.06427034 0.003151234
## MO 0.09006604 0.23244938 0.519601989
## P 0.02547473 0.09951363 -0.256730452
## K 0.11253380 -0.06878564 0.062450143
## Ca -0.14045951 0.13290536 0.010322383
## Mg 0.40150567 -0.18618697 0.015589698
## N_foliar -0.02351090 0.17846976 0.113015636
## K_foliar 0.21259886 0.20725596 0.132867942
## Ca_foliar 1.00000000 -0.10721922 -0.139750839
## Produtividade -0.10721922 1.00000000 0.054035615
## Qualidade -0.13975084 0.05403561 1.000000000resultado$P # p-valores## pH MO P K Ca
## pH NA 0.8104595570 0.87518309 0.8526913 0.6019715
## MO 0.8104596 NA 0.03240106 0.7003583 0.8040998
## P 0.8751831 0.0324010605 NA 0.6926212 0.4930032
## K 0.8526913 0.7003582635 0.69262119 NA 0.1016293
## Ca 0.6019715 0.8040998497 0.49300320 0.1016293 NA
## Mg 0.9600257 0.6260702346 0.55827465 0.4467852 0.5905257
## N_foliar 0.2448093 0.4508224305 0.35629195 0.3791997 0.7977347
## K_foliar 0.4992381 0.6306144908 0.45544018 0.7221212 0.6668274
## Ca_foliar 0.8991635 0.5562789830 0.86807220 0.4617311 0.3574212
## Produtividade 0.6748943 0.1243982983 0.51544034 0.6534490 0.3841158
## Qualidade 0.9836092 0.0002543963 0.08867998 0.6836130 0.9463446
## Mg N_foliar K_foliar Ca_foliar Produtividade
## pH 0.960025744 0.2448093 0.4992381 0.899163499 0.6748943
## MO 0.626070235 0.4508224 0.6306145 0.556278983 0.1243983
## P 0.558274646 0.3562920 0.4554402 0.868072196 0.5154403
## K 0.446785178 0.3791997 0.7221212 0.461731084 0.6534490
## Ca 0.590525680 0.7977347 0.6668274 0.357421218 0.3841158
## Mg NA 0.4918658 0.0508493 0.006263371 0.2207352
## N_foliar 0.491865832 NA 0.2718294 0.878162592 0.2408088
## K_foliar 0.050849300 0.2718294 NA 0.160881903 0.1719146
## Ca_foliar 0.006263371 0.8781626 0.1608819 NA 0.4832821
## Produtividade 0.220735201 0.2408088 0.1719146 0.483282115 NA
## Qualidade 0.919040759 0.4598031 0.3842509 0.359876443 0.7244338
## Qualidade
## pH 0.9836091763
## MO 0.0002543963
## P 0.0886799763
## K 0.6836130469
## Ca 0.9463446052
## Mg 0.9190407589
## N_foliar 0.4598031314
## K_foliar 0.3842508869
## Ca_foliar 0.3598764427
## Produtividade 0.7244338055
## Qualidade NA# 4) Heatmap (gráfico científico):
install.packages("corrplot")## Instalando pacote em 'C:/Users/diels/AppData/Local/R/win-library/4.6'
## (como 'lib' não foi especificado)## package 'corrplot' successfully unpacked and MD5 sums checked
##
## Os pacotes binários baixados estão em
## C:\Users\diels\AppData\Local\Temp\Rtmp8YBV6a\downloaded_packageslibrary(corrplot)## corrplot 0.95 loadedcorrplot(correlacao,
method = "color", ## diferenciação em gradiente
type = "upper",
tl.col = "black",
addCoef.col = "black", ## config
number.cex = 0.7)
# 5) Produtividade:
cor_prod <- correlacao[,"Produtividade"]
sort(cor_prod, decreasing = TRUE)## Produtividade MO K_foliar N_foliar Ca
## 1.00000000 0.23244938 0.20725596 0.17846976 0.13290536
## P Qualidade pH K Ca_foliar
## 0.09951363 0.05403561 -0.06427034 -0.06878564 -0.10721922
## Mg
## -0.18618697# 6) Qualidade:
cor_qualidade <- correlacao[,"Qualidade"]
sort(cor_qualidade, decreasing = TRUE)## Qualidade MO K_foliar N_foliar K
## 1.000000000 0.519601989 0.132867942 0.113015636 0.062450143
## Produtividade Mg Ca pH Ca_foliar
## 0.054035615 0.015589698 0.010322383 0.003151234 -0.139750839
## P
## -0.256730452=============================================================================
ANÁLISE DE CORRELAÇÃO CANÔNICA (CCA)
# 1) Gerando Dataset:
set.seed(555)
n <- 40 # número de lotes ou galpões
# Conjunto X (ambiente e manejo)
Temperatura <- rnorm(n, 30, 3) # °C
Umidade <- rnorm(n, 65, 10) # %
Densidade <- rnorm(n, 12, 2) # aves/m²
Ventilacao <- rnorm(n, 3, 0.8) # m/s
Amônia <- rnorm(n, 15, 5) # ppm
X <- data.frame(Temperatura, Umidade, Densidade, Ventilacao, Amônia)
# Conjunto Y (desempenho)
Ganho_Peso <- 2500 -
(Temperatura*20) -
(Amônia*10) +
(Ventilacao*50) +
rnorm(n, 0, 100)
Conversao_Alimentar <- 1.8 +
(Temperatura*0.02) +
(Amônia*0.03) -
(Ventilacao*0.05) +
rnorm(n, 0, 0.1)
Mortalidade <- 3 +
(Temperatura*0.2) +
(Amônia*0.3) +
rnorm(n, 0, 1)
Y <- data.frame(Ganho_Peso, Conversao_Alimentar, Mortalidade)# 2) Rodando a CCA:
cca_modelo <- cancor(X, Y)
cca_modelo## $cor
## [1] 0.9291258 0.4556859 0.2711185
##
## $xcoef
## [,1] [,2] [,3] [,4]
## Temperatura -0.030850322 0.03201707 -0.02756836 -0.0360477857
## Umidade 0.001440886 -0.01157135 -0.00758158 -0.0018553341
## Densidade -0.009713637 0.01607681 0.03499769 0.0243114349
## Ventilacao 0.043157781 -0.01893031 0.11776673 -0.1255244787
## Amônia -0.030467717 -0.01511599 0.01565405 0.0004618982
## [,5]
## Temperatura -0.018325785
## Umidade -0.004982992
## Densidade -0.060110046
## Ventilacao -0.026225359
## Amônia 0.005656036
##
## $ycoef
## [,1] [,2] [,3]
## Ganho_Peso 0.0003068096 -0.001211484 -0.0007420482
## Conversao_Alimentar -0.5870991464 -0.238037223 -1.2243536747
## Mortalidade -0.0309673089 -0.074383450 0.0965605121
##
## $xcenter
## Temperatura Umidade Densidade Ventilacao Amônia
## 29.882265 66.096829 11.960281 3.003791 15.294283
##
## $ycenter
## Ganho_Peso Conversao_Alimentar Mortalidade
## 1903.617407 2.689853 13.817638# 3) Correlações canônicas:
cca_modelo$cor## [1] 0.9291258 0.4556859 0.2711185# 4) Coeficientes canônicos:
cca_modelo$xcoef## [,1] [,2] [,3] [,4]
## Temperatura -0.030850322 0.03201707 -0.02756836 -0.0360477857
## Umidade 0.001440886 -0.01157135 -0.00758158 -0.0018553341
## Densidade -0.009713637 0.01607681 0.03499769 0.0243114349
## Ventilacao 0.043157781 -0.01893031 0.11776673 -0.1255244787
## Amônia -0.030467717 -0.01511599 0.01565405 0.0004618982
## [,5]
## Temperatura -0.018325785
## Umidade -0.004982992
## Densidade -0.060110046
## Ventilacao -0.026225359
## Amônia 0.005656036cca_modelo$ycoef## [,1] [,2] [,3]
## Ganho_Peso 0.0003068096 -0.001211484 -0.0007420482
## Conversao_Alimentar -0.5870991464 -0.238037223 -1.2243536747
## Mortalidade -0.0309673089 -0.074383450 0.0965605121# 5) Plot das variáveis canônicas:
U <- as.matrix(X) %*% cca_modelo$xcoef
V <- as.matrix(Y) %*% cca_modelo$ycoef
plot(U[,1], V[,1],
xlab = "Canônica Ambiental",
ylab = "Canônica Desempenho",
main = "CCA - Avicultura",
pch = 19)
abline(lm(V[,1] ~ U[,1]), col = "blue")
=============================================================================
ANÁLISE DE SIMILARIDADE (ANOSIM)
# 1) Criando Dataset:
set.seed(123) # Reprodutibilidade
# Tratamentos
tratamento = factor(rep(c("CONV", "PD", "ILP", "ORG"), each = 12))
# Simulação com base em padrões reais de campo
pH = c(
rnorm(12, 5.2, 0.2), # CONV (mais ácido)
rnorm(12, 5.6, 0.2), # PD
rnorm(12, 5.8, 0.2), # ILP
rnorm(12, 6.2, 0.2) # ORG (melhor correção)
)
MO = c(
rnorm(12, 18, 3),
rnorm(12, 25, 4),
rnorm(12, 28, 4),
rnorm(12, 35, 5)
)
P = c(
rnorm(12, 6, 2),
rnorm(12, 10, 3),
rnorm(12, 12, 3),
rnorm(12, 18, 4)
)
K = c(
rnorm(12, 0.15, 0.05),
rnorm(12, 0.25, 0.05),
rnorm(12, 0.30, 0.06),
rnorm(12, 0.40, 0.07)
)
Ca = c(
rnorm(12, 2.0, 0.5),
rnorm(12, 3.0, 0.6),
rnorm(12, 3.5, 0.7),
rnorm(12, 4.5, 0.8)
)
Mg = c(
rnorm(12, 0.8, 0.2),
rnorm(12, 1.2, 0.3),
rnorm(12, 1.4, 0.3),
rnorm(12, 1.8, 0.4)
)
CTC = Ca + Mg + K + rnorm(48, 2, 0.5)
# Data frame final
dados = data.frame(tratamento, pH, MO, P, K, Ca, Mg, CTC)
# Visualizar
head(dados)## tratamento pH MO P K Ca Mg
## 1 CONV 5.087905 20.33990 10.374666 0.06992319 2.047292 0.6422756
## 2 CONV 5.153965 17.74989 9.065221 0.12345467 1.552318 0.6995603
## 3 CONV 5.511742 18.75996 5.528599 0.07691222 1.344599 1.0992121
## 4 CONV 5.214102 17.91436 3.947158 0.18439584 2.998607 0.5725393
## 5 CONV 5.225858 17.87139 4.579187 0.25500545 2.300354 0.7641897
## 6 CONV 5.543013 22.10581 6.513767 0.08564848 1.374364 1.1804724
## CTC
## 1 4.649965
## 2 4.459366
## 3 5.104916
## 4 6.282632
## 5 5.892181
## 6 4.351751# 2) Uso de direto da ANOSIM:
library(vegan)## Carregando pacotes exigidos: permute## This is vegan 2.7-3variaveis = dados[,2:8]
distancia = vegdist(variaveis, method = "bray")
anosim_resultado = anosim(distancia, dados$tratamento)
summary(anosim_resultado)##
## Call:
## anosim(x = distancia, grouping = dados$tratamento)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.7108
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0490 0.0666 0.0807 0.1037
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 435.75 686.5 912.25 1128 864
## CONV 7 124.25 262.0 431.50 841 66
## ILP 1 76.25 180.5 268.25 423 66
## ORG 5 102.50 236.5 346.50 606 66
## PD 8 118.75 261.5 502.00 825 66=============================================================================
ANÁLISE DE REDUNDÂNCIA (RDA)
# 1) Criando Dataset:
set.seed(202)
n <- 36 # número de vacas/amostras
# Variáveis explicativas (nutrição e ambiente)
Proteina_Dieta <- rnorm(n, 16, 1.5) # %
Energia_Dieta <- rnorm(n, 68, 3) # NDT %
Fibra <- rnorm(n, 32, 4) # %
Consumo_MS <- rnorm(n, 18, 2) # kg/dia
Temperatura <- rnorm(n, 27, 2) # °C
Umidade <- rnorm(n, 70, 10) # %
# Variáveis resposta (produção leiteira)
Producao_Leite <- 20 +
(Proteina_Dieta*0.8) +
(Energia_Dieta*0.3) +
(Consumo_MS*1.2) -
(Temperatura*0.5) +
rnorm(n, 0, 2)
Gordura_Leite <- 3.5 +
(Fibra*0.05) -
(Energia_Dieta*0.02) +
rnorm(n, 0, 0.2)
Proteina_Leite <- 3.2 +
(Proteina_Dieta*0.1) +
rnorm(n, 0, 0.15)
# Dataframe
dados_leite <- data.frame(
Proteina_Dieta, Energia_Dieta, Fibra, Consumo_MS,
Temperatura, Umidade,
Producao_Leite, Gordura_Leite, Proteina_Leite
)
head(dados_leite)## Proteina_Dieta Energia_Dieta Fibra Consumo_MS Temperatura Umidade
## 1 14.30232 65.42386 32.11432 19.46783 30.09381 77.34628
## 2 15.33926 69.53735 31.63759 16.99656 28.80279 58.52731
## 3 15.49540 62.75073 31.60210 16.23006 26.56241 73.16835
## 4 14.72943 65.78132 29.62890 13.32423 26.99924 79.62614
## 5 15.78096 73.50795 27.64772 17.02677 26.77468 67.16852
## 6 13.85418 69.79082 22.56636 20.55297 28.11070 50.57001
## Producao_Leite Gordura_Leite Proteina_Leite
## 1 59.03745 3.960609 4.621453
## 2 58.51138 3.863942 4.711956
## 3 57.92708 4.036746 4.406387
## 4 57.24169 4.059189 4.712307
## 5 60.72601 3.647397 4.924965
## 6 58.50708 3.354951 4.486667# 2) Separando variáveis:
# Variáveis explicativas (X)
X <- dados_leite[, c("Proteina_Dieta", "Energia_Dieta", "Fibra",
"Consumo_MS", "Temperatura", "Umidade")]
# Variáveis resposta (Y)
Y <- dados_leite[, c("Producao_Leite", "Gordura_Leite", "Proteina_Leite")]# 3) Rodando RDA:
library(vegan)
rda_modelo <- rda(Y ~ ., data = X)
summary(rda_modelo)##
## Call:
## rda(formula = Y ~ Proteina_Dieta + Energia_Dieta + Fibra + Consumo_MS + Temperatura + Umidade, data = X)
##
## Partitioning of variance:
## Inertia Proportion
## Total 9.084 1.0000
## Constrained 4.831 0.5318
## Unconstrained 4.253 0.4682
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 PC1 PC2 PC3
## Eigenvalue 4.7531 0.051851 0.026121 4.1978 0.037043 0.018330
## Proportion Explained 0.5232 0.005708 0.002875 0.4621 0.004078 0.002018
## Cumulative Proportion 0.5232 0.528933 0.531808 0.9939 0.997982 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3
## Eigenvalue 4.7531 0.05185 0.026121
## Proportion Explained 0.9839 0.01073 0.005407
## Cumulative Proportion 0.9839 0.99459 1.000000# 4) Testes de significância:
anova(rda_modelo, permutations = 999)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Y ~ Proteina_Dieta + Energia_Dieta + Fibra + Consumo_MS + Temperatura + Umidade, data = X)
## Df Variance F Pr(>F)
## Model 6 4.8310 5.4901 0.001 ***
## Residual 29 4.2531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(rda_modelo, 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 ~ Proteina_Dieta + Energia_Dieta + Fibra + Consumo_MS + Temperatura + Umidade, data = X)
## Df Variance F Pr(>F)
## RDA1 1 4.7531 32.4088 0.002 **
## RDA2 1 0.0519 0.4023 0.821
## RDA3 1 0.0261 0.2088 0.821
## Residual 32 4.2531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(rda_modelo, 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 ~ Proteina_Dieta + Energia_Dieta + Fibra + Consumo_MS + Temperatura + Umidade, data = X)
## Df Variance F Pr(>F)
## Proteina_Dieta 1 1.0519 7.1727 0.010 **
## Energia_Dieta 1 0.1874 1.2779 0.273
## Fibra 1 0.1040 0.7089 0.410
## Consumo_MS 1 2.4124 16.4488 0.001 ***
## Temperatura 1 1.0579 7.2131 0.009 **
## Umidade 1 0.0175 0.1191 0.726
## Residual 29 4.2531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 5) Biplot:
plot(rda_modelo, scaling = 2)
=============================================================================
ANÁLISE DE CORRESPONDÊNCIA MÚLTIPLA (MCA)
# 1) Criando Dataset:
set.seed(101)
n <- 40
Sistema <- sample(c("Convencional", "Orgânico", "Agroecológico"), n, replace = TRUE)
Tecnologia <- sample(c("Baixa", "Média", "Alta"), n, replace = TRUE)
Irrigacao <- sample(c("Sim", "Não"), n, replace = TRUE)
Adubacao <- sample(c("Química", "Orgânica", "Mista"), n, replace = TRUE)
Assistencia <- sample(c("Sim", "Não"), n, replace = TRUE)
# Produtividade categorizada
Produtividade <- sample(c("Baixa", "Média", "Alta"), n, replace = TRUE)
dados_mca <- data.frame(
Sistema,
Tecnologia,
Irrigacao,
Adubacao,
Assistencia,
Produtividade
)
head(dados_mca)## Sistema Tecnologia Irrigacao Adubacao Assistencia Produtividade
## 1 Convencional Alta Sim Orgânica Sim Média
## 2 Convencional Média Sim Química Não Média
## 3 Orgânico Média Sim Orgânica Sim Baixa
## 4 Agroecológico Baixa Sim Química Sim Alta
## 5 Agroecológico Média Não Mista Sim Média
## 6 Convencional Baixa Não Orgânica Sim Baixa# 2) Carregando pacote:
install.packages("FactoMineR")## Instalando pacote em 'C:/Users/diels/AppData/Local/R/win-library/4.6'
## (como 'lib' não foi especificado)## package 'FactoMineR' successfully unpacked and MD5 sums checked
##
## Os pacotes binários baixados estão em
## C:\Users\diels\AppData\Local\Temp\Rtmp8YBV6a\downloaded_packagesinstall.packages("factoextra")## Instalando pacote em 'C:/Users/diels/AppData/Local/R/win-library/4.6'
## (como 'lib' não foi especificado)## package 'factoextra' successfully unpacked and MD5 sums checked
##
## Os pacotes binários baixados estão em
## C:\Users\diels\AppData\Local\Temp\Rtmp8YBV6a\downloaded_packageslibrary(FactoMineR)
library(factoextra)## Carregando pacotes exigidos: ggplot2## Welcome to factoextra!## Want to learn more? See two factoextra-related books at https://www.datanovia.com/en/product/practical-guide-to-principal-component-methods-in-r/# 2.1) Rodando ACM:
mca_modelo <- MCA(dados_mca, graph = FALSE)
summary(mca_modelo)##
## Call:
## MCA(X = dados_mca, graph = FALSE)
##
##
## Eigenvalues
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## Variance 0.333 0.306 0.232 0.185 0.168 0.130
## % of var. 19.966 18.381 13.900 11.077 10.053 7.776
## Cumulative % of var. 19.966 38.347 52.247 63.324 73.377 81.154
## Dim.7 Dim.8 Dim.9 Dim.10
## Variance 0.122 0.098 0.059 0.035
## % of var. 7.305 5.887 3.563 2.092
## Cumulative % of var. 88.458 94.345 97.908 100.000
##
## Individuals (the 10 first)
## Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3
## 1 | 0.745 4.165 0.216 | 0.713 4.144 0.198 | 1.069
## 2 | -0.415 1.293 0.100 | 0.137 0.154 0.011 | 0.719
## 3 | 0.856 5.499 0.433 | -0.618 3.115 0.226 | -0.043
## 4 | -0.155 0.181 0.017 | 0.401 1.310 0.110 | -0.301
## 5 | -0.071 0.037 0.003 | 0.644 3.385 0.273 | -0.611
## 6 | 0.465 1.626 0.130 | -0.504 2.074 0.152 | 0.185
## 7 | 0.509 1.944 0.158 | -0.483 1.904 0.143 | -0.086
## 8 | 0.178 0.237 0.020 | 0.695 3.943 0.300 | -0.604
## 9 | 0.427 1.370 0.079 | 0.447 1.632 0.087 | 0.853
## 10 | 0.385 1.116 0.056 | 0.848 5.875 0.270 | 0.867
## ctr cos2
## 1 12.323 0.446 |
## 2 5.579 0.302 |
## 3 0.020 0.001 |
## 4 0.975 0.062 |
## 5 4.031 0.246 |
## 6 0.368 0.020 |
## 7 0.079 0.004 |
## 8 3.941 0.227 |
## 9 7.848 0.315 |
## 10 8.109 0.282 |
##
## Categories (the 10 first)
## Dim.1 ctr cos2 v.test Dim.2 ctr cos2
## Agroecológico | -0.385 2.594 0.080 -1.763 | 0.920 16.112 0.456
## Convencional | -0.332 1.932 0.059 -1.521 | -0.099 0.186 0.005
## Orgânico | 0.836 10.503 0.300 3.418 | -0.958 14.972 0.393
## Tecnologia_Alta | 0.406 1.237 0.029 1.064 | 1.162 11.012 0.238
## Tecnologia_Baixa | 0.409 3.138 0.100 1.977 | -0.110 0.248 0.007
## Tecnologia_Média | -0.451 4.835 0.184 -2.678 | -0.280 2.023 0.071
## Irrigacao_Não | -0.495 6.448 0.271 -3.251 | -0.309 2.727 0.106
## Irrigacao_Sim | 0.547 7.127 0.271 3.251 | 0.342 3.015 0.106
## Mista | 0.166 0.554 0.018 0.848 | 0.473 4.867 0.149
## Orgânica | 1.036 13.439 0.358 3.735 | -0.373 1.891 0.046
## v.test Dim.3 ctr cos2 v.test
## Agroecológico 4.215 | -0.584 8.591 0.184 -2.677 |
## Convencional -0.453 | 0.878 19.401 0.415 4.022 |
## Orgânico -3.916 | -0.343 2.533 0.050 -1.401 |
## Tecnologia_Alta 3.048 | 1.509 24.563 0.402 3.958 |
## Tecnologia_Baixa -0.533 | -0.255 1.757 0.039 -1.235 |
## Tecnologia_Média -1.662 | -0.275 2.583 0.068 -1.633 |
## Irrigacao_Não -2.029 | -0.277 2.895 0.085 -1.818 |
## Irrigacao_Sim 2.029 | 0.306 3.200 0.085 1.818 |
## Mista 2.411 | -0.577 9.574 0.222 -2.941 |
## Orgânica -1.344 | 0.446 3.576 0.066 1.608 |
##
## Categorical variables (eta2)
## Dim.1 Dim.2 Dim.3
## Sistema | 0.300 0.575 0.424 |
## Tecnologia | 0.184 0.244 0.402 |
## Irrigacao | 0.271 0.106 0.085 |
## Adubacao | 0.582 0.151 0.223 |
## Assistencia | 0.346 0.079 0.189 |
## Produtividade | 0.313 0.684 0.067 |# 3) Variância explicada:
mca_modelo$eig## eigenvalue percentage of variance
## dim 1 0.33276214 19.965728
## dim 2 0.30635400 18.381240
## dim 3 0.23166455 13.899873
## dim 4 0.18462044 11.077226
## dim 5 0.16754953 10.052972
## dim 6 0.12960825 7.776495
## dim 7 0.12174310 7.304586
## dim 8 0.09811981 5.887189
## dim 9 0.05938338 3.563003
## dim 10 0.03486146 2.091688
## cumulative percentage of variance
## dim 1 19.96573
## dim 2 38.34697
## dim 3 52.24684
## dim 4 63.32407
## dim 5 73.37704
## dim 6 81.15353
## dim 7 88.45812
## dim 8 94.34531
## dim 9 97.90831
## dim 10 100.00000# 4) Biplot da ACM com agrupamentos
fviz_mca_biplot(
mca_modelo,
# Colorir por sistema produtivo
habillage = dados_mca$Sistema,
# Adicionar elipses dos grupos
addEllipses = TRUE,
# Evitar sobreposição dos textos
repel = TRUE,
# Mostrar apenas variáveis/categorias
label = "var",
# Tamanho dos pontos
pointsize = 2.5
) +
theme_classic() +
labs(
title = "Análise de Correspondência Múltipla (ACM)",
subtitle = "Agrupamentos associados ao perfil tecnológico e produtivo",
x = "Dimensão 1",
y = "Dimensão 2"
)
# 5) Contribuição das variáveis:
fviz_contrib(mca_modelo, choice = "var", axes = 1)
fviz_contrib(mca_modelo, choice = "var", axes = 2)
=============================================================================
ANÁLISE DE PORCENTAGEM DE SIMILARIDADE (SIMPER)
# 1) Criando Dataset:
set.seed(333)
sistema <- factor(rep(c("EXT", "SEMI", "INT"), each = 12))
Oxigenio <- c(rnorm(12, 7, 0.5), rnorm(12, 5.5, 0.6), rnorm(12, 4.5, 0.7))
pH <- c(rnorm(12, 7.0, 0.3), rnorm(12, 7.5, 0.3), rnorm(12, 8.0, 0.4))
Temperatura <- c(rnorm(12, 26, 1), rnorm(12, 28, 1), rnorm(12, 30, 1.5))
Amonia <- c(rnorm(12, 0.2, 0.05), rnorm(12, 0.5, 0.1), rnorm(12, 1.0, 0.2))
Turbidez <- c(rnorm(12, 10, 3), rnorm(12, 20, 5), rnorm(12, 35, 8))
dados <- data.frame(Oxigenio, pH, Temperatura, Amonia, Turbidez)# 2) Rondado o SIMPER:
install.packages("vegan")## Warning: o pacote 'vegan' está em uso e não será instaladolibrary(vegan)
simper_result <- simper(dados, sistema)
summary(simper_result)##
## Contrast: EXT_SEMI
##
## average sd ratio ava avb cumsum p
## Turbidez 0.08730 0.04920 1.77430 11.09900 20.94100 0.679 0.950
## Temperatura 0.02026 0.01014 1.99710 25.93300 28.20500 0.837 0.025 *
## Oxigenio 0.01264 0.00572 2.21070 6.91400 5.48200 0.935 0.006 **
## pH 0.00577 0.00312 1.84970 6.95200 7.59500 0.980 0.008 **
## Amonia 0.00257 0.00095 2.71240 0.21200 0.50600 1.000 0.997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Contrast: EXT_INT
##
## average sd ratio ava avb cumsum p
## Turbidez 0.18029 0.05000 3.60600 11.09900 34.51000 0.745 0.001
## Temperatura 0.02940 0.01327 2.21600 25.93300 29.73000 0.867 0.001
## Oxigenio 0.01817 0.00692 2.62600 6.91400 4.58000 0.942 0.001
## pH 0.00733 0.00419 1.74900 6.95200 7.88000 0.972 0.001
## Amonia 0.00670 0.00187 3.59200 0.21200 1.07000 1.000 0.001
##
## Turbidez ***
## Temperatura ***
## Oxigenio ***
## pH ***
## Amonia ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Contrast: SEMI_INT
##
## average sd ratio ava avb cumsum p
## Turbidez 0.09842 0.05670 1.73580 20.94100 34.51000 0.782 0.473
## Temperatura 0.01276 0.00938 1.35970 28.20500 29.73000 0.883 0.999
## Oxigenio 0.00711 0.00496 1.43400 5.48200 4.58000 0.940 1.000
## Amonia 0.00411 0.00176 2.33290 0.50600 1.07000 0.972 0.005 **
## pH 0.00350 0.00260 1.34620 7.59500 7.88000 1.000 0.997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999# 3) Comparação específica:
simper_result$EXT_INT## $species
## [1] "Oxigenio" "pH" "Temperatura" "Amonia" "Turbidez"
##
## $average
## Oxigenio pH Temperatura Amonia Turbidez
## 0.018172812 0.007333667 0.029395534 0.006703702 0.180292849
##
## $overall
## [1] 0.2418986
##
## $sd
## Oxigenio pH Temperatura Amonia Turbidez
## 0.006919327 0.004192578 0.013267531 0.001866212 0.050003963
##
## $ratio
## Oxigenio pH Temperatura Amonia Turbidez
## 2.626384 1.749202 2.215599 3.592144 3.605571
##
## $ava
## Oxigenio pH Temperatura Amonia Turbidez
## 6.914287 6.952327 25.932920 0.212474 11.099315
##
## $avb
## Oxigenio pH Temperatura Amonia Turbidez
## 4.578261 7.878826 29.725737 1.072654 34.505555
##
## $ord
## [1] 5 3 1 2 4
##
## $cusum
## Turbidez Temperatura Oxigenio pH Amonia
## 0.7453242 0.8668443 0.9419700 0.9722871 1.0000000
##
## $p
## Oxigenio pH Temperatura Amonia Turbidez
## 0.001 0.001 0.001 0.001 0.001# 4) Organizando a tabela:
resultado <- summary(simper_result)
resultado$EXT_INT## average sd ratio ava avb
## Turbidez 0.180292849 0.050003963 3.605571 11.099315 34.505555
## Temperatura 0.029395534 0.013267531 2.215599 25.932920 29.725737
## Oxigenio 0.018172812 0.006919327 2.626384 6.914287 4.578261
## pH 0.007333667 0.004192578 1.749202 6.952327 7.878826
## Amonia 0.006703702 0.001866212 3.592144 0.212474 1.072654
## cumsum p
## Turbidez 0.7453242 0.001
## Temperatura 0.8668443 0.001
## Oxigenio 0.9419700 0.001
## pH 0.9722871 0.001
## Amonia 1.0000000 0.001# 5) Plot
# Organizar resultados da SIMPER
simper_df <- as.data.frame(simper_result$EXT_INT)
# Adicionar nomes das variáveis
simper_df$Variavel <- rownames(simper_df)
# Ordenar da maior para menor contribuição
simper_df <- simper_df[order(simper_df$average, decreasing = TRUE), ]
# ------------------------------
# Gráfico robusto
# ------------------------------
barplot(
simper_df$average,
names.arg = simper_df$Variavel,
# Deixar nomes na horizontal
las = 1,
# Melhor espaçamento
cex.names = 0.9,
# Limites do eixo Y
ylim = c(0, max(simper_df$average) * 1.2),
# Títulos
main = "Análise SIMPER - Piscicultura",
sub = "Contribuição das variáveis para a dissimilaridade entre sistemas",
xlab = "Variáveis ambientais",
ylab = "Contribuição média (%)",
# Bordas das barras
border = "black"
)
# Adicionar valores acima das barras
text(
x = seq_along(simper_df$average),
y = simper_df$average,
labels = round(simper_df$average, 2),
pos = 3,
cex = 0.8
)
=============================================================================
HOMOGENEIDADE MULTIVARIADA DE DISPERSÕES DE GRUPOS (BETADISPER)
# 1) Criando Datset:
set.seed(444)
sistema <- factor(rep(c("CONF", "SEMI", "PASTO"), each = 15))
# Produção diária (L/vaca)
Producao <- c(
rnorm(15, 32, 2),
rnorm(15, 26, 3),
rnorm(15, 20, 4)
)
# Gordura (%)
Gordura <- c(
rnorm(15, 3.4, 0.2),
rnorm(15, 3.7, 0.3),
rnorm(15, 4.1, 0.4)
)
# Proteína (%)
Proteina <- c(
rnorm(15, 3.1, 0.1),
rnorm(15, 3.3, 0.15),
rnorm(15, 3.5, 0.2)
)
# CCS (x1000 células/mL)
CCS <- c(
rnorm(15, 250, 40),
rnorm(15, 350, 60),
rnorm(15, 500, 90)
)
# Nitrogênio ureico no leite
NUL <- c(
rnorm(15, 14, 2),
rnorm(15, 16, 2),
rnorm(15, 18, 3)
)
dados <- data.frame(
Producao,
Gordura,
Proteina,
CCS,
NUL
)# 2) Instalando e carregando o pacote 'vegan':
install.packages("vegan")## Warning: o pacote 'vegan' está em uso e não será instaladolibrary(vegan)# 3) Matriz de distância: usaremos a Eucliana.
distancia <- dist(dados, method = "euclidean")# 4) Rodando a BETADISPER:
bd <- betadisper(distancia, sistema)
bd##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = distancia, group = sistema)
##
## No. of Positive Eigenvalues: 5
## No. of Negative Eigenvalues: 0
##
## Average distance to median:
## CONF PASTO SEMI
## 29.58 38.60 64.73
##
## Eigenvalues for PCoA axes:
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5
## 5.802e+05 7.565e+02 2.410e+02 5.526e+00 9.423e-01# 5) Testando a Homogeneidade:
anova(bd)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 10001 5000.6 2.8943 0.06644 .
## Residuals 42 72564 1727.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 6) Teste por permutação:
permutest(bd)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 10001 5000.6 2.8943 999 0.062 .
## Residuals 42 72564 1727.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 7) Comparações múltiplas:
TukeyHSD(bd)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = distances ~ group, data = df)
##
## $group
## diff lwr upr p adj
## PASTO-CONF 9.017152 -27.856885 45.89119 0.8240254
## SEMI-CONF 35.153848 -1.720189 72.02788 0.0644541
## SEMI-PASTO 26.136696 -10.737341 63.01073 0.2089781plot(bd)
# 8) Gráfico
boxplot(bd,
xlab = "Sistema de produção",
ylab = "Distância ao centróide",
main = "Homogeneidade multivariada das dispersões")
=============================================================================