J'ai 17 ans (1995 à 2011) de données sur les certificats de décès liés aux décès par suicide pour un État aux États-Unis.Il y a beaucoup de mythologie au sujet des suicides et des mois / saisons, en grande partie contradictoires, et de la littérature I ' ve revu, je n'ai pas une idée claire des méthodes utilisées ou la confiance dans les résultats.
J'ai donc cherché à savoir si je pouvais déterminer si les suicides étaient plus ou moins susceptibles de se produire au cours d'un mois donné dans mon ensemble de données. Toutes mes analyses sont faites dans R.
Le nombre total de suicides dans les données est de 13 909.
Si vous regardez l'année avec le moins de suicides, ils se produisent sur 309/365 jours (85%). Si vous regardez l'année avec le plus de suicides, ils surviennent sur 339/365 jours (93%).
Il y a donc pas mal de jours chaque année sans suicide. Cependant, lorsqu'ils sont regroupés sur les 17 années, il y a des suicides tous les jours de l'année, y compris le 29 février (bien que seulement 5 lorsque la moyenne est de 38).
La simple addition du nombre de suicides chaque jour de l'année n'indique pas une saisonnalité claire (à mes yeux).
Au niveau mensuel, les suicides moyens par mois varient de:
(m = 65, sd = 7,4, à m = 72, sd = 11,1)
Ma première approche a été d'agréger l'ensemble de données par mois pour toutes les années et de faire un test du chi carré après avoir calculé les probabilités attendues pour l'hypothèse nulle, qu'il n'y avait pas de variance systématique du nombre de suicides par mois. J'ai calculé les probabilités pour chaque mois en tenant compte du nombre de jours (et en ajustant février pour les années bissextiles).
Les résultats du chi carré n'ont indiqué aucune variation significative par mois:
# So does the sample match expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131
L'image ci-dessous indique le nombre total de mois par mois. Les lignes rouges horizontales sont positionnées aux valeurs attendues pour février, 30 jours et 31 jours respectivement. Conformément au test du chi carré, aucun mois n'est en dehors de l'intervalle de confiance à 95% pour les dénombrements attendus.
Je pensais avoir fini jusqu'à ce que je commence à enquêter sur les données de séries chronologiques. Comme j'imagine que beaucoup de gens le font, j'ai commencé avec la méthode de décomposition saisonnière non paramétrique en utilisant la stl
fonction dans le package stats.
Pour créer les données de séries chronologiques, j'ai commencé avec les données mensuelles agrégées:
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
Et puis effectué la stl()
décomposition
# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)
À ce stade, je me suis inquiété car il me semble qu'il y a à la fois une composante saisonnière et une tendance. Après de nombreuses recherches sur Internet, j'ai décidé de suivre les instructions de Rob Hyndman et George Athanasopoulos comme indiqué dans leur texte en ligne "Prévision: principes et pratiques", spécifiquement pour appliquer un modèle saisonnier ARIMA.
J'ai utilisé adf.test()
et kpss.test()
pour évaluer la stationnarité et j'ai obtenu des résultats contradictoires. Ils ont tous deux rejeté l'hypothèse nulle (notant qu'ils testent l'hypothèse inverse).
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
J'ai ensuite utilisé l'algorithme dans le livre pour voir si je pouvais déterminer la quantité de différenciation qui devait être faite pour la tendance et la saison. J'ai fini avec nd = 1, ns = 0.
J'ai ensuite couru auto.arima
, qui a choisi un modèle qui avait à la fois une tendance et une composante saisonnière avec une constante de type "dérive".
# 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
Enfin, j'ai regardé les résidus de l'ajustement et si je comprends bien, puisque toutes les valeurs sont dans les limites de seuil, elles se comportent comme du bruit blanc et donc le modèle est assez raisonnable. J'ai exécuté un test de portemanteau comme décrit dans le texte, qui avait une valeur ap bien supérieure à 0,05, mais je ne suis pas sûr d'avoir les paramètres corrects.
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
Après être revenu en arrière et relire le chapitre sur la modélisation arima, je me rends compte maintenant que j'ai auto.arima
choisi de modéliser la tendance et la saison. Et je réalise également que les prévisions ne sont pas spécifiquement l'analyse que je devrais probablement faire. Je veux savoir si un mois spécifique (ou plus généralement une période de l'année) doit être signalé comme un mois à haut risque. Il semble que les outils de la littérature de prévision soient très pertinents, mais peut-être pas les meilleurs pour ma question. Toute contribution est très appréciée.
Je publie un lien vers un fichier csv qui contient les comptes quotidiens. Le fichier ressemble à ceci:
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
Le nombre est le nombre de suicides qui se sont produits ce jour-là. "t" est une séquence numérique de 1 au nombre total de jours dans le tableau (5533).
J'ai pris note des commentaires ci-dessous et pensé à deux choses liées à la modélisation du suicide et des saisons. Premièrement, en ce qui concerne ma question, les mois ne sont que des procurations pour marquer le changement de saison, je ne suis pas intéressé à savoir si un mois en particulier est différent des autres (c'est bien sûr une question intéressante, mais ce n'est pas ce à quoi je me suis fixé enquêter). Par conséquent, je pense qu'il est logique d' égaliser les mois en utilisant simplement les 28 premiers jours de tous les mois. Lorsque vous faites cela, vous obtenez un ajustement légèrement pire, que j'interprète comme une preuve supplémentaire d'un manque de saisonnalité. Dans la sortie ci-dessous, le premier ajustement est une reproduction d'une réponse ci-dessous en utilisant les mois avec leur nombre réel de jours, suivi d'un ensemble de données suicideByShortMonthoù le nombre de suicides a été calculé à partir des 28 premiers jours de tous les mois. Je suis intéressé par ce que les gens pensent de savoir si cet ajustement est une bonne idée, pas nécessaire ou nuisible?
> 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
La deuxième chose que j'ai examinée plus est la question de l'utilisation du mois comme proxy pour la saison. Un meilleur indicateur de la saison est peut-être le nombre d'heures de clarté qu'une zone reçoit. Ces données proviennent d'un État du nord qui a une variation substantielle de la lumière du jour. Voici un graphique de la lumière du jour de l'année 2002.
Lorsque j'utilise ces données plutôt que le mois de l'année, l'effet est toujours significatif, mais l'effet est très, très petit. La déviance résiduelle est beaucoup plus grande que les modèles ci-dessus. Si les heures de clarté sont un meilleur modèle pour les saisons et que l'ajustement n'est pas aussi bon, est-ce davantage une preuve d'un très faible effet saisonnier?
> 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
Je poste les heures du jour au cas où quelqu'un voudrait jouer avec ça. Notez que ce n'est pas une année bissextile, donc si vous voulez mettre les minutes pour les années bissextiles, extrapolez ou récupérez les données.
[ Modifier pour ajouter un tracé à partir de la réponse supprimée (j'espère que cela ne me dérange pas de déplacer le tracé dans la réponse supprimée ici à la question. Svannoy, si vous ne voulez pas que cela soit ajouté après tout, vous pouvez le revenir)]
la source
Réponses:
Et la régression de Poisson?
J'ai créé un bloc de données contenant vos données, plus un index
t
pour le temps (en mois) et une variablemonthdays
pour le nombre de jours de chaque mois.Cela ressemble donc à ceci:
Comparons maintenant un modèle avec un effet de temps et un effet de nombre de jours avec un modèle dans lequel nous ajoutons un effet de mois:
Voici le résumé du "petit" modèle:
Vous pouvez voir que les deux variables ont des effets marginaux largement significatifs. Regardez maintenant le modèle plus grand:
Eh bien, bien sûr, l'
monthdays
effet disparaît; il ne peut être estimé que grâce aux années bissextiles !! Le conserver dans le modèle (et prendre en compte les années bissextiles) permet d'utiliser les écarts résiduels pour comparer les deux modèles.Donc, pas d'effet mois (significatif)? Mais qu'en est-il d'un effet saisonnier? Nous pouvons essayer de capturer la saisonnalité en utilisant deux variables et :sin(2πtcos( 2 πt12) péché( 2 πt12)
Maintenant, comparez-le avec le modèle nul:
Donc, on peut sûrement dire que cela suggère un effet saisonnier ...
la source
Un test du chi carré est une bonne approche comme vue préliminaire à votre question.
La
stl
décomposition peut être trompeuse en tant qu'outil pour tester la présence de saisonnalité. Cette procédure parvient à renvoyer un modèle saisonnier stable même si un bruit blanc (signal aléatoire sans structure) est transmis en entrée. Essayez par exemple:Regarder les commandes choisies par une procédure de sélection automatique du modèle ARIMA est également un peu risqué car un modèle saisonnier ARIMA n'implique pas toujours la saisonnalité (pour plus de détails, voir cette discussion ). Dans ce cas, le modèle choisi génère des cycles saisonniers et le commentaire de @RichardHardy est raisonnable, cependant, un nouvel éclairage est nécessaire pour conclure que les suicides sont motivés par un modèle saisonnier.
Ci-dessous, je résume quelques résultats basés sur l'analyse de la série mensuelle que vous avez publiée. Il s'agit de la composante saisonnière estimée sur le modèle de base des séries chronologiques structurelles:
Un composant similaire a été extrait à l'aide du logiciel TRAMO-SEATS avec des options par défaut. Nous pouvons voir que la tendance saisonnière estimée n'est pas stable dans le temps et, par conséquent, ne soutient pas l'hypothèse d'une tendance récurrente du nombre de suicides au cours des mois au cours de la période d'échantillonnage. En exécutant le logiciel X-13ARIMA-SEATS avec des options par défaut, le test combiné de saisonnalité a conclu que la saisonnalité identifiable n'était pas présente.
Modifier (voir cette réponse et mon commentaire ci-dessous pour une définition de la saisonnalité identifiable ).
Compte tenu de la nature de vos données, il serait utile de compléter cette analyse basée sur des méthodes de séries chronologiques avec un modèle pour les données de comptage (par exemple le modèle de Poisson) et de tester la signification de la saisonnalité dans ce modèle. L'ajout des données de comptage des balises à votre question peut entraîner davantage de vues et de réponses potentielles dans ce sens.
la source
Comme indiqué dans mon commentaire, c'est un problème très intéressant. La détection de la saisonnalité n'est pas un simple exercice statistique. Une approche raisonnable serait de consulter la théorie et des experts tels que:
sur ce problème pour comprendre "pourquoi" il y aurait une saisonnalité pour compléter l'analyse des données. Pour en venir aux données, j'ai utilisé une excellente méthode de décomposition appelée modèle de composants non observés (UCM) qui est une forme de méthode de l'espace d'état. Voir aussi cet article très accessible de Koopman. Mon approche est similaire à @Javlacalle. Il fournit non seulement un outil pour décomposer les données de séries chronologiques, mais évalue également objectivement la présence ou l'absence de saisonnalité via des tests de signification. Je ne suis pas un grand fan des tests de signification sur des données non expérimentales, mais je ne connais aucune autre procédure permettant de tester votre hypothèse sur la présence / absence de saisonnalité sur des données de séries chronologiques.
Beaucoup ignorent mais une caractéristique très importante que l'on voudrait comprendre est le type de saisonnalité:
Pour une longue série de données temporelles comme la vôtre, il est possible que la saisonnalité ait changé au fil du temps. Encore une fois, l'UCM est la seule approche que je connaisse qui puisse détecter ces saisonnalités stochastiques / déterministes. UCM peut décomposer votre problème en "composants" suivants:
Vous pouvez également tester si le niveau, la pente, le cycle sont déterministes ou stochastiques. Veuillez le noter
level + slope = trend
. Ci-dessous, je présente une analyse de vos données en utilisant UCM. J'ai utilisé SAS pour faire l'analyse.Après plusieurs itérations prenant en compte différents composants et combinaisons, j'ai terminé avec un modèle parcimonieux de la forme suivante:
Il existe un niveau stochastique + une saisonnalité déterministe + quelques valeurs aberrantes et les données n'ont pas d'autres caractéristiques détectables.
Vous trouverez ci-dessous une analyse de la signification de divers composants. Notez que j'ai utilisé la trigonométrie (c'est-à-dire sin / cos dans la déclaration de saisonnalité dans PROC UCM) similaire à @Elvis et @Nick Cox. Vous pouvez également utiliser un codage factice dans UCM et lorsque j'ai testé, les deux ont donné des résultats similaires. Consultez cette documentation pour connaître les différences entre les deux façons de modéliser la saisonnalité dans SAS.
Comme indiqué ci-dessus, vous avez des valeurs aberrantes: deux impulsions et un changement de niveau en 2009 (la bulle économique / immobilière a-t-elle joué un rôle après 2009 ??), ce qui pourrait s'expliquer par une analyse approfondie plus approfondie. Une bonne caractéristique de l'utilisation
Proc UCM
est qu'elle fournit une excellente sortie graphique.Vous trouverez ci-dessous la saisonnalité et un graphique combiné de tendance et de saisonnalité. Tout ce qui reste est du bruit .
Un test de diagnostic plus important si vous souhaitez utiliser des valeurs de p et un test de signification consiste à vérifier si vos résidus sont sans schéma et normalement distribués, ce qui est satisfait dans le modèle ci-dessus à l'aide de l'UCM et comme indiqué ci-dessous dans les tracés de diagnostic résiduels comme acf / pacf et d'autres.
Conclusion : Sur la base de l'analyse des données utilisant l'UCM et des tests de signification, les données semblent avoir une saisonnalité et nous observons un nombre élevé de décès pendant les mois d'été de mai / juin / juillet et le plus bas pendant les mois d'hiver de décembre et février.
Considérations supplémentaires : veuillez également tenir compte de l'importance pratique de l'ampleur de la variation saisonnière. Pour réfuter les arguments contrefactuels, veuillez consulter des experts du domaine pour compléter et valider votre hypothèse.
Je ne dis nullement que c'est la seule approche pour résoudre ce problème. La fonctionnalité que j'aime dans UCM est qu'elle vous permet de modéliser explicitement toutes les fonctionnalités des séries temporelles et est également très visuelle.
la source
gretl
Pour une estimation visuelle initiale, le graphique suivant peut être utilisé. En traçant les données mensuelles avec la courbe de loess et son intervalle de confiance à 95%, il semble qu'il y ait une augmentation en milieu d'année avec un pic en juin. D'autres facteurs peuvent être à l'origine d'une large distribution des données, d'où la tendance saisonnière peut être masquée dans ce graphique de loess de données brutes. Les points de données ont été instables.
Modifier: le graphique suivant montre la courbe de Loess et l'intervalle de confiance pour le changement du nombre de cas par rapport au nombre du mois précédent:
Cela montre également qu'au cours des mois du premier semestre, le nombre de cas continue d'augmenter alors qu'il diminue au second semestre. Cela suggère également un pic au milieu de l'année. Cependant, les intervalles de confiance sont larges et vont de 0, c'est-à-dire aucun changement, tout au long de l'année, indiquant un manque de signification statistique.
La différence du nombre d'un mois peut être comparée à la moyenne des valeurs des 3 mois précédents:
Cela montre une nette augmentation des effectifs en mai et une baisse en octobre.
la source