set.seed(123)
# Variable respuesta: CONDUCTANCIA ESTOMÁTICA

ce = c( rnorm(n = 15, mean = 100, sd = 10),
        rnorm(n = 15, mean = 120, sd = 12))

# Factor 1: HORA DE EVALUACIÓN
hora = gl(2, 15, 30, labels = c(6, 12))

# Factor 2: ILUMINACIÓN
ilum = gl(3, 5, 30, c('IB', 'IM', 'IA'))

# Bloqueo: PARCELA
parc = gl(5, 1, 30, paste0('B', 1:5))

datos = data.frame(hora, ilum, parc, ce)
head(datos)
##   hora ilum parc        ce
## 1    6   IB   B1  94.39524
## 2    6   IB   B2  97.69823
## 3    6   IB   B3 115.58708
## 4    6   IB   B4 100.70508
## 5    6   IB   B5 101.29288
## 6    6   IM   B1 117.15065
### ANÁLISIS DESCRIPTIVO

library(ggplot2)

ggplot(datos)+
  aes(y=ce)+
  geom_boxplot()

ggplot(datos)+
  aes(hora, ce)+
  geom_boxplot()

ggplot(datos)+
  aes(ilum, ce)+
  geom_boxplot()

ggplot(datos)+
  aes(parc, ce)+
  geom_boxplot()

ggplot(datos)+
  aes(ilum, ce, fill=hora)+
  geom_boxplot()

ggplot(datos)+
  aes(parc, ce, fill=hora)+
  geom_col(position = 'dodge') +
  facet_wrap(~ilum)

MODELO

\[y_{ijk}= \mu + \tau_i + \beta_j + (\tau\beta)_{ij}+\delta_k +\epsilon_{ijk}\]