带有 glmmTMB 的 tab_model() 的奇怪输出

Question

当我使用 sjPlot 包的 tab_model() 函数和 glmmTMB 包的 glmmTMB 函数来拟合具有 beta 系列响应的广义线性混合模型时，我得到了一个奇怪的输出.截距和边际 R² 看起来很奇怪。

这是怎么回事？

df <- structure(list(date = structure(c(6L, 5L, 6L, 1L, 4L, 2L, 2L, 
2L, 2L, 4L, 6L, 1L, 6L, 6L, 2L, 2L, 4L, 4L, 5L, 1L), .Label = c("2021-03-17", 
"2021-04-07", "2021-04-13", "2021-04-27", "2021-05-11", "2021-05-27"
), class = "factor"), kettlehole = structure(c(4L, 6L, 6L, 4L, 
7L, 2L, 6L, 5L, 3L, 5L, 1L, 1L, 1L, 1L, 4L, 4L, 5L, 4L, 3L, 5L
), .Label = c("1189", "119", "1202", "149", "172", "2484", "552"
), class = "factor"), plot = structure(c(8L, 4L, 4L, 3L, 7L, 
8L, 1L, 3L, 6L, 4L, 4L, 3L, 6L, 1L, 2L, 7L, 5L, 8L, 1L, 1L), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8"), class = "factor"), treatment = structure(c(2L, 
2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
1L, 2L, 1L), .Label = c("a", "b"), class = "factor"), distance = structure(c(2L, 
2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 
2L, 1L, 1L), .Label = c("2", "5"), class = "factor"), soil_moisture_content = c(0.2173, 
0.1028, 0.148, 0.3852, 0.1535, 0.2618, 0.2295, 0.222, 0.3145, 
0.1482, 0.2442, 0.3225, 0.1715, 0.1598, 0.2358, 0.274, 0.1543, 
0.144, 0.128, 0.361), yield = c(0.518, 0.434, 0.35, 0.599, 0.594, 
0.73, 0.568, 0.442, 0.695, 0.73, 0.667, 0.49, 0.744, 0.56, 0.485, 
0.532, 0.668, 0.511, 0.555, 0.718), weed_coverage = c(0, 0.045, 
0.03, 0.002, 0.11, 0.003, 0.01, 0, 0.02, 0.002, 0, 0.008, 0, 
0.002, 0, 0.006, 0, 0, 0.02, 0.002)), row.names = c(NA, -20L), class = c("tbl_df", 
"tbl", "data.frame"))

library(sjPlot)
library(glmmTMB)

glmmTMB(yield ~ soil_moisture_content + weed_coverage + distance + treatment + (1/kettlehole/plot) + (1|date), family = "beta_family", data = df) -> modop

tab_model(modop)

编辑

所以这是我在 n=630 的实际数据集上使用的 tab_model() 结果的屏幕截图。我认为问题在于模型过拟合，正如 Ben 所提到的，需要通过消除不必要的预测变量来进行调整。

Answer 1

tl;dr 奇怪的拦截结果似乎是sjPlot::tab_model 中的一个错误，应该在sjPlot issues list 向维护人员报告——似乎tab_model 错误地将不应该的分散参数。但是，您的模型还有其他问题（可能是过度拟合），这些问题会破坏您的边际 R^2 值。

这里是一些合理的模拟数据，显示了tab_model() 的问题：

set.seed(101)
## rbeta() function parameterized by mean and shape
my_rbeta <- function(n, mu, shape0) {
  rbeta(n, shape1 = mu*shape0, shape2 = (1-mu)*shape0)
}
n <- 100; ng <- 10
dd <- data.frame(x = rnorm(n),
                 f = factor(rep(1:(n/ng), ng)))
dd <- transform(dd,
                y = my_rbeta(n,
                             mu = plogis(-1 + 2*x + rnorm(ng)[f]),
                             shape0 = 5))

m1 <- glmmTMB(y ~ x + (1|f), family = "beta_family", dd)
tab_model(m1)

sigma(m1)、print(m1)、summary(m1)的结果均同意估计的色散参数为5.56（接近其标称值5），同意confint(m1, "disp_")：

         2.5 %   97.5 % Estimate
sigma 4.068351 7.606602 5.562942

但是，tab_model() 报告：

显示两个问题：

• （主要）离散度报告为 exp(5.563) = 260.6 而不是 5.563，并且置信区间类似地（不正确地）取幂
• （次要）色散参数标记为(Intercept)，令人困惑（从技术上讲，它是色散模型的“截距”）

但是，R^2 值看起来很合理——我们会回到这个。

---

模型本身呢？

一个合理的经验法则（例如，请参阅 Harrell 回归建模策略）说，您通常应该针对每 10-20 个观察值设置大约 1 个参数。在固定效应和随机效应之间，您有 9 个参数（length(modop$fit$par) 或 nobs(modop) - df.residual(modop)）用于 20 次观察。

如果我们运行diagnose(modop)（注意我使用的是diagnose() 的固定/开发版本，您的结果可能会略有不同）给出：

diagnose(modop)
Unusually large coefficients (|x|>10):

theta_1|date.1 
     -11.77722 

zi 中的大负系数（零通胀的对数几率）， 分散或随机效应（对数标准差）建议 不必要的组件（在约束尺度上收敛到零）...

（如果您查看summary(modop)，您会看到date 随机效应的估计标准偏差为7e-6，比次大随机效应项小约4 个数量级......）

modop2 <- update(modop, . ~ . - (1|date))

diagnose(modop2) 说这个型号没问题。

然而，tab_model(modop2) 仍然给出了一个可疑的条件 R^2（1.038，即 >1）。直接运行performance::r2_nakagawa(modop2)（我相信这是tab_model()使用的底层机制）给出：

# R2 for Mixed Models
  Conditional R2: 1.038
     Marginal R2: 0.183

但有警告

1：1.5 的 mu 太接近于零，随机效应方差的估计可能不可靠。
2：模型的特定分布方差为负。结果不可靠。

我基本上会得出结论，这个数据集有点太小/模型太大而无法获得有用的 R^2 值。

FWIW 我有点担心tab_model() 报告此型号的N_plot = 8：它应该像summary(modop2) 一样报告N_{plot:kettlehole} = 18