GLM binomial negatif vs. transformasi log untuk data jumlah: peningkatan tingkat kesalahan Tipe I

18

Beberapa dari Anda mungkin telah membaca makalah yang bagus ini:

O'Hara RB, Kotze DJ (2010) Jangan log-transform data hitungan. Metode dalam Ekologi dan Evolusi 1: 118-122. klick .

Dalam bidang penelitian saya (ekotoksikologi), kita sedang berhadapan dengan eksperimen yang direplikasi dengan buruk dan GLM tidak banyak digunakan. Jadi saya melakukan simulasi serupa dengan O'Hara & Kotze (2010), tetapi meniru data ekotoksikologis.

Simulasi daya :

Saya mensimulasikan data dari desain faktorial dengan satu kelompok kontrol ( ) dan 5 kelompok perlakuan ( ). Kelimpahan dalam perawatan 1 identik dengan kontrol ( ), kelimpahan dalam perawatan 2-5 adalah setengah dari kelimpahan dalam kontrol ( ). Untuk simulasi saya memvariasikan ukuran sampel (3,6,9,12) dan kelimpahan pada kelompok kontrol (2, 4, 8, ..., 1024). Kelimpahan diambil dari distribusi binomial negatif dengan parameter dispersi tetap ( ). 100 set data dihasilkan dan dianalisis menggunakan GLM binomial negatif dan data transformasi log Gaussian GLM +.μ 1 - 5 μ 1 = μ c μ 2 - 5 = 0,5 μ c θ = 3,91μcμ1-5μ1=μcμ2-5=0,5μcθ=3.91

Hasilnya seperti yang diharapkan: GLM memiliki kekuatan yang lebih besar, terutama ketika tidak banyak sampel hewan. masukkan deskripsi gambar di sini Kode ada di sini.

Kesalahan Tipe I :

Selanjutnya saya melihat kesalahan tipe satu. Simulasi dilakukan seperti di atas, namun semua grup memiliki kelimpahan yang sama ( ).μc=μ1-5

Namun, hasilnya tidak seperti yang diharapkan: masukkan deskripsi gambar di sini GLM binomial negatif menunjukkan kesalahan Tipe-I yang lebih besar dibandingkan dengan transformasi LM +. Seperti yang diharapkan, perbedaan menghilang dengan meningkatnya ukuran sampel. Kode ada di sini.

Pertanyaan:

Mengapa ada peningkatan Tipe-I Kesalahan dibandingkan dengan transformasi lm +?

Jika kita memiliki data yang buruk (ukuran sampel kecil, kelimpahan rendah (banyak nol)), haruskah kita menggunakan transformasi lm +? Ukuran sampel yang kecil (2-4 per perawatan) adalah tipikal untuk eksperimen semacam itu dan tidak dapat ditingkatkan dengan mudah.

Meskipun, neg. tempat sampah. GLM dapat dibenarkan karena sesuai untuk data ini, transformasi lm + dapat mencegah kita dari kesalahan tipe 1.

EDI
sumber
1
Bukan jawaban untuk pertanyaan utama Anda, tetapi sesuatu untuk dicatat oleh pembaca: kecuali jika Anda membuat kesalahan tipe I yang sebenarnya setara untuk kedua prosedur, membandingkan daya tidak masuk akal; Saya selalu dapat membuat daya lebih tinggi untuk yang lebih rendah (dalam hal ini mengambil log dan cocok dengan yang normal) dengan mengangkat kesalahan tipe I. Di sisi lain, jika Anda menentukan situasi tertentu (ukuran sampel, kelimpahan), Anda bisa mendapatkan tingkat kesalahan tipe I (misalnya dengan simulasi), dan menghitung tingkat nominal yang akan diuji untuk mencapai tingkat kesalahan tipe I yang diinginkan. , jadi kekuatan mereka menjadi sebanding.
Glen_b -Reinstate Monica
Apakah nilai sumbu y di plot Anda dirata-rata di 100 set data?
shadowtalker
Saya harus mengklarifikasi komentar saya: dalam kasus di mana statistik secara inheren diskrit Anda tidak memiliki kontrol sempurna tingkat kesalahan tipe I, tetapi Anda umumnya dapat membuat tingkat kesalahan tipe I cukup dekat. Dalam situasi di mana Anda tidak bisa membuat mereka cukup dekat untuk dapat dibandingkan, satu-satunya cara untuk membuat mereka dapat dibandingkan adalah dengan tes acak.
Glen_b -Reinstate Monica
@ssdecontrol: Tidak, itu hanya proporsi kumpulan data (dari 100) di mana p <α
EDi
1
n

Jawaban:

17

Ini adalah masalah yang sangat menarik. Saya meninjau kode Anda dan tidak dapat menemukan kesalahan ketik yang jelas.

θθdrop1

Sebagian besar tes untuk model linier tidak mengharuskan Anda untuk menghitung ulang model di bawah hipotesis nol. Ini karena Anda dapat menghitung kemiringan geometrik (skor tes) dan memperkirakan lebar (Wald test) menggunakan estimasi parameter dan estimasi kovarians di bawah hipotesis alternatif saja.

