set.seed(123)
# Variable: Pesos
peso = c(runif(40, 130, 150),
       runif(40, 150, 180),
       runif(40, 80, 130))
# Variable: Compuestos Fenolicos
cf = c(rnorm(40, 20, 1),
       rnorm(40, 12, 1.2),
       rnorm(40, 16, 1.15))

# Factor 1: Tiempo de Coccion
tiemp = gl(3, 40, 120, c(0, 10, 12))

# Factor 2: Vierdades
varie = gl(2, 20, 120, c('var1', 'var2'))
df = data.frame(varie, tiemp, cf, peso)
head(df)
##   varie tiemp       cf     peso
## 1  var1     0 20.37964 135.7516
## 2  var1     0 19.49768 145.7661
## 3  var1     0 19.66679 138.1795
## 4  var1     0 18.98142 147.6603
## 5  var1     0 18.92821 148.8093
## 6  var1     0 20.30353 130.9111
View(df)
boxplot(df$cf ~ df$varie)

boxplot(df$cf ~ df$tiemp)
lattice::bwplot(df$cf ~ df$tiemp | df$varie)

library(collapsibleTree)
collapsibleTreeSummary(df,
                       hierarchy = c('varie', 'tiemp'),
                       collapsed = FALSE)
mod1 = aov(cf ~ tiemp*varie, df)
summary(mod1)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## tiemp         2 1399.9   700.0 602.599 <2e-16 ***
## varie         1    5.4     5.4   4.639 0.0334 *  
## tiemp:varie   2    0.1     0.0   0.033 0.9678    
## Residuals   114  132.4     1.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tapply(df$cf, df$varie, mean)
##     var1     var2 
## 15.74983 16.17367
tapply(df$cf, df$tiemp, mean)
##        0       10       12 
## 20.12759 11.76136 15.99629
modt = aov(cf+peso ~ tiemp*varie, df)
summary(modt)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## tiemp         2  59372   29686 271.832 <2e-16 ***
## varie         1    271     271   2.483  0.118    
## tiemp:varie   2     67      34   0.307  0.736    
## Residuals   114  12450     109                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(mod1, 'tiemp')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = cf ~ tiemp * varie, data = df)
## 
## $tiemp
##            diff       lwr       upr p adj
## 10-0  -8.366227 -8.938526 -7.793928     0
## 12-0  -4.131296 -4.703594 -3.558997     0
## 12-10  4.234931  3.662632  4.807230     0
TukeyHSD(modt, 'tiemp')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = cf + peso ~ tiemp * varie, data = df)
## 
## $tiemp
##            diff        lwr       upr p adj
## 10-0   13.53163   7.982541  19.08072 2e-07
## 12-0  -38.94101 -44.490101 -33.39192 0e+00
## 12-10 -52.47264 -58.021731 -46.92355 0e+00
res1 = mod1$residuals
# Prueba de normalidad en los residuales
shapiro.test(res1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.98475, p-value = 0.1949
trat = interaction(df$tiemp, df$varie)
bartlett.test(res1, trat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  res1 and trat
## Bartlett's K-squared = 7.1214, df = 5, p-value = 0.2118
plot(res1)

plot(mod1$fitted.values, res1)

rest = modt$residuals
# Prueba de normalidad en los residuales
shapiro.test(rest)
## 
##  Shapiro-Wilk normality test
## 
## data:  rest
## W = 0.99105, p-value = 0.6311
bartlett.test(rest, trat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rest and trat
## Bartlett's K-squared = 32.322, df = 5, p-value = 5.129e-06
plot(rest)

plot(mod1$fitted.values, rest)

library(ggplot2)

mdf = aggregate(list(cf = df$cf),
                list(varie = df$varie,
                     tiemp = df$tiemp),
                mean)
ggplot(mdf)+
  aes(x = tiemp, y = cf, group = varie, color = varie)+
  # geom_boxplot()+
  geom_line()

ggplot(mdf)+
  aes(x = varie, y = cf, group = tiemp, color = tiemp)+
  # geom_boxplot()+
  geom_line(linetype = 'dashed')

mod1b = aov(cf ~ tiemp + varie, df)
summary(mod1b)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## tiemp         2 1399.9   700.0 612.819 <2e-16 ***
## varie         1    5.4     5.4   4.718 0.0319 *  
## Residuals   116  132.5     1.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
shapiro.test(mod1b$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1b$residuals
## W = 0.98569, p-value = 0.2366
bartlett.test(mod1b$residuals, trat)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod1b$residuals and trat
## Bartlett's K-squared = 7.1214, df = 5, p-value = 0.2118
df$bloq = gl(2, 10, 120, c('biol','alim'))
head(df)
##   varie tiemp       cf     peso bloq
## 1  var1     0 20.37964 135.7516 biol
## 2  var1     0 19.49768 145.7661 biol
## 3  var1     0 19.66679 138.1795 biol
## 4  var1     0 18.98142 147.6603 biol
## 5  var1     0 18.92821 148.8093 biol
## 6  var1     0 20.30353 130.9111 biol
collapsibleTreeSummary(df,
                       hierarchy = c('bloq','varie', 'tiemp'),
                       collapsed = FALSE)
mod2 = aov(cf ~ bloq * tiemp * varie, df)
summary(mod2)
##                   Df Sum Sq Mean Sq F value Pr(>F)    
## bloq               1    0.1     0.1   0.066 0.7981    
## tiemp              2 1399.9   700.0 589.774 <2e-16 ***
## varie              1    5.4     5.4   4.541 0.0354 *  
## bloq:tiemp         2    0.8     0.4   0.351 0.7044    
## bloq:varie         1    0.1     0.1   0.055 0.8148    
## tiemp:varie        2    0.1     0.0   0.032 0.9684    
## bloq:tiemp:varie   2    3.3     1.6   1.375 0.2573    
## Residuals        108  128.2     1.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod2b = aov(cf ~ bloq + tiemp * varie, df)
summary(mod2b)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## bloq          1    0.1     0.1   0.067 0.7968    
## tiemp         2 1399.9   700.0 597.666 <2e-16 ***
## varie         1    5.4     5.4   4.601 0.0341 *  
## tiemp:varie   2    0.1     0.0   0.033 0.9680    
## Residuals   113  132.3     1.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lattice::bwplot(cf ~ tiemp | varie + bloq, df)