Unvoreingenommene Schätzung der Kovarianzmatrix für mehrfach zensierte Daten

Chemische Analysen von Umweltproben werden im Folgenden häufig an Meldegrenzen oder verschiedenen Nachweis- / Bestimmungsgrenzen zensiert. Letztere können variieren, normalerweise proportional zu den Werten anderer Variablen. Beispielsweise muss möglicherweise eine Probe mit einer hohen Konzentration einer Verbindung zur Analyse verdünnt werden, was zu einem proportionalen Aufpumpen der Zensurgrenzwerte für alle anderen Verbindungen führt, die zur gleichen Zeit in dieser Probe analysiert werden. Als ein anderes Beispiel kann manchmal das Vorhandensein einer Verbindung die Reaktion des Tests auf andere Verbindungen verändern (eine "Matrixinterferenz"); Wenn dies vom Labor festgestellt wird, werden die Meldegrenzen entsprechend erhöht.

Ich suche einen praktischen Weg, um die gesamte Varianz-Kovarianz-Matrix für solche Datensätze abzuschätzen, insbesondere wenn viele der Verbindungen zu mehr als 50% zensiert werden, was häufig der Fall ist. Ein herkömmliches Verteilungsmodell besteht darin, dass die Logarithmen der (wahren) Konzentrationen multinormal verteilt sind, und dies scheint in der Praxis gut zu passen, so dass eine Lösung für diese Situation nützlich wäre.

(Mit "praktisch" meine ich eine Methode, die zuverlässig in mindestens einer allgemein verfügbaren Softwareumgebung wie R, Python, SAS usw. so codiert werden kann, dass sie schnell genug ausgeführt werden kann, um iterative Neuberechnungen zu unterstützen, wie sie beispielsweise bei Mehrfachimputationen auftreten. und die einigermaßen stabil ist (weshalb ich eine BUGS-Implementierung nur ungern untersuchen möchte, obwohl Bayes'sche Lösungen im Allgemeinen erwünscht sind).

Vielen Dank im Voraus für Ihre Gedanken zu diesem Thema.

Ich habe das Problem der Matrixinterferenz noch nicht vollständig verinnerlicht, aber hier ist ein Ansatz. Lassen:

$Y$ ist ein Vektor, der die Konzentration aller Zielverbindungen in der unverdünnten Probe darstellt.

$Z$ ist der entsprechende Vektor in der verdünnten Probe.

$d$ ist der Verdünnungsfaktor, dh die Probe wird mit : 1 verdünnt .$d$

Unser Modell ist:

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

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

Dabei steht für den Fehler aufgrund von Verdünnungsfehlern.$\epsilon \sim N(0,\sigma^2\ I)$

Daraus folgt:

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

Bezeichne die obige Verteilung von mit .f Z ( . $Z$$f_Z(.)$

Sei die beobachtete Konzentration und der Schwellenwert des Testinstruments, unter dem es keine Verbindung nachweisen kann. Dann haben wir für die -Verbindung:τ i t $O$$\tau$$i^{th}$

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

Ohne Verlust der Allgemeinheit seien die ersten Verbindungen so, dass sie unterhalb der Schwelle liegen. Dann kann die Wahrscheinlichkeitsfunktion geschrieben werden als:$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 |-)}]$

woher

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

Bei der Schätzung kommt es dann darauf an, entweder die maximale Wahrscheinlichkeit oder Bayes'sche Ideen zu verwenden. Ich bin mir nicht sicher, wie gut das oben genannte funktioniert, aber ich hoffe, es gibt Ihnen einige Anregungen.

Eine weitere rechnerisch effizientere Option wäre die Anpassung der Kovarianzmatrix durch Momentananpassung unter Verwendung eines Modells, das als "dichomisiertes Gauß-Modell" bezeichnet wird, das eigentlich nur ein Gauß-Copula-Modell ist.

Ein kürzlich veröffentlichter Aufsatz von Macke et al. 2010 beschreibt ein Verfahren in geschlossener Form zur Anpassung dieses Modells, das nur die (zensierte) empirische Kovarianzmatrix und die Berechnung einiger bivariater normaler Wahrscheinlichkeiten umfasst. Dieselbe Gruppe (Bethge-Labor am MPI Tübingen) hat auch hybride diskrete / kontinuierliche Gauß-Modelle beschrieben, die wahrscheinlich das sind, was Sie hier wollen (dh da die Gauß-Wohnmobile nicht vollständig "dichotomisiert" sind - nur solche unterhalb der Schwelle).

Entscheidend ist, dass dies kein ML-Schätzer ist, und ich fürchte, ich weiß nicht, wie seine Verzerrungseigenschaften sind.

Wie viele Verbindungen enthält Ihre Probe? (Oder wie groß ist die betreffende Kovarianzmatrix?)

Alan Genz hat einen sehr guten Code in einer Vielzahl von Sprachen (R, Matlab, Fortran; siehe hier ) zur Berechnung von Integralen multivariater normaler Dichten über Hyperrechtecke (dh die Arten von Integralen, die Sie zur Bewertung der Wahrscheinlichkeit benötigen, wie von angegeben) user28).

Ich habe diese Funktionen ("ADAPT" und "QSIMVN") für Integrale bis zu 10-12 Dimensionen verwendet, und mehrere Funktionen auf dieser Seite bieten Integrale (und zugehörige Ableitungen, die Sie möglicherweise benötigen) für Probleme bis zur Dimension 100 an Sie wissen nicht, ob das für Ihre Zwecke ausreicht, aber in diesem Fall können Sie wahrscheinlich Schätzungen der maximalen Wahrscheinlichkeit anhand des Gradientenaufstiegs finden.