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")
#install.packages("emmeans")

library(readxl)
library(tidyverse)
library(skimr)
library(summarytools)
library(car)
library(agricolae)
library(emmeans)

2 La base de datos

Utilizaremos el dataset “VACAS” disponible en Agrarias Virtual.

# Con el argumento "sheet" le indicamos la hoja que debe cargar
VACUNAS<- read_excel("VACAS.xlsx", sheet = "tambo")
VACUNAS

## # A tibble: 30 × 4
##    TAMBO SUPLEMENTO REPETICION PROTEINA_BRUTA
##    <chr> <chr>           <dbl>          <dbl>
##  1 I     control             1           3.19
##  2 II    control             2           3.16
##  3 III   control             3           3.25
##  4 IV    control             4           3.48
##  5 V     control             5           3.25
##  6 VI    control             6           3.1 
##  7 I     S1                  1           3.03
##  8 II    S1                  2           3.07
##  9 III   S1                  3           3.23
## 10 IV    S1                  4           3.3 
## # ℹ 20 more rows

# Usamos la función skim del paquete "skimr"
str(VACUNAS)

## tibble [30 × 4] (S3: tbl_df/tbl/data.frame)
##  $ TAMBO         : chr [1:30] "I" "II" "III" "IV" ...
##  $ SUPLEMENTO    : chr [1:30] "control" "control" "control" "control" ...
##  $ REPETICION    : num [1:30] 1 2 3 4 5 6 1 2 3 4 ...
##  $ PROTEINA_BRUTA: num [1:30] 3.19 3.16 3.25 3.48 3.25 3.1 3.03 3.07 3.23 3.3 ...

skim(VACUNAS)

Data summary
Name	VACUNAS
Number of rows	30
Number of columns	4
_______________________
Column type frequency:
character	2
numeric	2
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
TAMBO	0	1	1	3	0	6	0
SUPLEMENTO	0	1	2	7	0	5	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
REPETICION	0	1	3.50	1.74	1.00	2.00	3.50	5.00	6.00	▇▃▃▃▃
PROTEINA_BRUTA	0	1	3.25	0.15	2.93	3.17	3.25	3.35	3.54	▂▃▇▅▃

Responder - Objetivo del estudio - ¿Cuál es la unidad experimental? - ¿Cuál es el factor en estudio y cuáles son sus niveles? - ¿Cuántas repeticiones tiene cada tratamiento?

3 Transformaciones de variables

Tambo -> factor
Suplemento -> factor
Repeticiones -> factor (es un identificador)
Proteína bruta -> numérica

VACUNAS <- VACUNAS %>% 
  mutate(
    TAMBO = factor(TAMBO),
    SUPLEMENTO = factor(SUPLEMENTO),
    REPETICION = factor(REPETICION)
  )

# COMPROBAMOS
str(VACUNAS)

## tibble [30 × 4] (S3: tbl_df/tbl/data.frame)
##  $ TAMBO         : Factor w/ 6 levels "I","II","III",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ SUPLEMENTO    : Factor w/ 5 levels "control","S1",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ REPETICION    : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ PROTEINA_BRUTA: num [1:30] 3.19 3.16 3.25 3.48 3.25 3.1 3.03 3.07 3.23 3.3 ...

skim(VACUNAS)

Data summary
Name	VACUNAS
Number of rows	30
Number of columns	4
_______________________
Column type frequency:
factor	3
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
TAMBO	1	FALSE	6	I: 5, II: 5, III: 5, IV: 5
SUPLEMENTO	1	FALSE	5	con: 6, S1: 6, S2: 6, S3: 6
REPETICION	1	FALSE	6	1: 5, 2: 5, 3: 5, 4: 5

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
PROTEINA_BRUTA	0	1	3.25	0.15	2.93	3.17	3.25	3.35	3.54	▂▃▇▅▃

4 Exploración de los datos

Medidas de resumen

# Opción 2: usando summarytools
VACUNAS %>% 
  group_by((SUPLEMENTO)) %>% 
  descr(PROTEINA_BRUTA)

