3.1 Pérdida de unidades experimentales

Comenzaremos por cargar la base de datos y eliminar 20 datos aleatoriamente y ver cómo queda el balance después de la eliminación

set.seed(213)
datos <- ToothGrowth
datos$dose <- as.factor(paste0("D", datos$dose))

datos <- datos[sample(1:60, 40), ]

table(datos$supp)

## 
## OJ VC 
## 21 19

table(datos$dose)

## 
## D0.5   D1   D2 
##   15   17    8

Como podemos apreciar el diseño quedó fuera de balance en varios de los niveles de los factores. Haremos primero ANOVA con suma de cuadrados tipo I:

modelo.1 <- lm(len ~ supp * dose, data = datos)
anova(modelo.1)

## Analysis of Variance Table
## 
## Response: len
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## supp       1  196.45  196.45 15.8233  0.000345 ***
## dose       2 1332.58  666.29 53.6665 3.027e-11 ***
## supp:dose  2   53.28   26.64  2.1457  0.132568    
## Residuals 34  422.12   12.42                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos ver que los resultados son muy parecidos a los del análisis balanceado, pero en este caso ya no podemos demostrar si la interacción es significativa.

Comparemos entonces los resultados con los que se obtienen utilizando suma de cuadrados tipo II:

library(car)

## Loading required package: carData

Anova(modelo.1, type = "II")

## Anova Table (Type II tests)
## 
## Response: len
##            Sum Sq Df F value    Pr(>F)    
## supp       219.10  1 17.6478 0.0001814 ***
## dose      1332.58  2 53.6665 3.027e-11 ***
## supp:dose   53.28  2  2.1457 0.1325679    
## Residuals  422.12 34                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aquí cambiaron un poco los niveles de significancia, pero lo que sí cambió mucho fueron las sumas de cuadrados de supp, pues estas ya no son las que se usan como referencia, ahora estas están calculadas en relación a dose. Veamos entonces cómo $SS-II$ es más conservador que $SS-I$ para estimar la suma de cuadrados.

Tenemos que la varianza total según $SS-I$ es:

# supp   dose    supp:dose residuales
196.45 + 1332.58 + 53.28 + 422.12

## [1] 2004.43

y que la varianza explicada por el modelo con $SS-II$ es:

# supp   dose   supp:dose residuales
219.10 + 1332.58 + 53.28 + 422.12

## [1] 2027.08

Con lo que podemos ver que la varianza total según $SS-II$ es mayor que $SS-I$, lo que altera los valores de $P$ para los diferentes factores, en especial supp que explicó una mayor varianza en comparación con dose en $SS-II$ que en $SS-I$.

3.2 Pérdida de niveles de un factor

Volvamos a cargar la base de datos para eliminar los datos de las combinaciones OJ-D1:

datos.2 <- ToothGrowth
datos.2$dose <- as.factor(paste0("D", datos.2$dose))
datos.2$len[datos.2$supp == "OJ" & datos.2$dose == "D1"] <- NA
datos.2 <- na.omit(datos.2)

Y hagamos un gráfico de cajas para verificar que falten esas combinaciones de tratamientos:

boxplot(len ~ supp + dose, data = datos.2)

Entonces haremos el análisis de este diseño desbalanceado, viendo las diferencias entre $SS$ I y $SS$ II:

modelo.2 <- lm(len ~ supp * dose, data = datos.2)
anova(modelo.2)

## Analysis of Variance Table
## 
## Response: len
##           Df  Sum Sq Mean Sq F value  Pr(>F)    
## supp       1   86.30   86.30  6.7601 0.01257 *  
## dose       2 2401.51 1200.76 94.0628 < 2e-16 ***
## supp:dose  1   71.02   71.02  5.5636 0.02274 *  
## Residuals 45  574.45   12.77                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(modelo.2, type = "II")

## Note: model has aliased coefficients
##       sums of squares computed by model comparison

## Anova Table (Type II tests)
## 
## Response: len
##            Sum Sq Df F value  Pr(>F)    
## supp        66.82  1  5.2346 0.02689 *  
## dose      2401.51  2 94.0628 < 2e-16 ***
## supp:dose   71.02  1  5.5636 0.02274 *  
## Residuals  574.45 45                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.3 Pérdida de observaciones y niveles de un factor

Eliminando los datos aleatoriamente:

datos.3 <- datos.2[sample(1:nrow(datos.2), 40), ]
table(datos.3$supp)

## 
## OJ VC 
## 16 24

table(datos.3$dose)

## 
## D0.5   D1   D2 
##   16    8   16

El análisis con $SS-I$:

modelo.3 <- lm(len ~ supp * dose, data = datos.3)
anova(modelo.3)

## Analysis of Variance Table
## 
## Response: len
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## supp       1   37.76   37.76  3.0742    0.0883 .  
## dose       2 1863.60  931.80 75.8554 1.887e-13 ***
## supp:dose  1   27.01   27.01  2.1989    0.1471    
## Residuals 35  429.94   12.28                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El análisis con $SS-II$:

Anova(modelo.3, type = "II")

## Note: model has aliased coefficients
##       sums of squares computed by model comparison

## Anova Table (Type II tests)
## 
## Response: len
##            Sum Sq Df F value    Pr(>F)    
## supp        30.42  1  2.4764    0.1246    
## dose      1863.60  2 75.8554 1.887e-13 ***
## supp:dose   27.01  1  2.1989    0.1471    
## Residuals  429.94 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos ver entonces que al haber perdido grados de libertad perdimos la capacidad de rechazar la hipótesis nula del tratamiento supp y la interacción supp : dose en ambos tipos de suma de cuadrados, pero el valor de $P$ para supp fue mayor en $SS-I$.

Ahora, finalmente veamos qué pasa si especificamos el modelo invirtiendo el orden de los factores:

modelo.4 <- lm(len ~ dose * supp, data = datos.3)
Anova(modelo.4, type = "II")

## Note: model has aliased coefficients
##       sums of squares computed by model comparison

## Anova Table (Type II tests)
## 
## Response: len
##            Sum Sq Df F value    Pr(>F)    
## dose      1863.60  2 75.8554 1.887e-13 ***
## supp        30.42  1  2.4764    0.1246    
## dose:supp   27.01  1  2.1989    0.1471    
## Residuals  429.94 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1