Imputation multiple pour les valeurs manquantes

Je voudrais utiliser l'imputation pour remplacer les valeurs manquantes dans mon ensemble de données sous certaines contraintes.

Par exemple, j'aimerais que la variable imputée x1soit supérieure ou égale à la somme de mes deux autres variables, disons x2et x3. Je veux aussi x3être imputé par ou 0ou >= 14et je veux x2être imputé par 0ou >= 16.

J'ai essayé de définir ces contraintes dans SPSS pour l'imputation multiple, mais dans SPSS je ne peux définir que des valeurs maximales et minimales. Existe-t-il un moyen de définir d'autres contraintes dans SPSS ou connaissez-vous un package R qui me permettrait de définir de telles contraintes pour l'imputation des valeurs manquantes?

Mes données sont les suivantes:

   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

Une solution consiste à écrire vos propres fonctions d'imputation personnalisées pour le micepackage. Le package est préparé pour cela et la configuration est étonnamment indolore.

D'abord, nous configurons les données comme suggéré:

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))

Ensuite, nous chargeons le micepackage et voyons quelles méthodes il choisit par défaut:

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

Le pmmreprésente l' appariement moyen prédictif - probablement l'algorithme d'imputation le plus populaire pour imputer des variables continues. Il calcule la valeur prédite à l'aide d'un modèle de régression et sélectionne les 5 éléments les plus proches de la valeur prédite (par distance euclidienne ). Ces éléments choisis sont appelés le pool de donateurs et la valeur finale est choisie au hasard dans ce pool de donateurs.

À partir de la matrice de prédiction, nous constatons que les méthodes obtiennent les variables qui sont intéressantes pour les restrictions. Notez que la ligne est la variable cible et la colonne les prédicteurs. Si x1 n'avait pas 1 dans la colonne x3, nous devrions ajouter ceci dans la matrice:imp_base$predictorMatrix["x1","x3"] <- 1

Passons maintenant à la partie amusante, générant les méthodes d'imputation. J'ai choisi une méthode plutôt grossière ici où je rejette toutes les valeurs si elles ne répondent pas aux critères. Cela peut entraîner un long temps de boucle et il peut être potentiellement plus efficace de conserver les imputations valides et de ne refaire que les autres, cela nécessiterait cependant un peu plus de réglages.

# 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)
  }

Une fois que nous avons fini de définir les méthodes, nous changeons simplement les méthodes précédentes. Si vous ne souhaitez modifier qu'une seule variable, vous pouvez simplement l'utiliser, imp_base$method["x2"] <- "pmm_x2"mais pour cet exemple, nous allons tout changer (la dénomination n'est pas nécessaire):

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)

Voyons maintenant le troisième ensemble de données imputé:

> 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, ça fait l'affaire. J'aime cette solution car vous pouvez vous superposer aux fonctions principales et ajouter simplement les restrictions que vous trouvez significatives.

Mise à jour

Afin d'appliquer les restrictions rigoureuses @ t0x1n mentionnées dans les commentaires, nous pouvons vouloir ajouter les capacités suivantes à la fonction wrapper:

1. Enregistrer les valeurs valides pendant les boucles afin que les données des exécutions précédentes partiellement réussies ne soient pas supprimées
2. Un mécanisme d'échappement pour éviter les boucles infinies
3. Gonflez le pool de donneurs après avoir essayé x fois sans trouver une correspondance appropriée (cela s'applique principalement à pmm)

Il en résulte une fonction wrapper légèrement plus compliquée:

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)
}

Notez que cela ne fonctionne pas très bien, probablement en raison du fait que l'ensemble de données suggéré ne respecte pas les contraintes pour tous les cas sans manquer. J'ai besoin d'augmenter la longueur de la boucle à 400-500 avant même qu'elle ne commence à se comporter. Je suppose que cela n'est pas intentionnel, votre imputation devrait imiter la façon dont les données réelles sont générées.

Optimisation

L'argument rycontient les valeurs non manquantes et nous pourrions éventuellement accélérer la boucle en supprimant les éléments que nous avons trouvés des imputations éligibles, mais comme je ne connais pas les fonctions internes, je m'en suis abstenu.

Je pense que la chose la plus importante lorsque vous avez de fortes contraintes qui prennent du temps à remplir est de paralléliser vos imputations ( voir ma réponse sur CrossValidated ). La plupart ont aujourd'hui des ordinateurs avec 4 à 8 cœurs et R n'en utilise qu'un par défaut. Le temps peut être (presque) divisé par deux en doublant le nombre de cœurs.

Paramètres manquants à l'imputation

Concernant le problème d' x2être manquant au moment de l'imputation - les souris n'alimentent jamais les valeurs manquantes dans le x- data.frame. La méthode des souris comprend le remplissage d'une valeur aléatoire au début. La partie chaîne de l'imputation limite l'impact de cette valeur initiale. Si vous regardez la micefonction-vous pouvez le trouver avant l'appel d'imputation (la mice:::samplerfonction-):

...
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))
  }
}
...

Le data.initpeut être fourni à la micefonction et le mice.imput.sample est une procédure d'échantillonnage de base.

Séquence de visites

Si la séquence de visites est importante, vous pouvez spécifier l'ordre dans lequel la micefonction-exécute les imputations. La valeur par défaut est de 1:ncol(data)mais vous pouvez définir visitSequencece que vous voulez.

