如何在R语言中获取GLM模型二次项的顶点值及相关统计量？

计算代码 #

# 拟合你的二项式GLM模型（确保DF数据集已加载）
GLM <- glm(Y ~ X1 + I(X1^2) + X2 + X3 + X4 + X5, family="binomial", data=DF)

# 提取X1和X1二次项的回归系数
coefs <- coef(GLM)
b1 <- coefs["X1"]
b2 <- coefs["I(X1^2)"]

# 计算二次项的顶点值（即曲线极值对应的X1取值）
vertex_x <- -b1 / (2 * b2)

# 提取系数的方差-协方差矩阵
vcov_mat <- vcov(GLM)
var_b1 <- vcov_mat["X1", "X1"]
var_b2 <- vcov_mat["I(X1^2)", "I(X1^2)"]
cov_b1b2 <- vcov_mat["X1", "I(X1^2)"]

# 用德尔塔法计算顶点值的标准误
se_vertex <- sqrt( (1/(2*b2))^2 * var_b1 + (b1/(2*b2^2))^2 * var_b2 + 2*(1/(2*b2))*(b1/(2*b2^2))*cov_b1b2 )

# 计算95%置信区间（基于正态分布近似）
ci_low <- vertex_x - 1.96 * se_vertex
ci_high <- vertex_x + 1.96 * se_vertex

# 整理输出结果
result <- data.frame(
  顶点值 = as.numeric(vertex_x),
  标准误 = as.numeric(se_vertex),
  95%置信下限 = as.numeric(ci_low),
  95%置信上限 = as.numeric(ci_high)
)

print(result)

代码说明 #

• 顶点值计算：基于二次函数极值点公式 -X1系数/(2*X1二次项系数)，对应曲线的峰值或谷值位置。
• 标准误计算：采用德尔塔法（Delta Method），通过系数的方差-协方差矩阵推导顶点值的方差，开平方后得到标准误，这是处理函数形式统计量的常用方法。
• 置信区间：使用1.96倍标准误计算95%置信区间，和Stata默认的正态近似逻辑一致。

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