绘图和 beta 回归输出不匹配

Question

在拟合一个似乎没有意义的 beta 回归模型后，我使用 tab_model() 得到了一个奇怪的输出。截距和系数似乎与图不对应。例如，soil_moisture 的 CI 似乎很高。另外，由于我必须回归线，所以不应该有两个截距吗？

我错过了什么？

这是我的数据：

df <- structure(list(weed_coverage = c(0.002, 0.01, 0.015, 0.01, 0.001, 
0.03, 0.006, 0.012, 0.03, 0.01, 0.002, 0.05, 0.004, 0.02, 0.02, 
0.006, 0.03, 0.01, 0.015, 0.01), soil_moisture = c(0.1125, 0.1343, 
0.1662, 0.3402, 0.2195, 0.1923, 0.2277, 0.2577, 0.148, 0.2715, 
0.104, 0.1495, 0.2788, 0.3477, 0.1835, 0.3175, 0.134, 0.3488, 
0.3477, 0.1097), distance = structure(c(2L, 2L, 1L, 2L, 1L, 1L, 
1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("2", 
"5"), class = "factor"), wc_pr = c(`1` = 0.0144421159096612, 
`2` = 0.0148173851570077, `3` = 0.0146600960637327, `4` = 0.0188698067335207, 
`5` = 0.0132256975788894, `6` = 0.0139395025623511, `7` = 0.0130176686719618, 
`8` = 0.0171297102952414, `9` = 0.0150581171360966, `10` = 0.0119600057973879, 
`11` = 0.0142983464947494, `12` = 0.0150847074475541, `13` = 0.0117921770637613, 
`14` = 0.019036340203918, `15` = 0.0141784723499839, `16` = 0.0109405819057696, 
`17` = 0.0148121562892363, `18` = 0.0190608859962305, `19` = 0.0103185336737258, 
`20` = 0.0163480105406738)), class = "data.frame", row.names = c(NA, 
-20L))

这是我的代码：

library(sjPlot)
library(betareg)
library(ggplot2)

betareg (weed_coverage ~ soil_moisture * distance, data = df) -> model_b # fit beta regression model

tab_model(model_b)

df %>% mutate(wc_pr= predict(model_b , type = "response")) -> df # create column with prediction values for weed_coverage

ggplot(df, aes(x = soil_moisture, y = weed_coverage, color = distance)) + # Plot the model
  geom_point(size = 2, shape = 21) +
  geom_line(aes(y = wc_pr, color = distance), data = df)+ 
  theme_bw()

tab_model(model_b)的输出：

这是我得到的情节：

Answer 1

好的，所以这里发生了很多事情，我将要经历这些事情。第一个是最简单的，就是澄清这个模型只有一个截距，也就是你的自变量为零时weed_coverage的估计值。因此，这些值并不总是可直接解释的。

图上的两条回归线是 soil_moisture 在两个不同距离值处的预测斜率。这是因为存在交互作用，所以soil_moisture 的直线斜率取决于distance 的水平。

现在，tab_model() 正确报告了 CI，但 sjPlot 包中的大多数函数的默认行为是使用 exp() 函数在对数尺度上转换值，以获得更接近概率值的值.如果您想要模型适合的日志空间中的原始模型值，您必须像这样指定transform=NULL：

tab_model(model_b,
          transform = NULL)

您现在会注意到模型报告的数字与您调用以下命令时相同：

summary(model_b)

Call:
betareg(formula = weed_coverage ~ soil_moisture * distance, data = df)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-2.2157 -0.5334 -0.0546  0.8254  1.7902 

Coefficients (mean model with logit link):
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -3.8822     0.7840  -4.952 7.35e-07 ***
soil_moisture            -1.9591     3.3122  -0.591    0.554    
distance5                -0.4752     0.9303  -0.511    0.610    
soil_moisture:distance5   3.1533     3.9372   0.801    0.423    

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)   
(phi)   102.30      34.72   2.946  0.00322 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 65.73 on 5 Df
Pseudo R-squared: 0.05632
Number of iterations: 399 (BFGS) + 6 (Fisher scoring) 

CI 非常大，但这是因为weed_coverage 的所有值都非常低。在 log-odds 空间中，小的百分比或比例变化在接近 0 或 1 时会变得非常大。在这个示例模型中，模型非常不确定，这就是概率尺度上的反向转换 CI 如此大的原因.还值得注意的是，这些估计值是针对weed_coverage 的变化，用于soil_moisture 的1 的变化。如果你缩放这个变量，模型不会改变，但数字会更容易在你的数据上下文中解释，因为不是估计回归线中的不确定性远远超出你的数据范围，一切都将报告一个soil_moisture 的 1 个标准差变化：

df$scaled_soil_moisture <- scale(df$soil_moisture)

model_b2 <- betareg(weed_coverage ~ scaled_soil_moisture * distance, data = df)

# Notice that I have let it back-transform the variable again.
tab_model(model_b2)

如果这里有什么不明白的地方请告诉我。

编辑：为了澄清您在评论中关于截距的问题，我用缩放自变量绘制了回归的两条回归线。当scaled_soil_moisture 等于 0 且距离等于 2 时，该模型中的截距对应于weed_coverage 的值（因为距离是一个二级因子，默认行为通常是创建一个虚拟变量，其中第一个由截距表示的水平）。可以看到，红线上x轴等于0的时候，y轴大约是0.01（已经向下取整）。这与非缩放模型形成对比，其中截距未在图表上显示，因为数据中没有零值。如果您将红线继续穿过数据到 x=0 的位置，您会看到它在 0.02 附近穿过 y 轴。

This page 解释说这是一种相当直观的方式。