8

我在使用 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')which 有效，但引入了太多的自由度，因为该模型过拟合了在可用年份较短的特定工作日内的假期。

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

4

3 回答 3

6

哇，很老的问题！

关于互动

虽然一周模式本身在一年中非常稳定，但其幅度也随季节而变化。

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

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

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

如果你不明白我为什么一直在谈论矩阵，那么你需要阅读 Simon 在 2006 年的书。虽然 GAM 估计解释现在已经相当过时，但 GAM 的构造/设置（如设计矩阵和惩罚矩阵）在第 1 章中解释4即使十年后也完全没有变化。

好的，让我再给你一个提示。用于构造/设计矩阵的行式 Kronecker 产品并不是一项新发明。teti

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

如果您有g15 个水平的因子和 10 个水平的另一个因子g2，它们的边际设计矩阵（对比后）有 4 列和 9 列，那么交互作用g1:g2将有 36 列。g1这样的设计矩阵，就是和的设计矩阵的逐行克罗内克积g2。

关于过拟合

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

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

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

方式2：继续使用GCV，手动增加平滑参数以获得更重的惩罚。gamma的参数在gam这里很有用。举个gamma = 1.5例子。

您的代码中有错字吗？

结的位置，应该是day_in_year = c(0, 365)？

于 2017-05-17T07:41:25.240 回答
6

您的模型没有多大意义，因为您声明存在：

  1. 季节性影响，
  2. 一周中的一天效果，
  3. 一种交互作用，使得星期几的效果随一年中的某一天而变化

这可以完全表示为张量积平滑。您在此处对其他答案的评论中提到的模型

y ~ f(一年中的一天) + g(一年中的一天) * h(一周中的一天)

如果您的意思*是主效应+交互，则只是对完整张量积的分解。在这种情况下，您拥有的模型是不可识别的-您有两次一年中的一天的功能。如果您的意思是:，那么您的模型没有星期几的主要影响，这似乎是不可取的。

我一直都适合这种形式的模型（仅适用于一年中的一年和一天）。我会通过以下方式解决这个问题：

gam(y ~ te(day_of_year, day_of_week, k = c(20, 6), bs = c("cc", "cc")),
    data = foo, method = "REML", knots = knots)

您还可以调整结定义。我倾向于使用以下内容：

knots <- list(day_of_year = c(0.5, 366.5),
              day_of_week = c(0.5, 7.5)

这不会有太大的区别，但你只是把边界结放在离数据更近一点的地方。

如果要分解效果，可以使用ti()平滑拟合模型

gam(y ~ ti(day_of_year, bs = "cc", k = 12) + ti(day_of_week, bs = "cc", k = 6) + 
      te(day_of_year, day_of_week, k = c(12, 6), bs = c("cc", "cc")),
    data = foo, method = "REML", knots = knots)

您可以调整 的值k以确定与您的数据和gam.check().

您还需要在模型中添加一个术语来处理假期。如果当天是假期，这将是应用调整的参数项 - 因此，创建一个因子holiday并将其添加到模型+ holiday中。您可以考虑更复杂的模型；如果一周有假期，则可能是一个因子索引，再加上该day_of_week组件的逐因子平滑，这样，如果一周是正常的一周，则您估计一个每周模式，如果该周包含假期，则估计第二个每周模式。

如果您向我们展示数据的示例/图，我可以不那么普遍地扩展或评论。

毫不奇怪，te()您安装的平滑度不能很好地适应假期；该模型假设周效应是平滑的，并且随着一年中的平滑一天效应平滑变化。假期不是从前一周或后一周的每周模式平稳偏离。假日效应并不能通过顺畅的关系很好地建模，需要其他东西来解释这种效应。

于 2018-01-20T17:59:21.280 回答
0

我的问题现在已经有几年了，我想添加最终对我最有效的解决方案。

首先，使用 mgcv 来拟合我的问题中描述的模型是不可能的。

有一段时间我使用了一个两阶段的过程。

model1 = gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc')
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)
# get_weekday_offset gets the coefficients for each weekday and normalizes them to have mean 0
data$weekday_offset = get_weekday_offset(model1)[data$day_in_week]

model2 = gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=52, bs='cc')
  + s(day_in_year, k=52, bs='cc', by=weekday_offset)
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)

第二个模型的形式是“y ~ f(年中的天) + g(年中的天) * h(周中的天)”，但h是一个固定函数，不适合第二个模型中的数据，而是基于最适合model1。这样做的明显缺点是 mgcv 需要运行两次，并且第一次运行的几乎所有信息（工作日的系数除外）都被丢弃。

我最终放弃这个模型的原因是除了振幅的变化之外，一年中形状确实发生了一些变化，所以我最终决定看起来像这样。

gam(
  y ~ s(day_in_year, k=52, bs='cc') 
  + s(day_in_year, k=10, bs='cc', by=as.factor(day_in_week)
  + as.factor(day_in_week)
  , knots=list(
    day_in_year=c(0, 366)
    , day_in_week=c(0,7)
  )
  , data = data
)

我放弃了s工作日的数据，因为我已经证明使用每周模式的最大自由度来证明数据是合理的，并且我在工作日和一年中的某一天（s(day_in_year, k=12, bs='cc', by=as.factor(day_in_week)）之间添加了交互效应，这估计了一条单独的季节性曲线对于每个工作日，但是由于它包含的自由度数量远低于季节性术语s(day_in_year, k=52, bs='cc')，因此我在尝试te术语时不会遇到问题。

于 2020-06-09T12:59:15.907 回答