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('1_', 1:4))

data = data.frame(gen, procedencia, diam_geom)
head(data)
##   gen procedencia diam_geom
## 1 g_1         1_1  1.743952
## 2 g_1         1_2  1.776982
## 3 g_1         1_3  1.955871
## 4 g_1         1_4  1.807051
## 5 g_2         1_1  2.015515
## 6 g_2         1_2  2.205808
library(collapsibleTree)
## Warning: package 'collapsibleTree' was built under R version 4.2.2
collapsibleTreeSummary(data, 
                       c('procedencia', 'gen'), collapsed = FALSE)

Análisis descriptivo

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
ggplot(data)+
  aes(gen, diam_geom)+
  geom_point(size = 3,
             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

Adevertencia No interpretar el p-value respecto a los bloques, únicamente se interpretan los correspondientes a los factores

En este caso p-value de los genótipo > 5% (7.62%) se rechaza la hipótesis nula, por lo tanto podemos asumir estadísticamente que los genótipos no son iguales

Eficiencia de bloqueo ¿Valió la pena bloquear? ¿Preocuparse por la procedencia de las semillas? Se utiliza H = 1.334 (F-value de los bloques), cuando H > 1 (vale la pena bloquear) y cuando H < 1 (NO vale la pena bloquear)

Conclusión = Si valió la pena bloquear al observar que H > 1

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

# 2. Prueba de normalidades residuales
shapiro.test(res_mod1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod1
## W = 0.97765, p-value = 0.9725
# 3. Probando homocedasticidad (Varianza 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 autocorrealcionados y es un problema para el ánalisis de varianza.

Conclusión Valió la pena bloquear Estadisticamente no difieren los genotipos *Se cumplen los supuestos del ANOVA

#3. Diseño factorial simple en boques generealizados del 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(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