El modelo en cuadrado greco-latino se puede considerar como una extensión del cuadrado latino en el que se incluye una tercera variable de control o variable de bloque. En este modelo, como en el diseño en cuadrado latino, todos los factores deben tener el mismo número de niveles K y el número de observaciones necesarias sigue siendo K^2. Este diseño es, por tanto, una fracción del diseño completo en bloques aleatorizados con un factor principal y 3 factores secundarios que requeriría K^4 observaciones. Los cuadrados greco-latinos se obtienen por superposición de dos cuadrados latinos del mismo orden y ortogonales entre sí, uno de los cuadrados con letras latinas el otro con letras griegas. Dos cuadrados reciben el nombre de ortogonales si, al superponerlos, cada letra latina y griega aparecen juntas una sola vez en el cuadrado resultante.
EJERCICIO
Un experimentador estudia los efectos que tienen cinco dosis diferentes de la fertilización utilizada en el maíz en el rendimiento del cultivo. Cada dosis se hace con un lote de fertilizante diferente que sólo alcanza para probar cinco dosis. Además, las fertilizaciones son preparadas para varias parcelas, y puede haber diferencias sustanciales en las características y propiedades en cada una de ellas; por tanto, al parecer, hay dos factores perturbadores que serán “calculados en promedio” en dicho diseño: los lotes de fertilizante y las parcelas. El diseño apropiado para este problema consiste en probar cada formulación exactamente una vez en cada una de las parcelas.
Para esto, se parte del diseño de un Cuadrado Latino, sin embargo, se supone que existe un factor adicional, los montajes de prueba, que podría ser importante. Este factor considera 5 montajes de prueba diferentes denotados por letras griegas.
Las letras latinas representan a dosis
Las filas representan el factor principal lote de fertilizante.
Las columnas representan el factor parcelas.
Las letras griegas representan el factor montajes de prueba.
Tamaño del experimento número total de observaciones: 25
library(agricolae)
trt1 = factor(c("A","B","C","D","E")) #Dosis
trt2 =factor(c("b","c","d","e","f")) #Montajes
Tabla <- design.graeco( trt1, trt2, seed = 1, serie =2)
print(Tabla$sketch)
## [,1] [,2] [,3] [,4] [,5]
## [1,] "E e" "B f" "A b" "D c" "C d"
## [2,] "B b" "A c" "D d" "C e" "E f"
## [3,] "A d" "D e" "C f" "E b" "B c"
## [4,] "D f" "C b" "E c" "B d" "A e"
## [5,] "C c" "E d" "B e" "A f" "D b"
# Datos del diseño
Dosis<- c("A","B","C","D","E",
"B","C","D","E","A",
"C","D","E","A","B",
"D","E","A","B","C",
"E","A","B","C","D")
Lote <- c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5)
Parcelagre <- c(1,5,4,2,3,2,3,5,4,1,3,4,1,5,2,4,2,3,1,5,5,1,2,3,4)
Montajes <- c("b","c","d","e","f",
"d","e","f","b","c",
"f","b","c","d","e",
"c","d","e","f","b",
"e","f","b","c","d")
Rendimientogreco <- c(3.1, 6.1, 4.1, 3.4, 4, 2.9, 6.3, 2.7, 3.3, 5.4, 6.9, 4.6, 5.6,2.3,2.3,3.1, 5.6, 3, 1.1, 6.5,4.2, 2, 2.1, 6.8, 4.8)
# Factores
Dosis <- factor(Dosis); Dosis
## [1] A B C D E B C D E A C D E A B D E A B C E A B C D
## Levels: A B C D E
Lote <- factor(Lote); Lote
## [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5
## Levels: 1 2 3 4 5
Parcelagre <- factor(Parcelagre); Parcelagre
## [1] 1 5 4 2 3 2 3 5 4 1 3 4 1 5 2 4 2 3 1 5 5 1 2 3 4
## Levels: 1 2 3 4 5
Montajes <- factor(Montajes); Montajes
## [1] b c d e f d e f b c f b c d e c d e f b e f b c d
## Levels: b c d e f
dataGrecol = data.frame(Rendimientogreco,Dosis,Lote,Parcelagre,Montajes)
dataGrecol
## Rendimientogreco Dosis Lote Parcelagre Montajes
## 1 3.1 A 1 1 b
## 2 6.1 B 1 5 c
## 3 4.1 C 1 4 d
## 4 3.4 D 1 2 e
## 5 4.0 E 1 3 f
## 6 2.9 B 2 2 d
## 7 6.3 C 2 3 e
## 8 2.7 D 2 5 f
## 9 3.3 E 2 4 b
## 10 5.4 A 2 1 c
## 11 6.9 C 3 3 f
## 12 4.6 D 3 4 b
## 13 5.6 E 3 1 c
## 14 2.3 A 3 5 d
## 15 2.3 B 3 2 e
## 16 3.1 D 4 4 c
## 17 5.6 E 4 2 d
## 18 3.0 A 4 3 e
## 19 1.1 B 4 1 f
## 20 6.5 C 4 5 b
## 21 4.2 E 5 5 e
## 22 2.0 A 5 1 f
## 23 2.1 B 5 2 b
## 24 6.8 C 5 3 c
## 25 4.8 D 5 4 d
# Modelo lineal y tabla ANOVA
Grecol <- lm(Rendimientogreco ~ Dosis + Lote + Parcelagre + Montajes)
AnovaGrecol <- aov(Grecol)
summary(AnovaGrecol)
## Df Sum Sq Mean Sq F value Pr(>F)
## Dosis 4 33.71 8.427 4.058 0.0437 *
## Lote 4 0.65 0.164 0.079 0.9867
## Parcelagre 4 1.92 0.479 0.231 0.9136
## Montajes 4 15.42 3.854 1.856 0.2117
## Residuals 8 16.61 2.077
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La dosis es significativa ya que el p valor es <0.05.
library(agricolae)
Metodo.LSD <- LSD.test( y=AnovaGrecol, trt= "Dosis", group=TRUE, console = T)
##
## Study: AnovaGrecol ~ "Dosis"
##
## LSD t Test for Rendimientogreco
##
## Mean Square Error: 2.076633
##
## Dosis, means and individual ( 95 %) CI
##
## Rendimientogreco std r LCL UCL Min Max
## A 3.16 1.3352902 5 1.673876 4.646124 2.0 5.4
## B 2.90 1.9026298 5 1.413876 4.386124 1.1 6.1
## C 6.12 1.1541230 5 4.633876 7.606124 4.1 6.9
## D 3.72 0.9311283 5 2.233876 5.206124 2.7 4.8
## E 4.54 1.0237187 5 3.053876 6.026124 3.3 5.6
##
## Alpha: 0.05 ; DF Error: 8
## Critical Value of t: 2.306004
##
## least Significant Difference: 2.101696
##
## Treatments with the same letter are not significantly different.
##
## Rendimientogreco groups
## C 6.12 a
## E 4.54 ab
## D 3.72 b
## A 3.16 b
## B 2.90 b
bar.group(x = Metodo.LSD$groups, horiz = T, col="orange",
xlab="Rendimiento",
ylab="Dosis",
xlim=c(0,8),
main="Prueba de comparación de medias por medio del método LDS")
Hay diferencias significativas entre las dosis C y D, A, B, el mayor rendimiento se presenta con la dosis de fertilizante C y el menor rendimiento con la dosis de fertilizante B.
#Comprobación de supuestos
shapiro.test(AnovaGrecol$residuals)
##
## Shapiro-Wilk normality test
##
## data: AnovaGrecol$residuals
## W = 0.94557, p-value = 0.1988
qqnorm(AnovaGrecol$residuals)
qqline(AnovaGrecol$residuals)
De acuerdo al shapiro.test los residuos presentan un comportamiento normal.
Tabla_greco <- read.csv("C:/Users/julia/OneDrive/Escritorio/Dat_greco.csv", sep=";")
names(Tabla_greco)
## [1] "Parcelagre" "Lote_fert" "Dosis" "Montajes"
## [5] "Rendimientogreco"
str(Tabla_greco)
## 'data.frame': 25 obs. of 5 variables:
## $ Parcelagre : int 1 5 4 2 3 2 3 5 4 1 ...
## $ Lote_fert : int 1 1 1 1 1 2 2 2 2 2 ...
## $ Dosis : chr "A" "B" "C" "D" ...
## $ Montajes : chr "b" "c" "d" "e" ...
## $ Rendimientogreco: num 3.1 6.1 4.1 3.4 4 2.9 6.3 2.7 3.3 5.4 ...
Datgr=data.frame(Tabla_greco)
Datgr
## Parcelagre Lote_fert Dosis Montajes Rendimientogreco
## 1 1 1 A b 3.1
## 2 5 1 B c 6.1
## 3 4 1 C d 4.1
## 4 2 1 D e 3.4
## 5 3 1 E f 4.0
## 6 2 2 B d 2.9
## 7 3 2 C e 6.3
## 8 5 2 D f 2.7
## 9 4 2 E b 3.3
## 10 1 2 A c 5.4
## 11 3 3 C f 6.9
## 12 4 3 D b 4.6
## 13 1 3 E c 5.6
## 14 5 3 A d 2.3
## 15 2 3 B e 2.3
## 16 4 4 D c 3.1
## 17 2 4 E d 5.6
## 18 3 4 A e 3.0
## 19 1 4 B f 1.1
## 20 5 4 C b 6.5
## 21 5 5 E e 4.2
## 22 1 5 A f 2.0
## 23 2 5 B b 2.1
## 24 3 5 C c 6.8
## 25 4 5 D d 4.8
Datgr$Parcelagre <- factor(Datgr$Parcelagre); Datgr$Parcelagre
## [1] 1 5 4 2 3 2 3 5 4 1 3 4 1 5 2 4 2 3 1 5 5 1 2 3 4
## Levels: 1 2 3 4 5
Datgr$Lote_fert <- factor(Datgr$Lote_fert); Datgr$Lote_fert
## [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5
## Levels: 1 2 3 4 5
Datgr$Dosis <- factor(Datgr$Dosis); Datgr$Dosis
## [1] A B C D E B C D E A C D E A B D E A B C E A B C D
## Levels: A B C D E
Datgr$Montajes <- factor(Datgr$Montajes); Datgr$Montajes
## [1] b c d e f d e f b c f b c d e c d e f b e f b c d
## Levels: b c d e f
Datgr$Rendimientogreco <- factor(Datgr$Rendimientogreco); Datgr$Rendimientogreco
## [1] 3.1 6.1 4.1 3.4 4 2.9 6.3 2.7 3.3 5.4 6.9 4.6 5.6 2.3 2.3 3.1 5.6 3 1.1
## [20] 6.5 4.2 2 2.1 6.8 4.8
## 22 Levels: 1.1 2 2.1 2.3 2.7 2.9 3 3.1 3.3 3.4 4 4.1 4.2 4.6 4.8 5.4 ... 6.9
bartlett.test(AnovaGrecol$residuals~Tabla_greco$Dosis)
##
## Bartlett test of homogeneity of variances
##
## data: AnovaGrecol$residuals by Tabla_greco$Dosis
## Bartlett's K-squared = 0.41776, df = 4, p-value = 0.981
No se rechaza la hipotesis nula, dado que el p valor es 0.981, se observa que hay homogeneidad entre los bloques y los tratamientos.
plot(AnovaGrecol$residuals)
No se rechaza la hipótesis nula, la dispersión de los datos indica independencia.
El diseño en bloques completos al azar trata de comparar tres fuentes de variabilidad: el factor de tratamientos, el factor de bloques y el error aleatorio. El adjetivo completo se refiere a que en cada bloque se prueban todos los tratamientos. La aleatorización se hace dentro de cada bloque.
Es uno de los diseños más utilizados en experimentación agrícola. Este diseño es utilizado cuando el material experimental, campo agrícola, invernadero etc. Presentan una fuente de variabilidad conocida, factible de evaluar y de deducir el error experimental. Con ello se logra disminuir el error experimental, lo que incrementa la precisión en la comparación entre tratamientos. Recibe el nombre de bloque completo al azar, por que el material experimental se fracciona en bloques o en estrato uniformes dentro de si pero diferente entre si
EJERCICIO
Abeto blanco, Abeto del Pirineo, es un árbol de gran belleza por la elegancia de sus formas y el exquisito perfume balsámico que destilan sus hojas y cortezas. Destilando hojas y madera se obtiene aceite de trementina muy utilizado en medicina contra torceduras y contusiones. En estos últimos años se ha observado que la producción de semillas ha descendido y con objeto de conseguir buenas producciones se proponen tres tratamientos. Se observa que árboles diferentes tienen distintas características naturales de reproducción, este efecto de las diferencias entre los árboles 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 árbol un bloque completo. En este caso se trata de un diseño en bloques completos aleatorizados. El objetivo del estudio es comparar los tres tratamientos, por lo que se trata de un factor con tres niveles. Sin embargo, al realizar la medición sobre los distintos abetos, es posible que estos influyan sobre el número se semillas observadas. Por ello, y al no ser directamente motivo de estudio, los abetos son un factor secundario que recibe el nombre de bloque.
Semillas <- read.csv("C:/Users/julia/OneDrive/Escritorio/Data_abet.csv", sep=";")
library(agricolae)
# 3 tratamientos (niveles) y 10 bloques (abetos)
trt<-c("T1","T2","T3")
r <- 10
# Aleatorización con la función design.rcbd de la librería agricolae
outdesign <-design.rcbd(trt,r,serie=2, seed =986,kinds ="Wichmann-Hill", first=TRUE,
continue=FALSE,randomization=TRUE)
book <-outdesign$book; book
## plots block trt
## 1 101 1 T1
## 2 102 1 T3
## 3 103 1 T2
## 4 201 2 T3
## 5 202 2 T2
## 6 203 2 T1
## 7 301 3 T2
## 8 302 3 T1
## 9 303 3 T3
## 10 401 4 T3
## 11 402 4 T2
## 12 403 4 T1
## 13 501 5 T3
## 14 502 5 T2
## 15 503 5 T1
## 16 601 6 T1
## 17 602 6 T2
## 18 603 6 T3
## 19 701 7 T3
## 20 702 7 T2
## 21 703 7 T1
## 22 801 8 T3
## 23 802 8 T2
## 24 803 8 T1
## 25 901 9 T1
## 26 902 9 T3
## 27 903 9 T2
## 28 1001 10 T3
## 29 1002 10 T2
## 30 1003 10 T1
fieldbook <- zigzag(outdesign)
print(outdesign$sketch)
## [,1] [,2] [,3]
## [1,] "T1" "T3" "T2"
## [2,] "T3" "T2" "T1"
## [3,] "T2" "T1" "T3"
## [4,] "T3" "T2" "T1"
## [5,] "T3" "T2" "T1"
## [6,] "T1" "T2" "T3"
## [7,] "T3" "T2" "T1"
## [8,] "T3" "T2" "T1"
## [9,] "T1" "T3" "T2"
## [10,] "T3" "T2" "T1"
print(matrix(fieldbook[,1],byrow=TRUE,ncol=3))
## [,1] [,2] [,3]
## [1,] 101 102 103
## [2,] 203 202 201
## [3,] 301 302 303
## [4,] 403 402 401
## [5,] 501 502 503
## [6,] 603 602 601
## [7,] 701 702 703
## [8,] 803 802 801
## [9,] 901 902 903
## [10,] 1003 1002 1001
# Para numeración continua de los plots
outdesign <-design.rcbd(trt,r,serie=0,continue=TRUE)
head(outdesign$book)
## plots block trt
## 1 1 1 T1
## 2 2 1 T3
## 3 3 1 T2
## 4 4 2 T1
## 5 5 2 T3
## 6 6 2 T2
Ahora, con el fin de hacer los cálculos necesarios se transforman las columnas tratamiento y bloques (abeto) en factores
Semillas$Tratamiento = factor(Semillas$Tratamiento); 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 3
Semillas$Abeto = factor(Semillas$Abeto); 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 10
Se realiza el análisis de varianza
mod = aov(y ~ Tratamiento + Abeto, data = Semillas); mod
## Call:
## aov(formula = y ~ Tratamiento + Abeto, data = Semillas)
##
## Terms:
## Tratamiento Abeto Residuals
## Sum of Squares 16.2 54.8 15.8
## Deg. of Freedom 2 9 18
##
## Residual standard error: 0.936898
## Estimated effects may be unbalanced
summary(mod) # TABLA ANOVA
## 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 ' ' 1
A partir de los resultados del ANOVA es posible observar que las medias son estadísticamente significativas, dado que los p valores son menores a 0.05, se rechaza la hipótesis nula y concluye que no todas las medias son iguales.
Ahora para el RCBD se debe validar el modelo propuesto, esto consiste en estudiar si las hipótesis básicas del modelo están o no en contradicción con los datos observados. Es decir, si se satisfacen los supuestos del modelo: Normalidad, Independencia, Homocedasticidad. Para ello se utilizan procedimientos gráficos y analíticos.
En este modelo se ha supuesto otra hipotesis adicional: Aditividad de los efectos de tratamiento y bloque (no existe interaccion entre tratamiento y bloque). Por lo que hay que contrastar la hipotesis de aditividad de los efectos de tratamiento y bloque.
#Hipótesis de aditividad entre los bloques y los tratamientos
library(daewr)
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
Tuk <- Tukey1df(Semillas); Tuk
## 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.7202
## NULL
Puesto que el p-valor (Pr > F) es 0.04 no rechazamos la hipótesis nula de no interacción, es decir, no hay interacción entre los tratamientos aplicados y los abetos, no se muestra agrupación por letras.
#COMPROBACIÓN DE SUPUESTOS
shapiro.test(mod$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.96415, p-value = 0.3935
Como se observa se tiene un p valor de 0.3935 que aceptaría la hipótesis de normalidad por ser mayor al 5% (nivel de significación usual).
qqnorm(mod$residuals)
qqline(mod$residuals)
En el gráfico se observa que no existe ninguna desviación marcada de la normalidad.
# Igualdad de varianzas del factor principal
bartlett.test(Semillas$y, Semillas$Tratamiento)
##
## Bartlett test of homogeneity of variances
##
## data: Semillas$y and Semillas$Tratamiento
## Bartlett's K-squared = 4.1729, df = 2, p-value = 0.1241
# Igualdad de varianzas para el factor bloque
bartlett.test(Semillas$y, Semillas$Abeto)
##
## Bartlett test of homogeneity of variances
##
## data: Semillas$y and Semillas$Abeto
## Bartlett's K-squared = 4.0723, df = 9, p-value = 0.9066
En ambos casos, con base a la interpretación clásica de p valor al ser mayor que 0.05 no se rechaza la hipótesis de igualdad de varianzas ni en el factor principal ni en el factor bloque.
layout(matrix(c(1,2,3,4),2,2))
plot(mod)
En el primer gráfico que representa los residuos frente a los valores ajustados se observa que no hay ninguna tendencia evidente. Por lo tanto, no hay evidencia para rechazar la hipótesis de independencia.
#COMPARACIONES MÚLTIPLES
library(agricolae)
duncanT =duncan.test(mod, "Tratamiento" , group = T); duncanT #Para el factor principal
## $statistics
## MSerror Df Mean CV
## 0.8777778 18 9.2 10.18367
##
## $parameters
## test name.t ntr alpha
## Duncan Tratamiento 3 0.05
##
## $duncan
## Table CriticalRange
## 2 2.971152 0.8802727
## 3 3.117384 0.9235973
##
## $means
## y std r Min Max Q25 Q50 Q75
## 1 8.3 1.337494 10 7 11 7.25 8 8.75
## 2 9.2 1.135292 10 8 12 9.00 9 9.00
## 3 10.1 2.183270 10 7 14 9.25 10 11.50
##
## $comparison
## NULL
##
## $groups
## y groups
## 3 10.1 a
## 2 9.2 b
## 1 8.3 c
##
## attr(,"class")
## [1] "group"
En el apartado “$groups” se concluye que los tres tratamientos difieren significativamente entre sí. Se observa que la concentración media del número de semillas es mayor con el Tratamiento3 (10.1) y menor con el Tratamiento1 (8.3)
library(agricolae)
duncanA = duncan.test(mod, "Abeto" , group = T); duncanA #Para el factor bloque
## $statistics
## MSerror Df Mean CV
## 0.8777778 18 9.2 10.18367
##
## $parameters
## test name.t ntr alpha
## Duncan Abeto 10 0.05
##
## $duncan
## Table CriticalRange
## 2 2.971152 1.607151
## 3 3.117384 1.686250
## 4 3.209655 1.736161
## 5 3.273593 1.770746
## 6 3.320327 1.796026
## 7 3.355651 1.815133
## 8 3.382941 1.829895
## 9 3.404326 1.841462
## 10 3.421226 1.850604
##
## $means
## y std r Min Max Q25 Q50 Q75
## 1 8.666667 1.5275252 3 7 10 8.0 9 9.5
## 10 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 2 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 3 10.000000 1.7320508 3 9 12 9.0 9 10.5
## 4 10.333333 1.5275252 3 9 12 9.5 10 11.0
## 5 12.333333 1.5275252 3 11 14 11.5 12 13.0
## 6 9.000000 1.0000000 3 8 10 8.5 9 9.5
## 7 7.333333 0.5773503 3 7 8 7.0 7 7.5
## 8 7.666667 0.5773503 3 7 8 7.5 8 8.0
## 9 8.666667 1.5275252 3 7 10 8.0 9 9.5
##
## $comparison
## NULL
##
## $groups
## y groups
## 5 12.333333 a
## 4 10.333333 b
## 3 10.000000 b
## 10 9.000000 bc
## 2 9.000000 bc
## 6 9.000000 bc
## 1 8.666667 bc
## 9 8.666667 bc
## 8 7.666667 c
## 7 7.333333 c
##
## attr(,"class")
## [1] "group"
Se observa que la prueba de Duncan ha agrupado los abetos 7, 8, 1, 9, 2, 6 y 10 en un mismo grupo; los abetos 1, 9 ,2, 6, 10, 3 y 4, en otro grupo y un tercero esta formado únicamente por el Abeto5. Inmediatamente se ve que por ejemplo el Abeto 5 difiere de todos los demás, siendo en este abeto donde se produce el mayor número de semillas (12.333) y el menor en el Abeto 7 (7.333)
Hemos estudiado que en el diseño en cuadrado latino se tiene que verificar que los tres factores tengan el mismo número de niveles, es decir que haya el mismo número de filas, de columnas y de letras latinas. Sin embargo, 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. En general un cuadrado de Youden es un diseño balanceado por bloques incompletos, simétrico. En este diseño cada tratamiento ocurre una vez en cada columna; la posición del tratamiento dentro de un bloque indica el nivel del factor secundario correspondiente a las columnas; y el número de réplicas de un tratamiento dado es igual al número de tratamientos por bloque.
Modelo lineal estadístico \[y_{ijk}=μ+τ_i+β_j+α_k+ϵ_{ijki}\] \[i=1,2..I\] \[j=1,2..J\] \[k=1,2..K\] \[I = J;K < I\] \[y~sus~condiciones~laterales\]
Diseño de Youden con la librería agricolae Uso
design.youden(trt, r, serie = 2, seed = 0, kinds = “Super-Duper”,first=TRUE ,randomization=TRUE)
Argumentos
trt: Tratamientos
r: Replicaciones o número de columnas
serie: Número del gráfico 1: 11,12; 2: 101,102; 3: 1001,1002
seed: semilla
kinds: Metodo para aleatorizar
first: TRUE o FALSE - aleatorizar rep 1
randomization: TRUE or FALSE = aleatorizar
EJERCICIO
Consideremos el ejercicio del investigador que quiere evaluar la productividad de cuatro variedades de aguacate, A, B, C y D. Para ello, decide realizar el ensayo en un terreno que posee un gradiente de pendiente de oriente a occidente y además, diferencias en la disponibilidad de Nitrógeno de norte a sur. Se seleccionan cuatro disponibilidades de nitrógeno, pero sólo dispone de tres gradientes de pendiente. Para controlar estas posibles fuentes de variabilidad, el investigador decide utilizar un diseño en cuadrado de Youden con cuatro filas, las cuatro disponibilidades de Nitrógeno (NI, N2, N3, N4), tres columnas, los tres gradientes de pendientes (P1, P2, P3) y cuatro letras latinas, las variedades de aguacates (A, B, C, D). Los datos corresponden a la producción en kg/parcela • Variable respuesta: Productividad. • Factor: Variedad de aguacate. Es un factor de efectos fijos ya que desde el principio se establecen los niveles concretos que se van a analizar. • Bloques: Disponibilidad de Nitrógeno y Pendiente, con 4 y 3 niveles, respectivamente y ambos de efectos fijos. • Tamaño del experimento Número total de observaciones: 12.
Aguacate <- read.csv("C:/Users/julia/OneDrive/Escritorio/Productividad_a.csv", sep=";")
library(agricolae)
Var <-c("A","B","C","D")
r<-3
# Aleatorización con design.youden de la librería agricolae
outdesign <-design.youden(Var,r,serie=2,seed=23)
youden <- outdesign$book
print(outdesign$sketch)
## [,1] [,2] [,3]
## [1,] "C" "D" "A"
## [2,] "B" "C" "D"
## [3,] "A" "B" "C"
## [4,] "D" "A" "B"
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 403
print(youden)
## plots row col Var
## 1 101 1 1 C
## 2 102 1 2 D
## 3 103 1 3 A
## 4 201 2 1 B
## 5 202 2 2 C
## 6 203 2 3 D
## 7 301 3 1 A
## 8 302 3 2 B
## 9 303 3 3 C
## 10 401 4 1 D
## 11 402 4 2 A
## 12 403 4 3 B
Youden = data.frame(Aguacate)
View(Youden)
Transformar tanto la columna de los tratamiento como la de los bloques en un factor para poder realizar los cálculos posteriores adecuadamente
Youden$Productividad <- factor(Youden$Productividad); Youden$Productividad
## [1] 756 720 689 596 855 780 750 975 899 950 870 880
## Levels: 596 689 720 750 756 780 855 870 880 899 950 975
Youden$Nitrogeno <- factor(Youden$Nitrogeno); Youden$Nitrogeno
## [1] N1 N1 N1 N2 N2 N2 N3 N3 N3 N4 N4 N4
## Levels: N1 N2 N3 N4
Youden$Pendiente <- factor(Youden$Pendiente); Youden$Pendiente
## [1] P1 P2 P3 P1 P2 P3 P1 P2 P3 P1 P2 P3
## Levels: P1 P2 P3
Youden$Variedad <- factor(Youden$Variedad); Youden$Variedad
## [1] D A C A B D C D B B C A
## Levels: A B C D
La función “BIBsize(t , k)” de la librería daewr 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
library(colorspace)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(daewr)
library(AlgDesign)
BIBsize(t = 4 , k = 3)
## Posible BIB design with b= 4 and r= 3 lambda= 2
El análisis se puede realizar de 2 formas.
# Factor principal: variedad
#Primero se introducen los factores bloques y después los tratamientos.
mod1 <- aov(Productividad~ Pendiente + Nitrogeno + Variedad, data = Aguacate); mod1
## Call:
## aov(formula = Productividad ~ Pendiente + Nitrogeno + Variedad,
## data = Aguacate)
##
## Terms:
## Pendiente Nitrogeno Variedad Residuals
## Sum of Squares 16952.00 73454.00 47805.25 3012.75
## Deg. of Freedom 2 3 3 3
##
## Residual standard error: 31.6899
## Estimated effects may be unbalanced
summary(mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pendiente 2 16952 8476 8.44 0.0586 .
## Nitrogeno 3 73454 24485 24.38 0.0131 *
## Variedad 3 47805 15935 15.87 0.0241 *
## Residuals 3 3013 1004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El p-valor, 0.0241, es menor que el nivel de significancia del 5%, se deduce que el factor principal: Variedades del aguacate, es significativo
# Factor bloque: pendiente
mod2 <- aov(Productividad~ Variedad + Nitrogeno + Pendiente, data = Aguacate); mod2
## Call:
## aov(formula = Productividad ~ Variedad + Nitrogeno + Pendiente,
## data = Aguacate)
##
## Terms:
## Variedad Nitrogeno Pendiente Residuals
## Sum of Squares 50344.67 70914.58 16952.00 3012.75
## Deg. of Freedom 3 3 2 3
##
## Residual standard error: 31.6899
## Estimated effects may be unbalanced
summary(mod2)
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 3 50345 16782 16.71 0.0224 *
## Nitrogeno 3 70915 23638 23.54 0.0138 *
## Pendiente 2 16952 8476 8.44 0.0586 .
## Residuals 3 3013 1004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El p-valor, 0.0586, es mayor que el nivel de significancia del 5%, se deduce que el Factor Bloque: Pendiente, no es significativo
# Factor bloque: nitrogeno
mod3 <- aov(Productividad~ Variedad + Pendiente + Nitrogeno, data = Aguacate ); mod3
## Call:
## aov(formula = Productividad ~ Variedad + Pendiente + Nitrogeno,
## data = Aguacate)
##
## Terms:
## Variedad Pendiente Nitrogeno Residuals
## Sum of Squares 50344.67 16952.00 70914.58 3012.75
## Deg. of Freedom 3 2 3 3
##
## Residual standard error: 31.6899
## Estimated effects may be unbalanced
summary(mod3)
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 3 50345 16782 16.71 0.0224 *
## Pendiente 2 16952 8476 8.44 0.0586 .
## Nitrogeno 3 70915 23638 23.54 0.0138 *
## Residuals 3 3013 1004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El p-valor es 0.0138; menor que el nivel de significancia del 5%, deducimos que Factor Bloque: Nitrógeno, es significativo.
library(car)
## Loading required package: carData
modComp <- lm(Productividad~ Variedad + Nitrogeno + Pendiente, data = Aguacate ); modComp
##
## Call:
## lm(formula = Productividad ~ Variedad + Nitrogeno + Pendiente,
## data = Aguacate)
##
## Coefficients:
## (Intercept) VariedadB VariedadC VariedadD NitrogenoN2 NitrogenoN3
## 634.38 133.12 -6.75 127.62 -24.63 108.62
## NitrogenoN4 PendienteP2 PendienteP3
## 176.50 92.00 49.00
car::Anova(modComp, type="III")
## Anova Table (Type III tests)
##
## Response: Productividad
## Sum Sq Df F value Pr(>F)
## (Intercept) 536576 1 534.3047 0.0001774 ***
## Variedad 47805 3 15.8676 0.0240651 *
## Nitrogeno 70915 3 23.5382 0.0137944 *
## Pendiente 16952 2 8.4401 0.0586204 .
## Residuals 3013 3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Los resultados obtenidos coinciden con los realizados primero a los bloques y después al tratamiento.
#COMPROBACIÓN DE SUPUESTOS
shapiro.test(modComp$residuals)
##
## Shapiro-Wilk normality test
##
## data: modComp$residuals
## W = 0.90534, p-value = 0.1859
Como se observa se tiene un p-valor de 0.1859 que aceptaría la hipótesis de normalidad por ser mayor al 5% (nivel de significación usual)
qqnorm(modComp$residuals)
qqline(modComp$residuals)
En el gráfico se observa que no existe ninguna desviación marcada de la normalidad.Sin embargo, la cantidad de réplicas puede respresentar un problema para el análisis adecuado, ya que se tienen datos atípicos.
bartlett.test(Aguacate$Productividad, Aguacate$Variedad)
##
## Bartlett test of homogeneity of variances
##
## data: Aguacate$Productividad and Aguacate$Variedad
## Bartlett's K-squared = 1.8021, df = 3, p-value = 0.6145
bartlett.test(Aguacate$Productividad, Aguacate$Pendiente)
##
## Bartlett test of homogeneity of variances
##
## data: Aguacate$Productividad and Aguacate$Pendiente
## Bartlett's K-squared = 0.49855, df = 2, p-value = 0.7794
bartlett.test(Aguacate$Productividad, Aguacate$Nitrogeno)
##
## Bartlett test of homogeneity of variances
##
## data: Aguacate$Productividad and Aguacate$Nitrogeno
## Bartlett's K-squared = 3.8699, df = 3, p-value = 0.2759
Los p-valores del factor tratamiento, Variedad de aguacate (0.6145), del factor bloque Pendiente (0.779) y del factor bloque Nitrógeno (0.2759) son mayores que 0.05, por lo tanto no se puede rechazar la hipótesis de homogeneidad de la varianza del tratamiento y de los bloques.
layout(matrix(c(1,2,3,4),2,2))
plot(modComp)
## hat values (leverages) are all = 0.75
## and there are no factor predictors; no plot no. 5
En el primer gráfico, que representa los residuos frente a los valores ajustados se observan datos que se pueden considerar atípicos, pero de manera general los datos en su conjunto no presentan una tendencia. Por lo tanto, no hay evidencia para rechazar la hipótesis de independencia.
#COMPARACIONES MÚLTIPLES
library(agricolae)
dunyouden=duncan.test(modComp, "Variedad" , group = T); dunyouden
## $statistics
## MSerror Df Mean CV
## 1004.25 3 810 3.912334
##
## $parameters
## test name.t ntr alpha
## Duncan Variedad 4 0.05
##
## $duncan
## Table CriticalRange
## 2 4.500659 82.34484
## 3 4.515652 82.61915
## 4 4.472854 81.83611
##
## $means
## Productividad std r Min Max Q25 Q50 Q75
## A 732.0000 142.37977 3 596 880 658.0 720 800.0
## B 901.3333 47.54296 3 855 950 877.0 899 924.5
## C 769.6667 92.08873 3 689 870 719.5 750 810.0
## D 837.0000 120.11245 3 756 975 768.0 780 877.5
##
## $comparison
## NULL
##
## $groups
## Productividad groups
## B 901.3333 a
## D 837.0000 ab
## C 769.6667 bc
## A 732.0000 c
##
## attr(,"class")
## [1] "group"
En la tabla se muestran los subgrupos formados de medias iguales al utilizar el método de Duncan. Hay tres subconjuntos que se diferencian entre sí. Por una parte, el formado por la variedad de aguacate B y D, el subgrupo formado por D y C y el formado por A y C. También se observa que la mayor productividad de aguacate es la del tipo B, con una producción de 901.33 Kg por parcela y la menor el tipo A, 732.00 kg por parcela
Este es un diseño experimental combinado que resulta útil cuando al estudiar simultáneamentevarios factores, alguno o algunos de ellos deben ser aplicados sobre unidades experimentales relativamente grandes, pudiéndose aplicar el otro o los otros en unidades experimentales menores, dentro de las unidades mayores. El caso más sencillo es aquél en el que se tienen sólo dos factores, asignando los niveles de uno de ellos a las unidades mayores y los niveles del otro a las subunidades. A las unidades experimentales mayores suele llamárseles parcelas grandes o parcelas principales y a las unidades experimentales menores se le llamasub parcelas o subunidades. El factor correspondiente a las parcelas principales puede asignarse a éstas utilizando cualquiera de los esquemas de aleatorización básicos: Completamente al Azar, en Bloques al Azar o en Cuadro Latino. El factor correspondiente a las subparcelas se asigna al azar dentro de cada parcela principal; en tal sentido, las parcelas principales son análogas a bloques, solo que por asignarse a éstas los niveles de un efecto fijo y por existir repeticiones de las mismas,es posible evaluar tanto los efectos principales del factor asignado a las mismas como su posible interacción con el otro factor.
Ejemplo: En un ensayo de arroz, con el fin de analizar el efecto del riego sobre el rendimiento se analizaron 2 láminas de riego (L0 y L1) distintas. Cada lámina (parcelas principales) se repitió tres veces en orden aleatorio. Luego, se dividió cada parcela en cuatro subparcelas para dar cabida a 4 variedades de arroz (A, B, C, D), las que fueron asignadas al azar dentro de cada parcela.
En este caso se tiene que: - Variable respuesta: Rendimiento - Factor: variedad con 4 niveles. - Bloqueo: lamina de riego con dos niveles y en 3 repeticiones - Modelo completo: numero total de observaciones 24.
El modelo de parcelas divididas (DPV) es el siguiente
\[y_{ijk}=\mu+\alpha_i+\eta_{k(i)}+\beta_j+(\alpha\beta)_{ij}+\epsilon_{k(ij)}\] Usando la libreria agricolae se puede generar el modelo así
lamina1<- c("L0","L1")
lamina1
## [1] "L0" "L1"
variedad1<- c("A","B","C","D")
variedad1
## [1] "A" "B" "C" "D"
library(agricolae)
dis1<-design.split(lamina1, variedad1,r=3, design=c("rcbd","crd","lsd"),serie = 2, seed = 0, kinds = "Super-Duper", first=TRUE,randomization=TRUE); dis1
## $parameters
## $parameters$design
## [1] "split"
##
## $parameters[[2]]
## [1] TRUE
##
## $parameters$trt1
## [1] "L0" "L1"
##
## $parameters$applied
## [1] "rcbd"
##
## $parameters$r
## [1] 3
##
## $parameters$serie
## [1] 2
##
## $parameters$seed
## [1] 405950935
##
## $parameters$kinds
## [1] "Super-Duper"
##
##
## $book
## plots splots block lamina1 variedad1
## 1 101 1 1 L0 A
## 2 101 2 1 L0 D
## 3 101 3 1 L0 C
## 4 101 4 1 L0 B
## 5 102 1 1 L1 D
## 6 102 2 1 L1 A
## 7 102 3 1 L1 B
## 8 102 4 1 L1 C
## 9 103 1 2 L0 B
## 10 103 2 2 L0 C
## 11 103 3 2 L0 D
## 12 103 4 2 L0 A
## 13 104 1 2 L1 A
## 14 104 2 2 L1 D
## 15 104 3 2 L1 B
## 16 104 4 2 L1 C
## 17 105 1 3 L1 B
## 18 105 2 3 L1 C
## 19 105 3 3 L1 A
## 20 105 4 3 L1 D
## 21 106 1 3 L0 C
## 22 106 2 3 L0 B
## 23 106 3 3 L0 D
## 24 106 4 3 L0 A
DPV<- dis1$book
View(DPV)
Retomando los datos del ejemplo
library(readxl)
datos <- read_excel("C:/Users/julia/OneDrive/Escritorio/Parcelas divididas.xlsx")
library(lattice)
bwplot (datos$Rendimiento~datos$Variedad|datos$Lamina)
Pra realizar los calculos se transforman las columnas de los tratamientos en factores
rendimiento_1<-(datos$Rendimiento)
rendimiento_1
## [1] 266.3 259.3 340.7 236.6 629.5 544.9 519.9 409.3 296.6 335.7 252.8 358.8
## [13] 311.7 639.0 477.0 445.4 350.2 390.5 327.2 299.9 624.5 516.4 585.7 585.7
lamina<- factor(datos$Lamina)
lamina
## [1] L0 L0 L0 L0 L1 L1 L1 L1 L0 L0 L0 L0 L1 L1 L1 L1 L0 L0 L0 L0 L1 L1 L1 L1
## Levels: L0 L1
repeticion<- factor(datos$Repeticion)
repeticion
## [1] r1 r1 r1 r1 r1 r1 r1 r1 r2 r2 r2 r2 r2 r2 r2 r2 r3 r3 r3 r3 r3 r3 r3 r3
## Levels: r1 r2 r3
variedad_1<- factor(datos$Variedad)
variedad_1
## [1] A B C D D B C A C D A B A D C B B D C A C A B D
## Levels: A B C D
df= data.frame(repeticion, lamina, variedad_1, rendimiento_1)
Se hace el análisis de varianza, recordando el modelo
\[y_{ijk}=\mu+\alpha_i+\eta_{k(i)}+\beta_j+(\alpha\beta)_{ij}+\epsilon_{k(ij)}\]
modelo_1<- lm(aov(rendimiento_1~repeticion+lamina*variedad_1+Error(repeticion/lamina), data=datos))
modelo1<-sp.plot(repeticion, lamina, variedad_1, rendimiento_1) #esta es una función especial de la librería agricolae que permite calcular el analisis de varianza especícifamente para parcelas divididas
##
## ANALYSIS SPLIT PLOT: rendimiento_1
## Class level information
##
## lamina : L0 L1
## variedad_1 : A B C D
## repeticion : r1 r2 r3
##
## Number of observations: 24
##
## Analysis of Variance Table
##
## Response: rendimiento_1
## Df Sum Sq Mean Sq F value Pr(>F)
## repeticion 2 22891 11446 2.2836 0.304548
## lamina 1 276147 276147 55.0952 0.017671 *
## Ea 2 10024 5012
## variedad_1 3 51101 17034 6.3792 0.007858 **
## lamina:variedad_1 3 18931 6310 2.3632 0.122478
## Eb 12 32042 2670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 17 %, cv(b) = 12.4 %, Mean = 416.8167
Se puede observar en los p valor que para las repeticiones no hay diferencias significativas, mientras que tanto para la lámina como para las variedades, independientemente, las diferencias entre las medias si son significativas lo que se puede apreciar con el gráfico mostrado anteriormente como las cajas de las variedades están más dispersas. Por tanto, el riego y la variedad usada influyen en el rendimiento del cultivo. Por otro lado, en la interacción lamina:variedad no muestra una diferencia significativa por lo que el comportamiento en la respuesta de las variedades entre las laminas es semejante.
Para continuar hay que hacer la #Comprobación de supuestos
shapiro.test(modelo_1$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo_1$residuals
## W = 0.98741, p-value = 0.9862
Dado que el p valor > 0.05 se puede concluir que hay normalidad en los residuos
bartlett.test(rendimiento_1, repeticion)
##
## Bartlett test of homogeneity of variances
##
## data: rendimiento_1 and repeticion
## Bartlett's K-squared = 0.21403, df = 2, p-value = 0.8985
bartlett.test(rendimiento_1, lamina)
##
## Bartlett test of homogeneity of variances
##
## data: rendimiento_1 and lamina
## Bartlett's K-squared = 5.0587, df = 1, p-value = 0.0245
En este caso se acepta la homogeneidad de varianzas para las repeticiones pero se rechaza para el caso de las laminas.
layout(matrix(c(1,2,3,4),2,2))
plot(modelo_1)
#Comparaciones múltiples
Tukey_A<-HSD.test(modelo_1, "variedad_1", group = TRUE);Tukey_A
## $statistics
## MSerror Df Mean CV MSD
## 2670.196 12 416.8167 12.39728 88.57409
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey variedad_1 4 4.19866 0.05
##
## $means
## rendimiento_1 std r Min Max Q25 Q50 Q75
## A 342.7333 101.3105 6 252.8 516.4 274.700 305.80 384.900
## B 424.0500 124.9362 6 259.3 585.7 352.350 402.10 520.025
## C 430.9833 129.9641 6 296.6 624.5 330.575 408.85 509.175
## D 469.5000 171.0079 6 236.6 639.0 349.400 488.10 618.550
##
## $comparison
## NULL
##
## $groups
## rendimiento_1 groups
## D 469.5000 a
## C 430.9833 ab
## B 424.0500 ab
## A 342.7333 b
##
## attr(,"class")
## [1] "group"
dunc<- duncan.test(modelo_1, "variedad_1", group = TRUE)
dunc
## $statistics
## MSerror Df Mean CV
## 2670.196 12 416.8167 12.39728
##
## $parameters
## test name.t ntr alpha
## Duncan variedad_1 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.081307 65.00262
## 3 3.225244 68.03909
## 4 3.312453 69.87884
##
## $means
## rendimiento_1 std r Min Max Q25 Q50 Q75
## A 342.7333 101.3105 6 252.8 516.4 274.700 305.80 384.900
## B 424.0500 124.9362 6 259.3 585.7 352.350 402.10 520.025
## C 430.9833 129.9641 6 296.6 624.5 330.575 408.85 509.175
## D 469.5000 171.0079 6 236.6 639.0 349.400 488.10 618.550
##
## $comparison
## NULL
##
## $groups
## rendimiento_1 groups
## D 469.5000 a
## C 430.9833 a
## B 424.0500 a
## A 342.7333 b
##
## attr(,"class")
## [1] "group"
El resultado de ambas pruebas (Duncan y Tukey) señalan a la variedad D como la de mayor rendimiento, tienden a agrupar las variedades C y B en un mismo grupo con rendimientos similares. Estas dos ultimas variedades para efectos de su uso en campo pueden ser intercambiables. Finalmente, la variedad A se muestra como la de menor rendimiento y que la diferencia en el rendimiento entre A y las demás variedades es apreciable, por lo que no se recomendaria su uso bajo las condiciones riego en la lamina L1.
Por otro lado en la librería “deobioreserch” se puede encontrar una función similar a la “split.plot” proporcionadas por la librería agricolae. Aquí un ejemplo:
library(doebioresearch)
##
## Attaching package: 'doebioresearch'
## The following object is masked from 'package:lattice':
##
## stripplot
dis2<- splitplot(datos[4],repeticion,lamina, variedad_1,1); dis2
## $Rendimiento
## $Rendimiento[[1]]
## $Rendimiento[[1]][[1]]
## Analysis of Variance Table
##
## Response: dependent.var
## Df Sum Sq Mean Sq F value Pr(>F)
## block 2 22891 11446 2.2836 0.304548
## main.plot 1 276147 276147 55.0952 0.017671 *
## Ea 2 10024 5012
## sub.plot 3 51101 17034 6.3792 0.007858 **
## main.plot:sub.plot 3 18931 6310 2.3632 0.122478
## Eb 12 32042 2670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Rendimiento[[1]][[2]]
## [1] "CV(a): 16.985 , CV(b) : 12.397"
##
## $Rendimiento[[1]][[3]]
## [1] "R Square 0.922"
##
## $Rendimiento[[1]][[4]]
##
## Shapiro-Wilk normality test
##
## data: model$residuals
## W = 0.98741, p-value = 0.9862
##
##
## $Rendimiento[[1]][[5]]
## [1] "Normality assumption is not violated"
##
## $Rendimiento[[1]][[6]]
## [1] "The means of one or more levels of main plot factor are not same, so go for multiple comparison test"
##
## $Rendimiento[[1]][[7]]
## $Rendimiento[[1]][[7]][[1]]
## MSerror Df Mean CV t.value LSD
## 5012.183 2 416.8167 16.98511 4.302653 124.3581
##
## $Rendimiento[[1]][[7]][[2]]
## dependent.var groups
## L1 524.0833 a
## L0 309.5500 b
##
##
## $Rendimiento[[1]][[8]]
## [1] "The means of one or more levels of sub plot factor are not same, so go for multiple comparison test"
##
## $Rendimiento[[1]][[9]]
## $Rendimiento[[1]][[9]][[1]]
## MSerror Df Mean CV t.value LSD
## 2670.196 12 416.8167 12.39728 2.178813 65.00262
##
## $Rendimiento[[1]][[9]][[2]]
## dependent.var groups
## D 469.5000 a
## C 430.9833 a
## B 424.0500 a
## A 342.7333 b
##
##
## $Rendimiento[[1]][[10]]
## [1] "All the interaction level means are same so dont go for any multiple comparison test"
##
## $Rendimiento[[1]][[11]]
## $Rendimiento[[1]][[11]][[1]]
## MSerror Df Mean CV t.value LSD
## 2670.196 12 416.8167 12.39728 2.178813 91.92759
##
## $Rendimiento[[1]][[11]][[2]]
## dependent.var groups
## L1:D 618.0667 a
## L1:C 540.4667 ab
## L1:B 525.3333 b
## L1:A 412.4667 c
## L0:B 322.7667 cd
## L0:C 321.5000 cd
## L0:D 320.9333 cd
## L0:A 273.0000 d
Los diseños en bloques divididos, “strip-plot”, “split-block” o “criss-cross”, constituyen un caso particular de los diseños “split-plot”, que a su vez forman parte de los experimentos factoriales, en el que los factores que intervienen no se combinan aletoriamente entre sí, sino que están subordinados unos a otros. En este caso los niveles de un primer factor A se asignan a franjas de parcelas a lo largo de los bloques en una dirección, mientras que los niveles del segundo factor B se asignan a franjas de parcelas orientadas perpendicularmente a las del primero. Debido a la orientación perpendicular de los niveles de los dos factores, éstos suelen denominarse factor horizontal (A) y factor vertical (B). Para cada bloque se realiza una aleatorización independiente de los niveles de los dos factores.
El modelo de bloques divididos (strip-plot) es el siguiente
\[y_{ijk}=\mu+\gamma_{k}+\alpha_i+\epsilon_{ik}+\beta_j+\epsilon_{jk}+(\alpha\beta)_{ij}+\epsilon_{k(ij)}\]
Ejemplo: Supóngase que tiene un experimento con dos niveles de irrigación (alta y moderada) y cuatro variedades de caña (v1, v2, v3, v4) en cuatro bloques, con lo cual se busca analizar el rendimiento en Kg/parcela de la caña. En este caso se tiene que: - Variable respuesta: Rendimiento - Factor (columnas): variedad de caña con 4 niveles. - Factor (filas): irrigación con 2 niveles. - Bloqueo: se organizó en 4 bloques - Modelo completo:numero total de observaciones 24.
Usando la libreria agricolae se puede generar el modelo así:
library(readxl)
franjas<- read_excel("C:/Users/julia/OneDrive/Escritorio/Franjas divididas.xlsx")
irr<- c("Alta","Moderada")
var<- c("v1","v2","v3","v4")
dis2<- design.strip(irr, var,r=3, serie = 2, seed = 0, kinds = "Super-Duper",randomization=TRUE); dis2
## $parameters
## $parameters$design
## [1] "strip"
##
## $parameters$trt1
## [1] "Alta" "Moderada"
##
## $parameters$trt2
## [1] "v1" "v2" "v3" "v4"
##
## $parameters$r
## [1] 3
##
## $parameters$serie
## [1] 2
##
## $parameters$seed
## [1] -366863141
##
## $parameters$kinds
## [1] "Super-Duper"
##
##
## $book
## plots block irr var
## 1 101 1 Alta v2
## 2 102 1 Alta v4
## 3 103 1 Alta v1
## 4 104 1 Alta v3
## 5 105 1 Moderada v2
## 6 106 1 Moderada v4
## 7 107 1 Moderada v1
## 8 108 1 Moderada v3
## 9 201 2 Alta v4
## 10 202 2 Alta v2
## 11 203 2 Alta v1
## 12 204 2 Alta v3
## 13 205 2 Moderada v4
## 14 206 2 Moderada v2
## 15 207 2 Moderada v1
## 16 208 2 Moderada v3
## 17 301 3 Moderada v3
## 18 302 3 Moderada v1
## 19 303 3 Moderada v4
## 20 304 3 Moderada v2
## 21 305 3 Alta v3
## 22 306 3 Alta v1
## 23 307 3 Alta v4
## 24 308 3 Alta v2
View(dis2$book)
Ahora se hacen los calculos con los datos del ejemplo.
irrigacion<- factor(franjas$Irrigación)
variedad_2<- factor(franjas$Variedad)
bloques<- factor(franjas$Bloques)
rendimiento_2 <- franjas$Rendimiento
modelostrip<- strip.plot(bloques, variedad_2, irrigacion, rendimiento_2)
##
## ANALYSIS STRIP PLOT: rendimiento_2
## Class level information
##
## variedad_2 : v1 v2 v3 v4
## irrigacion : Alta Moderada
## bloques : 1 2 3
##
## Number of observations: 24
##
## model Y: rendimiento_2 ~ bloques + variedad_2 + Ea + irrigacion + Eb + irrigacion:variedad_2 + Ec
##
## Analysis of Variance Table
##
## Response: rendimiento_2
## Df Sum Sq Mean Sq F value Pr(>F)
## bloques 2 213.05 106.53 10.1934 0.011757 *
## variedad_2 3 137.69 45.90 9.2452 0.011454 *
## Ea 6 29.79 4.96 0.4750 0.806552
## irrigacion 1 492.32 492.32 120.0659 0.008226 **
## Eb 2 8.20 4.10 0.3924 0.691599
## irrigacion:variedad_2 3 10.86 3.62 0.3464 0.793455
## Ec 6 62.70 10.45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 1.8 %, cv(b) = 1.7 %, cv(c) = 2.7 %, Mean = 121.6708
bwplot (rendimiento_2~variedad_2|irrigacion)
Dado el p valor <0,05 se puede observa en la tabla que tento el bloqueo como la variedad y la irrigación tienen efectos significativos en el rendimiento de la caña, sin embargo esto no ocurre con la interacción entre variedad e irrigación, no hay diferencias significativas en la interacción. Dado que el p valor de la irrigación es menor se puede suponer que este factor es de mayor importancia en este experimento.
No se pudo encontrar cuales eran los supuestos de este modelo.
#Comparaciones multiples
gla<-modelostrip$gl.a
glb<-modelostrip$gl.b
glc<-modelostrip$gl.c
ea<- modelostrip$Ea
eb<- modelostrip$Eb
ec<- modelostrip$Ec
out<-LSD.test(rendimiento_2, variedad_2, gla, ea, alpha = 0.05,); out
## $statistics
## MSerror Df Mean CV t.value LSD
## 4.964306 6 121.6708 1.831229 2.446912 3.147654
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none variedad_2 4 0.05
##
## $means
## rendimiento_2 std r LCL UCL Min Max Q25 Q50 Q75
## v1 119.7000 5.683309 6 117.4743 121.9257 111.2 128.2 118.20 118.70 122.200
## v2 125.5667 8.040066 6 123.3409 127.7924 117.2 138.3 120.70 122.70 130.025
## v3 119.7000 5.576737 6 117.4743 121.9257 113.2 128.2 115.70 119.20 122.700
## v4 121.7167 5.944886 6 119.4909 123.9424 113.3 128.8 117.55 123.05 125.550
##
## $comparison
## NULL
##
## $groups
## rendimiento_2 groups
## v2 125.5667 a
## v4 121.7167 b
## v1 119.7000 b
## v3 119.7000 b
##
## attr(,"class")
## [1] "group"
out2<-LSD.test(rendimiento_2, irrigacion, glb, eb, alpha = 0.05,); out2
## $statistics
## MSerror Df Mean CV t.value LSD
## 4.100417 2 121.6708 1.664284 4.302653 3.556925
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none irrigacion 2 0.05
##
## $means
## rendimiento_2 std r LCL UCL Min Max Q25 Q50
## Alta 126.2000 5.423015 12 123.6849 128.7151 118.2 138.3 122.950 125.3
## Moderada 117.1417 3.552069 12 114.6265 119.6568 111.2 123.2 114.725 117.2
## Q75
## Alta 128.35
## Moderada 119.45
##
## $comparison
## NULL
##
## $groups
## rendimiento_2 groups
## Alta 126.2000 a
## Moderada 117.1417 b
##
## attr(,"class")
## [1] "group"
out3<-LSD.test(rendimiento_2, irrigacion:variedad_2, glc, ec, alpha = 0.05,); out3
## $statistics
## MSerror Df Mean CV t.value LSD
## 10.45042 6 121.6708 2.656931 2.446912 6.458617
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD none irrigacion:variedad_2 8 0.05
##
## $means
## rendimiento_2 std r LCL UCL Min Max Q25 Q50
## Alta:v1 123.2000 5.000000 3 118.6331 127.7669 118.2 128.2 120.70 123.2
## Alta:v2 130.9333 8.136543 3 126.3664 135.5003 122.2 138.3 127.25 132.3
## Alta:v3 124.2000 3.605551 3 119.6331 128.7669 121.2 128.2 122.20 123.2
## Alta:v4 126.4667 2.081666 3 121.8997 131.0336 124.8 128.8 125.30 125.8
## Moderada:v1 116.2000 4.358899 3 111.6331 120.7669 111.2 119.2 114.70 118.2
## Moderada:v2 120.2000 3.000000 3 115.6331 124.7669 117.2 123.2 118.70 120.2
## Moderada:v3 115.2000 2.000000 3 110.6331 119.7669 113.2 117.2 114.20 115.2
## Moderada:v4 116.9667 4.041452 3 112.3997 121.5336 113.3 121.3 114.80 116.3
## Q75
## Alta:v1 125.7
## Alta:v2 135.3
## Alta:v3 125.7
## Alta:v4 127.3
## Moderada:v1 118.7
## Moderada:v2 121.7
## Moderada:v3 116.2
## Moderada:v4 118.8
##
## $comparison
## NULL
##
## $groups
## rendimiento_2 groups
## Alta:v2 130.9333 a
## Alta:v4 126.4667 ab
## Alta:v3 124.2000 b
## Alta:v1 123.2000 bc
## Moderada:v2 120.2000 bcd
## Moderada:v4 116.9667 cd
## Moderada:v1 116.2000 d
## Moderada:v3 115.2000 d
##
## attr(,"class")
## [1] "group"
Dados los resultados de las pruebas LDS se puede concluir que la variedad v2 y la irrigación alta son los factores que mayor afectan el rendimiento de la caña de azucar. Al igual que lo visto en la tabla ANOVA, en este test LDS la interaccione entre ambos factores no es muy significativa, o sea, no hay diferencias tan apreciables en la interacción pero si se puede asegurar que la interacción variedad 2 con irrigación alta es la que produce mayor rendimiento.
` Los diseños en bloques incompletos balanceados (DBIB) fueron introducidos por Yates (1936). Lo que caracteriza este arreglo del material experimental es lo siguiente: - Cada bloque contiene k unidades experimentales. - Hay más tratamientos que unidades experimentales en un bloque. - Cada tratamiento aparece exactamente en r bloques. - Cada par de tratamientos ocurre justo en el mismo número de bloques lambda veces La propiedad de que cada par de tratamientos aparezca junto lambda veces hace posible que cualquier par de tratamientos sea comparable con el mismo error estándar. Además, el balanceamiento facilita el análisis estadístico, ya que los totales de tratamiento se ajustan en una sola operación para el conjunto de bloques donde aparece el tratamiento i (i = 1, 2, . . . , t). En este tipo de diseño, el análisis estadístico se centra en la información intrabloque, en la cual para estimar el efecto de los tratamientos se considera inicialmente la estimación de las parcelas dentro del mismo bloque. Así, los efectos de tratamientos deben tener un proceso de ajuste. Con los tratamientos ajustados se lleva a cabo la estimación de los efectos de tratamientos.
El modelo bloques incompletos balanceados (BIB) es el siguiente
\[y_{ijk}=\mu+\tau_i+\beta_j+\epsilon_{ij}\]
Como un ejemplo del uso de este modelo se tiene:
library(agricolae)
# 4 tratamientos y k=3 como tamaño del bloque
trt<-c("A","B","C","D")
k<-3
outdesign<-design.bib(trt,k,serie=2,seed =41,kinds ="Super-Duper") # seed = 41
##
## Parameters BIB
## ==============
## Lambda : 2
## treatmeans : 4
## Block size : 3
## Blocks : 4
## Replication: 3
##
## Efficiency factor 0.8888889
##
## <<< Book >>>
print(outdesign$parameters)
## $design
## [1] "bib"
##
## $trt
## [1] "A" "B" "C" "D"
##
## $k
## [1] 3
##
## $serie
## [1] 2
##
## $seed
## [1] 41
##
## $kinds
## [1] "Super-Duper"
book<-outdesign$book
plots <-as.numeric(book[,1])
matrix(plots,byrow=TRUE,ncol=k)
## [,1] [,2] [,3]
## [1,] 101 102 103
## [2,] 201 202 203
## [3,] 301 302 303
## [4,] 401 402 403
print(outdesign$sketch)
## [,1] [,2] [,3]
## [1,] "D" "C" "A"
## [2,] "A" "D" "B"
## [3,] "A" "B" "C"
## [4,] "C" "B" "D"
Ejercicio: Supóngase que un ingeniero químico cree que el tiempo de reacción en un proceso químico es función del catalizador empleado. De hecho, cuatro catalizadores están siendo investigados (c1, c2, c3, c4). El procedimiento experimental consiste en seleccionar un lote de materia prima (l1, l2, l3, l4), cargar una planta piloto, aplicar cada catalizador a ensayos separados en dicha planta y observar el tiempo de reacción.
En este caso se tiene que: - Variable respuesta: Tiempo de reacción - Factor: catalizador con 4 niveles. - Factor: lote con 4 niveles. - Bloqueo: numero de repeticiones 3 - Modelo completo:numero total de observaciones 12.
\[y_{ijk}=\mu+\tau_i+\beta_j+\epsilon_{ij}\]
library(readxl)
BIB <- read_excel("C:/Users/julia/OneDrive/Escritorio/BIB.xlsx")
catalizador<-factor(BIB$Catalizador)
lote_2<- factor(BIB$Lote)
tiempo<- BIB$Tiempo
test1<-BIB.test(lote_2, catalizador, tiempo, test= c("lsd","tukey","duncan","waller","snk"),
alpha = 0.05, group = TRUE,console=TRUE)
##
## ANALYSIS BIB: tiempo
## Class level information
##
## Block: l1 l2 l4 l3
## Trt : c1 c2 c3 c4
##
## Number of observations: 12
##
## Analysis of Variance Table
##
## Response: tiempo
## Df Sum Sq Mean Sq F value Pr(>F)
## block.unadj 3 55.00 18.3333 28.205 0.001468 **
## trt.adj 3 22.75 7.5833 11.667 0.010739 *
## Residuals 5 3.25 0.6500
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## coefficient of variation: 1.1 %
## tiempo Means: 72.5
##
## catalizador, statistics
##
## tiempo mean.adj SE r std Min Max
## c1 72.66667 71.375 0.4868051 3 1.527525 71 74
## c2 71.33333 71.625 0.4868051 3 4.041452 67 75
## c3 72.00000 72.000 0.4868051 3 3.605551 68 75
## c4 74.00000 75.000 0.4868051 3 1.732051 72 75
##
## LSD test
## Std.diff : 0.698212
## Alpha : 0.05
## LSD : 1.794811
## Parameters BIB
## Lambda : 2
## treatmeans : 4
## Block size : 3
## Blocks : 4
## Replication: 3
##
## Efficiency factor 0.8888889
##
## <<< Book >>>
##
## Comparison between treatments means
## Difference pvalue sig.
## c1 - c2 -0.250 0.7350
## c1 - c3 -0.625 0.4118
## c1 - c4 -3.625 0.0034 **
## c2 - c3 -0.375 0.6142
## c2 - c4 -3.375 0.0048 **
## c3 - c4 -3.000 0.0078 **
##
## Treatments with the same letter are not significantly different.
##
## tiempo groups
## c4 75.000 a
## c3 72.000 b
## c2 71.625 b
## c1 71.375 b
test1
## $parameters
## lambda treatmeans blockSize blocks r alpha test
## 2 4 3 4 3 0.05 BIB
##
## $statistics
## Mean Efficiency CV
## 72.5 0.8888889 1.112036
##
## $comparison
## NULL
##
## $means
## tiempo mean.adj SE r std Min Max Q25 Q50 Q75
## c1 72.66667 71.375 0.4868051 3 1.527525 71 74 72.0 73 73.5
## c2 71.33333 71.625 0.4868051 3 4.041452 67 75 69.5 72 73.5
## c3 72.00000 72.000 0.4868051 3 3.605551 68 75 70.5 73 74.0
## c4 74.00000 75.000 0.4868051 3 1.732051 72 75 73.5 75 75.0
##
## $groups
## tiempo groups
## c4 75.000 a
## c3 72.000 b
## c2 71.625 b
## c1 71.375 b
##
## attr(,"class")
## [1] "group"
library(lattice)
boxplot (tiempo~catalizador)
boxplot (tiempo~lote_2)
Según el análisis de varianza realizado por el test, los p valor tanto del lote como del catalizador fueron <0,05 por lo que hay diferencias significativas entre las medias de ambos factores. Como se puede apreciar en la grafica hay mayor diferencia entre los lotes y en menor medidida hay diferencia entre los catalizadores.
Luego de esto se tiene que hacer la #comprobación de los supuestos
test1aov<- aov(tiempo~lote_2+catalizador, data=BIB)
summary(test1aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## lote_2 3 55.00 18.333 28.20 0.00147 **
## catalizador 3 22.75 7.583 11.67 0.01074 *
## Residuals 5 3.25 0.650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(test1aov$residuals)
##
## Shapiro-Wilk normality test
##
## data: test1aov$residuals
## W = 0.96945, p-value = 0.905
qqnorm(test1aov$residuals)
qqline(test1aov$residuals)
Seegún el p valor >0,05 se puede decir que hay normalidad en la distribución de los residuos
bartlett.test(tiempo, catalizador)
##
## Bartlett test of homogeneity of variances
##
## data: tiempo and catalizador
## Bartlett's K-squared = 2.2078, df = 3, p-value = 0.5304
bartlett.test(tiempo, lote_2)
##
## Bartlett test of homogeneity of variances
##
## data: tiempo and lote_2
## Bartlett's K-squared = 3.4979, df = 3, p-value = 0.321
De igual manera hay homogeneidad en las varianzas
layout(matrix(c(1,2,3,4),2,2))
plot(test1aov)
#Comparaciones múltiples
Tukey_A<-HSD.test(test1aov, "catalizador", group = TRUE);Tukey_A
## $statistics
## MSerror Df Mean CV MSD
## 0.65 5 72.5 1.112036 2.428998
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey catalizador 4 5.218325 0.05
##
## $means
## tiempo std r Min Max Q25 Q50 Q75
## c1 72.66667 1.527525 3 71 74 72.0 73 73.5
## c2 71.33333 4.041452 3 67 75 69.5 72 73.5
## c3 72.00000 3.605551 3 68 75 70.5 73 74.0
## c4 74.00000 1.732051 3 72 75 73.5 75 75.0
##
## $comparison
## NULL
##
## $groups
## tiempo groups
## c4 74.00000 a
## c1 72.66667 ab
## c3 72.00000 ab
## c2 71.33333 b
##
## attr(,"class")
## [1] "group"
Tukey_B<-HSD.test(test1aov, "lote_2", group = TRUE);Tukey_B
## $statistics
## MSerror Df Mean CV MSD
## 0.65 5 72.5 1.112036 2.428998
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey lote_2 4 5.218325 0.05
##
## $means
## tiempo std r Min Max Q25 Q50 Q75
## l1 73.66667 1.1547005 3 73 75 73.0 73 74.0
## l2 74.66667 0.5773503 3 74 75 74.5 75 75.0
## l3 69.00000 2.6457513 3 67 72 67.5 68 70.0
## l4 72.66667 2.0816660 3 71 75 71.5 72 73.5
##
## $comparison
## NULL
##
## $groups
## tiempo groups
## l2 74.66667 a
## l1 73.66667 a
## l4 72.66667 a
## l3 69.00000 b
##
## attr(,"class")
## [1] "group"
duncanBIB1<- duncan.test(test1aov, "catalizador", group = TRUE)
duncanBIB1
## $statistics
## MSerror Df Mean CV
## 0.65 5 72.5 1.112036
##
## $parameters
## test name.t ntr alpha
## Duncan catalizador 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.635351 1.692164
## 3 3.748500 1.744832
## 4 3.796455 1.767154
##
## $means
## tiempo std r Min Max Q25 Q50 Q75
## c1 72.66667 1.527525 3 71 74 72.0 73 73.5
## c2 71.33333 4.041452 3 67 75 69.5 72 73.5
## c3 72.00000 3.605551 3 68 75 70.5 73 74.0
## c4 74.00000 1.732051 3 72 75 73.5 75 75.0
##
## $comparison
## NULL
##
## $groups
## tiempo groups
## c4 74.00000 a
## c1 72.66667 ab
## c3 72.00000 b
## c2 71.33333 b
##
## attr(,"class")
## [1] "group"
duncanBIB2<- duncan.test(test1aov, "lote_2", group = TRUE)
duncanBIB2
## $statistics
## MSerror Df Mean CV
## 0.65 5 72.5 1.112036
##
## $parameters
## test name.t ntr alpha
## Duncan lote_2 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.635351 1.692164
## 3 3.748500 1.744832
## 4 3.796455 1.767154
##
## $means
## tiempo std r Min Max Q25 Q50 Q75
## l1 73.66667 1.1547005 3 73 75 73.0 73 74.0
## l2 74.66667 0.5773503 3 74 75 74.5 75 75.0
## l3 69.00000 2.6457513 3 67 72 67.5 68 70.0
## l4 72.66667 2.0816660 3 71 75 71.5 72 73.5
##
## $comparison
## NULL
##
## $groups
## tiempo groups
## l2 74.66667 a
## l1 73.66667 ab
## l4 72.66667 b
## l3 69.00000 c
##
## attr(,"class")
## [1] "group"
En el ejercicio no se especifica el tipo de catalisis que se está analizando (si se trata de acelerar o retardar la reaccción) por lo que la interpretación variará. En cualquier caso segun las pruebas de Tukey y Duncan, si se quiere hablar sobre acelerar la reacción por ejemplo, se obserca como escoger el lote 3 y el catalizador 2 disminuye considerablemnte el tiempo de reacción siendo la combinación más eficiente para este propósito. En el caso de los otros lotes de material y catalizadores (L1, L2, L3 yC1, C3, C4) las diferencias no son tan grandes entre ellas por lo que no son eficientes.
El diseño en bloques aleatorios es adecuado cuando una fuente de variabilidad extraña se elimina (control local) para poder comparar un conjunto de medias muestrales asociadas con los tratamientos. Una característica importante de este tipo de diseño es su balance, que se logra asignando el mismo número de observaciones a cada tratamiento dentro de cada bloque. La misma clase de balance puede lograrse en otros tipos de diseño más complicados, en los cuales es conveniente eliminar el efecto de varias fuentes extrañas de variabilidad (dos o más). El diseño en cuadrado latino (DCL) se usa para eliminar dos fuentes de variabilidad, es decir, permite hacer la formación de bloques sistemática en dos direcciones (en el sentido de las filas y columnas). Por lo tanto, las filas y columnas representan en realidad dos restricciones sobre la aleatorización. De esta forma, se llama cuadro latino a un arreglo experimental obtenido a partir de una matriz cuadrada t×t en la que aparecen t elementos diferentes dados de tal forma que cada fila y cada columna contenga una sola vez cada uno de los elementos en consideración. Cada una de las t^2 celdas resultantes contiene una de las t letras que corresponde a los tratamientos, cada letra ocurre una y solo una vez en cada fila y columna
Ejemplo: En Kenett y Zacks (2000) se presenta un experimento, en el que se probaron cuatro métodos distintos, A, B, C y D, para preparar mezclas de cemento. Consistieron los métodos de dos relaciones de cemento y agua, y dos duraciones de mezclado. Los cuatro métodos fueron controlados por cuatro lotes durante cuatro días. El cemento se coló en cilindros y se midió la resistencia a la compresión en kg/cm2, a los 7 días de almacenamiento en cámaras especiales con 200C de temperatura y 50 % de humedad relativa.
En este caso se tiene que: - Variable respuesta: compresión/resistencia - Factor: metodo de preparación con 4 niveles. - Bloqueo: lote con 4 bloques. - Modelo completo:numero total de observaciones 16.
El modelo de cuadrado latino es el siguiente
\[y_{ijk}=\mu+\beta_i+\gamma_j+\tau_k+\epsilon_{ij}\]
Se generar un nuevo modelo de cuadrado latino a partir del ejemplo a partir de la libreria agricale
metodo_1<- c("A","B","C","D")
library(agricolae)
dis3<- design.lsd(metodo_1, serie = 2, seed = 0, kinds = "Super-Duper",first=TRUE,randomization=TRUE)
dis3
## $parameters
## $parameters$design
## [1] "lsd"
##
## $parameters$trt
## [1] "A" "B" "C" "D"
##
## $parameters$r
## [1] 4
##
## $parameters$serie
## [1] 2
##
## $parameters$seed
## [1] -710785715
##
## $parameters$kinds
## [1] "Super-Duper"
##
## $parameters[[7]]
## [1] TRUE
##
##
## $sketch
## [,1] [,2] [,3] [,4]
## [1,] "C" "B" "A" "D"
## [2,] "A" "D" "C" "B"
## [3,] "D" "C" "B" "A"
## [4,] "B" "A" "D" "C"
##
## $book
## plots row col metodo_1
## 1 101 1 1 C
## 2 102 1 2 B
## 3 103 1 3 A
## 4 104 1 4 D
## 5 201 2 1 A
## 6 202 2 2 D
## 7 203 2 3 C
## 8 204 2 4 B
## 9 301 3 1 D
## 10 302 3 2 C
## 11 303 3 3 B
## 12 304 3 4 A
## 13 401 4 1 B
## 14 402 4 2 A
## 15 403 4 3 D
## 16 404 4 4 C
Ahora se realizan los calculos transformanto las columnas en factores
library(readxl)
DCL<- read_excel("C:/Users/julia/OneDrive/Escritorio/Cuadrado latino.xlsx")
dia<- factor(DCL$Dia)
dia
## [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
## Levels: 1 2 3 4
lote_1<- factor(DCL$Lote)
lote_1
## [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
## Levels: 1 2 3 4
metodo<- factor(DCL$Metodo)
metodo
## [1] A B C D B A D C C D A B D C B A
## Levels: A B C D
compresion<- DCL$Compresion
compresion
## [1] 303 299 290 290 280 321 313 282 275 315 319 300 304 293 295 305
Recordando el modelo: \[y_{ijk}=\mu+\beta_i+\gamma_j+\tau_k+\epsilon_{ij}\]
Se realiza el analisis de varianza
modelo_3<- lm(compresion~dia+lote_1+metodo, data=DCL)
modelo_3
##
## Call:
## lm(formula = compresion ~ dia + lote_1 + metodo, data = DCL)
##
## Coefficients:
## (Intercept) dia2 dia3 dia4 lote_12 lote_13
## 300.00 3.50 6.75 3.75 16.50 13.75
## lote_14 metodoB metodoC metodoD
## 3.75 -18.50 -27.00 -6.50
aovmodelo<- aov(modelo_3)
summary(aovmodelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## dia 3 91.5 30.5 0.685 0.59290
## lote_1 3 745.5 248.5 5.584 0.03591 *
## metodo 3 1750.0 583.3 13.109 0.00482 **
## Residuals 6 267.0 44.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(lattice)
boxplot (compresion~metodo)
boxplot (compresion~lote_1)
Por los p valor <0,05 obtenidos se puede observar que el metodo de preparación de cemento es el que tiene diferencias significativas entre las medias y en menor medida estas diferencias significativas también las presenta entre los lotes de control. En la gráfica se puede observar como varian los datos (y las medias) de los métodos de preparación y de igual manera los lotes.
Luego de esto se debe hacer la
#Comprobación de supuestos
shapiro.test(aovmodelo$residuals)
##
## Shapiro-Wilk normality test
##
## data: aovmodelo$residuals
## W = 0.97849, p-value = 0.9505
Por el p valor >0,05 se puede concluir que normalidad en los residuos.
Homogeneidad de varianzas
bartlett.test(compresion, metodo)
##
## Bartlett test of homogeneity of variances
##
## data: compresion and metodo
## Bartlett's K-squared = 0.3161, df = 3, p-value = 0.957
bartlett.test(compresion, lote_1)
##
## Bartlett test of homogeneity of variances
##
## data: compresion and lote_1
## Bartlett's K-squared = 0.40779, df = 3, p-value = 0.9386
Por el p valor >0,05 tanto para el lote como para el método hay homogeneidad en las varianzas.
layout(matrix(c(1,2,3,4),2,2))
plot(modelo_3)
#Comparaciones múltiples
Tukey_A<-HSD.test(aovmodelo, "metodo", group = TRUE);Tukey_A
## $statistics
## MSerror Df Mean CV MSD
## 44.5 6 299 2.231048 16.32886
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey metodo 4 4.895599 0.05
##
## $means
## compresion std r Min Max Q25 Q50 Q75
## A 312.0 9.309493 4 303 321 304.50 312.0 319.50
## B 293.5 9.255629 4 280 300 291.25 297.0 299.25
## C 285.0 8.124038 4 275 293 280.25 286.0 290.75
## D 305.5 11.387127 4 290 315 300.50 308.5 313.50
##
## $comparison
## NULL
##
## $groups
## compresion groups
## A 312.0 a
## D 305.5 ab
## B 293.5 bc
## C 285.0 c
##
## attr(,"class")
## [1] "group"
duncanlatino<- duncan.test(aovmodelo, "metodo", group = TRUE)
duncanlatino
## $statistics
## MSerror Df Mean CV
## 44.5 6 299 2.231048
##
## $parameters
## test name.t ntr alpha
## Duncan metodo 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.460456 11.54206
## 3 3.586498 11.96246
## 4 3.648934 12.17071
##
## $means
## compresion std r Min Max Q25 Q50 Q75
## A 312.0 9.309493 4 303 321 304.50 312.0 319.50
## B 293.5 9.255629 4 280 300 291.25 297.0 299.25
## C 285.0 8.124038 4 275 293 280.25 286.0 290.75
## D 305.5 11.387127 4 290 315 300.50 308.5 313.50
##
## $comparison
## NULL
##
## $groups
## compresion groups
## A 312.0 a
## D 305.5 a
## B 293.5 b
## C 285.0 b
##
## attr(,"class")
## [1] "group"
Tukey_B<-HSD.test(aovmodelo, "lote_1", group = TRUE);Tukey_B
## $statistics
## MSerror Df Mean CV MSD
## 44.5 6 299 2.231048 16.32886
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey lote_1 4 4.895599 0.05
##
## $means
## compresion std r Min Max Q25 Q50 Q75
## 1 290.50 15.15476 4 275 304 278.75 291.5 303.25
## 2 307.00 13.16561 4 293 321 297.50 307.0 316.50
## 3 304.25 13.93736 4 290 319 293.75 304.0 314.50
## 4 294.25 10.27538 4 282 305 288.00 295.0 301.25
##
## $comparison
## NULL
##
## $groups
## compresion groups
## 2 307.00 a
## 3 304.25 ab
## 4 294.25 ab
## 1 290.50 b
##
## attr(,"class")
## [1] "group"
duncanlatino2<- duncan.test(aovmodelo, "lote_1", group = TRUE)
duncanlatino2
## $statistics
## MSerror Df Mean CV
## 44.5 6 299 2.231048
##
## $parameters
## test name.t ntr alpha
## Duncan lote_1 4 0.05
##
## $duncan
## Table CriticalRange
## 2 3.460456 11.54206
## 3 3.586498 11.96246
## 4 3.648934 12.17071
##
## $means
## compresion std r Min Max Q25 Q50 Q75
## 1 290.50 15.15476 4 275 304 278.75 291.5 303.25
## 2 307.00 13.16561 4 293 321 297.50 307.0 316.50
## 3 304.25 13.93736 4 290 319 293.75 304.0 314.50
## 4 294.25 10.27538 4 282 305 288.00 295.0 301.25
##
## $comparison
## NULL
##
## $groups
## compresion groups
## 2 307.00 a
## 3 304.25 ab
## 4 294.25 bc
## 1 290.50 c
##
## attr(,"class")
## [1] "group"
Finalmente las pruebas de Tukey y Duncan agrupan a los metodos en dos principales grupos: a para metodos A y D, b para B y C, y comparando con la gráfica de métodos se puede concluir que el mejor método de preparación de cemento es el A seguido por el D mientras que B y C son menos eficientes, especialmente este último. Por otro lado para el lote de control ocurre de manera similar en la que se crean 2 grupos y sub grupos, en donde los lotes 2 y 3 son los que generan cemento (a) con mayor resistencia y los lotes 4 y 1 los de menos resistencia (b). Por tanto se puede concluir que el cemento con mayor resistencia es obtenido del método A a través de los lotes 2 o 3.
REFERENCIAS - Domínguez, A., Fernandez, R. & Trapero, A. (2010). Experimentación en agricultura. Andalucía. Consejería de Agricultura y Pesca. Andalucia, España. En: https://www.juntadeandalucia.es/export/drupaljda/1337160941EXPERIMENTACION.pdf
Gabriel J, Castro C, Valverde A, Indacochea B (2017). Diseños experimentales: Teoría y práctica para experimentos agropecuarios. Grupo COMPAS, Universidad Estatal del Sur de Manabí (UNESUM), Jipijapa, Ecuador. En: http://142.93.18.15:8080/jspui/bitstream/123456789/116/1/Dise%C3%B1os%20Experimentales%20-%20Enero%2030%2C%202017%20-%20VERSI%C3%93N%20FINAL.compressed.pdf
López, A., Melo, S. & Melo, O. (2020). Diseño de experimentos, métodos y aplicaciones. Coordinación de publicaciones-Facultad de Ciencias, Universidad Nacional de Colombia. Bogotá D.C, Colombia. En:
Mendiburu, F. (2020). Statistical Procedures for Agricultural Research. En: https://cran.r-project.org/web/packages/agricolae/agricolae.pdf
Universidad de Granada. (SF). Diseños en cuadrados greco-latinos. En: http://wpd.ugr.es/~bioestad/wp-content/uploads/GrecoLatinos1.pdf
Universidad de Granada. (SF).Diseños en cuadrados de Youden. En: http://wpd.ugr.es/~bioestad/wp-content/uploads/CuadradosYouden.pdf
Universidad de Granada. (SF). Práctica 7. Diseño estadístico de experimentos. En: https://wpd.ugr.es/~bioestad/guia-de-r/practica-7/#17
Vera, J. F. & Vera, J. A.(2018) Resumen de principios de diseños experimentales , Editorial Grupo Compás, Guayaquil Ecuador. En: https://repositorio.uteq.edu.ec/handle/43000/3764
Yepes, V. (2014). Diseño de experimentos por bloques completos al azar. Universitat Politècnica de València. En: https://victoryepes.blogs.upv.es/2014/06/30/diseno-de-experimentos-por-bloques-completos-al-azar/
Universidad Nacional de Colombia/Facultad de Ciencias Agrarias/Dpto. de Agronomía Diseño de experimentos/Parcial I/ 35%/ Apellidos y Nombres: Carlos Julián Garavito Arias & Leydy Viviana González Clavijo Cédula: 1000617957 & 1003519403 respectivamente
OPCIÓN 2 Parcial Diseño de experimentos
En ocasiones al realizar un experimento existen niveles de factores que resultan difíciles de cambiar en contraposición a otros que sí poseen dicha facilidad. En experimentos de parcelas divididas los factores difíciles de cambiar son variados menos y se les llaman factores de parcelas completas, mientras que a los fáciles de cambiar se les llama factores de subparcelas. Cuando en el experimento no se desea reiniciar todos los factores entre cada ejecución la mejor opción es escoger un diseño de parcelas divididas, que no significa omitir la aleatorización sino que la forma en que se aleatoriza es más estructurada.
Para el correcto análisis de este diseño hay que reconocer que al no reiniciar todo el factor parcela las observaciones obtenidas serán más relacionadas o similares entre ellas que cuando se reinician todos los factores. En parcelas divididas hay dos aleatorizaciones separadas: determinar el orden en que se corren las parcelas completas y las observaciones colectadas dentro de cada parcela completa. Además, otra razón para justificar el uso de este tipo de diseño es la reducción de costos del experimento al hacer cambios sobre los factores de las subparcelas en vez de las parcelas completas. .
La unidad experimental es conocida también como la unidad de replicación y se refiere a la entidad más pequeña que se asigna de manera independiente a todas las demás unidades a manera particular, está definida en términos de asignaciones de tratamiento independientes. Si a manera de ejemplo se tiene una situación en la que se quiere evaluar el efecto del agua de un arroyo contaminado en las lesiones de los peces y se instalan dos acuarios (control y contaminada) con 50 peces cada una y después de 30 días se capturan 10 peces de cada acuario para contar el número de heridas, las unidades experimentales no son los peces, para que fuera de esta manera se tendría que aplicar el tratamiento a cada pez de manera individual, lo que es poco práctico, aquí la unidad experimental es entonces el acuario con los 50 peces, donde se aplicó el tratamiento.
La unidad de observación o unidad de muestreo se define como la entidad física sobre la que se mide un resultado de interés en un experimento (en términos de medidas de resultado). Estas unidades de observación suelen estar contenidas en unidades experimentales. Con base en el ejemplo de los peces, los 10 sujetos sobre los que se midieron los efectos constituyen la unidad de muestreo
El diseño de experimentos tiene cuatro pilares principales: replicación, aleatorización, bloqueo y tamaño de las unidades experimentales. Para las investigaciones biológicas se sigue el precepto “el fracaso no es una opción” debido a los costos y tiempo de inversión para la realización de los experimentos. Una guía para diseñar experimentos exitosos contempla varios factores, entre ellos la forma y la escala de replicación, esto proporciona un mecanismo para estimar el error experimental, aumentar la precisión de un experimento, aumentar el alcance de la inferencia del experimento y afecta el control de error, puede ocurrir en múltiples niveles o escalas pero debe aplicarse al nivel de la unidad experimental.
Es posible que los tratamientos experimentales se repliquen a una escala insuficiente para su adecuado análisis, en cuyo caso se estarían confundiendo los tratamientos con unidades experimentales, a esto se le conoce como pseudorreplicación y tiene solución a partir de la replicación a la escala adecuada en el tiempo o en el espacio. Los pasos a seguir para la correcta formulación son: definir explícitamente la unidad experimental que forma el primer nivel de replicación contemplando la escala adecuada, se pueden diseñar niveles adicionales.
Deben considerarse el número de repeticiones y la distribución de réplicas entre las diversas formas de réplica, ya que el número de repeticiones repercute de manera directa, predecible, repetible y tangible sobre la precisión y la capacidad de detectar diferencias entre tratamientos, por lo tanto, con un manejo adecuado del experimento se evidencian diferencias significativas entre las medias de tratamiento con el nivel suficiente de replicación.
Existen situaciones bajo las cuales se marca una tendencia a invertir todos los recursos experimentales hacia múltiples tratamientos y no se contempla la observación independiente de estos ni tampoco la replicación, esto constituye una problemática, por la ineficacia de este método, aunque constituyen una opción viable para experimentos de tamaño reducido en granjas. Por su parte, los diseños aumentados representan una forma específica de diseño que puede manejar cientos o miles de tratamientos, la mayoría de los cuales no se replican. Respecto a la aleatorización como principio es aplicable a la correcta realización de experimentos en dos niveles. Primero, se debe definir la población y elegir una muestra aleatoria o representativa en el caso de que lo deseado sea representar una población más grande que la muestra, mientras que, cuando se desea una inferencia sólo para un pequeño número de niveles, cada uno de los cuales puede incluirse en el experimento, la opción suele ser tratar esto como un efecto fijo. El segundo aspecto de la aleatorización hace referencia a la asignación de tratamientos a las unidades experimentales con la misma probabilidad de aplicación. La aplicación más simple de aleatorización se da mediante el diseño completamente aleatorizado (CRD) que consta de un bloque.
El bloqueo, por su parte, se usa con el fin de tener precisión para crear grupos de unidades experimentales más homogéneos de lo que serían con un muestreo aleatorio de toda la población de unidades experimentales o por conveniencia, permitiendo diferentes tamaños de unidades experimentales cuando se requieran parcelas o áreas experimentales más grandes para la aplicación de un factor en comparación con otros. Un ejemplo es el diseño de bloques completos al azar (RCBD) como de bloque simple o el diseño de parcelas divididas como una restricción específica de aleatorización que se coloca en ciertos experimentos factoriales, con tratamientos de dos o más factores, entre otros diseños.
El tamaño de las unidades experimentales se rige por la ley de Smith que demuestra la relación empírica entre la varianza por unidad y el tamaño de la parcela en experimentos agronómicos de campo. Por lo general, los diseños más simples son los menos eficientes, es por eso que lo recomendable es buscar continuamente nuevos diseños que permitan diversos análisis y salir de lo conocido.
Artículo: De Mesquita, J., De Lima, A., Andrade, F., Da Silva, T., Ferreira, L., & De Oliveira, F. et al. (2020). Chlorophyll a fluorescence and development of zucchini plants under nitrogen and silicon fertilization. Agronomía Colombiana, 38(1), 45-52. https://revistas.unal.edu.co/index.php/agrocol/article/view/79172/74906
Respecto a la estructura factorial los autores lo plantean como: los tratamientos se distribuyeron en un esquema de parcelas divididas en un diseño de bloques completos al azar. Realmente la cantidad de pruebas/repeticiones que se usaron son demasiado pocas, contrario a lo que en literatura se recomienda para el uso de parcelas divididas, o sea que se trate de experimentos considerablemente amplios.
Los autores no especifican las razones, lo que representa un inconveniente para la comprensión del diseño experimental, sin embargo, es posible que la parcela principal se designa en función del nivel de silicio debido a que el efecto de este se evalúa en 2 niveles, por lo que es más sencillo el manejo experimental, por su parte, la dosis de nitrógeno se comprende en las subparcelas debido a que se cuenta con más niveles (cinco) que en el primer caso, por lo tanto, es más factible que sean estas las que se traten de esta forma. El diseño en parcelas divididas de la manera mencionada puede deberse, además, a una cuestión de costos.
En el análisis de los datos no se menciona la revisión para ninguno de los supuestos necesarios al realizar el análisis de varianza, los cuales son pertinentes para llevar a cabo el respectivo análisis bajo los parámetros establecidos, lo ideal es hacer mención a estos y sus respectivas comprobaciones.
Resume el análisis de varianza por los valores cuadrados medios, para clorofila a (Cla), clorofila b (Clb), clorofila total (tCl), relación clorofila a/b (Cla/Clb), fluorescencia inicial (F0) , máxima fluorescencia (Fm), fluorescencia variable (Fv) y rendimiento cuántico del fotosistema II (Fv/Fm) de calabacín (Cucurbita Pepo L.) bajo dosis de nitrógeno y fertilización foliar de silicio en Catolé do Rocha. Con todo lo anterior, proporciona un panorama completo para la interpretación de los resultados teniendo en cuenta cada variable; la interpretación de la tabla se puede hacer con facilidad, ya que se explica claramente cada componente. En esta tabla se evidencia que existen diferencias significativas de las medias para nitrógeno sobre el contenido de clorofila y fluorescencia, además, para la interacción nitrógeno - silicio para todas las variables, mientras que con la fertilización de silicio existe diferencia significativa para las variables de fluorescencia, a excepción de la variable Fv.
Los datos se sometieron a un análisis de varianza mediante la prueba F al 5% de probabilidad. Se evalúa la varianza de clorofila y fluorescencia de manera independiente, los resultados se presentan resumidos en una tabla, de acuerdo a los valores de cuadrado medio. La tabla de análisis no permite conocer si la clorofila se relaciona de alguna forma con la fluorescencia al presentar los datos de manera separada, no es posible apreciar directamente si existe un efecto de la fertilización y el nivel de silicio similar para ambas variables respuesta.
En la tabla 6 las medias seguidas de letras iguales en la misma columna y línea no difieren estadísticamente entre sí según la prueba de Tukey al 5% de probabilidad, este método permite probar las diferencias entre medias de tratamientos, la agrupación de acuerdo a las letras proporciona una descripción clara para la interpretación de los resultados y la división entre niveles de factor y de bloqueo facilitan la interpretación.
En la interacción entre la fertilización con nitrógeno y silicio sobre el crecimiento, el índice de clorofila y la fluorescencia de la clorofila de plantas de calabacín, los factores se mencionan de manera explícita, los autores hacen referencia directa a estos en la tabla de análisis de varianza, donde se aprecia que la interacción entre los factores afectó significativamente a todas las variables analizadas.
La parcela fue organizada en un único bloque dividido por el nivel de silicio (0 y 6 g/planta) con subparcelas por 5 niveles de nitrógeno (30, 60, 90, 120 y 150 kg/ha) con tres repeticiones y una planta por repetición. Por lo tanto, el factor de bloqueo es la presencia o no de silicio. En el artículo se hace referencia a un solo factor de bloqueo en dos niveles respecto al silicio de manera explícita.
El diseño experimental es balanceado, ya que presenta todos los datos necesarios para el desarrollo de los análisis pertinentes, las muestras de cada tratamiento se muestran el mismo número de veces.
La unidad experimental está claramente definida por la parcela dividida en un único bloque de acuerdo al nivel de silicio, con subparcelas en 5 niveles de nitrógeno.
Se utilizó el programa estadístico R, no se especifica la librería, esto es un problema en caso de que se quieran replicar y comprobar los análisis de los datos experimentales, por lo tanto, se dificultan las características de replicable y reproducible.
Muy pocas réplicas son utilizadas (tres repeticiones, una planta por cada una), por tratamiento, para que el experimento pueda ser concluyente sería ideal aumentar el número de repeticiones, posiblemente conservar las 30 parcelas experimentales pero disminuir el número de niveles para las dosis de nitrógeno.