Efectos del bloqueo experimental

dtos <- expand.grid(x = seq(7,84,7), y = seq(9,108,9))
plot(dtos)

seq(7,84,7)
##  [1]  7 14 21 28 35 42 49 56 63 70 77 84
set.seed(123)

dtos$rto <- c(c(rnorm(16, 3.2, 0.3), rnorm(16, 3.2, 0.3), rnorm(16, 3.3, 0.3)),
            c(rnorm(16, 3.0, 0.3), rnorm(16, 3.2, 0.3), rnorm(16, 3.5, 0.3)),
            c(rnorm(16, 3.1, 0.3), rnorm(16, 3.2, 0.3), rnorm(16, 3.5, 0.3))
  )

dtos$bloque <- gl(3, 48, 144)
dtos$fert <- gl(3, 16, 144, c('d0','d5','d10'))


dtos %>% 
  ggplot()+
  aes(x,y, fill = rto)+
  geom_tile()

dtos %>% 
  ggplot()+
  aes(x,y, fill = rto, col = bloque)+
  geom_tile(size=1)+
  geom_point(aes(shape = fert), size = 2)+
  theme(legend.position = 'top')

### Gráficos de interacción

interaction.plot(dtos$fert, dtos$bloque, dtos$rto)

interaction.plot(dtos$bloque, dtos$fert, dtos$rto)

mod1 <- aov(rto ~ bloque + fert, data = dtos)
summary(mod1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## bloque        2  0.066  0.0329   0.396    0.674    
## fert          2  3.401  1.7006  20.456 1.63e-08 ***
## Residuals   139 11.556  0.0831                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dtos %>% 
  ggplot() +
  aes(bloque, rto, fill = fert) +
  geom_boxplot()

\[H_0: Los~datos~sigue~una~distribución~normal\\ H_a:Los~datos~sigue~una~distribución~normal\]

res_mod1 <- mod1$residuals
shapiro.test(res_mod1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod1
## W = 0.9901, p-value = 0.4052

\[H_0: \sigma^2_1=\sigma^2_2=...=\sigma^2_k=\sigma^2\\ H_a:\sigma^2_i \neq \sigma^2_j~para~algún~i\neq j\]

bartlett.test(res_mod1, dtos$fert)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  res_mod1 and dtos$fert
## Bartlett's K-squared = 0.22047, df = 2, p-value = 0.8956
  • Como aparentemente el bloqueo no fue significativo (segun el p-valor), tenderia a eliminar los bloques del modelo. (No es correcto)
mod2 <- aov(rto ~ fert, dtos)
summary(mod2)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## fert          2  3.401  1.7006   20.63 1.38e-08 ***
## Residuals   141 11.621  0.0824                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Indice de Moran

# matriz de distancias
d = as.matrix(dist(dtos[,c('x','y')]))
di = 1/d
diag(di) = 0
W = di ## Matriz de pesos

# Modelo con bloques
Moran.I(mod1$residuals, W)
## $observed
## [1] -0.007791896
## 
## $expected
## [1] -0.006993007
## 
## $sd
## [1] 0.007941199
## 
## $p.value
## [1] 0.9198675
# Modelo sin bloques
Moran.I(mod2$residuals, W)
## $observed
## [1] -0.009434653
## 
## $expected
## [1] -0.006993007
## 
## $sd
## [1] 0.007943143
## 
## $p.value
## [1] 0.7585462

p valor> 5% no hay dependencia espacial

Eficiencia de bloqueo

summary(mod1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## bloque        2  0.066  0.0329   0.396    0.674    
## fert          2  3.401  1.7006  20.456 1.63e-08 ***
## Residuals   139 11.556  0.0831                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
H = summary(mod1)[[1]][1,3]/summary(mod1)[[1]][3,3]
H
## [1] 0.395927

\(H<1\) el bloqueo no es eficiente