Diseño Experimental

Cátedra de Bioestadística y Diseño Experimental. Equipo docente responsable. Año 2024.

Ing. Agr. Jorge Quiquinto, prof. Asociado, dedicación exclusiva.
Ing. Agr. Marta Leaño¹, prof. Adjunta, dedicación exclusiva (*).
Ing. Agr. Juan Manuel Solís², Jefe de Trabajos Prácticos, dedicación exclusiva (*).
Ing. Agr. Ivone Humacata, Jefe de Trabajos Prácticos, dedicación exclusiva.
Srta. Victoria López³, Ayudante de 2da, dedicación simple (*).
Sr. Daniel Vilca, Ayudante de 2da, dedicación simple.

(*) Equipo responsable del curso Diseño Experimental 2024.

Recomendación. No puede continuar con la clase si no tiene instaladas las siguientes librerías:

library(agricolae)

install.packages("agricolae") # Instalar

library(emmeans)

install.packages("emmeans") # Instalar

library(multcomp)

install.packages("multcomp") # Instalar

Recomiendo su instalación previa a la clase.

En caso de que su ordenador no le permita descagar las librerías mencionadas anteriormente puede trabajar desde RStudio Cloud (requiere internet durante la clase, tengalo en cuenta).

Para ello debe crearse un usuario (puede ser con su cuenta de google). Y luego debe crear un nuevo proyecto (New RStudio Project).

En la consola (console) ejecute las tres sentencias de install.packages("") (una a la vez) y listo. En clase se mostrará lo que se hace a continuación. Muchas gracias por leer las recomendaciones.

EJERCICIO 1

Se realizó un estudio sobre las diferencias en el canto posado de tres especies cercanamente emparentadas y simpátricas¹ del género Sturnella en el sur de la provincia de Buenos Aires, Argentina. El objetivo del estudio fue evaluar la importancia de estas diferencias en el mantenimiento del aislamiento reproductivo entre las especies. Una de las variables registradas fue la duración de la parte introductoria de la canción entera, medida en segundos. ¿Existen diferencias significativas entre las especies del género Sturnella? ¿Qué tratamientos son diferentes (estadísticamente)? Se obtuvieron los siguientes datos:

DATOS EJECICIO 1

S. loyca	S. superciliaris	S. defilippii
0.97	1.12	0.69
1.23	0.97	1.36
1.45	1.26	0.56
0.89	1.20	0.96
1.43	1.41	1.45
1.36	1.22	0.76
1.38	0.82	0.80
1.80	1.04	0.75
0.86	1.46	0.85
1.86	1.44	0.98

IMPORTAR DATOS A R

sturnella = data.frame(
  Tratamiento = factor(
    rep(
      c("loyca","superciliaris","defilippii"),
      times = 10
    )
  ),
  Duracion = c(
    0.97,   1.12,   0.69,
    1.23,   0.97,   1.36,
    1.45,   1.26,   0.56,
    0.89,   1.20,   0.96,
    1.43,   1.41,   1.45,
    1.36,   1.22,   0.76,
    1.38,   0.82,   0.80,
    1.80,   1.04,   0.75,
    0.86,   1.46,   0.85,
    1.86,   1.44,   0.98
  )
)
sturnella

sturnella

Tratamiento	Duracion
loyca	0.97
superciliaris	1.12
defilippii	0.69
loyca	1.23
superciliaris	0.97
defilippii	1.36
loyca	1.45
superciliaris	1.26
defilippii	0.56
loyca	0.89
superciliaris	1.20
defilippii	0.96
loyca	1.43
superciliaris	1.41
defilippii	1.45
loyca	1.36
superciliaris	1.22
defilippii	0.76
loyca	1.38
superciliaris	0.82
defilippii	0.80
loyca	1.80
superciliaris	1.04
defilippii	0.75
loyca	0.86
superciliaris	1.46
defilippii	0.85
loyca	1.86
superciliaris	1.44
defilippii	0.98

sturnella

str(sturnella)

'data.frame':   30 obs. of  2 variables:
 $ Tratamiento: Factor w/ 3 levels "defilippii","loyca",..: 2 3 1 2 3 1 2 3 1 2 ...
 $ Duracion   : num  0.97 1.12 0.69 1.23 0.97 1.36 1.45 1.26 0.56 0.89 ...

sturnella$Tratamiento

 [1] loyca         superciliaris defilippii    loyca         superciliaris
 [6] defilippii    loyca         superciliaris defilippii    loyca        
