Wielokrotna imputacja brakujących wartości

Chciałbym użyć imputacji do zastąpienia brakujących wartości w moim zbiorze danych z pewnymi ograniczeniami.

Na przykład chciałbym, aby zmienna przypisana x1była większa lub równa sumie moich dwóch innych zmiennych, powiedzmy x2i x3. Chcę też x3zostać przypisany przez jeden 0lub >= 14i chcę x2zostać przypisany przez jeden 0lub >= 16.

Próbowałem zdefiniować te ograniczenia w SPSS dla wielokrotnego imputacji, ale w SPSS mogę zdefiniować tylko wartości maksymalne i minimalne. Czy jest jakiś sposób na zdefiniowanie dalszych ograniczeń w SPSS lub znasz jakiś pakiet R, który pozwoliłby mi zdefiniować takie ograniczenia dla przypisania brakujących wartości?

Moje dane są następujące:

   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

Jednym z rozwiązań jest napisanie własnych niestandardowych funkcji imputacji dla micepakietu. Pakiet jest na to przygotowany, a konfiguracja zaskakująco bezbolesna.

Najpierw konfigurujemy dane zgodnie z sugestią:

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

Następnie ładujemy micepakiet i sprawdzamy, jakie metody wybierze domyślnie:

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

pmmOznacza predykcyjnej średniej dopasowywania - prawdopodobnie najbardziej popularnym algorytmem zaliczenia do przypisania zmiennych ciągłych. Oblicza przewidywaną wartość za pomocą modelu regresji i wybiera 5 elementów najbliższych przewidywanej wartości (według odległości euklidesowej ). Te wybrane elementy są nazywane pulą dawców, a ostateczna wartość jest wybierana losowo z tej puli dawców.

Z macierzy prognoz wynika, że ​​metody pobierają zmienne, które są interesujące dla ograniczeń. Zauważ, że wiersz jest zmienną docelową, a kolumna predyktorami. Gdyby x1 nie miał 1 w kolumnie x3, musielibyśmy dodać to do macierzy:imp_base$predictorMatrix["x1","x3"] <- 1

Teraz część zabawy, generowanie metod imputacji. Wybrałem tutaj dość prymitywną metodę, w której odrzucam wszystkie wartości, jeśli nie spełniają kryteriów. Może to skutkować długim czasem pętli i potencjalnie bardziej efektywne może być zachowanie prawidłowych imputacji i ponowne wykonanie tylko pozostałych, wymagałoby to jednak nieco drobniejszych poprawek.

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

Po zakończeniu definiowania metod, po prostu zmieniamy poprzednie metody. Jeśli chcesz zmienić tylko jedną zmienną, możesz po prostu użyć, imp_base$method["x2"] <- "pmm_x2"ale w tym przykładzie zmienimy wszystko (nazewnictwo nie jest konieczne):

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)

Teraz spójrzmy na trzeci przypisany zestaw danych:

> 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, to działa. Podoba mi się to rozwiązanie, ponieważ możesz korzystać z głównych funkcji i dodawać ograniczenia, które uważasz za istotne.

Aktualizacja

Aby wymusić rygorystyczne ograniczenia @ t0x1n wymienione w komentarzach, możemy chcieć dodać następujące funkcje do funkcji otoki:

1. Zapisz prawidłowe wartości podczas pętli, aby dane z poprzednich, częściowo udanych przebiegów nie zostały odrzucone
2. Mechanizm ucieczki w celu uniknięcia nieskończonych pętli
3. Napompuj pulę dawców po próbie x razy bez znalezienia odpowiedniego dopasowania (dotyczy to przede wszystkim pmm)

Powoduje to nieco bardziej skomplikowaną funkcję opakowania:

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

Zauważ, że to nie działa tak dobrze, najprawdopodobniej z powodu tego, że sugerowany zestaw danych nie spełnia ograniczeń dla wszystkich przypadków, nie tracąc ich. Muszę zwiększyć długość pętli do 400-500, zanim zacznie się ona zachowywać. Zakładam, że jest to niezamierzone, twoje przypisanie powinno naśladować sposób generowania rzeczywistych danych.

