1 Instalamos y cargamos paquetes

# Solo correr una vez si no estan instalados los siguientes paquetes:
# install.packages("readxl")
# install.packages("skimr")
# install.packages("car")
# install.packages("agricolae")
library(readxl)
library(tidyverse)
library(skimr)
library(car)
library(agricolae)

2 La base de datos

Utilizaremos el dataset “VACAS” disponible en Agrarias Virtual.

# Con el argumento "sheet" le indicamos la hoja que debe cargar
ALIMENTACION <- read_excel("VACAS.xlsx", sheet = "raciones")
ALIMENTACION
## # A tibble: 25 × 4
##    RACIONES  RAZA  EDAD  PESO
##    <chr>    <dbl> <dbl> <dbl>
##  1 R1           1     4    94
##  2 R1           2     2    70
##  3 R1           3     5    85
##  4 R1           4     1    90
##  5 R1           5     3    90
##  6 R2           1     3    70
##  7 R2           2     4    63
##  8 R2           3     2    60
##  9 R2           4     5    65
## 10 R2           5     1    64
## # ℹ 15 more rows
str(ALIMENTACION)
## tibble [25 × 4] (S3: tbl_df/tbl/data.frame)
##  $ RACIONES: chr [1:25] "R1" "R1" "R1" "R1" ...
##  $ RAZA    : num [1:25] 1 2 3 4 5 1 2 3 4 5 ...
##  $ EDAD    : num [1:25] 4 2 5 1 3 3 4 2 5 1 ...
##  $ PESO    : num [1:25] 94 70 85 90 90 70 63 60 65 64 ...

3 Transformación de variables

ALIMENTACION <- ALIMENTACION %>% 
  mutate(
    RACIONES = factor(RACIONES),
    RAZA = factor(RAZA),
    EDAD = factor(EDAD))
str(ALIMENTACION)
## tibble [25 × 4] (S3: tbl_df/tbl/data.frame)
##  $ RACIONES: Factor w/ 5 levels "R1","R2","R3",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ RAZA    : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ EDAD    : Factor w/ 5 levels "1","2","3","4",..: 4 2 5 1 3 3 4 2 5 1 ...
##  $ PESO    : num [1:25] 94 70 85 90 90 70 63 60 65 64 ...

4 Exploración de datos

ggplot(ALIMENTACION, aes(x = RACIONES, y = PESO, fill = RACIONES)) +
  geom_boxplot() +
  labs(title = "Peso de vacas según ración",
       x = "Ración", y = "Peso (kg)") +
  theme_minimal()

5 Análisis de la Varianza

  • Función con argumentos:

    aov(variable_respuesta ~ tratamiento + fila + columna, datos = base_de_datos)

modelo_dcl <- aov(PESO ~ RACIONES + RAZA + EDAD, data = ALIMENTACION)
summary(modelo_dcl)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## RACIONES     4 3154.8   788.7  60.825 7.24e-08 ***
## RAZA         4  348.8    87.2   6.725  0.00444 ** 
## EDAD         4  314.8    78.7   6.069  0.00657 ** 
## Residuals   12  155.6    13.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6 Supuestos del modelo

6.1 Normalidad

6.1.1 Gráfico Q-Q plot

plot(modelo_dcl, which = 2)

6.1.2 Prueba de Shapiro-Wilks

H0: SI se cumple el supuesto de Normalidad.

H1: NO se cumple el supuesto de Normalidad.