Karena binomial negatif tidak linier, saya pikir Anda harus menyesuaikan model nol.

EDIT:

Saya mengedit kode dan mendapatkan yang berikut: masukkan deskripsi gambar di sini

Kode yang diedit di sini: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

AdamO
sumber
Tapi saya pikir itu drop1() tidak cocok secara internal dengan model nol ...
Ben Bolker
4
glm.nbθdrop1logLikgetS3method('logLik', 'negbin'
ingin memberi +1 lagi tetapi saya tidak bisa. Bagus.
Ben Bolker
Terima kasih! Saya hanya melihat kode keduanya drop1()dan lrtest(). Anda benar, drop1.glmpenggunaan glm.fityang memberikan penyimpangan yang salah. Tidak menyadari bahwa kita tidak dapat menggunakan drop1()dengan glm.nb()!
EDi
Jadi skor tipikal dan tes Wald tidak valid dalam model binomial negatif?
shadowtalker
8

Makalah O'Hara dan Kotze (Metode dalam Ekologi dan Evolusi 1: 118-122) bukanlah titik awal yang baik untuk diskusi. Kekhawatiran saya yang paling serius adalah klaim dalam poin 4 ringkasan:

Kami menemukan bahwa transformasi berkinerja buruk, kecuali. . .. The quasi-Poisson dan model binomial negatif ... [menunjukkan] sedikit bias.

λθλ

λ

Kode R berikut mengilustrasikan poin:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

Atau coba

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

Skala estimasi parameter sangat penting!

λ

Perhatikan bahwa diagnostik standar berfungsi lebih baik pada skala log (x + c). Pilihan c mungkin tidak terlalu penting; sering 0,5 atau 1,0 masuk akal. Juga merupakan titik awal yang lebih baik untuk menyelidiki transformasi Box-Cox, atau varian Yeo-Johnson dari Box-Cox. [Yeo, I. dan Johnson, R. (2000)]. Lihat lebih lanjut halaman bantuan untuk powerTransform () dalam paket mobil R. Paket gamls R memungkinkan untuk mencocokkan jenis binomial negatif I (varietas umum) atau II, atau distribusi lain yang memodelkan dispersi serta rata-rata, dengan tautan transformasi daya sebesar 0 (= log, yaitu, tautan log) atau lebih . Fits mungkin tidak selalu menyatu.

Contoh: Deaths vs Base Damage Data untuk badai Atlantik bernama yang mencapai daratan AS. Data tersedia (nama hurricNamed ) dari rilis terbaru dari paket DAAG untuk R. Halaman bantuan untuk data memiliki detail.

Loglinear yang kuat vs fit binomial negatif

Grafik membandingkan garis pas yang diperoleh dengan menggunakan model linear kuat, dengan kurva diperoleh dengan mengubah kecocokan binomial negatif dengan tautan log ke skala log (hitung + 1) yang digunakan untuk sumbu y pada grafik. (Perhatikan bahwa kita harus menggunakan sesuatu yang mirip dengan skala log (hitung + c), dengan positif c, untuk menunjukkan titik dan "garis" yang pas dari kecocokan binomial negatif pada grafik yang sama.) Perhatikan bias besar yang terbukti untuk kecocokan binomial negatif pada skala log. Model linear yang kuat kurang bias pada skala ini, jika kita mengasumsikan distribusi binomial negatif untuk perhitungan. Model linier yang cocok akan tidak bias di bawah asumsi teori normal klasik. Saya menemukan bias yang mengejutkan ketika saya pertama kali membuat apa yang pada dasarnya adalah grafik di atas! Kurva akan lebih cocok dengan data, tetapi perbedaannya ada dalam batas-batas standar variabilitas statistik yang biasa. Model linier yang kuat mampu melakukan pekerjaan yang buruk untuk diperhitungkan pada skala rendah.

Catatan --- Studi dengan Data RNA-Seq: Perbandingan dua gaya model telah menarik untuk analisis data jumlah dari eksperimen ekspresi gen. Makalah berikut membandingkan penggunaan model linear yang kuat, bekerja dengan log (hitung +1), dengan penggunaan cocok binomial negatif (seperti pada tepi paket BioconductorR ). Sebagian besar hitungan, dalam aplikasi RNA-Seq yang terutama dalam pikiran, cukup besar sehingga model log-linear yang ditimbang sesuai sangat cocok untuk bekerja.

Hukum, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: bobot presisi membuka alat analisis model linier untuk jumlah baca RNA-seq. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB juga makalah terbaru:

Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, Wrobel N, Gharbi K, Simpson GG, Owen-Hughes T, Blaxter M, Barton GJ (2016). Berapa banyak ulangan biologis yang diperlukan dalam percobaan RNA-seq dan alat ekspresi diferensial mana yang harus Anda gunakan? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Sangat menarik bahwa model linier cocok menggunakan paket limma (seperti edgeR , dari kelompok WEHI) berdiri sangat baik (dalam arti menunjukkan sedikit bukti bias), relatif terhadap hasil dengan banyak ulangan, karena jumlah ulangan adalah berkurang.

Kode R untuk grafik di atas:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))    
John Maindonald
sumber
2
Terima kasih atas komentar Anda, Tn. Maindonald. Dalam dua tahun terakhir ada juga beberapa makalah lagi (lebih fokus pada pengujian hipotesis, kemudian bias): Ives 2015, Warton et al 2016, Szöcs 2015.
EDi
mungkin itu adalah titik awal yang baik untuk diskusi, bahkan jika titik khusus ini bermasalah? (Saya berpendapat lebih umum bahwa ini adalah alasan untuk tidak terlalu fokus pada bias, melainkan untuk mempertimbangkan sesuatu seperti RMSE ... [disclaimer, saya belum membaca ulang makalah ini belakangan ini, dan saya hanya membaca abstrak dari makalah Warton ...]
Ben Bolker
1
Poin Warton et al (2016), bahwa properti data harus menjadi dasar pilihan, sangat penting. Plot kuantil-kuantil adalah cara yang baik untuk membandingkan detail kecocokan. Khususnya kecocokan pada satu atau lain atau kedua ekstrem dapat menjadi penting untuk beberapa aplikasi. Model zero-inflated atau hurdle dapat menjadi penyempurnaan yang efektif untuk mendapatkan zero counts dengan benar. Pada bagian paling atas, salah satu model yang sedang dibahas dapat sangat dikompromikan. Warton et al do, patut dipuji, memiliki satu contoh. Saya ingin melihat perbandingan di berbagai set data ekologi.
John Maindonald
Tetapi bukankah dalam set data ekologis, spesies di bagian bawah (= spesies langka) menarik? Seharusnya tidak terlalu sulit untuk mengkompilasi beberapa dataset ekologis dan membandingkan ...
EDi
Sebenarnya, untuk low end dari kategori kerusakan itulah model binomial negatif tampaknya, untuk data kematian badai, menjadi paling tidak memuaskan. Paket gamls R memiliki fungsi yang memudahkan untuk membandingkan centile dari distribusi yang pas dengan centile data:
John Maindonald
6

Posting asli mencerminkan makalah Tony Ives: Ives (2015) . Jelas bahwa pengujian signifikansi memberikan hasil yang berbeda dengan estimasi parameter.

John Maindonald menjelaskan mengapa perkiraan itu bias, tetapi ketidaktahuannya tentang latar belakang itu menjengkelkan - dia mengkritik kita karena menunjukkan bahwa metode yang kita semua sepakat cacat cacat. Banyak ahli ekologi melakukan transformasi log secara membabi buta, dan kami berusaha menunjukkan masalah dengan melakukan itu.

Ada diskusi yang lebih bernuansa di sini: Warton (2016)

Ives, AR (2015), Untuk menguji signifikansi koefisien regresi, silakan dan data log-transform count. Metode Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. dan Ives, AR (2016), Tiga poin untuk dipertimbangkan ketika memilih tes LM atau GLM untuk data jumlah. Metode Ecol Evol. doi: 10.1111 / 2041-210X.12552

Bob O'Hara
sumber
Selamat Datang di CV. Meskipun membantu, respons ini sebagian besar merupakan jawaban jenis "hanya tautan". Tautan berubah dan tidak terhubung. Akan lebih bermanfaat bagi CV jika Anda menguraikan poin-poin utama di masing-masing.
Mike Hunter
Terima kasih atas tanggapannya. Saya pikir makalah Warton et al. koin keadaan diskusi saat ini.
EDi
Terima kasih, sama-sama! Saya telah mengambil kebebasan untuk menambahkan referensi secara penuh.
Scortchi
1
Tolong jelaskan poin-poin utama yang dibuat dalam referensi baru, dan di mana masuk akal untuk, juga menghubungkannya kembali ke pertanyaan awal. Ini adalah kontribusi yang berharga tetapi saat ini lebih dekat ke komentar pada jawaban lain daripada jawaban untuk pertanyaan (yang seharusnya memberikan konteks untuk tautan , misalnya). Beberapa kalimat tambahan konteks akan membantu posting secara substansial.
Glen_b -Reinstate Monica
3
Secara khusus, komentar saya membahas poin 4 dalam makalah O'Hara dan Kotze: "Kami menemukan bahwa transformasi berkinerja buruk, kecuali ... Model quasi-Poisson dan negatif binomial ... [menunjukkan] sedikit bias." Simulasi adalah komentar pada perbandingan antara rata-rata yang diharapkan pada skala y (hitungan), versus rata-rata yang diharapkan pada skala log (y + c), untuk distribusi kemiringan yang sangat positif, tidak lebih. Parameter binomial negatif lambda tidak bias pada skala y, sedangkan rata-rata log-normal tidak bias (di bawah normalitas pada skala itu) pada skala log (y + c).
John Maindonald