Metode yang baik untuk plot kerapatan variabel non-negatif dalam R?

36
plot(density(rexp(100))

Jelas semua kepadatan di sebelah kiri nol merupakan bias.

Saya ingin meringkas beberapa data untuk non-statistik, dan saya ingin menghindari pertanyaan tentang mengapa data non-negatif memiliki kepadatan di sebelah kiri nol. Plot-plot tersebut untuk pengecekan pengacakan; Saya ingin menunjukkan distribusi variabel berdasarkan kelompok perlakuan dan kontrol. Distribusi seringkali eksponensial-ish. Histogram rumit karena berbagai alasan.

Pencarian google cepat memberi saya pekerjaan oleh ahli statistik pada kernel non-negatif, misalnya: ini .

Tetapi apakah semua itu telah diterapkan dalam R? Dari metode yang diterapkan, apakah ada di antara mereka yang "terbaik" dalam beberapa cara untuk statistik deskriptif?

EDIT: bahkan jika fromperintah tersebut dapat menyelesaikan masalah saya saat ini, akan menyenangkan untuk mengetahui apakah ada orang yang telah mengimplementasikan kernel berdasarkan literatur pada estimasi kepadatan non-negatif

generic_user
sumber
3
Bukan apa yang Anda minta, tapi saya tidak akan menerapkan estimasi kepadatan kernel untuk sesuatu yang seharusnya eksponensial, terutama untuk presentasi kepada audiens non-statistik. Saya akan menggunakan plot kuantil-kuantil dan menjelaskan bahwa plot harus lurus jika distribusinya eksponensial.
Nick Cox
6
plot(density(rexp(100), from=0))?
Stéphane Laurent
4
Satu hal yang kadang-kadang saya lakukan dengan cukup sukses adalah mendapatkan kde pada log, dan kemudian mengubah estimasi kepadatan (tidak melupakan Jacobian). Kemungkinan lain adalah menggunakan estimasi kepadatan log-spline yang diatur sehingga ia tahu tentang batas tersebut.
Glen_b -Reinstate Monica
1
Saya membahas metode transformasi yang disebutkan oleh @Glen_b di stata-journal.com/sjpdf.html?articlenum=gr0003 (lihat hal.76-78). Nol mungkin diakomodasi dengan menggunakan log (x +1) daripada mencatat dan memodifikasi Jacobian.
Nick Cox

Jawaban:

21

Salah satu solusi, yang dipinjam dari pendekatan untuk pembobotan statistik spasial, adalah memotong kepadatan di sebelah kiri pada nol tetapi untuk menambah data yang paling dekat dengan nol. Idenya adalah bahwa setiap nilai adalah "menyebar" ke dalam kernel dari total area unit yang berpusat di x ; setiap bagian dari kernel yang akan masuk ke wilayah negatif dihapus dan kernel dinormalkan kembali ke area unit.xx

Misalnya, dengan kernel Gaussian , bobot renormalisasi adalahKh(y,x)=exp(12((yx)/h)2)/2π

w(x)=1/0K(y,x)dy=11Φx,h(0)

di mana adalah fungsi distribusi kumulatif dari varian Normal rata-rata x dan standar deviasi h . Formula yang sebanding tersedia untuk kernel lainnya.Φxh

Ini lebih sederhana - dan jauh lebih cepat dalam perhitungan - daripada mencoba mempersempit bandwidth di dekat . Sulit untuk meresepkan bagaimana bandwidth harus diubah mendekati 0 . Namun demikian, metode ini juga ad hoc : masih akan ada beberapa bias dekat 0 . Tampaknya berfungsi lebih baik daripada perkiraan kepadatan default. Berikut ini perbandingan menggunakan dataset largish:000

Angka

Biru menunjukkan kerapatan default sedangkan merah menunjukkan kerapatan yang disesuaikan untuk tepi pada . Distribusi mendasar yang sebenarnya dilacak sebagai garis putus-putus untuk referensi.0


Kode r

The densityfungsi dalam Rakan mengeluh bahwa jumlah bobot tidak kesatuan, karena ingin integral atas semua bilangan real menjadi kesatuan, sedangkan pendekatan ini membuat terpisahkan atas angka positif sama untuk persatuan. Sebagai cek, integral yang terakhir diperkirakan sebagai jumlah Riemann.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
whuber
sumber
21

Alternatifnya adalah pendekatan Kooperberg dan rekan, berdasarkan pada estimasi kepadatan menggunakan splines untuk memperkirakan kepadatan log data. Saya akan menunjukkan contoh menggunakan data dari jawaban @ whuber, yang akan memungkinkan untuk perbandingan pendekatan.

set.seed(17)
x <- rexp(1000)

Anda perlu menginstal paket logspline untuk ini; instal jika tidak:

install.packages("logspline")

Memuat paket dan memperkirakan kepadatan menggunakan logspline()fungsi:

require("logspline")
m <- logspline(x)

Berikut ini, saya menganggap bahwa objek ddari jawaban @ whuber ada di ruang kerja.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

Plot yang dihasilkan ditunjukkan di bawah ini, dengan kepadatan logspline ditunjukkan oleh garis merah

Kerapatan bawaan, terpotong, dan logspline

Selain itu, dukungan untuk kepadatan dapat ditentukan melalui argumen lbounddan ubound. Jika kita ingin mengasumsikan bahwa kerapatan adalah 0 di sebelah kiri 0 dan ada diskontinuitas pada 0, kita dapat menggunakan lbound = 0dalam panggilan ke logspline(), misalnya

m2 <- logspline(x, lbound = 0)

Menghasilkan estimasi kepadatan berikut (ditampilkan di sini dengan mpas logspline asli karena gambar sebelumnya sudah mulai sibuk).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

Plot yang dihasilkan ditunjukkan di bawah ini

Perbandingan perkiraan kepadatan logspline dengan dan tanpa batas bawah pada dukungan

xx=0x

Pasang kembali Monica - G. Simpson
sumber
1
01
@whuber Pertanyaan bagus. Saya sendiri baru saja menemukan pendekatan ini. Saya menduga pertanyaan yang bagus untuk ditanyakan di sini adalah, karena metode terpotong dan logspline hanya perkiraan kepadatan sebenarnya, apakah perbedaan dalam fit signifikan, secara statistik? Saya tidak yakin persis mengapa ia melakukannya dengan sangat baik di nol. Saya juga menghargai mengetahui mengapa.
Reinstate Monica - G. Simpson
@ GavinSimpson, Terima kasih atas jawaban yang bagus ini. Bisakah Anda mereproduksi plot terakhir dengan versi terbaru logspline? Bagi saya, kepadatan keduanya, versi terikat dan tidak terikat pergi ke nol pada x = 0.
cel
4

Untuk membandingkan distribusi berdasarkan kelompok (yang Anda katakan adalah tujuan dalam salah satu komentar Anda) mengapa tidak sesuatu yang lebih sederhana? Plot kotak paralel bekerja dengan baik jika N besar; plot strip paralel berfungsi jika N kecil (dan keduanya menunjukkan outlier dengan baik, yang menurut Anda merupakan masalah dalam data Anda).

Peter Flom - Pasang kembali Monica
sumber
1
Ya, terima kasih, itu berhasil. Tapi saya suka plot kepadatan. Mereka menunjukkan lebih banyak tentang data daripada plot kotak lakukan. Saya kira saya agak terkejut bahwa sepertinya belum ada yang diterapkan. Mungkin saya akan menerapkan salah satu dari hal ini sendiri suatu hari nanti. Orang mungkin akan menganggapnya berguna.
generic_user
1
Saya suka plot kerapatan juga; tetapi Anda harus mempertimbangkan audiens Anda.
Peter Flom - Pasang kembali Monica
1
Harus setuju dengan @PeterFlom untuk yang ini. Jangan terlalu rumit jika audiens Anda tidak berpengetahuan secara statistik. Anda juga dapat melakukan plot kotak komparatif / paralel dengan hamparan plot kupu-kupu di atasnya. Dengan cara itu ringkasan kotak-plot terlihat serta semua data.
doug.numbers
Saran bahwa orang yang berbeda memahami plot agregat secara berbeda tentu saja benar. Meskipun memahami apa plot kepadatan (dan memahami bahwa itu bukan probabilitas), saya tidak memiliki pemahaman tentang apa "plot kotak paralel" mungkin. Ini menyarankan plot koordinat paralel tetapi saya menduga itu tidak benar.
DWin
2

Sebagai komentar Stéphane Anda dapat menggunakan from = 0dan, selain itu, Anda dapat mewakili nilai Anda di bawah kurva kepadatan denganrug (x)

Aghila
sumber
4
Koreksi saya jika saya salah tetapi from=0tampak seolah-olah hanya menekan merencanakan nilai di bawah 0; itu tidak benar perhitungan untuk fakta bahwa beberapa distribusi telah dioleskan di bawah 0.
Nick Cox
1
Itu betul. Menggunakan fromperintah menghasilkan plot yang tampaknya memiliki puncak tepat nol. Tetapi jika Anda melihat histogram dengan nampan yang terus menerus lebih kecil, banyak data akan menunjukkan puncak AT nol. Ini fromhanyalah trik grafis.
generic_user
@NickCox Saya tidak yakin tetapi saya tidak berpikir from=0menekan apa pun. Itu hanya memulai "grid" di nol.
Stéphane Laurent
Perbedaannya adalah apakah estimasi kepadatan bukan nol untuk nilai negatif, bukan apakah diplot atau tidak. Peneliti dapat memutuskan untuk tidak khawatir tentang ini jika yang mereka inginkan hanyalah visualisasi.
Nick Cox
@NickCox Perintah density(rexp(100), from=0)ini tidak ada hubungannya dengan grafik
Stéphane Laurent