DE_21.04

Diseño factorial simple en bloques al azar

set.seed(123)

# Respuesta
diam_geom = c(
  rnorm(4, 1.8, 0.1),
  rnorm(4, 2.0, 0.12),
  rnorm(4, 1.9, 0.09)
)
# Factor
gen = gl(3, 4, 12, paste0('g_', 1:3))
#Bloqueo
procedencia = gl(4, 1, 12, paste0('l', 1:4))

data = data.frame(gen, procedencia, diam_geom)
head(data)

##   gen procedencia diam_geom
## 1 g_1          l1  1.743952
## 2 g_1          l2  1.776982
## 3 g_1          l3  1.955871
## 4 g_1          l4  1.807051
## 5 g_2          l1  2.015515
## 6 g_2          l2  2.205808

library(collapsibleTree)
collapsibleTree::collapsibleTreeSummary(data,
                      c('procedencia',
                        'gen',
                        'diam_geom'),
                      collapsed = FALSE)

Análisis descriptivo

library(ggplot2)

ggplot(data)+
  aes(gen, diam_geom)+
  geom_point(size=7,
             color='yellow')+
  facet_wrap(~procedencia)+
  theme_dark()

Análisis inferencial

\[H_0: \mu_{g_1}=\mu_{g_2}=\mu_{g_3}\]

mod1 = aov(diam_geom ~ procedencia + gen,
          data)
summary(mod1)

##             Df  Sum Sq Mean Sq F value Pr(>F)  
## procedencia  3 0.04373 0.01458   1.334 0.3483  
## gen          2 0.08908 0.04454   4.078 0.0762 .
## Residuals    6 0.06554 0.01092                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Advertencia no interpretar el p-value respecto a los bloques, unicamente se interpretan los correspondientes a los factores

En este caso p-value de los genotipos > 5% (7.62%) se rechaza la hipotesis nula, por lo tanto podemos asumir estadisticamente que los genotipos no son iguales

Eficiencia de bloqueo corresponde a la pregunta ¿Valió la pena el bloqueo? (tener en cuenta la procedencia) H = 1.334 (F-value de los bloques), cuando H>1 sugiere que si valio la pena bloquear.

mod2 = aov(diam_geom ~ procedencia + gen,
          data)
summary(mod2)

##             Df  Sum Sq Mean Sq F value Pr(>F)  
## procedencia  3 0.04373 0.01458   1.334 0.3483  
## gen          2 0.08908 0.04454   4.078 0.0762 .
## Residuals    6 0.06554 0.01092                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 1. Extraer residuales del modelo (con bloques)
res_mod1 = mod1$residuals

# 2. Probando normalidad
shapiro.test

## function (x) 
## {
##     DNAME <- deparse1(substitute(x))
##     stopifnot(is.numeric(x))
##     x <- sort(x[complete.cases(x)])
##     n <- length(x)
##     if (is.na(n) || n < 3L || n > 5000L) 
##         stop("sample size must be between 3 and 5000")
##     rng <- x[n] - x[1L]
##     if (rng == 0) 
##         stop("all 'x' values are identical")
##     if (rng < 1e-10) 
##         x <- x/rng
##     res <- .Call(C_SWilk, x)
##     RVAL <- list(statistic = c(W = res[1]), p.value = res[2], 
##         method = "Shapiro-Wilk normality test", data.name = DNAME)
##     class(RVAL) <- "htest"
##     return(RVAL)
## }
## <bytecode: 0x00000262f84451c8>
## <environment: namespace:stats>

# 3. Probando homocedasticidad (varianzas iguales)
bartlett.test(res_mod1, data$gen)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  res_mod1 and data$gen
## Bartlett's K-squared = 1.709, df = 2, p-value = 0.4255

En ambos casos los p_value de las pruebas fueron >5% por lo cual se cumplen los supuestos.

plot(data$diam_geom,
     res_mod1,
     pch=16)

Si los residuales mostraran un patron se dice que estan autocorrelacionados y es un ploblema para el análisis de varianza

Conclusion Valio la pena bloquear Estadisticamente no difieren los genotipos *Se cumplen los supuestos del anova

Bloques generalizados al azar

set.seed(123)
# Respuesta
diam_geom = c(
  rnorm(20, 1.8, 0.1),
  rnorm(20, 2.0, 0.12),
  rnorm(20, 1.9, 0.09)
)
# Factor
gen = gl(3, 20, 60, paste0('g_', 1:3))
# Bloqueo
procedencia = gl(4, 5, 60, paste0('l_',1:4))

data = data.frame(gen, procedencia, diam_geom)
head(data)

##   gen procedencia diam_geom
## 1 g_1         l_1  1.743952
## 2 g_1         l_1  1.776982
## 3 g_1         l_1  1.955871
## 4 g_1         l_1  1.807051
## 5 g_1         l_1  1.812929
## 6 g_1         l_2  1.971506

library(collapsibleTree)
collapsibleTree::collapsibleTreeSummary(data,
                      c('procedencia',
                        'gen',
                        'diam_geom'),
                      collapsed = FALSE)

library(lattice)

bwplot(diam_geom ~ gen | procedencia, data)

mod3 = aov(diam_geom ~ procedencia * gen,
           data)

summary(mod3)

##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## procedencia      3 0.0310 0.01034   1.135    0.344    
## gen              2 0.3233 0.16164  17.732 1.71e-06 ***
## procedencia:gen  6 0.0407 0.00678   0.744    0.617    
## Residuals       48 0.4376 0.00912                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1