## Descriptive Statistics  
## PROTEINA_BRUTA by (SUPLEMENTO)  
## Data Frame: VACUNAS  
## N: 30  
## 
##                     control       S1       S2       S3       S4
## ----------------- --------- -------- -------- -------- --------
##              Mean      3.24     3.15     3.15     3.34     3.38
##           Std.Dev      0.13     0.12     0.15     0.09     0.12
##               Min      3.10     3.03     2.93     3.22     3.20
##                Q1      3.16     3.05     3.06     3.28     3.33
##            Median      3.22     3.15     3.16     3.32     3.37
##                Q3      3.25     3.25     3.24     3.44     3.45
##               Max      3.48     3.30     3.33     3.45     3.54
##               MAD      0.07     0.15     0.13     0.10     0.09
##               IQR      0.08     0.19     0.18     0.14     0.10
##                CV      0.04     0.04     0.05     0.03     0.03
##          Skewness      0.80     0.06    -0.17     0.12    -0.10
##       SE.Skewness      0.85     0.85     0.85     0.85     0.85
##          Kurtosis     -0.86    -2.16    -1.77    -1.94    -1.37
##           N.Valid      6.00     6.00     6.00     6.00     6.00
##                 N      6.00     6.00     6.00     6.00     6.00
##         Pct.Valid    100.00   100.00   100.00   100.00   100.00

Visualización de datos

ggplot(VACUNAS, aes(x = SUPLEMENTO, y = PROTEINA_BRUTA, fill = SUPLEMENTO)) +
  geom_boxplot() +
  labs(title = "Proteína bruta en leche según suplemento",
       x = "Suplemento", y = "Proteína bruta (%)") +
  theme_minimal()

5 Análisis de la Varianza

Ajustamos un modelo ANAVA para evaluar diferencias entre tratamientos en un DBCA.
Vamos a usar la función aov() que es parte del paquete básico de R.
Función con argumentos:

aov(VARIABLE_RESPUESTA ~ FACTOR + BLOQUE, datos = BASE_DE_DATOS)

# Creamos el objeto modelo_dca
modelo_dbca <- aov(PROTEINA_BRUTA ~ SUPLEMENTO + TAMBO, data = VACUNAS)
summary(modelo_dbca)

##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## SUPLEMENTO   4 0.26035 0.06509   14.72 9.30e-06 ***
## TAMBO        5 0.28739 0.05748   13.00 1.03e-05 ***
## Residuals   20 0.08845 0.00442                     
## ---
## 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_dbca, which = 2)

6.1.2 Prueba de Shapiro-Wilks

H₀: SI se cumple el supuesto de Normalidad.

H₁: NO se cumple el supuesto de Normalidad.

