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

#analisis 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_{ikj}\] #analisis inferencial

# Análisis de varianza

mod1 = aov(ce ~ parc + hora + ilum + hora:ilum, datos)
summary(mod1)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## parc         4  262.4    65.6   0.462 0.76262   
## hora         1 1805.8  1805.8  12.725 0.00193 **
## ilum         2  231.5   115.8   0.816 0.45651   
## hora:ilum    2   74.7    37.4   0.263 0.77120   
## Residuals   20 2838.4   141.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#bloque>1: no hay necesidad de hacer este bloqueo #hora: p_valor<5%: rechazo H_0, causa diferencias estadisticas en la respuesta “si difiere la CE por hora” #ilum: p_valor>5%: No rechazo H_0 hay evidencia estaditica de que no causa difenrencias en CE por iluminación

#CONCLUSIÓN: la hora es el unico factor inluyente en la variable respuesta.

hist(mod1$residuals)# parece ser que hay un comportamiento de "distribución normal en los residuales"

resid = mod1$residuals 
plot(resid) #pero tienen un buen comportamiento

#revision de supuestos

#NORMALIDAD
shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.96166, p-value = 0.3413
#p_valor>5% No se rechaza la H_0: los supuestos cumplen ua distribuución normal 
datos$trt = interaction(datos$hora,
                        datos$ilum)
ggplot(datos) +
  aes(trt, ce) +
  geom_boxplot()

#Existen unas cajas más grandes a otras = Preocupante - Se procura que las varianzas sean lo más iguales posibles.

#HOMOCEDASTICIDAD
bartlett.test(mod1$residuals,datos$trt)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod1$residuals and datos$trt
## Bartlett's K-squared = 8.4276, df = 5, p-value = 0.1342
#p_valor<5% Rechazo la H_0: las varianzas de los datos no son similares entre ellas
tapply(datos$ce, datos$trt, var)
##      6.IB     12.IB      6.IM     12.IM      6.IA     12.IA 
##  65.77564 286.70940 135.37080  17.09574  40.89132 229.33821

#Esos tratamientos no son comparables ya que la variablidad es muy diferente, entonces ¿Que hacer? - Debe ser 4 veces más grande y no casi 20 veces más grande

#Trabajo con el logaritmo de la conductancia, ya que encoge los datos = LOGARITMO DE LA CONDUCTANCIA ESTOMÁTICA

tapply(datos$ce, datos$trt, var)
##      6.IB     12.IB      6.IM     12.IM      6.IA     12.IA 
##  65.77564 286.70940 135.37080  17.09574  40.89132 229.33821
datos$ce_log = log(datos$ce)
ggplot(datos) +
  aes(trt, ce_log) +
  geom_boxplot()

ce_log_var = tapply(datos$ce_log, datos$trt, var)

sort(ce_log_var)/min(ce_log_var)
##     12.IM      6.IA      6.IB      6.IM     12.IA     12.IB 
##  1.000000  2.824160  4.346150  9.559637 12.433759 15.473325

#Estandarización de resultados (posible solución *en otros casos) = modificación y genera problemas ya que se trabaja con una variable diferente

modlog = aov(ce_log ~ parc + hora + ilum + hora:ilum, datos)
summary(modlog)
##             Df  Sum Sq Mean Sq F value  Pr(>F)   
## parc         4 0.02122 0.00530   0.467 0.75939   
## hora         1 0.14615 0.14615  12.862 0.00185 **
## ilum         2 0.01722 0.00861   0.758 0.48172   
## hora:ilum    2 0.00357 0.00178   0.157 0.85574   
## Residuals   20 0.22726 0.01136                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
reslog = modlog$residuals
shapiro.test(reslog)
## 
##  Shapiro-Wilk normality test
## 
## data:  reslog
## W = 0.96786, p-value = 0.4824
TT= TukeyHSD(modlog, 'hora');TT
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ce_log ~ parc + hora + ilum + hora:ilum, data = datos)
## 
## $hora
##           diff        lwr       upr     p adj
## 12-6 0.1395957 0.05840189 0.2207896 0.0018458
plot(TT)

#No se necesita hacer la prueba de Tukey ya que se tienen solo dos, si existieran más, se puede hacer.