Estimasi unjuk kerja matriks kovarians untuk data yang disensor berlipat ganda

Analisis kimia terhadap sampel lingkungan sering disensor di bawah ini pada batas pelaporan atau berbagai batas deteksi / kuantisasi. Yang terakhir dapat bervariasi, biasanya sebanding dengan nilai-nilai variabel lain. Sebagai contoh, sampel dengan konsentrasi tinggi dari satu senyawa mungkin perlu diencerkan untuk analisis, menghasilkan inflasi proporsional batas sensor untuk semua senyawa lain yang dianalisis pada waktu yang sama dalam sampel tersebut. Sebagai contoh lain, kadang-kadang keberadaan senyawa dapat mengubah respons tes terhadap senyawa lain ("gangguan matriks"); saat ini terdeteksi oleh laboratorium, maka akan menaikkan batas pelaporannya.

Saya mencari cara praktis untuk memperkirakan seluruh matriks varians-kovarians untuk dataset tersebut, terutama ketika banyak senyawa mengalami lebih dari 50% sensor, yang sering terjadi. Model distribusi konvensional adalah bahwa logaritma konsentrasi (benar) terdistribusi secara multinormal, dan ini tampaknya cocok dalam praktiknya, sehingga solusi untuk situasi ini akan berguna.

(Dengan "praktis" yang saya maksud adalah metode yang dapat dipercaya dikodekan dalam setidaknya satu lingkungan perangkat lunak yang tersedia secara umum seperti R, Python, SAS, dll., Dengan cara yang dijalankan dengan cukup cepat untuk mendukung perhitungan ulang berulang seperti terjadi dalam beberapa imputasi, dan yang cukup stabil [itulah sebabnya saya enggan mengeksplorasi implementasi BUGS, meskipun solusi Bayesian secara umum diterima].)

Banyak terima kasih sebelumnya atas pemikiran Anda tentang masalah ini.

Saya belum sepenuhnya menginternalisasi masalah gangguan matriks tetapi di sini ada satu pendekatan. Membiarkan:

$Y$ menjadi vektor yang mewakili konsentrasi semua senyawa target dalam sampel yang tidak diencerkan.

$Z$ menjadi vektor yang sesuai dalam sampel encer.

menjadi faktor pengenceran yaitu, sampel diencerkan$d$$d$ : 1.

Model kami adalah:

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

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

di mana mewakili kesalahan karena kesalahan pengenceran.$\epsilon \sim N(0,\sigma^2\ I)$

Oleh karena itu, dapat disimpulkan bahwa:

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

Nyatakan distribusi oleh f Z ( . ) Di atas .$Z$$f_Z(.)$

Biarkan menjadi konsentrasi yang diamati dan τ mewakili ambang instrumen tes di bawah ini yang tidak dapat mendeteksi senyawa. Kemudian, untuk senyawa i t h kami memiliki:$O$$\tau$$i^{th}$

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

Tanpa kehilangan sifat umum, biarkan senyawa pertama sedemikian rupa sehingga berada di bawah ambang batas. Maka fungsi kemungkinan dapat ditulis sebagai:$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 |-)}]$

dimana

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

Estimasi kemudian adalah masalah menggunakan kemungkinan maksimum atau ide bayesian. Saya tidak yakin seberapa mudahnya penjelasan di atas, tetapi saya harap ini memberi Anda beberapa ide.

Pilihan lain yang lebih efisien secara komputasi adalah mencocokkan matriks kovarians dengan cara mencocokkan saat menggunakan model yang telah disebut "Gaussian dichomized", benar-benar hanya model kopula Gaussian.

Makalah terbaru dari Macke et al 2010 menjelaskan prosedur bentuk tertutup untuk menyesuaikan model ini yang hanya melibatkan matriks kovarians empiris (disensor) dan perhitungan beberapa probabilitas normal bivariat. Kelompok yang sama (laboratorium Bethge di MPI Tuebingen) juga menggambarkan model Gaussian diskrit / kontinyu hibrida yang mungkin Anda inginkan di sini (yaitu, karena RV Gaussian tidak sepenuhnya "dikotomisasi" - hanya yang di bawah ambang batas).

Secara kritis, ini bukan estimator ML, dan saya khawatir saya tidak tahu apa sifat biasnya.

Berapa banyak senyawa dalam sampel Anda? (Atau, seberapa besar matriks kovarian yang dimaksud?).

Alan Genz memiliki beberapa kode yang sangat bagus dalam berbagai bahasa (R, Matlab, Fortran; lihat di sini ) untuk menghitung integral dari kepadatan normal multivarian atas hiper-persegi panjang (yaitu, jenis integral yang Anda butuhkan untuk mengevaluasi kemungkinan, seperti dicatat oleh pengguna28).

Saya telah menggunakan fungsi-fungsi ini ("ADAPT" dan "QSIMVN") untuk integral hingga sekitar 10-12 dimensi, dan beberapa fungsi pada halaman tersebut mengiklankan integral (dan turunan terkait yang mungkin Anda perlukan) untuk masalah hingga dimensi 100. Saya tidak tidak tahu apakah itu dimensi yang cukup untuk keperluan Anda, tetapi jika demikian mungkin bisa memungkinkan Anda untuk menemukan perkiraan kemungkinan maksimum dengan kenaikan gradien.