Ratakan rangkaian waktu melingkar / periodik

9

Saya memiliki data kecelakaan kendaraan bermotor menurut jam sehari. Seperti yang Anda harapkan, mereka tinggi di tengah hari dan puncaknya pada jam sibuk. Geom_density default ggplot2 memuluskannya dengan baik

Subset dari data, untuk crash terkait drive-minum, tinggi di kedua ujung hari (malam hari dan pagi hari) dan tertinggi di ekstrem. Tapi geom_density default ggplot2 masih menurun di sisi kanan ekstrim.

Apa yang harus dilakukan? Tujuannya hanyalah visualisasi - tidak perlu (ada?) Untuk analisis statistik yang kuat.

Imgur

x <- structure(list(hour = c(14, 1, 1, 9, 2, 11, 20, 5, 22, 13, 21, 
                        2, 22, 10, 18, 0, 2, 1, 2, 15, 20, 23, 17, 3, 3, 16, 19, 23, 
                        3, 4, 4, 22, 2, 21, 20, 1, 19, 18, 17, 23, 23, 3, 11, 4, 23, 
                        4, 7, 2, 3, 19, 2, 18, 3, 17, 1, 9, 19, 23, 9, 6, 2, 1, 23, 21, 
                        22, 22, 22, 20, 1, 21, 6, 2, 22, 23, 19, 17, 19, 3, 22, 21, 4, 
                        10, 17, 23, 3, 7, 19, 16, 2, 23, 4, 5, 1, 20, 7, 21, 19, 2, 21)
               , count = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L))
          , .Names = c("hour", "count")
          , row.names = c(8L, 9L, 10L, 29L, 33L, 48L, 51L, 55L, 69L, 72L, 97L, 108L, 113L, 
                          118L, 126L, 140L, 150L, 171L, 177L, 184L, 202L, 230L, 236L, 240L, 
                          242L, 261L, 262L, 280L, 284L, 286L, 287L, 301L, 318L, 322L, 372L, 
                          380L, 385L, 432L, 448L, 462L, 463L, 495L, 539L, 557L, 563L, 566L, 
                          570L, 577L, 599L, 605L, 609L, 615L, 617L, 624L, 663L, 673L, 679L, 
                          682L, 707L, 730L, 733L, 746L, 754L, 757L, 762L, 781L, 793L, 815L, 
                          817L, 823L, 826L, 856L, 864L, 869L, 877L, 895L, 899L, 918L, 929L, 
                          937L, 962L, 963L, 978L, 980L, 981L, 995L, 1004L, 1005L, 1007L, 
                          1008L, 1012L, 1015L, 1020L, 1027L, 1055L, 1060L, 1078L, 1079L, 
                          1084L)
          , class = "data.frame")

ggplot(x, aes(hour)) + 
  geom_bar(binwidth = 1, position = "dodge", fill = "grey") +
  geom_density() + 
  aes(y = ..count..) +
  scale_x_continuous(breaks = seq(0,24,4))

Senang bagi siapa pun dengan kosa kata statistik yang lebih baik untuk mengedit pertanyaan ini, terutama judul dan tag.

r pdf kernel-smoothing splines nacnudus
sumber

Jawaban:

6

Untuk membuat kelancaran berkala (pada platform apa pun), cukup tambahkan datanya sendiri, haluskan daftar yang lebih panjang, dan potong ujungnya.

Berikut ini Rilustrasi:

y <- sqrt(table(factor(x[,"hour"], levels=0:23)))
y <- c(y,y,y)
x.mid <- 1:24; offset <- 24
plot(x.mid-1, y[x.mid+offset]^2, pch=19, xlab="Hour", ylab="Count")
y.smooth <- lowess(y, f=1/8)
lines(x.mid-1, y.smooth$y[x.mid+offset]^2, lwd=2, col="Blue")

