模型构建:如何构建有意义的游戏模型?(广义相加模型)

机器算法验证 r 模型选择 广义加法模型
2022-03-11 22:42:59

我看到有关于游戏的解释和构建的各种问题,这似乎说明了非统计学家处理这些问题的困难。不幸的是,从我阅读的任何线程或教程中,我都无法清楚地了解如何构建有意义的模型。

目前我正在研究有机农业对蜂群性能的影响。因此,我尝试将景观特征(例如半径 500m (bio.percent_b500) 内的有机农业百分比)与蜂群发展参数(例如蜂蜜储备)联系起来。我首先建立了一个基本的 gam 模型 (model0),只有一年中的一周作为解释变量,因为蜂箱中的蜂蜜量在一年中呈非线性变化。

library("gam")
library("mgcv")

model0 <- gam(honey.mean ~ s(week), data= my.data.frame) 
summary(model0)
plot(model0)

在此处输入图像描述

然后我试图包含一个包含有机农业百分比的平滑术语。然而,这失败了,我猜是因为超过 85% 的菌落在 500m 半径范围内没有有机田。

model1 <- gam(honey.mean ~ s(week) + s(bio.percent_b500),data = my.data.frame)
# Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : 
# A term has fewer unique covariate combinations than specified maximum 
# degrees of freedom

model2 = gam(honey.mean ~ s(week,bio.percent_b500) , data= my.data.frame)

然后我惊讶地发现,包含有机农业百分比和周数相互作用的模型奏效了。然而,我在一本德国统计书中读到,交互项不应包含在没有独立影响的模型中。作者提到了一个叫做“Marginalitätstheorem”(边际定理)的东西。因为我知道,从模型 1 中有机农业的平滑术语会导致问题,所以我只在一年中的一周中加入了一个额外的平滑术语。这个模型对我来说很直观,因为一年中的一周总是有效的。然而,有机农业的效果总是取决于一年中的时间。例如在夏天应该有更高的杂草花可用性。

model3 = gam(honey.mean ~ s(week) + s(week, bio.percent_b500) , data= my.data.frame)      

由于蜂巢中的蜂蜜储量可能取决于各种景观特征,因此我建立了模型,包括油菜的百分比 (osr.percent_b500)。

model4 = gam(honey.mean ~ s(week) + s(osr.percent_b500),data = my.data.frame)
vis.gam(model4, type = "response", plot.type = "persp")   
summary(model4)

model5 = gam(honey.mean ~ s(week,osr.percent_b500) + s(week,bio.percent_b500), data = my.data.frame)
summary(model5)

model6 = gam(honey.mean ~ s(week) + s(week,osr.percent_b500) + s(week,bio.percent_b500), data= my.data.frame)
summary(model6)

model7 = gam(honey.mean ~ s(week) + s(week,osr.percent_b500,bio.percent_b500), data= my.data.frame)
summary(model7)

由于上述原因,模型 0、3 和 6 在我看来最有意义。我不确定我是否应该考虑以另一种方式构建的模型,并通过 AIC 接受和比较它们。

AIC(model0,model2,model3,model4,model5,model6,model7)

AIC 值的比较确定模型 7 是最好的,因为它的模型自由度比模型 3 少。这让我再次感到惊讶,因为模型 7 包含比模型 3 更复杂的交互。

谁能给我关于如何构建有意义的游戏模型的建议?

1)交互项可以在没有独立项的情况下出现在(gam)模型中吗?

2)为什么更复杂的游戏交互平滑项会导致模型自由度降低?

3)上面提到的哪些模型是有意义的?

4)对于我正在尝试做的事情,是否有更好的替代广义加法模型的方法?

您可以在下面找到 my.data.frame:

