目录
1. 数据格式整理
（1）单组数据格式
a.cox回归数据提取
以R语言自带的数据集lung为例： 这是一份关于肺癌患者的生存数据。time是生存时间，以天为单位，status是生存状态，1代表删失，2代表死亡。
#加载系统自带的生存分析数据集
library(survival)
str(lung)
inst：机构代码
time：以天为单位的生存时间
status：删失状态1 = 删失，2 = 出现失效事件
age：岁
sex：性别，男= 1女= 2
ph.ecog：ECOG评分（0 =好，5 =死）
ph.karno：医师进行的Karnofsky评分（0 = 差，100 = 好）
pat.karno：患者自行进行的Karnofsky评分（0 = 差，100 = 好）
meal.cal：用餐时消耗的卡路里
wt.loss：最近六个月的体重减轻
使用cox回归探索危险因素。分类变量需要变为因子型，这样在进行回归时会自动进行哑变量设置。
lung$sex <- factor(lung$sex, labels = c("male","female"))
lung$ph.ecog <- factor(lung$ph.ecog, labels = c("asymptomatic", "symptomatic",'in bed <50%','in bed >50%'))
str(lung)拟合多因素Cox回归模型，这里我们只用sex/age/ph.ecog/ph.karno4个变量做演示：
（1）常规的cox分析如下：
fit.cox<-coxph(Surv(time,status)~sex+age+ph.karno+ph.ecog,data=lung)
#查看结果
summary(fit.cox)
(2) 现状需要对上述数据进行提取形成一个可以绘制森林图的表格(重点)
常见的格式：变量+HR+95%CI_low+95%CI_upper+P value
方法一：
思想：分别提取p value，HR，low CI，high CI
x <- summary(fit.cox) pvalue <- signif(as.matrix(x$coefficients)[,5],2)#把x$coefficients作为矩阵，提取第五列 HR <- signif(as.matrix(x$coefficients)[,2],2) low <- signif(x$conf.int[,3],2) high <- signif(x$conf.int[,4],2)x$coefficients；x$conf.int 分别是cox回归结果的第一个表和第二个表
将上述内容整合成为一个数据库格式
result_mul_cox <- data.frame(p.value=pvalue,HR=HR,low=low,high=high,stringsAsFactors = F) result_mul_cox
给数据框增加一列为HR（95%CI）
result_mul_cox$'HR(95%CI)' <- paste0(HR," (",low,"-",high,")",sep=" ") result_mul_cox
方法二（TableOne函数）：
第一步：提取变量HR+95%CI+95%CI
x <- summary(fit.cox) colnames(x$conf.int) multi1 <- as.data.frame(round(x$conf.int[,c(1,3,4)],2)) #提取出HR，95%CI，95%CI第二步：提取HR（95%CI）和pvalue
library(tableone) multi2 <- ShowRegTable(fit.cox,#fit.cox是你的cox回归模型 exp=T, #TRUE表示在表格中显示指数（HR） digits = 2, #小数点后保留的位数 pDigits = 3, #P值的小数点后显示的位数 printToggle = T, #TRUE表示打印结果表格 quote = F, #False表示不对列名添加引号 ciFun = confint) #表示使用confint函数计算置信区间
第三步：将两次提取的结果合并
result <- cbind(multi1,multi2) result
第四步：将行名转为第一列******
#将行名转为表格第一列 library(tibble) result <- rownames_to_column(result,var="Characteristics") result
b. COX回归数据的精修
（1）对第一列的数据名称进行修改
从表格中可以看到分类变量的名称是有点问题，由原变量的名称+分类变量名称组合而成，所以需要对原变量名称的字符串进行去除。
#修改变量名称
library(stringr)
result$Characteristics <- str_remove(result$Characteristics,"ph.ecog|sex")
result$Characteristics
（2）对数据的格式进行修改
#将第五行和第六行HR（95%CI）和P值改为字符型
str(result)
result[,5:6] <- lapply(result[,5:6],as.character)
str(result)#将2，3，4列转化为数值型，便于绘图
result[,2:4] <- lapply(result[,2:4],as.numeric)（3）创建新的result，补充空行，为了留给分类变量的名称
result <- rbind(
c("characteristic",NA,NA,NA,"HR (95%CI)","p"),
result[2,],
c("Sex",NA,NA,NA,NA,NA),
c("female",1.00,NA,NA,"1.00 (ref)",NA),
result[1,],
result[3,],
c("ph.ecog",NA,NA,NA,NA,NA),
c("asymptomatic",1.00,NA,NA,"1.00 (ref)",NA),
result[4,],
result[5,],
result[6,],
rep(NA,ncol(result)) #插入全为NA的空行
)
result
因为中间三列是数值型，所以不能给第一列标注名称！！！！！
C.总结每个变量的病例数
library(tableone)
str(lung)
myVars <- c("sex","age","ph.karno","ph.ecog") #所有需要统计的数据
catVars <- c("sex","ph.ecog") #所有的分类变量
table1 <- print(CreateTableOne(
vars=myVars,
data = lung,
factorVars = catVars),
showAllLevels=T) #设置showAlllevels参数可以显示所有的级别
插入基线表table1中的空行，使行数与变量result的行数一致
N <- rbind(c(NA,NA),
table1[4,],
c(NA,NA),
table1[2:3,],
table1[5,],
c(NA,NA),
table1[6:9,],
c(NA,NA));N
N <- N[,-1]
N <- as.data.frame(N)#注意对应相应的变量
head(N)
将result和N合并
result1 <- cbind(result[,1],N,result[,2:6])
#设置第一行名
result1[1,] <- c("characteristic","number (%)",NA,NA,NA,"HR (95%CI)","p") #因为中间三行为数值型不能命名
（2）多组数据格式
格式要求：
方法一：excel直接录入
方法二：使用上面的单组数据的方法，在列其他的HR
2.绘制单组森林图
（1）forestplot包
以前面的格式部分产生的result1为例，可以没有结果，可以把前面代码直接复制到R中跑一遍，生成result1
result1 library(forestplot) fig <- forestplot( result1[,c(1,2,6,7)],#需要在森林图中展示的列 mean=result1[,3], #均值列（HR） lower=result1[,4], #95%置信区间的下限数据列 upper=result1[,5], #95%置信区间的上限数据列 zero=1, #均值HR为1为参考线 boxsize=0.4, #方框的大小 graph.pos="right", #森林图在右侧 hrzl_lines=list( #水平线样式设置（三线表的线） "1"=gpar(lty=1,lwd=2), #均值线 "2"=gpar(lty=2), #下限和上限之间的虚线 "13"=gpar(lwd=2,lty=1,columns=c(1:4)) #下限和上限线 ), graphwidth=unit(.25,"npc"), #设置森林图宽度，整个图形的宽度 xlab="HR", #X轴标签 xticks=c(-0.8,1,3,6,10), #x轴刻度 is.summary=c(T,F,T,F,F,F,T,F,F,F,F), #判断是否为汇总行，不参与绘图的行设置为T，是否整行加粗 txt_gp=fpTxtGp( #设置文本样式 label = gpar(cex=0.8), #设置标签大小（正文内容字体大小） ticks = gpar(cex=1), #设置刻度大小 xlab = gpar(cex=0.9), #x轴标签大小 title = gpar(cex=1.2) #标题大小 ), lwd.zero=1, #HR=1时线条宽度，零刻度线得宽度 lwd.ci=2, #置信区间线得宽度 lwd.xaxis=2, #x轴得宽度 lty.ci=1.5, #置信区间得样式 ci.vertices=T, #是否显示置信区间得顶点 ci.vertices.height=0.2, #置信区间顶点的高度 clip=c(0.1,8), #X轴的截断范围 ineheight=unit(8,"mm"), #森林图内线的高度 line.margin=unit(8,"mm"), #森林图线的边距 colgap=unit(6,"mm"), #森林图中方框的间距 fn.ci_norm="fpDrawDiamondCI",#置信区间的形状 title="多因素COX回归森林图",#森林图的标题 col=fpColors( #设置颜色 box = "blue4", #方框的颜色 lines = "blue4", #线条颜色 zero = "black" #均值为0的颜色 ) ) fig
（2）forestploter包（重点）
a.格式要求
1.基本数据格式
重点：forestploter包格式不能包括第一行的变量名称
对result1进行修改，将第一行设置为列名
library(forestploter)
result2 <- result1
colnames(result2) <- result2[1,] #将第一行赋值为列名
#colnames(result2) <- c("characteristic","number","HR","low","high","HR(95%CI)","p" )
result2 <- result2[-1,]
result2 <- result2[-nrow(result2),]
标准格式如下：
前面可以通过加粗的方式修改分类变量的名称，这里可以通过分变量进行缩进体现
result2$characteristic <- ifelse(is.na(result2$number),
result2$characteristic,
paste0(" ",result2$characteristic))
2.把NA转为空字符串
result2$number <- ifelse(is.na(result2$number)," ",result2$number)
result2$`HR(95%CI)` <- ifelse(is.na(result2$`HR(95%CI)`)," ",result2$`HR(95%CI)`)
result2$p <- ifelse(is.na(result2$p)," ",result2$p)HR，low，high三行的NA不能随意修改为字符串，不然数据格式就被改变不是数值型，不能计算绘图。
3.计算标准误差SE，绘图时表现为正方形的大小
result2[,3:5] <- lapply(result2[,3:5],as.numeric)
str(result2)
result2$se <- -(log(result2$high)-log(result2$low))/1.964.为森林图添加空白列，为了产生一个绘图区间
result2$" " <- paste(rep(" ",12),collapse = " ")
b.设置主题
tm <- forest_theme(base_size = 10,
refline_col = "red4",
arrow_type = "closed",
footnote_col = "blue4")主题参数介绍：
所有主题参数都通过上述forest_theme( )函数调整
base_size=10 设置文本的基础大小
设置可信区间的外观：
ci_pch= 15 设置可信区间点的形状
ci_col="blue4" 设置可信区间的边框颜色
ci_fill="blue4" 设置可信区间的填充颜色
ci_alpha=0.8 设置可信区间的透明度
ci_lty=1, 设置可信区间的线条类型
ci_lwd=1.5 设置可信区间的线条宽度
ci_Theight=0.2 设置T字在可信区间末端的高度，默认NULL
设置x轴的粗细：
xaxis_lwd=0.6 设置x轴的粗细
xaxis_cex=1 设置x轴文本大小
设置参考线的外观：
refline_lwd=1 设置宽度
refline_lty="dashed" 设置线的类型
refline_col="grey20" 设置线的颜色
设置垂直线的外观：
vertline_lwd=1
vertline_lty="dashed"
vertline_col="grey20"
设置脚注的字体大小，字体样式和颜色：
footnote_cex=0.6 设置脚注字体大小
footnote_fontface="italic" 设置字体样式
footnote_col="red4" 设置脚注文本颜色
设置箭头格式：
arrow_type="closed" 设置箭头类型为闭合箭头
arrow_label_just="end" 设置箭头标签对齐方式为末尾
设置图例信息：
legend_name= "group" 设置图例标题
legend_value=c("trt1","trt2") 设置图例值
c.绘制森林图
p <- forest(result2[,c(1,2,9,6,7)], #选择要在森林图中显示的数据列，包括变量名称，样本数，绘图的空白列，HR（95%ci），p值
est = result2$HR,
lower = result2$low,
upper = result2$high,
sizes = result2$se,
ci_column = 3,
ref_line = 1,
arrow_lab = c("Low risk","High risk"),
xlim = c(-1,10),
ticks_at = c(-0.5,1,3,6,10),
theme = tm,
footnote = "XXXXXXX")
p
3.绘制多组森林图
以格式中提到的格式为例：
1.读取数据，设置缩进，设置绘图的空白区域
library(readxl)
library(forestploter)
slt <- read_excel("GRS森林图数据.xlsx",sheet=1)
str(slt)
slt$Var <- ifelse(is.na(slt$HR1),slt$Var,paste0(" ",slt$Var))
slt$HR1 <- ifelse(is.na(slt$HR1)," ",slt$HR1)
slt$HR2 <- ifelse(is.na(slt$HR2)," ",slt$HR2)
slt$HR3 <- ifelse(is.na(slt$HR3)," ",slt$HR3)
slt$plot1 <- paste0(rep(" ",20),collapse = " ")
slt$plot2 <- paste0(rep(" ",20),collapse = " ")
slt$plot3 <- paste0(rep(" ",20),collapse = " ")设置主题
tm <- forest_theme(
base_size = 10,
ci_pch = 18,
ci_col = "black",
ci_lty = 1,
ci_lwd =1.8,
ci_Theight = 0.2,
refline_lwd = 1,
refline_lty = "dashed",
refline_col = "grey20",
footnote_cex = 0.8,
footnote_fontface = "italic",
footnote_col = "red4"
)绘图函数，用list列表来描述每次需要描绘的结局：
p <- forest(slt[,c(1,2,17,7,18,12,19)],
est=list(slt$HR1_mean,
slt$HR2_mean,
slt$HR3_mean),
lower=list(slt$HR1_low,
slt$HR2_low,
slt$HR3_low),
upper=list(slt$HR1_high,
slt$HR2_high,
slt$HR3_high),
sizes = 0.6,
ci_column=c(3,5,7),
ref_line=c(1,1,1),
arrow_lab = c("Low risk","High risk"),
xlim = list(c(0.6,2),
c(0.6,2.2),
c(0.6,1.7)),
ticks_at = list(c(0.8,1,1.3,1.6),
c(0.8,1,1.3,1.6,2),
c(0.8,1,1.2,1.4)),
xlab = "HR",
theme=tm,
footnote = "
SSBs, sugar-sweetened beverages;
ASBs, artificial sweetened beverages;
PJs, pure fruit/vegetable juices;
HR, hazard ratio; CI, confidence interval."
)
p
对细节进行修饰，参考第五节部分内容细节调整
PP <- edit_plot(p,
part = "header",
gp=gpar(fontsize=12,
fontface="bold")) #标题栏放大
PP <- edit_plot(PP,
row=c(1,6,11),
gp=gpar(fontface="bold")) #变量名称加粗
PP <- add_border(PP,part = "header",where = "bottom") #增加边际线
PP <- add_border(PP,part = "header",where = "top")
PP <- edit_plot(PP,
row = c(5,10),
col=2,
gp=gpar(fontface="bold"))#统计学意义的HR加粗
PP <- edit_plot(PP,
row = c(5,10),
col=3,
which = "ci",
gp=gpar(col="green4",
lwd=2))
PP <- edit_plot(PP,
row = c(4,5,9,10),
col=4,
gp=gpar(fontface="bold"))
PP <- edit_plot(PP,
row = c(4,5,9,10),
col=5,
which = "ci",
gp=gpar(col="green4",
lwd=2))
PP <- edit_plot(PP,
row = c(3,8,13,15),
col=6,
gp=gpar(fontface="bold"))
PP <- edit_plot(PP,
row = c(3,8,13,15),
col=7,
which = "ci",
gp=gpar(col="green4",
lwd=2))
PP <- edit_plot(PP,
which = "background",
gp=gpar(fill="white"));PP#设置背景色
4.双变量森林图
5.常用的细节调整
以下内容只能修改forestplot包产生的森林图：
基本介绍：
使用edit_plot函数进行修改，基本修改格式为
edit_plot（p，
row=，
which=，
gp=gpar（））
p表示需要编辑的图片名称
row需要修改的行
which需要修改的内容：包括text文本，background（背景），ci(置信区间)三个内容
gpar（）中可以修改颜色col=，字体fontface=，填充颜色fill=
a.改变某行颜色
PP <- edit_plot(p,row=5,gp=gpar(col="red4",fontface="italic"))
PP
b.改变置信区间CI条的颜色
将7到10行的第三列置信区间改为绿色
PP <- edit_plot(PP,row=c(7:10),col=3,
which="ci",
gp=gpar(col="green4"));PP
c.改变文本字体，字号
PP <- edit_plot(PP,
row=c(2,6),
gp=gpar(fontface="bold"));PP
d.改变某行背景颜色
PP <- edit_plot(PP,
row=5,
which = "background",
gp=gpar(fill="yellow"));PP
e.顶部插入文本（森林图标题）
col设置标题在位置，在哪几列之间
PP <- insert_text(PP,
text="COX回归分析森林图",
col=1:5,
part="header",
gp=gpar(fontface="bold"));PP 
f.列名添加横线（三线表的横线）
PP <- add_border(PP,
part = "header",
where="bottom");PP
PP <- add_border(PP,
row = 10,
where="bottom");PP 
g.某行插入文本
just表示文本的位置，可以设置为center，left，right
PP <- insert_text(PP,
text="Smoking",
row=5,
just = "left",
gp=gpar(col="blue4",
cex=1,
fontface="italic"));PP
PP <- insert_text(PP,
text="80(58.56)",
row=5,
col=2,
just = "center",
gp=gpar(col="blue4",
cex=1,
fontface="italic"));PP
h.添加文本
PP <- add_text(p,"HR (95%CI)",col=2,part = "header",just = "left");PPI.修改列名****（重要）
因为列名重复造成的命名错误
plot <- edit_plot(
plot,
row = 1,
col = c(3,5), #修改第一行，第三列和第五列
part = "header",#列名，如果修改其他部分part设置为body
which="text",#设置为修改文本
label="Placebo\r\n(N=2841)"#填写修改的内容 \r\n表示换行
)
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
