Est-ce une méthode appropriée pour tester les effets saisonniers des données sur le nombre de suicides?

24

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).

entrez la description de l'image ici

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. entrez la description de l'image ici

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 stlfonction 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")

entrez la description de l'image ici

     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)

entrez la description de l'image ici

À 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

entrez la description de l'image ici

> 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))

entrez la description de l'image ici

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.arimachoisi 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

daily_suicide_data.csv

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.

entrez la description de l'image ici

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.

state.daylight.2002.csv

[ 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)]

entrez la description de l'image ici

svannoy
la source
1
L'expression "l'un de nos 50 États" implique que tous les lecteurs appartiennent aux États-Unis. Manifestement, de nombreux extraterrestres se cachent ici aussi.
Nick Cox
1
Est-ce à partir d'un ensemble de données public? Pourriez-vous rendre les données hebdomadaires ou même quotidiennes disponibles?
Elvis
1
@Elvis - J'ai publié un lien vers les données de comptage quotidiennes. Les données proviennent de certificats de décès qui sont des «archives publiques» mais nécessitent un processus pour les obtenir; cependant, les données de comptage agrégées ne le font pas. PS - J'ai essayé le lien moi-même et cela a fonctionné, mais je n'ai pas posté dans un dossier public de cette manière auparavant, alors faites-moi savoir si le lien ne fonctionne pas.
svannoy
1
Étant donné que vos données sont des chiffres, je m'attends à ce que la variance soit liée à la moyenne. Les modèles de séries chronologiques habituels ne tiennent pas compte de cela (cependant, vous pouvez essayer de dire une transformation , peut-être un Freeman-Tukey , par exemple), ou vous pouvez regarder un modèle de série chronologique conçu pour les données de comptage. (Si vous ne le faites pas, ce ne sera peut-être pas un gros problème car le nombre ne dépasse que deux ou deux.)
Glen_b -Reinstate Monica
1
prévisionniste - parce que la répartition des données en nombre est liée à la moyenne. par exemple, considérons un nombre de Poisson , avec une moyenne . Dans ce cas, , donc à mesure que la moyenne augmente, la variance (et donc aussi l'écart) augmente. [En effet, vous trouvez souvent que la propagation est liée à la moyenne d'une manière ou d'une autre pour presque toutes les données qui ont une limite supérieure ou inférieure; il y a des raisons claires pour lesquelles cela pourrait être prévu dans ces cas.]μ t Var ( y t ) = μ tytμtVar(yt)=μt
Glen_b -Reinstate Monica

Réponses:

13

Et la régression de Poisson?

J'ai créé un bloc de données contenant vos données, plus un index tpour le temps (en mois) et une variable monthdayspour le nombre de jours de chaque mois.

T <- read.table("suicide.txt", header=TRUE)
U <- data.frame( year = as.numeric(rep(rownames(T),each=12)), 
         month = rep(colnames(T),nrow(T)), 
         t = seq(0, length = nrow(T)*ncol(T)), 
         suicides = as.vector(t(T)))
U$monthdays <- c(31,28,31,30,31,30,31,31,30,31,30,31)
U$monthdays[ !(U$year %% 4) & U$month == "Feb" ] <- 29

Cela ressemble donc à ceci:

> head(U,14)
   year month  t suicides monthdays
1  1995   Jan  0       62        31
2  1995   Feb  1       47        28
3  1995   Mar  2       55        31
4  1995   Apr  3       74        30
5  1995   May  4       71        31
6  1995   Jun  5       70        30
7  1995   Jul  6       67        31
8  1995   Aug  7       69        31
9  1995   Sep  8       61        30
10 1995   Oct  9       76        31
11 1995   Nov 10       68        30
12 1995   Dec 11       68        31
13 1996   Jan 12       64        31
14 1996   Feb 13       69        29

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:

> a0 <- glm( suicides ~ t + monthdays, family="poisson", data = U )
> a1 <- glm( suicides ~ t + monthdays + month, family="poisson", data = U )

Voici le résumé du "petit" modèle:

> summary(a0)

Call:
glm(formula = suicides ~ t + monthdays, family = "poisson", data = U)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-2.7163  -0.6865  -0.1161   0.6363   3.2104

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8060135  0.3259116   8.610  < 2e-16 ***
t           0.0013650  0.0001443   9.461  < 2e-16 ***
monthdays   0.0418509  0.0106874   3.916 9.01e-05 ***
---
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: 196.64  on 201  degrees of freedom
AIC: 1437.2

Number of Fisher Scoring iterations: 4

Vous pouvez voir que les deux variables ont des effets marginaux largement significatifs. Regardez maintenant le modèle plus grand:

> summary(a1)

Call:
glm(formula = suicides ~ t + monthdays + month, family = "poisson",
    data = U)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-2.56164  -0.72363  -0.05581   0.58897   3.09423

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.4559200  2.1586699   0.674    0.500
t            0.0013810  0.0001446   9.550   <2e-16 ***
monthdays    0.0869293  0.0719304   1.209    0.227
monthAug    -0.0845759  0.0832327  -1.016    0.310
monthDec    -0.1094669  0.0833577  -1.313    0.189
monthFeb     0.0657800  0.1331944   0.494    0.621
monthJan    -0.0372652  0.0830087  -0.449    0.653
monthJul    -0.0125179  0.0828694  -0.151    0.880
monthJun     0.0452746  0.0414287   1.093    0.274
monthMar    -0.0638177  0.0831378  -0.768    0.443
monthMay    -0.0146418  0.0828840  -0.177    0.860
monthNov    -0.0381897  0.0422365  -0.904    0.366
monthOct    -0.0463416  0.0830329  -0.558    0.577
monthSep     0.0070567  0.0417829   0.169    0.866
---
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: 182.72  on 190  degrees of freedom
AIC: 1445.3

Number of Fisher Scoring iterations: 4

Eh bien, bien sûr, l' monthdayseffet 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.

> anova(a0, a1, test="Chisq")
Analysis of Deviance Table

Model 1: suicides ~ t + monthdays
Model 2: suicides ~ t + monthdays + month
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       201     196.65
2       190     182.72 11   13.928    0.237

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)sin(2πt12)

