Cuadrados de Youden en R.

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 se tenga 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 seria un diseño en cuadrado latino incompleto. Estos diseños fueron estudiados por W.J. Youden (1937) y se conocen con el nombre de cuadrados de Youden.

Usos.

design.youden(trt, r, serie = 2, seed = 0, kinds = “Super-Duper”,first=TRUE,randomization=TRUE)

Argumentos.

trt = Tratamientos.

r = Replicaciones o número de columna

serie = Gráfico de números, 1: 11,12; 2: 101,102; 3: 1001,1002

seed = Semilla

kinds = Metodo para aleatorizar

first TRUE or FALSE - Aleatoriazar la repetición 1

randomization = TRUE or FALSE - Aleatoriazar

Detalles.

kinds <- c(“Wichmann-Hill”, “Marsaglia-Multicarry”, “Super-Duper”, “Mersenne-Twister”, “KnuthTAOCP”, “user-supplied”, “Knuth-TAOCP-2002”, “default” )

Valor.

parameters = Parametros de diseño. sketch = Boceto de diseño. book = Libro de campo.

Modelo estadístico.

\[y_{ijk} = \mu + \tau_i + \beta_j + \gamma_k + \epsilon_{ijk} \\ i= 1,2...I \\ j= 1,2...J \\ k = 1,2...K \\ con \\I=J \\ K<1\]

Ejemplo de aleatorizacion.

library(agricolae)
variedad<-c("Pastusa Suprema","Rubí","Criolla","Parda Pastusa")
r <-3
outdesign <-design.youden(variedad,r,serie=4,seed=0)
print(outdesign$sketch)
##      [,1]              [,2]              [,3]             
## [1,] "Pastusa Suprema" "Criolla"         "Rubí"           
## [2,] "Criolla"         "Parda Pastusa"   "Pastusa Suprema"
## [3,] "Rubí"            "Pastusa Suprema" "Parda Pastusa"  
## [4,] "Parda Pastusa"   "Rubí"            "Criolla"
book <- outdesign$book
print(book) # field book.
##    plots row col        variedad
## 1  10001   1   1 Pastusa Suprema
## 2  10002   1   2         Criolla
## 3  10003   1   3            Rubí
## 4  20001   2   1         Criolla
## 5  20002   2   2   Parda Pastusa
## 6  20003   2   3 Pastusa Suprema
## 7  30001   3   1            Rubí
## 8  30002   3   2 Pastusa Suprema
## 9  30003   3   3   Parda Pastusa
## 10 40001   4   1   Parda Pastusa
## 11 40002   4   2            Rubí
## 12 40003   4   3         Criolla
print(matrix(as.numeric(book[,1]),byrow = TRUE, ncol = r))
##       [,1]  [,2]  [,3]
## [1,] 10001 10002 10003
## [2,] 20001 20002 20003
## [3,] 30001 30002 30003
## [4,] 40001 40002 40003
### Serpentine enumeration.
book <- zigzag(outdesign)
print(matrix(as.numeric(book[,1]),byrow = TRUE, ncol = r))
##       [,1]  [,2]  [,3]
## [1,] 10001 10002 10003
## [2,] 20003 20002 20001
## [3,] 30001 30002 30003
## [4,] 40003 40002 40001

Ejemplo

Se quiere evaluar la productividad de cuatro variedades de papa, 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 se cuenta con 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 papa (A, B, C, D). Los datos corresponden a la producción en kg/parcela.

knitr::include_graphics("/cloud/project/Captura.PNG")

Figura 1.Tabla de datos del Ejercicio

\[H_o = \mu_A = \mu_B = \mu_C = \mu_D \\ H_a = \mu_i \neq \mu_j \\ Algún \\i\neq j\]

library(readxl)
excel_data <- read_excel("excel.data.xlsx")
data <- data.frame(excel_data)
data
##    Productividad Nitrogeno Pendiente Variedad
## 1            756        N1        P1        D
## 2            720        N1        P2        A
## 3            689        N1        P3        C
## 4            596        N2        P1        A
## 5            855        N2        P2        B
## 6            780        N2        P3        D
## 7            750        N3        P1        C
## 8            975        N3        P2        D
## 9            899        N3        P3        B
## 10           950        N4        P1        B
## 11           870        N4        P2        C
## 12           880        N4        P3        A
data$Nitrogeno <- factor(data$Nitrogeno)
data$Nitrogeno
##  [1] N1 N1 N1 N2 N2 N2 N3 N3 N3 N4 N4 N4
## Levels: N1 N2 N3 N4
data$Pendiente <- factor(data$Pendiente) 
data$Pendiente
##  [1] P1 P2 P3 P1 P2 P3 P1 P2 P3 P1 P2 P3
## Levels: P1 P2 P3
data$Variedad <- factor(data$Variedad)
data$Variedad
##  [1] D A C A B D C D B B C A
## Levels: A B C D

