我可以用 GAMM 估计 R 中随时间变化的季节性影响吗?

can I estimate a time varying seasonal effect in R with GAMM?

我想使用广义加性模型来研究 R 中的 time-series 数据。我的数据是每月的,我想估计季节性效应和更长的 运行 趋势效应。我关注了 Gavin Simpson 的一些有用帖子 here and here:

我的数据是这样的:

我的 github page:

上有完整的数据集

我试图指定一个具有平滑季节性和趋势项的广义加性模型,如下所示:

    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 方法分析了 these 数据:

使用 STL 范例,季节性影响的波动或平滑程度似乎是一个偏好或选择的问题。如果我能让数据告诉我随机误差和变化的季节性峰值之间的区别,我会更愿意。 GAMS 似乎更适合这个目标,因为它们更容易进行统计模型 fitting-type 练习……但我想知道 R 包中是否有一个参数用于拟合 gams,允许随时间变化的季节性影响。

答案是:是的,可以为您感兴趣的问题制定一个 GAM 模型。如果我们假设模型的趋势和季节性成分平滑地相互作用,我们就有一个平滑等价于连续-连续相互作用。可以使用两个边缘平滑的张量积在 GAM 中拟合这种相互作用:

  1. 季节性循环平滑,
  2. 长期趋势平稳

顺便说一句,我还有关于这些的更多博文:

阅读这些内容以获取更多详细信息,但基本方面要符合以下模型:

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

第二个模型的替代方法是对两个主效应和交互作用进行类似方差分析的分解。在上面的修改模型中,所有三个分量都绑定在单个张量积平滑中,即模型的 te() 部分。

类方差分析分解变体将使用

进行拟合
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() 一起使用。