Estructura del Diseño Experimental

Con el Diseño en Cuadro Grecolatino (DCGL) se controlan tres factores de bloque, además del factor de tratamiento. Se llama cuadro grecolatino porque los cuatro factores involucrados se prueban en la misma cantidad de niveles, de aquí que se pueda escribir como un cuadro, como se muestra en la Tabla 1. Se utilizan letras para denotar los tratamientos, y letras griegas para nombrar los niveles del tercer factor de bloque, mientras que el primer y segundo factor de bloque se ubican en los renglones y en las columnas (Pulido & Vara Salazar, 2012):
Tabla 1. Estructura del Diseño en Cuadro Grecolatino
1 2 3 4
A\(\alpha\) B\(\beta\) C\(\gamma\) D\(\delta\)
B\(\delta\) A\(\gamma\) C\(\beta\) D\(\alpha\)
C\(\beta\) D\(\alpha\) A\(\delta\) B\(\gamma\)
D\(\gamma\) C\(\delta\) B\(\alpha\) A\(\beta\)
El modelo estadístico que describe las mediciones en un cuadro grecolatino está dado por:
\[y_{ijk}=\mu+\tau_i+\beta_j+\delta_k+\varphi_l+\varepsilon_{ijkl}\]
Donde \(y_{ijkl}\) es la observación o respuesta que se encuentra en el tratamiento i (i-ésima letra latina), en el renglón j, en la columna k y en la l-ésima letra griega, que son los niveles del tercer factor de bloque; el término \(\varepsilon_{ijkl}\) representa el error aleatoio atribuible a la medición \(y_{ijkl}\).

##Caso Práctico

Para ejemplificar la metodología del DCGL tomaremos como base el ejercicio 25 de la página 111 del libro de Análsis y Diseño de Experimentos (Pulido & Vara Salazar, 2012), en el que se considera que un investigador está interesado en el efecto del porcentaje de lisina y del porcentaje de de proteína en la producción de vacas lecheras. Se consideran siete niveles de cada factor:
  • Porcentaje de lisina: 0.0(A),0.1(B),0.2(C),0.3(D),0.4(E),0.5(F),0.6(G).
  • Porcentaje de proteína: 2(\(\alpha\)),4(\(\beta\)),6(\(\chi\)),8(\(\delta\)),10(\(\varepsilon\)),12(\(\varphi\)),14(\(\gamma\)).
Para el estudio, se seleccionan siete vacas al azar, a las cuales se les da un seguimiento de siete periodos de tres meses. Los datos en galones de leche fueron los siguientes:
Tabla 2. Resultados del experimento.
Vaca/Periodo 1 2 3 4 5 6 7
1 304(A\(\alpha\)) 436(B\(\varepsilon\)) 350(C\(\beta\)) 504(D\(\varphi\)) 417(E\(\chi\)) 519(F\(\gamma\)) 432(G\(\delta\))
2 381(B\(\beta\)) 505(C\(\varphi\)) 425(D\(\chi\)) 564(E\(\gamma\)) 494(F\(\delta\)) 350(G\(\alpha\)) 413(GA\(\varepsilon\))
3 432(C\(\chi\)) 566(D\(\gamma\)) 479(E\(\delta\)) 357(F\(\alpha\)) 461(G\(\varepsilon\)) 340(A\(\beta\)) 502(B\(\varphi\))
4 442(D\(\delta\)) 372(E\(\alpha\)) 536(F\(\varepsilon\)) 366(G\(\beta\)) 495(A\(\varphi\)) 425(B\(\chi\)) 507(C\(\gamma\))
5 496(E\(\varepsilon\)) 449(F\(\beta\)) 493(G\(\varphi\)) 345(A\(\chi\)) 509(B\(\gamma\)) 481(C\(\delta\)) 380(D\(\alpha\))
6 534(F\(\varphi\)) 421(G\(\chi\)) 352(A\(\gamma\)) 427(B\(\delta\)) 346(C\(\alpha\)) 478(D\(\varepsilon\)) 397(E\(\beta\))
7 543(G\(\gamma\)) 386(A\(\delta\)) 435(B\(\alpha\)) 485(C\(\varepsilon\)) 406(D\(\beta\)) 554(E\(\varphi\)) 410(F\(\chi\))
  1. Analice este experimento. ¿Qué factores tienen efecto en la producción de leche?
  2. Interprete los resultados usando gráfico de medias.
  3. ¿Cómo puede explicar la falta de efectos en vacas y periodo?
  4. ¿Qué porcentajes de lisina y proteína dan los mejores resultados?
  5. Verifique los supuestos del modelo.

