如何检验 GAMM 中随机效应的统计显着性？
How to test the statistical significance of a random effect in GAMM?
我现在正在使用包 mgcv 在 R 中构建一个 GAMM，我的问题是：
- 首先，我们如何知道随机效应是否具有统计显着性？
- 其次，我们如何提取模型中的随机截距值？
- 三、
gamm中的"offset"是什么意思？我查看了 R 帮助，但我仍然对函数中的 "offset" 术语感到困惑？感谢您的帮助。
例子取自书本Generalized additive models: an introduction with R
library(mgcv)
library(gamair)
data(sole)
sole$off <- log(sole$a.1-sole$a.0)
sole$a<-(sole$a.1+sole$a.0)/2
solr<-sole
solr$t<-solr$t-mean(sole$t)
solr$t<-solr$t/var(sole$t)^0.5
solr$la<-solr$la-mean(sole$la)
solr$lo<-solr$lo-mean(sole$lo)
solr$station <- factor(with(solr,paste(-la,-lo,-t,sep="")))
som <- gamm(eggs~te(lo,la,t,bs=c("tp","tp"),k=c(25,5),d=c(2,1))
+s(t,k=5,by=a)+offset(off), family=quasipoisson,
data=solr,random=list(station=~1))
请注意，对于此模型，通过家庭 tw 与 gam() 和 bam() 使用 Tweedie 响应可能更有意义，不能与 [=16] 一起使用=].事实上，Simon Wood 和 Matteo Fasiolo 将这些数据与位置尺度 Tweedie GAM 进行拟合（其中他们使用单独的线性预测器 [模型] 对 Tweedie 分布的均值、方差和功率参数进行建模）。
按照@BenBolker 的建议：我什至懒得专门测试这个模型中的随机效应，而且通常我不关心它是否显着。这取决于我正在研究的问题或假设。由于我希望包含在模型中的数据中存在一些聚类，而不管重要性如何，我经常希望它出现在模型中。
但是，我不相信（广义）似然比检验 (GLRT) 的理论不适用于准似然的使用在这种情况下 . Simon Wood 在他的 GAMS 教科书第 2 版附录 A 中介绍了推导，表明如果我们用对数准似然替换对数似然，则先前推导的最大似然估计结果（包括 GLRT 的结果）仍然成立。至少 Simon 似乎在争论，这表明我在下面提到的测试的解释以及在 summary.gam() 中针对随机效应实施的测试的解释与基于适当可能性的测试一样可靠。
除非我真的需要，否则我会用 gam() 或 bam() 和 gamm4()（后者来自 gamm4 package), before gamm(), 特别是对于非高斯模型，因为 gamm() 函数必须将此模型拟合为使用惩罚拟似然的混合效应模型，这不一定是最好的方法估计这些模型。
library(mgcv)
library(gamair)
devtools::install_github('gavinsimpson/gratia')
library(gratia)
data(sole)
sole$off <- log(sole$a.1-sole$a.0)
sole$a<-(sole$a.1+sole$a.0)/2
solr <- sole
solr$t <- solr$t-mean(sole$t)
solr$t <- solr$t/var(sole$t)^0.5
solr$la <- solr$la-mean(sole$la)
solr$lo <- solr$lo-mean(sole$lo)
solr$station <- factor(with(solr,paste(-la,-lo,-t,sep="")))
som <- gam(eggs ~ te(lo, la, t, bs = c('tp','tp'), k = c(25, 5), d = c(2,1)) +
s(t, k = 5, by = a) + s(station, bs = 're') + offset(off),
family = quasipoisson, data = solr, method = 'REML')
然后 summary(som) 按照@BenBolker 的建议给出了基于似然比检验的检验，但是为了在参数 space.[=30= 的边界上检验参考分布被修正]
> summary(som)
Family: quasipoisson
Link function: log
Formula:
eggs ~ te(lo, la, t, bs = c("tp", "tp"), k = c(25, 5), d = c(2,
1)) + s(t, k = 5, by = a) + s(station, bs = "re") + offset(off)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4016 0.3061 -11.11 <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
te(lo,la,t) 56.025 65.456 2.547 4.62e-10 ***
s(t):a 4.535 4.886 54.790 < 2e-16 ***
s(station) 128.563 388.000 1.175 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.833 Deviance explained = 88%
-REML = -7.9014 Scale est. = 0.58148 n = 1575
我无法使用 gamm() 让没有随机效应的模型收敛，所以我无法测试随机效应项，甚至在尝试 [=25 的多模型形式时遇到错误=].
如果你想获得随机效果，使用 gam() 模型，你可以使用我的 gratia 包（希望几天后在 CRAN 上，但可以如上所示从 github 安装），然后：
> evaluate_smooth(som, 's(station)')
# A tibble: 394 x 5
smooth by_variable station est se
<chr> <fct> <chr> <dbl> <dbl>
1 s(statio… NA -0.0004304761904734280.419685714285714-… -0.0396 2.55
2 s(statio… NA -0.0004304761904734280.6586857142857140… 1.48 1.20
3 s(statio… NA -0.0004304761904734281.15968571428571-1… -0.00606 2.63
4 s(statio… NA -0.0004304761904734281.176685714285710.… -0.0767 2.48
5 s(statio… NA -0.002430476190475870.9096857142857141.… -0.00654 2.63
6 s(statio… NA -0.01243047619047390.4106857142857140.0… -0.802 1.61
7 s(statio… NA -0.0154304761904740.631685714285714-0.4… -0.138 2.35
8 s(statio… NA -0.02043047619047660.375685714285714-0.… -0.426 1.94
9 s(statio… NA -0.02543047619047911.14668571428571-0.4… -0.0333 2.57
10 s(statio… NA -0.02743047619047450.875685714285714-0.… -0.0673 2.49
# … with 384 more rows
并且您想要 est 列。
偏移量是模型中的一个术语，其固定效应为 1。在本例中，它用于标准化计数响应，以便您可以对每个对象进行同类比较；在这种情况下，它被用来整合在这个样本中发现的鸡蛋的年龄。阅读第Simon 的 GAM 书第 2 版的 143，以了解更多有关此模型所做的工作以及偏移量的含义。
更一般地说，假设您用两张网对河流进行采样；一张网的面积是另一张的两倍。你更有可能在更大的网中捕获更多的东西，因此由于更大的采样努力，更大的网中的计数会更高——你用更大的网扫了更多的河流（假设你采样了相同数量的水）时间）。为确保考虑到这种工作量差异，您可以在模型中包含一个偏移量。偏移量将为（对于具有 log link 的泊松模型）offset(log(net_area))。我们必须在 link 尺度上包括偏移量，因此 log()。现在我们建模的是单位面积网的数量。
我现在正在使用包 mgcv 在 R 中构建一个 GAMM，我的问题是：
- 首先，我们如何知道随机效应是否具有统计显着性？
- 其次，我们如何提取模型中的随机截距值？
- 三、
gamm中的"offset"是什么意思？我查看了 R 帮助，但我仍然对函数中的 "offset" 术语感到困惑？感谢您的帮助。
例子取自书本Generalized additive models: an introduction with R
library(mgcv)
library(gamair)
data(sole)
sole$off <- log(sole$a.1-sole$a.0)
sole$a<-(sole$a.1+sole$a.0)/2
solr<-sole
solr$t<-solr$t-mean(sole$t)
solr$t<-solr$t/var(sole$t)^0.5
solr$la<-solr$la-mean(sole$la)
solr$lo<-solr$lo-mean(sole$lo)
solr$station <- factor(with(solr,paste(-la,-lo,-t,sep="")))
som <- gamm(eggs~te(lo,la,t,bs=c("tp","tp"),k=c(25,5),d=c(2,1))
+s(t,k=5,by=a)+offset(off), family=quasipoisson,
data=solr,random=list(station=~1))
请注意，对于此模型，通过家庭 tw 与 gam() 和 bam() 使用 Tweedie 响应可能更有意义，不能与 [=16] 一起使用=].事实上，Simon Wood 和 Matteo Fasiolo 将这些数据与位置尺度 Tweedie GAM 进行拟合（其中他们使用单独的线性预测器 [模型] 对 Tweedie 分布的均值、方差和功率参数进行建模）。
按照@BenBolker 的建议：我什至懒得专门测试这个模型中的随机效应，而且通常我不关心它是否显着。这取决于我正在研究的问题或假设。由于我希望包含在模型中的数据中存在一些聚类，而不管重要性如何，我经常希望它出现在模型中。
但是，我不相信（广义）似然比检验 (GLRT) 的理论不适用于准似然的使用在这种情况下 . Simon Wood 在他的 GAMS 教科书第 2 版附录 A 中介绍了推导，表明如果我们用对数准似然替换对数似然，则先前推导的最大似然估计结果（包括 GLRT 的结果）仍然成立。至少 Simon 似乎在争论，这表明我在下面提到的测试的解释以及在 summary.gam() 中针对随机效应实施的测试的解释与基于适当可能性的测试一样可靠。
除非我真的需要，否则我会用 gam() 或 bam() 和 gamm4()（后者来自 gamm4 package), before gamm(), 特别是对于非高斯模型，因为 gamm() 函数必须将此模型拟合为使用惩罚拟似然的混合效应模型，这不一定是最好的方法估计这些模型。
library(mgcv)
library(gamair)
devtools::install_github('gavinsimpson/gratia')
library(gratia)
data(sole)
sole$off <- log(sole$a.1-sole$a.0)
sole$a<-(sole$a.1+sole$a.0)/2
solr <- sole
solr$t <- solr$t-mean(sole$t)
solr$t <- solr$t/var(sole$t)^0.5
solr$la <- solr$la-mean(sole$la)
solr$lo <- solr$lo-mean(sole$lo)
solr$station <- factor(with(solr,paste(-la,-lo,-t,sep="")))
som <- gam(eggs ~ te(lo, la, t, bs = c('tp','tp'), k = c(25, 5), d = c(2,1)) +
s(t, k = 5, by = a) + s(station, bs = 're') + offset(off),
family = quasipoisson, data = solr, method = 'REML')
然后 summary(som) 按照@BenBolker 的建议给出了基于似然比检验的检验，但是为了在参数 space.[=30= 的边界上检验参考分布被修正]
> summary(som)
Family: quasipoisson
Link function: log
Formula:
eggs ~ te(lo, la, t, bs = c("tp", "tp"), k = c(25, 5), d = c(2,
1)) + s(t, k = 5, by = a) + s(station, bs = "re") + offset(off)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4016 0.3061 -11.11 <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
te(lo,la,t) 56.025 65.456 2.547 4.62e-10 ***
s(t):a 4.535 4.886 54.790 < 2e-16 ***
s(station) 128.563 388.000 1.175 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.833 Deviance explained = 88%
-REML = -7.9014 Scale est. = 0.58148 n = 1575
我无法使用 gamm() 让没有随机效应的模型收敛，所以我无法测试随机效应项，甚至在尝试 [=25 的多模型形式时遇到错误=].
如果你想获得随机效果，使用 gam() 模型，你可以使用我的 gratia 包（希望几天后在 CRAN 上，但可以如上所示从 github 安装），然后：
> evaluate_smooth(som, 's(station)')
# A tibble: 394 x 5
smooth by_variable station est se
<chr> <fct> <chr> <dbl> <dbl>
1 s(statio… NA -0.0004304761904734280.419685714285714-… -0.0396 2.55
2 s(statio… NA -0.0004304761904734280.6586857142857140… 1.48 1.20
3 s(statio… NA -0.0004304761904734281.15968571428571-1… -0.00606 2.63
4 s(statio… NA -0.0004304761904734281.176685714285710.… -0.0767 2.48
5 s(statio… NA -0.002430476190475870.9096857142857141.… -0.00654 2.63
6 s(statio… NA -0.01243047619047390.4106857142857140.0… -0.802 1.61
7 s(statio… NA -0.0154304761904740.631685714285714-0.4… -0.138 2.35
8 s(statio… NA -0.02043047619047660.375685714285714-0.… -0.426 1.94
9 s(statio… NA -0.02543047619047911.14668571428571-0.4… -0.0333 2.57
10 s(statio… NA -0.02743047619047450.875685714285714-0.… -0.0673 2.49
# … with 384 more rows
并且您想要 est 列。
偏移量是模型中的一个术语，其固定效应为 1。在本例中，它用于标准化计数响应，以便您可以对每个对象进行同类比较；在这种情况下，它被用来整合在这个样本中发现的鸡蛋的年龄。阅读第Simon 的 GAM 书第 2 版的 143，以了解更多有关此模型所做的工作以及偏移量的含义。
更一般地说，假设您用两张网对河流进行采样；一张网的面积是另一张的两倍。你更有可能在更大的网中捕获更多的东西，因此由于更大的采样努力，更大的网中的计数会更高——你用更大的网扫了更多的河流（假设你采样了相同数量的水）时间）。为确保考虑到这种工作量差异，您可以在模型中包含一个偏移量。偏移量将为（对于具有 log link 的泊松模型）offset(log(net_area))。我们必须在 link 尺度上包括偏移量，因此 log()。现在我们建模的是单位面积网的数量。