[11] superciliaris defilippii    loyca         superciliaris defilippii   
[16] loyca         superciliaris defilippii    loyca         superciliaris
[21] defilippii    loyca         superciliaris defilippii    loyca        
[26] superciliaris defilippii    loyca         superciliaris defilippii   
Levels: defilippii loyca superciliaris

Relación Objeto$Atributos

Objeto que almacena atributos$Atributos
sturnella$Tratamiento
sturnella$Duracion

GRÁFICO DE CAJAS

Tratamiento	Media	Desvio	Mediana	n
defilippii	0.92	0.29	0.82	10
loyca	1.32	0.35	1.37	10
superciliaris	1.19	0.21	1.21	10

CONSTRUCCIÓN DEL MODELO Y ANAVA

modelo = lm(Duracion ~ Tratamiento, data = sturnella)
anova(modelo)

Analysis of Variance Table

Response: Duracion
            Df  Sum Sq Mean Sq F value  Pr(>F)  
Tratamiento  2 0.86525 0.43262  5.2586 0.01178 *
Residuals   27 2.22129 0.08227                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Existen diferencias significativas en la duración de la parte introductoria de la canción entera entre las tres especies del género Sturnella?

Hipótesis

\[ H_0: \mu_1=\mu_2=\mu_3 \]
\[ H_1: \text{Al menos una de las medias es diferente} \]

¿Pero cuál?

SUPUESTOS

#Test de Shapiro-Wilk para evaluar Normalidad
residuos = resid(modelo)
shapiro.test(residuos)


    Shapiro-Wilk normality test

data:  residuos
W = 0.96114, p-value = 0.3311

#Test de Bartlett para evaluar Homogeneidad de Varianzas
bartlett.test(residuos, sturnella$Tratamiento)


    Bartlett test of homogeneity of variances

data:  residuos and sturnella$Tratamiento
Bartlett's K-squared = 2.0236, df = 2, p-value = 0.3636

¿El modelo es válido?

CONTINUACION DEL ANAVA

Comparaciones a posteriori

Se realizan solamente luego de encontrar efectos significativos en el ANAVA

\[\begin{gather*} \text{p-valor }< \alpha \end{gather*}\]

Test de Scheffé

Generalmente se lo aplica en comparaciones a posteriori
Prueba todo y cualquier contraste
Condición de ortogonalidad de los contrastes no necesaria
Poco poderoso
Útil cuando se compara más de dos medias (posteriori)

Código

library(agricolae)
scheffe = scheffe.test(
  y = modelo,
  trt = "Tratamiento",
  alpha = 0.05,
  group = TRUE)

scheffe.test()

Requiere:

Importar datos (sturnella)
Modelo (modelo = lm(Duracion ~ Tratamiento, data = sturnella))

Scheffé

scheffe$groups

              Duracion groups
loyca            1.323      a
superciliaris    1.194     ab
defilippii       0.916      b

Letras distintas indican diferencias significativas de acuerdo a los resultados de Scheffé

Gráfico Scheffé

library(agricolae)
bar.group(x = scheffe$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,2),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="Sturnella",
      ylab="Duración")

Test de Tukey

Comparación a posteriori
Más usado
Compara todo y cualquier contraste entre dos medias

\[\begin{gather*} HSD = q * \sqrt{\dfrac{CMEE}{repeticiones}} \text{ donde } q_{(\alpha;n_{trat};gl_{EE})} \end{gather*}\]

Código

library(agricolae)
Tukey= HSD.test(
  y = modelo,
  trt = "Tratamiento",
  alpha = 0.05,
  group = TRUE,
  unbalanced = FALSE)

HSD.test()

Requiere:

Importar datos (sturnella)
Modelo (modelo = lm(Duracion ~ Tratamiento, data = sturnella))

Tukey

Tukey$groups

              Duracion groups
loyca            1.323      a
superciliaris    1.194     ab
defilippii       0.916      b

Letras distintas indican diferencias significativas de acuerdo a los resultados de Tukey

Gráfico Tukey

library(agricolae)
bar.group(x = Tukey$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,2),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="Sturnella",
      ylab="Duración")

HSD (Diferencia Honestamente Significativa)

#Diferencia Mínina Significativa
MSD = Tukey$statistics
MSD

	MSerror	Df	Mean	CV	MSD
	0.08227	27	1.144333	25.06503	0.3180427

\[\begin{gather*} HSD = q * \sqrt{\dfrac{CMEE}{repeticiones}} \text{ donde } q_{(\alpha;n_{trat};gl_{EE})} \end{gather*}\]