Analisis del experimento

Para analizar este experimento, es importante reacomodar los datos en un dataset de manera adecuada para que el código pueda ejecutarse, para esto, utilizaremos la libreria readxl, misma que permite la lectura de datos desde un archivo de Excel. 
library(readxl)
datos=read_excel("dataset.xlsx")
attach(datos)
print(datos)
## # A tibble: 49 x 5
##    Galones Lisina  Vaca Periodo Proteina
##      <dbl> <chr>  <dbl>   <dbl> <chr>   
##  1     304 a          1       1 a       
##  2     381 b          2       1 b       
##  3     432 c          3       1 c       
##  4     442 d          4       1 d       
##  5     496 e          5       1 e       
##  6     534 f          6       1 f       
##  7     543 g          7       1 g       
##  8     436 b          1       2 e       
##  9     505 c          2       2 f       
## 10     566 d          3       2 g       
## # … with 39 more rows
En este caso, como es de observarse, se cambiaron los datos del dataset en el caso de la columna que contiene las letras griegas, para que los datos pudieran ser procesados adecuadamente.
Tomando como base el modelo estadístico propuesto al inicio de este manuscrito, se programa la siguiente etapa del análisis, escribiendo como factores las columnas correspondientes al Porcentaje de Lisina, las Vacas, los Periodos y el Porcentaje de Proteína. 
f_lisina=factor(Lisina)
f_vacas=factor(Vaca)
f_periodo=factor(Periodo)
f_proteina=factor(Proteina)
modelo=lm(Galones~(f_lisina+f_vacas+f_periodo+f_proteina))
summary(modelo)
## 
## Call:
## lm(formula = Galones ~ (f_lisina + f_vacas + f_periodo + f_proteina))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.082 -13.796  -0.082  11.204  56.776 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 281.7959    22.4284  12.564 4.80e-12 ***
## f_lisinab    68.5714    16.7839   4.086 0.000424 ***
## f_lisinac    67.2857    16.7839   4.009 0.000515 ***
## f_lisinad    80.8571    16.7839   4.818 6.60e-05 ***
## f_lisinae    92.0000    16.7839   5.481 1.23e-05 ***
## f_lisinaf    94.8571    16.7839   5.652 8.07e-06 ***
## f_lisinag    61.5714    16.7839   3.668 0.001212 ** 
## f_vacas2     24.2857    16.7839   1.447 0.160843    
## f_vacas3     25.0000    16.7839   1.490 0.149375    
## f_vacas4     25.8571    16.7839   1.541 0.136499    
## f_vacas5     27.2857    16.7839   1.626 0.117073    
## f_vacas6     -1.0000    16.7839  -0.060 0.952983    
## f_vacas7     36.7143    16.7839   2.187 0.038684 *  
## f_periodo2    0.4286    16.7839   0.026 0.979840    
## f_periodo3   -8.8571    16.7839  -0.528 0.602542    
## f_periodo4  -12.0000    16.7839  -0.715 0.481525    
## f_periodo5   -0.5714    16.7839  -0.034 0.973122    
## f_periodo6    2.1429    16.7839   0.128 0.899471    
## f_periodo7  -13.0000    16.7839  -0.775 0.446169    
## f_proteinab  20.7143    16.7839   1.234 0.229087    
## f_proteinac  47.2857    16.7839   2.817 0.009537 ** 
## f_proteinad  85.2857    16.7839   5.081 3.38e-05 ***
## f_proteinae 108.7143    16.7839   6.477 1.07e-06 ***
## f_proteinaf 149.0000    16.7839   8.878 4.77e-09 ***
## f_proteinag 145.1429    16.7839   8.648 7.74e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.4 on 24 degrees of freedom
## Multiple R-squared:  0.8938, Adjusted R-squared:  0.7876 
## F-statistic: 8.417 on 24 and 24 DF,  p-value: 8.969e-07
En el caso de este modelo estadísitico, puede concluirse que la variabilidad de los galones de leche producidos por las vacas sometidas a las pruebas está explicada en un 89.38%, sin embargo dado que el modelo matemático tiene agregados factores de bloque, es conveniente observar el Coeficiente de Determinación Ajustado, el cual arroja como resultado que la variabilidad de los galones de leche producida por las vacas sometidas al experimento está explicada en un 78.76% por los factores involucrados en el análisis, situación que según (Pulido & Vara Salazar, 2012), lleva a la conclusión de que el modelo es apto para predicción dado que el Coeficiente de DEterminación Ajustado es mayor a 70%.
Una vez establecida la medida de ajuste del modelo, procederemos a revisar la gráfica de cajas del modelo, para poder verificar, en terminos de la media de cada factor, cual es su influencia en la variable de respuesta. 
boxplot(Galones~f_lisina, col="red", xlab= "Porcentaje de Lisina", ylab = "Galones de leche")