Optymalizacja

Argument ryzawiera brakujące wartości i prawdopodobnie moglibyśmy przyspieszyć pętlę, usuwając elementy, dla których znaleźliśmy kwalifikujące się przypisania, ale ponieważ nie jestem zaznajomiony z funkcjami wewnętrznymi, powstrzymałem się od tego.

Myślę, że najważniejszą rzeczą, kiedy masz silne ograniczenia, które wymagają pełnego wypełnienia, jest równoległe przypisanie imputacji ( zobacz moją odpowiedź na temat CrossValidated ). Większość ma dziś komputery z 4-8 rdzeniami, a R domyślnie używa tylko jednego z nich. Czas można (prawie) skrócić na pół, podwajając liczbę rdzeni.

Brakujące parametry przy imputacji

Jeśli chodzi o problem x2zaginięcia w momencie przypisania, myszy tak naprawdę nigdy nie wprowadzają brakujących wartości do x- data.frame. Metoda myszy obejmuje wypełnienie losowej wartości na początku. Część łańcuchowa imputacji ogranicza wpływ tej wartości początkowej. Jeśli spojrzysz na mice-funkcję, możesz to znaleźć przed wywołaniem imputacji ( mice:::sampler-funkcja):

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

data.initMoże być dostarczany do micefunkcji i mice.imput.sample to podstawowa procedura pobierania próbek.

Sekwencja odwiedzin

Jeśli kolejność odwiedzin jest ważna, możesz określić kolejność, w której micefunkcja uruchamia imputacje. Domyślnie jest z, 1:ncol(data)ale możesz ustawić dowolną wartość visitSequence.

Najbliższe, co mogłem znaleźć, to wcześniejsze włączenie informacji Amelii . Zobacz rozdział 4.7 winiety , w szczególności 4.7.2:

Priory na poziomie obserwacji

Naukowcy często dysponują dodatkowymi wcześniejszymi informacjami na temat brakujących wartości danych na podstawie wcześniejszych badań, konsensusu akademickiego lub osobistych doświadczeń. Amelia może uwzględnić te informacje, aby uzyskać znacznie lepsze przypisania. Algorytm Amelii pozwala użytkownikom uwzględniać pouczające priory bayesowskie o pojedynczych brakujących komórkach danych zamiast bardziej ogólnych parametrów modelu, z których wiele nie ma bezpośredniego znaczenia.

Włączenie priorytetów następuje po podstawowej analizie bayesowskiej, w której imputacja okazuje się średnią ważoną imputacji opartej na modelu i średniej wcześniejszej, gdzie wagi są funkcjami względnej siły danych i wcześniejszej: kiedy model bardzo dobrze przewiduje , przypisanie obniży wagę wcześniejszego i odwrotnie (Honaker i King, 2010).

Priorytety dotyczące poszczególnych obserwacji powinny opisywać przekonanie analityka dotyczące rozmieszczenia brakującej komórki danych. Może to przybrać formę średniej i odchylenia standardowego lub przedziału kondycji. Na przykład moglibyśmy wiedzieć, że stawki tariów z 1986 r. W Tajlandii wynoszą około 40%, ale mamy wątpliwości co do dokładnej wartości. Nasze wcześniejsze przekonanie o rozkładzie brakującej komórki danych koncentruje się zatem na 40 ze standardowym odchyleniem, które odzwierciedla ilość niepewności, jaką mamy odnośnie do naszego wcześniejszego przekonania.

Aby wprowadzić priors, musisz zbudować macierz priors z czterema lub pięcioma kolumnami. Każdy rząd macierzy reprezentuje pierwszeństwo jednej obserwacji lub jednej zmiennej. W dowolnym wierszu wpis w pierwszej kolumnie jest rzędem obserwacji, a wpis w drugiej kolumnie jest kolumną obserwacji. W czterokolumnowej macierzy priorytetów trzecia i czwarta kolumna są średnią i odchyleniem standardowym wcześniejszego rozkładu brakującej wartości.

