Estimation impartiale de la matrice de covariance pour les données à censure multiple

22

Les analyses chimiques d'échantillons environnementaux sont souvent censurées ci-dessous aux limites de déclaration ou à diverses limites de détection / quantification. Ces dernières peuvent varier, généralement proportionnellement aux valeurs des autres variables. Par exemple, un échantillon avec une concentration élevée d'un composé pourrait devoir être dilué pour l'analyse, ce qui entraînerait une inflation proportionnelle des limites de censure pour tous les autres composés analysés en même temps dans cet échantillon. Comme autre exemple, parfois la présence d'un composé peut modifier la réponse du test à d'autres composés (une «interférence matricielle»); lorsque cela est détecté par le laboratoire, il gonfle ses limites de déclaration en conséquence.

Je cherche un moyen pratique d'estimer la matrice variance-covariance pour de tels ensembles de données, en particulier lorsque de nombreux composés subissent une censure supérieure à 50%, ce qui est souvent le cas. Un modèle de distribution conventionnel est que les logarithmes des (vraies) concentrations sont distribués multinormalement, et cela semble bien cadrer dans la pratique, donc une solution à cette situation serait utile.

(Par "pratique", j'entends une méthode qui peut être codée de manière fiable dans au moins un environnement logiciel généralement disponible comme R, Python, SAS, etc., d'une manière qui s'exécute assez rapidement pour prendre en charge des recalculs itératifs tels que ceux qui se produisent lors d'une imputation multiple, et qui est raisonnablement stable [c'est pourquoi je suis réticent à explorer une implémentation de BUGS, bien que les solutions bayésiennes en général soient les bienvenues].)

Merci d'avance pour vos réflexions à ce sujet.

whuber
la source
Pour que je comprenne bien le problème de la censure: lorsque vous diluez un échantillon, la concentration d'un composé tombe si bas que l'instrument de test peut ne pas détecter sa présence. Est-ce une reformulation précise du problème de censure?
Oui, c'est exact: la dilution d'un facteur D augmente également toutes les limites de détection d'un facteur D. (Le problème des interférences matricielles est plus difficile à quantifier et la situation générale est extrêmement complexe. Pour simplifier, le modèle conventionnel est qu'une suite de tests sur un échantillon donne un vecteur (x [1], ..., x [k ]) où les x [i] sont soit des nombres réels soit des intervalles de réels, généralement avec le point d'extrémité gauche à -infinity; un intervalle identifie un ensemble dans lequel la valeur vraie est supposée se trouver.)
whuber
Pourquoi les limites de détection augmenteraient-elles? Ne sont-elles pas une caractéristique de l'instrument de test plutôt que celle de l'échantillon testé?
Par exemple, supposons que la limite de détection d'un instrument soit de 1 microgramme par litre (ug / L). Un échantillon est dilué à 10: 1 (avec une grande précision, donc nous ne nous inquiétons pas des erreurs ici) et l'instrument lit "<1"; c'est-à-dire non détectable pour l'échantillon dilué. Le laboratoire en déduit que la concentration dans l'échantillon est inférieure à 10 * 1 = 10 ug / L et la signale comme telle; c'est-à-dire "<10".
whuber
1
@amoeba Je vois que j'aurais dû expliquer ces choses dans la question elle-même. Les réponses sont: PCA; la dimensionnalité variera de 3 à quelques centaines; la taille des échantillons dépasse toujours largement la dimensionnalité mais les taux de censure peuvent être très élevés (pouvoir gérer jusqu'à 50% est nécessaire et jusqu'à 95% est souhaitable).
whuber

Réponses:

3

Je n'ai pas complètement intériorisé le problème des interférences matricielles mais voici une approche. Laisser:

être un vecteur qui représente la concentration de tous les composés cibles dans l'échantillon non dilué.Y

est le vecteur correspondant dans l'échantillon dilué.Z

est le facteur de dilution, c'est-à-dire que l'échantillon est dilué d : 1.dd

Notre modèle est:

YN(μ,Σ)

Z=Yd+ϵ

représente l'erreur due aux erreurs de dilution.ϵN(0,σ2 I)

Il s'ensuit donc que:

ZN(μd,Σ+σ2 I)

Notons la distribution ci-dessus de par f Z ( . ) .ZfZ(.)

Soit les concentrations observées et τ représente le seuil de l'instrument de test en dessous duquel il ne peut pas détecter un composé. Ensuite, pour le composé i t h , nous avons:Oτith

Oi=ZiI(Zi>τ)+0I(Ziτ)

Sans perte de généralité, que les premiers composés soient tels qu'ils soient inférieurs au seuil. La fonction de vraisemblance peut alors s'écrire:k

L(O1,...Ok,Ok+1,...On|)=[i=1i=kPr(Ziτ)][i=k+1i=nf(Oi|)]

f(Oi|)=jifZ(Oi|)I(Oi>τ)

L'estimation consiste alors à utiliser soit le maximum de vraisemblance, soit les idées bayésiennes. Je ne sais pas à quel point ce qui précède est traitable, mais j'espère que cela vous donnera quelques idées.


la source
Merci beaucoup pour cette pensée. En effet, il s'agit d'une approche standard et bien documentée de la censure multiple. Une difficulté réside dans son caractère intraitable: ces intégrales sont notoirement difficiles à calculer. Il y a aussi un problème de modélisation qui se cache ici: la valeur de d est généralement positivement corrélée avec Y , comme l'indique le premier paragraphe de ma description.
whuber
2

Une autre option plus efficace en termes de calcul serait d'ajuster la matrice de covariance en faisant correspondre les moments en utilisant un modèle qui a été appelé "gaussien dichomisé", en réalité juste un modèle de copule gaussien.

Un article récent de Macke et al 2010 décrit une procédure de forme fermée pour ajuster ce modèle qui implique uniquement la matrice de covariance empirique (censurée) et le calcul de quelques probabilités normales bivariées. Le même groupe (laboratoire de Bethge au MPI Tuebingen) a également décrit des modèles gaussiens discrets / continus hybrides qui sont probablement ce que vous voulez ici (c'est-à-dire, puisque les VR gaussiens ne sont pas complètement "dichotomisés" - seulement ceux en dessous du seuil).

Surtout, ce n'est pas un estimateur ML, et je crains de ne pas savoir quelles sont ses propriétés de biais.

jpillow
la source
@jp Merci: je vais examiner cela. (Cela peut prendre un certain temps ...)
whuber
1

Combien de composés y a-t-il dans votre échantillon? (Ou, quelle est la taille de la matrice de covariance en question?).

Alan Genz a un très bon code dans une variété de langages (R, Matlab, Fortran; voir ici ) pour calculer les intégrales de densités normales multivariées sur des hyper-rectangles (c'est-à-dire les types d'intégrales dont vous avez besoin pour évaluer la probabilité, comme indiqué par utilisateur28).

J'ai utilisé ces fonctions ("ADAPT" et "QSIMVN") pour les intégrales jusqu'à environ 10-12 dimensions, et plusieurs fonctions sur cette page annoncent les intégrales (et les dérivés associés dont vous pourriez avoir besoin) pour les problèmes jusqu'à la dimension 100. Je ne Je ne sais pas si c'est assez de dimensions pour vos besoins, mais si c'est le cas, cela pourrait probablement vous permettre de trouver des estimations du maximum de vraisemblance par ascension du gradient.

jpillow
la source
Oh, désolé - je suis nouveau ici et je n'ai pas remarqué depuis combien de temps cela a été publié - probablement trop tard pour être d'une grande aide!
jpillow
@jp Il s'agit d'un problème important et permanent, de sorte que le temps écoulé entre la question et la réponse n'a que peu d'importance. Merci d'avoir répondu!
whuber