En este gráfico se puede observar como el nivel A de Lisina genera una menor producción promedio de leche, mientras que, en lo general, para los demas niveles se mantiene la media y la variabilidad en niveles razonablemente constantes. 
boxplot(Galones~f_vacas, col="skyblue", xlab = "Identificador de Vacas", ylab="Galones de leche")

Para el caso del factor vaca, es posible concluir mediante la gráfica de cajas que, en lo general la producción debido al factor de bloque analizado se mantiene en niveles constantes en promedio. 
boxplot(Galones~f_periodo, xlab="Identificador de Periodo", ylab = "Galones de leche")

Para el caso del factor de bloque Periodo, es posible concluir mediante la gráfica de cajas que, en lo general la producción debido al factor considerado se mantiene en niveles constantes en promedio. 
boxplot(Galones~f_proteina, col="orange", xlab="Porcentaje de Proteína", ylab="Galones de leche")

Para el caso del factor proteína, es posible concluir mediante la gráfica de cajas que, la variación en los niveles de proteína suministrados a las unidades experimentales generan cambios aparentemente significativos en la producción de leche.

Para establecer conclusiones con un nivel de confianza definido, en este caso, de 95%, procederemos a generar la Tabla ANOVA,para la cual se probarán las siguientes hipótesis estadísticas:

