Bezstronna ocena macierzy kowariancji dla wielokrotnie cenzurowanych danych

Analizy chemiczne próbek środowiskowych są często cenzurowane poniżej limitów sprawozdawczych lub różnych limitów wykrywalności / ilościowych. Te ostatnie mogą się różnić, zwykle proporcjonalnie do wartości innych zmiennych. Na przykład, próbka o wysokim stężeniu jednego związku może wymagać rozcieńczenia do analizy, co spowoduje proporcjonalne zawyżenie limitów cenzury dla wszystkich innych związków analizowanych jednocześnie w tej próbce. Jako inny przykład, czasami obecność związku może zmienić odpowiedź testu na inne związki („interferencja matrycy”); kiedy laboratorium wykryje to, odpowiednio zwiększy swoje limity raportowania.

Szukam praktycznego sposobu oszacowania całej macierzy wariancji-kowariancji dla takich zestawów danych, zwłaszcza gdy wiele związków doświadcza ponad 50% cenzury, co często ma miejsce. Konwencjonalny model dystrybucji polega na tym, że logarytmy (prawdziwych) stężeń są rozkładane wielonormalnie, co wydaje się dobrze pasować w praktyce, więc przydatne byłoby rozwiązanie tej sytuacji.

(Przez „praktyczny” rozumiem metodę, którą można niezawodnie zakodować w co najmniej jednym ogólnie dostępnym środowisku oprogramowania, takim jak R, Python, SAS itp., W sposób, który wykonuje się wystarczająco szybko, aby obsługiwać iteracyjne ponowne obliczenia, takie jak wielokrotne przypisywanie, i który jest dość stabilny [dlatego niechętnie badam implementację BŁĘDU, chociaż rozwiązania bayesowskie są ogólnie mile widziane])

Z góry dziękuję za przemyślenia na ten temat.

Nie w pełni zinternalizowałem problem interferencji macierzy, ale oto jedno podejście. Pozwolić:

$Y$ jest wektorem, który reprezentuje stężenie wszystkich docelowych związków w nierozcieńczonej próbce.

$Z$ oznacza odpowiedni wektor w rozcieńczonej próbce.

$d$ jest współczynnikiem rozcieńczenia, tzn. próbka jest rozcieńczana : 1.$d$

Nasz model to:

$Y \sim N(\mu,\Sigma)$

$Z = \frac{Y}{d} + \epsilon$

gdzie reprezentuje błąd wynikający z błędów rozcieńczania.$\epsilon \sim N(0,\sigma^2\ I)$

Wynika stąd, że:

$Z \sim N(\frac{\mu}{d}, \Sigma + \sigma^2\ I)$

Oznacz powyższy rozkład przez .$Z$$f_Z(.)$

Niech będzie obserwowanymi stężeniami, a reprezentuje próg przyrządu testowego, poniżej którego nie może wykryć związku. Następnie dla związku mamy:$O$$\tau$$i^{th}$

$O_i = Z_i I(Z_i > \tau) + 0 I(Z_i \le \tau)$

Bez utraty ogólności niech pierwsze związków będzie takie, aby były poniżej progu. Następnie funkcję prawdopodobieństwa można zapisać jako:$k$

$L(O_1, ... O_k, O_{k+1},...O_n |- ) = [\prod_{i=1}^{i=k}{Pr(Z_i \le \tau)}] [\prod_{i=k+1}^{i=n}{f(O_i |-)}]$

gdzie

$f(O_i |-) = \int_{j\neq i}{f_Z(O_i|-) I(O_i > \tau)}$

Oszacowanie polega zatem na wykorzystaniu albo maksymalnego prawdopodobieństwa, albo pomysłów bayesowskich. Nie jestem pewien, na ile powyższe jest wykonalne, ale mam nadzieję, że dostarczy ci kilku pomysłów.

Inną bardziej wydajną obliczeniowo opcją byłoby dopasowanie macierzy kowariancji poprzez dopasowanie momentu za pomocą modelu, który został nazwany „dychomizowanym gaussowskim”, tak naprawdę tylko modelem kopuły Gaussa.

Niedawny artykuł Macke i in. 2010 opisuje procedurę zamkniętej formy dopasowania tego modelu, która obejmuje tylko (ocenzurowaną) empiryczną macierz kowariancji i obliczenie niektórych dwuwymiarowych normalnych prawdopodobieństw. Ta sama grupa (laboratorium Bethge'a z MPI Tuebingen) opisała również hybrydowe dyskretne / ciągłe modele gaussowskie, które prawdopodobnie są tutaj potrzebne (tj. Ponieważ Gaussowskie RV nie są w pełni „dychotomizowane” - tylko te poniżej progu).

Krytycznie nie jest to oszacowanie ML i obawiam się, że nie wiem, jakie są jego właściwości uprzedzające.

Ile związków jest w twojej próbce? (Lub, jak duża jest omawiana macierz kowariancji?).

Alan Genz ma bardzo ładny kod w różnych językach (R, Matlab, Fortran; patrz tutaj ) do obliczania całek wielowymiarowych normalnych gęstości w hiperprostokątach (tj. Rodzajów całek, których potrzebujesz do oceny prawdopodobieństwa, jak zauważono przez użytkownik 28).

Użyłem tych funkcji („ADAPT” i „QSIMVN”) dla całek o wielkości do około 10-12 wymiarów, a kilka funkcji na tej stronie reklamuje całki (i powiązane pochodne, których możesz potrzebować) w przypadku problemów do wymiaru 100. Nie nie wiem, czy jest to wystarczająca liczba wymiarów dla twoich celów, ale jeśli tak, to prawdopodobnie pozwoli ci znaleźć szacunki maksymalnego prawdopodobieństwa na podstawie wzrostu gradientu.