tab_model() 与 glmmTMB 的奇怪输出

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，并且置信区间同样（错误地）取幂
(minor) 色散参数被标记为 (Intercept)，这令人困惑（从技术上讲它是色散模型的“截距”）

但是，R^2 值看起来很合理 — 我们稍后再谈。

模型本身呢？

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

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

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

theta_1|date.1 
     -11.77722

Large negative coefficients in zi (log-odds of zero-inflation), dispersion, or random effects (log-standard deviations) suggest unnecessary components (converging to zero on the constrained scale) ...

（如果您查看 summary(modop)，您会发现 date 随机效应的估计标准差为 7e-6，比 [=88= 小约 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: mu of 1.5 is too close to zero, estimate of random effect variances may be unreliable.
2: Model's distribution-specific variance is negative. Results are not reliable.

我基本上可以得出结论，这个数据集有点太大 small/model 太大而无法获得有用的 R^2 值。

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

tab_model() 与 glmmTMB 的奇怪输出

Weird output of tab_model() with glmmTMB

r

mixed-models

sjplot

glmmtmb