为什么 R 的 glm 的协方差矩阵与 minitab 的概率分析的协方差矩阵不同

Question

使用下面 R 脚本中的数据集。 R 二项式概率 glm vcov 函数的输出与 Minitab 的概率分析方差协方差表的输出不同。 vcov 的 R 文档说它返回方差协方差表。 拟合的模型系数相同（至 3 dp），但协方差矩阵完全不同。如何从 R 获得与从 Minitab 获得相同的矩阵输出，反之亦然？

我原以为 R 和 Minitab 中的表格是相同的。

我询问了 Minitab 支持人员，我已向他们支付了年度订阅费用，但得到了这个毫无帮助的答复：

“感谢您的询问。我们无法评论其他应用程序中使用的方法。但是，您会在下面的页面中找到使用的公式。

https://support.minitab.com/en-us/minitab/21/help-and-how-to/statistical-modeling/reliability/how-to/probit-analysis/methods-and-formulas/equation/

它指向包含方程 πj = c + (1 − c) g(β0 + xjβ) 的页面，没有其他内容 - 我期待有关奇异值分解或类似内容的信息。

我只订阅 Minitab 以获得支持，但现在它似乎没有什么帮助，所以我可能只使用 R。

来自R

df<-data.frame(stimulus = c(0.615,0.634,0.655,0.675),
   success = c(0,3,3,2),
   failure = c(5,4,3,0))
fm <- glm(cbind(success,failure)~stimulus,data=df,
    family=binomial("probit"))
coef(fm)

(Intercept)    stimulus 
  -29.68637    45.89187 

vcov(fm)
            (Intercept)  stimulus
(Intercept)    156.9548 -244.3060
stimulus      -244.3060  380.5229 

来自 Minitab

相同的系数但不同的协方差矩阵

回归表

可变系数标准误差 Z P
常数 -29.6864 13.1351 -2.26 0.024
刺激 45.8920 20.4461 2.24 0.025

自然
回复 0

矩阵VCCO1

172.531 -268.482
-268.482 418.044

Answer 1

“带着两块手表的人永远不知道现在是几点。”

就其价值而言，

glmmTMB

— 使用与

glm

完全不同的拟合过程 — 给出了接近 Minitab 的协方差矩阵。

m_glm <- glm(cbind(success,failure)~stimulus,data=df,
      family=binomial("probit"))
library(glmmTMB)
m_tmb <- glmmTMB(cbind(success,failure)~stimulus,data=df,
      family=binomial("probit"))
vcov(m_tmb)
Conditional model:
            (Intercept)  stimulus
(Intercept)    172.5314 -268.4821
stimulus      -268.4821  418.0434

有更多关于如何计算协方差矩阵的信息

glm

here ...

也就是说，您还应该小心使用标准误差 (

sqrt(diag(vcov(...)))

) 作为 GLM 参数不确定性的度量，因为它们仅对 large 数据集可靠。

glm

和

glmmTMB

的剖面置信区间非常相似。

> confint.default(m_glm)
                2.5 %    97.5 %
(Intercept) -54.24111 -5.131622
stimulus      7.65886 84.124881
> confint(m_glm)
Waiting for profiling to be done...
                2.5 %   97.5 %
(Intercept) -59.70676 -6.66262
stimulus     10.01528 92.56096

> confint(m_tmb, method = "wald", estimate = FALSE) ## equiv. confint.default()
                 2.5 %    97.5 %
(Intercept) -55.430789 -3.942062
stimulus      5.818323 85.965601

> confint(m_tmb, method = "profile", estimate = FALSE) ## equiv. confint()
                2.5 %    97.5 %
(Intercept) -59.70642 -6.660542
stimulus     10.01219 92.560959