Ejemplo

$d_1 = \overline{x}_{superciliaris}-\overline{x}_{defilippii}$

$d_1=0.407$

q

#Valor q
q = Tukey$parameters
q

	test	name.t	ntr	StudentizedRange	alpha
	Tukey	Tratamiento	3	3.506426	0.05

qtukey(p = 0.05,
       nmeans = q$ntr,
       df = MSD$Df,
       lower.tail = FALSE)

[1] 3.506426

Test de Duncan

Comparación a posteriori
Compara contrastes entre dos medias
Considera el número de medias que abarca el contraste
Medias ordenadas de mayor a menor
Más sensible que el de Tukey
Sí aumenta el n° de tratamientos, disminuye el grado de confiabilidad

Código

library(agricolae)
duncan = duncan.test(
  y = modelo,
  trt = "Tratamiento",
  alpha = 0.05,
  group = TRUE)

duncan.test()

Requiere:

Importar datos (sturnella)
Modelo (modelo = lm(Duracion ~ Tratamiento, data = sturnella))

Duncan

duncan$groups

              Duracion groups
loyca            1.323      a
superciliaris    1.194      a
defilippii       0.916      b

Letras distintas indican diferencias significativas de acuerdo a los resultados de Duncan

Gráfico Duncan

library(agricolae)
bar.group(x = duncan$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,2),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="Sturnella",
      ylab="Duración")

EJERCICIO 2

Se sabe que el dióxido de carbono tiene un efecto crítico en el crecimiento microbiológico. Cantidades pequeñas de $CO_2$ estimulan el crecimiento de muchos microorganismos, mientras que altas concentraciones inhiben el crecimiento de la mayor parte de ellos. Este último efecto se utiliza comercialmente cuando se almacenan productos alimenticios perecederos. Se realizó un estudio para investigar el efecto de $CO_2$ sobre la tasa de crecimiento de Pseudomonas fragi, un corruptor de alimentos. Se administró $CO_2$ a cinco presiones atmosféricas diferentes (0.0, 0.083, 0.29, 0.50, 0.86). La respuesta anotada fue el cambio porcentual en la masa celular después de un tiempo de crecimiento de una hora. ¿Qué tratamientos son diferentes (estadísticamente)?

DATOS EJERCICIO 2

Se obtuvieron los siguientes datos:

0.0	0.083	0.29	0.5	0.86
62.6	50.9	45.5	29.5	24.9
59.6	44.3	41.1	22.8	17.2
64.5	47.5	29.8	19.2	7.8
59.3	49.5	38.3	20.6	10.5
58.6	48.5	40.2	29.2	17.8
64.6	50.4	38.5	24.1	22.1
50.9	35.2	30.2	22.6	22.6
56.2	49.9	27.0	32.7	16.8
52.3	42.6	40.0	24.4	15.9
62.8	41.6	33.9	29.6	8.8

IMPORTAR DATOS A R

datos = data.frame(
  CO2 = factor(
    rep(
      c("0.0","0.083","0.29","0.50","0.86"),
      times = 10
    )
  ),
  respuesta = c(
    62.6,   50.9,   45.5,   29.5,   24.9,
    59.6,   44.3,   41.1,   22.8,   17.2,
    64.5,   47.5,   29.8,   19.2,   7.8,
    59.3,   49.5,   38.3,   20.6,   10.5,
    58.6,   48.5,   40.2,   29.2,   17.8,
    64.6,   50.4,   38.5,   24.1,   22.1,
    50.9,   35.2,   30.2,   22.6,   22.6,
    56.2,   49.9,   27.0,   32.7,   16.8,
    52.3,   42.6,   40.0,   24.4,   15.9,
    62.8,   41.6,   33.9,   29.6,   8.8
  )
)
datos

datos

CO2	respuesta
0.0	62.6
0.083	50.9
0.29	45.5
0.50	29.5
0.86	24.9
0.0	59.6
0.083	44.3
0.29	41.1
0.50	22.8
0.86	17.2
0.0	64.5
0.083	47.5
0.29	29.8
0.50	19.2
0.86	7.8
0.0	59.3
0.083	49.5
0.29	38.3
0.50	20.6
0.86	10.5
0.0	58.6
0.083	48.5
0.29	40.2
0.50	29.2
0.86	17.8
0.0	64.6
0.083	50.4
0.29	38.5
0.50	24.1
0.86	22.1
0.0	50.9
0.083	35.2
0.29	30.2
0.50	22.6
0.86	22.6
0.0	56.2
0.083	49.9
0.29	27.0
0.50	32.7
0.86	16.8
0.0	52.3
0.083	42.6
0.29	40.0
0.50	24.4
0.86	15.9
0.0	62.8
0.083	41.6
0.29	33.9
0.50	29.6
0.86	8.8

