Imputación múltiple para valores perdidos

Me gustaría usar la imputación para reemplazar los valores faltantes en mi conjunto de datos bajo ciertas restricciones.

Por ejemplo, me gustaría que la variable imputada x1sea ​​mayor o igual a la suma de mis otras dos variables, digamos x2y x3. También quiero x3ser imputado por cualquiera 0o >= 14y quiero x2ser imputado por cualquiera 0o >= 16.

Intenté definir estas restricciones en SPSS para la imputación múltiple, pero en SPSS solo puedo definir valores máximos y mínimos. ¿Hay alguna forma de definir más restricciones en SPSS o conoces algún paquete R que me permita definir tales restricciones para la imputación de valores perdidos?

Mis datos son los siguientes:

   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)
   dat=data.frame(x1=x1, x2=x2, x3=x3)
   > dat
       x1 x2 x3
   1   21  0  0
   2   50 NA  0
   3   31 18  0
   4   15  0  0
   5   36 19  0
   6   82  0 54
   7   14 NA  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  NA NA NA
   15  NA NA NA
   16  NA NA NA
   17  50 22 NA
   18  52 NA  0
   19  26  0  0
   20  24  0  0

Una solución es escribir sus propias funciones de imputación personalizadas para el micepaquete. El paquete está preparado para esto y la configuración sorprendentemente sin dolor.

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 micepaquete y vemos qué métodos elige de forma predeterminada:

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 pmmrepresenta predictivo coincidencia media - probablemente el algoritmo de imputación más populares para imputar las variables continuas. Calcula el valor predicho utilizando un modelo de regresión y selecciona los 5 elementos más cercanos al valor predicho (por distancia euclidiana ). Estos elementos elegidos se denominan grupo de donantes y el valor final se elige al azar de este grupo de donantes.

De la matriz de predicción encontramos que los métodos obtienen las variables que son interesantes para las restricciones. Tenga en cuenta que la fila es la variable objetivo y la columna los predictores. Si x1 no tuviera 1 en la columna x3, tendríamos que agregar esto en la matriz:imp_base$predictorMatrix["x1","x3"] <- 1

Ahora a la parte divertida, generando los métodos de imputación. Elegí un método bastante crudo aquí donde descarto todos los valores si no cumplen con los criterios. Esto puede resultar en un tiempo de bucle prolongado y puede ser potencialmente más eficiente mantener las imputaciones válidas y solo rehacer las restantes, aunque requeriría un poco más de ajustes.

# 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 hayamos terminado de definir los métodos, simplemente cambiaremos los métodos anteriores. Si solo desea cambiar una sola variable, simplemente puede usarla, imp_base$method["x2"] <- "pmm_x2"pero para este ejemplo cambiaremos todas (la nomenclatura 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)

Ahora echemos un vistazo al tercer conjunto de datos imputado:

> 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

Ok, eso hace el trabajo. Me gusta esta solución, ya que puedes aprovechar las funciones principales y solo agregar las restricciones que consideres significativas.

Actualizar

Para hacer cumplir las rigurosas restricciones @ t0x1n mencionadas en los comentarios, es posible que deseemos agregar las siguientes habilidades a la función de envoltura:

1. Guarde valores válidos durante los bucles para que no se descarten los datos de ejecuciones parcialmente exitosas anteriores
2. Un mecanismo de escape para evitar bucles infinitos.
3. Infle el grupo de donantes después de intentar x veces sin encontrar una coincidencia adecuada (esto se aplica principalmente a pmm)

Esto da como resultado una función de envoltura 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)
}

Tenga en cuenta que esto no funciona tan bien, probablemente debido a que el conjunto de datos sugerido falla sin restricciones en todos los casos. Necesito aumentar la longitud del bucle a 400-500 incluso antes de que comience a comportarse. Supongo que esto no es intencional, su imputación debe imitar cómo se generan los datos reales.

Mejoramiento

El argumento rycontiene los valores que no faltan y posiblemente podríamos acelerar el ciclo eliminando los elementos que hemos encontrado imputaciones elegibles, pero como no estoy familiarizado con las funciones internas, me he abstenido de esto.

Creo que lo más importante cuando tienes fuertes restricciones que tardan en completarse es paralelizar tus imputaciones ( consulta mi respuesta en CrossValidated ). La mayoría tiene hoy computadoras con 4-8 núcleos y R solo usa una de ellas de manera predeterminada. El tiempo puede dividirse (casi) por la mitad duplicando el número de núcleos.

Faltan parámetros en la imputación

