如何测试GAMM中随机效应的统计显著性？

Question

我现在正在使用R中的包mgcv来构建GAMM模型，我的问题如下：

• 首先，我们如何知道随机效应是否具有统计学意义？
• 其次，我们如何从模型中提取随机截距值？
• 第三，gamm中的"offset"是什么意思？我已经查看了R帮助文档，但我仍然对该函数中的"offset"术语感到困惑。感谢任何帮助。

本例摘自书籍《广义加性模型：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))

- Yang Yang

---

som 和 str(som) 可能会显示它所包含的一些信息。 - Henry

@Henry 是的，som和str(som)可以显示一些信息，但它们无法展示随机效应的统计学意义。 - Yang Yang

从 http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-significance-of-random-effects 开始，对完整模型和简化模型进行似然比检验。 - Ben Bolker

@BenBolker您好，感谢您的帮助。您能否正式回答一下问题，这样我就可以接受您的答案了吗？同时，其他人也可以受益。谢谢。 - Yang Yang

Answer 1

请注意，对于此模型，使用gam()和bam()通过tw家族使用Tweedie响应可能更有意义，而不能使用gamm()。事实上，Simon Wood和Matteo Fasiolo使用位置尺度Tweedie GAM（其中他们单独用线性预测器[模型]建模Tweedie分布的平均值，方差和功率参数）适配这些数据。

根据@BenBolker的建议：我甚至不会特别在此模型中测试随机效应，而且我经常不关心它是否显着。这取决于我正在处理的问题或假设。通常，由于数据中存在一些聚类，我希望将其包含在模型中，而不管其显着性如何。

然而，我并不确信在这种情况下推广极大似然比检验（GLRT）理论是否适用于准似然方法。 Simon Wood在他有关GAMS的教材第二版的附录A中提出导数，显示如果我们用准似然函数替换对数似然函数，则先前推导出来的最大似然估计结果（其中包括GLRT的结果）仍然成立。至少Simon似乎是这样认为的，这表明我之前提到的测试的解释以及summary.gam()中实现的随机效应是可靠的，就像它基于适当的似然一样。

除非我真的需要，否则我会使用gam()或bam()来拟合此模型，然后再使用gamm()（后者来自 gamm4包），特别是对于非高斯模型，因为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所建议的那样，但参考分布已经校正，以便在参数空间边界进行测试。

> 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()模型，但无法使没有随机效应的模型收敛，因此无法测试随机效应项，甚至在尝试anova()的多模型形式时遇到错误。

如果您想获取随机效应，可以使用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页，了解更多关于此模型正在做什么以及偏移量的含义。

更一般地说，假设您用两个网捕捉河流；一个网的面积是另一个网的两倍。您更有可能在较大的网中捕获更多东西，因此较大网的计数将更高，因为采样工作更多-使用较大的网扫过了更多的河流（假设您采样相同的时间）。为确保您考虑到这种努力差异，您可以在模型中包含偏移量。对于具有对数链接的泊松模型，偏移量将为offset(log(net_area))。我们必须在链接比例上包含偏移量，因此使用log()。现在我们正在建模每个单位面积网的计数。

- Gavin Simpson

---

你好，感谢你的帮助。我想知道 est 是什么意思？我该如何使用你的软件包来解释随机效应的重要性呢？谢谢。 - Yang Yang

“est”列是每个站点的估计效应，“se”列是该效应的标准误。我已经展示了如何测试随机效应-查看summary()输出中的s（station）行。如果您想讨论“est”列中的单个效应，对于GAM（M），我不知道除了简单地在“est”上形成95％置信区间（从“est”减去和加上2倍的“se”）之外还有其他方法，因为非参数引导法对平滑器来说不可取。 - Gavin Simpson

@GavinSimpson，你是否担心对准似然模型进行似然比检验？它是欠散的，因此很难转换为标准的似然模型（例如，负二项式不是一个选项；我们需要广义伽玛或COM-Poisson或其他一些奇特的东西...）。我可能会更努力地说服OP，也许他们不需要对随机效应进行显著性检验... - Ben Bolker

@BenBolker 是的；我更喜欢用其他方式建模 - 我在最后一句话中的评论被推到了一边，因为我回答了其他部分。Simon Wood和Matteo Fasiolo在最近的工作中实际上使用Tweedie分布对这些数据进行建模（包括通过位置比例Tweedie进行建模，但目前mgcv仅适用于他们扩展的Fellner-Schall方法进行平滑选择，这不是默认值）。我个人的偏好是只包括随机效应而不费心测试它。 - Gavin Simpson

@BenBolker 我想看看Simon对此的看法，他在使用使用对数准似然的GLRT方面不太含蓄。我已经在答案中添加了一些内容，希望能解决其中的一些问题(?) - Gavin Simpson