Diseño de bloques aleatorizados

library(faraway)
## Warning: package 'faraway' was built under R version 4.2.3
head(rabbit)
##   treat gain block
## 1     f 42.2    b1
## 2     b 32.6    b1
## 3     c 35.2    b1
## 4     c 40.9    b2
## 5     a 40.1    b2
## 6     b 38.1    b2
?rabbit
## starting httpd help server ... done

A nutritionist studied the effects of six diets, on weight gain of domestic rabbits. From past experience with sizes of litters, it was felt that only 3 uniform rabbits could be selected from each available litter. There were ten litters available forming blocks of size three.

Visualización del modelo

matrix(rabbit$treat, nrow = 3)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] "f"  "c"  "c"  "a"  "e"  "b"  "d"  "a"  "d"  "f"  
## [2,] "b"  "a"  "f"  "e"  "c"  "f"  "a"  "e"  "b"  "d"  
## [3,] "c"  "b"  "d"  "c"  "d"  "e"  "b"  "f"  "e"  "a"
matrix(rabbit$block, nrow=3)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"
## [2,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"
## [3,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"
matrix(rabbit$gain, nrow=3)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 42.2 40.9 34.6 44.9 32.0 37.3 37.9 44.0 27.5  41.7
## [2,] 32.6 40.1 34.3 40.8 40.9 42.8 45.2 38.5 30.6  42.3
## [3,] 35.2 38.1 37.5 43.9 37.3 40.5 40.6 51.9 20.6  37.3

Tabla de contingencia

xtabs(gain~treat+block, data=rabbit)
##      block
## treat   b1  b10   b2   b3   b4   b5   b6   b7   b8   b9
##     a  0.0 37.3 40.1  0.0 44.9  0.0  0.0 45.2 44.0  0.0
##     b 32.6  0.0 38.1  0.0  0.0  0.0 37.3 40.6  0.0 30.6
##     c 35.2  0.0 40.9 34.6 43.9 40.9  0.0  0.0  0.0  0.0
##     d  0.0 42.3  0.0 37.5  0.0 37.3  0.0 37.9  0.0 27.5
##     e  0.0  0.0  0.0  0.0 40.8 32.0 40.5  0.0 38.5 20.6
##     f 42.2 41.7  0.0 34.3  0.0  0.0 42.8  0.0 51.9  0.0
plot(gain~treat+block, rabbit, col = 050)

ANOVA

lm.rabbit <- lm(gain~block+treat, data=rabbit)
ANOVA1<-anova(lm.rabbit)
ANOVA1
## Analysis of Variance Table
## 
## Response: gain
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## block      9 730.39  81.154  8.0738 0.0002454 ***
## treat      5 158.73  31.745  3.1583 0.0381655 *  
## Residuals 15 150.77  10.052                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Block: El estadístico F asciende a 8.07 con un p valor de 0.0002454 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Por tanto, los valores medios de gain son diferentes por tipo de treat.

Treat: El estadístico F asciende a 3.15 con un p valor de 0.03 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Por tanto, los valores medios de gain son diferentes por tipo de treat.

Observamos que tanto el bloque como el tratamiento son significativos.

plot(ANOVA1)

ANOVA1$`Pr(>F)`
## [1] 0.0002454165 0.0381654787           NA

¿Qué tabla es adecuada para probar el efecto del tratamiento o el efecto de bloque?

summary(lm.rabbit)
## 
## Call:
## lm(formula = gain ~ block + treat, data = rabbit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9583 -1.6146 -0.6083  1.9396  4.3028 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 36.01389    2.58863  13.912 5.59e-10 ***
## blockb10     3.29722    2.79604   1.179  0.25667    
## blockb2      4.13333    2.69433   1.534  0.14583    
## blockb3     -1.80278    2.69433  -0.669  0.51360    
## blockb4      8.79444    2.79604   3.145  0.00667 ** 
## blockb5      2.30556    2.79604   0.825  0.42253    
## blockb6      5.40833    2.69433   2.007  0.06309 .  
## blockb7      5.77778    2.79604   2.066  0.05651 .  
## blockb8      9.42778    2.79604   3.372  0.00419 ** 
## blockb9     -7.48056    2.79604  -2.675  0.01729 *  
## treatb      -1.74167    2.24182  -0.777  0.44930    
## treatc       0.40000    2.24182   0.178  0.86078    
## treatd       0.06667    2.24182   0.030  0.97667    
## treate      -5.22500    2.24182  -2.331  0.03413 *  
## treatf       3.30000    2.24182   1.472  0.16169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.17 on 15 degrees of freedom
## Multiple R-squared:  0.855,  Adjusted R-squared:  0.7197 
## F-statistic: 6.318 on 14 and 15 DF,  p-value: 0.000518
        Estimate Std. Error t value Pr(>|t|)

