cox回归亚组分析中计算交互效应的三种方法

1.原理

   交互作用：某一因素的真实效应（单独效应）随着另一因素水平的改变而改变。当两种或两种以上暴露因素同时存在时所致的效应不等于它们单个作用相联合的效应时，则称因素之间存在交互作用。

① 因素A的效应在因素B的不同水平上存在差异，则认为因素A、B之间存在交互作用。

② 因素A、B的联合效应不等于两因素独立效应之和或之积。

    统计建模中一般线性模型交互项反映的是因素间相加交互作用（additive interactions，INTA），而logistic和Cox等广义线性模型则反映因素间统计学上相乘交互作用（multiplicative interaction，INTM）

    相乘交互作用的定义：假设研究多风险因素中交互作用的两暴露因素为AB，则OR00表示AB均无暴露，即OR00＝1；OR10表示A暴露、B无暴露，OR01表示A无暴露、B暴露，OR11表示A、B均暴露。则相乘交互作用INTM＝ORA×B＝OR11/（OR10×OR01）

2.示例

下载所需要的包

install.packages('dplyr')
install.packages('survival')
install.packages('Publish')
install.packages('jstable')

用survival包的lung数据做示例

library(dplyr)
library(survival)

survival::lung
#    inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
#1      3  306      2  74   1       1       90       100     1175      NA
#2      3  455      2  68   1       0       90        90     1225      15
#3      3 1010      1  56   1       0       90        90       NA      15
#4      5  210      2  57   1       1       90        60     1150      11
#5      1  883      2  60   1       0      100        90       NA       0
#6     12 1022      1  74   1       1       50        80      513       0
#7      7  310      2  68   2       2       70        60      384      10
#8     11  361      2  71   2       2       60        80      538       1
#9      1  218      2  53   1       1       70        80      825      16
#10     7  166      2  61   1       2       70        70      271      34
#…………
#inst       病人检测机构代码
#time       生存时间（天）
#status     删失状态 1=删失，2=死亡
#ph.ecog    医生评定的ECOG体力状态评分，因子数字越小健康状况越好
#ph.karno   医生评定的karnofsky功能状态评分，越高健康状况越好
#pat.karno  病人评定的karnofsky功能状态评分
#meal.cal   饮食能量摄入
#wt.loss    最近六个月体重减轻

设置变量类型

#将sex，ph.age设为因子，将age四分后设为因子
data_int <- survival::lung
data_int$sex <- as.factor(data_int$sex)
data_int$ph.ecog <- as.factor(data_int$ph.ecog)
qq <- quantile(data_int$age)
data_int$AGE <- 1
data_int<-within(data_int,{
  AGE[age<qq[2]]<-1
  AGE[age>=qq[2]&age<qq[3]]<-2
  AGE[age>=qq[3]&age<qq[4]]<-3
  AGE[age>=qq[4]]<-4 })
data_int$AGE <- as.factor(data_int$AGE)
data_int <- data_int %>% dplyr::select(time,status,AGE,sex,ph.ecog) %>% na.omit()

1.常规方法anova

#以status~ph.ecog的sex亚组和AGE亚组为例，3种方法都用似然比检验
#1 常规方法anova
fit1 <- coxph(Surv(time,status)~ph.ecog+AGE+sex,data_int)
fit2 <- coxph(Surv(time,status)~ph.ecog*sex+AGE,data_int)
anova(fit1,fit2,test="Chisq")
#Analysis of Deviance Table
# Cox model: response is  Surv(time, status)
# Model 1: ~ ph.ecog + AGE + sex
# Model 2: ~ ph.ecog * sex + AGE
#   loglik  Chisq Df Pr(>|Chi|)
#1 -728.46                     
#2 -728.35 0.2196  2      0.896
fit3 <- coxph(Surv(time,status)~ph.ecog*AGE+sex,data_int)
anova(fit1,fit3,test="Chisq")
#Analysis of Deviance Table
# Cox model: response is  Surv(time, status)
# Model 1: ~ ph.ecog + AGE + sex
# Model 2: ~ ph.ecog * AGE + sex
#   loglik  Chisq Df Pr(>|Chi|)  
#1 -728.46                       
#2 -720.98 14.966  6    0.02052 *
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#优势：好理解，准确
#劣势：代码行数多，结果需整理

2.Publish::subgroupAnalysis()

