MICE (Multivariate imputation by chained equations) es un método de imputación múltiple que se utiliza para reemplazar los valores de datos que faltan en un conjunto de datos bajo ciertos supuestos sobre el mecanismo de falta de datos, el planteamiento de ecuaciones encadenadas es muy flexible y puede manejar variables de distintos tipos (por ejemplo, continuas o binarias).
Funciona creando múltiples versiones de un conjunto de datos, donde cada versión contiene diferentes imputaciones (valores estimados) para los datos faltantes. Luego, se realiza el análisis estadístico en cada conjunto de datos imputado y los resultados se combinan para obtener una estimación más robusta y precisa, para cada variable con datos faltantes, se crea un modelo predictivo utilizando las otras variables del conjunto de datos como predictores.
Para el presente proyecto se realizara este metodo para la base de datos Framinghan, la cual corresponde a una base pequeña del Framinghan Heart Study.
fram<- read.csv("framingham_dataset.csv", header = TRUE)
#Para poder imputar los datos tomaremos la libreria "mice" de R
library(mice)
#Observamos los datos faltantes de nuestra base
library(VIM)
aggr(fram)
Hacemos la imputacion para toda la base, la desventaja de este metodo es que genera una base de datos nueva con los datos imputados y se genera una nueva por cada iteracion ocupando recursos computacionales, sin embargo, MICE permite trabajar con conjuntos de datos que tienen valores faltantes, lo que es común en muchos conjuntos de datos del mundo real.
fram_imp <- mice(fram) #Por defecto itera 5 veces y la documentacion recomienda de 5 a 10, en este caso lo dejaremos en 5
##
## iter imp variable
## 1 1 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 1 2 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 1 3 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 1 4 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 1 5 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 2 1 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 2 2 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 2 3 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 2 4 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 2 5 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 3 1 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 3 2 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 3 3 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 3 4 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 3 5 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 4 1 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 4 2 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 4 3 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 4 4 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 4 5 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 5 1 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 5 2 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 5 3 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 5 4 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
## 5 5 totchol cigpday bmi bpmeds heartrte glucose educ hdlc ldlc
#completamos la data
fram_comp <- complete(fram_imp)
aggr(fram_comp)
Usuaremos las mismas variables del ejercicio anterior (Tabla 2 Datos Framingham) pero ahora con la base de datos imputados:
#definimos las variables:
age <- fram_comp$age
sex <- factor(fram_comp$sex,levels = c(1,2) , labels = c("Male","Female"))
glucose <- fram_comp$glucose #mg/dL
diabetes <- factor(fram_comp$diabetes,levels = c(0,1) , labels = c("no","si"))
stroke <- factor(fram_comp$stroke,levels = c(0,1) , labels = c("no","si"))
cigarettes <- fram_comp$cigpday #cigarettes per day
hyperten <- factor(fram_comp$hyperten,levels = c(0,1) , labels = c("no","si"))
cholesterol <- fram_comp$totchol #Total cholesterol (mg,dL)
bmi <- fram_comp$bmi #Body mass index
rate <- fram_comp$heartrte # Heart rate (beats per minute)
Para ver la diferencia entre los datos tomaremos la variable de glucosa la cual tiene una de las variables con mas datos faltantes y la compararemos con la base de datos nueva
plot(fram$age,fram$glucose,main = "Scatter plot: Age vs Glucose whit NA", xlab = "Age (Years)", ylab = "Glucose (mg/dL)")
plot(age,glucose,main = "Scatter plot: Age vs Cholesterol data complete", xlab = "Age (Years)", ylab = "Glucose (mg/dL)")
Observamos una grafica mas densa en los datos imputados que en la base original, ahora proseguimos con la tabla 2 con datos imputados
library(tidyverse)
library(broom)
fchd <- factor(fram_comp$mi_fchd , levels = c(0,1) , labels = c("no","si"))
modelo <- glm(fchd ~ age + sex + glucose + diabetes + hyperten + cigarettes + stroke + rate + cholesterol, data = fram, family = binomial(link = "logit"))
summary(modelo)
##
## Call:
## glm(formula = fchd ~ age + sex + glucose + diabetes + hyperten +
## cigarettes + stroke + rate + cholesterol, family = binomial(link = "logit"),
## data = fram)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3098951 0.2990460 -14.412 < 2e-16 ***
## age 0.0273002 0.0033062 8.257 < 2e-16 ***
## sex -1.2811627 0.0639924 -20.021 < 2e-16 ***
## glucose 0.0032734 0.0011190 2.925 0.00344 **
## diabetes 0.7956830 0.1305418 6.095 1.09e-09 ***
## hyperten 0.6285881 0.0787245 7.985 1.41e-15 ***
## cigarettes 0.0092015 0.0023181 3.969 7.21e-05 ***
## stroke 0.3582672 0.0879392 4.074 4.62e-05 ***
## rate 0.0028080 0.0023516 1.194 0.23244
## cholesterol 0.0072562 0.0006333 11.458 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8742.8 on 10186 degrees of freedom
## Residual deviance: 7841.7 on 10177 degrees of freedom
## (1440 observations deleted due to missingness)
## AIC: 7861.7
##
## Number of Fisher Scoring iterations: 5
Viendo la comparacion de ambas tablas observamos que cambian los coeficientes y disminuyen los errores estandar, dandonos asi un nuevo modelo de prediccion probabilistico dandole a este mas precision
tabla_2 <- tidy(modelo, conf.int = TRUE, exponentiate = TRUE) %>%
mutate(across(c(estimate, conf.low, conf.high), ~round(., 2))) %>%
mutate(p.value = ifelse(p.value < 0.001, "<0.001", round(p.value, 3))) %>%
select(term, estimate, std.error, p.value, conf.low, conf.high)
knitr::kable(tabla_2,
col.names = c("Variable", "OR", "Error Estándar", "p-valor", "IC 95% (bajo)", "IC 95% (alto)"),
caption = " Tabla 2: Modelo de regresión logística para predecir enfermedad de miocardio")
Variable | OR | Error Estándar | p-valor | IC 95% (bajo) | IC 95% (alto) |
---|---|---|---|---|---|
(Intercept) | 0.01 | 0.2990460 | <0.001 | 0.01 | 0.02 |
age | 1.03 | 0.0033062 | <0.001 | 1.02 | 1.03 |
sex | 0.28 | 0.0639924 | <0.001 | 0.24 | 0.31 |
glucose | 1.00 | 0.0011190 | 0.003 | 1.00 | 1.01 |
diabetes | 2.22 | 0.1305418 | <0.001 | 1.71 | 2.86 |
hyperten | 1.87 | 0.0787245 | <0.001 | 1.61 | 2.19 |
cigarettes | 1.01 | 0.0023181 | <0.001 | 1.00 | 1.01 |
stroke | 1.43 | 0.0879392 | <0.001 | 1.20 | 1.70 |
rate | 1.00 | 0.0023516 | 0.232 | 1.00 | 1.01 |
cholesterol | 1.01 | 0.0006333 | <0.001 | 1.01 | 1.01 |
Las OR de los datos originales comparados con los datos imputados no difieren, solo en la variable hyperten que aumenta 3 milesimas, por lo tanto el nuevo modelo probabilistico quedaria de la siguiente forma:
\[ logit(L) = -4.312 + 0.027*age -1.281*sex + 0.003*glucose + 0.792*diabetes + 0.682*hyperten + 0.009*cigpday\]
\[ + 0.360*stroke + 0.002*rate + 0.007*cholesterol \]