Несмещенная оценка ковариационной матрицы для данных с множественной цензурой

22

Химические анализы проб окружающей среды часто подвергаются цензуре ниже пределов отчетности или различных пределов обнаружения / количественного определения. Последние могут варьироваться, как правило, пропорционально значениям других переменных. Например, для анализа может потребоваться разведение образца с высокой концентрацией одного соединения, что приведет к пропорциональному раздуванию пределов цензуры для всех других соединений, анализируемых одновременно в этом образце. В качестве другого примера, иногда присутствие соединения может изменить реакцию теста на другие соединения («матричное вмешательство»); когда лаборатория обнаружит это, она соответственно увеличит свои пределы отчетности.

Я ищу практический способ оценки всей дисперсионно-ковариационной матрицы для таких наборов данных, особенно когда многие из соединений подвергаются цензуре более чем на 50%, что часто имеет место. Традиционная модель распределения состоит в том, что логарифмы (истинных) концентраций распределены по нескольким нормам, и это, по-видимому, хорошо подходит на практике, поэтому решение для этой ситуации было бы полезно.

(Под «практическим» я подразумеваю метод, который можно надежно кодировать, по крайней мере, в одной общедоступной программной среде, такой как R, Python, SAS и т. Д., Способом, который выполняется достаточно быстро для поддержки итеративных пересчетов, таких как многократное вменение, и который достаточно стабилен [именно поэтому я неохотно исследую реализацию BUGS, хотя байесовские решения в целом приветствуются].)

Заранее большое спасибо за ваши мысли по этому вопросу.

Whuber
источник
Именно поэтому я правильно понимаю проблему цензуры: когда вы разбавляете образец, концентрация соединения падает настолько низко, что тестирующий прибор может не обнаружить его присутствие. Это точная перефразировка проблемы цензуры?
Да, это правильно: разбавление с коэффициентом D также увеличивает все пределы обнаружения с коэффициентом D. (Проблему матричных помех количественно оценить сложнее, а общая ситуация чрезвычайно сложна. Чтобы упростить это, обычная модель состоит в том, что набор тестов для одной выборки дает вектор (x [1], ..., x [k ]) где x [i] - либо действительные числа, либо интервалы действительных значений, обычно с левой конечной точкой в ​​-infinity; интервал идентифицирует набор, в котором предполагается, что истинное значение лежит.)
whuber
Почему пределы обнаружения возрастают? Не являются ли они особенностью испытательного прибора, а не образца, который тестируется?
В качестве примера предположим, что предел обнаружения прибора составляет 1 микрограмм на литр (мкг / л). Образец разбавляется 10: 1 (с большой точностью, поэтому мы не беспокоимся об ошибке здесь), и прибор читает «<1»; то есть не обнаруживаемый, для разведенного образца. Лаборатория делает вывод, что концентрация в образце составляет менее 10 * 1 = 10 мкг / л, и сообщает об этом как таковую; то есть как <10.
whuber
1
@amoeba Я вижу, я должен был объяснить эти вещи в самом вопросе. Ответы: PCA; размерность будет варьироваться от 3 до нескольких сотен; размеры выборки всегда значительно превышают размерность, но уровни цензуры могут быть очень высокими (требуется возможность обрабатывать до 50% и желательно до 95%).
whuber

Ответы:

3

Я не полностью усвоил проблему матричных помех, но здесь есть один подход. Позволять:

будет вектором, который представляет концентрацию всех целевых соединений в неразбавленном образце.Y

Z

dd

Наша модель это:

YN(μ,Σ)

Z=Yd+ϵ

ϵN(0,σ2 я)

Следовательно, следует, что:

Z~N(μd,Σ+σ2 я)

ZеZ(,)

ОτяTчас

Оязнак равноZяя(Zя>τ)+0я(Zяτ)

К

L(О1,,,,ОК,ОК+1,,,,ОN|-)знак равно[Πязнак равно1язнак равноКпр(Zяτ)][Πязнак равноК+1язнак равноNе(Оя|-)]

где

е(Оя|-)знак равноJяеZ(Оя|-)я(Оя>τ)

Таким образом, оценка заключается в использовании либо максимальной вероятности, либо байесовских идей. Я не уверен, насколько податливы вышесказанное, но я надеюсь, что это даст вам некоторые идеи.


источник
Большое спасибо за эту мысль. Действительно, это стандартный и хорошо документированный подход к множественной цензуре. Одна трудность заключается в его неразрешимости: эти интегралы, как известно, трудно вычислить. Здесь также скрывается проблема моделирования: значение d обычно положительно коррелирует с Y , как подразумевается в первом абзаце моего описания.
whuber
2

Другой более эффективный с точки зрения вычислений вариант - подгонка ковариационной матрицы путем сопоставления моментов с использованием модели, которая называется «дихомизированный гауссовский», на самом деле просто модель гауссовой связки.

В недавней статье Macke et al 2010 описывается процедура закрытой формы для подгонки этой модели, которая включает только (цензурированную) эмпирическую ковариационную матрицу и вычисление некоторых двумерных нормальных вероятностей. Та же группа (лаборатория Бетге в MPI Tuebingen) также описала гибридные дискретные / непрерывные гауссовские модели, которые, вероятно, вам здесь нужны (т. Е. Поскольку гауссовые RV не полностью «дихотомизированы» - только те, которые ниже порогового значения).

Критически, это не оценка ML, и я боюсь, что я не знаю, каковы ее свойства смещения.

jpillow
источник
@jp Спасибо: я посмотрю на это. (Это может занять некоторое время ...)
whuber
1

Сколько соединений в вашем образце? (Или насколько велика рассматриваемая ковариационная матрица?).

У Alan Genz есть несколько очень хороших кодов на разных языках (R, Matlab, Fortran; см. Здесь ) для вычисления интегралов многомерных нормальных плотностей по гипер-прямоугольникам (т. Е. Видов интегралов, которые вам нужны для оценки вероятности, как отмечено user28).

Я использовал эти функции («ADAPT» и «QSIMVN») для интегралов примерно до 10–12 измерений, и несколько функций на этой странице объявляют интегралы (и связанные с ними производные, которые могут вам понадобиться) для задач до измерения 100. Я надеваю Не знаю, достаточно ли измерений для ваших целей, но если это так, то это, вероятно, позволит вам найти максимальные вероятностные оценки по градиентному всплытию.

jpillow
источник
Ой, извините - я новичок здесь и не заметил, как давно это было опубликовано - возможно, слишком поздно, чтобы оказать большую помощь!
jpillow
@jp Это постоянная важная проблема, поэтому прошедшее время между вопросом и ответом не имеет большого значения. Спасибо за ответ!
whuber