structure(list(year = c(2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L), apiary = c(4L, 8L, 8L, 8L, 18L, 18L, 18L, 19L, 
19L, 19L, 23L, 23L, 23L, 23L, 34L, 34L, 34L, 45L, 45L, 45L, 46L, 
46L, 46L, 49L, 49L, 49L, 3L, 3L, 3L, 3L, 9L, 9L, 9L, 9L, 14L, 
14L, 14L, 14L, 17L, 17L, 17L, 17L, 20L, 20L, 20L, 28L, 28L, 28L, 
28L, 31L, 31L, 31L, 31L, 33L, 33L, 33L, 33L, 33L, 35L, 35L, 35L, 
44L, 44L, 44L, 44L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
12L, 12L, 26L, 26L, 26L, 26L, 26L, 30L, 30L, 30L, 30L, 30L, 32L, 
32L, 32L, 32L, 32L, 37L, 37L, 37L, 37L, 37L, 42L, 42L, 42L, 42L, 
42L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 48L, 48L, 48L, 48L, 48L, 
50L, 50L, 50L, 50L, 1L, 1L, 1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 6L, 
7L, 7L, 7L, 7L, 7L, 7L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 24L, 
24L, 24L, 27L, 27L, 27L, 27L, 27L, 27L, 36L, 36L, 36L, 36L, 36L, 
40L, 40L, 40L, 40L, 40L, 41L, 41L, 41L, 41L, 41L, 43L, 43L, 43L, 
43L, 43L, 43L, 43L, 2L, 2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 
10L, 10L, 10L, 10L, 10L, 10L, 13L, 13L, 13L, 13L, 13L, 15L, 15L, 
15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 21L, 21L, 21L, 
21L, 21L, 21L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 29L, 29L, 29L, 
29L, 29L, 29L, 29L, 39L, 39L, 39L, 39L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 18L, 18L, 18L, 18L, 18L, 18L, 
18L, 19L, 19L, 19L, 19L, 19L, 19L, 23L, 23L, 23L, 23L, 23L, 23L, 
23L, 34L, 34L, 34L, 34L, 34L, 34L, 38L, 38L, 38L, 38L, 38L, 38L, 
38L, 45L, 45L, 45L, 45L, 45L, 45L, 46L, 46L, 46L, 46L, 46L, 46L, 
46L, 49L, 49L, 49L, 49L, 49L, 49L), week = c(26L, 24L, 26L, 28L, 
23L, 28L, 31L, 23L, 24L, 28L, 24L, 26L, 28L, 29L, 23L, 26L, 28L, 
24L, 26L, 29L, 23L, 28L, 29L, 23L, 28L, 31L, 18L, 20L, 22L, 32L, 
18L, 20L, 30L, 32L, 16L, 22L, 26L, 32L, 16L, 18L, 24L, 28L, 16L, 
24L, 32L, 16L, 24L, 28L, 30L, 18L, 20L, 22L, 26L, 16L, 20L, 22L, 
26L, 30L, 16L, 24L, 28L, 18L, 26L, 28L, 32L, 20L, 21L, 33L, 35L, 
39L, 21L, 25L, 27L, 29L, 31L, 35L, 21L, 25L, 27L, 31L, 35L, 21L, 
23L, 29L, 35L, 39L, 17L, 27L, 33L, 35L, 39L, 17L, 20L, 27L, 35L, 
39L, 17L, 21L, 23L, 25L, 35L, 17L, 20L, 21L, 25L, 27L, 31L, 33L, 
17L, 21L, 23L, 29L, 39L, 20L, 31L, 33L, 39L, 19L, 21L, 23L, 29L, 
37L, 19L, 21L, 23L, 29L, 33L, 39L, 17L, 19L, 25L, 29L, 31L, 35L, 
19L, 33L, 37L, 39L, 15L, 19L, 23L, 35L, 37L, 39L, 15L, 17L, 21L, 
29L, 33L, 35L, 17L, 23L, 25L, 29L, 39L, 17L, 19L, 21L, 29L, 35L, 
17L, 19L, 21L, 25L, 39L, 15L, 19L, 27L, 31L, 33L, 37L, 39L, 13L, 
23L, 27L, 33L, 35L, 39L, 23L, 25L, 27L, 31L, 37L, 13L, 15L, 19L, 
23L, 29L, 37L, 29L, 33L, 35L, 37L, 39L, 13L, 21L, 25L, 27L, 29L, 
35L, 23L, 29L, 31L, 35L, 37L, 39L, 15L, 19L, 21L, 27L, 33L, 39L, 
13L, 15L, 23L, 27L, 29L, 35L, 39L, 13L, 15L, 23L, 27L, 29L, 31L, 
35L, 13L, 31L, 35L, 37L, 16L, 20L, 26L, 38L, 40L, 42L, 44L, 16L, 
24L, 32L, 34L, 38L, 40L, 44L, 18L, 20L, 24L, 34L, 38L, 42L, 44L, 
24L, 28L, 32L, 40L, 42L, 44L, 16L, 20L, 26L, 38L, 40L, 42L, 44L, 
18L, 20L, 22L, 32L, 38L, 44L, 16L, 20L, 22L, 28L, 30L, 34L, 38L, 
18L, 20L, 22L, 28L, 32L, 44L, 16L, 22L, 24L, 28L, 32L, 34L, 38L, 
22L, 28L, 32L, 34L, 38L, 40L), bio.percent_b500 = c(0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.13, 16.13, 16.13, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.73, 15.73, 15.73, 
15.73, 15.73, 0, 0, 0, 0, 0, 0.75, 0.75, 0.75, 0.75, 0.75, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 2.14, 2.14, 2.14, 2.14, 2.14, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.69, 13.69, 13.69, 13.69, 
13.69, 13.69, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 6.47, 6.47, 6.47, 6.47, 6.47, 6.47, 6.47, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.68, 5.68, 
5.68, 5.68, 5.68, 5.68, 5.68, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 44.93, 44.93, 44.93, 44.93, 44.93, 
44.93, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0), osr.percent_b500 = c(10.12, 1.51, 1.51, 
1.51, 0, 0, 0, 4.85, 4.85, 4.85, 0, 0, 0, 0, 8.94, 8.94, 8.94, 
0, 0, 0, 0, 0, 0, 1.2, 1.2, 1.2, 6.41, 6.41, 6.41, 6.41, 0, 0, 
0, 0, 8.27, 8.27, 8.27, 8.27, 4.67, 4.67, 4.67, 4.67, 7.2, 7.2, 
7.2, 5.84, 5.84, 5.84, 5.84, 20.51, 20.51, 20.51, 20.51, 10.22, 
10.22, 10.22, 10.22, 10.22, 9.85, 9.85, 9.85, 0.02, 0.02, 0.02, 
0.02, 14.33, 14.33, 14.33, 14.33, 14.33, 21.6, 21.6, 21.6, 21.6, 
21.6, 21.6, 0, 0, 0, 0, 0, 6.1, 6.1, 6.1, 6.1, 6.1, 3.18, 3.18, 
3.18, 3.18, 3.18, 5.45, 5.45, 5.45, 5.45, 5.45, 0, 0, 0, 0, 0, 
22.65, 22.65, 22.65, 22.65, 22.65, 22.65, 22.65, 0.52, 0.52, 
0.52, 0.52, 0.52, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
5.59, 5.59, 5.59, 5.59, 5.59, 5.59, 7.41, 7.41, 7.41, 7.41, 4.13, 
4.13, 4.13, 4.13, 4.13, 4.13, 21.77, 21.77, 21.77, 21.77, 21.77, 
21.77, 3.58, 3.58, 3.58, 3.58, 3.58, 7.09, 7.09, 7.09, 7.09, 
7.09, 18.35, 18.35, 18.35, 18.35, 18.35, 0.78, 0.78, 0.78, 0.78, 
0.78, 0.78, 0.78, 0.41, 0.41, 0.41, 0.41, 0.41, 0.41, 12.2, 12.2, 
12.2, 12.2, 12.2, 0.26, 0.26, 0.26, 0.26, 0.26, 0.26, 7.57, 7.57, 
7.57, 7.57, 7.57, 12.8, 12.8, 12.8, 12.8, 12.8, 12.8, 34.1, 34.1, 
34.1, 34.1, 34.1, 34.1, 18.33, 18.33, 18.33, 18.33, 18.33, 18.33, 
12.44, 12.44, 12.44, 12.44, 12.44, 12.44, 12.44, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1.97, 1.97, 1.97, 1.97, 1.97, 1.97, 1.97, 
18.06, 18.06, 18.06, 18.06, 18.06, 18.06, 18.06, 0, 0, 0, 0, 
0, 0, 0, 16.76, 16.76, 16.76, 16.76, 16.76, 16.76, 0, 0, 0, 0, 
0, 0, 0, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 0, 0, 0, 0, 0, 0, 
0, 5.28, 5.28, 5.28, 5.28, 5.28, 5.28, 7.99, 7.99, 7.99, 7.99, 
7.99, 7.99, 7.99, 18.09, 18.09, 18.09, 18.09, 18.09, 18.09), 
    honey.mean = c(2.48, 3.99666666666667, 2.36, 2.94, 3.42, 
    3.71, 4.09, 2.12, 3.92, 4.145, 6.27, 6.92, 9.16, 6.75, 6.8, 
    1.07, 6.06, 1.7, 3.4, 5.805, 4.45, 4.19, 13.61, 3.695, 2.86, 
    8.32, 7.67, 6.81, 3.68, 14.335, 2.78, 3.62, 19.035, 12.77, 
    5.81, 3.05, 10.22, 10.44, 4.43, 8.64, 2.4, 16.41, 2.9, 7.175, 
    15.735, 3.16, 1.49, 5.48, 18.95, 6.885, 4.46, 7.9, 0.68, 
    1.4, 2.5, 8.12, 3.09, 14.72, 5.85, 1.885, 16.44, 8.055, 6.68, 
    8.58, 24.7, 8.135, 8.43, 26.08, 16.83, 9.72, 5.24, 5.65, 
    5.19, 7.35, 17.25, 8.82, 14.95, 12.05, 7.3, 62.4, 16.68, 
    1, 10.65, 10.28, 19.65, 17.26, 6.64, 9.94, 65.15, 12.07, 
    20.62, 7.7, 6.31, 1.68, 20.97, 23.825, 6.5, 6.14, 4.22, 2.47, 
    17.97, 2.61, 3.17, 3.24, 0.57, 0.54, 33.07, 49.8, 9.1, 8.41, 
    7.29, 10.61, 19.67, 3.09, 37.125, 24.99, 18.62, 24.15, 17.96, 
    16.61, 28.86, 7.74, 18.95, 18.45, 15.56, 48.35, 16.045, 8.37, 
    23.47, 5.44, 1.8, 64.27, 17.08, 20.62, 18.465, 18.255, 16.5, 
    23.17, 7.49, 12.55, 7.45, 16.72, 23.29, 7.965, 9.83, 15.39, 
    11.19, 35.85, 16.755, 18.8, 19.51, 10.39, 14.02, 32.82, 12.9466666666667, 
    14.68, 15.79, 12.8, 40.37, 22.27, 14.63, 16.9, 6.65, 2.42, 
    18.24, 9.3, 23.08, 17.94, 57.78, 24.34, 20.06, 18.2, 3.99, 
    6.465, 2.93, 25.98, 19.87, 17.25, 13.21, 9.07, 5.21, 9.48, 
    11.825, 7.58, 3.41, 12.56, 13.58, 22.17, 19.43, 11.7, 36.5, 
    18, 12.675, 5.8, 7.72, 4.41, 1.96, 2.83, 12.04, 17.24, 15.77, 
    17.655, 40.15, 21.87, 17.42, 19.16, 8.91, 5.41, 19.91, 9.65, 
    43.54, 17.72, 2.85, 3.41, 7.4, 7.38, 13.73, 14.16, 20.25, 
    2.77, 5.93, 11.185, 2.36, 12.62, 30.24, 13.97, 9.11, 13.985, 
    12.54, 11.13, 1.54, 8.91, 1.3, 4.03, 9.2, 8.86, 9.12, 1.11, 
    7.83, 17.985, 0.86, 14.5, 4.17, 5.18, 5.76, 6.22, 3.79, 17.18, 
    15.83, 11.195, 9.99, 12.395, 7.42, 26.15, 18.29, 15.955, 
    14.76, 2.18, 4.41, 3.53, 11.77, 10.1, 12.81, 20.25, 4.9, 
    10.43, 0.84, 8.81, 19.59, 24.94, 1.42, 6.57, 11.38, 1.92, 
    6.97, 19.31, 17.885, 8.07, 11.25, 6.05, 5.55, 30.23, 9.82, 
    4.8, 4.94, 3.835, 2.54, 21.73, 20.84, 19.02, 5.62, 0.72, 
    23.335, 10.745, 10.43, 7.34)), .Names = c("year", "apiary", 
"week", "bio.percent_b500", "osr.percent_b500", "honey.mean"), row.names = c(NA, 
296L), class = "data.frame")
1个回答

