Estimação imparcial da matriz de covariância para dados multiplamente censurados

22

As análises químicas de amostras ambientais são frequentemente censuradas abaixo nos limites de relatório ou em vários limites de detecção / quantificação. Este último pode variar, geralmente na proporção dos valores de outras variáveis. Por exemplo, uma amostra com uma alta concentração de um composto pode precisar ser diluída para análise, resultando em inflação proporcional dos limites de censura para todos os outros compostos analisados ​​ao mesmo tempo nessa amostra. Como outro exemplo, às vezes a presença de um composto pode alterar a resposta do teste a outros compostos (uma "interferência da matriz"); quando isso for detectado pelo laboratório, aumentará seus limites de relatório de acordo.

Estou procurando uma maneira prática de estimar toda a matriz de variância-covariância para esses conjuntos de dados, especialmente quando muitos dos compostos experimentam mais de 50% de censura, o que geralmente é o caso. Um modelo distributivo convencional é que os logaritmos das concentrações (verdadeiras) são multinormalmente distribuídos, e isso parece se encaixar bem na prática; portanto, uma solução para essa situação seria útil.

(Por "prático", refiro-me a um método que pode ser codificado de maneira confiável em pelo menos um ambiente de software geralmente disponível, como R, Python, SAS etc., de uma maneira que seja executada com rapidez suficiente para suportar recálculos iterativos, como ocorrem em imputações múltiplas, e que é razoavelmente estável [é por isso que reluto em explorar uma implementação do BUGS, embora as soluções bayesianas em geral sejam bem-vindas].)

Muito obrigado antecipadamente por seus pensamentos sobre este assunto.

whuber
fonte
Apenas para entender corretamente a questão da censura: Quando você dilui uma amostra, a concentração de um composto cai tão baixo que o instrumento de teste pode falhar em detectar sua presença. É uma reformulação precisa do problema de censura?
Sim, isso está correto: a diluição por um fator de D também aumenta todos os limites de detecção por um fator de D. (A questão da interferência da matriz é mais difícil de quantificar e a situação geral é extremamente complexa. Para simplificar isso, o modelo convencional é que um conjunto de testes em uma amostra produz um vetor (x [1], ..., x [k ]) onde x [i] é um número real ou é um intervalo de reais, normalmente com o ponto final esquerdo no infinito; um intervalo identifica um conjunto no qual o valor verdadeiro é assumido.)
whuber
Por que os limites de detecção aumentariam? Elas não são uma característica do instrumento de teste e não a amostra que está sendo testada?
Como exemplo, suponha que o limite de detecção de um instrumento seja de 1 micrograma por litro (ug / L). Uma amostra é diluída 10: 1 (com grande precisão, para não nos preocuparmos com erros aqui) e o instrumento lê "<1"; isto é, não detectável, para a amostra diluída. O laboratório deduz que a concentração na amostra é menor que 10 * 1 = 10 µg / L e a informa como tal; isto é, como "<10".
whuber
1
@amoeba Vejo que eu deveria ter explicado essas coisas na própria pergunta. As respostas são: PCA; a dimensionalidade varia de 3 a algumas centenas; os tamanhos das amostras sempre excedem muito a dimensionalidade, mas as taxas de censura podem ser muito altas (é possível lidar com até 50% e é desejável até 95%).
whuber

Respostas:

3

Ainda não internalizei completamente a questão da interferência matricial, mas aqui está uma abordagem. Deixei:

é um vetor que representa a concentração de todos os compostos alvo na amostra não diluída.Y

é o vetor correspondente na amostra diluída.Z

é o fator de diluição, ou seja, a amostra é diluída d : 1.dd

Nosso modelo é:

YN(μ,Σ)

Z=Yd+ϵ

onde representa o erro devido a erros de diluição.ϵN(0,σ2 I)

Portanto, segue-se que:

ZN(μd,Σ+σ2 I)

Denote a distribuição acima de por f Z ( . ) .ZfZ(.)

Seja as concentrações observadas e τ represente o limiar do instrumento de teste abaixo do qual ele não pode detectar um composto. Então, para o i t h composto, temos:Oτith

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

Sem perda de generalidade, deixe o primeiro k compostos sejam tais que estejam abaixo do limiar. Em seguida, a função de probabilidade pode ser escrita como:

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

Onde

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

A estimativa é, então, uma questão de usar a probabilidade máxima ou idéias bayesianas. Não sei ao certo o quão tratável é o acima, mas espero que lhe dê algumas idéias.


fonte
Muito obrigado por esse pensamento. De fato, essa é uma abordagem padrão e bem documentada para a censura múltipla. Uma dificuldade está na sua intratabilidade: essas integrais são notoriamente difíceis de calcular. Também há um problema de modelagem oculto aqui: o valor de d geralmente está correlacionado positivamente com Y , como está implícito no primeiro parágrafo da minha descrição.
whuber
2

Outra opção mais computacionalmente eficiente seria ajustar a matriz de covariância pela correspondência de momentos usando um modelo que foi chamado de "gaussiano dichomizado", na verdade apenas um modelo de cópula gaussiana.

Um artigo recente de Macke et al. 2010 descreve um procedimento de formulário fechado para o ajuste deste modelo, que envolve apenas a matriz de covariância empírica (censurada) e o cálculo de algumas probabilidades normais bivariadas. O mesmo grupo (laboratório de Bethge da MPI Tuebingen) também descreveu modelos gaussianos discretos / contínuos híbridos que provavelmente são o que você deseja aqui (ou seja, como os RVs gaussianos não são totalmente "dicotomizados" - apenas aqueles abaixo do limite).

Criticamente, esse não é um estimador de ML, e temo não saber quais são suas propriedades de polarização.

jpillow
fonte
@ jp Obrigado: vou olhar para isso. (Pode demorar um pouco ...)
whuber
1

Quantos compostos existem na sua amostra? (Ou qual é o tamanho da matriz de covariância em questão?).

Alan Genz tem um código muito interessante em várias linguagens (R, Matlab, Fortran; veja aqui ) para calcular integrais de densidades normais multivariadas sobre hiper retângulos (ou seja, os tipos de integrais necessários para avaliar a probabilidade, conforme observado por usuário28).

Eu usei essas funções ("ADAPT" e "QSIMVN") para integrais de até 10 a 12 dimensões, e várias funções nessa página anunciam integrais (e derivados associados que você pode precisar) para problemas até a dimensão 100. Eu não você sabe se essas dimensões são suficientes para seus propósitos, mas, se for o caso, presumivelmente, poderá encontrar estimativas de probabilidade máxima por subida gradiente.

jpillow
fonte
Oh, desculpe, eu sou novo por aqui e não percebi há quanto tempo isso foi publicado - provavelmente muito tarde para ser de muita ajuda!
jpillow
@jp Este é um problema importante em andamento, portanto o tempo decorrido entre a pergunta e a resposta é de pouca importância. Obrigado por responder!
whuber