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)
library(collapsibleTree)
## Warning: package 'collapsibleTree' was built under R version 4.2.2
collapsibleTreeSummary(data,
                       c('procedencia',
                         'gen',
                         'diam_geom'),
                       collapsed = FALSE)

Diámetro geometrico.

\[Q_g = \sqrt[n_d]{d_1 * d_2 * d_3 ...}\]

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.2
ggplot(data)+
  aes(gen, diam_geom)+
  geom_point(size=7,
             color='yellow')+
  facet_wrap(~procedencia)+
  theme_dark()

Parece ser que el genotipo 2 es el que presente mejor compartamiento en las localidades 1, 2 y 3, en la 4, es el genotipo 3.

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

Aviso: No interpretar p valor respecto a los bloques, unicamente se interpretan los p-value correspondientes a los fatores, como el de gen.

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

Eficiencia de bloqueo: Corresponde a la pregunta ¿Valió la peba el bloqueo? (Tener en cuenta la procedencia) H = 1.334: Eficiencia de bloqueo (F-value de os bloques). Cuando H > 1, Si valió la pena bloquear.

Revisión de supuestos:

# 1. Extraer residuales del modelo (con bloques).

res_mod1 = mod1$residuals

# 2. Prueba de normalidad. 
shapiro.test(res_mod1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod1
## W = 0.97765, p-value = 0.9725
# 3. Prueba homocedaticidad (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)

Los residulaes que tienen patró se dice que están autocorrelacionado y representa un problema para el análisis de varianza.

Ejemplo: https://fvela.files.wordpress.com/2010/10/heterocedasticidad_notas.pdf

Conclusión:

3. Diseño factorial simple en 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 ('1_' , 1:4))

data = data.frame(gen, procedencia, diam_geom)
head(data)
library(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