Una solución es escribir sus propias funciones de imputación personalizadas para el mice paquete. El paquete está preparado para ello y la configuración es sorprendentemente sencilla.
Primero configuramos los datos como se sugiere:
dat=data.frame(x1=c(21, 50, 31, 15, 36, 82, 14, 14, 19, 18, 16, 36, 583, NA,NA,NA, 50, 52, 26, 24),
x2=c(0, NA, 18,0, 19, 0, NA, 0, 0, 0, 0, 0, 0,NA,NA, NA, 22, NA, 0, 0),
x3=c(0, 0, 0, 0, 0, 54, 0 ,0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0))
A continuación, cargamos el mice y ver qué métodos elige por defecto:
library(mice)
# Do a non-imputation
imp_base <- mice(dat, m=0, maxit = 0)
# Find the methods that mice chooses
imp_base$method
# Returns: "pmm" "pmm" "pmm"
# Look at the imputation matrix
imp_base$predictorMatrix
# Returns:
# x1 x2 x3
#x1 0 1 1
#x2 1 0 1
#x3 1 1 0
El pmm significa coincidencia de medias predictivas - probablemente el algoritmo de imputación más popular para imputar variables continuas. Calcula el valor predicho utilizando un modelo de regresión y elige los 5 elementos más cercanos al valor predicho (por Distancia euclidiana ). Estos elementos elegidos se denominan pool de donantes y el valor final se elige al azar de este pool de donantes.
A partir de la matriz de predicción encontramos que los métodos obtienen las variables pasadas que son de interés para las restricciones. Nótese que la fila es la variable objetivo y la columna los predictores. Si x1 no tuviera 1 en la columna x3 tendríamos que añadirlo en la matriz: imp_base$predictorMatrix["x1","x3"] <- 1
Ahora pasamos a la parte divertida, la generación de los métodos de imputación. Aquí he elegido un método bastante burdo en el que descarto todos los valores si no cumplen los criterios. Esto puede dar lugar a un largo tiempo de bucle y puede ser potencialmente más eficiente para mantener las imputaciones válidas y sólo rehacer los restantes, que requeriría un poco más de ajuste sin embargo.
# Generate our custom methods
mice.impute.pmm_x1 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
max_sum <- sum(max(x[,"x2"], na.rm=TRUE),
max(x[,"x3"], na.rm=TRUE))
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals < max_sum)){
break
}
}
return(vals)
}
mice.impute.pmm_x2 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals == 0 | vals >= 14)){
break
}
}
return(vals)
}
mice.impute.pmm_x3 <-
function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "",
...)
{
repeat{
vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...)
if (all(vals == 0 | vals >= 16)){
break
}
}
return(vals)
}
Una vez que hemos terminado de definir los métodos, simplemente cambiamos los métodos anteriores. Si sólo quieres cambiar una sola variable entonces puedes usar simplemente imp_base$method["x2"] <- "pmm_x2" pero para este ejemplo lo cambiaremos todo (la denominación no es necesaria):
imp_base$method <- c(x1 = "pmm_x1", x2 = "pmm_x2", x3 = "pmm_x3")
# The predictor matrix is not really necessary for this example
# but I use it just to illustrate in case you would like to
# modify it
imp_ds <-
mice(dat,
method = imp_base$method,
predictorMatrix = imp_base$predictorMatrix)
Veamos ahora el tercer conjunto de datos imputados:
> complete(imp_ds, action = 3)
x1 x2 x3
1 21 0 0
2 50 19 0
3 31 18 0
4 15 0 0
5 36 19 0
6 82 0 54
7 14 0 0
8 14 0 0
9 19 0 0
10 18 0 0
11 16 0 0
12 36 0 0
13 583 0 0
14 50 22 0
15 52 19 0
16 14 0 0
17 50 22 0
18 52 0 0
19 26 0 0
20 24 0 0
Bien, eso hace el trabajo. Me gusta esta solución, ya que se puede añadir a las funciones principales y añadir las restricciones que se consideren significativas.
Actualización
Para hacer cumplir las rigurosas restricciones que @t0x1n mencionó en los comentarios, podríamos añadir las siguientes capacidades a la función envolvente:
- Guardar los valores válidos durante los bucles para que no se descarten los datos de las ejecuciones anteriores, parcialmente exitosas
- Un mecanismo de escape para evitar los bucles infinitos
- Inflar la reserva de donantes después de intentar x veces sin encontrar una coincidencia adecuada (esto se aplica principalmente a pmm)
Esto da lugar a una función envolvente ligeramente más complicada:
mice.impute.pmm_x1_adv <- function (y, ry,
x, donors = 5,
type = 1, ridge = 1e-05,
version = "", ...) {
# The mice:::remove.lindep may remove the parts required for
# the test - in those cases we should escape the test
if (!all(c("x2", "x3") %in% colnames(x))){
warning("Could not enforce pmm_x1 due to missing column(s):",
c("x2", "x3")[!c("x2", "x3") %in% colnames(x)])
return(mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
version = "", ...))
}
# Select those missing
max_vals <- rowSums(x[!ry, c("x2", "x3")])
# We will keep saving the valid values in the valid_vals
valid_vals <- rep(NA, length.out = sum(!ry))
# We need a counter in order to avoid an eternal loop
# and for inflating the donor pool if no match is found
cntr <- 0
repeat{
# We should be prepared to increase the donor pool, otherwise
# the criteria may become imposs
donor_inflation <- floor(cntr/10)
vals <- mice.impute.pmm(y, ry, x,
donors = min(5 + donor_inflation, sum(ry)),
type = 1, ridge = 1e-05,
version = "", ...)
# Our criteria check
correct <- vals < max_vals
if (all(!is.na(valid_vals) |
correct)){
valid_vals[correct] <-
vals[correct]
break
}else if (any(is.na(valid_vals) &
correct)){
# Save the new valid values
valid_vals[correct] <-
vals[correct]
}
# An emergency exit to avoid endless loop
cntr <- cntr + 1
if (cntr > 200){
warning("Could not completely enforce constraints for ",
sum(is.na(valid_vals)),
" out of ",
length(valid_vals),
" missing elements")
if (all(is.na(valid_vals))){
valid_vals <- vals
}else{
valid_vals[is.na(valid_vals)] <-
vals[is.na(valid_vals)]
}
break
}
}
return(valid_vals)
}
Obsérvese que esto no funciona tan bien, probablemente debido a que el conjunto de datos sugerido no cumple con las restricciones para todo casos sin faltar. Tengo que aumentar la longitud del bucle a 400-500 antes de que empiece a comportarse. Supongo que esto es involuntario, su imputación debe imitar cómo se generan los datos reales.
Optimización
El argumento ry contiene los valores no ausentes y posiblemente podríamos acelerar el bucle eliminando los elementos que hemos encontrado imputables, pero como no estoy familiarizado con las funciones internas me he abstenido de hacerlo.
Creo que lo más importante cuando se tienen fuertes restricciones que tardan en completarse es paralelizar las imputaciones ( ver mi respuesta en CrossValidated ). La mayoría tiene hoy en día ordenadores con 4-8 núcleos y R sólo utiliza uno de ellos por defecto. El tiempo puede ser (casi) cortado a la mitad duplicando el número de núcleos.
Parámetros ausentes en la imputación
En cuanto al problema de x2 que faltan en el momento de la imputación - los ratones en realidad nunca introducen los valores que faltan en el x - data.frame . El ratones incluye la introducción de un valor aleatorio al inicio. La parte de la imputación en cadena limita el impacto de este valor inicial. Si se observa el mice -se puede encontrar esto antes de la llamada de imputación (el mice:::sampler -función):
...
if (method[j] != "") {
for (i in 1:m) {
if (nmis[j] < nrow(data)) {
if (is.null(data.init)) {
imp[[j]][, i] <- mice.impute.sample(y,
ry, ...)
}
else {
imp[[j]][, i] <- data.init[!ry, j]
}
}
else imp[[j]][, i] <- rnorm(nrow(data))
}
}
...
El data.init puede ser suministrado al mice y la función mice.imput.sample es un procedimiento básico de muestreo.
Secuencia de visitas
Si la secuencia de visitas es importante, puede especificar el orden en el que los mice -función ejecuta las imputaciones. Por defecto es de 1:ncol(data) pero se puede establecer el visitSequence para ser lo que quieras.
