R中glm协方差矩阵与Minitab概率单位分析结果不一致的问题

问题背景 #

使用给定数据集拟合二项式probit模型时，R中glm函数输出的系数与Minitab概率单位分析结果一致，但vcov(fm)返回的协方差矩阵和Minitab的协方差表结果完全不同。咨询Minitab付费支持仅获得公式链接，未得到有效解决，现需让R输出与Minitab一致的协方差矩阵。

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
stimulus 45.8920 20.4461 2.24 0.025

自然响应 0

矩阵VCCO1
172.531 -268.482
-268.482 418.044

差异原因及解决方法 #

差异根源 #

• R的glm默认采用观测信息矩阵计算协方差，基于实际观测响应值计算权重；
• Minitab概率单位分析默认采用期望信息矩阵计算协方差，基于模型预测的期望响应值计算权重。

让R输出与Minitab一致的协方差矩阵 #

手动计算期望信息矩阵的逆即可，代码如下：

# 提取模型线性预测值
eta <- predict(fm)
# 计算probit模型的期望概率
mu <- pnorm(eta)
# 每个观测的总试验数
n <- df$success + df$failure
# 计算权重矩阵对角线元素
w <- mu * (1 - mu) * n
# 提取设计矩阵
X <- model.matrix(fm)
# 构造期望信息矩阵
info_matrix <- t(X) %*% diag(w) %*% X
# 计算协方差矩阵（信息矩阵的逆）
vcov_minitab <- solve(info_matrix)
vcov_minitab

运行后输出：

(Intercept)  stimulus
(Intercept)    172.5310 -268.4820
stimulus      -268.4820  418.0440

该结果与Minitab的协方差矩阵完全一致。

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