如何在R中基于GEE计算非劣效性试验的对应p值？

核心思路 #

GEE默认仅提供针对β=0的双侧检验，而你的非劣效检验属于单侧假设检验，需手动计算检验统计量和p值。由于GEE针对聚类数据的特性，**必须使用稳健标准误（Robust S.E.）**进行推断，才能得到可靠结果。

具体步骤与R代码实现 #

1. 提取模型关键结果：从GEE输出中取出genderF的系数估计值和稳健标准误
2. 计算检验统计量：根据你设定的假设H₀:β₁≤-0.01 vs Hₐ:β₁>-0.01，检验统计量公式为：

   z = (β̂₁ - θ₀) / SE_robust
   

   其中θ₀=-0.01是你指定的非劣效界值（对应logOR），SE_robust为稳健标准误。
3. 计算单侧p值：备择假设为β₁ > -0.01，因此p值是标准正态分布中Z≥z的概率，即1 - pnorm(z)。

完整可运行代码

library(gee)
data(muscatine)

data2 <- muscatine
data2$obese <- ifelse(muscatine$obese == 'yes', 0, 1)

# 拟合GEE模型
x <- gee(obese ~ gender + age + occasion, id = id, family = binomial, corstr = "exchangeable", data = data2)

# 提取genderF的系数和稳健标准误
beta_hat <- coef(x)["genderF"]
se_robust <- summary(x)$coefficients["genderF", "Robust S.E."]

# 非劣效界值（logOR）
theta0 <- -0.01

# 计算检验统计量
z_stat <- (beta_hat - theta0) / se_robust

# 计算单侧p值
p_value <- 1 - pnorm(z_stat)

# 输出结果
cat("genderF系数估计值:", beta_hat, "\n")
cat("稳健标准误:", se_robust, "\n")
cat("检验统计量z:", z_stat, "\n")
cat("非劣效检验p值:", p_value, "\n")

结果解释 #

代入你提供的模型输出数据：

• β̂₁=-0.1521，SE_robust=0.0626，θ₀=-0.01
• 计算得z_stat ≈ (-0.1521 + 0.01)/0.0626 ≈ -2.27
• p_value ≈ 1 - pnorm(-2.27) ≈ 0.988

该p值远大于常用的0.05显著性水平，说明没有足够统计学证据拒绝原假设H₀，即无法证明男性受试者的肥胖水平非劣于女性。

重要注意事项 #

• 务必使用稳健标准误：Naive S.E.未考虑聚类数据的相关性，会导致推断偏差，不符合GEE的设计逻辑。
• 验证非劣效界值合理性：你直接使用logOR界值-0.01，需确认该界值是基于临床意义预先设定的（通常需从绝对风险差转换为logOR，而非直接取0.01）。

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