1 单因素COX和多因素COX

clinical.txt

library(survival)
rt=read.table("clinical.txt",header=T,sep="\t",check.names=F,row.names=1)
#单因素cox
uniTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
  cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
  coxSummary = summary(cox)
  uniTab=rbind(uniTab,
               cbind(id=i,
                     HR=coxSummary$conf.int[,"exp(coef)"],
                     HR.95L=coxSummary$conf.int[,"lower .95"],
                     HR.95H=coxSummary$conf.int[,"upper .95"],
                     pvalue=coxSummary$coefficients[,"Pr(>|z|)"])
  )
}
write.table(uniTab,file="uniCox.txt",sep="\t",row.names=F,quote=F)



#多因素cox
multiCox=coxph(Surv(futime, fustat) ~ ., data = rt)
multiCoxSum=summary(multiCox)
multiTab=data.frame()
multiTab=cbind(
  HR=multiCoxSum$conf.int[,"exp(coef)"],
  HR.95L=multiCoxSum$conf.int[,"lower .95"],
  HR.95H=multiCoxSum$conf.int[,"upper .95"],
  pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"])
multiTab=cbind(id=row.names(multiTab),multiTab)
write.table(multiTab,file="multiCox.txt",sep="\t",row.names=F,quote=F)

uniCox.txt

2 绘制森林图

R-forestplot包| HR结果绘制森林图
R语言 | forestplot包绘制森林图

2.1 载入数据

#载入R包
library(forestplot)
#就从这个网站下载数据
#数据来源：https://www.r-bloggers.com/forest-plot-with-horizontal-bands/
data <- read.csv("ForestPlotData.csv", stringsAsFactors=FALSE)
#查看数据
head(data)

ForestPlotData.csv

2.2 简单森林图绘制

## 构建tabletext，更改列名称，展示更多信息
np <- ifelse(!is.na(data$Count), paste(data$Count," (",data$Percent,")",sep=""), NA)

## The rest of the columns in the table.
tabletext <- cbind(c("Subgroup","\n",data$Variable),
                   c("No. of Patients (%)","\n",np),
                   c("4-Yr Cum. Event Rate\n PCI","\n",data$PCI.Group),
                   c("4-Yr Cum. Event Rate\n Medical Therapy","\n",data$Medical.Therapy.Group),
                   c("P Value","\n",data$P.Value))



##绘制森林图
forestplot(labeltext=tabletext, graph.pos=3,
           mean=c(NA,NA,data$Point.Estimate),
           lower=c(NA,NA,data$Low), upper=c(NA,NA,data$High),
           boxsize=0.5)

2.3 优化森林图

## 定义亚组，方便后面线条区分
subgps <- c(4,5,8,9,12,13,16,17,20,21,24,25,28,29,32,33)
data$Variable[subgps] <- paste("  ",data$Variable[subgps])

forestplot(labeltext=tabletext,
           graph.pos=3, #为Pvalue箱线图所在的位置
           mean=c(NA,NA,data$Point.Estimate),
           lower=c(NA,NA,data$Low), upper=c(NA,NA,data$High),
           #定义标题
           title="Hazard Ratio Plot",
           ##定义x轴
           xlab="    <---PCI Better---   ---Medical Therapy Better--->",
           ##根据亚组的位置，设置线型，宽度造成“区块感”
            hrzl_lines=list("3" = gpar(lwd=1, col="#99999922"),
                            "7" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
                            "15" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
                            "23" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
                            "31" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922")),
           #fpTxtGp函数中的cex参数设置各个组件的大小
           txt_gp=fpTxtGp(label=gpar(cex=1.25),
                          ticks=gpar(cex=1.1),
                          xlab=gpar(cex = 1.2),
                          title=gpar(cex = 1.2)),
           ##fpColors函数设置颜色
           col=fpColors(box="#1c61b6", lines="#1c61b6", zero = "gray50"),
           #箱线图中基准线的位置
           zero=1,
           cex=0.9, lineheight = "auto",
           colgap=unit(8,"mm"),
           #箱子大小，线的宽度
           lwd.ci=2, boxsize=0.5,
           #箱线图两端添加小竖线，高度
           ci.vertices=TRUE, ci.vertices.height = 0.4)

image.png