PARCIAL
Luis Fernando Rosero Males
2023-12-14
Crear un vector de paquetes
pacotes <- c("agricolae", "performance", "ggplot2",
"outliers","car", "readxl","tidyverse",
"RColorBrewer","nortest", "stats",
"agridat", "performance")Script para instalar y cargar librerias y dependencias
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
instalador <- pacotes[!pacotes %in% installed.packages()]
for(i in 1:length(instalador)) {
install.packages(instalador, dependencies = T)
break()}
sapply(pacotes, require, character = T)
} else {
sapply(pacotes, require, character = T)
}## Loading required package: agricolae## Loading required package: performance## Loading required package: ggplot2## Loading required package: outliers## Loading required package: car## Loading required package: carData## Loading required package: readxl## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: RColorBrewer
##
## Loading required package: nortest
##
## Loading required package: agridat## Warning: package 'agridat' was built under R version 4.3.2## agricolae performance ggplot2 outliers car readxl
## TRUE TRUE TRUE TRUE TRUE TRUE
## tidyverse RColorBrewer nortest stats agridat performance
## TRUE TRUE TRUE TRUE TRUE TRUEDATOS
library(agricolae)
f <- system.file("external/ssp.csv", package="agricolae")
ssp<-read.csv(f)
view (ssp)
str(ssp)## 'data.frame': 135 obs. of 5 variables:
## $ block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ nitrogen : int 0 0 0 50 50 50 80 80 80 110 ...
## $ management: chr "m1" "m2" "m3" "m1" ...
## $ variety : int 1 1 1 1 1 1 1 1 1 1 ...
## $ yield : num 3.32 3.77 4.66 3.19 3.62 ...summary(ssp)## block nitrogen management variety yield
## Min. :1 Min. : 0 Length:135 Min. :1 Min. : 3.132
## 1st Qu.:1 1st Qu.: 50 Class :character 1st Qu.:1 1st Qu.: 5.301
## Median :2 Median : 80 Mode :character Median :2 Median : 6.509
## Mean :2 Mean : 76 Mean :2 Mean : 6.554
## 3rd Qu.:3 3rd Qu.:110 3rd Qu.:3 3rd Qu.: 7.644
## Max. :3 Max. :140 Max. :3 Max. :10.360Valores atípicos
ssp%>%
select_if(is.numeric)%>%
outlier()## block nitrogen variety yield
## 3.00 0.00 3.00 10.36##ANOVA DE 3 FACTORES PARA DCA: para lo cual se verifica la estructura de la base de datos
str (ssp)## 'data.frame': 135 obs. of 5 variables:
## $ block : int 1 1 1 1 1 1 1 1 1 1 ...
## $ nitrogen : int 0 0 0 50 50 50 80 80 80 110 ...
## $ management: chr "m1" "m2" "m3" "m1" ...
## $ variety : int 1 1 1 1 1 1 1 1 1 1 ...
## $ yield : num 3.32 3.77 4.66 3.19 3.62 ...ssp$block <- as.factor(ssp$block)
ssp$nitrogen <- as.factor(ssp$nitrogen)
ssp$variety <- as.factor(ssp$variety)
str (ssp)## 'data.frame': 135 obs. of 5 variables:
## $ block : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ nitrogen : Factor w/ 5 levels "0","50","80",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ management: chr "m1" "m2" "m3" "m1" ...
## $ variety : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ yield : num 3.32 3.77 4.66 3.19 3.62 ...model_DCA <- aov(yield ~ management * variety * nitrogen, data = ssp)
summary(model_DCA)## Df Sum Sq Mean Sq F value Pr(>F)
## management 2 42.94 21.47 48.120 6.13e-15 ***
## variety 2 206.01 103.01 230.886 < 2e-16 ***
## nitrogen 4 61.64 15.41 34.542 < 2e-16 ***
## management:variety 4 3.85 0.96 2.158 0.080061 .
## management:nitrogen 8 1.10 0.14 0.309 0.960824
## variety:nitrogen 8 14.14 1.77 3.963 0.000468 ***
## management:variety:nitrogen 16 3.70 0.23 0.518 0.931339
## Residuals 90 40.15 0.45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cv.model(model_DCA)## [1] 10.19059Se verifica la normalidad mediante una prueba de Shapiro-Wilk
shapiro.test(residuals(model_DCA))##
## Shapiro-Wilk normality test
##
## data: residuals(model_DCA)
## W = 0.99198, p-value = 0.6403Homogeneidad de varianzas ——————-
Calcular los residuos del modelo ANOVA
residuos <- residuals(model_DCA)Agrupar los residuos por los niveles de las interacciones
grupo <- paste(ssp$variety, ssp$nitrogen, sep = "_")
grupos_residuos <- split(residuos, grupo)Prueba de homogeneidad de varianzas con fligner.test
fligner.test(grupos_residuos)##
## Fligner-Killeen test of homogeneity of variances
##
## data: grupos_residuos
## Fligner-Killeen:med chi-squared = 17.786, df = 14, p-value = 0.2167Gráfico
residuos_estandarizados <- rstandard(model_DCA)Obtener los valores ajustados (fitted values)
valores_ajustados <- fitted(model_DCA)Graficar residuos estandarizados vs. valores ajustados
plot(valores_ajustados, residuos_estandarizados, pch = 14 ,
col = "black",
xlab = "Valores Ajustados", ylab = "Residuos Estandarizados",
main = "Gráfico de Residuos Est. vs. Valores Ajustados")
Otra forma de verificar los supuestos:
plot(model_DCA)
check_model(model_DCA)
EVALUAR LOS EFECTOS COMPUESTOS POR LA DIF. SIGNIF. EN LA INTERACCIÓN
resultados_lsd<-LSD.test(model_DCA,c("nitrogen","variety"),
console=T)##
## Study: model_DCA ~ c("nitrogen", "variety")
##
## LSD t Test for yield
##
## Mean Square Error: 0.4461352
##
## nitrogen:variety, means and individual ( 95 %) CI
##
## yield std r se LCL UCL Min Max Q25 Q50
## 0:1 4.513111 0.8176433 9 0.2226445 4.070789 4.955433 3.320 5.915 3.864 4.507
## 0:2 5.162889 0.7955654 9 0.2226445 4.720567 5.605211 4.166 6.573 4.815 5.096
## 0:3 6.478111 1.0853742 9 0.2226445 6.035789 6.920433 5.244 8.020 5.536 6.462
## 110:1 5.444667 0.7656558 9 0.2226445 5.002344 5.886989 4.246 6.829 4.863 5.345
## 110:2 6.924556 0.7410263 9 0.2226445 6.482233 7.366878 5.779 7.856 6.209 6.992
## 110:3 8.442889 0.9350460 9 0.2226445 8.000567 8.885211 6.902 9.660 7.778 8.966
## 140:1 5.077778 1.1663266 9 0.2226445 4.635456 5.520100 3.132 7.309 4.375 5.217
## 140:2 7.288444 0.7236999 9 0.2226445 6.846122 7.730767 6.573 8.950 6.860 6.974
## 140:3 9.335556 0.7091084 9 0.2226445 8.893233 9.777878 8.032 10.360 9.224 9.314
## 50:1 4.763667 0.8726879 9 0.2226445 4.321344 5.205989 3.188 6.046 4.752 4.809
## 50:2 6.016222 0.9620228 9 0.2226445 5.573900 6.458544 4.478 7.442 5.390 5.925
## 50:3 7.881111 1.1232814 9 0.2226445 7.438789 8.323433 6.546 9.942 7.092 7.646
## 80:1 5.834889 0.7497605 9 0.2226445 5.392567 6.277211 4.422 7.106 5.468 5.788
## 80:2 6.588556 0.7006203 9 0.2226445 6.146233 7.030878 5.442 7.991 6.398 6.533
## 80:3 8.563778 0.7680380 9 0.2226445 8.121456 9.006100 6.698 9.320 8.514 8.650
## Q75
## 0:1 4.875
## 0:2 5.495
## 0:3 7.442
## 110:1 5.869
## 110:2 7.565
## 110:3 9.080
## 140:1 5.389
## 140:2 7.422
## 140:3 9.712
## 50:1 5.232
## 50:2 6.780
## 50:3 8.592
## 80:1 6.215
## 80:2 6.914
## 80:3 9.112
##
## Alpha: 0.05 ; DF Error: 90
## Critical Value of t: 1.986675
##
## least Significant Difference: 0.625538
##
## Treatments with the same letter are not significantly different.
##
## yield groups
## 140:3 9.335556 a
## 80:3 8.563778 b
## 110:3 8.442889 bc
## 50:3 7.881111 cd
## 140:2 7.288444 de
## 110:2 6.924556 ef
## 80:2 6.588556 fg
## 0:3 6.478111 fg
## 50:2 6.016222 gh
## 80:1 5.834889 h
## 110:1 5.444667 hi
## 0:2 5.162889 ij
## 140:1 5.077778 ijk
## 50:1 4.763667 jk
## 0:1 4.513111 kCrear un data frame con los resultados de la prueba LSD
resultados_lsd_df <- data.frame(
nitrogen_variety = rownames(resultados_lsd$groups),
yield = resultados_lsd$groups$yield,
groups = resultados_lsd$groups$groups
)Separar la columna nitrogen_variety en nitrogen y variety
resultados_lsd_df <- separate(resultados_lsd_df, nitrogen_variety, into = c("nitrogen", "variety"), sep = ":")
resultados_lsd_df$nitrogen<-
factor(resultados_lsd_df$nitrogen,
levels = c("0",
"50", "80", "110", "140"))Crear un gráfico de barras agrupadas
ggplot(resultados_lsd_df, aes(x = variety, y = yield, fill = nitrogen, label = groups)) +
geom_bar(position = "dodge", stat = "identity") +
geom_text(position = position_dodge(width = 0.9),
vjust = -0.5,
aes(group = nitrogen),
size = 4) +
labs(title = "Comparación de Medias Significativas",
x = "Variedad", y = "Producción") +
scale_fill_brewer(palette = "Spectral") +
theme_minimal() +
theme(legend.position = "top") 
Crea un gráfico de dispersión tridimensional
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_ly(ssp, x = ~nitrogen, y = ~management, z = ~variety,
marker = list(color = ~yield, colorscale = "YlGnBu", size = 5),
type = "scatter3d") %>%
layout(scene = list(
xaxis = list(title = "Nitrogen"),
yaxis = list(title = "Management"),
zaxis = list(title = "Variety"),
bgcolor = "black" # Color oscuro para el fondo
))## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode