你知道log-binomial GLM吗？OR≠RR
在队列研究中，常以相对危险度RR（包括危险度比和率比）来说明某因素与结局事件的关联强度；在病例对照中，我们常用OR值来说明某因素与结局事件的关联强度。虽然两者代表的意义相同，但计算公式不同（OR=ad/bc，RR=Ie/I0）。那么，两者之间又有着什么样的关系呢？由于病例对照研究无法得到事件的发生率，也就无法计算RR值，所以通过计算病例组与对照组的暴露比值之比—OR值—来近似估计RR值。既然是估计，那么估计的精度或准度是否受到一些因素的影响呢？有研究表明，当事件的发生率超过15%时（也有研究人员说是超过10%），用OR值去估计或代替RR值便不再精确，或者说偏差较大。
现实情况下，在队列研究中我们也会用二元logistic回归和OR值来说明某因素与结局事件的关联。然而，当发病率超过15%时，我们不建议仍使用logistic回归/模型，而应该改用log binomial 广义线性模型 （Generalized linear models, GLM），因为该模型可以得到我们想要的RR值，而不是OR值。在这里简单的解释一下GLM，其之所以叫广义，是因为他可以对多种分布类型的因变量进行建模/回归分析，主要是通过连接函数 link和分布函数family来完成，而我们常说的二元logistic模型在GLM中对应的连接函数和分布函数分别是logit和binomial（logit binomial GLM）。关于GLM的详细知识，请大家自行百度补充，这里就不多说了。接下来，我们通过公式计算和GLM拟合的方式来说明OR和RR的差异是怎么随着发病率变化的。
我们用的是来自SEER数据库的数据，看T分期为T1b的患者发生淋巴结转移的风险是T1a的多少倍。当结局事件发生率>15%时，我们以T1a为未暴露组，T1b为暴露组， 分别用公式和R语言的广义线性模型来计算T1b组患者发生结局事件的OR值/RR值。 OR值和RR值的计算公式如下：OR=ad/bc，RR=Ie/I0
****** 导入数据 ******
dat1<-read.csv("data1.csv",sep = ",")
dat2<-read.csv("data2.csv",sep = ",")
dat3<-read.csv("data3.csv",sep = ",")
****** 结局事件的频数和发生率 ******
table(dat1$LNM)
## 0 1
## 394 111
prop.table(table(dat1$LNM))
## 0 1
## 0.780198 0.219802
****** 利用公式计算OR值和RR值 ******
table(dat1$T_stage,dat1$LNM)
## 0 1
## T1a 162 14
## T1b 232 97
OR1=(97*162)/(232*14)
## [1] 4.838054
RR1=(97/(232+97))/(14/(14+162))
## [1] 3.70647
通过OR值和RR值的计算公式，我们得到了在事件发生率为21.98%时的OR值和RR值分别为4.838054和3.70647。接下来我们通过GLM来计算这两个值。
****** 利用广义线性模型计算OR值和RR值 ******
**** 单因素广义线性模型 ****
logit_T<-glm(LNM ~ T_stage,family = binomial("logit"),data = dat1)
exp(logit_T$coefficients)
## T_stageT1b
## 4.83805419
log_T<-glm(LNM ~ T_stage,family = binomial("log"),data = dat1)
exp(log_T$coefficients)
## T_stageT1b
## 3.70646982
从上面的结果我们可以看到，通过logit-binomial GLM和log-binomial GLM得到的OR值和RR值与手动计算的一致，而OR值（4.838）和RR值(3.706)相差1.132。同样，多因素模型给出的结果也有很大差异。
**** 多因素广义线性模型 ****
logit_fit<-glm(LNM ~ tumor_size + Grade + LN_number,family = binomial("logit"),data = dat1)
summary(logit_fit) ## logit-binomial GLM
## Call:
## glm(formula = LNM ~ tumor_size + Grade + LN_number, family = binomial("logit"), data = dat1)
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.1901 0.6310 -6.640 3.13e-11 ***
## tumor_size1- 1.3045 0.5165 2.526 0.01154 *
## tumor_size2- 1.5155 0.5198 2.915 0.00355 **
## tumor_size3- 1.7954 0.5528 3.248 0.00116 **
## tumor_size4- 2.1820 0.5121 4.261 2.04e-05 ***
## GradeⅡ 0.8727 0.4740 1.841 0.06562 .
## GradeⅢ 1.3398 0.4634 2.891 0.00384 **
## LN_number15- 0.7304 0.2361 3.093 0.00198 **
log_fit<-glm(LNM ~ tumor_size + Grade + LN_number,family = binomial("log"),data = dat1)
summary(log_fit) ## log-binomial GLM
## Call:
## glm(formula = LNM ~ tumor_size + Grade + LN_number, family = binomial("log"), data = dat1)
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9135 0.5602 -6.986 2.83e-12 ***
## tumor_size1- 1.1655 0.4687 2.487 0.012896 *
## tumor_size2- 1.3179 0.4682 2.815 0.004875 **
## tumor_size3- 1.4891 0.4818 3.091 0.001995 **
## tumor_size4- 1.7381 0.4552 3.818 0.000134 ***
## GradeⅡ 0.6994 0.4081 1.714 0.086534 .
## GradeⅢ 1.0494 0.3973 2.641 0.008255 **
## LN_number15- 0.5058 0.1748 2.894 0.003803 **
当结局事件的发生率<15%时，结果是怎样的呢？我们接下来仍以公式和拟合GLM的形式来分别计算OR值和RR值。在原数据的基础上，我们人为地调低了结局事件的频数，其他内容与前面的数据一致。
****** 利用公式计算OR值和RR值 ******
table(dat2$LNM) ## 结局事件的频数
## 0 1
## 435 70
prop.table(table(dat2$LNM)) ## 结局事件的发生率
## 0 1
## 0.8613861 0.1386139
table(dat2$T_stage,dat2$LNM) ## 暴露组和对照组结局事件的频数
## 0 1
## T1a 164 12
## T1b 271 58
prop.table(table(dat2$T_stage,dat2$LNM),1) ## 暴露组和对照组结局事件的发生率
## 0 1
## T1a 0.93181818 0.06818182
## T1b 0.82370821 0.17629179
OR2=(58*164)/(12*271)
## [1] 2.924969
RR2=(58/(58+271))/(12/(12+164)) = 0.17629179/0.06818182
## [1] 2.585613
****** 利用广义线性模型计算OR值和RR值 ******
logit_T<-glm(LNM ~ T_stage,family = binomial("logit"),data = dat2)
exp(logit_T$coefficients)
## T_stageT1b
## 2.92496925
log_T<-glm(LNM ~ T_stage,family = binomial("log"),data = dat2)
exp(log_T$coefficients)
## T_stageT1b
## 2.58561297
从上面的结果我们得知，OR=2.9250，RR=2.5856，差值为0.3394，差值变小了。接下来，我们接着调整事件的频数，看看OR值和RR值的差值是否会继续变小。
****** 利用公式计算OR值和RR值 ******
table(dat3$LNM) ## 结局事件的频数
## 0 1
## 458 47
prop.table(table(dat3$LNM)) ## 结局事件的发生率
## 0 1
## 0.90693069 0.09306931
table(dat3$T_stage,dat3$LNM) ## 暴露组和对照组结局事件的频数
## 0 1
## T1a 164 12
## T1b 294 35
prop.table(table(dat3$T_stage,dat3$LNM),1) ## 暴露组和对照组结局事件的发生率
## 0 1
## T1a 0.93181818 0.06818182
## T1b 0.89361702 0.10638298
OR3=(35*164)/(12*294)
## [1] 1.626984
RR3=(35/(35+294))/(12/(12+164))
## [1] 1.560284
****** 利用广义线性模型计算OR值和RR值 ******
logit_T<-glm(LNM ~ T_stage,family = binomial("logit"),data = dat3)
exp(logit_T$coefficients)
## T_stageT1b
## 1.62698413
log_T<-glm(LNM ~ T_stage,family = binomial("log"),data = dat3)
exp(log_T$coefficients)
## T_stageT1b
## 1.56028369
通过上面的结果我们可以看到，随着事件发生率的较低，OR值和RR值的差值（0.0667）在变小，也就是说当事件发生的率较低时用OR值近似估计RR值比较合适。那么，我们来粗糙地下个结论，当事件发生率较大时（>10%或15%），建议适用log-binomial GLM，以精确估计研究因素与结果的统计学关联。
最后留下一个问题。通过上面一系列的模拟，我们可以看到使用log binomial GLM得到的RR值与公式计算所得的RR值是一样的。那么，在队列研究中，是不是可以直接用log binomial GLM，无论结局事件发生率为多大？欢迎各位留言讨论。尊重数据，助力科研，我一直在路上！
数据分析与可视化最后编辑于 2022-06-08 · 浏览 2487
