Estimación imparcial de la matriz de covarianza para datos censurados múltiples

Los análisis químicos de las muestras ambientales a menudo se censuran a continuación en los límites de notificación o en varios límites de detección / cuantificación. Este último puede variar, generalmente en proporción a los valores de otras variables. Por ejemplo, una muestra con una alta concentración de un compuesto puede necesitar diluirse para el análisis, lo que resulta en una inflación proporcional de los límites de censura para todos los demás compuestos analizados al mismo tiempo en esa muestra. Como otro ejemplo, a veces la presencia de un compuesto puede alterar la respuesta de la prueba a otros compuestos (una "interferencia de matriz"); cuando el laboratorio lo detecte, aumentará sus límites de notificación en consecuencia.

Estoy buscando una forma práctica de estimar la matriz de varianza-covarianza completa para dichos conjuntos de datos, especialmente cuando muchos de los compuestos experimentan más del 50% de censura, que a menudo es el caso. Un modelo de distribución convencional es que los logaritmos de las concentraciones (verdaderas) están distribuidos de manera multinormal, y esto parece encajar bien en la práctica, por lo que sería útil una solución para esta situación.

(Por "práctico" me refiero a un método que puede codificarse de manera confiable en al menos un entorno de software generalmente disponible como R, Python, SAS, etc., de una manera que se ejecute lo suficientemente rápido como para admitir recálculos iterativos como los que ocurren en la imputación múltiple, y que es razonablemente estable [es por eso que soy reacio a explorar una implementación de ERRORES, aunque las soluciones bayesianas en general son bienvenidas].

Muchas gracias de antemano por sus pensamientos sobre este asunto.

No he internalizado por completo el tema de la interferencia de la matriz, pero aquí hay un enfoque. Dejar:

sea ​​un vector que represente la concentración de todos los compuestos objetivo en la muestra sin diluir.$Y$

$Z$

$d$$d$

Nuestro modelo es:

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

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

$\epsilon \sim N(0,\sigma^2\ I)$

Por lo tanto, se deduce que:

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

$Z$$f_Z(.)$

$O$$\tau$$i^{th}$

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

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

dónde

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

La estimación es, entonces, una cuestión de utilizar la máxima probabilidad o ideas bayesianas. No estoy seguro de cuán manejable es lo anterior, pero espero que te dé algunas ideas.

Otra opción más eficiente desde el punto de vista computacional sería ajustar la matriz de covarianza por coincidencia de momentos utilizando un modelo que se ha denominado "Gaussiano dicomizado", realmente solo un modelo de cópula gaussiana.

Un artículo reciente de Macke et al 2010 describe un procedimiento de forma cerrada para ajustar este modelo que involucra solo la matriz de covarianza empírica (censurada) y el cálculo de algunas probabilidades normales bivariadas. El mismo grupo (laboratorio de Bethge en MPI Tuebingen) también ha descrito modelos gaussianos híbridos discretos / continuos que probablemente son lo que desea aquí (es decir, ya que los RV gaussianos no están completamente "dicotomizados", solo aquellos por debajo del umbral).

Críticamente, este no es un estimador de ML, y me temo que no sé cuáles son sus propiedades de sesgo.

¿Cuántos compuestos hay en su muestra? (O, ¿qué tan grande es la matriz de covarianza en cuestión?).

Alan Genz tiene un código muy agradable en una variedad de lenguajes (R, Matlab, Fortran; ver aquí ) para calcular integrales de densidades normales multivariadas sobre hiper-rectángulos (es decir, los tipos de integrales que necesita para evaluar la probabilidad, como lo señala usuario28).

He usado estas funciones ("ADAPT" y "QSIMVN") para integrales de hasta 10-12 dimensiones, y varias funciones en esa página anuncian integrales (y derivados asociados que pueda necesitar) para problemas hasta la dimensión 100. No lo uso No sé si eso es suficiente para sus propósitos, pero si es así, presumiblemente podría permitirle encontrar estimaciones de máxima probabilidad por gradiente de ascenso.