shapiro.test(modelo_dbca$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_dbca$residuals
## W = 0.97055, p-value = 0.5544

6.2 Homocedasticidad

6.2.1 Gráfico de residuos vs predichos

plot(modelo_dbca, which = 1) # Residos vs predichos

6.2.2 Test de Levene

H₀: SI cumple el supuesto de Homocedasticidad.

H₁: NO se cumple el supuesto de Homocedasticidad

# Para esta prueba usamos una función de la librería "car"
leveneTest(PROTEINA_BRUTA ~ SUPLEMENTO, data = VACUNAS)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  0.5269 0.7169
##       25

6.3 Aditividad ( No interacción Bloque - Tratamiento)

# Usamos el paquete emmeans
emmip(modelo_dbca, TAMBO ~ SUPLEMENTO)

7 Pruebas post hoc

7.1 Test de tukey

# Usaremos una función del paquete agricolae
TUKEY_DBCA <- HSD.test(modelo_dbca, "SUPLEMENTO")
TUKEY_DBCA

## $statistics
##       MSerror Df     Mean       CV       MSD
##   0.004422333 20 3.250667 2.045753 0.1148897
## 
## $parameters
##    test     name.t ntr StudentizedRange alpha
##   Tukey SUPLEMENTO   5         4.231857  0.05
## 
## $means
##         PROTEINA_BRUTA        std r         se  Min  Max    Q25   Q50    Q75
## control       3.238333 0.13136463 6 0.02714877 3.10 3.48 3.1675 3.220 3.2500
## S1            3.155000 0.11794066 6 0.02714877 3.03 3.30 3.0550 3.150 3.2450
## S2            3.146667 0.14827902 6 0.02714877 2.93 3.33 3.0650 3.160 3.2400
## S3            3.336667 0.09352362 6 0.02714877 3.22 3.45 3.2800 3.315 3.4175
## S4            3.376667 0.11518102 6 0.02714877 3.20 3.54 3.3350 3.370 3.4350
## 
## $comparison
## NULL
## 
## $groups
##         PROTEINA_BRUTA groups
## S4            3.376667      a
## S3            3.336667     ab
## control       3.238333     bc
## S1            3.155000      c
## S2            3.146667      c
## 
## attr(,"class")
## [1] "group"

plot(TUKEY_DBCA)

7.2 Test de Fisher (LSD)

# Usamos librería agricolae
FISHER_DBCA <- LSD.test(modelo_dbca, "SUPLEMENTO")  
FISHER_DBCA

## $statistics
##       MSerror Df     Mean       CV  t.value       LSD
##   0.004422333 20 3.250667 2.045753 2.085963 0.0800888
## 
## $parameters
##         test p.ajusted     name.t ntr alpha
##   Fisher-LSD      none SUPLEMENTO   5  0.05
## 
## $means
##         PROTEINA_BRUTA        std r         se      LCL      UCL  Min  Max
## control       3.238333 0.13136463 6 0.02714877 3.181702 3.294965 3.10 3.48
## S1            3.155000 0.11794066 6 0.02714877 3.098369 3.211631 3.03 3.30
## S2            3.146667 0.14827902 6 0.02714877 3.090035 3.203298 2.93 3.33
## S3            3.336667 0.09352362 6 0.02714877 3.280035 3.393298 3.22 3.45
## S4            3.376667 0.11518102 6 0.02714877 3.320035 3.433298 3.20 3.54
##            Q25   Q50    Q75
## control 3.1675 3.220 3.2500
## S1      3.0550 3.150 3.2450
## S2      3.0650 3.160 3.2400
## S3      3.2800 3.315 3.4175
## S4      3.3350 3.370 3.4350
## 
## $comparison
## NULL
## 
## $groups
##         PROTEINA_BRUTA groups
## S4            3.376667      a
## S3            3.336667      a
## control       3.238333      b
## S1            3.155000      c
## S2            3.146667      c
## 
## attr(,"class")
## [1] "group"

plot(FISHER_DBCA)

8 Asignación de tratamientos con el paquete `agricolae`

Vamos a asignar los tratamientos seleccionados a las U.E. o parcelas usando el paquete agricolae de R.

Recordemos: en un DBCA los tratamientos se asignan al azar dentro de cada bloque, asegurando que en cada bloque estén presentes todos los tratamientos.

# Vector de tratamientos (5 tratamientos)
tratamientos <- c("control", "S1", "S2", "S3", "S4")

# Número de repeticiones = número de bloques (4 bloques)
tambos <- 6

# Asignación aleatoria con agricolae (una sola vez)
# design.rcbd: Randomized Complete Block Design
diseno_dbca <- design.rcbd(tratamientos, r = tambos, serie = 1, seed = 123)

# Mostrar la tabla ("field book") con la asignación aleatoria de tratamientos
diseno_dbca$book

##    plots block tratamientos
## 1     11     1      control
## 2     12     1           S2
## 3     13     1           S3
## 4     14     1           S4
## 5     15     1           S1
## 6     21     2           S1
## 7     22     2      control
## 8     23     2           S4
## 9     24     2           S2
## 10    25     2           S3
## 11    31     3      control
## 12    32     3           S3
## 13    33     3           S1
## 14    34     3           S4
## 15    35     3           S2
## 16    41     4           S4
## 17    42     4           S2
## 18    43     4           S3
## 19    44     4           S1
## 20    45     4      control
## 21    51     5           S3
## 22    52     5      control
## 23    53     5           S2
## 24    54     5           S4
## 25    55     5           S1
## 26    61     6           S2
## 27    62     6           S1
## 28    63     6           S3
## 29    64     6      control
## 30    65     6           S4

# plots = parcelas experimentales
# block = el bloque al que pertenece la unidad experimental

Argumentos

trat: número de tratamientosa realizar
rep_dca: número de repeticiones por tratamiento
Define el número de dígitos que tendrá el código de parcela o unidad experimental:
- serie = 0 → no agrega dígitos (parcela = 1, 2, 3 …).
- serie = 1 → agrega un dígito delante (parcela = 11, 12, 13 …).
- serie = 2 → agrega dos dígitos (parcela = 101, 102, …).
seed: Número de semilla aleatoria para reproducibilidad. Si se fija (ejemplo seed = 123), cada vez que corra la función se obtendrá la misma aleatorización.
randomization = TRUE: Indica que los tratamientos se asignan aleatoriamente a las unidades experimentales.

Si deseo que sea totalmente aleatoria…

dbca_random <- design.rcbd(tratamientos, r = tambos, serie = 1, randomization = TRUE) 
dbca_random$book

##    plots block tratamientos
## 1     11     1           S3
## 2     12     1      control
## 3     13     1           S2
## 4     14     1           S1
## 5     15     1           S4
## 6     21     2           S2
## 7     22     2           S1
## 8     23     2      control
## 9     24     2           S4
## 10    25     2           S3
## 11    31     3           S3
## 12    32     3      control
## 13    33     3           S1
## 14    34     3           S4
## 15    35     3           S2
## 16    41     4           S4
## 17    42     4           S2
## 18    43     4           S1
## 19    44     4           S3
## 20    45     4      control
## 21    51     5           S3
## 22    52     5           S4
## 23    53     5      control
## 24    54     5           S1
## 25    55     5           S2
## 26    61     6           S4
## 27    62     6      control
## 28    63     6           S2
## 29    64     6           S1
## 30    65     6           S3

dbca_random$sketch # Esquema de distribución del ensayo

##      [,1] [,2]      [,3]      [,4] [,5]     
## [1,] "S3" "control" "S2"      "S1" "S4"     
## [2,] "S2" "S1"      "control" "S4" "S3"     
## [3,] "S3" "control" "S1"      "S4" "S2"     
## [4,] "S4" "S2"      "S1"      "S3" "control"
## [5,] "S3" "S4"      "control" "S1" "S2"     
## [6,] "S4" "control" "S2"      "S1" "S3"

# en las filas del esquema se encuentran los bloques
# en las columnas del esquema los tratamientos

Diseño en Bloques Completamente al Azar

contacto: sbustos@agrarias.unca.edu.ar

Facultad de Ciencias Agrarias - Universidad Nacional de Catamarca

1 Instalamos y cargamos paquetes

2 La base de datos

3 Transformaciones de variables

4 Exploración de los datos

5 Análisis de la Varianza

6 Supuestos del modelo

6.1 Normalidad

6.1.1 Gráfico Q-Q plot

6.1.2 Prueba de Shapiro-Wilks

6.2 Homocedasticidad

6.2.1 Gráfico de residuos vs predichos

6.2.2 Test de Levene

6.3 Aditividad ( No interacción Bloque - Tratamiento)

7 Pruebas post hoc

7.1 Test de tukey

7.2 Test de Fisher (LSD)

8 Asignación de tratamientos con el paquete `agricolae`

Diseño en Bloques Completamente al Azar

contacto: sbustos@agrarias.unca.edu.ar

Facultad de Ciencias Agrarias - Universidad Nacional de Catamarca

1 Instalamos y cargamos paquetes

2 La base de datos

3 Transformaciones de variables

4 Exploración de los datos

5 Análisis de la Varianza

6 Supuestos del modelo

6.1 Normalidad

6.1.1 Gráfico Q-Q plot

6.1.2 Prueba de Shapiro-Wilks

6.2 Homocedasticidad

6.2.1 Gráfico de residuos vs predichos

6.2.2 Test de Levene

6.3 Aditividad ( No interacción Bloque - Tratamiento)

7 Pruebas post hoc

7.1 Test de tukey

7.2 Test de Fisher (LSD)

8 Asignación de tratamientos con el paquete agricolae

8 Asignación de tratamientos con el paquete `agricolae`