Vicente Mallama

2023-12-13

Instalar paquetes y sus dependencias

Crear un vector de paquetes

pacotes <- c("agricolae", "performance", "ggplot2",
             "outliers","car", "readxl","tidyverse", 
             "RColorBrewer","nortest", "stats", 
             "agridat", "performance")

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.2     ✔ 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

Importar datos

library(agricolae)
f <- system.file("external/ssp.csv", package="agricolae")
df<-read.csv(f)
head(df)
##   block nitrogen management variety yield
## 1     1        0         m1       1 3.320
## 2     1        0         m2       1 3.766
## 3     1        0         m3       1 4.660
## 4     1       50         m1       1 3.188
## 5     1       50         m2       1 3.625
## 6     1       50         m3       1 5.232

Convertir a factores

df$block <- factor(df$block)
df$nitrogen <- factor(df$nitrogen)
df$variety <- factor(df$variety)

Verificar los cambios

str(df)
## '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 ...

Valores atípicos

df%>%
  select_if(is.numeric)%>%
  outlier() #desconsideramos porque la base de datos es pequeña
## yield 
## 10.36

ANOVA Trifactorial para un DCA

Modelo de ANOVA

modelo_anova <- aov(yield ~ management * variety * nitrogen, data = df)

Resumen del modelo

summary(modelo_anova)
##                             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 ' ' 1

Normalidad PRUEBAS

shapiro.test(residuals(modelo_anova)) #shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_anova)
## W = 0.99198, p-value = 0.6403
# es mayor a 0,05, tiene una distribución normal
ad.test(residuals(modelo_anova)) #darling
## 
##  Anderson-Darling normality test
## 
## data:  residuals(modelo_anova)
## A = 0.30438, p-value = 0.5657
# es mayor a 0,05, tiene una distribución normal

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$nitrogen, df$management, df$variety, sep = "_")
grupos_residuos <- split(residuos, grupo)

Realizar la prueba de homogeneidad de varianzas con fligner.test

se hace para verificar los supuestos, si son homogeneas o heterogeneas, puedo ver si estoy cometiendo algun tipo de error.

bartlett.test(grupos_residuos) #acepto h0, hay homogeneidad en las varianzas
## 
##  Bartlett test of homogeneity of variances
## 
## data:  grupos_residuos
## Bartlett's K-squared = 49.016, df = 44, p-value = 0.2789
fligner.test(grupos_residuos) #acepto h0, hay homogeneidad en las varianzas
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  grupos_residuos
## Fligner-Killeen:med chi-squared = 18.798, df = 44, p-value = 0.9997

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

plot(valores_ajustados, residuos_estandarizados, pch = 12, 
     col = "blue",
     xlab = "Valores Ajustados", ylab = "Residuos Estandarizados",
     main = "Gráfico de Residuos Est. vs. Valores Ajustados")
abline(h = 0, lty = 2, col = "red")  # Línea de referencia en y = 0

Otra forma de verificar los supuestos:

plot(modelo_anova)

check_model(modelo_anova)

Evaluar los efectos compuestos porque dieron dif. Signif. En la interacción

resultados_lsd<-LSD.test(modelo_anova,c("nitrogen","variety"),
                         console=T)
## 
## Study: modelo_anova ~ 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      k
library(ggplot2)

Crear 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"))

Gráfico

Crear un gráfico de barras agrupadas con paleta pastel

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),    # Ajusta según el ancho de las barras
            vjust = -0.5,                              # Ajusta la posición vertical del texto
            aes(group = nitrogen),                      # Agrupa por DosisN para que se alinee correctamente
            size = 4) +                                # Aumenta el tamaño de la letra
  labs(title = "Comparación de Medias Significativas",
       x = "Variedad", y = "Rendimiento (tn/ha)") +
  scale_fill_brewer(palette = "Greens") +  # Paleta de colores pastel
  theme_minimal() +
  theme(legend.position = "top")  # Mover la leyenda a la parte superior

# Instala y carga la biblioteca plotly si no lo has hecho aún
# install.packages("plotly")
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':
## 
##     layout

Gráfico de dispersión tridimensional

plot_ly(df, x = ~nitrogen, y = ~yield, z = ~variety, 
        marker = list(color = ~yield, colorscale = "Greens", size = 5),
        type = "scatter3d") %>%
  layout(scene = list(xaxis = list(title = "Nitrogeno"),
                      yaxis = list(title = "Rendimiento"),
                      zaxis = list(title = "Variedad")))
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Gracias