第一个错误的原因是bio.percent_b500没有k- 1 个唯一值。如果您k将此样条线设置为较低的值,则模型将适合。为什么二维薄板版本有效,我认为 IIRC 是因为它计算k未提供时的默认值的方式,它必须从数据中执行此操作。所以对于model1,设置k = 9有效:

> model1 <- gam(honey.mean ~ s(week) + s(bio.percent_b500, k = 9), data = df)

在 mgcv 中,模型s(x1, x2) 包含主效应交互作用。请注意,s()假设变量中的缩放比例相似,所以我怀疑您是否真的想要s()- 尝试te()使用张量积平滑。如果您想分解成主要影响和交互,请尝试:

model2 <- gam(honey.mean ~ ti(week) + ti(bio.percent_b500) + 
                ti(week, bio.percent_b500), data = df)

summary()它表明不需要交互:

> summary(model2)

Family: gaussian 
Link function: identity 

Formula:
honey.mean ~ ti(week) + ti(bio.percent_b500) + ti(week, bio.percent_b500)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.2910     0.5263   23.35   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                            edf Ref.df      F p-value    
ti(week)                  3.753  3.963 22.572  <2e-16 ***
ti(bio.percent_b500)      3.833  3.974  2.250  0.0461 *  
ti(week,bio.percent_b500) 1.246  1.448  1.299  0.2036    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.25   Deviance explained = 27.3%
GCV = 84.697  Scale est. = 81.884    n = 296

