1 Instalamos y cargamos paquetes

# Solo correr una vez si no estan instalados los siguientes paquetes:

#install.packages("skimr")
# install.packages("car")
# install.packages("agricolae")
library(tidyverse)
library(skimr)
library(summarytools)
library(car)
library(agricolae)

2 La base de datos

Utilizaremos el dataset PlantGrowth disponible en R.

# Exploramos al base de datos
str(PlantGrowth)
## 'data.frame':    30 obs. of  2 variables:
##  $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
##  $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...
glimpse(PlantGrowth)
## Rows: 30
## Columns: 2
## $ weight <dbl> 4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14, 4.8…
## $ group  <fct> ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, trt…
# Usamos la función skim del paquete "skimr"
skim(PlantGrowth)
Data summary
Name PlantGrowth
Number of rows 30
Number of columns 2
_______________________
Column type frequency:
factor 1
numeric 1
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
group 0 1 FALSE 3 ctr: 10, trt: 10, trt: 10

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
weight 0 1 5.07 0.7 3.59 4.55 5.15 5.53 6.31 ▂▇▇▇▆

Interpretación:

  1. skim_variable: Nombre de la variable que se está describiendo. En este caso: group.

  2. n_missing: Número de valores faltantes (NA). Aquí es 0, es decir no hay datos perdidos.

  3. complete_rate: Proporción de valores completos (no faltantes). Un valor de 1 significa que el 100% de los datos están completos.

  4. ordered: Indica si el factor es ordenado (factor ordinal) o no. Aquí aparece FALSE, o sea que group es un factor nominal (no hay un orden natural entre los niveles).

  5. n_unique: Número de categorías únicas (niveles del factor). En este caso es 3, entonces group tiene tres niveles distintos.

  6. top_counts: Muestra los niveles más frecuentes y sus frecuencias. ctr: 10, trt1: 10, trt2: 10 significa que el nivel ctr aparece 10 veces, lo mismo que trt1 y trt2.


Variables presentes en la base de datos:

weight: peso seco de la planta (en gramos).

group: tratamientos aplicados (control, trt1, trt2)


Responder:

  • ¿Cuál es el objetivo del estudio?
?PlantGrowth
  • ¿Cuál es la unidad experimental?
  • ¿Cuál es el factor en estudio y cuáles son sus niveles?
  • ¿Qué representa la variable de respuesta?
  • ¿Cuántas repeticiones tiene cada tratamiento?

3 Exploración de los datos

  1. Medidas de resumen
# Opción 1: usando los paquetes del tidyverse
PlantGrowth %>%
  group_by(group) %>%
  summarise(
    media = mean(weight),
    sd = sd(weight),
    n = n())
## # A tibble: 3 × 4
##   group media    sd     n
##   <fct> <dbl> <dbl> <int>
## 1 ctrl   5.03 0.583    10
## 2 trt1   4.66 0.794    10
## 3 trt2   5.53 0.443    10
# Opción 2: usando summarytools
PlantGrowth %>% 
  group_by((group)) %>% 
  descr(weight)
## Descriptive Statistics  
## weight by (group)  
## Data Frame: PlantGrowth  
## N: 30  
## 
##                       ctrl     trt1     trt2
## ----------------- -------- -------- --------
##              Mean     5.03     4.66     5.53
##           Std.Dev     0.58     0.79     0.44
##               Min     4.17     3.59     4.92
##                Q1     4.53     4.17     5.26
##            Median     5.15     4.55     5.44
##                Q3     5.33     4.89     5.80
##               Max     6.11     6.03     6.31
##               MAD     0.72     0.53     0.36
##               IQR     0.74     0.66     0.47
##                CV     0.12     0.17     0.08
##          Skewness     0.23     0.47     0.48
##       SE.Skewness     0.69     0.69     0.69
##          Kurtosis    -1.12    -1.10    -1.16
##           N.Valid    10.00    10.00    10.00
##                 N    10.00    10.00    10.00
##         Pct.Valid   100.00   100.00   100.00
  1. Visualización de datos
ggplot(PlantGrowth, aes(x = group, y = weight, fill = group)) +
  geom_boxplot() +
  labs(title = "Peso de plantas según tratamiento",
       x = "Tratamientos", y = "Peso (g)") +
  theme_minimal()

4 Renombrar variables del dataset

bd_crecimiento_plantas <- PlantGrowth
bd_crecimiento_plantas <- bd_crecimiento_plantas %>% 
  rename(peso = weight,
         tratamiento = group)