GRÁFICO DE CAJAS

CO2	Media	Desvio	Mediana	n
0.0	59.14	4.80	59.45	10
0.083	46.04	5.05	48.00	10
0.29	36.45	5.93	38.40	10
0.50	25.47	4.48	24.25	10
0.86	16.44	5.89	17.00	10

CONSTRUCCIÓN DEL MODELO Y ANAVA

modelo.ej2 = lm(respuesta ~ CO2, data = datos)
anova(modelo.ej2)

Analysis of Variance Table

Response: respuesta
          Df Sum Sq Mean Sq F value    Pr(>F)    
CO2        4  11274 2818.58  101.63 < 2.2e-16 ***
Residuals 45   1248   27.73                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Existen diferencias significativas en el cambio porcentual en la masa celular entre las cinco presiones atmosféricas de $CO_2$?

Hipótesis

\[ H_0: \mu_1=\mu_2=\mu_3=\mu_4=\mu_5 \]
\[ H_1: \text{Al menos una de las medias es diferente} \]

¿Pero cuál?

SUPUESTOS

#Test de Shapiro-Wilk para evaluar Normalidad
residuos.ej2 = resid(modelo.ej2)
shapiro.test(residuos.ej2)


    Shapiro-Wilk normality test

data:  residuos.ej2
W = 0.9627, p-value = 0.1153

#Test de Bartlett para evaluar Homogeneidad de Varianzas
bartlett.test(residuos.ej2,datos$CO2)


    Bartlett test of homogeneity of variances

data:  residuos.ej2 and datos$CO2
Bartlett's K-squared = 1.0701, df = 4, p-value = 0.899

¿El modelo es válido?

COMPARACIONES A priori

	0.0	0.083	0.29	0.50	0.86
0.0 vs. 0.083,0.29,0.50 y 0.86	4	-1	-1	-1	-1
0.083 vs. 0.29,0.50 y 0.86	0	3	-1	-1	-1
0.29 vs. 0.50 y 0.86	0	0	2	-1	-1
0.50 vs. 0.86	0	0	0	1	-1

CONSTRASTES Y CONSTRASTES ORTOGONALES

c1 = c(4,-1,-1,-1,-1)
c2 = c(0,3,-1,-1,-1)
c3 = c(0,0,2,-1,-1)
c4 = c(0,0,0,1,-1)
sum(c1)    #Contraste

[1] 0

sum(c1*c2) #Constrastes ortogonales

[1] 0

#¿c2:c4 son contrastes?
#¿se cumple la ortogonalidad entre el resto de los contrastes?

Código

library(emmeans)
medias.modelo=emmeans(object = modelo.ej2,
                      specs = ~CO2)
medias.modelo

CO2	emmean	SE	df	lower.CL	upper.CL
0.0	59.14	1.665358	45	55.7858	62.4942
0.083	46.04	1.665358	45	42.6858	49.3942
0.29	36.45	1.665358	45	33.0958	39.8042
0.50	25.47	1.665358	45	22.1158	28.8242
0.86	16.44	1.665358	45	13.0858	19.7942

emmeans()

Requiere:

Importar datos (datos)
Modelo (modelo.ej2 = lm(respuesta ~ CO2, data = datos))

Código

contrastes = list(c1, c2, c3,c4)
library(emmeans)
contrast(object = medias.modelo,
         contrastes)

 contrast             estimate   SE df t.ratio p.value
 c(4, -1, -1, -1, -1)   112.16 7.45 45  15.060  <.0001
 c(0, 3, -1, -1, -1)     59.76 5.77 45  10.359  <.0001
 c(0, 0, 2, -1, -1)      30.99 4.08 45   7.597  <.0001
 c(0, 0, 0, 1, -1)        9.03 2.36 45   3.834  0.0004

El contraste difiere significativamente sí y solo sí el p-valor es < 0.05

Test de Dunnett

Comparación a posteriori
Comparación de a pares
Control o testigo

Código

class(datos$CO2)

[1] "factor"

levels(datos$CO2)

[1] "0.0"   "0.083" "0.29"  "0.50"  "0.86"

