R语言处理非治疗因素在两组变量不均衡的方法（结果为二分类变量）
2020-10-25 本文已影响0人
灵活胖子的进步之路
library(survival)
## Warning: package 'survival' was built under R version 3.6.3
library(compareGroups)
## Warning: package 'compareGroups' was built under R version 3.6.3
data(colon)#加载R内置数据集colon
str(colon)#显示coloon数据集结构
## 'data.frame': 1858 obs. of 16 variables:
## $ id : num 1 1 2 2 3 3 4 4 5 5 ...
## $ study : num 1 1 1 1 1 1 1 1 1 1 ...
## $ rx : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
## $ sex : num 1 1 1 1 0 0 0 0 1 1 ...
## $ age : num 43 43 63 63 71 71 66 66 69 69 ...
## $ obstruct: num 0 0 0 0 0 0 1 1 0 0 ...
## $ perfor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ adhere : num 0 0 0 0 1 1 0 0 0 0 ...
## $ nodes : num 5 5 1 1 7 7 6 6 22 22 ...
## $ status : num 1 1 0 0 1 1 1 1 1 1 ...
## $ differ : num 2 2 2 2 2 2 2 2 2 2 ...
## $ extent : num 3 3 3 3 2 2 3 3 3 3 ...
## $ surg : num 0 0 0 0 0 0 1 1 1 1 ...
## $ node4 : num 1 1 0 0 1 1 1 1 1 1 ...
## $ time : num 1521 968 3087 3087 963 ...
## $ etype : num 2 1 2 1 2 1 2 1 2 1 ...
newdata<-na.omit(colon)#删除缺失值
str(newdata)
## 'data.frame': 1776 obs. of 16 variables:
## $ id : num 1 1 2 2 3 3 4 4 5 5 ...
## $ study : num 1 1 1 1 1 1 1 1 1 1 ...
## $ rx : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
## $ sex : num 1 1 1 1 0 0 0 0 1 1 ...
## $ age : num 43 43 63 63 71 71 66 66 69 69 ...
## $ obstruct: num 0 0 0 0 0 0 1 1 0 0 ...
## $ perfor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ adhere : num 0 0 0 0 1 1 0 0 0 0 ...
## $ nodes : num 5 5 1 1 7 7 6 6 22 22 ...
## $ status : num 1 1 0 0 1 1 1 1 1 1 ...
## $ differ : num 2 2 2 2 2 2 2 2 2 2 ...
## $ extent : num 3 3 3 3 2 2 3 3 3 3 ...
## $ surg : num 0 0 0 0 0 0 1 1 1 1 ...
## $ node4 : num 1 1 0 0 1 1 1 1 1 1 ...
## $ time : num 1521 968 3087 3087 963 ...
## $ etype : num 2 1 2 1 2 1 2 1 2 1 ...
## - attr(*, "na.action")= 'omit' Named int 127 128 165 166 179 180 187 188 197 198 ...
## ..- attr(*, "names")= chr "127" "128" "165" "166" ...
descrTable(rx~ ., data = newdata)#对比不同术后化疗方案患者的基线信息情况
## Warning in cor(x, y): 标准差为零
##
## --------Summary descriptives table by 'rx'---------
##
## ______________________________________________________
## Obs Lev Lev+5FU p.overall
## N=610 N=588 N=578
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
## id 459 (269) 480 (271) 461 (268) 0.354
## study 1.00 (0.00) 1.00 (0.00) 1.00 (0.00) 0.355
## sex 0.53 (0.50) 0.56 (0.50) 0.47 (0.50) 0.007
## age 59.5 (12.0) 60.1 (11.7) 59.9 (12.0) 0.653
## obstruct 0.20 (0.40) 0.20 (0.40) 0.18 (0.38) 0.473
## perfor 0.03 (0.17) 0.03 (0.18) 0.03 (0.16) 0.810
## adhere 0.15 (0.35) 0.15 (0.36) 0.13 (0.34) 0.553
## nodes 3.83 (3.75) 3.71 (3.60) 3.44 (3.23) 0.143
## status 0.55 (0.50) 0.53 (0.50) 0.40 (0.49) <0.001
## differ 2.08 (0.50) 2.02 (0.52) 2.09 (0.51) 0.033
## extent 2.89 (0.50) 2.89 (0.43) 2.87 (0.50) 0.508
## surg 0.30 (0.46) 0.26 (0.44) 0.25 (0.43) 0.166
## node4 0.28 (0.45) 0.28 (0.45) 0.24 (0.43) 0.331
## time 1439 (931) 1479 (965) 1716 (922) <0.001
## etype 1.50 (0.50) 1.50 (0.50) 1.50 (0.50) .
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
newdata$treat<-ifelse(newdata$rx=="Obs","obs","chemo")#根据化疗方案分为观察组及化疗组
newdata$status<-factor(newdata$status,labels=c("alive","dead"))#生存状态因子化变为二分变量并赋值
newdata$treat<-as.factor(newdata$treat)#治疗方式因子化
descrTable(treat~ .-rx,show.ratio =TRUE,data = newdata)#不同治疗方式患者的基线信息差异
## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## Variables 'study' have been removed since some errors occurred
##
## --------Summary descriptives table by 'treat'---------
##
## ____________________________________________________________________
## chemo obs OR p.ratio p.overall
## N=1166 N=610
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
## id 471 (269) 459 (269) 1.00 [1.00;1.00] 0.381 0.381
## sex 0.51 (0.50) 0.53 (0.50) 1.06 [0.87;1.29] 0.548 0.548
## age 60.0 (11.8) 59.5 (12.0) 1.00 [0.99;1.00] 0.377 0.380
## obstruct 0.19 (0.39) 0.20 (0.40) 1.11 [0.87;1.42] 0.408 0.413
## perfor 0.03 (0.17) 0.03 (0.17) 0.95 [0.54;1.70] 0.873 0.873
## adhere 0.14 (0.35) 0.15 (0.35) 1.04 [0.79;1.38] 0.768 0.769
## nodes 3.57 (3.42) 3.83 (3.75) 1.02 [0.99;1.05] 0.145 0.156
## status: 0.001
## alive 626 (53.7%) 274 (44.9%) Ref. Ref.
## dead 540 (46.3%) 336 (55.1%) 1.42 [1.17;1.73] <0.001
## differ 2.05 (0.51) 2.08 (0.50) 1.12 [0.93;1.36] 0.232 0.229
## extent 2.88 (0.46) 2.89 (0.50) 1.05 [0.86;1.29] 0.619 0.629
## surg 0.25 (0.44) 0.30 (0.46) 1.23 [0.99;1.53] 0.063 0.067
## node4 0.26 (0.44) 0.28 (0.45) 1.09 [0.87;1.36] 0.457 0.460
## time 1597 (951) 1439 (931) 1.00 [1.00;1.00] 0.001 0.001
## etype 1.50 (0.50) 1.50 (0.50) 1.00 [0.82;1.22] 1.000 1.000
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
可以看到两者在surg有差异，为后续分析演示，同时纳入nodes两个变量，考虑这两个变量在组件有差异
descrTable(status~ .-rx,show.ratio =TRUE,data = newdata)
## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## Variables 'study' have been removed since some errors occurred
##
## --------Summary descriptives table by 'status'---------
##
## ____________________________________________________________________
## alive dead OR p.ratio p.overall
## N=900 N=876
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
## id 470 (264) 463 (275) 1.00 [1.00;1.00] 0.603 0.603
## sex 0.53 (0.50) 0.51 (0.50) 0.93 [0.77;1.12] 0.460 0.460
## age 60.0 (11.6) 59.6 (12.3) 1.00 [0.99;1.00] 0.416 0.416
## obstruct 0.17 (0.38) 0.21 (0.41) 1.27 [1.00;1.61] 0.050 0.050
## perfor 0.02 (0.15) 0.04 (0.19) 1.51 [0.87;2.63] 0.141 0.139
## adhere 0.12 (0.32) 0.17 (0.38) 1.58 [1.21;2.06] 0.001 0.001
## nodes 2.75 (2.46) 4.60 (4.18) 1.21 [1.16;1.25] <0.001 <0.001
## differ 2.03 (0.49) 2.09 (0.53) 1.27 [1.06;1.53] 0.010 0.010
## extent 2.81 (0.53) 2.96 (0.41) 2.02 [1.63;2.51] <0.001 <0.001
## surg 0.24 (0.43) 0.30 (0.46) 1.38 [1.12;1.71] 0.003 0.003
## node4 0.15 (0.36) 0.38 (0.49) 3.37 [2.69;4.23] <0.001 <0.001
## time 2323 (447) 741 (586) 1.00 [1.00;1.00] <0.001 0.000
## etype 1.51 (0.50) 1.49 (0.50) 0.93 [0.77;1.12] 0.448 0.448
## treat: 0.001
## chemo 626 (69.6%) 540 (61.6%) Ref. Ref.
## obs 274 (30.4%) 336 (38.4%) 1.42 [1.17;1.73] <0.001
## ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ
#以status为结局变量，观察个自变量对结果的影响，可以看到，nodes和surg都是由影响的。所以不同治疗措施的生存结果差异可能是两个因素影响的了
psModel<-glm(status~adhere+differ,family=binomial(link="logit"),data=newdata)
#以上述挑选的两个自变量为因变量，结局变量为因变量建立二元回归方程
newdata$ps=predict(psModel,type="response")
#计算倾向性评分并在数据集内添加一列为PS的列，内容为评分
newdata$IPTW<-ifelse(newdata$status==1,1/newdata$ps,1/(1-newdata$ps))
#计算各个患者的加权值（IPTW）
fit<-glm(status~treat,family=binomial(link="logit"),data=newdata)
summary(fit)
##
## Call:
## glm(formula = status ~ treat, family = binomial(link = "logit"),
## data = newdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.265 -1.115 -1.115 1.241 1.241
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.14778 0.05873 -2.516 0.011861 *
## treatobs 0.35176 0.10037 3.505 0.000457 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2461.7 on 1775 degrees of freedom
## Residual deviance: 2449.4 on 1774 degrees of freedom
## AIC: 2453.4
##
## Number of Fisher Scoring iterations: 3
#在不矫正的情况下观察治疗方式对于生存的影响，P值为0.
fit.IPTW<-glm(status~treat,weights = ps,family=binomial(link="logit"),data=newdata)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
summary(fit.IPTW)
##
## Call:
## glm(formula = status ~ treat, family = binomial(link = "logit"),
## data = newdata, weights = ps)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0116 -0.7744 -0.7291 0.8484 0.9801
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.12652 0.08363 -1.513 0.1303
## treatobs 0.34432 0.14285 2.410 0.0159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1214.4 on 1775 degrees of freedom
## Residual deviance: 1208.5 on 1774 degrees of freedom
## AIC: 679.98
##
## Number of Fisher Scoring iterations: 3
#可以看到，加权后的P值为0.0159，虽然也有意义，但是在校准后其变小了，说明加权后其效果变小了
