DISEÑO CUADRADO GRECOLATINO

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.

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.

DISEÑO DE EXPERIMENTOS POR BLOQUES COMPLETOS AL AZAR

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)

DISEÑO CUADRADO LATINO INCOMPLETO (CUADRADO DE YOUDEN)

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.

  1. El análisis evaluando primero el efecto del tratamiento y después el de los bloques utilizando tres funciones
# 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.

  1. El análisis evaluando tanto los tratamientos como los bloques ejecutando solo una función.
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

PARCELAS DIVIDIDAS

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

BLOQUES O FRANJAS DIVIDIDAS

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.

BLOQUES INCOMPLETOS BALANCEADOS

` 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.

CUADRADO LATINO

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

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

8) Realice un resumen con la nota que aparece en las siguientes direcciones sobre:

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.

9) Seleccionar un artículo científico de una revista de agronomía donde se haya utilizado un diseño en parcelas divididas. Hacer las críticas constructivas sobre:

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.