Estimasi unjuk kerja matriks kovarians untuk data yang disensor berlipat ganda

22

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.

whuber
sumber
Persis seperti yang saya pahami dengan benar masalah sensor: Ketika Anda mencairkan sampel, konsentrasi suatu senyawa turun sangat rendah sehingga instrumen uji bisa gagal mendeteksi keberadaannya. Apakah itu pengungkapan ulang masalah sensor yang akurat?
Ya, itu benar: pengenceran oleh faktor D meningkatkan semua batas deteksi oleh faktor D juga. (Masalah gangguan matriks lebih sulit untuk dikuantifikasi dan situasi umumnya sangat kompleks. Untuk menyederhanakan ini, model konvensional adalah bahwa serangkaian pengujian pada satu sampel menghasilkan vektor (x [1], ..., x [k ]) di mana x [i] adalah bilangan real atau interval real, biasanya dengan titik akhir kiri pada-infinity; suatu interval mengidentifikasi set di mana nilai sebenarnya diasumsikan terletak.)
whuber
Mengapa batas deteksi naik? Apakah itu bukan fitur instrumen tes daripada sampel yang sedang diuji?
Sebagai contoh, misalkan batas deteksi instrumen adalah 1 mikrogram per Liter (ug / L). Sampel diencerkan 10: 1 (dengan presisi tinggi, jadi kami tidak khawatir tentang kesalahan di sini) dan instrumen bertuliskan "<1"; yaitu, tidak terdeteksi, untuk sampel encer. Laboratorium menyimpulkan bahwa konsentrasi dalam sampel kurang dari 10 * 1 = 10 ug / L dan melaporkannya; yaitu, sebagai "<10".
whuber
1
@amoeba saya mengerti saya harus menjelaskan hal-hal itu dalam pertanyaan itu sendiri. Jawabannya adalah: PCA; dimensi akan bervariasi dari 3 hingga beberapa ratus; ukuran sampel selalu jauh melebihi dimensi tetapi tingkat sensor mungkin sangat tinggi (mampu menangani hingga 50% diperlukan dan hingga 95% diinginkan).
whuber

Jawaban:

3

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 ddd : 1.

Model kami adalah:

YN(μ,Σ)

Z=Yd+ϵ

di mana mewakili kesalahan karena kesalahan pengenceran.ϵN(0,σ2 I)

Oleh karena itu, dapat disimpulkan bahwa:

ZN(μd,Σ+σ2 I)

Nyatakan distribusi oleh f Z ( . ) Di atas .ZfZ(.)

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τith

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

Tanpa kehilangan sifat umum, biarkan senyawa pertama sedemikian rupa sehingga berada di bawah ambang batas. Maka fungsi kemungkinan dapat ditulis sebagai:k

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

dimana

f(HAIsaya|-)=jsayafZ(HAIsaya|-)saya(HAIsaya>τ)

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.


sumber
Terima kasih banyak atas pemikiran ini. Memang, ini adalah pendekatan standar dan didokumentasikan dengan baik untuk sensor ganda. Satu kesulitan terletak pada sifatnya yang tidak bisa dipraktikkan: integral-integral itu terkenal sulit untuk dihitung. Ada masalah pemodelan yang mengintai di sini, juga: nilai d biasanya berkorelasi positif dengan Y , seperti yang tersirat pada paragraf pertama deskripsi saya.
whuber
2

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.

jpillow
sumber
@ jp Terima kasih: Saya akan melihat ini. (Mungkin perlu waktu ...)
whuber
1

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.

jpillow
sumber
Oh, maaf — saya baru di sini dan tidak menyadari sudah berapa lama ini diposting — mungkin sudah terlambat untuk banyak membantu!
jpillow
@ jp Ini adalah masalah penting yang sedang berlangsung, jadi waktu yang berlalu antara pertanyaan dan jawaban adalah konsekuensi kecil. Terima kasih untuk balasannya!
whuber