shapiro.test(modelo_dcl$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_dcl$residuals
## W = 0.96177, p-value = 0.4508

6.2 Homocedasticidad

6.2.1 Gráfico de residuos vs predichos

plot(modelo_dcl, which = 1)

6.2.2 Test de Levene

H0: SI cumple el supuesto de Homocedasticidad.

H1: NO se cumple el supuesto de Homocedasticidad

leveneTest(PESO ~ RACIONES, data = ALIMENTACION)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  1.0468 0.4081
##       20

7 Pruebas post hoc

7.1 Test de tukey

TUKEY_DCL <- HSD.test(modelo_dcl, "RACIONES")
TUKEY_DCL
## $statistics
##    MSerror Df Mean       CV      MSD
##   12.96667 12   68 5.295479 7.259139
## 
## $parameters
##    test   name.t ntr StudentizedRange alpha
##   Tukey RACIONES   5          4.50771  0.05
## 
## $means
##    PESO      std r       se Min Max Q25 Q50 Q75
## R1 85.8 9.391486 5 1.610383  70  94  85  90  90
## R2 64.4 3.646917 5 1.610383  60  70  63  64  65
## R3 58.8 4.969909 5 1.610383  53  66  56  58  61
## R4 75.6 8.294577 5 1.610383  67  85  68  75  83
## R5 55.4 3.130495 5 1.610383  50  58  56  56  57
## 
## $comparison
## NULL
## 
## $groups
##    PESO groups
## R1 85.8      a
## R4 75.6      b
## R2 64.4      c
## R3 58.8     cd
## R5 55.4      d
## 
## attr(,"class")
## [1] "group"
plot(TUKEY_DCL)

7.2 Test de Fisher (LSD)

FISHER_DCL <- LSD.test(modelo_dcl, "RACIONES")
FISHER_DCL
## $statistics
##    MSerror Df Mean       CV  t.value      LSD
##   12.96667 12   68 5.295479 2.178813 4.962084
## 
## $parameters
##         test p.ajusted   name.t ntr alpha
##   Fisher-LSD      none RACIONES   5  0.05
## 
## $means
##    PESO      std r       se      LCL      UCL Min Max Q25 Q50 Q75
## R1 85.8 9.391486 5 1.610383 82.29128 89.30872  70  94  85  90  90
## R2 64.4 3.646917 5 1.610383 60.89128 67.90872  60  70  63  64  65
## R3 58.8 4.969909 5 1.610383 55.29128 62.30872  53  66  56  58  61
## R4 75.6 8.294577 5 1.610383 72.09128 79.10872  67  85  68  75  83
## R5 55.4 3.130495 5 1.610383 51.89128 58.90872  50  58  56  56  57
## 
## $comparison
## NULL
## 
## $groups
##    PESO groups
## R1 85.8      a
## R4 75.6      b
## R2 64.4      c
## R3 58.8      d
## R5 55.4      d
## 
## attr(,"class")
## [1] "group"
plot(FISHER_DCL)

8 Asignación de tratamientos con el paquete agricolae

# Vector de tratamientos
tratamientos_dcl <- c("R1", "R2", "R3", "R4", "R5")

# Generar el cuadrado latino
diseno_dcl <- design.lsd(tratamientos_dcl, serie = 1, seed = 123)

# Mostrar asignación de tratamientos a las unidades experimentales
# plot = parcela o unidad experimental 
# fila = row = RAZA
# Columna = col = EDAD
diseno_dcl$book
##    plots row col tratamientos_dcl
## 1     11   1   1               R1
## 2     12   1   2               R4
## 3     13   1   3               R2
## 4     14   1   4               R5
## 5     15   1   5               R3
## 6     21   2   1               R4
## 7     22   2   2               R2
## 8     23   2   3               R5
## 9     24   2   4               R3
## 10    25   2   5               R1
## 11    31   3   1               R2
## 12    32   3   2               R5
## 13    33   3   3               R3
## 14    34   3   4               R1
## 15    35   3   5               R4
## 16    41   4   1               R5
## 17    42   4   2               R3
## 18    43   4   3               R1
## 19    44   4   4               R4
## 20    45   4   5               R2
## 21    51   5   1               R3
## 22    52   5   2               R1
## 23    53   5   3               R4
## 24    54   5   4               R2
## 25    55   5   5               R5
# Mostrar un esquema de campo
diseno_dcl$sketch
##      [,1] [,2] [,3] [,4] [,5]
## [1,] "R1" "R4" "R2" "R5" "R3"
## [2,] "R4" "R2" "R5" "R3" "R1"
## [3,] "R2" "R5" "R3" "R1" "R4"
## [4,] "R5" "R3" "R1" "R4" "R2"
## [5,] "R3" "R1" "R4" "R2" "R5"

Si deseo que sea totalmente aleatoria…

dcl_random <- design.lsd(tratamientos_dcl, serie = 2, randomization = TRUE) 
dcl_random$book
##    plots row col tratamientos_dcl
## 1    101   1   1               R1
## 2    102   1   2               R3
## 3    103   1   3               R2
## 4    104   1   4               R5
## 5    105   1   5               R4
## 6    201   2   1               R4
## 7    202   2   2               R1
## 8    203   2   3               R5
## 9    204   2   4               R3
## 10   205   2   5               R2
## 11   301   3   1               R3
## 12   302   3   2               R5
## 13   303   3   3               R4
## 14   304   3   4               R2
## 15   305   3   5               R1
## 16   401   4   1               R5
## 17   402   4   2               R2
## 18   403   4   3               R1
## 19   404   4   4               R4
## 20   405   4   5               R3
## 21   501   5   1               R2
## 22   502   5   2               R4
## 23   503   5   3               R3
## 24   504   5   4               R1
## 25   505   5   5               R5
dcl_random$sketch # Esquema de distribución del ensayo
##      [,1] [,2] [,3] [,4] [,5]
## [1,] "R1" "R3" "R2" "R5" "R4"
## [2,] "R4" "R1" "R5" "R3" "R2"
## [3,] "R3" "R5" "R4" "R2" "R1"
## [4,] "R5" "R2" "R1" "R4" "R3"
## [5,] "R2" "R4" "R3" "R1" "R5"
# en las filas del esquema se encuentran los bloques
# en las columnas del esquema los tratamientos