如何对loglinear模型执行power analysis（功效分析）

基础说明 #

你当前拟合的Admit * Dept * Gender为包含三向交互项的饱和对数线性泊松模型，这类模型的功效分析默认采用模拟法实现，操作步骤如下：

操作步骤 #

步骤1：准备依赖工具与基准模型 #

# 首次使用先安装依赖包
install.packages("simr")
# 加载所需包
library(simr)
library(MASS)

# 处理示例数据集
UCBAdf <- as.data.frame(UCBAdmissions)
# 拟合基准对数线性模型
glmFitT <- glm(Freq ~ Admit * Dept * Gender, family=poisson(link="log"), data=UCBAdf)

步骤2：确定待检验效应 #

首先查看模型系数，确认你需要做功效分析的目标效应对应的名称：

# 输出所有系数的命名
coef(glmFitT)

你可以选择检验单系数、某类交互项的整体显著性，比如三向交互Admit:Dept:Gender、Dept和Gender的两向交互Dept:Gender等。

步骤3：执行功效计算 #

使用模拟法计算指定效应的检验功效，默认显著水平为0.05：

# 示例：计算三向交互项的检验功效，nsim为模拟次数，数值越大结果越稳定
power_res <- powerSim(
  fit = glmFitT,
  test = fixed("Admit:Dept:Gender", method = "fdr"), # 指定检验的效应和多重检验校正方法
  nsim = 1000
)
# 输出功效结果
print(power_res)

输出结果中的power值即为当前样本量、当前效应量下，能检测到目标效应显著的概率，行业通用标准为功效≥0.8时检验效能达标。

步骤4：生成不同样本量下的功效曲线 #

如果需要评估样本量变化对功效的影响，可以生成功效曲线：

# breaks参数对应原样本量的放大倍数，此处设置为原样本的1/2/3/4倍
power_curve_res <- powerCurve(
  fit = glmFitT,
  test = fixed("Admit:Dept:Gender"),
  along = "Freq",
  nsim = 200,
  breaks = c(1,2,3,4)
)
# 绘制功效曲线
plot(power_curve_res)

注意事项 #

• 如果要检测自定义的最小有意义效应，可以通过fixef(glmFitT)["目标效应名"] <- 预期效应值修改基准模型的系数后再运行功效分析
• 正式分析时建议将nsim参数设置为1000及以上，降低随机模拟带来的结果波动
• 同时检验多个效应时，需在test参数中指定对应的多重检验校正方法，避免一类错误膨胀

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