Imputación múltiple para una regresión simple
Andrés Gutiérrez
Octubre de 2015
Introducción
Suponga 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\). Bajo esta perspectiva, suponga que la variable dependiente \(Y\) sólo pudo ser observada para un conjunto de individuos de tamaño \(n_1\), mientras que para los \(n_0\) individuos restantes, no existen datos para la variable respuesta. 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.
Vamos a simular un conjunto de \(n = 500\) datos con una pendiente \(\beta = 10\) y con una dispersión de \(\sigma = 2\). A su vez, el conjunto de datos tendrá \(n_0 = 200\) valores faltantes en la variable respuesta. La siguiente función genera la estructura de los datos
generar <- function(n = 500, n_0 = 200, beta = 10, sigma = 2){
x <- runif(n)
mu <- beta * x
y <- mu + rnorm(n, mean = 0, sd = sigma)
datos <- data.frame(x = x, y = y)
faltantes <- sample(n, n_0)
datos$faltantes <- "No"
datos$faltantes[faltantes] <- "Si"
datos$y.per <- y
datos$y.per[faltantes] <- NA
return(datos)
}
datos <- generar()
head(datos)## x y faltantes y.per
## 1 0.7275323 5.854265 No 5.854265
## 2 0.8481436 8.448366 Si NA
## 3 0.1442695 4.635120 No 4.635120
## 4 0.8738171 8.380226 No 8.380226
## 5 0.4505050 2.537086 No 2.537086
## 6 0.6925783 8.789090 Si NAPara este ejemplo, los gráficos que relacionan las variables (con y sin valores faltantes) son:
Ahora, con el 40% de valores faltantes, es necesario imputar los datos faltantes. Para esto, utilizaremos la técnica de imputación múltiple propuesta por Rubin (1987). La idea consiste en generar \(M > 1\) conjuntos de valores para los datos faltantes. Al final, el valor imputado corresponderá al promedio de esos \(M\) valores. Por tanto el modelo final de imputación (para los valores faltantes) toma la siguiente forma:
\[\dot{y}_i = \dot{\beta} x_{i - (missing)}+ \dot{\varepsilon_i}\]
El símbolo \(\dot{}\) representa la naturaleza estocástica de la imputación. Por ejemplo, aunque \(\dot{\beta}\) representa un valor diferente a \(\beta = 10\), es un valor que viene de una distribución de probabilidad. Hay varias maneras de realizar la imputación:
- Ingenua: \(\dot{y}_i = \beta x_{i - (missing)}\). Esta clase de imputación carece de aleatoriedad y por tanto, la varianza de \(\beta\) va a ser subestimada.
- Bootstrap: \(\dot{y}_i = \dot{\beta} x_{i - (missing)}+ \dot{\varepsilon_i}\). Se seleccionan \(m\) muestras bootstrap, y para cada una se estiman los parámetros \(\beta\) y \(\sigma\) para generar \(\dot{y}_i\). Al final se promedian los \(m\) valores y se imputa el valor faltante.
- Bayesiana: \(\dot{y}_i = \dot{\beta} x_{i - (missing)}+ \dot{\varepsilon_i}\). Se definen las distribuciones posteriores de \(\beta\) y \(\sigma\) para generar \(M\) valores de estos parámetros y por tanto \(M\) valores de \(\dot{y}_i\). Al final se promedian los \(M\) valores y se imputa el valor faltante.
El valor agregado de la imputación múltiple realmente está en la estimación de los errores estándar. No tener en cuenta la naturaleza estocástica de los valores imputados arroja estimaciones de la varianza mucho menores. Esto puede ser uno de los errores más graves que un estadístico puede cometer, puesto que afecta la cobertura nominal de los intervalos de confianza y a su vez influye en las pruebas de hipótesis y en el cálculo de los valores p.
Por ejemplo, si el interés es la estimación de la pendiente de la regresión 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\]
Por otro lado, 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\]
Imputación ingenua
Una función que realiza esta imputación es la siguiente:
im.ingenua <- function(datos){
Ind <- is.na(datos$y.per)
na.omit(datos)
beta1 <- lm(y ~ 0 + x, data = datos, na.action = na.omit)$coeff
datos$y.per[Ind] <- datos$x[Ind] * beta1
model.input <- lm(y.per ~ 0 + x, data = datos)
beta.input <- model.input$coeff
beta.sd <- summary(model.input)$coeff[2]
result <- list(new = datos, beta = beta.input, sd = beta.sd)
}Al aplicar la función sobre el conjunto de datos creado, se obtienen las siguientes salidas:
datos <- generar()
im.ingenua(datos)$beta## x
## 10.09468im.ingenua(datos)$sd## [1] 0.1187957head(im.ingenua(datos)$new)## x y faltantes y.per
## 1 0.7865446 11.301002 No 11.301002
## 2 0.6821044 4.740224 Si 6.867376
## 3 0.1734049 1.994628 No 1.994628
## 4 0.5161599 3.141765 Si 5.196659
## 5 0.6095868 2.103371 No 2.103371
## 6 0.7082076 6.043318 No 6.043318Nótese que todos los valores faltantes imputados caen, como era de esperarse, sobre la línea de regresión. Esta es una característica que implica una subestimación de la varianza de la variable respuesta.
Imputación Bootstrap
Una función que realiza esta imputación es la siguiente:
im.bootstrap <- function(datos, M = 15){
library(dplyr)
n <- nrow(datos)
datos1 <- na.omit(datos)
n1 <- nrow(datos1)
n0 <- n - n1
Ind <- is.na(datos$y.per)
faltantes.boot <- NULL
beta1 <- NULL
sigma1 <- NULL
for (m in 1:M){
datos.m <- dplyr::sample_n(datos1, n1, replace = TRUE)
model1 <- lm(y ~ 0 + x, data = datos.m)
beta <- model1$coeff
sigma <- sqrt(anova(model1)[["Mean Sq"]][2])
faltantes.boot <- rnorm(n0, datos$x[Ind] * beta, sd = sigma)
datos$y.per[Ind] <- faltantes.boot
model.input <- lm(y.per ~ 0 + x, data = datos)
beta1[m] <- model.input$coeff
sigma1[m] <- summary(model.input)$coeff[2]
}
beta.input <- mean(beta1)
u.bar <- mean(sigma1 ^ 2)
B <- var(beta1)
beta.sd <- sqrt(u.bar + B + B/M)
result <- list(new = datos, beta = beta.input, sd = beta.sd)
}Al aplicar la función sobre el conjunto de datos creado, se obtienen las siguientes salidas:
datos <- generar()
im.bootstrap(datos)$beta## [1] 10.25154im.bootstrap(datos)$sd## [1] 0.2345388head(im.bootstrap(datos)$new)## x y faltantes y.per
## 1 0.97794174 9.683813 Si 12.992767
## 2 0.44932770 5.823474 No 5.823474
## 3 0.66077397 5.651702 No 5.651702
## 4 0.01300767 -1.497382 No -1.497382
## 5 0.28542032 2.388563 No 2.388563
## 6 0.13981969 2.988707 No 2.988707Nótese que existe una buena dispersión en los valores imputados.
Imputación Bayesiana
Bajo distribuciones previas no informativas, la distribución posterior de \(\sigma^2\) es: \[\sigma^2| y, x \sim \frac{\sum_{i = 1}^{n_1} (y_i - \hat{\beta} x_i)^2}{\chi ^2_{n_1-1}}\]
con \(\hat{\beta} = \frac{\sum_{i = 1}^{n_1} x_i y_i}{\sum_{i = 1}^{n_1} x_i^2}\). La distribución posterior de \(\beta\) es: \[\beta | \sigma^2, y, x \sim Normal(\hat{\beta}, \frac{\sigma^2}{\sum_{i = 1}^{n_1} x_i^2} )\]
Una función que realiza esta imputación es la siguiente:
im.bayes <- function(datos, M = 15){
library(dplyr)
n <- nrow(datos)
datos1 <- na.omit(datos)
n1 <- nrow(datos1)
n0 <- n - n1
Ind <- is.na(datos$y.per)
model1 <- lm(y ~ 0 + x, data = datos1)
beta <- model1$coeff
sigma <- sqrt(anova(model1)[["Mean Sq"]][2])
faltantes.bayes <- NULL
beta1 <- NULL
sigma1 <- NULL
for (m in 1:M) {
g <- rchisq(1, n1 - 1)
sigma.pos <- (n1 - 1) * sigma ^ 2 / g
var.pos <- sigma.pos / sum(datos1$x ^ 2)
beta.pos <- rnorm(1, beta, sd = sqrt(var.pos))
faltantes.bayes <- rnorm(n0, datos$x[Ind] * beta.pos, sd = sqrt(sigma.pos))
datos$y.per[Ind] <- faltantes.bayes
model.input <- lm(y.per ~ 0 + x, data = datos)
beta1[m] <- model.input$coeff
sigma1[m] <- summary(model.input)$coeff[2]
}
beta.input <- mean(beta1)
u.bar <- mean(sigma1 ^ 2)
B <- var(beta1)
beta.sd <- sqrt(u.bar + B + B/M)
result <- list(new = datos, beta = beta.input, sd = beta.sd)
}Al aplicar la función sobre el conjunto de datos creado, se obtienen las siguientes salidas:
datos <- generar()
im.bayes(datos)$beta## [1] 9.97836im.bayes(datos)$sd## [1] 0.207771head(im.bayes(datos)$new)## x y faltantes y.per
## 1 0.4622356 7.177236 No 7.1772358
## 2 0.5092983 2.707868 No 2.7078680
## 3 0.4112130 2.346065 Si 4.5994068
## 4 0.5185065 5.894138 Si 7.4832260
## 5 0.3391639 2.691695 Si -0.2335784
## 6 0.7489785 4.703660 No 4.7036598Nótese que existe una buena dispersión en los valores imputados.
Simulación de Monte Carlo
La siguiente función agrupa todas las anteriores funciones de imputación y arroja resultados sobre la esperanza de las imputaciones, así como sus respectivos errores estándar.
simu <- function(){
im.betas <- NULL
im.sd <- NULL
datos <- generar()
im.betas[1] <- im.ingenua(datos)$beta
im.betas[2] <- im.bootstrap(datos)$beta
im.betas[3] <- im.bayes(datos)$beta
im.betas[4] <- im.ingenua(datos)$sd
im.betas[5] <- im.bootstrap(datos)$sd
im.betas[6] <- im.bayes(datos)$sd
return(im.betas)
}
resultados <- t(replicate(n = 1000, simu()))
resultados <- as.data.frame(resultados)
colnames(resultados) <- c("beta.ing", "beta.boot", "beta.bayes", "sd.ing", "sd.boot", "sd.bayes")
head(resultados)## beta.ing beta.boot beta.bayes sd.ing sd.boot sd.bayes
## 1 10.056772 10.115601 10.171105 0.1215380 0.1947824 0.1950999
## 2 10.160699 10.241731 10.133356 0.1355619 0.2299487 0.2122155
## 3 10.112279 10.048669 10.079439 0.1169640 0.2330632 0.1882248
## 4 9.994243 9.986130 9.992585 0.1192448 0.1726904 0.2034506
## 5 9.911859 9.937705 9.910568 0.1203331 0.2011541 0.2016907
## 6 9.816548 9.870926 9.836582 0.1224530 0.2258406 0.2037236resultados$long.ing <- (resultados[,1] + 1.96 * resultados[,4]) -
(resultados[,1] - 1.96 * resultados[,4])
resultados$long.boot <- (resultados[,2] + 1.96 * resultados[,5]) -
(resultados[,2] - 1.96 * resultados[,5])
resultados$long.bayes <- (resultados[,3] + 1.96 * resultados[,6]) -
(resultados[,3] - 1.96 * resultados[,6])
resultados$cob.ing <- (resultados[,1] - 1.96 * resultados[,4]) < 10 &
10 < (resultados[,1] + 1.96 * resultados[,4])
resultados$cob.boot <- (resultados[,2] - 1.96 * resultados[,5]) < 10 &
10 < (resultados[,2] + 1.96 * resultados[,5])
resultados$cob.bayes <- (resultados[,3] - 1.96 * resultados[,6]) < 10 &
10 < (resultados[,3] + 1.96 * resultados[,6])
colMeans(resultados)## beta.ing beta.boot beta.bayes sd.ing sd.boot sd.bayes
## 10.0023061 10.0025197 10.0031333 0.1200185 0.2015411 0.2018985
## long.ing long.boot long.bayes cob.ing cob.boot cob.bayes
## 0.4704725 0.7900411 0.7914423 0.8340000 0.9500000 0.9410000Nótese que inclusive la imputación ingenua arroja estimaciones insesgadas. Sin emabrgo, el error estándar de esta es el más bajo. Por lo tanto, la longitud de los intervalos que utilizan imputación múltiple es un poco mayor y como resultado la cobertura del intervalo de confianza ingenuo no corresponde a un nivel de confianza del 0.95.