rm(list=ls())
gc() #清列表清內(nèi)存
#單因素和多因素COX回歸
setwd("C:/AW/COX")
library(survival)
data=read.csv("UCSC.csv",header=T) #載入需要做COX的數(shù)據(jù)，包括存活狀態(tài)，時間，和對存活有影響的變量
unimodel=coxph(Surv(OSTIME, OS) ~ AGE, data)? #對AGE分組做單因素COX
print(unimodel)? # 單因素結(jié)果，exp(coef)即HR, 保存exp(coef)和p值
exp(confint(unimodel)) #得到HR的95%置信區(qū)間，保存
summary(coxmodel) #詳細(xì)版本的單因素結(jié)果
coxmodel=coxph(Surv(OSTIME, OS) ~ AGE? + T + STAGE + M + N , data) #多因素分析
print(coxmodel)? #多因素結(jié)果
exp(confint(coxmodel)) #多因素置信區(qū)間
#做森林圖
cox=read.csv("COX.csv",header=T) #森林圖需要HR值，和HR的95%置信區(qū)間上下限
cox$group = ifelse(cox$lower >1|cox$upper <1,"red3", "black")? #對有意義的變量標(biāo)紅
library(ggplot2)
ggplot(data=cox)+
? aes(x=hr,y=reorder(X,hr))+ #根據(jù)HR值對縱坐標(biāo)的變量排序
? geom_errorbarh(aes(xmax=upper,xmin=lower),color=cox$group,height=0,size=0.8)+ #畫出HR值上下限
? geom_point(size=1.5,shape=15,color=cox$group)+ #畫出HR值
? geom_vline(xintercept=1,linetype="dashed",size=0.6,color="darkblue")+ #在1處標(biāo)藍(lán)色虛線
? coord_trans(x="log2")+ #對X軸取log
? scale_x_continuous(limits=c(0.2,10),breaks=c(0.5,1,2,3,5,10))+ #規(guī)定X軸上下限和需要顯示的坐標(biāo)
? labs(x="Hazard ratio",y="")+
? theme_bw()+
? theme(panel.grid.minor = element_blank()) #去除多余雜線