\[H_o:\tau_i=\tau_j=0\] \[H_1:\tau_i\neq\tau_j\neq0\]
Para el caso del factor proteína, es posible concluir mediante la gráfica de cajas que, la variación en los niveles de proteína suministrados a las unidades experimentales generan cambios aparentemente significativos en la producción de leche .
En esta hipótesis estadística se hace referencia, en el caso de la hipótesis nula, a que el cambio en los niveles de lisina no produce efectos significativos en la producción de leche medida en galones por parte de las unidades experimentales.
En el caso de los factores de bloque, sólo se establecerán las concluisiones de manera parcial debido a la falte de aleatoriedad de los mismos.
La Tabla ANOVA se obtiene mediante la siguiente secuencia de comandos: 
anova=aov(modelo)
summary(anova)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## f_lisina     6  42784    7131   7.232 0.000171 ***
## f_vacas      6   8754    1459   1.480 0.227194    
## f_periodo    6   1761     293   0.298 0.931964    
## f_proteina   6 145880   24313  24.660 3.77e-09 ***
## Residuals   24  23663     986                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tal y como se estableció de manera visual mediante las gráficas de caja y bigote, la evidencia muestral obtenida sugiere, con un 95% de confianza que el cambio en los niveles de lisina produce cambios significativos en la producción de leche por parte de las unidades experimentales, en el caso del factor de bloque vacas, se puede concluir que no existen cambios significativos en la producción de leche por efectos inducidos por las diferencias naturales de los individuos involucrados en el experimento, para el caso del factor de bloque denominado periodo, se puede concluir que no existen cambios significativos en los niveles de producción de leche por efectos o ruido estadístico inducido por el paso del tiempo, sin embargo, para el factor de bloque denominado proteína, se concluye que sí existen cambios significativos en la producción de leche por parte de las unidades experimentales debido a las variaciones de la proteína suministrada durante el experimento, situación que hace sugerir que será necesario realizar un nuevo experimento, considerando a la proteína como un factor de tratamiento para verificar cómo, tanto por separado como en interacción, los factores de lisina y proteína afectan la producción de leche.
Ahora prodeceremos a realizar las pruebas de rango múltiple para el factor de tratamiento, que corresponde a los niveles de lisina y para el factor de bloque que corresponde a los niveles de proteína suministrados, para el caso del factor de bloque denominado vacas y periodo, dado que la evidencia muestral arrojó como resultado que no tienen efectos significativos, no se considerarán en las pruebas de rango múltiple.
En el caso de las pruebas de rango múltiple para el factor de tratamiento denominado lisina, aplicaremos la siguiente secuencia de comandos: 
library(agricolae)
rm_lisina=duncan.test(y=Galones,trt=f_lisina,MSerror = 986, DFerror = 24, alpha = 0.05,group = TRUE)
print(rm_lisina)
## $statistics
##   MSerror Df     Mean      CV
##       986 24 442.8776 7.09014
## 
## $parameters
##     test   name.t ntr alpha
##   Duncan f_lisina   7  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.918793      34.64119
## 3 3.065610      36.38366
## 4 3.159874      37.50241
## 5 3.226454      38.29261
## 6 3.276155      38.88247
## 7 3.314602      39.33878
## 
## $means
##    Galones      std r Min Max   Q25 Q50   Q75
## a 376.4286 62.77701 7 304 495 342.5 352 399.5
## b 445.0000 45.36151 7 381 509 426.0 435 469.0
## c 443.7143 69.90878 7 346 507 391.0 481 495.0
## d 457.2857 63.65196 7 380 566 415.5 442 491.0
## e 468.4286 75.68984 7 372 564 407.0 479 525.0
## f 471.2857 68.58988 7 357 536 429.5 494 526.5
## g 438.0000 68.10776 7 350 543 393.5 432 477.0
## 
## $comparison
## NULL
## 
## $groups
##    Galones groups
## f 471.2857      a
## e 468.4286      a
## d 457.2857      a
## b 445.0000      a
## c 443.7143      a
## g 438.0000      a
## a 376.4286      b
## 
## attr(,"class")
## [1] "group"
plot.group(rm_lisina,main="Grupos para el factor Lisina")

Los resultados de la prueba de rango múltiple arrojan como resultado que, en lo general, todos los niveles de lisina generan en promedio la misma producción de leche, salvo en el caso del nivel A, que corresponde al nivel en el que no se suministra lisina a las unidades experimentales.
En el caso de las pruebas de rango múltiple para el factor de bloque denominado proteína, aplicaremos la siguiente secuencia de comandos: 
rm_proteina=duncan.test(y=Galones,trt=f_proteina,MSerror = 986, DFerror = 24, alpha = 0.05,group = TRUE)
print(rm_proteina)
## $statistics
##   MSerror Df     Mean      CV
##       986 24 442.8776 7.09014
## 
## $parameters
##     test     name.t ntr alpha
##   Duncan f_proteina   7  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.918793      34.64119
## 3 3.065610      36.38366
## 4 3.159874      37.50241
## 5 3.226454      38.29261
## 6 3.276155      38.88247
## 7 3.314602      39.33878
## 
## $means
##    Galones      std r Min Max   Q25 Q50   Q75
## a 363.4286 39.84912 7 304 435 348.0 357 376.0
## b 384.1429 37.19959 7 340 449 358.0 381 401.5
## c 410.7143 29.79214 7 345 432 413.5 421 425.0
## d 448.7143 38.16506 7 386 494 429.5 442 480.0
## e 472.1429 40.36264 7 413 536 448.5 478 490.5
## f 512.4286 22.76589 7 493 554 498.5 504 519.5
## g 508.5714 73.23673 7 352 566 508.0 519 553.5
## 
## $comparison
## NULL
## 
## $groups
##    Galones groups
## f 512.4286      a
## g 508.5714      a
## e 472.1429      b
## d 448.7143      b
## c 410.7143      c
## b 384.1429     cd
## a 363.4286      d
## 
## attr(,"class")
## [1] "group"
plot.group(rm_proteina,main="Grupos para el factor de bloque Proteína")