(Karena ini adalah jumlah yang saya pilih untuk menghaluskan akar kuadrat mereka; mereka dikonversi kembali ke jumlah untuk plot.) Rentang dalam lowesstelah menyusut jauh dari standarnya f=2/3karena (a) kita sekarang memproses array tiga kali lebih lama, yang seharusnya menyebabkan kita mengurangi menjadi , dan (b) Saya ingin smooth lokal yang cukup sehingga tidak ada efek titik akhir yang cukup besar muncul di sepertiga tengah.f2/9

Ini telah melakukan pekerjaan yang cukup bagus dengan data ini. Secara khusus, anomali pada jam 0 telah dihaluskan.

Merencanakan

whuber
sumber
Ini menjawab kebutuhan saya akan visualisasi sederhana, tetapi karena ketertarikan, apakah ini sedikit omong kosong? Apakah menggunakan sesuatu dari tautan Nick menghindari efek titik akhir?
nacnudus
1
Ini persis sama dengan metode yang saya gunakan selama lebar jendela dipilih dengan hati-hati, seperti yang dilakukan @whuber. Tetapi perangkat lunak R sudah tersedia untuk melakukan apa yang saya lakukan. (Awalnya saya mendelegasikan tugas untuk menemukannya ke ahli R, tetapi mereka tidak menyadarinya.)
Nick Cox
3
Saya tidak melihatnya sebagai kluge: teknik ini didasarkan pada definisi periodisitas. Ini berfungsi untuk semua kelancaran lokal. (Ini tidak akan berfungsi untuk kelancaran global, tapi itu bukan masalah, karena kebanyakan smoothers global berasal dari metode periodik inheren seperti Fourier series.) @Nick One tidak harus sangat berhati-hati: ketika menggunakan pemul maksimum setengah-lebar , kita hanya perlu menempelkan nilai nilai terakhir dari urutan ke awal dan ke ujung, tetapi tidak ada salahnya memperluas urutan secara konservatif dengan lebih banyak - hanya saja kurang efisien . kk−1k−1
whuber
1
@whuber sangat. Saya hanya menyinggung kebenaran bahwa apa yang Anda tambahkan sebagai salinan depan dan belakang dari data aktual harus konsisten dengan seberapa banyak Anda memuluskan.
Nick Cox
7

Saya tidak menggunakan R secara rutin dan saya tidak pernah menggunakannya ggplot, tetapi ada cerita sederhana di sini, atau jadi saya kira.

Waktu dalam sehari adalah variabel lingkaran atau periodik. Dalam data Anda, Anda memiliki jam 0 (1) 23 yang membungkus, sehingga 23 diikuti oleh 0. Namun, ggplottidak tahu itu, setidaknya dari informasi yang Anda berikan. Sejauh menyangkut mungkin ada nilai-nilai pada -1, -2, dll atau pada 24, 25, dll dan beberapa kemungkinan kemungkinan diperhalus melampaui batas data yang diamati, dan memang di luar batas data yang mungkin.

Ini akan terjadi untuk data utama Anda juga, tetapi tidak begitu terlihat.

Jika Anda menginginkan estimasi kepadatan kernel untuk data tersebut, Anda memerlukan rutin yang cukup pintar untuk menangani variabel periodik atau melingkar seperti itu dengan benar. "Benar" berarti bahwa rutin merapikan pada ruang melingkar, mengakui bahwa 0 mengikuti 23. Dalam beberapa hal perataan distribusi seperti itu lebih mudah daripada kasus biasa, karena tidak ada masalah batas (karena tidak ada batas). Orang lain harus dapat memberi saran tentang fungsi yang digunakan dalam R.

Jenis data ini berada di antara deret waktu periodik dan statistik sirkuler.

Data yang disajikan memiliki 99 observasi. Untuk itu histogram berfungsi dengan baik, walaupun saya dapat melihat bahwa Anda mungkin ingin sedikit memperhalusnya.

masukkan deskripsi gambar di sini

