使用 mgcv gam 的负二项式混合模型中固定效应的意义答案

【问题标题】：Significance of fixed effects in negative binomial mixed model using mgcv gam使用 mgcv gam 的负二项式混合模型中固定效应的意义
【发布时间】：2018-10-19 20:10:21
【问题描述】：

我正在使用 mgcv 包中的 gam 来分析包含 24 个条目的数据集：

ran  f1     f2   y
1   3000    5   545
1   3000    10  1045
1   10000   5   536
1   10000   10  770
2   3000    5   842
2   3000    10  2042
2   10000   5   615
2   10000   10  1361
3   3000    5   328
3   3000    10  1028
3   10000   5   262
3   10000   10  722
4   3000    5   349
4   3000    10  665
4   10000   5   255
4   10000   10  470
5   3000    5   680
5   3000    10  1510
5   10000   5   499
5   10000   10  1422
6   3000    5   628
6   3000    10  2062
6   10000   5   499
6   10000   10  2158

数据有两种固定效应（f1 和f2）和一种随机效应（ran）。相关数据为y。因为依赖数据y 表示计数并且过度分散，所以我使用的是负二项式模型。

gam 模型及其summary 输出如下：

library(mgcv)
summary(gam(y ~ f1 * f2 + s(ran, bs = "re"), data = df2, family = nb, method = "REML"))

Family: Negative Binomial(27.376) 
Link function: log 

Formula:
y ~ f1 * f2 + s(ran, bs = "re")

Parametric coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.500e+00  3.137e-01  17.533  < 2e-16 ***
f1          -3.421e-05  3.619e-05  -0.945    0.345    
f2           1.760e-01  3.355e-02   5.247 1.55e-07 ***
f1:f2        2.665e-07  4.554e-06   0.059    0.953    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
         edf Ref.df Chi.sq p-value    
s(ran) 4.726      5  85.66  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.866   Deviance explained = 93.6%
-REML = 185.96  Scale est. = 1         n = 24

来自summary 的 Wald 测试对 f2 (P = 1.55e-07) 具有非常高的意义。然而，当我通过使用方差分析比较两个不同的模型来测试f2 的重要性时，我得到了截然不同的结果：

anova(gam(y ~ f1 * f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
gam(y ~ f1 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
test="Chisq")

Analysis of Deviance Table

Model 1: y ~ f1 * f2 + s(ran, bs = "re")
Model 2: y ~ f1 + s(ran, bs = "re")
  Resid. Df Resid. Dev      Df Deviance Pr(>Chi)
1    14.843     18.340                          
2    16.652     21.529 -1.8091   -3.188   0.1752

f2 不再重要。模型从 REML 更改为 ML，建议用于评估固定效应。

如果交互作用被保留，使用方差分析 f2 仍然是微不足道的：

anova(gam(y ~ f1 + f2 + f1:f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
gam(y ~ f1 + f1:f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
test="Chisq")
Analysis of Deviance Table

Model 1: y ~ f1 + f2 + f1:f2 + s(ran, bs = "re")
Model 2: y ~ f1 + f1:f2 + s(ran, bs = "re")
  Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
1    14.843     18.340                           
2    15.645     19.194 -0.80159 -0.85391   0.2855

如果您能就这些方法中的哪一种更合适提供建议，我将不胜感激。非常感谢！

【问题讨论】：

标签： anova mixed-models gam mgcv

【解决方案1】：

?anova.gam 的警告部分说：

如果模型 a 和 b 仅在没有未惩罚的组件（例如随机效应）方面有所不同，则来自 anova(a,b) 的 p 值是不可靠的，并且通常太低了。

我猜 p 值是不可靠的，但在这种情况下，您会观察到与预期相反的情况 - p 值要大得多。

但是，我认为您没有比较正确的模型。除非您知道自己在做什么，否则在将模型与交互进行比较时应遵守边际原则。

因此，我会将具有 f1 和 f2 主效应的模型与包含这些主效应和它们的相互作用的模型进行比较。

型号1：y ~ f1 * f2 + s(ran, bs = "re")
型号2：y ~ f1 + f2 + s(ran, bs = "re")

除非您没有告诉我们您的模型是如何设置的，否则您不应在不包含在高阶项中找到的低阶项的情况下包含高阶项。例如，您有 f1 + f1:f2 和 f2 在二阶项中找到，但在模型中未作为一阶项找到。

【讨论】：

感谢您的回答 Gavin，也感谢您将 mgcv::gam 之火的礼物带给大众的所有努力。我阅读了边际原则，并决定最好使用 multcomp::glht 在 f2 的两个级别和 f1 的两个级别上对 f1 进行事后测试。这样，f1 有点显着（两个水平上 P ~ 0.04），f2 高度显着（P ~ 1e-15 两个水平），交互不显着（P = 0.96）。
另外，我接受了你的回答，但没有足够的 Stack Overflow 能力来显示绿色小箭头。
@Bob。我认为您添加复选标记的能力是一次限制，而不是重复限制。很确定你现在可以做到。（但也许也不会添加赞成票。）
非常感谢您的意见，@42-（我以前曾有幸享受并从中受益）。正如你所建议的，我现在可以添加漂亮的绿色复选标记（但不是赞成票），并且确实这样做了。再次感谢！