R语言动态广义相加模型GAM张量积交互项、傅立叶项、谐波回归分析季节性时间序列航空乘客数据
全文链接：https://tecdat.cn/?p=36497
原文出处：拓端数据部落公众号
季节性在真实的时间序列中是非常常见的。许多系列以周期性、规律性的方式变化。例如，冰淇淋销售在温暖的假期月份往往更高，而候鸟数量围绕年度迁徙周期强烈波动。由于季节性非常普遍，已经开发了许多时间序列和预测方法专门用于处理这个特征。
本文的目的是强调一种策略，用于在动态广义相加模型中捕捉季节性和时变季节性模式。
绘制季节性时间序列
我们将使用AirPassengers数据集，这是一个从1949年到1960年每月国际航空乘客的时间序列（以千为单位）。
从这张图中可以清晰地看出，乘客数量存在强烈的季节性波动，但这种模式也似乎在随时间变化。一种常见的检查可能季节性模式的方法是计算和绘制STL分解（使用Loess的季节性和趋势分解），该方法利用Loess平滑来迭代地平滑和移除数据的某些成分（即可能是非线性的长期趋势和可能随时间变化的季节性）。
autoplot()可以通过修改t.window和s.window参数来调整趋势和季节性成分被平滑的周期，但这些对于本文来说并不是很重要。关键信息是，似乎存在一个长期的非线性趋势和一个随时间变化的循环季节性模式，这种变化不仅体现在幅度上，可能还在形状上略有不同。接下来，我们将介绍一些模型。首先，我们将时间序列转换为适合使用{mgcv}和{mvgam}包进行建模的data.frame()对象。这个函数还会自动将数据拆分为训练集和测试集，以便轻松评估不同模型的预测能力。
在将数据转换为正确的“长”格式之后，我们可以轻松地使用{
mvgam}包的S3函数绘制时间序列的一些特征。
r复制代码
plot_mvgam_s自相关函数（ACF）图清晰地显示了强烈的自相关性和季节性证据，这正是我们希望在我们的模型中捕捉到的。
使用张量积交互项拟合mvgam模型
现在我们可以拟合一个模型来尝试捕捉这两个组成部分，即趋势和季节性模式。
mod1 <- mvgam(y ~
te(season, time,
为了更深入地了解模型中的时间与时节交互平滑函数，我们可以方便地使用gratia::draw()函数从模型底层的gam对象（存储在返回模型对象的mgcv_model槽位中）进行查看。
尽管此模型的样本外预测表现尚可，但不确定性区间可能稍显宽泛。如果你花些时间理解样条如何超出训练数据进行外推，这并不令人意外。在实际应用中，可能需要进一步调整模型或考虑其他方法来缩小预测的不确定性区间。
该模型表现尚可，但出于以下几个原因可能不是最优选择：
- 我们使用样条函数来捕捉时间维度，但众所周知，样条函数在外推时通常表现不佳。
- 该模型假设季节性模式随时间样条函数的变化而变化的速率相同，而这可能并非总是最合适的模型。
使用傅里叶项来捕捉时变季节性
另一种捕捉时变季节性的方法是使用傅里叶变换。傅里叶级数通过不同频率的正弦和余弦项集合，通过加权来近似各种周期性函数。首先，我们使用forecast::fourier()函数来计算正弦和余弦函数，该函数会自动识别AirPassengers数据的月周期性，以计算正确频率的傅里叶项。在这里，我们将阶数K设置为4，这给我们提供了总共8个正弦和余弦预测器，这些预测器可以作为我们动态模型中的回归变量。这种模型通常被称为谐波回归，因为连续的傅里叶项可以被认为是前一个傅里叶项的谐波。
dplyr::glimpse(fou_terms)观察这些傅里叶项（针对前24个时间点）有助于我们理解它们如何被视为一种特定的基函数扩展。
接下来，我们可以将这些新的预测变量添加到训练和测试数据集中。重要的是，我们也需要为样本外测试集添加这些项，因为这将用于帮助我们从模型中计算预测值。
然后，我们可以拟合一个mvgam()模型，该模型包括平滑的非线性时间趋势和时变季节性，但使用傅里叶项代替上面使用的循环样条。这首先通过使用时间的单调样条来实现，确保时间趋势永远不会下降，这在当前情况下似乎是合适的。时变季节性通过允许傅里叶项的系数随时间变化来实现，这通过使用s()包装器中的'by'参数并利用单独的样条来实现。这实际上设置了时间的独立样条，然后与傅里叶项进行交互，为每个傅里叶项制定了一个时变系数模型。
mod2 <- mvgam(y ~)
# 时间单调平滑，以确保趋势
# 要么增加要么变平，但不会
# 下降
s(time, bs = 'moi', k = 10) +
# 时变傅立叶系数，以捕捉
# 可能变化的季节性
s(time, by = s1_12, k = 5) + # 时间变化的傅立叶系数，以捕捉 # 可能变化的季节性。当然，该模型仍然使用时间样条来模拟趋势，但单调性限制确保了预测结果至少在一定程度上是可控的。基于样本内拟合的效果，我们更倾向于选择第二个模型，并可以通过近似留一交叉验证（Approximate Leave-One-Out Cross-Validation）轻松地比较两个模型的性能。在这里，我们更倾向于具有更高期望对数预测密度（ELPD）的模型。
但预测结果如何呢？这里我们绘制了两个模型的样本外预测结果，结果再次表明第二个模型略胜一筹
查看平滑局部效应图，了解某些傅立叶系数的估计值随时间的变化情况
绘制谐波回归的趋势和季节成分
绘制该模型的图需要更多的工作，因为我们需要使用该模型的目标预测来直观地显示两个主要效应（即非线性趋势和随时间变化的季节性）。
首先是趋势图，在这个模型中使用参数可以很容易地计算出模型的期望值（即结果标度上的预测值，但忽略了负二项采样分布的超分散参数）。
现在是时变季节模式，需要从总预测值中减去条件趋势预测值
season_preds <- with_season - agg_r_season现在，我们将这些效果绘制在一起
这幅图清楚地显示了模型是如何估算出这两个组成部分的，以及季节模式在形状和幅度上随时间的变化情况
时变周期性
上例说明了如何利用傅立叶项来捕捉时变季节性。但如果周期性本身随时间变化呢？换句话说，如果一个序列开始时每 8 周出现一次有规律的周期性变化，但随后周期性窗口不断缩小或扩大，我们该如何建立模型呢？下面是一个非常简单的示例，说明我们如何扩展傅立叶策略，以对这种情况进行建模
首先，我定义了构建单调递增或递减系数的函数。这将有助于模拟时变周期性
x <- sort(as.vector(scale(times)))
exp(k * x) / (1 + exp(k * x))现在进行模拟，包括构建两个不同周期的傅里叶级数（本例中为 "k = 12 "和 "k = 6"），每个周期有三个傅里叶项。然后，我模拟这些项的时变系数，迫使周期为 12 的项随时间下降为零，而周期为 6 的项随时间上升。这些系数加上模拟的傅立叶项本身，就可以模拟出一个周期性随时间从 ~12 缩小到 ~6 的时间序列。
绘制随时间变化的傅立叶系数图
现在计算线性预测因子并绘制曲线
最后，添加一些高斯随机漫步噪声，绘制模拟时间序列图
现在我们可以构建用于建模的数据；这需要我们仔细考虑问题，并将两个周期性的傅立叶级数作为预测因子纳入其中
现在我们可以拟合动态 GAM，它使用状态空间表示法来捕捉随时间变化的季节性、周期性和随机漫步自相关过程
和之前一样，绘制时变傅立叶系数图
draw(trend_mgcv_model)和以前一样，我们的模型成功地捕捉到了这些项的时变系数。请注意，许多 12 期项（即 s1_12 和 c3_12）估计会随时间推移而下降，而许多 6 期项（即 s2_6 和 c3_6）估计会随时间推移而上升。令人鼓舞的是，该模型在处理随机漫步自相关过程和两个方差项（过程误差和观测误差）的同时，还能发现这些问题。
我们还可以使用后验期望绘制季节性时变图，就像上面第一个例子那样。但要注意的是，该图的预测不包括随机漫步部分。
现在，我们可以提取估算出的 RW 趋势，并减去季节性预测，这样就可以看到模型估算出的潜在非线性趋势了
treneds <- hinst(mod3, type = 'expected')$hiasts[[1]] -
sepreds将二者绘制在一起。
参考文献
以下论文和资源提供了有关动态谐波回归和使用 GAM 进行时间序列建模/预测的有用资料
Autin, M. A., & Edwards, D. 河口水质数据的非参数谐波回归 Environmetrics 21.6 (2010)： 588-605.
Karunarathna, K. A. N. K., Wells, K., & Clark, N. J. (2024). 用分层动态广义加法模型模拟沙漠啮齿动物物种对环境变化的非线性响应. Ecological Modelling, 490, 110648.
Young, P. C., Pedregal, D. J., & Tych, W. (1999). Dynamic harmonic regression. Journal of Forecasting, 18, 369-394.
