我能用GAMM在R中估计时间变化的季节效应吗？

Question

我希望在R中使用广义可加模型来研究时间序列数据。我的数据是按月计算的，我想估计季节效应和长期趋势效应。我已经参考了Gavin Simpson的一些有用文章这里和这里。

我的数据长这样：

我在我的 GitHub 页面 上有完整的数据集:

我尝试使用平滑的季节性和趋势项来指定广义可加模型，如下所示:

    df <- read.csv('trips.csv')
    head(df)
    # A tibble: 276 × 2
     date ntrips
   <date>  <int>
    1  1994-01-01    157
    2  1994-02-01    169
    3  1994-03-01    195
    4  1994-04-01    124
    5  1994-05-01    169

    #add a time column
    trips <- tbl_df(trips) %>% mutate(time=as.numeric(date))

    mod1 <- gamm(ntrips~s(month,bs="cc",k=12) + s(time),data=trips)

我按照以下方式提取了季节效应的估计值：

    pred <- predict(mod1$gam,newdata=trips,type="terms")
    seas <- data.frame(s=pred[,1],date=trips$date)
    ggplot(seas,aes(x=date,y=s)) + geom_line()

以下包含此图：

我的问题是：在原始数据中，季节性高峰每年会有些移动。在我指定的极其简单的GAM中，季节效应是恒定的。有没有一种方法可以适应时间变化的季节性与GAM？

我使用Cleveland等人的STL方法分析了这些数据：

使用STL范例，允许季节性效应的波动或平滑似乎是个人偏好或选择的问题。我希望能够让数据告诉我随机误差和移动季节峰值之间的区别。 GAMS似乎更适合这个目标，因为它们更容易用于统计模型拟合类型的练习……但我想知道R软件包中是否有适用于拟合gams的参数，允许时间变化的季节性影响。

- aaronmams

Answer 1

答案是：可以针对您感兴趣的问题制定GAM模型。如果我们假设模型的趋势和季节成分相互平滑地交互作用，那么我们就有了一个连续-连续交互的平滑等效物。这种交互可以使用两个边际平滑的张量积来适应GAM：

1. 季节周期平滑；
2. 长期趋势平滑。

顺便说一下，我还有关于这些的进一步博客文章：

• https://www.fromthebottomoftheheap.net/2015/11/21/climate-change-and-spline-interactions/ 和
• https://www.fromthebottomoftheheap.net/2015/11/23/are-some-seasons-warming-more-than-others/

阅读这些内容以获取更多详细信息，但基本要素是拟合以下模型：

## fix the knots are just passed the ends of the month numbers
## this allows Dec and Jan to have their own estimates
knots <- list(month = c(0.5, 12.5))

## original model, fixed seasonal component
m1 <- gam(ntrips ~ s(month, bs="cc", k=12) + s(time), data = trips,
          knots = knots)

## modified model with
m2a <- gam(ntrips ~ te(month, time, bs = c("cc","tp"), k = c(12, 10)), data = trips,
          knots = knots))

第二个模型的替代方案是对两个主要影响及其交互作用进行类似于ANOVA的分解。在上面修改后的模型中，所有三个组成部分都包含在单个张量积平滑函数中，即模型的te()部分。

类似于ANOVA的分解变体将使用以下方法进行拟合：

m2b <- gam(ntrips ~ ti(month, bs = 'cc', k = 12) +
             ti(time, bs = 'tp', k = 10) +
             ti(month, time, bs = c("cc","tp")), data = trips,
           knots = knots)

第三个ti()是平滑交互，与季节性和长期趋势的主要平滑效应分开。

我用gam()展示了这些拟合结果，但如果您需要在模型残差中包括ARMA过程，也可以使用gamm()。

- Gavin Simpson