#2 Publish包
library(Publish)
sub_cox <-subgroupAnalysis(fit1,data_int,treatment="ph.ecog",subgroups=c('sex','AGE'))
sub_cox
# subgroups  level sample_0 sample_1 sample_2 sample_3 event_0 event_1 event_2 event_3 #time_0 time_1 time_2
#      <char> <fctr>    <int>    <int>    <int>    <int>   <num>   <num>   <num>   <num>  
#<num>  <num>  <num>
#1:       sex      1       36       71       29        1      64     125      57       2  #13184  19414   6299
#2:       sex      2       27       42       21       NA      36      70      37      NA   #8984  16118   5405
#3:       AGE      1       14       30        5       NA      19      52      10      NA   #4471   9738   1221
#4:       AGE      2       21       22       13       NA      32      38      26      NA   #8499   7176   2125
#5:       AGE      3       13       33        9       NA      22      55      17      NA   #4823   8907   3662
#6:       AGE      4       15       28       23        1      27      50      41       2   #4375   9711   4696
#   time_3 HazardRatio pinteraction subgroup  Lower confint
#    <num>       <num>       <char>   <char> <char>  <char>
#1:    118   3.1014850       0.8960     <NA>     NA   NA-NA
#2:     NA          NA       0.8960     <NA>     NA   NA-NA
#3:     NA          NA       0.0205     <NA>     NA   NA-NA
#4:     NA   0.8538482       0.0205     <NA>     NA   NA-NA
#5:     NA   2.1016129       0.0205     <NA>     NA   NA-NA
#6:    118          NA       0.0205     <NA>     NA   NA-NA
#优势：可一行代码跑所有协变量的交互作用

3.jstable::TableSubgroupMultiCox()

#3 jstable包
library(jstable)
cox_sub_plot <- TableSubgroupMultiCox(formula = Surv(time,status) ~ ph.ecog,
                                      var_subgroups = c("AGE","sex"),
                                      var_cov = c("AGE","sex"),
                                      data=data_int)
cox_sub_plot
#        Variable Count Percent    Levels Point Estimate Lower Upper   KM P value P for interaction
#ph.ecog  Overall   227     100 ph.ecog=0      Reference             92.3                      <NA>
#1                              ph.ecog=1           1.53  1.03  2.28 93.9   0.037              <NA>
#2                              ph.ecog=2           2.66  1.69  4.19  100  <0.001              <NA>
#3                              ph.ecog=3           6.88   0.9 52.44  100   0.063              <NA>
#4            AGE  <NA>    <NA>      <NA>           <NA>  <NA>  <NA> <NA>    <NA>             0.021
#5              1    49    21.6 ph.ecog=0      Reference              100                      <NA>
#6                              ph.ecog=1           1.87  0.69  5.04 90.3   0.215              <NA>
#7                              ph.ecog=2           3.11  0.86 11.26  100   0.085              <NA>
#8                              ph.ecog=3           <NA>  <NA>  <NA> <NA>    <NA>              <NA>
#9              2    56    24.7 ph.ecog=0      Reference             87.8                      <NA>
#10                             ph.ecog=1           2.09  0.91  4.81  100   0.082              <NA>
#11                             ph.ecog=2           9.52  3.78 23.98  100  <0.001              <NA>
#12                             ph.ecog=3           <NA>  <NA>  <NA> <NA>    <NA>              <NA>
#13             3    55    24.2 ph.ecog=0      Reference             80.8                      <NA>
#14                             ph.ecog=1           1.56   0.7  3.47 92.9   0.275              <NA>
#15                             ph.ecog=2           1.11  0.43   2.9  100   0.828              <NA>
#16                             ph.ecog=3           <NA>  <NA>  <NA> <NA>    <NA>              <NA>
#17             4    67    29.5 ph.ecog=0      Reference              100                      <NA>
#18                             ph.ecog=1           0.98  0.48     2 95.2   0.948              <NA>
#19                             ph.ecog=2           2.78  1.23  6.32 87.9   0.014              <NA>
#20                             ph.ecog=3           4.96   0.6 40.97  100   0.137              <NA>
#21           sex  <NA>    <NA>      <NA>           <NA>  <NA>  <NA> <NA>    <NA>             0.896
#22             1   137    60.4 ph.ecog=0      Reference             94.3                      <NA>
#23                             ph.ecog=1           1.38  0.86  2.22 94.8   0.184              <NA>
#24                             ph.ecog=2           2.32  1.35  3.99  100   0.002              <NA>
#25                             ph.ecog=3           5.33  0.69 41.27  100   0.109              <NA>
#26             2    90    39.6 ph.ecog=0      Reference             74.2                      <NA>
#27                             ph.ecog=1           1.62  0.76  3.48 93.6   0.215              <NA>
#28                             ph.ecog=2           4.57  1.82 11.51  100   0.001              <NA>
#29                             ph.ecog=3           <NA>  <NA>  <NA> <NA>    <NA>              <NA>
#优势：可一行代码同时实现亚组分析和交互作用分析

总的来说，用jstable::TableSubgroupMultiCox()显然是在做分析想要做出森林图时更快速方便的方式；