运行广义估计方程后如何计算效应大小和置信区间？

Question

我是一名博士生，没有统计学背景。

我目前正在进行一项随机对照试验，我需要使用 GEE 模型进行统计分析。

我有两组（干预组和对照组），我每个月（7个月内）收集整个赛季受伤的受试者（受伤与否，即受伤发生率），最后我想比较整个赛季的两组。我读过两篇使用 GEE 的论文，其研究设计与我的研究相似（Harøy 等人，2017 年：https://doi.org/10.1136/bjsports-2017-098937 和 Åkerlund 等人，2022 年：https： //doi.org/10.1007/s00167-021-06644-2），虽然他们使用SPSS进行统计分析，但也不知道如何在那里进行GEE。

我开始阅读“gee”包中的文档并最终运行此代码：

migee<-gee(Injury ~ GRUPO * Tiempo, family = poisson(link = "log"), data = GEE_trial, corstr = "exchangeable", id = id)

“Injury”是一个二分变量，表示球员是否受伤 “GRUPO”是群体，要么是控制，要么是干预 “Tiempo”是每个月 GEE_Trial 是数据框 “id”是分配给每个主题的编号

我使用了对数泊松和可交换相关结构，因为我在上述论文中读到他们使用了这种格式。虽然我不确定这是否正确执行。

运行分析和“摘要”功能后，出现了各种系数（请注意，我只使用三个月而不是整个赛季的数据进行分析）以及标准误差和 z 分数。 系数：

                 Estimate Naive S.E.    Naive z Robust S.E.   Robust z
(Intercept)    -2.1400662  0.6492832 -3.2960444   0.6642112 -3.2219666
GRUPO1          1.1415373  0.7362179  1.5505427   0.7289733  1.5659522
Tiempo2         0.6931472  0.7952062  0.8716571   0.7952062  0.8716571
Tiempo3         0.6931472  0.7952062  0.8716571   0.7952062  0.8716571
GRUPO1:Tiempo2 -1.9459101  1.0836834 -1.7956445   1.0816408 -1.7990354
GRUPO1:Tiempo3 -1.5404450  1.0167826 -1.5150191   1.0016388 -1.5379247          

我读到，通过对“GRUPO1”系数求幂，我可以获得优势比，因为这是我需要在论文中呈现的效应大小。然而，在上述论文中，他们提出了优势比（或比率比 (Åkerlund 2022) 和置信区间。

首先，我不知道我是否正确地进行了分析。 其次，如果正确，我如何计算该效应大小的置信区间（优势比或比率）？ 第三，如果不正确，我该如何正确地进行分析？

感谢您阅读到这里，如果有必要，我愿意使用其他软件，我想要的只是完成我的博士学位并进行正确的分析。

Answer 1

我们无法帮助解答“这是否正确”，因为这是基于当前的科学问题，只有您和主题专家才能确定。如果您对统计学及其在分析中的作用很幼稚，我建议您咨询您正在攻读博士学位的机构内的统计学家来指导您。

就客观编码问题而言，我不相信

gee

对象可以使用典型函数来收集置信区间（即

confint()

。

但是，您可以直接使用模型的结果计算 Wald 置信区间。

您没有提供数据，因此使用此可重现的数据：

set.seed(123)
n <- 1e3
GEE_trial <- data.frame(Injury = sample(0:1, n, replace = TRUE),
                 id = factor(LETTERS[1:5]),
                 GRUPO =  factor(sample(0:1, n, replace = TRUE)),
                 Tiempo = factor(sample(1:3,  n, replace = TRUE)))

您可以建模 (

midge

)、提取结果 (

mdlresults

)、计算置信区间并将结果组织在数据框中 (

finalres

)：

# model
migee <- gee::gee(Injury ~ GRUPO * Tiempo,
                  family = poisson(link = "log"),
                  data = GEE_trial,
                  corstr = "exchangeable",
                  id = id)

# Wald confidence interval
ci <- 0.95
xx <- qnorm((1 + ci) / 2)

# Wald confidence interval
ci <- 0.95
xx <- qnorm((1 + ci) / 2)

# final results
finalres <- data.frame(param = rownames(mdlresults),
                       est = mdlresults$Estimate,
                       lwr =  mdlresults$Estimate - xx * mdlresults$Robust.S.E.,
                       upr = mdlresults$Estimate + xx * mdlresults$Robust.S.E.)

输出

# > finalres
#            param         est         lwr         upr
# 1    (Intercept) -0.79094492 -0.96054071 -0.62134914
# 2         GRUPO1  0.14851284 -0.07233009  0.36935577
# 3        Tiempo2  0.07898841 -0.15237451  0.31035133
# 4        Tiempo3  0.12337459 -0.09494036  0.34168953
# 5 GRUPO1:Tiempo2 -0.24411386 -0.56936113  0.08113342
# 6 GRUPO1:Tiempo3 -0.12646163 -0.42464761  0.17172435