> a2 <- glm( suicides ~ t + monthdays + cos(2*pi*t/12) + sin(2*pi*t/12),
             family="poisson", data = U )
> summary(a2)

Call:
glm(formula = suicides ~ t + monthdays + cos(2 * pi * t/12) +
    sin(2 * pi * t/12), family = "poisson", data = U)

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.8676170  0.3381954   8.479  < 2e-16 ***
t                   0.0013711  0.0001444   9.493  < 2e-16 ***
monthdays           0.0397990  0.0110877   3.589 0.000331 ***
cos(2 * pi * t/12) -0.0245589  0.0122658  -2.002 0.045261 *
sin(2 * pi * t/12)  0.0172967  0.0121591   1.423 0.154874
---
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

Maintenant, comparez-le avec le modèle nul:

> anova(a0, a2, test="Chisq")
Analysis of Deviance Table

Model 1: suicides ~ t + monthdays
Model 2: suicides ~ t + monthdays + cos(2 * pi * t/12) + sin(2 * pi *
    t/12)
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       201     196.65                   
2       199     190.38  2   6.2698   0.0435 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Donc, on peut sûrement dire que cela suggère un effet saisonnier ...

Elvis
la source
2
J'aime beaucoup l'idée de régression de Poisson (+1), je pense que c'est une hypothèse de modélisation beaucoup plus éduquée mais " sûrement dire " sur la base d'une valeur de 0,044? p
J'amorcerais
2
Je suis totalement d'accord, c'est ce que j'impliquais :) "On peut sûrement dire que cela suggère un effet"; et non "Il y a un effet"! Ce que je pense intéressant, c'est cette transformation trigonométrique, elle est très naturelle et je ne comprends pas pourquoi elle n'est pas vue plus souvent. Mais ce n'est qu'un point de départ ... l'amorçage, l'évaluation du modèle ... beaucoup à faire.
Elvis
1
Pas de problème! Mon mal alors, je n'ai pas pu détecter la blague. :)
usεr11852 dit Réintégrer Monic le
2
+1. Poisson rencontre Fourier ... Je pense que les économistes et certains autres mettent l'accent sur les variables indicatrices parce que la saisonnalité peut être pointue, mais l'approche trigonométrique aide souvent.
Nick Cox
2
Effectivement. Une revue de tutoriel que j'ai écrite est accessible sur stata-journal.com/article.html?article=st0116
Nick Cox le
8