Con respecto al problema de x2faltar en el momento de la imputación, los ratones nunca introducen valores perdidos en el x- data.frame. El método de los ratones incluye completar algún valor aleatorio al inicio. La parte de la cadena de la imputación limita el impacto de este valor inicial. Si observa la micefunción, puede encontrarla antes de la llamada de imputación (la mice:::samplerfunció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.initpuede suministrarse a la micefunción y mouse.imput.sample es un procedimiento de muestreo básico.

Secuencia de visita

Si la secuencia de visitas es importante, puede especificar el orden en que la función mice-funciona las imputaciones. El valor predeterminado es from 1:ncol(data)pero puedes configurarlo visitSequencepara que sea lo que quieras.

Lo más cercano que pude encontrar es la inclusión de información previa de Amelia . Vea el capítulo 4.7 en la viñeta , específicamente 4.7.2:

Nivel previo de observación

Los investigadores a menudo tienen información previa adicional sobre valores de datos faltantes basados ​​en investigaciones previas, consenso académico o experiencia personal. Amelia puede incorporar esta información para producir imputaciones muy mejoradas. El algoritmo Amelia permite a los usuarios incluir información previa bayesiana informativa sobre celdas de datos faltantes individuales en lugar de los parámetros más generales del modelo, muchos de los cuales tienen poco significado directo.

La incorporación de anteriores sigue el análisis bayesiano básico donde la imputación resulta ser un promedio ponderado de la imputación basada en el modelo y la media anterior, donde los pesos son funciones de la fuerza relativa de los datos y anteriores: cuando el modelo predice muy bien , la imputación bajará el peso del anterior y viceversa (Honaker y King, 2010).

Los antecedentes sobre observaciones individuales deben describir la creencia del analista sobre la distribución de la celda de datos que falta. Esto puede tomar la forma de una media y una desviación estándar o un intervalo de condensación. Por ejemplo, podríamos saber que las tasas de tari de 1986 en Tailandia rondan el 40%, pero tenemos cierta incertidumbre sobre el valor exacto. Nuestra creencia previa sobre la distribución de la celda de datos faltantes, entonces, se centra en 40 con una desviación estándar que refleja la cantidad de incertidumbre que tenemos sobre nuestra creencia previa.

Para ingresar priors debe construir una matriz de priors con cuatro o cinco columnas. Cada fila de la matriz representa un previo en una observación o en una variable. En cualquier fila, la entrada en la primera columna es la fila de la observación y la entrada es la segunda columna es la columna de la observación. En la matriz de cuatro columnas anteriores, las columnas tercera y cuarta son la media y la desviación estándar de la distribución previa del valor faltante.

Entonces, aunque generalmente no podrá decir algo como x1<x2+x3, podría recorrer su conjunto de datos y agregar un nivel de observación anterior para cada caso relevante. También se pueden aplicar límites constantes (como configurar x1, x2 y x3 para que no sean negativos). Por ejemplo:

priors = matrix(NA, nrow=0, ncol=5);
for (i in seq(1, length(data))) 
{
    x1 = data$x1[i];
    x2 = data$x2[i];
    x3 = data$x3[i];

    if (is.na(x1) && !is.na(x2) && !is.na(x3))
    {
        priors = rbind(priors, c(i, 1, 0, x2+x3, 0.999999))
    }
}

amelia(data, m=1, bound = rbind(c(1, 0, Inf), c(2, 0, Inf), c(3, 0, Inf)), pr = priors);

Las restricciones son probablemente más fáciles de implementar en la media predictiva que coincide con la imputación múltiple. Esto supone que hay un número significativo de observaciones con variables de restricción no faltantes que cumplen con las restricciones. Estoy pensando en implementar esto en la función del Hmiscpaquete R. aregImputeEs posible que desee volver en un mes más o menos. Será importante especificar la distancia máxima desde el objetivo que puede ser la observación de un donante, porque las restricciones empujarán a los donantes más lejos del donante ideal sin restricciones.

Creo que el Ameliapaquete (Amelia II) actualmente tiene el soporte más completo para especificar restricciones de rango de valores de datos. Sin embargo, el problema es que Ameliasupone que los datos son multivariados normales.

Si en su caso no se aplica el supuesto de normalidad multivariante, es posible que desee verificar el micepaquete, que implementa la imputación múltiple (MI) a través de ecuaciones encadenadas . Este paquete no asume la normalidad multivariante . También tiene una función que podría ser suficiente para especificar restricciones , pero no estoy seguro de hasta qué punto. La función se llama squeeze(). Puede leer sobre esto en la documentación: http://cran.r-project.org/web/packages/mice/mice.pdf . Un beneficio adicional micees su flexibilidad en términos de permitir la especificación de funciones de imputación definidas por el usuario y una selección más amplia de algoritmos. Aquí hay un tutorial sobre cómo realizar MI, usando mice:http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm .

Según tengo entendido, el Hmiscpaquete del Dr. Harrell , usando el mismo enfoque de ecuaciones encadenadas ( coincidencia de predicción media ), probablemente admite datos no normales (con la excepción del normpmmmétodo). Tal vez ya haya implementado la funcionalidad de especificación de restricciones según su respuesta anterior. No lo he usado aregImpute(), así que no puedo decir mucho más al respecto (lo he usado Ameliay mice, pero definitivamente no soy un experto en estadística, solo trato de aprender todo lo que puedo).

Finalmente, puede encontrar interesante lo siguiente, un poco anticuado, pero aún así agradable, descripción general de enfoques, métodos y software para la imputación múltiple de datos con valores faltantes: http://www.ncbi.nlm.nih.gov/pmc/articles / PMC1839993 . Estoy seguro de que hay documentos generales más recientes sobre MI, pero eso es todo lo que sé en este momento. Espero que esto sea de alguna ayuda.

Si entiendo su pregunta correctamente, me parece que ya sabe qué valores deben tomar las variables que faltan, sujetas a algunas restricciones. No estoy muy familiarizado con SPSS, pero en RI creo que puede escribir una función para hacer eso (lo que no debería ser demasiado difícil dependiendo de su experiencia, debería decir). No conozco ningún paquete que funcione con tales restricciones.