包括泊松模型中的偏移量是候选模型的两倍答案

【问题标题】：Including offset in poisson models double number of candidate models包括泊松模型中的偏移量是候选模型的两倍
【发布时间】：2020-06-12 21:36:10
【问题描述】：

修改后的问题和代码 (2020-07-10)

随着分析的进行，这个问题已经进行了多次修改。我在顶部包含了最新信息，但是对于那些在其他步骤中遇到问题的人，请参阅下面的信息。

回顾：我正在对几个连续变量（A、B 和 C）、一个分类变量 (YN) 和一个随机截距 (ID) 之间的关系进行建模。我已经包含了一个偏移量 (T)，因为计数数据是在不同持续时间的行为观察期间收集的。还有一个基于随机截距的零通胀术语。我正在比较该模型的所有子集以确定最适合的模型。

目前，当我尝试比较所有子集时，疏通函数会抛出错误。

#simulated data
df <- data.frame(
  count = sample(size = 300, c(0:6), replace = TRUE, prob = c(0.7, 0.15, 0.05, 0.04, 0.03, 0.02, 0.01)),
  A = sample(size = 300, x = 24:65, replace = TRUE),
  B = runif(n = 300, min = 7, max = 19),
  C = rbeta(n = 300, shape1 = 1, shape2 = 25)*100,
  YN = sample(size = 300, c("Y", "N"), replace = TRUE, prob = c(0.2, 0.8)),
  T = runif(n = 300, min = 10, max = 20),
  ID = sample(size = 300, letters, replace = TRUE)
)
  
#standardize explanatory variables
dfSTD <- stdize(df, 
                omit.cols = c("count", "YN", "T", "ID"), 
                center = TRUE, scale = TRUE) 

#specify full model with standardized covariates for analysis
full <- stdizeFit(glmmTMB(count ~ z.A + z.B + z.C + YN + offset(log(T)) + (1|ID),
                          zi = ~ (1|ID),
                          family = poisson(),
                          data = dfSTD,
                          na.action = "na.fail"),
                  data = dfSTD, 
                  which = "formula", 
                  evaluate = FALSE)

#compare all possible model subsets
all <- dredge(full, 
              beta = "sd", 
              fixed = "cond(offset(log(T)))")

Error in nobs.default(global.model) : no 'nobs' method is available

我正在使用 R 版本 3.6.1。

我正在对几个连续变量（A、B 和 C）、一个分类变量 (YN) 和一个随机截距 (ID) 之间的关系进行建模。我已经包含了一个偏移量 (T)，因为计数数据是在不同持续时间的行为观察期间收集的。还有一个基于随机截距的零通胀术语。

我正在比较该模型的所有子集以确定最适合的模型。

我的代码如下所示：

full <- glmmTMB(count ~ A + B + C + YN + offset(log(T)) + (1|ID), 
                zi = ~ (1|ID),
                family = poisson(),
                data = df, 
                na.action = "na.fail")
View(dredge(full, beta = sd))

但是当我查看表格时，模型似乎是重复的，即有偏移和没有偏移。该表包括 32 个模型，当我省略偏移量时，只有 16 个模型。加上系数非常相似。（显示此表的子集）

我已尝试更改我的代码，以便始终包含如下所示的偏移量，但这会导致警告。

    full <- glmmTMB(count ~ A + B + C + YN + (1|ID), 
                    offset = log(T), 
                    zi = ~ (1|ID),
                    family = poisson(),
                    data = df, 
                    na.action = "na.fail")
