Diseño Factorial Completo en arreglo Completamente al Azar (DFCCA)

set.seed(123)
# Respuesta: Volumen de raices
vol = c(rnorm(40, 6, 1), 
        rnorm(20, 8, 1.3))
# Factor 1: variedad de papa (2 niveles)
variedad = gl(2, 30, 60, c('Colombia', 'Pastusa'))
# Factor 2: ferilizante foliar (3 niveles)
ferti = gl(3, 10, 60, c('0cc','5cc','10cc'))

df = data.frame(vol, variedad, ferti)
head(df)
##        vol variedad ferti
## 1 5.439524 Colombia   0cc
## 2 5.769823 Colombia   0cc
## 3 7.558708 Colombia   0cc
## 4 6.070508 Colombia   0cc
## 5 6.129288 Colombia   0cc
## 6 7.715065 Colombia   0cc
medias = tapply(vol, list(variedad, ferti), mean); medias
##               0cc      5cc     10cc
## Colombia 6.074626 6.208622 5.575441
## Pastusa  6.322045 7.988670 8.288192
desv = tapply(vol, list(variedad, ferti), sd); desv
##                0cc      5cc      10cc
## Colombia 0.9537841 1.038073 0.9308092
## Pastusa  0.5273024 1.407274 1.1133786
cv = 100*desv/medias; cv
##                0cc      5cc     10cc
## Colombia 15.701117 16.71987 16.69481
## Pastusa   8.340694 17.61587 13.43331
library(lattice)

bwplot(vol~ ferti | variedad)

xyplot(vol~ ferti | variedad)

Modelo:

Hipotesis: \[H_O:(tau*beta)_ij=0\]

Nota: Las dos siguientes solo son interpretables si no ahy interaccion

mod1 = aov(vol ~ variedad*ferti, df)
summary(mod1)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## variedad        1  37.45   37.45  35.375 2.07e-07 ***
## ferti           2   9.18    4.59   4.334  0.01797 *  
## variedad:ferti  2  15.49    7.75   7.318  0.00154 ** 
## Residuals      54  57.17    1.06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ifelse(mod1$p.value<0.05,"Hay interaccion","No hay interaccion")
## logical(0)

Como el p-value de la interaccion es menor a 0.05, entonces se rechaza la hipotesis de No interaccion

interaction.plot(df$ferti,
                 df$variedad,
                 df$vol)

df_mean = aggregate(list(mvol = df$vol),
          list(variedad = df$variedad,
               ferti = df$ferti),
          mean)
library(ggplot2)

ggplot(df_mean)+
  aes(ferti, mvol, color = variedad)+
  geom_line(aes(group = variedad), size = 1)+
  geom_point(size = 2)

ggplot(df_mean)+
  aes(variedad, mvol, color = ferti)+
  geom_line(aes(group = ferti), size = 1)+
  geom_point(size = 2)

¿Es esta interacción remediable?

df2 = df[df$ferti != '0cc',]
head(df2)
##         vol variedad ferti
## 11 7.224082 Colombia   5cc
## 12 6.359814 Colombia   5cc
## 13 6.400771 Colombia   5cc
## 14 6.110683 Colombia   5cc
## 15 5.444159 Colombia   5cc
## 16 7.786913 Colombia   5cc
mod2 = aov(vol ~ variedad*ferti, df2)
summary(mod2)
##                Df Sum Sq Mean Sq F value  Pr(>F)    
## variedad        1  50.46   50.46  39.088 3.2e-07 ***
## ferti           1   0.28    0.28   0.216   0.645    
## variedad:ferti  1   2.17    2.17   1.685   0.203    
## Residuals      36  46.48    1.29                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df_mean2 = aggregate(list(mvol = df2$vol),
          list(variedad = df2$variedad,
               ferti = df2$ferti),
          mean)