Tak więc, chociaż ogólnie nie będziesz w stanie powiedzieć czegoś takiego x1<x2+x3, możesz zapętlić swój zestaw danych i dodać poziom obserwacji przed każdym odpowiednim przypadkiem. Można również zastosować stałe granice (takie jak ustawienie wartości x1, x2 i x3, aby były nieujemne). Na przykład:

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

Ograniczenia są prawdopodobnie łatwiejsze do wdrożenia w predykcyjnym dopasowaniu średnim wielokrotnej imputacji. Zakłada się, że istnieje znaczna liczba obserwacji z brakującymi zmiennymi ograniczającymi, które spełniają ograniczenia. Myślę o zaimplementowaniu tego w funkcji Hmiscpakietu R. aregImputeMożesz sprawdzić ponownie za około miesiąc. Ważne będzie określenie maksymalnej odległości od celu, jaką może być obserwacja dawcy, ponieważ ograniczenia popchną dawców dalej od idealnego, nieograniczonego dawcy.

Uważam, że Ameliapakiet (Amelia II) ma obecnie najbardziej kompleksowe wsparcie dla określania ograniczeń zakresu wartości danych. Problem polega jednak na tym Amelia, że dane są normalne na wielu odmianach.

Jeśli w twoim przypadku nie obowiązuje założenie wielowymiarowej normalności, możesz sprawdzić micepakiet, który implementuje wielokrotną imputację (MI) za pomocą równań łańcuchowych . Ten pakiet nie zakłada normalności wielowymiarowej . Ma również funkcję, która może wystarczyć do określenia ograniczeń , ale nie jestem pewien, do jakiego stopnia. Funkcja jest wywoływana squeeze(). Możesz o tym przeczytać w dokumentacji: http://cran.r-project.org/web/packages/mice/mice.pdf . Dodatkową zaletą micejest jego elastyczność w zakresie umożliwiania specyfikacji zdefiniowanych przez użytkownika funkcji imputacji i szerszego wyboru algorytmów. Oto samouczek na temat wykonywania MI przy użyciu mice:http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm .

O ile rozumiem, Hmiscpakiet dr. Harrella , wykorzystujący takie same łańcuchowe równania ( predykcyjne dopasowywanie średnich ), prawdopodobnie obsługuje dane nienormalne (z wyjątkiem normpmmmetody). Być może zaimplementował już funkcjonalność specyfikacji ograniczeń zgodnie z powyższą odpowiedzią. Nie korzystałem aregImpute(), więc nie mogę powiedzieć więcej na ten temat (korzystałem Ameliai mice, ale zdecydowanie nie jestem ekspertem w dziedzinie statystyki, staram się tylko nauczyć jak najwięcej).

Na koniec możesz znaleźć interesujące, nieco przestarzałe, ale wciąż ładne, omówienie podejść, metod i oprogramowania do wielokrotnego przypisywania danych z brakującymi wartościami: http://www.ncbi.nlm.nih.gov/pmc/articles / PMC1839993 . Jestem pewien, że są nowsze artykuły przeglądowe na temat MI, ale to wszystko, o czym wiem obecnie. Mam nadzieję, że jest to nieco pomocne.

Jeśli dobrze rozumiem twoje pytanie, wydaje mi się, że już wiesz, jakie wartości powinny przyjąć brakujące zmienne z zastrzeżeniem pewnych ograniczeń. Nie jestem zbyt biegły w SPSS, ale w RI myślę, że możesz napisać funkcję, aby to zrobić (co nie powinno być zbyt trudne w zależności od twojego doświadczenia, powinienem powiedzieć). Nie znam żadnego pakietu, który działałby z takimi ograniczeniami.