R语言计算glmmTMB模型x-intercept时遇非数值参数错误

我用glmmTMB拟合了负二项模型，模型运行正常，但计算x-intercept时触发错误：

Error in b[1]/b[2] : non-numeric argument to binary operator

模型代码及摘要 #

mod <-  glmmTMB(total_count ~  mean_temp + (1|month), family = nbinom1, data = df)

summary(mod)

Family: nbinom1  ( log )
Formula:          total_count ~ mean_temp + (1 | month)
Data: df

     AIC      BIC   logLik deviance df.resid 
   251.2    258.9   -121.6    243.2       46 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 month  (Intercept) 3.415    1.848   
Number of obs: 50, groups:  month, 10

Dispersion parameter for nbinom1 family (): 43.8 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.58553    2.01287   0.291    0.771
mean_temp    0.06466    0.12167   0.531    0.595

计算x-intercept的代码（负二项模型需加入阈值C） #

C <- 0.01

b <- fixef(mod)

x_intercept <- log(C) - b[1]/b[2]

尝试将b转为numeric时，又出现：

Error: 'list' object cannot be coerced to type 'double'

复现数据 #

df <- structure(
  list(
    month = structure(
      c(
        1L,
        1L,
        1L,
        1L,
        1L,
        2L,
        2L,
        2L,
        2L,
        2L,
        3L,
        3L,
        3L,
        3L,
        3L,
        4L,
        4L,
        4L,
        4L,
        4L,
        5L,
        5L,
        5L,
        5L,
        5L,
        6L,
        6L,
        6L,
        6L,
        6L,
        7L,
        7L,
        7L,
        7L,
        7L,
        8L,
        8L,
        8L,
        8L,
        8L,
        9L,
        9L,
        9L,
        9L,
        9L,
        10L,
        10L,
        10L,
        10L,
        10L
      ),
      .Label = c(
        "April",
        "August",
        "December",
        "July",
        "June",
        "March",
        "May",
        "November",
        "October",
        "September"
      ),
      class = "factor"
    ),
    total_count = c(
      0L,
      0L,
      0L,
      0L,
      0L,
      29L,
      44L,
      10L,
      5L,
      2L,
      2L,
      43L,
      1L,
      0L,
      4L,
      10L,
      2L,
      4L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      3L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      1L,
      0L,
      0L,
      0L,
      0L,
      0L,
      0L,
      41L,
      7L,
      47L,
      1L,
      131L,
      81L,
      82L,
      409L,
      204L
    ),
    mean_temp = c(
      11.35,
      11.35,
      11.35,
      11.35,
      11.35,
      18.82,
      18.82,
      18.82,
      18.82,
      18.82,
      4.66,
      4.66,
      4.66,
      4.66,
      4.66,
      20.49,
      20.49,
      20.49,
      20.49,
      20.49,
      20.78,
      20.78,
      20.78,
      20.78,
      20.78,
      10.31,
      10.31,
      10.31,
      10.31,
      10.31,
      20.59,
      20.59,
      20.59,
      16.3,
      16.3,
      11.66,
      11.66,
      11.66,
      11.66,
      11.66,
      14.14,
      14.14,
      14.14,
      14.14,
      14.14,
      17.95,
      17.95,
      17.95,
      17.95,
      17.95
    )
  ),
  class = c("tbl_df",
            "tbl", "data.frame"),
  row.names = c(NA,-50L)
)

问题原因 #

fixef(mod)返回的是列表而非数值向量，因为glmmTMB会将固定效应按模型组件（条件模型、零膨胀模型等）分组存储。直接用b[1]/b[2]会尝试对列表元素做除法，自然触发非数值参数错误。

解决方法 #

需要提取条件模型的固定效应向量，有两种方式：

方式1：直接指定组件（推荐） #

b <- fixef(mod)$cond  # 提取条件模型的固定效应

方式2：转为向量（仅当模型只有条件模型时适用） #

b <- unlist(fixef(mod))  # 将列表转为数值向量

修正后的x-intercept计算代码（注意括号优先级，先计算log(C)-b[1]再除以b[2]）：

C <- 0.01

b <- fixef(mod)$cond
x_intercept <- (log(C) - b[1])/b[2]

验证结果 #

运行修正后的代码，得到x-intercept值：

> x_intercept
[1] -92.04755

内容的提问来源于stack exchange，提问作者Ahsk