La chose la plus proche que j'ai pu trouver est l' inclusion d'informations préalables d'Amelia . Voir le chapitre 4.7 de la vignette , en particulier 4.7.2:

Prieurs de niveau observation

Les chercheurs ont souvent des informations préalables supplémentaires sur les valeurs de données manquantes basées sur des recherches antérieures, un consensus académique ou une expérience personnelle. Amelia peut intégrer ces informations pour produire des imputations considérablement améliorées. L'algorithme Amelia permet aux utilisateurs d'inclure des a priori bayésiens informatifs sur les cellules de données manquantes individuelles au lieu des paramètres de modèle plus généraux, dont beaucoup ont peu de signification directe.

L'incorporation des a priori suit une analyse bayésienne de base où l'imputation se révèle être une moyenne pondérée de l'imputation basée sur le modèle et la moyenne antérieure, où les poids sont fonction de la force relative des données et a priori: lorsque le modèle prédit très bien , l'imputation réduira la pondération de l'a priori et vice versa (Honaker et King, 2010).

Les priorités des observations individuelles doivent décrire la croyance de l'analyste à propos de la distribution de la cellule de données manquante. Cela peut prendre la forme d'une moyenne et d'un écart-type ou d'un intervalle de condence. Par exemple, nous pouvons savoir que les taux de tari de 1986 en Thaïlande sont d'environ 40%, mais nous avons une certaine incertitude quant à la valeur exacte. Notre croyance antérieure quant à la distribution de la cellule de données manquante se concentre donc sur 40 avec un écart-type qui reflète la quantité d'incertitude que nous avons sur notre croyance antérieure.

Pour saisir des priorités, vous devez créer une matrice de priorités avec quatre ou cinq colonnes. Chaque ligne de la matrice représente un a priori sur une observation ou une variable. Dans n'importe quelle ligne, l'entrée dans la première colonne est la ligne de l'observation et l'entrée est la deuxième colonne est la colonne de l'observation. Dans la matrice a priori à quatre colonnes, les troisième et quatrième colonnes sont la moyenne et l'écart-type de la distribution antérieure de la valeur manquante.

Ainsi, bien que vous ne puissiez généralement pas dire quelque chose comme cela x1<x2+x3, vous pouvez parcourir votre ensemble de données et ajouter un niveau d'observation avant pour chaque cas pertinent. Des limites constantes peuvent également être appliquées (par exemple, en définissant x1, x2 et x3 pour qu'elles soient non négatives). Par exemple:

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);

Les contraintes sont probablement plus faciles à mettre en œuvre dans l'imputation multiple à appariement moyen prédictif. Cela suppose qu'il existe un nombre significatif d'observations avec des variables contraignantes non manquantes qui satisfont aux contraintes. Je pense à implémenter cela dans la fonction de Hmiscpackage R. aregImputeVous voudrez peut-être revenir dans un mois environ. Il sera important de spécifier la distance maximale de la cible que peut être une observation de donneur, car les contraintes pousseront les donneurs plus loin du donneur sans contrainte idéal.

Je pense que le Ameliapackage (Amelia II) a actuellement le support le plus complet pour spécifier les contraintes de plage de valeurs de données. Cependant, le problème est que cela Ameliasuppose que les données sont normales à plusieurs variables.

Si dans votre cas, l'hypothèse de normalité multivariée ne s'applique pas, vous pouvez vérifier le micepackage, qui implémente l' imputation multiple (MI) via des équations chaînées . Ce package n'a pas l'hypothèse de normalité multivariée . Il a également une fonction qui pourrait être suffisante pour spécifier des contraintes , mais je ne sais pas dans quelle mesure. La fonction est appelée squeeze(). Vous pouvez en lire plus dans la documentation: http://cran.r-project.org/web/packages/mice/mice.pdf . Un avantage supplémentaire de micesa flexibilité en termes de permettre la spécification de fonctions d'imputation définies par l'utilisateur et une sélection plus large d'algorithmes. Voici un tutoriel sur l'exécution de MI, en utilisant mice:http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm .

Autant que je sache, le Hmiscprogiciel du Dr Harrell , utilisant la même approche d' équations chaînées ( correspondance moyenne prédictive ), prend probablement en charge les données non normales (à l'exception de la normpmmméthode). Peut-être qu'il a déjà implémenté la fonctionnalité de spécification des contraintes selon sa réponse ci-dessus. Je n'ai pas utilisé aregImpute(), donc je ne peux pas en dire plus à ce sujet (j'ai utilisé Ameliaet mice, mais je ne suis certainement pas un expert en statistiques, j'essaie simplement d'apprendre autant que possible).

Enfin, vous trouverez peut-être intéressant ce qui suit, un peu daté, mais toujours agréable, un aperçu des approches, méthodes et logiciels pour l'imputation multiple de données avec des valeurs manquantes: http://www.ncbi.nlm.nih.gov/pmc/articles / PMC1839993 . Je suis sûr qu'il existe de nouveaux articles sur l'IM, mais c'est tout ce que je sais à l'heure actuelle. J'espère que cela est quelque peu utile.

Si je comprends bien votre question, il me semble que vous savez déjà quelles valeurs les variables manquantes doivent prendre sous réserve de certaines contraintes. Je ne suis pas très familier avec SPSS mais en RI je pense que vous pouvez écrire une fonction pour faire cela (ce qui ne devrait pas être trop difficile en fonction de votre expérience je devrais dire). Je ne connais aucun paquet qui fonctionne avec de telles contraintes.