Instalar paquetes y sus dependencias ——————
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 TRUE
Datos
?agridat
## starting httpd help server ... done
data("gomez.nitrogen")
?gomez.nitrogen
str(gomez.nitrogen)
## 'data.frame': 96 obs. of 4 variables:
## $ trt : Factor w/ 8 levels "T1","T2","T3",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ nitro: num 3.26 3.84 3.5 3.43 3.43 3.68 2.97 3.11 1.88 2.36 ...
## $ rep : Factor w/ 4 levels "R1","R2","R3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ stage: Factor w/ 3 levels "P1","P2","P3": 1 1 1 1 1 1 1 1 2 2 ...
df <- gomez.nitrogen
str(df)
## 'data.frame': 96 obs. of 4 variables:
## $ trt : Factor w/ 8 levels "T1","T2","T3",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ nitro: num 3.26 3.84 3.5 3.43 3.43 3.68 2.97 3.11 1.88 2.36 ...
## $ rep : Factor w/ 4 levels "R1","R2","R3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ stage: Factor w/ 3 levels "P1","P2","P3": 1 1 1 1 1 1 1 1 2 2 ...
df
## trt nitro rep stage
## 1 T1 3.26 R1 P1
## 2 T2 3.84 R1 P1
## 3 T3 3.50 R1 P1
## 4 T4 3.43 R1 P1
## 5 T5 3.43 R1 P1
## 6 T6 3.68 R1 P1
## 7 T7 2.97 R1 P1
## 8 T8 3.11 R1 P1
## 9 T1 1.88 R1 P2
## 10 T2 2.36 R1 P2
## 11 T3 2.20 R1 P2
## 12 T4 2.32 R1 P2
## 13 T5 1.98 R1 P2
## 14 T6 2.01 R1 P2
## 15 T7 2.66 R1 P2
## 16 T8 2.53 R1 P2
## 17 T1 1.40 R1 P3
## 18 T2 1.53 R1 P3
## 19 T3 1.33 R1 P3
## 20 T4 1.61 R1 P3
## 21 T5 1.11 R1 P3
## 22 T6 1.26 R1 P3
## 23 T7 1.87 R1 P3
## 24 T8 1.76 R1 P3
## 25 T1 2.98 R2 P1
## 26 T2 3.74 R2 P1
## 27 T3 3.49 R2 P1
## 28 T4 3.45 R2 P1
## 29 T5 3.24 R2 P1
## 30 T6 3.24 R2 P1
## 31 T7 2.90 R2 P1
## 32 T8 3.04 R2 P1
## 33 T1 1.74 R2 P2
## 34 T2 2.14 R2 P2
## 35 T3 2.28 R2 P2
## 36 T4 2.33 R2 P2
## 37 T5 1.70 R2 P2
## 38 T6 2.33 R2 P2
## 39 T7 2.74 R2 P2
## 40 T8 2.22 R2 P2
## 41 T1 1.24 R2 P3
## 42 T2 1.21 R2 P3
## 43 T3 1.54 R2 P3
## 44 T4 1.33 R2 P3
## 45 T5 1.25 R2 P3
## 46 T6 1.44 R2 P3
## 47 T7 1.81 R2 P3
## 48 T8 1.28 R2 P3
## 49 T1 2.78 R3 P1
## 50 T2 3.09 R3 P1
## 51 T3 3.03 R3 P1
## 52 T4 2.81 R3 P1
## 53 T5 3.45 R3 P1
## 54 T6 2.84 R3 P1
## 55 T7 2.92 R3 P1
## 56 T8 3.20 R3 P1
## 57 T1 1.76 R3 P2
## 58 T2 1.75 R3 P2
## 59 T3 2.48 R3 P2
## 60 T4 2.16 R3 P2
## 61 T5 1.78 R3 P2
## 62 T6 2.22 R3 P2
## 63 T7 2.67 R3 P2
## 64 T8 2.61 R3 P2
## 65 T1 1.44 R3 P3
## 66 T2 1.28 R3 P3
## 67 T3 1.46 R3 P3
## 68 T4 1.40 R3 P3
## 69 T5 1.39 R3 P3
## 70 T6 1.12 R3 P3
## 71 T7 1.31 R3 P3
## 72 T8 1.23 R3 P3
## 73 T1 2.77 R4 P1
## 74 T2 3.36 R4 P1
## 75 T3 3.36 R4 P1
## 76 T4 3.32 R4 P1
## 77 T5 3.09 R4 P1
## 78 T6 2.91 R4 P1
## 79 T7 2.42 R4 P1
## 80 T8 2.81 R4 P1
## 81 T1 2.00 R4 P2
## 82 T2 1.57 R4 P2
## 83 T3 2.47 R4 P2
## 84 T4 1.99 R4 P2
## 85 T5 1.74 R4 P2
## 86 T6 2.00 R4 P2
## 87 T7 2.98 R4 P2
## 88 T8 2.22 R4 P2
## 89 T1 1.25 R4 P3
## 90 T2 1.17 R4 P3
## 91 T3 1.41 R4 P3
## 92 T4 1.12 R4 P3
## 93 T5 1.20 R4 P3
## 94 T6 1.24 R4 P3
## 95 T7 1.56 R4 P3
## 96 T8 1.29 R4 P3
Verificar Valores atípicos
df%>%
select_if(is.numeric)%>%
outlier()
## nitro
## 3.84
ANOVA bifactorial para un DCA ———-
Realizar el análisis de varianza (ANOVA)
modelo_anova <- aov(nitro ~ trt + stage + trt:stage,
data = df)
summary(modelo_anova)
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 7 1.27 0.181 3.875 0.00124 **
## stage 2 52.04 26.021 557.602 < 2e-16 ***
## trt:stage 14 3.57 0.255 5.459 5.62e-07 ***
## Residuals 72 3.36 0.047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Homogeneidad de varianzas ——————-
Calcular los residuos del modelo ANOVA
residuos <- residuals(modelo_anova)
Agrupar los residuos por los niveles de las interacciones
grupo <- paste(df$trt, df$nitro, sep = "_")
grupos_residuos <- split(residuos, grupo)
Realizar la 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 = NaN, df = 94, p-value = NA
Gráfico
residuos_estandarizados <- rstandard(modelo_anova)
Obtener los valores ajustados (fitted values)
valores_ajustados <- fitted(modelo_anova)
Crear un gráfico de residuos estandarizados vs. valores
ajustados
ggplot
## function (data = NULL, mapping = aes(), ..., environment = parent.frame())
## {
## UseMethod("ggplot")
## }
## <bytecode: 0x000001ff857e37d8>
## <environment: namespace:ggplot2>
plot(valores_ajustados, residuos_estandarizados, pch = 11 ,
col = "blue",
xlab = "Valores Ajustados", ylab = "Residuos Estandarizados",
main = "Gráfico de Residuos Est. vs. Valores Ajustados")
plot(valores_ajustados, residuos_estandarizados, pch = 11 ,
col = "blue",
xlab = "Valores Ajustados", ylab = "Residuos Estandarizados",
main = "Gráfico de Residuos Est. vs. Valores Ajustados")
abline(h = 0, lty = 2, col = "red")
Obtener lsmeans (medias ajustadas) para la interacción
EVALUAR LOS EFECTOS SIMPLES COMO “Tratamiento”
LSD.test(modelo_anova,c("trt"),
console=TRUE)
##
## Study: modelo_anova ~ c("trt")
##
## LSD t Test for nitro
##
## Mean Square Error: 0.04666667
##
## trt, means and individual ( 95 %) CI
##
## nitro std r se LCL UCL Min Max Q25 Q50
## T1 2.041667 0.7186962 12 0.06236096 1.917352 2.165981 1.24 3.26 1.4300 1.820
## T2 2.253333 1.0058767 12 0.06236096 2.129019 2.377648 1.17 3.84 1.4675 1.945
## T3 2.379167 0.8271139 12 0.06236096 2.254852 2.503481 1.33 3.50 1.5200 2.375
## T4 2.272500 0.8326969 12 0.06236096 2.148186 2.396814 1.12 3.45 1.5575 2.240
## T5 2.113333 0.9190938 12 0.06236096 1.989019 2.237648 1.11 3.45 1.3550 1.760
## T6 2.190833 0.8435473 12 0.06236096 2.066519 2.315148 1.12 3.68 1.3950 2.115
## T7 2.400833 0.6000675 12 0.06236096 2.276519 2.525148 1.31 2.98 1.8550 2.665
## T8 2.275000 0.7339247 12 0.06236096 2.150686 2.399314 1.23 3.20 1.6425 2.375
## Q75
## T1 2.7725
## T2 3.1575
## T3 3.1125
## T4 2.9375
## T5 3.1275
## T6 2.8575
## T7 2.9050
## T8 2.8675
##
## Alpha: 0.05 ; DF Error: 72
## Critical Value of t: 1.993464
##
## least Significant Difference: 0.175807
##
## Treatments with the same letter are not significantly different.
##
## nitro groups
## T7 2.400833 a
## T3 2.379167 a
## T8 2.275000 ab
## T4 2.272500 ab
## T2 2.253333 ab
## T6 2.190833 bc
## T5 2.113333 bc
## T1 2.041667 c
EVALUAR LOS EFECTOS SIMPLES COMO “Stage”
LSD.test(modelo_anova,c("stage"),
console=TRUE)
##
## Study: modelo_anova ~ c("stage")
##
## LSD t Test for nitro
##
## Mean Square Error: 0.04666667
##
## stage, means and individual ( 95 %) CI
##
## nitro std r se LCL UCL Min Max Q25 Q50
## P1 3.170625 0.3225597 32 0.03818813 3.094498 3.246752 2.42 3.84 2.9175 3.155
## P2 2.181875 0.3496767 32 0.03818813 2.105748 2.258002 1.57 2.98 1.9550 2.210
## P3 1.370000 0.1948035 32 0.03818813 1.293873 1.446127 1.11 1.87 1.2400 1.320
## Q75
## P1 3.4300
## P2 2.3875
## P3 1.4450
##
## Alpha: 0.05 ; DF Error: 72
## Critical Value of t: 1.993464
##
## least Significant Difference: 0.1076593
##
## Treatments with the same letter are not significantly different.
##
## nitro groups
## P1 3.170625 a
## P2 2.181875 b
## P3 1.370000 c
EVALUAR LOS EFECTOS COMPUESTOS PORQUE DIERON DIF. SIGNIF. EN LA
INTERACCIÓN
resultados_lsd<-LSD.test(modelo_anova,c("trt","stage"),
console=T)
##
## Study: modelo_anova ~ c("trt", "stage")
##
## LSD t Test for nitro
##
## Mean Square Error: 0.04666667
##
## trt:stage, means and individual ( 95 %) CI
##
## nitro std r se LCL UCL Min Max Q25 Q50
## T1:P1 2.9475 0.22969182 4 0.1080123 2.732181 3.162819 2.77 3.26 2.7775 2.880
## T1:P2 1.8450 0.12041595 4 0.1080123 1.629681 2.060319 1.74 2.00 1.7550 1.820
## T1:P3 1.3325 0.10242884 4 0.1080123 1.117181 1.547819 1.24 1.44 1.2475 1.325
## T2:P1 3.5075 0.34673477 4 0.1080123 3.292181 3.722819 3.09 3.84 3.2925 3.550
## T2:P2 1.9550 0.35986108 4 0.1080123 1.739681 2.170319 1.57 2.36 1.7050 1.945
## T2:P3 1.2975 0.16152915 4 0.1080123 1.082181 1.512819 1.17 1.53 1.2000 1.245
## T3:P1 3.3450 0.21946906 4 0.1080123 3.129681 3.560319 3.03 3.50 3.2775 3.425
## T3:P2 2.3575 0.13961256 4 0.1080123 2.142181 2.572819 2.20 2.48 2.2600 2.375
## T3:P3 1.4350 0.08812869 4 0.1080123 1.219681 1.650319 1.33 1.54 1.3900 1.435
## T4:P1 3.2525 0.30048572 4 0.1080123 3.037181 3.467819 2.81 3.45 3.1925 3.375
## T4:P2 2.2000 0.16020820 4 0.1080123 1.984681 2.415319 1.99 2.33 2.1175 2.240
## T4:P3 1.3650 0.20207259 4 0.1080123 1.149681 1.580319 1.12 1.61 1.2775 1.365
## T5:P1 3.3025 0.17036725 4 0.1080123 3.087181 3.517819 3.09 3.45 3.2025 3.335
## T5:P2 1.8000 0.12436505 4 0.1080123 1.584681 2.015319 1.70 1.98 1.7300 1.760
## T5:P3 1.2375 0.11701140 4 0.1080123 1.022181 1.452819 1.11 1.39 1.1775 1.225
## T6:P1 3.1675 0.38361222 4 0.1080123 2.952181 3.382819 2.84 3.68 2.8925 3.075
## T6:P2 2.1400 0.16227549 4 0.1080123 1.924681 2.355319 2.00 2.33 2.0075 2.115
## T6:P3 1.2650 0.13203535 4 0.1080123 1.049681 1.480319 1.12 1.44 1.2100 1.250
## T7:P1 2.8025 0.25669372 4 0.1080123 2.587181 3.017819 2.42 2.97 2.7800 2.910
## T7:P2 2.7625 0.14930394 4 0.1080123 2.547181 2.977819 2.66 2.98 2.6675 2.705
## T7:P3 1.6375 0.25630386 4 0.1080123 1.422181 1.852819 1.31 1.87 1.4975 1.685
## T8:P1 3.0400 0.16673332 4 0.1080123 2.824681 3.255319 2.81 3.20 2.9825 3.075
## T8:P2 2.3950 0.20469489 4 0.1080123 2.179681 2.610319 2.22 2.61 2.2200 2.375
## T8:P3 1.3900 0.24805913 4 0.1080123 1.174681 1.605319 1.23 1.76 1.2675 1.285
## Q75
## T1:P1 3.0500
## T1:P2 1.9100
## T1:P3 1.4100
## T2:P1 3.7650
## T2:P2 2.1950
## T2:P3 1.3425
## T3:P1 3.4925
## T3:P2 2.4725
## T3:P3 1.4800
## T4:P1 3.4350
## T4:P2 2.3225
## T4:P3 1.4525
## T5:P1 3.4350
## T5:P2 1.8300
## T5:P3 1.2850
## T6:P1 3.3500
## T6:P2 2.2475
## T6:P3 1.3050
## T7:P1 2.9325
## T7:P2 2.8000
## T7:P3 1.8250
## T8:P1 3.1325
## T8:P2 2.5500
## T8:P3 1.4075
##
## Alpha: 0.05 ; DF Error: 72
## Critical Value of t: 1.993464
##
## least Significant Difference: 0.3045066
##
## Treatments with the same letter are not significantly different.
##
## nitro groups
## T2:P1 3.5075 a
## T3:P1 3.3450 ab
## T5:P1 3.3025 abc
## T4:P1 3.2525 abc
## T6:P1 3.1675 bcd
## T8:P1 3.0400 cde
## T1:P1 2.9475 de
## T7:P1 2.8025 e
## T7:P2 2.7625 e
## T8:P2 2.3950 f
## T3:P2 2.3575 f
## T4:P2 2.2000 fg
## T6:P2 2.1400 fgh
## T2:P2 1.9550 ghi
## T1:P2 1.8450 hij
## T5:P2 1.8000 ij
## T7:P3 1.6375 jk
## T3:P3 1.4350 kl
## T8:P3 1.3900 kl
## T4:P3 1.3650 kl
## T1:P3 1.3325 l
## T2:P3 1.2975 l
## T6:P3 1.2650 l
## T5:P3 1.2375 l
Crear un data frame con los resultados de la prueba LSD
resultados_lsd_df <- data.frame(
trt_stage = rownames(resultados_lsd$groups),
nitro = resultados_lsd$groups$nitro,
groups = resultados_lsd$groups$groups
)
Separar la columna trt_stage en trt y stage
resultados_lsd_df <- separate(resultados_lsd_df,
trt_stage, into = c("trt", "stage"),
sep = ":")
resultados_lsd_df$trt<-
factor(resultados_lsd_df$trt,
levels = c("T1",
"T2", "T3", "T4","T5","T6","T7","T8"))
resultados_lsd_df
## trt stage nitro groups
## 1 T2 P1 3.5075 a
## 2 T3 P1 3.3450 ab
## 3 T5 P1 3.3025 abc
## 4 T4 P1 3.2525 abc
## 5 T6 P1 3.1675 bcd
## 6 T8 P1 3.0400 cde
## 7 T1 P1 2.9475 de
## 8 T7 P1 2.8025 e
## 9 T7 P2 2.7625 e
## 10 T8 P2 2.3950 f
## 11 T3 P2 2.3575 f
## 12 T4 P2 2.2000 fg
## 13 T6 P2 2.1400 fgh
## 14 T2 P2 1.9550 ghi
## 15 T1 P2 1.8450 hij
## 16 T5 P2 1.8000 ij
## 17 T7 P3 1.6375 jk
## 18 T3 P3 1.4350 kl
## 19 T8 P3 1.3900 kl
## 20 T4 P3 1.3650 kl
## 21 T1 P3 1.3325 l
## 22 T2 P3 1.2975 l
## 23 T6 P3 1.2650 l
## 24 T5 P3 1.2375 l
Crear un gráfico de barras agrupadas
library(ggplot2)
ggplot(resultados_lsd_df, aes(x = trt, y = nitro, fill = stage, label = groups)) +
geom_bar(position = "dodge", stat = "identity") +
labs(title = "Comparación de Medias Significativas",
x = "Tratamiento", y = "Cant. de Nitrogeno") +
scale_fill_manual(values = c("P1" = "blue", "P2" = "green", "P3" = "yellow")) +
theme_minimal() +
theme(legend.position = "top")
Visualizar todas las paletas disponibles
display.brewer.all()
Crear un gráfico de barras agrupadas con paleta pastel
ggplot(resultados_lsd_df, aes(x = trt, y = nitro, fill = stage, label = groups)) +
geom_bar(position = "dodge", stat = "identity") +
labs(title = "Comparación de Medias Significativas",
x = "Tratamiento", y = "Cant, de Nitrógeno") +
scale_fill_brewer(palette = "Oranges") +
theme_minimal() +
theme(legend.position = "top")
Colocando los nombres
Crear un gráfico de barras agrupadas con paleta pastel
Aumentar el tamaño de letra
library(ggplot2)
ggplot(resultados_lsd_df, aes(x = trt, y = nitro, fill = stage, label = groups)) +
geom_bar(position = "dodge", stat = "identity") +
geom_text(position = position_dodge(width = 0.9),
vjust = -0.5,
aes(group = stage),
size = 4) +
labs(title = "Comparación de Medias Significativas",
x = "Tratamiento", y = "Cant, de Nitrógeno") +
scale_fill_brewer(palette = "Oranges") +
theme_minimal() +
theme(legend.position = "top")
