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

Question

我正在使用 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

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

Answer 1

?anova.gam 的警告部分说：

If models a and b differ only in terms with no un-penalized components (such as random effects) then p values from anova(a,b) are unreliable, and usually much too low.

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

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

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

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

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

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

Significance of fixed effects in negative binomial mixed model using mgcv gam

anova

gam

mixed-models

mgcv