Clase
Stefania Hernández Azuero
2023-05-26
En la clase de hoy se analizará los diseños desbalanceados y cómo crear un análisis de varianza adecuado para cada tipo de diseño.
#Análisis de varianza desbalanceado
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
acido = gl(30, 40, 120, c('C0', 'C1', 'C2'))
datos = data.frame(acido, porc_germ)
head(datos)## acido porc_germ
## 1 C0 56.63715
## 2 C0 58.61894
## 3 C0 69.35225
## 4 C0 60.42305
## 5 C0 60.77573
## 6 C0 70.29039table(datos$acido)##
## C0 C1 C2
## 40 40 40datos_des = datos[-c(50, 111, 120),]
table(datos_des$acido)##
## C0 C1 C2
## 40 39 38En esta parte, se crean los datos de un diseño balanceado y uno desbalanceado con un sólo factor. Así, en el diseño balanceado todos los niveles cuentan con 40 datos. En cambio, en el diseño desbalanceado en el nivel ‘C1’ falta un (1) dato y en el diseño ‘C2’ faltan dos (2) datos.
#Análisis de varianza balanceado
mod1 = aov(porc_germ~acido, datos)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## acido 2 7835 3918 98.15 <2e-16 ***
## Residuals 117 4670 40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(datos$porc_germ ~ datos$acido)
Usamos la función ‘aov’ como usualmente y procedemos a obtener el
boxplot, en el cuál evidenciamos la existencia de un dato atípico en el
nivel ‘C1’.
#Análisis de varianza balanceado (Con datos desbalanceados)
mod2 = aov(porc_germ~acido, datos_des)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## acido 2 7898 3949 98.39 <2e-16 ***
## Residuals 114 4576 40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En este ejemplo, usamos la función ‘aov’ para obtener el análisis de varianza del experimento desbalanceado. Cabe destacar que ‘Aov’no está diseñado para modelos desbalanceado sino balanceados. Por lo tanto, a función que se debería usar para experimentos desbalanceados es ’Anova’ con la librería ‘car’.
mod3 = lm(porc_germ~acido, datos_des)
library(car)## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3mod3_res = Anova(mod3, type = 'II')
mod3_res## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## acido 7898.3 2 98.392 < 2.2e-16 ***
## Residuals 4575.6 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No obstante, como se evidencia en este caso a pesar de usamos la función ‘Anova’ los valores que se nos entregan en esta tabla son bastante similares a los obtenidos en el análisis anterior (mod2) con la función ‘Aov’. Esto se debe a que cuando el diseño es de un solo factor, se puede usar ‘Aov’ o ‘Anova’ y no va a cambiar significativamente el resultado. No obstante, cuando se empieza a aumentar la cantidad de factores, ahí si se debería usar únicamente la función ‘Anova’. Más adelante se explicará la importancia de usar siempre el ‘Type = ’II’’ en la función.
#Ejemplo Bloques (desbalanceado)
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
bloq = gl(3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl(4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ)
datos_des = datos[-c(50, 111, 120),]
datos_des## acido bloq porc_germ
## 1 C0 B0 56.63715
## 2 C0 B0 58.61894
## 3 C0 B0 69.35225
## 4 C0 B0 60.42305
## 5 C0 B0 60.77573
## 6 C0 B0 70.29039
## 7 C0 B0 62.76550
## 8 C0 B0 52.40963
## 9 C0 B0 55.87888
## 10 C0 B0 57.32603
## 11 C1 B0 67.34449
## 12 C1 B0 62.15888
## 13 C1 B0 62.40463
## 14 C1 B0 60.66410
## 15 C1 B0 56.66495
## 16 C1 B0 70.72148
## 17 C1 B0 62.98710
## 18 C1 B0 48.20030
## 19 C1 B0 64.20814
## 20 C1 B0 57.16325
## 21 C2 B0 53.59306
## 22 C2 B0 58.69215
## 23 C2 B0 53.84397
## 24 C2 B0 55.62665
## 25 C2 B0 56.24976
## 26 C2 B0 49.87984
## 27 C2 B0 65.02672
## 28 C2 B0 60.92024
## 29 C2 B0 53.17118
## 30 C2 B0 67.52289
## 31 C3 B0 62.55879
## 32 C3 B0 58.22957
## 33 C3 B0 65.37075
## 34 C3 B0 65.26880
## 35 C3 B0 64.92949
## 36 C3 B0 64.13184
## 37 C3 B0 63.32351
## 38 C3 B0 59.62853
## 39 C3 B0 58.16422
## 40 C3 B0 57.71717
## 41 C0 B1 65.13705
## 42 C0 B1 68.54458
## 43 C0 B1 61.14223
## 44 C0 B1 85.18269
## 45 C0 B1 78.45573
## 46 C0 B1 62.13824
## 47 C0 B1 67.17981
## 48 C0 B1 66.73341
## 49 C0 B1 75.45976
## 51 C1 B1 71.77323
## 52 C1 B1 69.80017
## 53 C1 B1 69.69991
## 54 C1 B1 79.58022
## 55 C1 B1 68.41960
## 56 C1 B1 80.61529
## 57 C1 B1 59.15873
## 58 C1 B1 74.09230
## 59 C1 B1 70.86698
## 60 C1 B1 71.51159
## 61 C2 B1 72.65748
## 62 C2 B1 66.48374
## 63 C2 B1 67.66755
## 64 C2 B1 62.86997
## 65 C2 B1 62.49746
## 66 C2 B1 72.12470
## 67 C2 B1 73.13747
## 68 C2 B1 70.37103
## 69 C2 B1 76.45587
## 70 C2 B1 84.35059
## 71 C3 B1 66.56278
## 72 C3 B1 53.83582
## 73 C3 B1 77.04017
## 74 C3 B1 65.03559
## 75 C3 B1 65.18394
## 76 C3 B1 77.17900
## 77 C3 B1 68.00659
## 78 C3 B1 61.45498
## 79 C3 B1 71.26912
## 80 C3 B1 69.02776
## 81 C0 B2 80.04611
## 82 C0 B2 83.08224
## 83 C0 B2 77.03472
## 84 C0 B2 85.15501
## 85 C0 B2 78.23611
## 86 C0 B2 82.65426
## 87 C0 B2 88.77471
## 88 C0 B2 83.48145
## 89 C0 B2 77.39255
## 90 C0 B2 89.19046
## 91 C1 B2 87.94803
## 92 C1 B2 84.38718
## 93 C1 B2 81.90985
## 94 C1 B2 74.97675
## 95 C1 B2 90.88522
## 96 C1 B2 75.19792
## 97 C1 B2 97.49866
## 98 C1 B2 92.26089
## 99 C1 B2 78.11440
## 100 C1 B2 71.78863
## 101 C2 B2 74.31675
## 102 C2 B2 82.05507
## 103 C2 B2 78.02646
## 104 C2 B2 77.21966
## 105 C2 B2 72.38705
## 106 C2 B2 79.63978
## 107 C2 B2 73.72076
## 108 C2 B2 66.65646
## 109 C2 B2 76.95819
## 110 C2 B2 87.35197
## 112 C3 B2 84.86371
## 113 C3 B2 67.05694
## 114 C3 B2 79.55550
## 115 C3 B2 84.15526
## 116 C3 B2 82.40923
## 117 C3 B2 80.84541
## 118 C3 B2 74.87435
## 119 C3 B2 73.20237table(datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 10 10 10
## B1 9 10 10 10
## B2 10 10 10 8Ahora, procedemos a crear los datos desbalanceados para un diseño factorial simple en bloques. Para este caso todos los bloques con los respectivos niveles tienen de a diez (10) datos cada uno, a excepción del Bloque 1, nivel ‘C0’ que le falta un (1) datos y el bloque 2, nivel ‘C3’ que le faltan dos (2) datos.
#Corriendo el diseño con Aov
mod1 = aov(porc_germ ~ bloq *acido, datos_des)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## bloq 2 7898 3949 102.076 <2e-16 ***
## acido 3 248 83 2.133 0.100
## bloq:acido 6 266 44 1.145 0.342
## Residuals 105 4062 39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Corriendo el diseño con lm y Anova
mod2 = lm(porc_germ ~ bloq * acido, datos_des)
mod2_res = Anova(mod2, type = 'II')
mod2_res## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como se evidencia los análisis obtenidos con ‘Aov’ y ‘Anova’ (mod1 y mod2, respectivamente) son bastantes similares y sus valores no difieren significativamente. Esto puede deberse a que solo hay tres datos que están faltando, para ello eliminaremos más datos y veremos si ya hay diferencias significativas en los análisis.
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
bloq = gl(3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl(4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ)
datos_des = datos[-sample(120, 5),]
datos_des## acido bloq porc_germ
## 1 C0 B0 56.63715
## 2 C0 B0 58.61894
## 3 C0 B0 69.35225
## 4 C0 B0 60.42305
## 5 C0 B0 60.77573
## 6 C0 B0 70.29039
## 7 C0 B0 62.76550
## 8 C0 B0 52.40963
## 9 C0 B0 55.87888
## 10 C0 B0 57.32603
## 11 C1 B0 67.34449
## 12 C1 B0 62.15888
## 13 C1 B0 62.40463
## 15 C1 B0 56.66495
## 17 C1 B0 62.98710
## 18 C1 B0 48.20030
## 19 C1 B0 64.20814
## 20 C1 B0 57.16325
## 21 C2 B0 53.59306
## 22 C2 B0 58.69215
## 23 C2 B0 53.84397
## 24 C2 B0 55.62665
## 25 C2 B0 56.24976
## 26 C2 B0 49.87984
## 27 C2 B0 65.02672
## 28 C2 B0 60.92024
## 29 C2 B0 53.17118
## 30 C2 B0 67.52289
## 31 C3 B0 62.55879
## 32 C3 B0 58.22957
## 34 C3 B0 65.26880
## 35 C3 B0 64.92949
## 36 C3 B0 64.13184
## 37 C3 B0 63.32351
## 38 C3 B0 59.62853
## 39 C3 B0 58.16422
## 40 C3 B0 57.71717
## 41 C0 B1 65.13705
## 42 C0 B1 68.54458
## 43 C0 B1 61.14223
## 44 C0 B1 85.18269
## 45 C0 B1 78.45573
## 46 C0 B1 62.13824
## 47 C0 B1 67.17981
## 48 C0 B1 66.73341
## 49 C0 B1 75.45976
## 50 C0 B1 69.41642
## 51 C1 B1 71.77323
## 52 C1 B1 69.80017
## 53 C1 B1 69.69991
## 54 C1 B1 79.58022
## 55 C1 B1 68.41960
## 56 C1 B1 80.61529
## 57 C1 B1 59.15873
## 58 C1 B1 74.09230
## 59 C1 B1 70.86698
## 60 C1 B1 71.51159
## 61 C2 B1 72.65748
## 62 C2 B1 66.48374
## 63 C2 B1 67.66755
## 64 C2 B1 62.86997
## 65 C2 B1 62.49746
## 66 C2 B1 72.12470
## 67 C2 B1 73.13747
## 68 C2 B1 70.37103
## 69 C2 B1 76.45587
## 70 C2 B1 84.35059
## 71 C3 B1 66.56278
## 72 C3 B1 53.83582
## 73 C3 B1 77.04017
## 74 C3 B1 65.03559
## 75 C3 B1 65.18394
## 76 C3 B1 77.17900
## 77 C3 B1 68.00659
## 78 C3 B1 61.45498
## 79 C3 B1 71.26912
## 80 C3 B1 69.02776
## 81 C0 B2 80.04611
## 82 C0 B2 83.08224
## 83 C0 B2 77.03472
## 84 C0 B2 85.15501
## 85 C0 B2 78.23611
## 86 C0 B2 82.65426
## 88 C0 B2 83.48145
## 89 C0 B2 77.39255
## 90 C0 B2 89.19046
## 91 C1 B2 87.94803
## 92 C1 B2 84.38718
## 93 C1 B2 81.90985
## 94 C1 B2 74.97675
## 95 C1 B2 90.88522
## 96 C1 B2 75.19792
## 97 C1 B2 97.49866
## 98 C1 B2 92.26089
## 99 C1 B2 78.11440
## 100 C1 B2 71.78863
## 101 C2 B2 74.31675
## 102 C2 B2 82.05507
## 103 C2 B2 78.02646
## 104 C2 B2 77.21966
## 105 C2 B2 72.38705
## 106 C2 B2 79.63978
## 107 C2 B2 73.72076
## 108 C2 B2 66.65646
## 109 C2 B2 76.95819
## 110 C2 B2 87.35197
## 111 C3 B2 75.39722
## 112 C3 B2 84.86371
## 113 C3 B2 67.05694
## 114 C3 B2 79.55550
## 115 C3 B2 84.15526
## 116 C3 B2 82.40923
## 118 C3 B2 74.87435
## 119 C3 B2 73.20237
## 120 C3 B2 71.80697table(datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 8 10 9
## B1 10 10 10 10
## B2 9 10 10 9Se crea una tabla de datos, dónde están faltando cinco datos en el diseño. Estos datos hace falta en el Bloque 0, nivel ‘C1’ y ‘C3’ y en el Bloque 3, nivel ‘C0’ y ‘C3’.
mod1 = aov(porc_germ ~ bloq *acido, datos_des)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## bloq 2 7478 3739 97.858 <2e-16 ***
## acido 3 238 79 2.073 0.108
## bloq:acido 6 276 46 1.203 0.311
## Residuals 103 3936 38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod2 = lm(porc_germ ~ bloq * acido, datos_des)
mod2_res = Anova(mod2, type = 'II')
mod2_res## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7399.0 2 96.8208 <2e-16 ***
## acido 237.7 3 2.0734 0.1083
## bloq:acido 275.8 6 1.2030 0.3108
## Residuals 3935.6 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A pesar de haber aumentado la cantidad de datos faltantes, se observa que los valores obtenidos en el análisis de varianza con la función ‘Aov’ y ‘Anova’ no son significativamente diferentes, especialmente en los P-valores.
#Ajustando desbalanceado con diferente orden
mod3 = lm(porc_germ ~ bloq + acido + bloq:acido, datos_des)
Anova(mod3, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7399.0 2 96.8208 <2e-16 ***
## acido 237.7 3 2.0734 0.1083
## bloq:acido 275.8 6 1.2030 0.3108
## Residuals 3935.6 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ acido + bloq + bloq:acido, datos_des)
Anova(mod3, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## acido 237.7 3 2.0734 0.1083
## bloq 7399.0 2 96.8208 <2e-16 ***
## acido:bloq 275.8 6 1.2030 0.3108
## Residuals 3935.6 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ bloq:acido + acido + bloq, datos_des)
Anova(mod3, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## acido 237.7 3 2.0734 0.1083
## bloq 7399.0 2 96.8208 <2e-16 ***
## bloq:acido 275.8 6 1.2030 0.3108
## Residuals 3935.6 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ bloq:acido + bloq + acido, datos_des)
Anova(mod3, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7399.0 2 96.8208 <2e-16 ***
## acido 237.7 3 2.0734 0.1083
## bloq:acido 275.8 6 1.2030 0.3108
## Residuals 3935.6 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En esta parte, estamos evaluando que pasaría si colocasemos el factor, bloques e interacción del experimentos en distintos ordenes a la hora de crear el análisis de varianza. Como se observa, a pesar de que se cambió el órden de los datos en cada Anova, los resultados son idénticos en cada tabla. Esto se debe al ‘Type=’II’’, el cuál arregla la suma de cuadrados. Tradicionalmente, siempre se ha usado el Tipo I (para experimentos balanceados) y el tipo III (para experimentos desbalanceados), pero el Tipo II nos permite en estos casos obtener resultados certeros sin importar el órden de los componentes en el diseño.
Lo anteriormente dicho, se puede verificar en el artículo “ANOVA for unbalanced data: Use type II instead of type III sums of squares.”
#Ahora, con una covariable
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
diam_med = sort(rnorm(120, 12, 1.3))
bloq = gl(3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl(4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ, diam_med)
datos_des = datos[-sample(120, 5),]
table(datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 9 9 10
## B1 10 10 9 9
## B2 10 9 10 10Para calcular diámetro medio de algo que no es circular, se usa la media geométrica. Por ejemplo, raíz de la multiplicación de las medidas.
#Diseño Factorial Simple en Bloques completamente Generalizados y al Azar, Desbalanceado y con la técnica Análisis de Covarianza.
FS –> Acido BC –> 3Bloques GA –> Rep. en los bloques D –> Repetición dif. Por esto ya no sería Anova sino Ancova
mod1 = lm(porc_germ ~ diam_med + bloq + acido + bloq:acido,
datos_des)
Anova(mod1, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## diam_med 6.0 1 0.1574 0.6924
## bloq 1069.4 2 14.0706 4.013e-06 ***
## acido 151.8 3 1.3315 0.2683
## bloq:acido 235.4 6 1.0325 0.4087
## Residuals 3876.1 102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Procedemos a crear el Ancova para este experimento desbalanceado que ahora incluye una covariable.
#Ahora, con un dato N/A
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
diam_med = sort(rnorm(120, 12, 1.3))
bloq = gl(3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl(4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ, diam_med)
datos_des = datos
datos_des[sample(120, 5), 'porc_germ'] = NA
table(datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 10 10 10
## B1 10 10 10 10
## B2 10 10 10 10Se crean nuevamente los datos, pero esta vez el desbalanceamiento lo denotaremos con NA en las filas que se hayan perdido.
tapply(datos_des$porc_germ,
datos_des$acido,
mean)## C0 C1 C2 C3
## 70.96384 NA NA NAComo se evidencia, a la hora de obtener la media de los niveles, sólo se nos entrega la del nivel ‘C0’ y el resto de niveles tiene como media NA. Esto se debe a que R no puede dar la media cuando están los valores NA en las filas de los datos perdidos. Para ellos, tendremos que eliminar o hacer que R ignore estos NA.
tapply(datos_des$porc_germ,
datos_des$acido,
mean, na.rm = TRUE)## C0 C1 C2 C3
## 70.96384 72.04659 68.44023 69.07068El problema anterior se soluciona al utilizar la función ‘na.rm = True’. Esta función se usa para excluir valores faltantes en los datos cuando están codificados como NA, ya con esta función podemos obtener las medias. Como otra opción está imputar los datos faltantes, y así también podría arreglarse el problema del desbalanceo.
#Punto 2. ¿Qué pasa en un análisis de varianza cuando aparece NA en los datos faltantes en vez de eliminar la fila?
Cuando se desea correr un Anova pero hay datos faltantes, no habría problema siempre cuando sean pocos y halla al menos un valor en cada fila para cada conjunto de datos y así se pueda ajustar el modelo. No obstante, cuando los valores faltantes son NA, se pueden uso de diversas técnicas y funciones para omitirlos y obtener tablas de análisis de varianza acertadas. Algunas de las formas para solucionar este problema que encontré son:
- Listwise Deletion (eliminación por lista)
Aquí se tiene un ejemplo de unos datos, que en la variable Y tiene varios NA valores. Para este método, creamos un Anova con la función ‘lm’ y lo corremos, a través de la Listwise Deletion sólo se tomará en cuenta los valores Y observados y se omitirán los NA. Cabe destacar que esta es una configuración estándar en la mayoría de programas estadísticos. Lo anteriormente dicho se verifica en la tabla ‘mod1a’ dónde en la parte inferior vemos el mensaje ‘(31 observations deleted due to missingness)’, por lo tanto, los valores NA fueron omitidos correctamente.
library(readxl)
dat = read_excel("C:/Users/stefa/Downloads/dat.xlsx")
dat$Y <- as.numeric(gsub("\\.", "", dat$Y))## Warning: NAs introducidos por coerciónmod1a <- lm(Y ~ 1 + group, data = dat)
summary(mod1a)##
## Call:
## lm(formula = Y ~ 1 + group, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.526e+16 -1.969e+16 -5.756e+15 2.186e+16 3.481e+16
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.709e+16 1.255e+16 2.159 0.0487 *
## groupB 1.927e+16 1.699e+16 1.134 0.2758
## groupC 1.008e+16 1.699e+16 0.593 0.5627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.806e+16 on 14 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.0842, Adjusted R-squared: -0.04663
## F-statistic: 0.6436 on 2 and 14 DF, p-value: 0.5403- Multiple imputation (Imputación múltiple)
Este método es uno de los más confiables que podemos usar, ya que reemplaza los valores faltantes a través del modelo de imputación. Asimismo, este método tiene como beneficio que nos permite tener en cuenta covariables que son relevantes para los datos que faltan aunque no tenga una importancia significativa en el Anova. En este método se tiene la siguiente serie de pasos para proceder:
- Los datos que faltan se imputan varias veces.
- Se analizan por separado los conjuntos de datos imputados.
- Se forma un conjunto final de estimaciones e inferencias.
Ejemplo (con datos usados para el método de eliminación por lista)
#1. Se cargan las librerías necesarias para la imputación múltiple
library(mice)## Warning: package 'mice' was built under R version 4.2.3##
## Attaching package: 'mice'## The following object is masked from 'package:stats':
##
## filter## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(mitml)## Warning: package 'mitml' was built under R version 4.2.3## *** This is beta software. Please report any bugs!
## *** See the NEWS file for recent changes.#2. Se corre la imputación múltiple y se determina la cantidad de imputaciones a realizar ('m=')
imp <- mice(data = dat, method = "norm", m = 10)##
## iter imp variable
## 1 1 Y
## 1 2 Y
## 1 3 Y
## 1 4 Y
## 1 5 Y
## 1 6 Y
## 1 7 Y
## 1 8 Y
## 1 9 Y
## 1 10 Y
## 2 1 Y
## 2 2 Y
## 2 3 Y
## 2 4 Y
## 2 5 Y
## 2 6 Y
## 2 7 Y
## 2 8 Y
## 2 9 Y
## 2 10 Y
## 3 1 Y
## 3 2 Y
## 3 3 Y
## 3 4 Y
## 3 5 Y
## 3 6 Y
## 3 7 Y
## 3 8 Y
## 3 9 Y
## 3 10 Y
## 4 1 Y
## 4 2 Y
## 4 3 Y
## 4 4 Y
## 4 5 Y
## 4 6 Y
## 4 7 Y
## 4 8 Y
## 4 9 Y
## 4 10 Y
## 5 1 Y
## 5 2 Y
## 5 3 Y
## 5 4 Y
## 5 5 Y
## 5 6 Y
## 5 7 Y
## 5 8 Y
## 5 9 Y
## 5 10 Y## Warning: Number of logged events: 1#3. Los datos faltantes ya han sido reemplazados y por ende, se prosigue a crear un conjunto de datos completos.
implist <- mids2mitml.list(imp)#4. Ahora, se proceder a crear el ANOVA y estimar los parámetros.
fit2 <- with(implist, lm(Y ~ 1 + group))
testEstimates(fit2)##
## Call:
##
## testEstimates(model = fit2)
##
## Final parameter estimates and inferences obtained from 10 imputed data sets.
##
## Estimate Std.Error t.value df P(>|t|) RIV FMI
## (Intercept) 2.677e+16 1.172e+16 2.284 220.544 0.023 0.253 0.209
## groupB 1.721e+16 1.592e+16 1.081 489.436 0.280 0.157 0.139
## groupC 1.216e+16 1.622e+16 0.750 324.372 0.454 0.200 0.172
##
## Unadjusted hypothesis test as appropriate in larger samples.#Veamos el modelo con los coeficientes correspondientes a las diferencias entre grupos.
fit2.reduced <- with(implist, lm(Y ~ 1))
testModels(fit2, fit2.reduced, method = "D1")##
## Call:
##
## testModels(model = fit2, null.model = fit2.reduced, method = "D1")
##
## Model comparison calculated from 10 imputed data sets.
## Combination method: D1
##
## F.value df1 df2 P(>F) RIV
## 0.600 2 452.496 0.549 0.191
##
## Unadjusted hypothesis test as appropriate in larger samples.En conclusión, existen varios métodos eliminiar valores faltantes (NA) en R. Estos dos me parecieron bastante sencillos y su complejidad puede variar de acuerdo al tipo de diseño. Si se corre el ANOVA con valores NA puede que no se obtengan valores poco confiables (siempre y cuando sean pocos los valores faltantes), pero en mi opinión, es mejor utilizar algún método para eliminar o reemplazar dichos valores NA y así tomar decisiones acertadas y obtener análisis de varianza correctos.