bd_crecimiento_plantas 
##    peso tratamiento
## 1  4.17        ctrl
## 2  5.58        ctrl
## 3  5.18        ctrl
## 4  6.11        ctrl
## 5  4.50        ctrl
## 6  4.61        ctrl
## 7  5.17        ctrl
## 8  4.53        ctrl
## 9  5.33        ctrl
## 10 5.14        ctrl
## 11 4.81        trt1
## 12 4.17        trt1
## 13 4.41        trt1
## 14 3.59        trt1
## 15 5.87        trt1
## 16 3.83        trt1
## 17 6.03        trt1
## 18 4.89        trt1
## 19 4.32        trt1
## 20 4.69        trt1
## 21 6.31        trt2
## 22 5.12        trt2
## 23 5.54        trt2
## 24 5.50        trt2
## 25 5.37        trt2
## 26 5.29        trt2
## 27 4.92        trt2
## 28 6.15        trt2
## 29 5.80        trt2
## 30 5.26        trt2

5 Análisis de la Varianza

  • Ajustamos un modelo ANAVA para evaluar diferencias entre tratamientos en un DCA.

  • Vamos a usar la función aov() que es parte del paquete básico de R.

  • El ANAVA solo es válido si se cumplen los supuesto.

  • Función con argumentos:

    aov(VARIABLE_RESPUESTA ~ FACTOR, datos = BASE_DE_DATOS)

# Creamos el objeto modelo_dca
modelo_dca <- aov(peso ~ tratamiento, data = bd_crecimiento_plantas)
summary(modelo_dca)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## tratamiento  2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6 Supuestos del modelo

Trabajamos con los residuos y predichos del modelo para realizar las pruebas de los errores.

residuos_dca <- modelo_dca$residuals # residuos = observado - predicho
predichos_dca <- modelo_dca$fitted.values # predichos
tabla <- data.frame(
  Observado = bd_crecimiento_plantas$peso,
  Predichos  = predichos_dca,
  Residual  = residuos_dca,
  Tratamiento = bd_crecimiento_plantas$tratamiento)

print(tabla)
##    Observado Predichos Residual Tratamiento
## 1       4.17     5.032   -0.862        ctrl
## 2       5.58     5.032    0.548        ctrl
## 3       5.18     5.032    0.148        ctrl
## 4       6.11     5.032    1.078        ctrl
## 5       4.50     5.032   -0.532        ctrl
## 6       4.61     5.032   -0.422        ctrl
## 7       5.17     5.032    0.138        ctrl
## 8       4.53     5.032   -0.502        ctrl
## 9       5.33     5.032    0.298        ctrl
## 10      5.14     5.032    0.108        ctrl
## 11      4.81     4.661    0.149        trt1
## 12      4.17     4.661   -0.491        trt1
## 13      4.41     4.661   -0.251        trt1
## 14      3.59     4.661   -1.071        trt1
## 15      5.87     4.661    1.209        trt1
## 16      3.83     4.661   -0.831        trt1
## 17      6.03     4.661    1.369        trt1
## 18      4.89     4.661    0.229        trt1
## 19      4.32     4.661   -0.341        trt1
## 20      4.69     4.661    0.029        trt1
## 21      6.31     5.526    0.784        trt2
## 22      5.12     5.526   -0.406        trt2
## 23      5.54     5.526    0.014        trt2
## 24      5.50     5.526   -0.026        trt2
## 25      5.37     5.526   -0.156        trt2
## 26      5.29     5.526   -0.236        trt2
## 27      4.92     5.526   -0.606        trt2
## 28      6.15     5.526    0.624        trt2
## 29      5.80     5.526    0.274        trt2
## 30      5.26     5.526   -0.266        trt2

6.1 Análisis gráfico de los residuos

hist(residuos_dca, main = "Histograma de residuos")

boxplot(residuos_dca, main = "Boxplot de residuos")

6.2 Normalidad

6.2.1 Gráfico Q-Q plot

plot(modelo_dca, which = 2)

6.2.2 Prueba de Shapiro-Wilks

H0: SI se cumple el supuesto de Normalidad.

H1: NO se cumple el supuesto de Normalidad.

