泊松回归呢？

我创建了一个数据框，其中包含您的数据、t时间索引（以月为单位）和monthdays每个月天数的变量。

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

所以它看起来像这样：

> 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

现在让我们将具有时间效应和天数效应的模型与添加了月份效应的模型进行比较：

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

以下是“小”模型的总结：

> 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

您可以看到这两个变量具有很大的边际效应。现在看看更大的模型：

> 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

好吧，monthdays效果当然消失了；只有闰年才能估计！！将其保留在模型中（并考虑闰年）允许使用剩余偏差来比较两个模型。

> 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

那么，没有（显着的）月份效应？但是季节性影响呢？我们可以尝试使用两个变量和来捕捉季节性：$\cos\left( {2\pi t \over 12}\right)$$\sin\left( {2\pi t \over 12}\right)$

> 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

现在将其与空模型进行比较：

> 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

所以，可以肯定地说，这表明了季节性效应......

卡方检验是一种很好的方法，可以作为您问题的初步观点。

stl作为测试季节性存在的工具，分解可能会产生误导。即使将白噪声（没有结构的随机信号）作为输入传递，此过程也设法返回稳定的季节性模式。尝试例如：

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

查看由自动 ARIMA 模型选择程序选择的订单也有点冒险，因为季节性 ARIMA 模型并不总是涉及季节性（有关详细信息，请参阅此讨论）。在这种情况下，所选模型会产生季节性周期，@RichardHardy 的评论是合理的，但是，需要进一步了解才能得出自杀是由季节性模式驱动的结论。

下面，我根据对您发布的月度系列的分析总结了一些结果。这是根据基本结构时间序列模型估计的季节性分量：

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)

使用带有默认选项的软件 TRAMO-SEATS 提取了一个类似的组件。我们可以看到，估计的季节性模式随着时间的推移并不稳定，因此不支持样本期间跨月自杀人数反复出现模式的假设。使用默认选项运行软件 X-13ARIMA-SEATS，季节性组合测试得出的结论是不存在可识别的季节性。

编辑（ 有关可识别季节性的定义，请参见下面的此答案和我的评论）。

鉴于您的数据的性质，值得使用计数数据模型（例如泊松模型）来补充基于时间序列方法的分析，并测试该模型中季节性的重要性。将标签计数数据添加到您的问题可能会在这个方向上产生更多视图和潜在答案。

正如我在评论中指出的，这是一个非常有趣的问题。检测季节性不仅仅是一项统计工作。一个合理的方法是咨询理论和专家，例如：

• 心理学家
• 精神科医生
• 社会学家

对这个问题要了解“为什么”会有季节性来补充数据分析。谈到数据，我使用了一种出色的分解方法，称为未观察组件模型(UCM)，它是状态空间方法的一种形式。另请参阅 Koopman 的这篇非常易于理解的文章。我的方法类似于@Javlacalle。它不仅提供了分解时间序列数据的工具，而且还通过显着性检验客观地评估了季节性的存在与否。我不是对非实验数据进行显着性测试的忠实粉丝，但我不知道有任何其他程序可以测试您对时间序列数据是否存在季节性的假设。

许多人忽略了一个非常重要的特征，人们想要了解的是季节性的类型：

1. 随机 - 随机变化且难以预测
2. 确定性 - 不会改变，完全可以预测。您可以使用虚拟或三角函数（正弦/余弦等）来建模

对于像您这样的冗长时间序列数据，季节性可能会随着时间而改变。同样，UCM 是我所知道的唯一可以检测这些随机/确定性季节性的方法。UCM 可以将您的问题分解为以下“组件”：

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

您还可以测试水平、斜率、周期是确定性的还是随机的。请注意level + slope = trend. 下面我将使用 UCM 对您的数据进行一些分析。我使用SAS进行分析。

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;

在考虑了不同的组件和组合的几次迭代之后，我以以下形式的简约模型结束：

存在随机水平 + 确定性季节性 + 一些异常值，并且数据没有任何其他可检测的特征。

以下是各种成分的显着性分析。请注意，我使用了类似于@Elvis 和@Nick Cox 的三角函数（即 PROC UCM 中季节性语句中的 sin/cos）。您也可以在 UCM 中使用虚拟编码，当我测试时，两者都给出了相似的结果。有关在 SAS 中建模季节性的两种方法之间的差异，请参阅此文档。

如上所示，您有异常值：2009 年的两个脉冲和一个水平转变（经济/房地产泡沫在 2009 年之后发挥了作用吗？？）这可以通过进一步的深入分析来解释。使用的一个很好的特点Proc UCM是它提供了出色的图形输出。

以下是季节性以及趋势和季节性的组合图。剩下的就是噪音。

如果要使用 p 值和显着性测试，更重要的诊断测试是检查残差是否无模式且正态分布，这在使用 UCM 的上述模型中得到满足，如下面的残差诊断图中所示，如 acf/pacf和别的。

结论：基于使用 UCM 的数据分析和显着性检验，数据似乎具有季节性，我们看到 5 月/6 月/7 月的夏季死亡人数较多，而冬季 12 月和 2 月的死亡人数最低。

其他注意事项：还请考虑季节性变化幅度的实际意义。要否定反事实论点，请咨询领域专家以进一步补充和验证您的假设。

我绝不是说这是解决这个问题的唯一方法。我喜欢 UCM 的特点是它允许您显式地对所有时间序列特征进行建模，并且还具有高度的可视性。

对于初始视觉估计，可以使用下图。用黄土曲线及其 95% 的置信区间绘制月度数据，似乎有一个在 6 月达到峰值的年中上涨。其他因素可能导致数据分布广泛，因此季节性趋势可能在此原始数据黄土图中被掩盖。数据点已经抖动。

编辑：下图显示黄土曲线和上个月病例数变化的置信区间：

这也表明，在上半年的几个月里，病例数一直在上升，而在下半年则在下降。这也暗示了年中的一个高峰。然而，置信区间很宽并且跨越0，即全年没有变化，表明缺乏统计学意义。

一个月数字的差异可以与前 3 个月的平均值进行比较：

这表明 5 月份数量明显增加，10 月份数量下降。