r - 是否可以在 mgcv gam 模型中包含两个平滑项的乘积

Question

我使用 gam 对时间序列数据的季节性建模取得了巨大成功。除了季节性变化之外，我的最新模型清楚地显示了每周模式。虽然每周模式本身在一年中非常稳定，但其幅度也随季节变化。所以理想情况下，我想将我的数据建模为:

y ~ f(day in year) + g(day in year) * h(day in week)

哪里 f , g和 h是 mgcv 中的循环平滑函数

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc'):s(day_in_week, k=5, bs='cc')
  , knots=list(
    day_in_year=c(0, 356)
    , day_in_week=c(0,7)
  )
  , data = data
)

不幸的是，这不起作用并抛出错误 NA/NaN argument .我尝试使用 te(day_in_year, day_in_week, k=c(52, 5), bs='cc')这是有效的，但引入了太多的自由度，因为该模型过度拟合了在可用年份较短的特定工作日内的假期。

是否可以按照我尝试的方式指定模型？

Answer 1

哇，一个很老的问题!

关于互动

While the weekly pattern itself is very stable over the year, its amplitude also varies with the seasons.

张量积样条基的使用 te是正确的交互方式，虽然更合适的构造函数是 ti .你说的是 te返回许多参数。当然。您有 k = 52第一个边距，和 k = 5第二个，那么你最终会得到 52 * 5 - 1这个张量项的系数。但这只是创建交互的方式。

请注意，在 mgcv GAM 公式， :或 *仅对参数项之间的交互有效。平滑之间的交互由 te 处理或 ti .

如果这不是您所希望的，那么您期望“产品”是什么？两个边际设计矩阵的 Hadamard 乘积？那么这有什么意义呢？顺便说一下，Hadamard 乘积需要两个相同维度的矩阵。但是，您的两个边距的列数不同。

如果你不明白为什么我一直在谈论矩阵，那么你需要阅读 Simon 在 2006 年的书。虽然 GAM 估计解释现在已经过时了，GAM 的构造/设置(如设计矩阵和惩罚矩阵)在章节中解释4 即使在十年之后也根本没有改变。

好的，让我再给你一个提示。 行明智 Kronecker 产品 用于建筑 te/ ti设计矩阵并不是一项新发明。

平稳期 s(x)很像一个因子变量 g ，就好像它们看起来是一个单一的变量，它们被构造为具有许多列的设计矩阵。对于 g它是一个虚拟矩阵，而对于 f(x)它是一个基矩阵。 因此，两个平滑函数之间的交互作用的构建方式与构建两个因素之间的交互作用的方式相同。

如果你有一个因素 g1 5 个级别，还有另一个因素 g2 10个级别，他们的边际设计矩阵(对比后)有4列和9列，那么交互作用 g1:g2将有 36 列。这样的设计矩阵，就是 g1的设计矩阵的行式Kronecker积。和 g2 .

关于过拟合

正如你所说，你只有几年的数据，也许是 2 或 3 年？在这种情况下，您的模型已通过使用 k = 52 过度参数化。为 day_in_year .尝试将其简化为例如 k = 30 .

如果过度拟合仍然很明显，这里有一些方法可以解决它。

方式一:您正在使用 GCV 进行平滑度选择。试试 method = "REML" . GCV 总是倾向于过度拟合数据。

方式二:继续使用 GCV，手动增加平滑参数以获得更重的惩罚。 gamma gam 的参数在这里很有用。试试 gamma = 1.5例如。

你的代码有错别字？

结点位置，应该是 day_in_year = c(0, 365) ?