Un test du chi carré est une bonne approche comme vue préliminaire à votre question.

La stldé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:

plot(stl(ts(rnorm(144), frequency=12), s.window="periodic"))

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:

require(stsm)
m <- stsm.model(model = "llm+seas", y = x)
fit <- stsmFit(m, stsm.method = "maxlik.td.scoring")
plot(tsSmooth(fit)$states[,2], ylab = "")
mtext(text = "seasonal component", side = 3, adj = 0, font = 2)

composante saisonnière estimée

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.

javlacalle
la source
Merci @javiacalle, je vais enquêter sur les méthodes que vous proposez. Puis-je vous interroger sur votre conclusion du graphique que vous avez publié, est-ce le fait que l'amplitude augmente au fil des années qui est à la base de votre commentaire, "nous pouvons voir que le modèle saisonnier estimé n'est pas stable dans le temps", ou est-ce que cela inclure l'observation plus subtile que la forme de chaque pic est légèrement différente? Je suppose que le premier, mais nous savons où les hypothèses nous mènent.
svannoy
2
χ2
@svannoy La principale conclusion basée sur les méthodes de séries chronologiques utilisées dans ma réponse est qu'il n'y a pas de tendance saisonnière claire dans les données de l'échantillon. Bien que les cycles saisonniers expliquent une partie de la variabilité des données, un modèle saisonnier ne peut pas être identifié de manière fiable car il est masqué par un degré élevé de fluctuations irrégulières (cela pourrait également être vérifié en affichant la fonction de gain du modèle ARIMA choisi montré dans la question) .
javlacalle
@oDDsKooL J'ai également fait le test du chi carré le jour de la semaine, samedi / dimanche sont un peu en deçà des attentes et lundi / mardi sont juste au-dessus ....
svannoy
6

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:

  • Psychologue
  • Psychiatre
  • Sociologue

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é:

  1. Stochastique - changements aléatoires et difficiles à prévoir
  2. Déterministe - ne change pas, parfaitement prévisible. Vous pouvez utiliser un mannequin ou une trigonométrie (sin / cos etc.,) pour modéliser

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:

Time Series Data = level + Slope + Seasonality + Cycle + Causal + Error(Noise).

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.

data input;
format date mmddyy10.;
date = intnx( 'month', '1jan1995'd, _n_-1 );
      input deaths@@;
datalines;
62    47  55  74  71  70  67  69  61  76  68  68
64    69  68  53  72  73  62  63  64  72  55  61
71    61  64  63  60  64  67  50  48  49  59  72
67    54  72  69  78  45  59  53  48  65  64  44
69    64  65  58  73  83  70  73  58  75  71  58
60    54  67  59  54  69  62  60  58  61  68  56
67    60  54  57  51  61  67  63  55  70  54  55
65    68  65  72  79  72  64  70  59  66  63  66
69    50  59  67  73  77  64  66  71  68  59  69
68    61  66  62  69  84  73  62  71  64  59  70
67    53  76  65  77  68  65  60  68  71  60  79
65    54  65  68  69  68  81  64  69  71  67  67
77    63  61  78  73  69  92  68  72  61  65  77
67    73  81  73  66  63  96  71  75  74  81  63
80    68  76  65  82  69  74  88  80  86  78  76
80    77  82  80  77  70  81  89  91  82  71  73
93    64  87  75  101 89  87  78  106 84  64  71
;
run;