(PEMBARUAN) Ini masalah selera dan penilaian, tetapi saya akan menganggap kurva mulus Anda secara drastis melampaui.

Di sini sebagai sampel adalah estimasi kepadatan biweight. Saya menggunakan program Stata saya sendiri untuk data sirkuler dalam derajat dengan konversi ad hoc 15 * (jam + 0,5) tetapi kepadatan dinyatakan per jam. Sebaliknya ini agak kurang meyakinkan, tetapi Anda dapat menyesuaikan pilihan Anda.

masukkan deskripsi gambar di sini

Nick Cox
sumber
1
Setuju bahwa itu oversmoothed, tapi itulah prinsip yang saya maksudkan. Beberapa googling dari vocab bermanfaat Anda (sirkuler, periodik) mengungkap secara mengejutkan sedikit minat pada masalah semacam ini, tetapi saya akan menunggu sedikit lebih lama bagi siapa saja untuk berpadu dengan saran R.
nacnudus
3
cran.r-project.org/web/packages/circular/circular.pdf
Nick Cox
5

Melakukan Tukey's 4253H, dua kali pada tiga salinan gabungan jumlah mentah dan kemudian mengambil set tengah nilai dihaluskan memberikan banyak gambaran yang sama seperti kerendahan gumpalan pada akar kuadrat dari hitungan.
masukkan deskripsi gambar di sini

Ray Koopman
sumber
2
+1 Saya lebih suka smoothers Tukey dan senang melihat contoh dari satu acara di sini.
whuber
1
Resep yang tepat ini dirancang oleh Paul F. Velleman, tetapi tidak diragukan lagi di bawah bimbingan Tukey. "42" memotong artefak tangga.
Nick Cox
2

Selain itu, dan sebagai alternatif yang lebih kompleks, dari apa yang telah disarankan, Anda mungkin ingin melihat splines periodik. Anda dapat menemukan alat untuk menyesuaikan mereka dalam paket R splinesaand mgcv. Keuntungan yang saya lihat dari pendekatan yang sudah disarankan adalah bahwa Anda dapat menghitung derajat kebebasan yang sesuai, yang tidak jelas dengan metode 'tiga salinan'.

F. Tusell
sumber
1
(+1) Beberapa komentar: Pertama, "tiga salinan" adalah aplikasi tertentu, bukan aturan umum. Kedua, saya percaya perhitungan DF sama mudahnya: jumlah data tetap sama dan satu mengurangi jumlah parameter yang digunakan dalam pemasangan spline.
whuber
@whuber: tidak jelas bagi saya bagaimana melakukan bit terakhir (bagaimana menghitung parameter yang digunakan pas spline jika Anda memasukkannya ke "tiga salinan").
F. Tusell
1
Bagian penyalinan tidak mengubah jumlah data, jadi yang terpenting dalam memperkirakan DF adalah menghitung parameter yang digunakan oleh splines.
whuber
1

Masih pendekatan lain, splines periodik (seperti yang disarankan dalam jawaban oleh F.Tusell), tetapi di sini kami juga menunjukkan implementasi dalam R. Kami akan menggunakan Poisson glm agar sesuai dengan jumlah histogram, menghasilkan histogram berikut dengan halus:

masukkan deskripsi gambar di sini

Kode yang digunakan (dimulai dengan objek data yang xdiberikan dalam pertanyaan):

library(pbs) # basis for periodic spline

x.tab <- with(x, table(factor(hour,levels=as.character(0:23))))
x.df <- data.frame(time=0:23, count=as.vector(x.tab))
mod.hist <- with(x.df, glm(count ~ pbs::pbs(time, df=4, Boundary.knots=c(0,24)), family=poisson))
pred <- predict(mod.hist, type="response", newdata=data.frame(time=0:24))

with(x.df, {plot(time, count,type="h",col="blue", main="Histogram") ; lines(time, pred[1:24], col="red")} )
kjetil b halvorsen
sumber