前面我們講過了單因素和多因素cox回歸分析，那么怎么樣將結果以森林圖的形式來展示呢？
森林圖簡介
森林圖（forest plot），從定義上講，它一般是在平面直角坐標系中，以一條垂直于X軸的無效線（通常坐標X=1或0）為中心，用若干條平行于X軸的線段，來表示每個研究的效應量大小及其95%可信區(qū)間，并用一個棱形來表示多個研究合并的效應量及可信區(qū)間，它是Meta分析中最常用的結果綜合表達形式，現在也廣泛應用在biomarker此類研究中。
森林圖的科研用途
提到森林圖，很多人的第一反應就是Meta分析。實際上，除了Meta分析，森林圖還有很多用處。森林圖可以直觀的反映出效應量（例如RR、OR、HR或者WMD）大小及其95% CI，這些效應量指標通常都是通過采用多因素回歸分析所得，因此我們同樣可以把森林圖借鑒過來，用于展示單因素或者多因素回歸分析的結果。總結來說，森林圖的科研用途主要用于Meta和臨床實驗。
臨床實驗普通分析，常規(guī)森林圖
下圖就是常規(guī)Cox回歸結果的森林圖展示，主要體現了變量、病人數量、P值和HR值。比如: ph.ecog變量位于無效線（即中間的那條豎線）右側，說明ph.ecog有助于死亡。森林圖在常規(guī)情況下事件結局是"生/死"這種兩分類，但有時候事件結局是"有效/無效"、"治療/未治療"等等其他二分類情況，評估事件是好事還是壞事。比如生存（生：0；死：1），位于無效線左側的變量，說明這些變量不利于事件發(fā)生，是保護因素；位于無效線右側的變量，說明這些變量有助于事件發(fā)生，是危險因素；當與無效線相交時，說明這些變量與事件發(fā)生之間關系不強！在整體數據上，用來評估這些變量因素對事件結局的影響！
小編在下面這篇文章中
看到了如下的森林圖，
今天小編就帶大家一起來重現這張圖，我們還是用單因素和多因素cox回歸分析中提到的lung這套數據來舉例。小編用三種不同的方法來實現這張圖。
第一種，我們用最原始的plot函數，lines函數從底層來實現。后邊兩種方法，我們用現成的R包來實現。
#加載這兩個R包
library("survival")
library("survminer")
?
#加載肺癌這套數據
data("lung")
?
###########################################
#批量單因素cox回歸分析
############################################
#假設我們要對如下5個特征做單因素cox回歸分析
covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")
#分別對每一個變量，構建生存分析的公式
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~', x)))
#對每一個特征做cox回歸分析
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
?
#提取HR，95%置信區(qū)間和p值
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
#獲取p值
p.value<-signif(x$wald["pvalue"], digits=2)
#獲取HR
HR <-signif(x$coef[2], digits=2);
#獲取95%置信區(qū)間
HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res<-c(p.value,HR)
names(res)<-c("p.value","HR (95% CI for HR)")
return(res)
})
#轉換成數據框，并轉置
res <- t(as.data.frame(univ_results, check.names = FALSE))
res <-as.data.frame(res,stringsAsFactors=F)
?
#############################################################
#對HR (95% CI for HR)做處理，得到HR和low .95和high .95
#當然也可以改計算univ_results這一步的代碼，不要將HR和CI貼起來
############################################################
HR=gsub("[\\(\\)]","",res$`HR (95% CI for HR)`)
HR=gsub("-"," ",HR)
HR=as.data.frame(do.call(cbind,strsplit(HR," ")),stringsAsFactors=F)
names(HR)=rownames(res)
?
#################################
#開始繪圖，直接保存到pdf文件中
#################################
pdf(file="univariate_forest.pdf",width=7)
#左邊和右邊邊距稍微留多一點來寫變量名稱，pvalue和HR
par(mar=c(5,6,4,13))
#先用小方塊畫出HR
plot(as.numeric(HR[1,]),1:dim(HR)[2],
pch=15,cex=2,col="blue",bty='n',yaxt='n',ylab=NA,xlab="Hazard Ratio",
xlim=range(as.numeric(unlist(HR)))
)
#添加中線
abline(v=1,col="grey",lwd=2,lty=2)
?
for(i in 1:ncol(HR)){
x=as.numeric(HR[2:3,i])
#循環(huán)畫出CI
lines(x,c(i,i),col="blue")
#添加變量名
text(0.2,i,rownames(res)[i],xpd=T,adj = c(0,0))
#添加p值
text(2.1,i,as.numeric(res[i,1]),xpd=T,adj = c(0,0))
#添加HR和CI
text(2.7,i,as.character(res[i,2]),xpd=T,adj = c(0,0))
}
#添加標題
text(2.1,ncol(HR)+0.5,"pvalue",xpd=T,adj = c(0,0))
text(2.7,ncol(HR)+0.5,"HR(CI)",xpd=T,adj = c(0,0))
dev.off()
會得到下面這張圖，是不是跟文章中的長的很像，這可是小編純手工打造的。
多因素的森林圖，如果你理解了上面單因素的繪圖的思路和原理，應該也不難，大家可以自己練練手。
第二種方法，使用survivalAnalysis包來實現
這個包不僅可以畫forest圖，還可以計算cox回歸的結果。
先來看單因素cox分析的結果和forest圖
再來看看多因素cox分析的結果和forest圖
第三種方法，使用ggforest函數來實現
完整代碼參考