ods graphics on;
 proc ucm data = input plots=all; 
      id date interval = month; 
      model deaths ; 
      irregular ; 
      level checkbreak; 
      season length = 12 type=trig var = 0 noest; * Note I have used trigonometry to model seasonality;
   run;

   ods graphics off;

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.

entrez la description de l'image ici

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.

entrez la description de l'image ici

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 UCMest 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 .

entrez la description de l'image ici entrez la description de l'image ici

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.

entrez la description de l'image ici

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.

prévisionniste
la source
Merci pour cette réponse et pour les commentaires intéressants. Je ne connais pas UCM, cela semble très intéressant, je vais essayer de garder cela à l'esprit ...
Elvis
(+1) Analyse intéressante. Je serais toujours prudent quant à la conclusion de la présence d'un schéma saisonnier déterministe significatif, mais vos résultats appellent à un examen plus approfondi de cette possibilité. Le test de Canova et Hansen pour la stabilité saisonnière pourrait être utile, il est décrit par exemple ici .
javlacalle
gretlπ/6
1
+1. Beaucoup de commentaires intéressants et utiles. À votre liste de psychologue, psychiatre, sociologue pourrait s'ajouter météorologue / climatologue. Une telle personne voudrait ajouter qu'il n'y a pas deux années identiques dans les cycles de précipitations et de température. J'aurais deviné grossièrement à plus de dépression en hiver (des journées plus courtes, etc.), mais tant pis pour une estimation étant donné certaines données.
Nick Cox
Merci @forecaster, cela ajoute beaucoup à mon apprentissage. Je suis psychologue, diplômée en santé publique. J'ajouterais l'épidémiologiste à votre liste. Comme je l'ai mentionné au début, il y a beaucoup de mythologie (aka théorisation) sur les tendances saisonnières et le suicide. On peut argumenter fortement en faveur des tendances saisonnières, dans n'importe quelle direction, nous avons donc besoin d'analyses quantitatives pour (dé) confirmer. Du point de vue de la santé publique, si nous découvrions de fortes discontinuités, nous pourrions cibler les interventions. Je ne vois pas cela dans ces données. Du point de vue de la théorie du suicide, la confirmation de tendances même infimes pourrait éclairer le développement de la théorie.
svannoy
1

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.

entrez la description de l'image ici

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:

entrez la description de l'image ici

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:

entrez la description de l'image ici

Cela montre une nette augmentation des effectifs en mai et une baisse en octobre.

rnso
la source
(-1) Il existe déjà trois réponses de qualité à cette question. Votre réponse ne répond pas non plus à la question publiée - vous pouvez la poster en tant que commentaire . Vous ne fournissez pas de réponse sur la façon dont ces données pourraient être analysées.
Tim
J'avais précédemment posté un commentaire ici (voir ci-dessous la question), mais je ne peux pas poster de chiffre dans les commentaires.
rnso
Bien que je comprenne l'éditorial ici, je dirai que @rnso a fourni un joli graphique qui illustre bien la composante saisonnière potentielle et aurait dû faire partie de mon message d'origine.
svannoy
Je comprends cela et je suis d'accord, mais ce n'est toujours pas une réponse mais plutôt un commentaire ou une amélioration. @rnso aurait pu vous suggérer via un commentaire que vous pourriez regarder ou inclure un tel complot.
Tim
@Glen_b, @ Tim: J'ai ajouté une autre intrigue qui pourrait être utile et que je ne peux pas mettre dans un commentaire.
rnso