ggplot(df_mean2)+
  aes(ferti, mvol, color = variedad)+
  geom_line(aes(group = variedad), size = 1)+
  geom_point(size = 2)

Al no tener interacción se puede hacer comparación de medias

TukeyHSD(mod2, 'variedad')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = vol ~ variedad * ferti, data = df2)
## 
## $variedad
##                      diff      lwr      upr p adj
## Pastusa-Colombia 2.246399 1.517693 2.975106 3e-07
res = residuals(mod2)
trt = interaction(df2$variedad, df2$ferti)

shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.97349, p-value = 0.4608
bartlett.test(df2$vol, trt)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  df2$vol and trt
## Bartlett's K-squared = 1.6537, df = 3, p-value = 0.6473
df_mean2
##   variedad ferti     mvol
## 1 Colombia   5cc 6.208622
## 2  Pastusa   5cc 7.988670
## 3 Colombia  10cc 5.575441
## 4  Pastusa  10cc 8.288192
# aumento porcentual del rendimiento entre dosis fertilizante para variedad Pastusa
100 * (df_mean2$mvol[4]-df_mean2$mvol[2])/df_mean2$mvol[2]
## [1] 3.749335

Factorial Completo en Bloques Completos al Azar

set.seed(123)
# Respuesta: Volumen de raices
vol = c(rnorm(20, 6, 1), 
        rnorm(10, 8, 1.3))
# Factor 1: variedad de papa (2 niveles)
variedad = gl(2, 15, 30, c('Colombia', 'Pastusa'))
# Factor 2: ferilizante foliar (3 niveles)
ferti = gl(3, 5, 30, c('0cc','5cc','10cc'))
# Bloque: Nivel de alguna caracteristica del suelo
bloq = gl(5, 1, 30, c('--','-','+-','+','++'))

df = data.frame(vol, variedad, ferti, bloq)
head(df)
##        vol variedad ferti bloq
## 1 5.439524 Colombia   0cc   --
## 2 5.769823 Colombia   0cc    -
## 3 7.558708 Colombia   0cc   +-
## 4 6.070508 Colombia   0cc    +
## 5 6.129288 Colombia   0cc   ++
## 6 7.715065 Colombia   5cc   --
library(collapsibleTree)

collapsibleTree(df, c('bloq', 'variedad','ferti'))
table(df$bloq, df$ferti, df$variedad)
## , ,  = Colombia
## 
##     
##      0cc 5cc 10cc
##   --   1   1    1
##   -    1   1    1
##   +-   1   1    1
##   +    1   1    1
##   ++   1   1    1
## 
## , ,  = Pastusa
## 
##     
##      0cc 5cc 10cc
##   --   1   1    1
##   -    1   1    1
##   +-   1   1    1
##   +    1   1    1
##   ++   1   1    1
ggplot(df)+
  aes(ferti, vol, color=variedad)+
  geom_point(size = 2)+
  geom_line(aes(group=variedad))+
  facet_wrap(.~bloq, ncol = 5)

mod3 = aov(vol ~ bloq + variedad*ferti, df)
summary(mod3)
##                Df Sum Sq Mean Sq F value Pr(>F)  
## bloq            4  2.256   0.564   0.418 0.7934  
## variedad        1  5.412   5.412   4.015 0.0588 .
## ferti           2  4.384   2.192   1.626 0.2216  
## variedad:ferti  2  3.522   1.761   1.307 0.2929  
## Residuals      20 26.958   1.348                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(mod3, 'variedad')
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = vol ~ bloq + variedad * ferti, data = df)
## 
## $variedad
##                       diff         lwr      upr     p adj
## Pastusa-Colombia 0.8494454 -0.03487277 1.733764 0.0588366
med_variedad = tapply(df$vol, df$variedad, mean)
med_variedad
## Colombia  Pastusa 
## 6.152384 7.001830
ggplot(df)+
  aes(vol, fill = variedad)+
  geom_density(alpha = 0.3)