Diseño desbalanceado FS-BG-CA

Modelo

\[y_{ijk} = \mu + \tau_i + \beta_{j} + tau\beta_{ij} + \epsilon_{ijk}\]

\(i=1,2, \dots,a:\text{num tratamientos}\) \(a\) numero de niveles del factor \(j=1,2,\dots,\text{num tratamientos}\) \(b\) numero de bloques \(k=1,2, \dots,r_i\) \(r_i\) repeticion de cada tratamiento

set.seed(123)

aceite = c(
  rnorm(12, 10, 0.8),
  rnorm(12, 11,0.78),
  rnorm(11, 9, 0.70),
  
  rnorm(12, 10, 0.8),
  rnorm(12, 11,0.78),
  rnorm(11, 9, 0.70)
)

bloque = gl(2,35,70, c('b1', 'b2'))

metodo = rep(rep(c('T1', 'T2', 'T3'), c(12,12,11)), 2)

datos = data.frame(metodo, bloque, aceite)
head(datos)

##   metodo bloque    aceite
## 1     T1     b1  9.551619
## 2     T1     b1  9.815858
## 3     T1     b1 11.246967
## 4     T1     b1 10.056407
## 5     T1     b1 10.103430
## 6     T1     b1 11.372052

table(datos$metodo, datos$bloque)

##     
##      b1 b2
##   T1 12 12
##   T2 12 12
##   T3 11 11

Corriendo como si fuera balanceado (no aplica para este caso)

mod1 = aov(aceite ~ bloque * metodo, datos)
summary (mod1)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## bloque         1   0.05   0.054   0.107    0.745    
## metodo         2  41.24  20.622  40.834 3.72e-12 ***
## bloque:metodo  2   0.65   0.323   0.640    0.531    
## Residuals     64  32.32   0.505                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Corriendo desbalanceado

ADVERTENCIA: aunque los dos den lo mismo en este caso, no siempre es asi, cuando hay más factores puede cambiar, entonces la siguiente manera es la que se debe utilizar

mod2 = anova (lm(aceite ~ bloque * metodo, datos))
mod2

## Analysis of Variance Table
## 
## Response: aceite
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## bloque         1  0.054  0.0541  0.1071    0.7446    
## metodo         2 41.243 20.6217 40.8342 3.716e-12 ***
## bloque:metodo  2  0.646  0.3231  0.6399    0.5307    
## Residuals     64 32.321  0.5050                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cuadrado Latino

lote <- c(rep("lote1",1), 
            rep("lote2",1), 
            rep("lote3",1), 
            rep("lote4",1), 
            rep("lote5",1)
            )
genotipo <- c(
  rep("gA",5), 
  rep("gB",5), 
  rep("gC",5), 
  rep("gD",5), 
  rep("gE",5)
  )
prov <- c("A","E","C","B","D",
          "C","B","A","D","E", 
          "B","C","D","E","A",
          "D","A","E","C","B",
          "E","D","B","A","C")
biom <- c(42,45,41,56,47, 47,
          54,46,52,49, 55,52,
          57,49,45, 51,44,47,
          50,54, 44,50,48,43,
          46)
 
data <- data.frame(
  lote, genotipo, prov, biom
  )

head(data)

##    lote genotipo prov biom
## 1 lote1       gA    A   42
## 2 lote2       gA    E   45
## 3 lote3       gA    C   41
## 4 lote4       gA    B   56
## 5 lote5       gA    D   47
## 6 lote1       gB    C   47

Graficos descriptivos

library(lattice)

bwplot(biom ~ genotipo | prov + lote,
       data)

Modelo

\[y = \mu + \tau_1 + \beta_j + \delta_k + \epsilon_{ijk}\]

\(i=1, \dots,p\) \(j=1, \dots,p\) \(k=1, \dots,p\)

tbl = matrix(data$prov, 5)
colnames(tbl) = unique(data$genotipo)
rownames(tbl) = unique(data$lote)
tbl

##       gA  gB  gC  gD  gE 
## lote1 "A" "C" "B" "D" "E"
## lote2 "E" "B" "C" "A" "D"
## lote3 "C" "A" "D" "E" "B"
## lote4 "B" "D" "E" "C" "A"
## lote5 "D" "E" "A" "B" "C"

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

mod <- lm(biom ~ lote + genotipo + prov,
          data)
anova(mod)

## Analysis of Variance Table
## 
## Response: biom
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## lote       4  17.76   4.440  0.7967 0.549839    
## genotipo   4 109.36  27.340  4.9055 0.014105 *  
## prov       4 286.16  71.540 12.8361 0.000271 ***
## Residuals 12  66.88   5.573                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Observaciones:

El valor del genotipo es 0.01 < 5% por lo tanto rechazo la hiótesis nula.
La eficiencia de bloqueo por lote fue de 0.7 < 1

bwplot(biom ~ genotipo | prov,
       data)

library(ggplot2)
ggplot(data)+
  aes(genotipo,
      biom,
      fill=prov)+
  geom_col(position = 'dodge')

Revision de supuestos

res_mod = mod$residuals
  
# 1. Normalidad
shapiro.test(res_mod)

## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod
## W = 0.97691, p-value = 0.8178

# se cumplio el supuesto de normalidad

#2. Igualdad de varianzas
bartlett.test(res_mod,
              genotipo)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  res_mod and genotipo
## Bartlett's K-squared = 5.9223, df = 4, p-value = 0.205

# se cumple el supuesto (varianzas iguales)

install.packages("TukeyC")

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)

library(TukeyC)

tt = TukeyC(mod, 'genotipo')
plot(tt)

DE_CLASE 28 ABRIL 2023

Valentina Melo

2023-04-28