uni*_*ue2 7 regression r gam mgcv
我使用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')哪种方法有效,但引入了太多的自由度,因为模型过度适应假期,这些假期在特定工作日内可用的短数量.
是否可以按照我想要的方式指定模型？
哇,这是一个很老的问题!
虽然每周模式本身在一年中非常稳定,但其幅度也随季节而变化.
使用张量积样条基te是交互的正确方法,尽管更合适的构造函数是ti.你说过这te会给你带来很多参数.当然.你有k = 52第一个边距,而k = 5第二个,你最终52 * 5 - 1得到这个张量项的系数.但这只是创建交互的方式.
请注意,在mgcvGAM公式中,:或*仅对参数项之间的交互有效.平滑之间的相互作用由te或处理ti.
如果这不是您所希望的,那么您对"产品"的期望是什么？两个边缘设计矩阵的Hadamard积？那有什么意义呢？顺便说一下,Hadamard产品需要两个相同尺寸的矩阵.但是,您的两个边距没有相同的列数.
如果你不明白为什么我一直在谈论矩阵,那么你需要在2006年阅读Simon的书.尽管GAM估计现在已经过时了,但GAM的构造/设置(如设计矩阵和惩罚矩阵)在第即使在十年之后,4也完全没有变化.
好的,让我再给你一个提示.用于构造/ 设计矩阵的行式Kronecker产品不是新发明.teti
平滑项s(x)非常类似于因子变量g,就好像它们看起来是单个变量一样,它们被构造成具有许多列的设计矩阵.因为g它是一个虚拟矩阵,而f(x)它是一个基矩阵.因此,两个平滑函数之间的相互作用的构造方式与构造两个因子之间的相互作用的方式相同.
如果你有一个g15级的因子,另一个g210级的因子,它们的边际设计矩阵(在对比之后)有4列和9列,那么交互g1:g2将有36列.这样的设计矩阵,只是行和Kronecker产品的设计矩阵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)吗？
您声明的模型并没有多大意义：
这可以完全表示为张量积平滑。您在此处对其他答案的评论中提到的模型
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()您安装的平滑度不能很好地适应假期。该模型假设周效应是平滑的，并且随着年中日平滑的变化而变化。假期并不是上周或下周每周模式的平稳变化。该节日效应没有很好地平滑关系建模和别的东西需要考虑到这种效果。