(Intercept) 36.01389 2.58863 13.912 5.59e-10 * blockb4 8.79444 2.79604 3.145 0.00667 blockb6 5.40833 2.69433 2.007 0.06309 .
blockb7 5.77778 2.79604 2.066 0.05651 .
blockb8 9.42778 2.79604 3.372 0.00419 ** blockb9 -7.48056 2.79604 -2.675 0.01729 *
treate -5.22500 2.24182 -2.331 0.03413 *

Éstos son los tratamientos y bloques individualmente significativos.

El R2 es 0.855; por tanto, los tratamientos y bloque explican el 85% de la variabilidad del Gain (R2 < 75%; por tanto, el modelo es adecuado).

plot(lm.rabbit)

residuos<-lm.rabbit$residuals
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.93647, p-value = 0.07311

Normalidad.

Se verifican los supuestos del modelo: Linealidad, independencia, homocedasticidad, ausencia de atípicos influyentes. Por lo tanto, el modelo es válido.

lm.rabbitr <- lm(gain~treat, data=rabbit)
anova(lm.rabbitr)
## Analysis of Variance Table
## 
## Response: gain
##           Df Sum Sq Mean Sq F value Pr(>F)
## treat      5 293.38  58.676  1.8864 0.1342
## Residuals 24 746.51  31.104

El treat por sí solo no es significativo, pues su p valor es 0.13>0.05.

summary(lm.rabbitr)
## 
## Call:
## lm(formula = gain ~ treat, data = rabbit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.880  -3.050   1.200   2.825   9.320 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   42.300      2.494  16.960 7.28e-15 ***
## treatb        -6.460      3.527  -1.831   0.0795 .  
## treatc        -3.200      3.527  -0.907   0.3733    
## treatd        -5.800      3.527  -1.644   0.1131    
## treate        -7.820      3.527  -2.217   0.0363 *  
## treatf         0.280      3.527   0.079   0.9374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.577 on 24 degrees of freedom
## Multiple R-squared:  0.2821, Adjusted R-squared:  0.1326 
## F-statistic: 1.886 on 5 and 24 DF,  p-value: 0.1342

Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.300 2.494 16.960 7.28e-15 ** treatb -6.460 3.527 -1.831 0.0795 .
treate -7.820 3.527 -2.217 0.0363

El tratamiento b y e son individualmente significativos al 10%.

El estadístico F de significación conjunta para los tratamientos asciende a 1.88 con p valor 0.1342>0.05. Por tanto, los tratamientos no son conjuntamente significativos.

¿Cuál es la dfierencia entre lm.rabbit y lm.rabbitr? El modelo restringido no incluye los tratamientos –> comparando los modelos con un contraste Restringido vs no restringido, podemos evaluar la significación de incluir los bloques en el modelo.

anova(lm.rabbit, lm.rabbitr)
## Analysis of Variance Table
## 
## Model 1: gain ~ block + treat
## Model 2: gain ~ treat
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     15 150.77                                  
## 2     24 746.51 -9   -595.74 6.5854 0.0007602 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p valor es 0.0007 < 0.05. Tenemos evidencia de que incluir los bloques aporta información estadísticamente significativa al modelo.

Concluimos el análisis de bloques es bueno para este caso. Hay que notar que hay significancia entre los modelos.

INFO DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 59 × 3
##    package  loadedversion source        
##    <chr>    <chr>         <chr>         
##  1 boot     1.3-28        CRAN (R 4.2.2)
##  2 bslib    0.4.2         CRAN (R 4.2.2)
##  3 cachem   1.0.6         CRAN (R 4.2.2)
##  4 callr    3.7.3         CRAN (R 4.2.3)
##  5 cli      3.6.0         CRAN (R 4.2.2)
##  6 crayon   1.5.2         CRAN (R 4.2.3)
##  7 devtools 2.4.5         CRAN (R 4.2.3)
##  8 digest   0.6.31        CRAN (R 4.2.2)
##  9 ellipsis 0.3.2         CRAN (R 4.2.2)
## 10 evaluate 0.20          CRAN (R 4.2.2)
## # ℹ 49 more rows