10.8 Imputación múltiple.
La imputación múltiple consiste en crear múltiples copias del conjunto de datos, donde los valores faltantes en cada copia son imputados utilizando modelos estadísticos. Estos modelos se basan en las relaciones entre las variables en el conjunto de datos y se utilizan para estimar los valores faltantes de manera plausible. Luego de que se han creado múltiples copias completas del conjunto de datos, se realizan análisis separados en cada copia para generar resultados. Los resultados de cada análisis se combinan para obtener un único resultado final que refleje la incertidumbre causada por la imputación de los valores faltantes.
La imputación múltiple es una técnica poderosa para lidiar con los datos faltantes, ya que proporciona resultados más precisos y menos sesgados en comparación con otros métodos que simplemente eliminan las observaciones con valores faltantes. Sin embargo, la imputación múltiple es un proceso computacionalmente intensivo y requiere un conocimiento sólido de la teoría estadística para su implementación efectiva. Esta técnica de imputación fue propuesta por Rubin (1987) y consite en en generar \(M > 1\) conjuntos de valores para los datos faltantes.
En este sentido, suponga entonces que existe un conjunto de \(n\) datos que relaciona dos variables \(X\), \(Y\), a través del siguiente modelo de regresión simple:
\[y_i = \beta x_i + \varepsilon_i\]
Para todo individuo \(i = 1, \ldots, n.\), de tal manera que los errores tienen distribución normal con \(E(\varepsilon) = 0\) y \(Var(\varepsilon) = \sigma ^2\). Suponga que \(Y_{Obs}\) es un vector conteniendo los valores observados para un conjunto de individuos de tamaño \(n_1\), \(Y_{NoObs}\) es el vector de los valores no observados de tamaño \(n_0\), es decir, \(n_1 + n_0 = n\). Además, se asume que sí fue posible observar los valores de la covariable \(X\) para todos los individuos en la muestra.
Al final, el valor imputado corresponderá al promedio de esos \(M\) valores. Dado que el interés es la estimación de la pendiente de la regresión lineal simple \(\beta\), entonces la esperanza estimada al utilizar la metodología de imputación múltiple está dada por:
\[E(\hat{\beta} | Y_{obs}) = E(E(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs})\]
Esta expresión es estimada por el promedio de las \(M\) estimaciones puntuales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dado por:
\[\bar{\hat{\beta}} = \frac{1}{M} \sum_{m = 1} ^ M \hat{\beta}_m\]
La varianza estimada al utilizar la metodología de imputación múltiple está dada por la siguiente expresión:
\[ V(\hat{\beta} | Y_{obs}) = E(V(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs}) + V(E(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs}) \]
La primera parte de la anterior expresión se estima como el promedio de las varianzas muestrales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dado por:
\[\bar{U} = \frac{1}{M} \sum_{m = 1} ^ M Var(\beta)\]
El segundo término se estima como la varianza muestral de las \(M\) estimaciones puntuales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dada por:
\[B = \frac{1}{M-1} \sum_{m = 1} ^ M (\hat{\beta}_m - \bar{\hat{\beta}})\]
Es necesario tener en cuenta un factor de corrección (puesto que \(M\) es finito). Por tanto, la estimación del segundo término viene dada por la siguiente expresión:
\[ (1 + \frac{1}{M}) B \]
Por tanto, la varianza estimada es igual a: \[\hat{V}(\hat{\beta} | Y_{obs}) = \bar{U} + (1 + \frac{1}{M}) B\]
Se ejemplificará la técnica de imputación múltiple para los datos de la encuesta que se utiliza de ejemplo en este texto:
encuesta$Income_imp <- encuesta$Income_missing
encuesta$Employment_imp <- encuesta$Employment_missin
encuesta_obs <- filter(encuesta,!is.na(Income_missing))
encuesta_no_obs <- filter(encuesta, is.na(Income_missing))
n0 <- nrow(encuesta_no_obs)
n1 <- nrow(encuesta_obs)Inicialmente, se extraen los datos a imputar y se calculan los tamaños de los datos observados y no observados.
M <- 10
set.seed(1234)
for (ii in 1:M) {
vp <- paste0("Income_vp_", ii)
vp2 <- paste0("Employment_vp_", ii)
encuesta_temp <-
encuesta_obs %>%
sample_n(size = n1, replace = TRUE)
mod <- lm(Income ~ Zone + Sex + Expenditure, data = encuesta_temp)
encuesta_no_obs[[vp]] <- predict(mod, encuesta_no_obs)
encuesta_obs[[vp]] <- encuesta_obs$Income
}El código anterior ajusta un modelo de regresión lineal con la variable Income como respuesta y las variables Zone, Sex y Expenditure como covariables. Este modelo se utiliza para predecir los valores de la variable Income en la base de datos encuesta_no_obs. Una vez corrido el código anterior, se seleccionan las variables de ingresos y sus valores plausibles (diferentes conjuntos de imputaciones) como se muestra a continuación:
| Income | Income_vp_1 | Income_vp_2 | Income_vp_3 |
|---|---|---|---|
| 409.87 | 550.2195 | 566.0432 | 567.8296 |
| 409.87 | 561.1194 | 529.3457 | 541.8310 |
| 90.92 | 210.5682 | 225.7715 | 163.9913 |
| 90.92 | 221.4682 | 189.0739 | 137.9927 |
| 90.92 | 210.5682 | 225.7715 | 163.9913 |
| 135.33 | 222.6937 | 237.9191 | 178.4083 |
| 135.33 | 222.6937 | 237.9191 | 178.4083 |
| 1539.75 | 784.8579 | 801.1103 | 846.8098 |
| 336.00 | 507.8955 | 472.7913 | 439.7467 |
| 685.48 | 593.0111 | 558.0623 | 540.9473 |
A continuación, se grafica la distribución de los ingresos y los valores plausibles, observándose que las distribuciones son muy similares:
encuesta <- bind_rows(encuesta_obs, encuesta_no_obs)
dat_plot10 <- tidyr::gather(
encuesta %>%
dplyr::select(Zone, Sex, matches("Income_vp_")),
key = "Caso",
value = "Income2",-Zone,-Sex
)
p1 <- ggplot(dat_plot10, aes(x = Income2, col = Caso)) +
geom_density(alpha = 0.2) +
theme_bw() +
theme(legend.position = "bottom") + geom_density(data = encuesta ,
aes(x = Income),
col = "black",
size = 1.2)
p1
Con los valores plausibles enocntrados anteriormente, se procede a definir el diseño muestral utilizado en este ejemplo y así poder hacer la estimación de los parámetros. A continuación, se define el diseño muestral:
library(srvyr)
diseno <- encuesta %>%
as_survey_design(
strata = Stratum,
ids = PSU,
weights = wk,
nest = T
)Con el diseño anterior, se estiman los ingresos medios para cada valor plausible junto ocn su varianza, como se muestra a continuación:
estimacion_vp <- diseno %>%
summarise(
vp1 = survey_mean(Income_vp_1, vartype = c("var")),
vp2 = survey_mean(Income_vp_2, vartype = c("var")),
vp3 = survey_mean(Income_vp_3, vartype = c("var")),
vp4 = survey_mean(Income_vp_4, vartype = c("var")),
vp5 = survey_mean(Income_vp_5, vartype = c("var")),
vp6 = survey_mean(Income_vp_6, vartype = c("var")),
vp7 = survey_mean(Income_vp_7, vartype = c("var")),
vp8 = survey_mean(Income_vp_8, vartype = c("var")),
vp9 = survey_mean(Income_vp_9, vartype = c("var")),
vp10 = survey_mean(Income_vp_10, vartype = c("var"))
)
estimacion_vp| vp1 | vp1_var | vp2 | vp2_var | vp3 | vp3_var | vp4 | vp4_var | vp5 | vp5_var | vp6 | vp6_var | vp7 | vp7_var | vp8 | vp8_var | vp9 | vp9_var | vp10 | vp10_var |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 619.3525 | 845.706 | 617.6058 | 844.5521 | 617.4322 | 867.5133 | 617.748 | 856.7732 | 619.9729 | 857.805 | 617.1098 | 852.6605 | 618.3378 | 860.7981 | 619.2779 | 867.9485 | 616.9401 | 850.131 | 615.7853 | 870.0434 |
A continuación se presentan los datos anteriores discriminado por promedio y varianza:
require(tidyr)
(
estimacion_vp %<>%
tidyr::gather() %>%
separate(key, c("vp", "estimacion")) %>%
mutate(estimacion = ifelse(is.na(estimacion), "promedio", "var")) %>%
spread(estimacion, value) %>%
mutate(vp = 1:10)
)| vp | promedio | var |
|---|---|---|
| 1 | 619.3525 | 845.7060 |
| 2 | 615.7853 | 870.0434 |
| 3 | 617.6058 | 844.5521 |
| 4 | 617.4322 | 867.5133 |
| 5 | 617.7480 | 856.7732 |
| 6 | 619.9729 | 857.8050 |
| 7 | 617.1098 | 852.6605 |
| 8 | 618.3378 | 860.7981 |
| 9 | 619.2779 | 867.9485 |
| 10 | 616.9401 | 850.1310 |
Por último, para obtener la estimación de la media y su varianza utilizando la imputación múltiple, se realizan los siguientes cálculos que se derivan de las expresiones matemáticas antes mostradas:
## [1] 857.3931## [1] 1.645758| Media_vp | Media_vp_se |
|---|---|
| 617.9562 | 29.30005 |
estimacion_var_vp <- diseno %>%
summarise_at(vars(matches("Income_vp")),
survey_var, vartype = "var")Por otro lado, otro parámetro de interés es la varianza de los ingresos. Este parámetro permite medir la variabilidad de los ingresos de los ciudadanos de la base de datos de ejemplo. La forma de estimarla es la misma que para el promedio de los ingresos y se utilizarán los mismo códigos mostrados anteriormente, cambiando el parámetro a estimar:
(
estimacion_var_vp %<>%
tidyr::gather() %>%
separate(key, c("A", "B", "vp", "estimacion")) %>%
mutate(
estimacion = ifelse(is.na(estimacion), "promedio", "var"),
A = NULL,
B = NULL,
vp = as.numeric(vp)
) %>%
spread(estimacion, value)
)| vp | promedio | var |
|---|---|---|
| 1 | 262689.8 | 3074674460 |
| 2 | 263092.0 | 3079895581 |
| 3 | 274370.1 | 3237991671 |
| 4 | 269127.4 | 3164842634 |
| 5 | 270450.0 | 3165285304 |
| 6 | 264992.6 | 3106582542 |
| 7 | 270916.0 | 3175865575 |
| 8 | 276069.5 | 3251836316 |
| 9 | 265060.8 | 3111426636 |
| 10 | 275462.0 | 3258214990 |
Por último, se utilizan las ecuaciones mostradas anteriormente:
Media_var_vp <- mean(estimacion_var_vp$promedio)
Ubar <- mean(estimacion_var_vp$var)
B <- var(estimacion_var_vp$promedio)
var_var_vp <- Ubar + (1 + 1 / M) * B
resultado$var_vp <- Media_var_vp
resultado$var_vp_se <- sqrt(var_var_vp)
cbind(Media_var_vp, var_var_vp)| Media_var_vp | var_var_vp |
|---|---|
| 269223 | 3191037330 |