Saya memiliki 17 tahun (1995 hingga 2011) data sertifikat kematian yang terkait dengan kematian bunuh diri untuk sebuah negara bagian di AS. Ada banyak mitologi di luar sana tentang bunuh diri dan bulan / musim, banyak di antaranya kontradiktif, dan literatur saya. Sudah diulas, saya tidak mendapatkan pemahaman yang jelas tentang metode yang digunakan atau kepercayaan pada hasil.
Jadi saya telah menetapkan untuk melihat apakah saya dapat menentukan apakah bunuh diri lebih atau kurang mungkin terjadi pada bulan tertentu dalam kumpulan data saya. Semua analisis saya dilakukan dalam R.
Jumlah total kasus bunuh diri dalam data adalah 13.909.
Jika Anda melihat tahun dengan bunuh diri paling sedikit, itu terjadi pada 309/365 hari (85%). Jika Anda melihat tahun dengan bunuh diri terbanyak, mereka terjadi pada 339/365 hari (93%).
Jadi ada cukup banyak hari setiap tahun tanpa bunuh diri. Namun, ketika dikumpulkan di seluruh 17 tahun, ada bunuh diri pada setiap hari sepanjang tahun, termasuk 29 Februari (meskipun hanya 5 ketika rata-rata 38).
Cukup menambahkan jumlah bunuh diri pada setiap hari dalam setahun tidak menunjukkan musim yang jelas (bagi saya).
Diagregasi pada tingkat bulanan, rata-rata bunuh diri per bulan berkisar dari:
(m = 65, sd = 7.4, hingga m = 72, sd = 11.1)
Pendekatan pertama saya adalah untuk mengumpulkan kumpulan data berdasarkan bulan untuk semua tahun dan melakukan uji chi-square setelah menghitung probabilitas yang diharapkan untuk hipotesis nol, bahwa tidak ada varians sistematis dalam jumlah bunuh diri berdasarkan bulan. Saya menghitung probabilitas untuk setiap bulan dengan mempertimbangkan jumlah hari (dan menyesuaikan Februari untuk tahun kabisat).
Hasil chi-square menunjukkan tidak ada variasi yang signifikan berdasarkan bulan:
# So does the sample match expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131
Gambar di bawah ini menunjukkan jumlah total per bulan. Garis merah horizontal diposisikan pada nilai yang diharapkan untuk bulan Februari, 30 hari, dan 31 hari masing-masing. Konsisten dengan uji chi-square, tidak ada bulan di luar interval kepercayaan 95% untuk jumlah yang diharapkan.
Saya pikir saya sudah selesai sampai saya mulai menyelidiki data deret waktu. Seperti yang saya bayangkan banyak orang lakukan, saya mulai dengan metode dekomposisi musiman non-parametrik menggunakan stl
fungsi dalam paket statistik.
Untuk membuat data deret waktu, saya mulai dengan data bulanan teragregasi:
suicideByMonthTs <- ts(suicideByMonth$monthlySuicideCount, start=c(1995, 1), end=c(2011, 12), frequency=12)
# Plot the monthly suicide count, note the trend, but seasonality?
plot(suicideByMonthTs, xlab="Year",
ylab="Annual monthly suicides")
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1995 62 47 55 74 71 70 67 69 61 76 68 68
1996 64 69 68 53 72 73 62 63 64 72 55 61
1997 71 61 64 63 60 64 67 50 48 49 59 72
1998 67 54 72 69 78 45 59 53 48 65 64 44
1999 69 64 65 58 73 83 70 73 58 75 71 58
2000 60 54 67 59 54 69 62 60 58 61 68 56
2001 67 60 54 57 51 61 67 63 55 70 54 55
2002 65 68 65 72 79 72 64 70 59 66 63 66
2003 69 50 59 67 73 77 64 66 71 68 59 69
2004 68 61 66 62 69 84 73 62 71 64 59 70
2005 67 53 76 65 77 68 65 60 68 71 60 79
2006 65 54 65 68 69 68 81 64 69 71 67 67
2007 77 63 61 78 73 69 92 68 72 61 65 77
2008 67 73 81 73 66 63 96 71 75 74 81 63
2009 80 68 76 65 82 69 74 88 80 86 78 76
2010 80 77 82 80 77 70 81 89 91 82 71 73
2011 93 64 87 75 101 89 87 78 106 84 64 71
Dan kemudian dilakukan stl()
dekomposisi
# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)
Pada titik ini saya menjadi khawatir karena saya melihat ada komponen musiman dan tren. Setelah banyak riset internet, saya memutuskan untuk mengikuti instruksi Rob Hyndman dan George Athanasopoulos sebagaimana tercantum dalam teks on-line mereka "Peramalan: prinsip dan praktik", khususnya untuk menerapkan model ARIMA musiman.
Saya menggunakan adf.test()
dan kpss.test()
untuk menilai stasioneritas dan mendapatkan hasil yang bertentangan. Mereka berdua menolak hipotesis nol (mencatat bahwa mereka menguji hipotesis sebaliknya).
adfResults <- adf.test(suicideByMonthTs, alternative = "stationary") # The p < .05 value
adfResults
Augmented Dickey-Fuller Test
data: suicideByMonthTs
Dickey-Fuller = -4.5033, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
kpssResults <- kpss.test(suicideByMonthTs)
kpssResults
KPSS Test for Level Stationarity
data: suicideByMonthTs
KPSS Level = 2.9954, Truncation lag parameter = 3, p-value = 0.01
Saya kemudian menggunakan algoritma dalam buku ini untuk melihat apakah saya dapat menentukan jumlah perbedaan yang perlu dilakukan untuk tren dan musim. Saya berakhir dengan nd = 1, ns = 0.
Saya kemudian berlari auto.arima
, yang memilih model yang memiliki tren dan komponen musiman bersama dengan konstanta tipe "drift".
# Extract the best model, it takes time as I've turned off the shortcuts (results differ with it on)
bestFit <- auto.arima(suicideByMonthTs, stepwise=FALSE, approximation=FALSE)
plot(theForecast <- forecast(bestFit, h=12))
theForecast
> summary(bestFit)
Series: suicideByMonthFromMonthTs
ARIMA(0,1,1)(1,0,1)[12] with drift
Coefficients:
ma1 sar1 sma1 drift
-0.9299 0.8930 -0.7728 0.0921
s.e. 0.0278 0.1123 0.1621 0.0700
sigma^2 estimated as 64.95: log likelihood=-709.55
AIC=1429.1 AICc=1429.4 BIC=1445.67
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.2753657 8.01942 6.32144 -1.045278 9.512259 0.707026 0.03813434
Akhirnya, saya melihat residu dari fit dan jika saya memahami ini dengan benar, karena semua nilai berada dalam batas ambang, mereka berperilaku seperti white noise dan dengan demikian modelnya cukup masuk akal. Saya menjalankan tes portmanteau seperti yang dijelaskan dalam teks, yang memiliki nilai ap jauh di atas 0,05, tetapi saya tidak yakin bahwa saya memiliki parameter yang benar.
Acf(residuals(bestFit))
Box.test(residuals(bestFit), lag=12, fitdf=4, type="Ljung")
Box-Ljung test
data: residuals(bestFit)
X-squared = 7.5201, df = 8, p-value = 0.4817
Setelah kembali dan membaca bab tentang pemodelan arima lagi, saya menyadari sekarang bahwa auto.arima
memang memilih untuk memodelkan tren dan musim. Dan saya juga menyadari bahwa peramalan tidak secara khusus analisis yang mungkin harus saya lakukan. Saya ingin tahu apakah bulan tertentu (atau lebih umum waktu dalam setahun) harus ditandai sebagai bulan berisiko tinggi. Tampaknya alat dalam literatur perkiraan sangat relevan, tetapi mungkin bukan yang terbaik untuk pertanyaan saya. Setiap dan semua input sangat dihargai.
Saya memposting tautan ke file csv yang berisi jumlah harian. File terlihat seperti ini:
head(suicideByDay)
date year month day_of_month t count
1 1995-01-01 1995 01 01 1 2
2 1995-01-03 1995 01 03 2 1
3 1995-01-04 1995 01 04 3 3
4 1995-01-05 1995 01 05 4 2
5 1995-01-06 1995 01 06 5 3
6 1995-01-07 1995 01 07 6 2
Hitung adalah jumlah bunuh diri yang terjadi pada hari itu. "t" adalah urutan numerik dari 1 hingga total jumlah hari dalam tabel (5533).
Saya telah mencatat komentar di bawah ini dan memikirkan dua hal yang berkaitan dengan pemodelan bunuh diri dan musim. Pertama, sehubungan dengan pertanyaan saya, bulan hanyalah proksi untuk menandai pergantian musim, saya tidak tertarik pada apakah bulan tertentu berbeda dari yang lain (yang tentu saja merupakan pertanyaan yang menarik, tapi bukan itu yang saya maksudkan untuk menyelidiki). Oleh karena itu, saya pikir masuk akal untuk menyamakan bulan dengan hanya menggunakan 28 hari pertama dari semua bulan. Ketika Anda melakukan ini, Anda mendapatkan kecocokan sedikit lebih buruk, yang saya tafsirkan sebagai lebih banyak bukti terhadap kurangnya musim. Pada output di bawah, fit pertama adalah reproduksi dari jawaban di bawah ini menggunakan bulan dengan jumlah hari sebenarnya, diikuti oleh set data bunuh diriByShortMonthdi mana jumlah bunuh diri dihitung dari 28 hari pertama dari semua bulan. Saya tertarik dengan pendapat orang tentang apakah penyesuaian ini adalah ide yang baik, tidak perlu, atau berbahaya?
> summary(seasonFit)
Call:
glm(formula = count ~ t + days_in_month + cos(2 * pi * t/12) +
sin(2 * pi * t/12), family = "poisson", data = suicideByMonth)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4782 -0.7095 -0.0544 0.6471 3.2236
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8662459 0.3382020 8.475 < 2e-16 ***
t 0.0013711 0.0001444 9.493 < 2e-16 ***
days_in_month 0.0397990 0.0110877 3.589 0.000331 ***
cos(2 * pi * t/12) -0.0299170 0.0120295 -2.487 0.012884 *
sin(2 * pi * t/12) 0.0026999 0.0123930 0.218 0.827541
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 302.67 on 203 degrees of freedom
Residual deviance: 190.37 on 199 degrees of freedom
AIC: 1434.9
Number of Fisher Scoring iterations: 4
> summary(shortSeasonFit)
Call:
glm(formula = shortMonthCount ~ t + cos(2 * pi * t/12) + sin(2 *
pi * t/12), family = "poisson", data = suicideByShortMonth)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.2414 -0.7588 -0.0710 0.7170 3.3074
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.0022084 0.0182211 219.647 <2e-16 ***
t 0.0013738 0.0001501 9.153 <2e-16 ***
cos(2 * pi * t/12) -0.0281767 0.0124693 -2.260 0.0238 *
sin(2 * pi * t/12) 0.0143912 0.0124712 1.154 0.2485
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 295.41 on 203 degrees of freedom
Residual deviance: 205.30 on 200 degrees of freedom
AIC: 1432
Number of Fisher Scoring iterations: 4
Hal kedua yang saya lihat lebih dalam adalah masalah menggunakan bulan sebagai proxy untuk musim. Mungkin indikator musim yang lebih baik adalah jumlah jam siang hari yang diterima suatu daerah. Data ini berasal dari negara bagian utara yang memiliki variasi substansial di siang hari. Di bawah ini adalah grafik siang hari dari tahun 2002.
Ketika saya menggunakan data ini daripada bulan dalam setahun, efeknya masih signifikan, tetapi efeknya sangat, sangat kecil. Penyimpangan residual jauh lebih besar daripada model di atas. Jika siang hari adalah model yang lebih baik untuk musim, dan kecocokannya tidak sebaik, apakah ini lebih banyak bukti efek musiman yang sangat kecil?
> summary(daylightFit)
Call:
glm(formula = aggregatedDailyCount ~ t + daylightMinutes, family = "poisson",
data = aggregatedDailyNoLeap)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.0003 -0.6684 -0.0407 0.5930 3.8269
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.545e+00 4.759e-02 74.493 <2e-16 ***
t -5.230e-05 8.216e-05 -0.637 0.5244
daylightMinutes 1.418e-04 5.720e-05 2.479 0.0132 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 380.22 on 364 degrees of freedom
Residual deviance: 373.01 on 362 degrees of freedom
AIC: 2375
Number of Fisher Scoring iterations: 4
Saya memposting jam siang kalau-kalau ada yang ingin bermain-main dengan ini. Catatan, ini bukan tahun kabisat, jadi jika Anda ingin memasukkan menit untuk tahun kabisat, ekstrapolasi atau ambil data.
[ Sunting untuk menambahkan plot dari jawaban yang dihapus (mudah-mudahan juga tidak keberatan saya memindahkan plot dalam jawaban yang dihapus di sini ke pertanyaan. Svannoy, jika Anda tidak ingin ini ditambahkan, Anda dapat mengembalikannya)]
sumber
Jawaban:
Bagaimana dengan regresi Poisson?
Saya membuat kerangka data yang berisi data Anda, ditambah indeks
t
untuk waktu (dalam bulan) dan variabelmonthdays
untuk jumlah hari di setiap bulan.Jadi terlihat seperti ini:
Sekarang mari kita bandingkan model dengan efek waktu dan efek beberapa hari dengan model di mana kita menambahkan efek bulan:
Berikut ini ringkasan model "kecil":
Anda dapat melihat bahwa kedua variabel memiliki efek marginal yang sangat besar. Sekarang lihat model yang lebih besar:
Yah, tentu saja
monthdays
efeknya menghilang; itu bisa diperkirakan hanya berkat tahun kabisat !! Menyimpannya dalam model (dan memperhitungkan tahun kabisat) memungkinkan untuk menggunakan penyimpangan residual untuk membandingkan kedua model.Sekarang bandingkan dengan model nol:
Jadi, orang dapat mengatakan bahwa ini menunjukkan efek musiman ...
sumber
Tes chi-square adalah pendekatan yang baik sebagai pandangan awal untuk pertanyaan Anda.
The
stl
dekomposisi dapat menyesatkan sebagai alat untuk menguji keberadaan musiman. Prosedur ini berhasil mengembalikan pola musiman yang stabil bahkan jika white noise (sinyal acak tanpa struktur) dilewatkan sebagai input. Coba misalnya:Melihat pesanan yang dipilih oleh prosedur pemilihan model ARIMA otomatis juga agak berisiko karena model ARIMA musiman tidak selalu melibatkan musiman (untuk detail lihat diskusi ini ). Dalam kasus ini, model yang dipilih menghasilkan siklus musiman dan komentar oleh @RichardHardy masuk akal, namun, wawasan lebih lanjut diperlukan untuk menyimpulkan bahwa bunuh diri didorong oleh pola musiman.
Di bawah, saya merangkum beberapa hasil berdasarkan analisis dari seri bulanan yang Anda posting. Ini adalah komponen musiman yang diestimasi berdasarkan model deret waktu struktural dasar:
Komponen serupa diekstraksi menggunakan perangkat lunak TRAMO-SEATS dengan opsi default. Kita dapat melihat bahwa pola musiman yang diperkirakan tidak stabil sepanjang waktu dan, karenanya, tidak mendukung hipotesis pola berulang dalam jumlah bunuh diri di bulan selama periode sampel. Menjalankan perangkat lunak X-13ARIMA-SEATS dengan opsi default, tes gabungan untuk musiman menyimpulkan bahwa musiman yang dapat diidentifikasi tidak ada.
Sunting (lihat jawaban ini dan komentar saya di bawah ini untuk definisi musiman yang dapat diidentifikasi ).
Mengingat sifat data Anda, ada baiknya melengkapi analisis ini berdasarkan metode deret waktu dengan model untuk menghitung data (misalnya model Poisson) dan menguji signifikansi musiman pada model itu. Menambahkan data jumlah tag ke pertanyaan Anda dapat menghasilkan lebih banyak tampilan dan kemungkinan jawaban dalam arah ini.
sumber
Seperti disebutkan dalam komentar saya, ini adalah masalah yang sangat menarik. Mendeteksi musiman bukanlah latihan statistik saja. Pendekatan yang masuk akal adalah berkonsultasi dengan teori dan pakar seperti:
pada masalah ini untuk memahami "mengapa" akan ada musiman untuk melengkapi analisis data. Datang ke data, saya menggunakan metode dekomposisi yang sangat baik yang disebut unobserved components model (UCM) yang merupakan bentuk metode state space. Lihat juga artikel yang sangat mudah diakses oleh Koopman ini. Pendekatan saya mirip dengan @Javlacalle. Itu tidak hanya menyediakan alat untuk menguraikan data deret waktu tetapi juga secara obyektif menilai ada atau tidaknya musiman melalui pengujian signifikansi. Saya bukan penggemar pengujian signifikansi pada data non-eksperimental, tetapi saya tidak tahu prosedur lain yang dapat Anda uji hipotesis Anda tentang ada / tidaknya musiman pada data time series.
Banyak yang mengabaikan tetapi fitur yang sangat penting yang ingin dipahami adalah jenis musiman:
Untuk data seri waktu yang panjang seperti milik Anda, ada kemungkinan bahwa musiman mungkin berubah seiring waktu. Sekali lagi UCM adalah satu-satunya pendekatan yang saya tahu yang dapat mendeteksi musiman stokastik / deterministik ini. UCM dapat menguraikan masalah Anda menjadi "komponen" berikut:
Anda juga dapat menguji apakah level, kemiringan, siklus adalah deterministik atau stokastik. Harap perhatikan itu
level + slope = trend
. Di bawah ini saya menyajikan beberapa analisis pada data Anda menggunakan UCM. Saya menggunakan SAS untuk melakukan analisis.Setelah beberapa iterasi mempertimbangkan komponen dan kombinasi yang berbeda, saya mengakhiri dengan model pelit dari bentuk berikut:
Ada tingkat stokastik + musiman deterministik + beberapa outlier dan data tidak memiliki fitur lain yang terdeteksi.
Di bawah ini adalah analisis signifikansi berbagai komponen. Perhatikan bahwa saya menggunakan trigonometri (yaitu sin / cos dalam pernyataan musiman di PROC UCM) yang mirip dengan @Elvis dan @Nick Cox. Anda juga bisa menggunakan dummy coding di UCM dan ketika saya menguji keduanya memberikan hasil yang sama. Lihat dokumentasi ini untuk perbedaan antara dua cara untuk memodelkan musiman di SAS.
Seperti yang ditunjukkan di atas, Anda memiliki outlier: dua pulsa dan satu tingkat pergeseran pada tahun 2009 (Apakah gelembung ekonomi / perumahan berperan setelah 2009 ??) yang dapat dijelaskan dengan analisis menyelam lebih dalam. Fitur yang baik untuk digunakan
Proc UCM
adalah memberikan output grafis yang luar biasa.Di bawah ini adalah musiman dan tren gabungan dan plot musiman. Apa pun yang tersisa adalah kebisingan .
Tes diagnostik yang lebih penting jika Anda ingin menggunakan nilai p dan pengujian signifikansi adalah untuk memeriksa apakah residu Anda kurang pola dan terdistribusi normal yang dipenuhi dalam model di atas menggunakan UCM dan seperti yang ditunjukkan di bawah ini dalam plot diagnostik residual seperti acf / pacf dan lain-lain.
Kesimpulan : Berdasarkan analisis data menggunakan UCM dan pengujian signifikansi data tampaknya memiliki musiman dan kami melihat jumlah kematian yang tinggi di bulan-bulan musim panas Mei / Juni / Juli dan terendah di bulan-bulan musim dingin di bulan Desember dan Februari.
Pertimbangan Tambahan : Harap pertimbangkan juga arti praktis dari besarnya variasi musiman. Untuk meniadakan argumen kontrafaktual, silakan berkonsultasi dengan pakar domain untuk lebih melengkapi dan memvalidasi hipotesis Anda.
Saya sama sekali tidak mengatakan bahwa ini adalah satu-satunya pendekatan untuk menyelesaikan masalah ini. Fitur yang saya sukai tentang UCM adalah fitur ini memungkinkan Anda untuk memodelkan semua fitur deret waktu secara eksplisit dan juga sangat visual.
sumber
gretl
Untuk estimasi visual awal, grafik berikut dapat digunakan. Merencanakan data bulanan dengan kurva loess dan interval kepercayaan 95%, tampak bahwa ada kenaikan pertengahan tahun memuncak pada Juni. Faktor-faktor lain mungkin menyebabkan data memiliki distribusi yang luas, maka tren musiman mungkin tertutupi dalam plot loess data mentah ini. Poin data telah gugup.
Sunting: Plot berikut menunjukkan kurva loess dan interval kepercayaan untuk perubahan jumlah kasus dari jumlah di bulan sebelumnya:
Ini juga menunjukkan bahwa selama bulan-bulan di paruh pertama tahun ini, jumlah kasus terus meningkat, sementara mereka turun di paruh kedua tahun ini. Ini juga menunjukkan puncaknya di pertengahan tahun. Namun, interval kepercayaan yang luas dan melampaui 0, yaitu tidak ada perubahan, sepanjang tahun, menunjukkan kurangnya signifikansi statistik.
Perbedaan angka satu bulan dapat dibandingkan dengan rata-rata nilai 3 bulan sebelumnya:
Ini menunjukkan peningkatan yang jelas dalam jumlah di bulan Mei dan penurunan di bulan Oktober.
sumber