Diseños de la libreria agricolae
Kahomy Avila 1000854408, Carlos Burgos 1001089617
10/11/2020
DISEÑO CUADRADO LATINO
Los diseños en cuadrados latinos se emplean cuando es necesario controlar dos fuentes de variabilidad. En estos diseños el número de niveles del factor principal es igual al número de niveles de las dos variables de bloque y además presenta las siguientes caracteristicas:
Cada uno de los factores tiene el mismo número de niveles, K.
Cada nivel del factor principal aparece una vez en cada fila y una vez en cada columna.
No hay interacción entre los factores.
Modelo estadistico
\[ y_{ijk}=\mu+\beta_i+\tau_j+\lambda_k+\epsilon_{ijk}\\ Codiciones~~laterales \]
Ejemplo con libreria agricolae
Se esta estudia el rendimiento presentado por 5 variedades hibridas de maiz (a,b,c,d,e), se consideran 5 hileras y 5 columnas en una parcela. se consideran 5 niveles de cada uno de estos factores. - La variable respuesta seria el rendimiento. - Las filas y las columnas son los bloques, ambos con 5 niveles. - La variedad es el factor. tiene 5 niveles. - El numero de observaciones es 25.
set.seed(20)
hilera = factor( c(rep("1",5),
rep("2",5),
rep("3",5),
rep("4",5),
rep("5"),5))
columna = factor( c(rep("1",5),
rep("2",5),
rep("3",5),
rep("4",5),
rep("5"),5))
Variedad = factor(c("a","b","c","d","e"))
rendimiento = (c(rnorm(n = 25,mean = 18,sd = 2)))Acontinuacion se realizara la aleatorizacion de los tratamientos mediante el uso de la funcion design.lsd de la libreria agricolae.
library(agricolae)## Warning: package 'agricolae' was built under R version 4.0.3diseñolatin = design.lsd(trt = Variedad,serie = 5,seed = 20)Una vez realizada la aleatorizacion procedemos a extraer el sketch y el book para tener el dataframe aleatorizado y la matriz del diseño.
book = diseñolatin$book
sk= diseñolatin$sketchAqui podemos ver la matriz del diseño sin la variable respuesta observada.
colnames(sk)=paste0("Columna", 1:5)
rownames(sk)=paste0("Hilera", 1:5)
sk## Columna1 Columna2 Columna3 Columna4 Columna5
## Hilera1 "b" "c" "e" "d" "a"
## Hilera2 "d" "e" "b" "a" "c"
## Hilera3 "c" "d" "a" "e" "b"
## Hilera4 "e" "a" "c" "b" "d"
## Hilera5 "a" "b" "d" "c" "e"Podemos observar que la letra latina no se repite ni en su misma fila ni en su misma columna.
Acontinuacion procedemos a realizar el analisis de varianzas, pero antes se deben ingresar los datos de la variable respuesta con el dataframe previamente aleatorizado.
latin= data.frame(book, rendimiento)
lm1 = lm(latin$rendimiento ~ latin$Variedad + latin$col + latin$row)
anovalatin = aov(lm1)
summary(anovalatin)## Df Sum Sq Mean Sq F value Pr(>F)
## latin$Variedad 4 1.19 0.297 0.056 0.993
## latin$col 4 12.24 3.060 0.578 0.684
## latin$row 4 15.26 3.816 0.720 0.594
## Residuals 12 63.56 5.297Observando que los p valores son mayores que el nivel de significación del 5% podemos decir que ningún efecto es significativo.
Prueba de algunos supuestos
shapiro.test(anovalatin$res)##
## Shapiro-Wilk normality test
##
## data: anovalatin$res
## W = 0.97218, p-value = 0.7007Al realizar shapiro test y obtener un p valor mayor a 0.05 en los residuales del anova podemos evidencia que se acepta la hipotesis de la normalidad de los residuos.
Acontinuacion podemos ver de manera grafica una prueba de normalidad con qqplot donde se evidencia que los datos se ajustan a la linea de normalidad.
qqnorm(anovalatin$residuals, col="blue")
qqline(anovalatin$residuals)
Evidenciamos tambien la distribucion normal mendiante el siguiente histograma.
hist(anovalatin$residuals,col = "lightblue", xlab= "Residuos" , ylab= "Densidad",freq = F)
lines(density(anovalatin$residuals),col="red",lwd=3)
DISEÑO CUADRADO GRECOLATINO
El diseño en cuadrado grecolatino se puede ver como una extensión del diseño cuadrado latino en el que se añade una tercera variable bloque. En este modelo como en el diseño en cuadrado latino, todos los factores deben tener el mismo número de niveles.
Modelo estadistico
\[ y_{ijkl}=\mu+\theta_i+\tau_j+\omega_k+\psi_i +\epsilon_{ijkl}\\ Codiciones~~laterales \] ### Ejemplo con la libreria agricolae
En la sintesis de un determinado producto agroquímico se está interesado en comparar 4 procedimientos. En esta sintesis pueden influir la temperatura, presión y tipo de catalizador utilizado, decidiéndose realizar un experimento en cuadrado greco-latino. Para ello, se consideran 4 niveles de cada uno de estos factores.
La variable respuesta seria la cantidad de agroquimico obtenida.
Los procedimientos, temperaturas y presion son los factores bloque, cada uno con cuatro niveles.
El factor principal es el catalizador, que tiene 4 niveles.
El numero de observaciones es 16.
set.seed(20)
procedimientos = factor( c(rep("p1",4),
rep("p2",4),
rep("p3",4),
rep("p4",4)))
temperatura = factor(c(rep("t1",4),
rep("t2",4),
rep("t3",4),
rep("t4",4)))
presion = factor(c("a","b","c","d"))
catalizador = factor(c("alpha", "beta", "gama", "lamda"))
cantidad = (c(rnorm(n = 16,mean = 10)))Acontinuacion se realizara la aleatorizacion de los tratamientos mediante el uso de la funcion design.lsd de la libreria agricolae.
diseñogreek <- design.graeco(trt1 = presion,trt2 = catalizador,serie = 4,seed = 20)Una vez realizada la aleatorizacion procedemos a extraer el sketch y el book para tener el dataframe aleatorizado y la matriz del diseño.
dis1 <-diseñogreek$book
sk1 <- diseñogreek$sketchAqui podemos ver la matriz del diseño sin la variable respuesta observada.
colnames(sk1)=paste0("Temperatura", 1:4)
rownames(sk1)=paste0("Procedimiento", 1:4)
sk1## Temperatura1 Temperatura2 Temperatura3 Temperatura4
## Procedimiento1 "b beta" "a alpha" "c gama" "d lamda"
## Procedimiento2 "a gama" "b lamda" "d beta" "c alpha"
## Procedimiento3 "c lamda" "d gama" "b alpha" "a beta"
## Procedimiento4 "d alpha" "c beta" "a lamda" "b gama"Podemos observar que ni la letra latina, ni la griega se repiten en su misma fila ni en su misma columna.
Acontinuacion procedemos a realizar el analisis de varianzas, pero antes se deben ingresar los datos de la variable respuesta con el dataframe previamente aleatorizado.
datgreek = data.frame(dis1,cantidad)
datgreek## plots row col presion catalizador cantidad
## 1 10001 1 1 b beta 8.559857
## 2 10002 1 2 a alpha 8.627640
## 3 10003 1 3 c gama 10.921764
## 4 10004 1 4 d lamda 10.592630
## 5 20001 2 1 a gama 10.118595
## 6 20002 2 2 b lamda 10.463420
## 7 20003 2 3 d beta 8.839000
## 8 20004 2 4 c alpha 10.214091
## 9 30001 3 1 c lamda 10.198672
## 10 30002 3 2 d gama 10.472681
## 11 30003 3 3 b alpha 8.165710
## 12 30004 3 4 a beta 8.902460
## 13 40001 4 1 d alpha 10.398615
## 14 40002 4 2 c beta 9.534995
## 15 40003 4 3 a lamda 10.553053
## 16 40004 4 4 b gama 9.128063Dado que los tratamientos aleatorizados anteriormente tienen nombre row y col se procede a extraer la aleatorizacion bajo otros nombres que representen temperatura y procedimiento. Y a finalmente poner esto en un dataframe completo.
proced = datgreek$row
temp = datgreek$col
datafinal= data.frame(datgreek$plots, proced, temp, datgreek$presion,datgreek$catalizador,datgreek$cantidad)
datafinal## datgreek.plots proced temp datgreek.presion datgreek.catalizador
## 1 10001 1 1 b beta
## 2 10002 1 2 a alpha
## 3 10003 1 3 c gama
## 4 10004 1 4 d lamda
## 5 20001 2 1 a gama
## 6 20002 2 2 b lamda
## 7 20003 2 3 d beta
## 8 20004 2 4 c alpha
## 9 30001 3 1 c lamda
## 10 30002 3 2 d gama
## 11 30003 3 3 b alpha
## 12 30004 3 4 a beta
## 13 40001 4 1 d alpha
## 14 40002 4 2 c beta
## 15 40003 4 3 a lamda
## 16 40004 4 4 b gama
## datgreek.cantidad
## 1 8.559857
## 2 8.627640
## 3 10.921764
## 4 10.592630
## 5 10.118595
## 6 10.463420
## 7 8.839000
## 8 10.214091
## 9 10.198672
## 10 10.472681
## 11 8.165710
## 12 8.902460
## 13 10.398615
## 14 9.534995
## 15 10.553053
## 16 9.128063Acontinuacion procedemos a realizar el analisis de varianzas.
modeloGreco= lm(datafinal$datgreek.cantidad ~ datafinal$datgreek.presion + datafinal$proced + datafinal$temp + datafinal$datgreek.catalizador, datafinal)
anovaGreco=aov(modeloGreco)
summary(anovaGreco)## Df Sum Sq Mean Sq F value Pr(>F)
## datafinal$datgreek.presion 3 3.251 1.0837 1.499 0.374
## datafinal$proced 3 0.609 0.2029 0.281 0.838
## datafinal$temp 3 0.090 0.0299 0.041 0.987
## datafinal$datgreek.catalizador 3 5.776 1.9252 2.663 0.221
## Residuals 3 2.169 0.7228Observando los valores de los p valores 0.374 , 0.838 , 0.987, 0.221 son mayores que el nivel de significación del 5% por lo que podemos decir que ningún efecto es significativo.
Prueba de algunos supuestos
shapiro.test(anovaGreco$res) ##
## Shapiro-Wilk normality test
##
## data: anovaGreco$res
## W = 0.85127, p-value = 0.01422Al realizar shapiro test y obtener un p valor de 0.01 en los residuales del anova podemos evidencia que se rechaza la hipotesis de la normalidad de los residuos.
summary(anovaGreco$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.53236 -0.22888 0.04575 0.00000 0.27463 0.44086Se evidencia media= 0, lo que supone la normalidad de los datos.
DISEÑO EN BLOQUES COMPLETOS ALEATORIZADOS
La palabra bloque en el nombre hace referencia al hecho de que se ha agrupado a las unidades experimentales en función de alguna variable extraña, lo aleatorizado hace referencia al hecho de que los tratamientos se asignan aleatoriamente dentro de los bloques y la completes implica que se utiliza cada tratamiento exactamente una vez dentro de cada bloque. Estos diseños se caracterizan tambien porque los efectos bloque y tratamiento son aditivos; es decir no hay interacción entre los bloques y los tratamientos. Los diseños en bloques completos aleatorizados Todos los tratamientos se prueban en cada bloque exactamente vez.
Modelo estadistico
\[ y_{ij}=\mu+\tau_i+\beta+\epsilon_{ij} ,\\ i=1,...,3;\\ j=1,...10.\\En \ general\\ i=1,...I;\\ j=1,...,J \]
Ejemplo (tomado de pagina de estadistica de la Universidad de Granada, practica 7 del curso de diseño de experimentos a cargo de Ana Maria Lara Porras)
El Abeto Abies alba es un arbol importante del cual se extraen compuestos para perfumeria. En estos ultimos años se ha observado que la producción de semillas de estos ha disminuido y con el objetivo de lograr buenas producciones se proponen tres tratamientos. Se observa que arboles diferentes tienen distintas caracteristicas naturales de reproduccion, este efecto de las diferencias entre los arboles se debe de controlar y este control se realiza mediante bloques. En el experimento se utilizan 10 abetos, dentro de cada abeto se seleccionan tres ramas semejantes. Cada rama recibe exactamente uno de los tres tratamientos que son asignados aleatoriamente. Constituyendo cada arbol un bloque completo. No hay ningun otro factor que pueda afectar de forma significativa a los resultados
La variable respuesta es el numero de semillas
El factor es el Tratamiento que tiene tres niveles.
El Bloque es el Abeto que tiene diez niveles. Es un factor de efectos fijos ya que viene decidido que niveles concretos se van a utilizar.
Los tres tratamientos se prueban en cada bloque exactamente una vez.
El Tamaño del experimento es de 30 observaciones.
Acontinuacion visualizamos los datos donde y representa el numero de semillas.
library(data.table)## Warning: package 'data.table' was built under R version 4.0.3library(curl)## Warning: package 'curl' was built under R version 4.0.3datos_semillas = fread("https://wpd.ugr.es/~bioestad/wp-content/uploads/supuesto3.txt")
datos_semillas## y Tratamiento Abeto
## 1: 7 1 1
## 2: 9 2 1
## 3: 10 3 1
## 4: 8 1 2
## 5: 9 2 2
## 6: 10 3 2
## 7: 9 1 3
## 8: 9 2 3
## 9: 12 3 3
## 10: 10 1 4
## 11: 9 2 4
## 12: 12 3 4
## 13: 11 1 5
## 14: 12 2 5
## 15: 14 3 5
## 16: 8 1 6
## 17: 10 2 6
## 18: 9 3 6
## 19: 7 1 7
## 20: 8 2 7
## 21: 7 3 7
## 22: 8 1 8
## 23: 8 2 8
## 24: 7 3 8
## 25: 7 1 9
## 26: 9 2 9
## 27: 10 3 9
## 28: 8 1 10
## 29: 9 2 10
## 30: 10 3 10
## y Tratamiento AbetoAcontinuacion imputamos los datos del dataframe la carateristica factor para su correcto uso mas adelante.
datos_semillas$Tratamiento = factor(datos_semillas$Tratamiento)
datos_semillas$Tratamiento## [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## Levels: 1 2 3datos_semillas$Abeto = factor(datos_semillas$Abeto)
datos_semillas$Abeto## [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9
## [26] 9 9 10 10 10
## Levels: 1 2 3 4 5 6 7 8 9 10observamos los 10 niveles mencionados en el planteamiento.
Para calcular la tabla ANOVA primero hacemos uso de la función aov de la forma como lo establece el modelo estadistico como se ve acontinuacion.
mod = aov(y ~ Tratamiento + Abeto, data = datos_semillas)
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 2 16.2 8.100 9.228 0.00174 **
## Abeto 9 54.8 6.089 6.937 0.00026 ***
## Residuals 18 15.8 0.878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tratamiento muestra un p-valor de 0.001, es menor que el nivel de significación del 5%, por lo que se rechaza la Hipótesis nula de igualdad de tratamientos.
Abeto muestra un p-valor de 0.0002, menor que el nivel de significación del 5%, por lo que se rechaza la Hipótesis nula de igualdad de bloques. Y adicionalmente muestra un valor grande de F de los bloques de 6.937 lo que implica que el factor bloque tiene un efecto grande.
Verificacion de algunos supuestos
La interaccion entre el factor bloque y los tratamientos se estudiara mediante el Test de Tukey. Urilizando la libreria daewr la cual es utilizada en el libro “Design and Analysis of Experiments with R”
library(daewr)## Warning: package 'daewr' was built under R version 4.0.3## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.designTukey1df(datos_semillas)## Source df SS MS F Pr>F
## A 2 16.2 8.1
## B 9 54.8 6.0889
## Error 18 15.8 71.1
## NonAdditivity 1 3.5573 3.5573 4.94 0.0401
## Residual 17 12.2427 0.7202analizando el tukey test podemos observar que no hay interacción entre los tratamientos aplicados y los abetos.
Hipótesis de Normalidad
shapiro.test(mod$residuals)##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.96415, p-value = 0.3935Podemos observar un p-valor de 0.3935 que aceptaría la hipótesis de normalidad por ser mayor al 5%.
Graficamente lo podemos ver mediante el siguiente qqplot
qqnorm(mod$residuals, col="blue")
qqline(mod$residuals)
Hipotesis de Homogeneidad de Varianzas
bartlett.test(datos_semillas$y, datos_semillas$Tratamiento)##
## Bartlett test of homogeneity of variances
##
## data: datos_semillas$y and datos_semillas$Tratamiento
## Bartlett's K-squared = 4.1729, df = 2, p-value = 0.1241El p-valor de 0.1241 es mayor al nivel significacion por lo que no podemos rechazar la hipótesis de igualdad de varianzas en el factor principal.
bartlett.test(datos_semillas$y, datos_semillas$Abeto)##
## Bartlett test of homogeneity of variances
##
## data: datos_semillas$y and datos_semillas$Abeto
## Bartlett's K-squared = 4.0723, df = 9, p-value = 0.9066El p-valor es mayor que 0.05 por lo que no podemos rechazar la hipótesis de igualdad de varianzas en el factor bloque.
Se rechaza la Hipótesis nula de igualdad de tratamientos. Por lo que en conclusion se puede decir que existen diferencias significativas en el numero de semillas entre los tres tratamientos.
DISEÑO EN BLOQUES INCOMPLETOS ALEATORIZADOS
El diseño de bloques incompletos aleatorizados se da cuando no se puede realizar todos los tratamientos en cada bloque.
Modelo estadistico
\[ y_{ij}=\mu+\tau_i+\epsilon_{ij}\] ### Ejemplo
Se evalua la efectividad en el retraso del crecimiento de colonias de hongos utilizando cuatro soluciones diferentes para lavar las botas a la entrada de un invernadero. El análisis se realiza en el laboratorio y sólo se pueden realizar seis pruebas en un mismo dia. Como los días son una fuente de variabilidad potencial, el investigador decide utilizar un diseño aleatorizado por bloques, pero al recopilar las observaciones durante seis dias no ha sido posible aplicar todos los tratamientos en cada dia, sino que sslo se han podido aplicar dos de las cuatro soluciones cada dia. El objetivo principal es estudiar la efectividad en el retraso del crecimiento de colonia de hongos utilizando cuatro soluciones, por lo que se trata de un factor con cuatro niveles. Sin embrago, como los días son una fuente de variabilidad potencial, consideramos un factor bloque con seis niveles.
La variable respuesta es el numero de colonia de hongos.
El Factor es las Soluciones que tiene cuatro niveles.
El factores de bloque son los días que tiene seis niveles.
El Modelo incompleto ya que Todos los tratamientos no se prueban en cada bloque.
El Tamaño del experimento es de 12 observaciones.
no_hongos = c(12,24,31,21,20,21,19,18,19,15,19,47)
soluciones = c(1,1,1,2,2,2,3,3,3,4,4,4)
dias = c(1,2,3,1,5,6,3,4,6,2,4,5)
datos_hongos = data.frame(no_hongos,soluciones,dias)
datos_hongos## no_hongos soluciones dias
## 1 12 1 1
## 2 24 1 2
## 3 31 1 3
## 4 21 2 1
## 5 20 2 5
## 6 21 2 6
## 7 19 3 3
## 8 18 3 4
## 9 19 3 6
## 10 15 4 2
## 11 19 4 4
## 12 47 4 5library(daewr)
library(AlgDesign)## Warning: package 'AlgDesign' was built under R version 4.0.3Usaremos la función “BIBsize(t , k)” de la librería daewr la cual nos permite saber si el diseño puede realizarse. Calcula los parámetros del diseño donde
t = número de niveles del factor tratamiento.
k = número de tratamientos por bloque.
BIBsize(t = 4 , k = 2)## Posible BIB design with b= 6 and r= 3 lambda= 1Para poder evaluar el efecto de los tratamientos la suma de cuadrados de tratamientos deben ajustarse por bloques, por lo tanto primero se introducen los bloques y después los tratamientos.
Acontinuacion calculamos el ANOVA haciendo uso de la función aov.
mod3 <- aov(no_hongos ~ dias + soluciones, data = datos_hongos )
summary(mod3)## Df Sum Sq Mean Sq F value Pr(>F)
## dias 1 80.3 80.26 0.876 0.374
## soluciones 1 2.6 2.63 0.029 0.869
## Residuals 9 824.8 91.64Al observar el Analisis de varianzas soluciones tiene un p valor mayor al nivel de significacion del 5%, se rechaza la hipotesis nula que los tratamientos sean iguales y por lo tanto podemos concluir que la solucion no influye en el retraso del crecimiento de las colonias de hongos.
Para evaluar el efecto de los bloques la suma de cuadrados de bloques debe ajustarse por los tratamientos, por lo tanto primero se introducen los tratamientos y después los bloques.
mod4 <- aov(no_hongos ~ soluciones + dias, data = datos_hongos )
summary(mod4)## Df Sum Sq Mean Sq F value Pr(>F)
## soluciones 1 21.6 21.60 0.236 0.639
## dias 1 61.3 61.29 0.669 0.435
## Residuals 9 824.8 91.64Al observar el Analisis de varianzas soluciones tiene un p valor mayor al nivel de significacion del 5%, se rechaza la hipotesis nula que los tratamientos sean iguales y por lo tanto podemos concluir que los dias en los que se realiza la prueba no influye en el retraso del crecimiento de las colonias de hongos
DISEÑO CUADRADOS DE YOUDEN
puede suceder que el número de niveles disponibles de uno de los factores de control sea menor que el número de tratamientos, en este caso estariamos ante un diseño en cuadrado latino incompleto. Estos diseños fueron estudiados por W.J. Youden y se conocen con el nombre de cuadrados de Youden.En general un cuadrado de Youden es un diseño balanceado por bloques incompletos
Para aleatorizar con la libreria agricolae:
\[design.youden(trt, r, serie = 2, seed = 0, kinds = “Super-Duper”,first=TRUE ,randomization=TRUE)\]
Donde
trt: Tratamientos
r: replicas o número de columnas
serie: Número del gráfico 1: 11,12; 2: 101,102; 3: 1001,1002
seed: semilla
kinds: Método para aleatorizar
first: TRUE o FALSE - aleatorizar rep 1
randomization: TRUE or FALSE = aleatorizar
library(agricolae)
varieties<-c("perricholi","yungay","maria bonita","tomasa")
r<-3
outdesign <-design.youden(varieties,r,serie=2,seed=23)
youden <- outdesign$book
print(outdesign$sketch)## [,1] [,2] [,3]
## [1,] "maria bonita" "tomasa" "perricholi"
## [2,] "yungay" "maria bonita" "tomasa"
## [3,] "perricholi" "yungay" "maria bonita"
## [4,] "tomasa" "perricholi" "yungay"plots <-as.numeric(youden[,1])
print(matrix(plots,byrow=TRUE,ncol=r))## [,1] [,2] [,3]
## [1,] 101 102 103
## [2,] 201 202 203
## [3,] 301 302 303
## [4,] 401 402 403print(youden) # field book.## plots row col varieties
## 1 101 1 1 maria bonita
## 2 102 1 2 tomasa
## 3 103 1 3 perricholi
## 4 201 2 1 yungay
## 5 202 2 2 maria bonita
## 6 203 2 3 tomasa
## 7 301 3 1 perricholi
## 8 302 3 2 yungay
## 9 303 3 3 maria bonita
## 10 401 4 1 tomasa
## 11 402 4 2 perricholi
## 12 403 4 3 yungayEjemplo
En este experimento se desea conocer el rendimiento de la semilla de trigo en kg/ha, en el que se está interesado en estudiar 4 tipos de semillas y se desea eliminar estadísticamente el efecto del tipo de insecticida y abono. Se dispone de 3 tipos de abono. Para realizar este experimento se decidió utilizar un cuadrado de Youden con 4 filas, los tipos de insecticidas (i1, i2, i3 ,i4), 3 columnas, los tipos de abono (a1, a2, a3) y 4 tipos de semillas (A, B, C, D);(ya los datos se encuntran aleatorizados).
library(readxl)## Warning: package 'readxl' was built under R version 4.0.3Datos <- read_excel("C:/Users/USER/Desktop/Kahomy/Datos.xlsx")d1 = data.frame(Datos); d1## rendimiento insecticida abono semilla
## 1 26 d1 a1 A
## 2 18 d2 a1 B
## 3 19 d3 a1 C
## 4 21 d4 a1 D
## 5 25 d1 a2 B
## 6 15 d2 a2 C
## 7 25 d3 a2 D
## 8 25 d4 a2 A
## 9 16 d1 a3 C
## 10 21 d2 a3 D
## 11 24 d3 a3 A
## 12 20 d4 a3 BEl modelo estadistico para este cuadro de Youden es el siguiente
\[ y_{ijk}=\mu +\tau_i+\beta_j+\gamma_k+\epsilon_{ijk}\\ i = 1,2,3,4\\ j = 1,2,3,4\\ k = 1,2,3 \\ Condiciones~~laterales\] Realizar el análisis de varianzas para cada factor (principal y de bloqueo)
d1$insecticida = factor(d1$insecticida); d1$insecticida## [1] d1 d2 d3 d4 d1 d2 d3 d4 d1 d2 d3 d4
## Levels: d1 d2 d3 d4d1$abono = factor(d1$abono); d1$abono## [1] a1 a1 a1 a1 a2 a2 a2 a2 a3 a3 a3 a3
## Levels: a1 a2 a3d1$semilla = factor(d1$semilla); d1$semilla## [1] A B C D B C D A C D A B
## Levels: A B C Dm1 = aov (rendimiento~insecticida+abono+semilla, data = d1)
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## insecticida 3 42.92 14.31 5.202 0.1044
## abono 2 10.50 5.25 1.909 0.2919
## semilla 3 94.58 31.53 11.465 0.0376 *
## Residuals 3 8.25 2.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2 = aov (rendimiento~semilla+abono+insecticida, data = d1)
summary(m2)## Df Sum Sq Mean Sq F value Pr(>F)
## semilla 3 108.92 36.31 13.202 0.031 *
## abono 2 10.50 5.25 1.909 0.292
## insecticida 3 28.58 9.53 3.465 0.167
## Residuals 3 8.25 2.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m3 = aov (rendimiento~semilla+ insecticida+abono, data = d1)
summary(m3)## Df Sum Sq Mean Sq F value Pr(>F)
## semilla 3 108.92 36.31 13.202 0.031 *
## insecticida 3 28.58 9.53 3.465 0.167
## abono 2 10.50 5.25 1.909 0.292
## Residuals 3 8.25 2.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mediante el ANOVA realizado para los factores de bloqueo como el abono y el insecticida, se puede decir que los efectos de la semilla sobre los tratamientos son significatvos, independiente del factor que se use como bloqueo.
ANOVA del modelo
mod=lm(aov(rendimiento~semilla+insecticida+abono, data=d1))
summary(mod)##
## Call:
## lm(formula = aov(rendimiento ~ semilla + insecticida + abono,
## data = d1))
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10 11
## 0.375 -1.000 1.375 -0.750 0.625 -0.375 -0.625 0.375 -1.000 1.375 -0.750
## 12
## 0.375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.625 1.436 17.843 0.000384 ***
## semillaB -2.750 1.436 -1.915 0.151401
## semillaC -7.875 1.436 -5.483 0.011929 *
## semillaD -1.375 1.436 -0.957 0.408983
## insecticidad2 -3.875 1.436 -2.698 0.073898 .
## insecticidad3 -0.125 1.436 -0.087 0.936125
## insecticidad4 -2.500 1.436 -1.741 0.180095
## abonoa2 1.500 1.173 1.279 0.290795
## abonoa3 -0.750 1.173 -0.640 0.567924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.658 on 3 degrees of freedom
## Multiple R-squared: 0.9472, Adjusted R-squared: 0.8064
## F-statistic: 6.727 on 8 and 3 DF, p-value: 0.07233Mediante el ANOVA del modelo y los ANOVA de los factores, podemos confirmar que los efectos de la semilla sobre los tratamientos fueron significativos, mientras que el abono y el insecticida no fueron de relevancia, por lo que fue bueno bloquear con ellos.
Revisión de supuestos
Normalidad de residuales
qqnorm(mod$residuals, col="purple")
qqline(mod$residuals)
hist(mod$residuals, col="lightgreen", main = "Histograma de los Residuos", freq = F, xlab="Residuos",ylab="Densidad")
lines(density(mod$residuals), col="red", lwd=3)
norm = shapiro.test(mod$residuals)
ifelse(norm$p.value>0.05, "Los residuales son normales", "Los residuales NO son normales" )## [1] "Los residuales son normales"Varianzas de los trataminetos estadísticamente iguales
boxplot(mod$residuals~d1$semilla, border = c("blue", "green", "red", "orange"), col= "white", xlab = "Semilla", ylab = "Residuales")
vari = bartlett.test(mod$residuals, d1$semilla)
ifelse(vari$p.value>0.05, "Los tratamientos tiene una varianza igual", "Los tratamientos NO tiene una varianza igual")## [1] "Los tratamientos tiene una varianza igual"dt = duncan.test(mod, "semilla");dt## $statistics
## MSerror Df Mean CV
## 2.75 3 21.25 7.803823
##
## $parameters
## test name.t ntr alpha
## Duncan semilla 4 0.05
##
## $duncan
## Table CriticalRange
## 2 4.500659 4.309053
## 3 4.515652 4.323407
## 4 4.472854 4.282431
##
## $means
## rendimiento std r Min Max Q25 Q50 Q75
## A 25.00000 1.000000 3 24 26 24.5 25 25.5
## B 21.00000 3.605551 3 18 25 19.0 20 22.5
## C 16.66667 2.081666 3 15 19 15.5 16 17.5
## D 22.33333 2.309401 3 21 25 21.0 21 23.0
##
## $comparison
## NULL
##
## $groups
## rendimiento groups
## A 25.00000 a
## D 22.33333 a
## B 21.00000 a
## C 16.66667 b
##
## attr(,"class")
## [1] "group"Con la prueba Duncan test determinamos que con la semilla A se obtuvo el mejor rendimineto con 25 kg/ha, y el menor rendimiento con la semilla C 16.7 kg/ha, por esto la mejor semilla para la producción de trigo es la semilla A.
DISEÑO BLOQUES DIVIDIDOS (STRIP)
Diseño en parcelas en franja o bloques divididos, se utilizan cuando hay dos tipos de tratamientos (factores) y se aplican por separado en parcelas grandes, llamado franjas, en una dirección vertical y horizontal del bloque, cada boque constituye una repetición. Para generar un diseño de bloques divididos con agricolae, se requiere:
\[design.strip(trt1, trt2,r, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE) \] Donde:
trt1 = filas tratamientos
trt2 = columnas tratamientos
r = repeticiones
serie number plot, 1: 11,12; 2: 101,102; 3: 1001,1002
seed = semilla para aleatorización
kinds = método de aleatorización
randomization = TRUE or FALSE - aleatorización
Para el análisis de varianza se debe usar
\[strip.plot(BLOCK, COL,ROW,Y)\] Donde:
Block = bloqueo
COL = columna de tratamiento
ROW = Fila de tratamiento
Y = variable respuesta
Ejemplo
Este ejemplo se compone de tres variedades y cuatro niveles de fertilizantes diferentes. Donde las variedades se mantienen en la parcela principal y los niveles de fertilizante en las subparcelas, mediante el análisis de los datos, indique si las variedades y los niveles de fertilización influyen en en rendimiento. (los datos ya se encuentran aleatorizados)
library(agricolae)
library(readxl)
strip <- read_excel("C:/Users/USER/Desktop/Kahomy/strip.xlsx")
strip## # A tibble: 48 x 4
## bloque variedad fertilizante rendimiento
## <dbl> <chr> <chr> <dbl>
## 1 1 V1 F1 10.2
## 2 1 V1 F2 11.1
## 3 1 V1 F3 6.8
## 4 1 V1 F4 5.3
## 5 1 V2 F1 8
## 6 1 V2 F2 9.7
## 7 1 V2 F3 8.6
## 8 1 V2 F4 3.4
## 9 1 V3 F1 2
## 10 1 V3 F2 10.9
## # ... with 38 more rowsdebemos transformar las variables fertilizante y variedad en factor, para su análisis
strip$variedad = factor(strip$variedad); strip$variedad## [1] V1 V1 V1 V1 V2 V2 V2 V2 V3 V3 V3 V3 V1 V1 V1 V1 V2 V2 V2 V2 V3 V3 V3 V3 V1
## [26] V1 V1 V1 V2 V2 V2 V2 V3 V3 V3 V3 V1 V1 V1 V1 V2 V2 V2 V2 V3 V3 V3 V3
## Levels: V1 V2 V3strip$fertilizante = factor(strip$fertilizante);strip$fertilizante## [1] F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1
## [26] F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4 F1 F2 F3 F4
## Levels: F1 F2 F3 F4modelo
\[ y_{ijk} = \mu+\alpha_i+\delta_k+\eta_{ik}+\beta_j+\alpha\beta_{ij}+ \epsilon_{ijk}\\ i = 1,2,...I\\ j = 1,2,...J\\ k = 1,2,...K \\ Condiciones~~laterales \]
El análisis de varianza de un diseño de parcela de franjas se divide en tres partes:
- El análisis de factor horizontal
- El análisis de factor vertical
- El análisis de interacción Para ellos agricolae usa la función strip.plot(), que ajusta el modelo de análisis de varianza.
model = strip.plot(BLOCK = strip$bloque,
COL = strip$variedad,
ROW = strip$fertilizante,
Y = strip$rendimiento)##
## ANALYSIS STRIP PLOT: strip$rendimiento
## Class level information
##
## strip$variedad : V1 V2 V3
## strip$fertilizante : F1 F2 F3 F4
## strip$bloque : 1 2 3 4
##
## Number of observations: 48
##
## model Y: strip$rendimiento ~ strip$bloque + strip$variedad + Ea + strip$fertilizante + Eb + strip$fertilizante:strip$variedad + Ec
##
## Analysis of Variance Table
##
## Response: strip$rendimiento
## Df Sum Sq Mean Sq F value Pr(>F)
## strip$bloque 3 13.692 4.564 2.4086 0.1007122
## strip$variedad 2 163.007 81.503 57.0397 0.0001248 ***
## Ea 6 8.573 1.429 0.7541 0.6144958
## strip$fertilizante 3 152.685 50.895 17.1086 0.0004638 ***
## Eb 9 26.773 2.975 1.5700 0.1983665
## strip$fertilizante:strip$variedad 6 40.320 6.720 3.5465 0.0170071 *
## Ec 18 34.107 1.895
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 16 %, cv(b) = 23 %, cv(c) = 18.4 %, Mean = 7.491667Siendo el rendimiento la variable respuesta, vemos que el p-valor para la variedad y el fertilizante es menor a 0.05, indicando así que estos afectan signifiativamente al rendimineto.
DISEÑO FACTORIAL
Un diseño factorial es aquél en el que se investigan todas las posibles combinaciones de los niveles de los factores en cada ensayo completo. En este caso se dicen que están cruzados, apareciendo el concepto de interacción.Se supone la existencia de repeticiones del experimento en cada una de las posibles combinaciones de los niveles del factor correspondiente.
Para aleatorizar
\[design.ab(trt, r, serie = 2, design=c(“rcbd”,“crd”,“lsd”), seed = 0, kinds = “Super-Duper”,first=TRUE,randomization=TRUE)\]
Donde:
trt: n niveles del factor
r: replicaciones o bloques
serie: número del gráfico, 1: 11,12; 2: 101,102; 3: 1001,1002
design: tipo de diseño
seed: semilla
kinds: Método de aleatorizar
first: TRUE o FALSE = aleatoriza rep 1
randomization: TRUE o FALSE = aleatorizar
Ejemplo
Con la libreria agricolae se generan y aleatorizan los datos para un factorial 3x2 con 4 bloques El factor A es definido por 3 sustancias de crecimiento (A) a 2 niveles de concentración (B) aplicados en un experimento, con tres repeticones para cad uno, para evaluar la propagación vegetativa de un cultivo sobre medios artificiales. La formación de callos se medirá a la cuarta semana.
library(agricolae)
library(readxl)
factorial <- read_excel("C:/Users/USER/Desktop/Kahomy/factorial.xlsx")Se transforman los vectores en factores, para su análisis
factorial$A = factor(factorial$A)
factorial$B = factor(factorial$B)
factorial$bloque = factor(factorial$bloque)Modelo
\[y_{ijk} = \mu+\alpha_i+\delta_k+\beta_j+\alpha\beta_{ij}+ \epsilon_{ijk}\\ i = 1,2\\ j = 1,2\\ k = 1,2,3,4\\ Condiciones~~laterales \]
modelof <- aov(Y~bloque+A+B+A:B,data=factorial)
anova(modelof) ## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## bloque 3 73.125 24.3750 6.0811 0.006429 **
## A 1 7.042 7.0417 1.7568 0.204864
## B 2 38.583 19.2917 4.8129 0.024281 *
## A:B 2 2.083 1.0417 0.2599 0.774549
## Residuals 15 60.125 4.0083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La respuesta, formación de callos, por medio del análisis de varianzas vemos que el facor bloque y nivel de concentración (B), afectan significativamente la formación de callos.
Revisión de supuestos
Normalidad de los residuales
hist(modelof$residuals, col=c("lightblue", "yellow", "purple"), main = "Histograma de los Residuos", freq = F, xlab="Residuos",ylab="Densidad")
lines(density(modelof$residuals), col="black", lwd=3)
qqnorm(modelof$residuals, col="orange4")
qqline(modelof$residuals)
ps = shapiro.test(modelof$residuals)
ifelse(ps$p.value>0.05, "Los residuales son normales", "Los residuales NO son normales" )## [1] "Los residuales son normales"Igualdad de varianzas
boxplot(modelof$residuals~factorial$bloque, border = c("blue", "green", "red", "orange"), col= "white", ylab = "Residuales")
varif = bartlett.test(modelof$residuals, factorial$bloque)
ifelse(varif$p.value>0.05, "Los tratamientos tiene una varianza igual", "Los tratamientos NO tiene una varianza igual")## [1] "Los tratamientos tiene una varianza igual"dtf = duncan.test(modelof, "bloque");dtf## $statistics
## MSerror Df Mean CV
## 4.008333 15 6.958333 28.77244
##
## $parameters
## test name.t ntr alpha
## Duncan bloque 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.014325 2.463748
## 3 3.159826 2.582673
## 4 3.250248 2.656579
##
## $means
## Y std r Min Max Q25 Q50 Q75
## 1 4.500000 2.345208 6 2 8 2.50 4.5 5.75
## 2 6.500000 1.224745 6 5 8 6.00 6.0 7.50
## 3 7.500000 2.664583 6 4 12 6.25 7.5 8.00
## 4 9.333333 2.732520 6 6 14 8.00 9.0 10.00
##
## $comparison
## NULL
##
## $groups
## Y groups
## 4 9.333333 a
## 3 7.500000 ab
## 2 6.500000 bc
## 1 4.500000 c
##
## attr(,"class")
## [1] "group"Con la prueba Duncan test determinamos que con el bloque 4 se obtuvo la mayor formación de callo y la menor formación con la semilla el bloque 1.
Punto 8. Resúmenes
- El uso de los diseños en parcelas divididas:
http://207.67.83.164/quality-progress/2007/10/laboratory/when-should-you-consider-a- split-plot-design.html
Los diseños en parcelas divididas tienen su origen en entornos agrícolas, es por esto que gran parte de la nomenclatura se refiere a parcelas de tierra. Estos diseños resultan de gran utilidad cuando se tienen situaciones donde algunos factores tienen niveles más difíciles de cambiar (factores de parcela completa) que otros (factores de subparcela), Como por ejemplo el ejemplo que se menciona en el informe donde un experimento de campo agrícola en el que varias parcelas de tierra adyacentes reciben el mismo fertilizante (el factor de parcela completa) y luego el tipo de semilla en particular (factor de subparcela) se asigna al azar y se aplica a parcelas dentro de ese grupo más grande. Allí el número de parcelas completas va a ser igual a el número de veces que se restablecen los factores de la parcela completa; mientras que, dado que se restablecen los factores de la subparcela después de cada ejecución, el número total de experimentos va a ser igual al número de subparcelas.
Es preciso para un buena ejecución de un experimento bajo el diseño en parcelas tener especial cuidado con el análisis realizado, pues puede darse el caso que se ha realizado un experimento de parcela dividida, pero al momento de analizar se tratan las corridas como independientes (como si fuera un diseño completamente al azar) y esto supone un error que podría conducir a falsas conclusiones sobre los efectos de los factores. A este efecto se le conoce como una parcela dividida inadvertida. También resulta necesario al analizar este diseño tener en cuenta que al no reiniciar los factores de parcela completa entre cada ejecución, las observaciones obtenidas estarán más relacionadas entre sí que cuando se reinician todos los factores.
El uso de diseños en parcela dividida resulta apropiado cuando no se quiere reiniciar todos los factores entre cada ejecución, esto no implica necesariamente no realizar una aleatorización, si no por el contrario que esta sea más estructurada. Adicionalmente, otra ventaja de usar diseño en parcelas divididas es que suponen una reducción de costos del experimento debido a que como se mencionaba anteriormente no se realizan tantos cambios sobre los factores de parcela completa y los subfactores. Esto a la vez supone una ventaja indirecta, y es que al reducir costos en la ejecución del experimento, esos recursos ahorrados se pueden invertir en instrumentos de mayor precisión para el levantamiento de los datos en las observaciones o en un mejor análisis de los datos.
- Sobre lo que significa unidad experimental y unidad de observación
Short communication: On recognizing the proper experimental unit in animal studies in the dairy sciences (https://www.sciencedirect.com/science/article/pii/S002203021630621X)
https://online.stat.psu.edu/stat502/lesson/6/6.1-0
La unidad experimental o unidad de replicación es la entidad más pequeña que se asigna independientemente de todas las demás unidades de un tratamiento particular, en general en estadística se dice que no presentan diferencias significativas entre ellas y que las inferencias que se hagan sobre ellas son válidas y confiables, indiferenciadamente del tratamiento que se le asigno a cada una.
Pensemos en que se quiere saber el efecto de un tratamiento para el control de la mastitis en las vacas, este tratamiento es inyectado individualmente, donde cada vaca es una unidad experimental, y no importa que estas estén juntas o interaccionando, como en dos corrales, en donde se distribuyan todas vacas, pero que pasaría si en vez de un tratamiento que contrala la mastitis se quiera ver el efecto de dos dietas diferentes frente la ganancia de masa corporal, de modo que se alimentaran de la misma a todas las vacas que se encuentren en el mismo corral.
La unidad de observación también conocida como unidad de muestreo se define como la entidad física sobre la que se mide un resultado de interés en un experimento. Por lo que en el ejemplo anterior la unidad de observación seria cada vaca, y la unidad experimental seria cada corral cada uno con n vacas.
Pero hay que tener cuidado con la asignación de unidad experimental y unidad de observación, esto depende de lo que se quiere estudiar, es decir de lo que se va a medir, en el ejemplo anterior si ahora se midieron individualmente la producción de leche con respecto a la dieta implantada, se encuentra un sesgo o desajuste natural entre la unidad independiente o experimental (corral) y la unidad de observación (vaca), en este caso se presenta una estructura de diseño anidada, es decir, donde se debe tener en cuenta la correlación (falta de independencia) que puedan tener las vacas como individuo con respecto a la toma de datos sobre el coral como conjunto.
Se denomina submuestra, psudorreplicas o replicas técnicas cuando las unidades de observación se encuentran anidadas en la unidad experimental, es importante determinar la jerarquía de los datos y distinguir entre el error de muestreo (unidades de observación) y el error experimental (unidades experimentales).
Podemos inferir que las unidades experimentales se definen en términos de asignaciones de tratamiento independientes, mientras que las unidades de observación se definen en términos de medidas de resultados, donde para obtener un ANOVA correcto los valores promedio de las unidades de muestreo deben calcularse para cada unidad experimental.
- Guía para diseñar experimentos exitosos
https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2013.0114
Diseñar experimentos une a las matemáticas con la biología, donde lo fundamental son las preguntas y / o hipótesis, seguido del ciclo de retroalimentación en donde el investigador puede crear nuevas hipótesis científicas y modificaciones del diseño para eventos futuros. El diseño experimental tiene cuatro pilares y cada uno de estos pilares tiene un impacto profundo en el diseño experimental, el análisis de datos y las conclusiones que resultan de la conducta experimenta, la replicación, aleatorización, el bloqueo y las unidades experimentales son estos pilares.
REPLICACIÓN: un mecanismo para estimar el error experimental, necesario para pruebas de hipótesis validas, estimar el ruido en las estimaciones de los efectos del tratamiento; aumenta la precisión del experimento, esto basado en la formula clásica del error estándar; la inferencia del experimento tiene un mayor alcance, pues tiene un rango más alto de observaciones; el investigador tiene el control sobre el error, permitiendo así una mayor potencia en el experimento. El hacer repeticiones tiene un efecto directo en el experimento, es el tema más importante y estas pueden darse en cuatro niveles dentro del experimento: en la unidad experimental, en el experimento completo, en muestreo dentro de las unidades experimentales o en medidas repetidas. Diseños aumentados implican cientos o miles de tratamientos no repetidos organizados en bloques pequeños (fincas pequeñas). Aumentando el grado de complejidad, el diseño se vuelve más eficiente, los diseños entre más simples son menos grandes.
ALEATORIZACIÓN: implica que cada tratamiento se aplique a cada unidad experimental y su uso principal es para minimizar el efecto de otras variables que no son las de interés, aumenta la posibilidad de que el supuesto de independencia de los errores se cumpla, es decir que se encuentren distribuidos independientemente.
BLOQUEO: se utiliza con uno o ambos propósitos, por precisión o por conveniencia, cuando existe un factor que es difícil de aleatorizar porque requiere de unidades experimentales más grandes que los factores de la subparcela, por ello se coloca en la parcela principal. Existen diseños en bloque incompletos, cuadrado de celosía, latino, los diseños alfa, etc., que ofrecen flexibilidad para variar el número de repeticiones, de tratamientos y el tamaño del bloque. Existen diseños en bloques incompletos y parcialmente balanceados.
El fracaso no es una opción en diseño de experimentos, ya que la investigación es muy costosa, tanto monetaria como emocionalmente, es por ello que esto es un mantra para los investigadores, es por ello que se debe convertir en un experto a la hora de diseñar un experimento y rodearse de personas aptas para lograrlo y nunca rendirse.
Bibliografia
Lawson, J. (2014). Design and analysis of experiments with r. CRC Press.
Lara Porras, A. (s.f) . Diseño de Experimentos. Universidad de Granada. Tomado de http://wpd.ugr.es/~bioestad/wp-content/uploads/Contenidos1.pdf
Práctica 7 | Estadística. (s. f.). Universidad de Granada. Recuperado de https://wpd.ugr.es/~bioestad/guia-de-r/practica-7/
Montgomery, D. C. (2005). Diseño y análisis de experimentos. Limusa Wiley.
CAPITULO IV Factoriales reducido. F.de Mendiburu. https://tarwi.lamolina.edu.pe/~fmendiburu/index-filer/academic/design/Factoriales.pdf
Diseños en cuadrados de Youden. http://wpd.ugr.es/~bioestad/wp-content/uploads/CuadradosYouden.pdf
Strip plot analysis using R. AGRON Stats Lectures. http://agroninfotech.blogspot.com/2018/06/strip-plot-analysis-using-r.html#example-data-set
LA ESTADÍSTICA: UNA ORQUESTA HECHA INSTRUMENTO. Test de Duncan. Jaume Llopis Pérez.https://jllopisperez.com/2013/01/28/test-de-duncan/