Warning messages:
1: In (function (start, objective, gradient = NULL, hessian = NULL,  :
  NA/NaN function evaluation
repeated 9 times

编辑 2020/06/16

当偏移量作为一个术语包含时，我已经用完整的模型摘要更新了这篇文章：

 Family: poisson  ( log )
Formula:          count ~ A + B + C + YN + offset(log(T)) +  
    (1 | ID)
Zero inflation:      ~(1 | ID)
Data: df

     AIC      BIC   logLik deviance df.resid 
   286.4    317.4   -135.2    270.4      344 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 ID   (Intercept) 0.1006   0.3172  
Number of obs: 352, groups:  name, 64

Zero-inflation model:
 Groups Name        Variance  Std.Dev. 
 ID   (Intercept) 2.089e-08 0.0001445
Number of obs: 352, groups:  name, 64

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -4.582134   0.382926 -11.966  < 2e-16 ***
A              0.003517   0.170809   0.021  0.98357    
B              0.356596   0.162703   2.192  0.02840 *  
C              0.594737   0.095198   6.247 4.17e-10 ***
YN             -1.397989   0.538510  -2.596  0.00943 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4744     0.7323  -0.648    0.517

当偏移量作为公式之外的单独参数包含时：

Family: poisson  ( log )
Formula:          count ~ A + B + C + YN + (1 | ID) #NOTE: offset isn't listed here despite being included when I specified the model above
Zero inflation:      ~(1 | ID)
Data: df
 Offset: log(duration)

     AIC      BIC   logLik deviance df.resid 
   286.4    317.4   -135.2    270.4      344 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 ID   (Intercept) 0.1006   0.3172  
Number of obs: 352, groups:  name, 64

Zero-inflation model:
 Groups Name        Variance  Std.Dev. 
 ID   (Intercept) 2.089e-08 0.0001445
Number of obs: 352, groups:  name, 64

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -4.582134   0.382926 -11.966  < 2e-16 ***
A              0.003517   0.170809   0.021  0.98357    
B              0.356596   0.162703   2.192  0.02840 *  
C              0.594737   0.095198   6.247 4.17e-10 ***
YN             -1.397989   0.538510  -2.596  0.00943 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4744     0.7323  -0.648    0.517

编辑 2020/06/24

修改后的疏浚功能解决了在所有模型中保留偏移变量的问题。

View(dredge(full, 
            beta = "sd"), #note sd was not in quotes above!
            fixed = "cond(offset(log(time)))") #retains offset in all models

但是，这导致了一个新问题：

Warning message:
In dredge(full, beta = "sd", fixed = "cond(offset(log(time)))") :
  do not know how to standardize coefficients of 'glmmTMB', argument 'beta' ignored

@Kamil Barton 在下面建议我需要使用罢工和 stdizeFit 来解决此警告消息。但到目前为止，如果没有其他错误消息，我无法让此代码正常工作。此外，我担心标准化我的响应变量（如下所示）是错误的方法，并且由于标准化时间变量会导致负值，因此无法再将其记录在模型中（产生 NA 值）。

无论如何，这是我尝试使用与帮助文件中几乎相同的方法来实施卡米尔建议的方法：

dfSTD <- stdize(select(df, counts, A, B, C, YN, T, ID),
                center = TRUE, scale = TRUE) 

full <- stdizeFit(glmmTMB(z.counts ~ z.A + z.B + z.C + z.YN + 
                              offset(log(z.T)) + (1|ID),
                          zi = ~ (1|ID),
                          family = poisson(),
                          data = dfSTD,
                          na.action = "na.fail"),
                   which = "formula", evaluate = FALSE)

Error in stdizeFit(glmmTMB(z.counts ~ z.A + z.B + z.C + z.YN + : 
  argument "data" is missing, with no default

任何关于解决这些最新错误的想法，无论是卡米尔建议的方法还是我上面所做的，都将不胜感激。

编辑 2020/07/10

我已经能够使用以下代码块解决上述问题：

dfSTD <- stdize(df, omit.cols = c("counts", "YN", "T", "ID"), 
                center = TRUE, scale = TRUE) 

full <- stdizeFit(glmmTMB(counts ~ z.A + z.B + z.C + YN + 
                                  offset(log(T)) + (1|ID),
                          zi = ~ (1|ID),
                          family = poisson(),
                          data = dfSTD,
                          na.action = "na.fail"),
                  data = dfSTD, which = "formula", evaluate = FALSE)

但是，现在我在使用疏通功能时，又出现了一个错误。

allEVCounts <- dredge(full, 
                      beta = "sd", 
                      fixed = "cond(offset(log(T)))")

Error in nobs.default(global.model) : no 'nobs' method is available

在确定如何以及在何处解决“nobs”方法的任何帮助将不胜感激。

【问题讨论】：

也许您可以提供标准的摘要输出？这可能比自编译的 Excel 表格更有帮助。
@Daniel 这就是您使用 View() 和 dredge() 时表格的样子。您是否正在寻找单个模型的摘要？该表显示了带有和不带有偏移量的模型的相似性，但实际上每个模型中都应该包含偏移量。
“原始”摘要可能有助于确定它是否与模型拟合或dredge() 函数有关。我不确定dregde() 在这里对 glmmTMB 模型做了什么。
@Daniel 我已经添加了上面的摘要。希望这可以为您提供更清晰的信息。我的目标是简单地比较这些模型的所有子集，但所有模型都必须包含偏移项。
可以给我们一个minimal reproducible example吗？

标签： r regression poisson mumin glmmtmb

【解决方案1】：

dredge 在公式中包含偏移项时，其处理方式与固定模型项相同。要在所有模型中保留偏移量，请添加参数fixed = "<offset term name>"，在本例中为offset(log(T))：

dredge(full, beta = "sd", fixed = "cond(offset(log(T)))")

请注意，在您的代码中，beta = sd（不带引号）被视为beta = "none"，因为您传递了一个函数 sd，而不是字符串"sd"，这不是一个可接受的值.

关于“NA/NaN 函数评估”的警告来自glmmTMB，某些子模型可能不会与偏移量收敛。设置options(warn = 1)（立即警告消息）并使用dredge 和trace = TRUE 找出哪些。

【讨论】：

Barton 感谢您的回答，感谢您了解我的 sd 术语。但是，当我在引号中包含 sd 时，会收到以下警告消息： In dredge(full, beta = "sd", fixed = "cond(offset(log(T)))") ：不知道如何标准化系数'glmmTMB'，参数 'beta' 被忽略。有解决这个问题的想法吗？
@tnt 需要实现函数std.coef 来标准化模型系数。默认不工作，glmmTMB 有三种系数。
我很困惑。在运行任何模型之前，我总是缩放系数：(x-mean)/sd。这不适用于 glmmTMB 吗？如何使用 std.coef？我试过 A2
@tnt 目前std.coef 不适用于glmmTMB 对象。但是，您可以标准化数据（例如使用stdize），然后在模型上使用stdizeFit。
这是我现在所处的位置：1) 标准化系数：df % mutate(Ac = stdize(A, center = TRUE, scale = TRUE), Bc = stdize(B , center = TRUE, scale = TRUE), Cc = stdize(C, center = TRUE, scale = TRUE))

【解决方案2】：

来自Kamil Bartoń 的解决方案是：在stdizeFit 中设置evaluate=TRUE，否则它会产生调用而不是拟合模型对象，这不是您（和dredge）所期望的。

【讨论】：

您的答案可以通过额外的支持信息得到改进。请edit 添加更多详细信息，例如引用或文档，以便其他人可以确认您的答案是正确的。你可以找到更多关于如何写好答案的信息in the help center。