El análisis de este modelo lo podemos realizar en R de dos formas: 1. Realizaremos el análisis evaluando primero el efecto del tratamiento y después el de los bloques utilizando tres funciones.

• Factor principal: Variedad

Para evaluar el efecto de los tratamientos, la suma de cuadrados de tratamientos debe ajustarse por bloques, por lo tanto primero se introducen los bloques y después los tratamientos.

mod1 <- aov(Productividad~  Pendiente + Nitrogeno + Variedad, data = data)
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 significación del 5%, deducimos que el factor principal: Variedades de papa es significativo/** Hay evidencia estadistica de que en la variedad hay interacción**

mod_1 <- aov(Productividad~  Pendiente * Nitrogeno + Variedad, data = data)
summary(mod_1)
##                     Df Sum Sq Mean Sq
## Pendiente            2  16952    8476
## Nitrogeno            3  73454   24485
## Variedad             3  47805   15935
## Pendiente:Nitrogeno  3   3013    1004

• Factor Bloque: Pendiente

Para evaluar el efecto del primero de los bloques, la suma de cuadrados de bloques debe ajustarse por los tratamientos, por lo tanto primero se introducen los tratamientos y después los bloques:

mod2 <- aov(Productividad~ Variedad +  Nitrogeno + Pendiente, data = data )
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 significación del 5%, deducimos que el Factor Bloque: Pendiente no es significativo / Hay evidencia estadistica de que en el factor pendiente no hay interacción.

• Factor Bloque: Nitrogeno.

Para evaluar el efecto del segundo bloque, la suma de cuadrados de bloques debe ajustarse también por los tratamientos, por lo tanto primero se introducen los tratamientos y después los bloques:

mod3 <- aov(Productividad~ Variedad +  Pendiente + Nitrogeno, data = data )
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 significación del 5%, deducimos que Factor Bloque: Nitrógeno es significativo / Hay evidencia estadistica de que en el factor nitrogeno no hay interacción

  1. El análisis evaluando tanto los tratamientos como los bloques ejecutando solo una función.
mod4 <- lm(Productividad~ Variedad +  Nitrogeno + Pendiente, data = data )
summary(mod4)
## 
## Call:
## lm(formula = Productividad ~ Variedad + Nitrogeno + Pendiente, 
##     data = data)
## 
## Residuals:
##       1       2       3       4       5       6       7       8       9      10 
##  -6.000  -6.375  12.375 -13.750  20.125  -6.375  13.750  12.375 -26.125   6.000 
##      11      12 
## -26.125  20.125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   634.38      27.44  23.115 0.000177 ***
## VariedadB     133.12      27.44   4.851 0.016722 *  
## VariedadC      -6.75      27.44  -0.246 0.821586    
## VariedadD     127.62      27.44   4.650 0.018753 *  
## NitrogenoN2   -24.62      27.44  -0.897 0.435695    
## NitrogenoN3   108.62      27.44   3.958 0.028792 *  
## NitrogenoN4   176.50      27.44   6.431 0.007621 ** 
## PendienteP2    92.00      22.41   4.106 0.026155 *  
## PendienteP3    49.00      22.41   2.187 0.116614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.69 on 3 degrees of freedom
## Multiple R-squared:  0.9787, Adjusted R-squared:  0.9218 
## F-statistic:  17.2 on 8 and 3 DF,  p-value: 0.01967

Los resultados obtenidos coinciden con los realizados primero a los bloques y después al tratamiento.

car::Anova(mod4, 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

Homocedasticidad

bartlett.test(data$Productividad, data$Variedad)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$Productividad and data$Variedad
## Bartlett's K-squared = 1.8021, df = 3, p-value = 0.6145
bartlett.test(data$Productividad, data$Pendiente)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$Productividad and data$Pendiente
## Bartlett's K-squared = 0.49855, df = 2, p-value = 0.7794
bartlett.test(data$Productividad, data$Nitrogeno)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$Productividad and data$Nitrogeno
## Bartlett's K-squared = 3.8699, df = 3, p-value = 0.2759

Los p-valores del factor tratamiento, Variedad de papa (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.

library(agricolae)
(duncan=duncan.test(mod4, "Variedad" , group = T))
## $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 papa 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 papa es la del tipo B, con una producción de 901.3333 Kg por parcela y la menor el tipo A, 732.0000 kg por parcela.

Diseño Carolina

Análisis estadístico de los diseños genéticos de Carolina I, II y III. Esta clase de diseños fueron elaborados por Comstock y Robinson (1948) en Carolina del Norte. Según la técnica de apareamiento entre progenitores se reconocen tres métodos:

Diseño I. Es también denominado diseño anidado o jerárquico. Bajo este esquema, cada macho es apareado con un grupo de hembras, con restricción que cada hembra sólo participa en una sola cruza.

Diseño II. Mencionado, como diseño factorial o cruzado. Esta técnica de apareamiento consiste en cruzar un grupo de progenitores machos con un conjunto de hembras en todas las combinaciones posibles. La única restricción es que unos progenitores actúan como machos y otros sólo como hembras.

Diseño III. Este fue desarrollado con la finalidad de estimar el grado de dominacia de los genes que controlan los caracteres en estudio.

Uso

carolina (modelo, datos)

Argumentos

modelo Constante marco de datos de datos

Detalles

modelo = 1,2 y 3 es I, II y III ver carolina1,2 y 3.

Valor

modelo de análisis del modelo (I, II o III) de diseño caroline y varianza y varianza aditiva de la interacción masculina, femenina y masculina.

Ejemplo

library(agricolae)
data(DC) ## DC, son datos para el análisis del diseño genético carolina I, II y III
carolina1 <- DC$carolina1
# str(carolina1)
output<-carolina(model=1,carolina1)
## Response(y):  yield 
## 
## Analysis of Variance Table
## 
## Response: y
##                             Df  Sum Sq Mean Sq F value    Pr(>F)    
## set                          1  0.5339  0.5339  7.2120 0.0099144 ** 
## set:replication              2  2.9894  1.4947 20.1914 4.335e-07 ***
## set:male                     4 22.1711  5.5428 74.8743 < 2.2e-16 ***
## set:male:female              6  4.8250  0.8042 10.8630 1.311e-07 ***
## set:replication:male:female 10  3.2072  0.3207  4.3325 0.0002462 ***
## Residuals                   48  3.5533  0.0740                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## CV: 8.286715     Mean: 3.283333
output[][-1]
## $var.m
## [1] 0.3948843
## 
## $var.f
## [1] 0.08057407
## 
## $var.A
## [1] 1.579537
## 
## $var.D
## [1] -1.257241
library(VCA)
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
da <-DC$carolina1
datos <- as.data.frame(da)
model <-anovaVCA(yield~set+rep:set+male:set+female:male:set+rep:female:male:set,datos)
model
## 
## 
## Result Variance Component Analysis:
## -----------------------------------
## 
##   Name                DF       SS        MS       VC       %Total    SD      
## 1 total               8.044979                    0.69694  100       0.834829
## 2 set                 1        0.533889  0.533889 0*       0*        0*      
## 3 set:rep             2        2.989444  1.494722 0.065222 9.358372  0.255386
## 4 set:male            4        22.171111 5.542778 0.394884 56.659736 0.628398
## 5 set:male:female     6        4.825     0.804167 0.080574 11.561124 0.283856
## 6 set:rep:male:female 10       3.207222  0.320722 0.082231 11.798936 0.28676 
## 7 error               48       3.553333  0.074028 0.074028 10.621832 0.27208 
##   CV[%]    
## 1 25.42627 
## 2 0*       
## 3 7.778267 
## 4 19.139031
## 5 8.645352 
## 6 8.733817 
## 7 8.286715 
## 
## Mean: 3.283333 (N = 72) 
## 
## Experimental Design: balanced  |  Method: ANOVA | * VC set to 0 | adapted MS used for total DF
#Varianza de machos *4
0.394884*4
## [1] 1.579536
#Varianza de Machos-Varianza de Hembras
(0.080574-0.394884)*4
## [1] -1.25724

Diseño Alpha

Estos diseños son generados por los arreglos alfa, y las tablas son rectangulares s por k (con s bloques y k <s columnas). El número de tratamientos debe ser igual a s * k y todas las unidades experimentales r * s * k (r repeticiones ).

Uso

design.alpha (trt, k, r, serie = 2, seed = 0, types = “Super-Duper”, randomization = TRUE)

Argumentos

trt = Tratamientos

k = bloque de tamaño

r = Replicaciones

Serie = Numero del grafico 1:11, 12; 2: 101,102; 3: 1001,1002

seed = semilla

Kind = método de para aleatorizar

randomization = TRUE OR FALSE- aleatorizar

library(agricolae)
trt<-1:30
t <- length(trt)
# size block k
k<-3
# Blocks s
s<-t/k
# replications r
r <- 2
outdesign<- design.alpha(trt,k,r,serie=2)
## 
## Alpha Design (0,1) - Serie  I 
## 
## Parameters Alpha Design
## =======================
## Treatmeans : 30
## Block size : 3
## Blocks     : 10
## Replication: 2 
## 
## Efficiency factor
## (E ) 0.6170213 
## 
## <<< Book >>>
book<-outdesign$book
plots<-book[,1]
dim(plots)<-c(k,s,r)
for (i in 1:r) print(t(plots[,,i]))
##       [,1] [,2] [,3]
##  [1,]  101  102  103
##  [2,]  104  105  106
##  [3,]  107  108  109
##  [4,]  110  111  112
##  [5,]  113  114  115
##  [6,]  116  117  118
##  [7,]  119  120  121
##  [8,]  122  123  124
##  [9,]  125  126  127
## [10,]  128  129  130
##       [,1] [,2] [,3]
##  [1,]  201  202  203
##  [2,]  204  205  206
##  [3,]  207  208  209
##  [4,]  210  211  212
##  [5,]  213  214  215
##  [6,]  216  217  218
##  [7,]  219  220  221
##  [8,]  222  223  224
##  [9,]  225  226  227
## [10,]  228  229  230
outdesign$sketch
## $rep1
##       [,1] [,2] [,3]
##  [1,] "16" "22" "15"
##  [2,] "25" "21" "9" 
##  [3,] "4"  "17" "10"
##  [4,] "24" "28" "19"
##  [5,] "26" "30" "14"
##  [6,] "29" "18" "23"
##  [7,] "11" "13" "12"
##  [8,] "27" "5"  "3" 
##  [9,] "7"  "2"  "6" 
## [10,] "20" "1"  "8" 
## 
## $rep2
##       [,1] [,2] [,3]
##  [1,] "16" "5"  "20"
##  [2,] "14" "24" "25"
##  [3,] "3"  "28" "8" 
##  [4,] "6"  "10" "12"
##  [5,] "19" "26" "1" 
##  [6,] "29" "30" "21"
##  [7,] "2"  "18" "9" 
##  [8,] "17" "11" "15"
##  [9,] "7"  "23" "4" 
## [10,] "27" "22" "13"

Diseño lattice

Diseños de lattice SIMPLE y TRIPLE. Aleatoriza tratamientos en lattice k x k.

Uso

design.lattice (trt, r = 3, serie = 2, seed = 0, types = “Super-Duper”, randomization = TRUE)

Argumentos

trt = tratamientos r — r = 2 (simple) o r = 3 (triple) serie = numero de grafico, 1:11, 12; 2: 101,102; 3: 1001,1002 seed = semilla Kind = método de para aleatorizar randomization = método para aleatorizar aleatorización = TRUE OR FALSE- aleatorizar

Ejemplo

library(agricolae)
# triple lattice
trt<-LETTERS[1:9]
outdesign<-design.lattice(trt,r=3,serie=2) # triple lattice design ( 9 trt)
## 
## Lattice design,  triple   3 x 3 
## 
## Efficiency factor
## (E ) 0.7272727 
## 
## <<< Book >>>
# simple lattice
trt<-1:100
outdesign<-design.lattice(trt,r=2,serie=3) # simple lattice design, 10x10
## 
## Lattice design,  simple   10 x 10 
## 
## Efficiency factor
## (E ) 0.8461538 
## 
## <<< Book >>>

Referencias

-Mendiburu, F. 2020. Package ‘agricolae’

-Mendiburu, F. 2020. agricolae tutorial (Version 1.3-3)

-Morales, I. 2001. COMPARACIÓN DE DOS DISEÑOS GENETICOS EN LA ESTIMACIÓN DE LOS COMPONENTES DE VARIANZA EN UNA POBLACIÓN DE MAÍZ.

-Roblero, E. 2017. Diseno Genetico Carolina del Norte.

-Universidad de Granada. 2020. Disponible en: https://wpd.ugr.es/~bioestad/guia-de-r/practica-7/#17