Los resultados de la prueba de rango múltiple arrojan como resultado que, los niveles del factor de bloque denominado proteína que generan una producción promedio más alta son los niveles \(\varphi\) y \(\gamma\), en orden descendente, los niveles que seguirían son \(\varepsilon\) y \(\delta\), con la misma producción promedio entre ellos, para el caso del nivel \(\chi\), pertenece a dos grupos, los cuales no están claramente diferenciados junto con el factor \(\beta\), el misma factor \(\beta\) se incluye junto con el factor \(\alpha\) en el último grupo con la menor producción de leche.
Una vez obtenidos los resultados del experimento se está en condiciones de concluir que para obtener los mejores resultados en la producción de leche, debe proporsionarse al menos un 0.1% de Lisina, dado que, todos los niveles probados tienen en promedio el mismo efecto, salvo en el caso de el nivel en el que no se suministra lisina, y por lo menos un 12% de proteína.

Verificación de los supuestos del modelo.

Los supuestos del modelo son:
1. Normalidad de residuos
2. Homocedasticidad de residuos.

Normalidad de Residuos

PAra el caso de la Prueba de Normalidad, procederemos a aplicar la Prueba de Shapiro-Wilks, considerando las siguientes hipótesis estadísticas:

\[ H_0: {x \in N(\mu=0,\sigma^2=Constante)} \] \[ H_1: {x \not\in N(\mu=0,\sigma^2=Constante)} \] mediante la sigueinte secuencia de comandos:. 
normalidad=shapiro.test(resid(modelo))
print(normalidad)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(modelo)
## W = 0.98663, p-value = 0.8467
#Gráfica de Probabiliad Normal
qqnorm(resid(modelo), main= "Gráfica de Probabilidad para los Residuales del Modelo", xlab="Cuantiles Teoricos", ylab = "Cuantiles de muestra")
qqline(resid(modelo))

Como es observarse, se acepta la hipótesis nula de la Prueba de Shapiro Wilks, y se concluye que los residuales provienen de una problación con Distribución de Probabilidad Normal con \(\mu=0\), conclusion que es confirmada por la Gráfica de Probabilidad Normal.

Prueba de Homocedasticidad

La Prueba de Homocedasticidad, o de Varianzas Iguales, se realiza mediante la Prueba de Bartlett, misma que para el caso particular, se aplicará al factor de tratamiento y al factor de bloque denominados proteína.
La Prueba de Bartlett tiene las siguientes hipótesis:
\[ H_0:\sigma^2_i=\sigma^2_j=Constante \] \[ H_1:\sigma^2_i\neq\sigma^2_j\neq Constante \]
La Prueba de Bartlett se realiza con la siguiente línea de comandos: 
homocedasticidad_lisina=bartlett.test(resid(modelo)~Lisina,data=datos)
print(homocedasticidad_lisina)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid(modelo) by Lisina
## Bartlett's K-squared = 6.4864, df = 6, p-value = 0.371
homocedasticidad_proteina=bartlett.test(resid(modelo)~Proteina,data=datos)
print(homocedasticidad_proteina)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid(modelo) by Proteina
## Bartlett's K-squared = 5.7725, df = 6, p-value = 0.4492
Para este caso, se aceptan las hipótesis nulas de ambas pruebas, por lo que se concluye que el modelo matemático propuesto es adecuado para estimar la producción de leche medida en galones por las unidades experimentales, en términos de la cantidad de lisina y proteína suministradas durante el experimento.

Bibliografía

Pulido, H. G., & Vara Salazar, R. de la. (2012). Analisis y diseño de experimentos (3rd ed.). McGraw Hill.