library(multcomp)
Dunnett = glht(model = modelo.ej2,
              linfct = mcp(CO2="Dunnett"))

summary(Dunnett)

Atención: Observar salida summary(Dunnett) en R.

El contraste difiere significativamente sí y solo sí el p-valor es < 0.05

glht()

Requiere:

Importar datos (datos)
Modelo (modelo.ej2 = lm(respuesta ~ CO2, data = datos))

Gráfico Dunnett

plot(Dunnett,
     cex.axis = 0.72,
     col = "#8B0000")

Test LSD

Código

Comparación a posteriori

library(agricolae)
LSD = LSD.test(
  y = modelo.ej2,
  trt = "CO2",
  alpha = 0.05,
  p.adj= "none",
  group = TRUE)

LSD.test()

Requiere:

Importar datos (datos)
Modelo (modelo.ej2 = lm(respuesta ~ CO2, data = datos))

Código

LSD$groups

      respuesta groups
0.0       59.14      a
0.083     46.04      b
0.29      36.45      c
0.50      25.47      d
0.86      16.44      e

Letras distintas indican diferencias significativas de acuerdo a los resultados de LSD

Gráfico LSD

library(agricolae)
bar.group(x = LSD$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,66),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="CO2",
      ylab="Variable respuesta")

Especificar control

datos$CO2 = factor(datos$CO2, 
                   levels = c("0.86","0.50","0.29","0.083","0.0"),
                   ordered = TRUE)
modelo.ej2.2 = lm(respuesta ~ CO2, data = datos)
library(multcomp)
Dunnett.0.86 = glht(model = modelo.ej2.2,
                    linfct = mcp(CO2 = "Dunnett"))
summary(Dunnett.0.86)

Atención: Observar salida summary(Dunnett.0.86) en R.

El contraste difiere significativamente sí y solo sí el p-valor es < 0.05

Gráfico Dunnett.0.86

plot(Dunnett.0.86,
     cex.axis = 0.70,
     col = "#8B0000")

EJERCICIO DBCA (CLASE 2)

DATOS = data.frame(
  trat = factor(rep(c("A","B","C"),
                    each = 5)),
  bl = factor(rep(c("I","II","III","IV","V"),
                  times = 3)),
  as = c(9.04,8.91,9.68,10.36,10.2,
         10.32,10.04,10.8,11.11,11.42,
         10.04,9.79,10.44,11.1,10.76)
)
DATOS

DATOS DBCA

trat	bl	as
A	I	9.04
A	II	8.91
A	III	9.68
A	IV	10.36
A	V	10.20
B	I	10.32
B	II	10.04
B	III	10.80
B	IV	11.11
B	V	11.42
C	I	10.04
C	II	9.79
C	III	10.44
C	IV	11.10
C	V	10.76

Test de Tukey

modelo = lm(as ~ trat + bl, DATOS)
library(agricolae)
Tukey= HSD.test(
  y = modelo,
  trt = "trat",
  alpha = 0.05,
  group = TRUE,
  unbalanced = FALSE)

Tukey

Tukey$groups

      as groups
B 10.738      a
C 10.426      b
A  9.638      c

Letras distintas indican diferencias significativas de acuerdo a los resultados de Tukey

Gráfico Tukey

library(agricolae)
bar.group(x = Tukey$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,15),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="Especies",
      ylab="Concentración de As en ppm")

Test de Duncan

library(agricolae)
duncan = duncan.test(
  y = modelo,
  trt = "trat",
  alpha = 0.05,
  group = TRUE)

duncan.test()

Requiere:

Importar datos (DATOS)
Modelo (modelo = lm(as ~ trat + bl, data = DATOS))

Duncan

duncan$groups

      as groups
B 10.738      a
C 10.426      b
A  9.638      c

Letras distintas indican diferencias significativas de acuerdo a los resultados de Duncan

Gráfico Duncan

library(agricolae)
bar.group(x = duncan$groups,
          horiz = FALSE,
          decreasing = TRUE,
          ylim=c(0,15),
          border="#710F99",
          col = "#E9B2FF",
          las = 1)
title(cex.main=1,
      main="Comparación entre \nlas medias de los tratamientos",
      xlab="Especies",
      ylab="Concetración de As en ppm")

Bibliografía

Garibaldi, L. A., Oddi, F. J., Aristimuño, F. J., Behnisch, A. N. Modelos estadísticos en lenguaje R. 2019.
Quinteros, H. Diseño Experimental. 2013.
Glosario (link de interés)

Muchas gracias por su atención