你从广义似然比检验中得到类似的信息”

> model1 <- gam(honey.mean ~ ti(week) + ti(bio.percent_b500, k = 9), data = df)
> anova(model1, model2, test = "LRT")
Analysis of Deviance Table

Model 1: honey.mean ~ ti(week) + ti(bio.percent_b500, k = 9)
Model 2: honey.mean ~ ti(week) + ti(bio.percent_b500) + ti(week, bio.percent_b500)
  Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
1    285.65      23363                           
2    286.17      23433 -0.51811  -69.675   0.1831

这里ti()的项是张量积交互项,其中模型中其他项的主要/边际效应被删除。你并不严格需要ti()我认为model1te()甚至s()应该工作,但mgcv中的示例现在都使用这种形式,所以我会这样做。

model3对我来说没有多大意义,特别是如果你使用张量产品;s(week)术语包含在 中s(week, bio.percent_b500),但正如我提到的双变量s()术语假设各向同性,因此这可能过度限制了组件,从而为 进来并解释某些东西week留下了空间。s(week)一般来说,您不应该这样做,因为您正确使用了双变量术语。我不清楚您是否可以使用 500m 变量获得真正的双变量项。

问题:

第一季度

我怀疑您是否想要没有主要/边际影响的互动。你的model2包括两者。

第二季度

就您设想的函数集而言,该模型可能更复杂,但您需要考虑为每个平滑项生成的基数,另外mgcv使用带有惩罚的样条曲线来进行平滑度选择,因此您最终可能会得到一个位于乍一看更复杂,但平滑项已被缩小,因为惩罚使得它们在拟合后使用更少的自由度。

第三季度

从技术上讲,只有model1真正正确指定。我认为您不想假设model2. 我认为您会遇到与model3,model6model7一般的拟合问题;如果您正确设置了张量积的基础,则不需要s(week)这些模型中的单独术语。