Wielokrotna imputacja brakujących wartości

13

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
Róża
źródło
Zmieniłem 0 or 16 or >= 16na, 0 or >= 16ponieważ >=16zawiera wartość 16. Mam nadzieję, że to nie zepsuło twojego znaczenia. To samo dotyczy0 or 14 or >= 14
Alexis,

Odpowiedzi:

16

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.

Max Gordon
źródło
+1 To świetne rzeczy, dokładnie to, co miałem na myśli (zobacz mój komentarz do odpowiedzi Franka), a na pewno kandydat numer 1 na nagrodę. Kilka rzeczy o mnie kłopot pmm_x1jednak: (1) Biorąc pod uwagę sumę maksymalną każdej możliwej kombinacji x2i x3od całego zestawu danych jest znacznie Stringer niż oryginalny ograniczeń. Prawidłowe rzeczą byłoby sprawdzić, czy dla każdego wiersza , x1 < x2 + x3. Oczywiście im więcej masz rzędów, tym mniejsza jest szansa na spełnienie takiego ograniczenia (ponieważ jeden zły rząd rujnuje wszystko) i tym dłużej może potencjalnie uzyskać pętlę.
t0x1n
(2) jeśli oba x1i x2brakuje, można przypisać wartość x1, dla których ograniczenia są utrzymywane (powiedzmy 50), ale gdy x2zostanie przypisane są łamane (powiedzmy to być przypisana do 55). Czy istnieje sposób przypisywania „w poziomie”, a nie w pionie? W ten sposób możemy przypisać jeden wiersz x1, x2oraz x3i po prostu ponownie przypisać go do momentu, że specyficzny rząd podlega ograniczeniom. To powinno być wystarczająco szybkie, a gdy to zrobimy, możemy przejść do następnego rzędu. Oczywiście, jeśli MI ma charakter „pionowy”, nie mamy szczęścia. W takim razie może podejście, o którym wspominał Aleksandr?
t0x1n
Fajne rozwiązanie, +1! Może być szczególnie przydatny, ponieważ obecnie używam micepakietu. Dzięki za udostępnienie.
Aleksandr Blekh
1
@ t0x1n Zaktualizowałem odpowiedź o bardziej zaawansowaną funkcję opakowania zgodnie z Twoimi komentarzami. Jeśli chcesz nurkować głębiej, zalecamy zabawę z nim debug(), aby zobaczyć, jak mice.impute.pmmi jego rodzeństwo działa pod maską.
Max Gordon
1
@ t0x1n: Chyba - sprawdź swoje przypisane wartości. Jeśli wydają się nierealne, możesz wybrać moje podejście, aby przypisać tylko te, które są mniej centralne dla modelu. W moim przypadku zdecydowałem się wykluczyć te, które nie są kontrolowane promieniami rentgenowskimi, ponieważ są one w centrum badań, a przypisania nie zapewniają klinicznie wiarygodnych wartości (wydłużenie nogi po złamaniu). Nie jestem do końca z tego zadowolony, ale wydaje się to rozsądnym kompromisem.
Max Gordon
8

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);
t0x1n
źródło
5

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.

Frank Harrell
źródło
Też chciałbym to mieć. Powiedzmy, że potrzebuję tylko najbardziej podstawowych ograniczeń między zmiennymi x<y<z.
t0x1n
Wybacz mi moją niewiedzę, jeśli jestem daleko, ale miałem wrażenie, że wiele technik imputacji wymaga wyciągania wartości z odpowiedniego rozkładu. Czy nie powinno być prostą sprawą zastosowanie próbkowania odrzucającego? np. kontynuować rysowanie, dopóki nie x1<x2zostanie spełnione określone ograniczenie (np. )?
t0x1n
To właśnie mogę zrobić z aregImputefunkcją R z predykcyjnym dopasowaniem średnim. Ale co, jeśli żadna obserwacja dawcy (w pobliżu dopasowań prognoz) nie spełnia ograniczeń przypisywanych obserwacji celu, chociaż oczywiście musiały one spełniać ograniczenia zbioru zmiennych dawcy?
Frank Harrell,
W takim przypadku może wziąć bezpośrednio przewidywaną wartość? Czy taka próbka polega tylko na regresji (bez fazy PMM)?
t0x1n
Imputacja regresji jest nieco bardziej prawdopodobna w przypadku wartości przypisywanych, które są niespójne z resztą zapisu podmiotu. Więc nie sądzę, że jest to powód do unikania PMM.
Frank Harrell,
4

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.

Aleksandr Blekh
źródło
1
Ten miły komentarz sprawia, że ​​myślę, że predykcyjne porównywanie średnich, które zastępuje braki faktycznie obserwowanymi wartościami, może już zawierać pewne ograniczenia, jeśli wszystkie obserwowane dane spełniają te ograniczenia. Byłbym wdzięczny za przemyślenie tego. Nie wdrożyłem jeszcze żadnych specjalnych ograniczeń w aregImpute.
Frank Harrell,
1
Masz rację. Właśnie zdałem sobie sprawę, że obserwacje dawcy dostarczają wartości, które są zgodne z ich innymi zmiennymi, ale nie z innymi zmiennymi w zmiennej docelowej.
Frank Harrell,
1
Czy oprócz założeń dotyczących dystrybucji dokonanych przez Amelię, byłeś w stanie sprecyzować ograniczenia bardziej szczegółowo, niż wykazałem w mojej odpowiedzi? Problem squeezepolega na tym, że jego granice są stałe, więc nie można podać niczego podobnego x1<x2. Wydaje się również, że jest on wywoływany na domniemanym wektorze wyników, który moim zdaniem jest za późno. Wydaje mi się, że granice należy brać pod uwagę podczas procesu przypisywania, więc mają one większe znaczenie niż korekta po fakcie.
t0x1n
1
@ t0x1n: Niestety nie miałem okazji sprecyzować ograniczeń Amelia, ponieważ przeszedłem z niego na mice, gdy tylko moje testy potwierdziły, że moje dane nie są normalne na wielu odmianach. Jednak ja niedawno biegł tym bardzo ładnym zestawem slajdów prezentacji na temat (MI metod i oprogramowania): statistik.lmu.de/~fkreuter/imputation_sose2011/downloads/... . Jeśli dobrze zrozumiałem, opisuje potencjalne rozwiązanie problemu ograniczeń (patrz strona 50 w PDF - nie slajd nr 50!). Mam nadzieję że to pomoże.
Aleksandr Blekh
1
@ t0x1n: W rzeczywistości rozwiązanie opisano na stronach 50 i 51.
Aleksandr Blekh
0

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.

ThinkStatsme
źródło