shapiro.test(modelo_dca$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_dca$residuals
## W = 0.96607, p-value = 0.4379

6.3 Homocedasticidad

6.3.1 Gráfico de residuos vs predichos

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

6.3.2 Test de Levene

H0: SI cumple el supuesto de Homocedasticidad.

H1: NO se cumple el supuesto de Homocedasticidad

# Para esta prueba usamos una función de la librería "car"
leveneTest(peso ~ tratamiento, data = bd_crecimiento_plantas)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.1192 0.3412
##       27

7 Pruebas post hoc

7.1 Test de tukey

# Usaremos una función del paquete agricolae
TUKEY <- HSD.test(modelo_dca, "tratamiento")
TUKEY
## $statistics
##     MSerror Df  Mean       CV       MSD
##   0.3885959 27 5.073 12.28809 0.6912161
## 
## $parameters
##    test      name.t ntr StudentizedRange alpha
##   Tukey tratamiento   3         3.506426  0.05
## 
## $means
##       peso       std  r        se  Min  Max    Q25   Q50    Q75
## ctrl 5.032 0.5830914 10 0.1971284 4.17 6.11 4.5500 5.155 5.2925
## trt1 4.661 0.7936757 10 0.1971284 3.59 6.03 4.2075 4.550 4.8700
## trt2 5.526 0.4425733 10 0.1971284 4.92 6.31 5.2675 5.435 5.7350
## 
## $comparison
## NULL
## 
## $groups
##       peso groups
## trt2 5.526      a
## ctrl 5.032     ab
## trt1 4.661      b
## 
## attr(,"class")
## [1] "group"
plot(TUKEY)

7.2 Test de Fisher (LSD)

# Usamos librería agricolae
FISHER <- LSD.test(modelo_dca, "tratamiento")  
FISHER
## $statistics
##     MSerror Df  Mean       CV  t.value       LSD
##   0.3885959 27 5.073 12.28809 2.051831 0.5720126
## 
## $parameters
##         test p.ajusted      name.t ntr alpha
##   Fisher-LSD      none tratamiento   3  0.05
## 
## $means
##       peso       std  r        se      LCL      UCL  Min  Max    Q25   Q50
## ctrl 5.032 0.5830914 10 0.1971284 4.627526 5.436474 4.17 6.11 4.5500 5.155
## trt1 4.661 0.7936757 10 0.1971284 4.256526 5.065474 3.59 6.03 4.2075 4.550
## trt2 5.526 0.4425733 10 0.1971284 5.121526 5.930474 4.92 6.31 5.2675 5.435
##         Q75
## ctrl 5.2925
## trt1 4.8700
## trt2 5.7350
## 
## $comparison
## NULL
## 
## $groups
##       peso groups
## trt2 5.526      a
## ctrl 5.032     ab
## trt1 4.661      b
## 
## attr(,"class")
## [1] "group"
plot(FISHER)

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 DCA los tratamientos se asignan al azar a las unidades experimentales, sin considerar bloques ni restricciones adicionales. Esto permitira asegurarnos el cumplimiento del supuesto de independencia.

# Creamos el vector de los tratamientos
trat <- c("Control", "trat1", "trat2")

# Vector de repeticiones por trratamiento
rep_dca <- 10

# design.crd: Completely Randomized Design
# La función design.crd() permite generar la asignación aleatoria de tratamientos a parcelas o unidades exp.

dca <- design.crd(trat, rep_dca, serie = 0, seed = 123) 

# Mostar la tabla ("field book") con la asignación aleatoria de tratamientos.
# plots = parcelas experimentales (en inglés: plots = “parcelas” o “lotes” de campo).

dca$book 
##    plots  r    trat
## 1      1  1   trat2
## 2      2  1 Control
## 3      3  2 Control
## 4      4  3 Control
## 5      5  4 Control
## 6      6  5 Control
## 7      7  2   trat2
## 8      8  1   trat1
## 9      9  3   trat2
## 10    10  6 Control
## 11    11  4   trat2
## 12    12  2   trat1
## 13    13  7 Control
## 14    14  5   trat2
## 15    15  6   trat2
## 16    16  8 Control
## 17    17  7   trat2
## 18    18  3   trat1
## 19    19  8   trat2
## 20    20  9   trat2
## 21    21  4   trat1
## 22    22  5   trat1
## 23    23  9 Control
## 24    24 10 Control
## 25    25  6   trat1
## 26    26  7   trat1
## 27    27  8   trat1
## 28    28  9   trat1
## 29    29 10   trat1
## 30    30 10   trat2

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…

dca_random <- design.crd(trat, rep_dca, serie = 0, randomization = TRUE) 
dca_random$book
##    plots  r    trat
## 1      1  1   trat1
## 2      2  2   trat1
## 3      3  1 Control
## 4      4  1   trat2
## 5      5  3   trat1
## 6      6  2 Control
## 7      7  4   trat1
## 8      8  2   trat2
## 9      9  5   trat1
## 10    10  3   trat2
## 11    11  3 Control
## 12    12  6   trat1
## 13    13  4 Control
## 14    14  5 Control
## 15    15  6 Control
## 16    16  7 Control
## 17    17  4   trat2
## 18    18  5   trat2
## 19    19  6   trat2
## 20    20  7   trat1
## 21    21  8 Control
## 22    22  8   trat1
## 23    23  9 Control
## 24    24  7   trat2
## 25    25  8   trat2
## 26    26 10 Control
## 27    27  9   trat2
## 28    28  9   trat1
## 29    29 10   trat1
## 30    30 10   trat2