锹形虫
2023/02/05阅读：22主题：嫩青
一些帮助练习的轮子
来看一些用于练习的“轮子”
“分析过程中常需要出一些图表报告，但已经被定义的函数并不总能达到要求，我们试试“重复造轮子”来加以练习吧~
1. 输出些制表数据看看？
举个常用的 cox 回归和逻辑回归的例子，建模过程中总需要不断调整模型，然后重新获取 summary 结果，抑或是一次性写很多可能的模型都报告一下用于后续讨论。在后续制表的时候，则需要各个模型的参数。如果能一次将所有模型的参数都按要求做成数据框就好了，这个功能自己攒攒就能写，但我也提供了一个例子。用于提示在 R 中的很多结果展示都是以数据对象的形式出现的，对于 summary 结果，我们可以用 coef(), 当然也可以用 $coefficient 来获得我们想要的数据，在 R 中只要格式满足，就可以调用。
这里准备了逻辑回归和 cox 回归的模型参数获取代码，额外提示是，我们可以将所有获取的数据直接做 rbind() 合并然后定义成数据框，即可直接复制粘贴到 Excel 中开始制表了。请参考下述代码：
# 逻辑回归分析参数输出
log_info <- function(fit)
{
a <- as.data.frame(coef(summary(fit)))
P.value <- format.pval(a$`Pr(>|z|)`[c(2)], eps = 0.001)
OR <- round(exp(a$Estimate[c(2)]),2)
b <- as.data.frame(confint(fit))[c(2),]
CI95 <- paste('(',round(exp(b$`2.5 %`),2),', ',round(exp(b$`97.5 %`),2), ')',sep = '')
dat <- data.frame(P.value = c(P.value), OR = c(OR), CI95 = c(CI95))
return(dat)
}
# Cox 回归分析参数输出
cox_info <- function(fit)
{
a <- as.data.frame(summary(fit)$coefficients)
P.value <- format.pval(a$`Pr(>|z|)`[1], eps = 0.001)
HR <- round(a$`exp(coef)`,2)[1]
b <- as.data.frame(summary(fit)$conf.int)[1,]
CI95 <- paste('(',round((b$`lower .95`),2),', ',round((b$`upper .95`),2), ')',sep = '')
dat <- data.frame(P.value = c(P.value), HR = c(HR), CI95 = c(CI95))
return(dat)
}
2. 画个森林图康康？
森林图在分析中也十分常见，很多模型下想要关注或者对比不同模型的主要研究因素的效应大小，森林图是一种不错的展示形式。下面是一个该图的例子：
虽然网络上也提供了很多函数可以直接输出这样的图，比如 forestplot() 函数，但真正使用的时候往往又遇到这样那样的问题。比如，想让误差线别那么长？想让标目加加粗？想换个字体？这些已定义函数可能就鞭长莫及了。但还好森林图称之为森林图，最起码的它是一张图。因此，在 ggplot2 中它就可以很好地被绘制。这里我提供了一个多 cox 回归模型的森林图绘制方法，它只关注了主要研究因素的效应。在这里你会看到 （1. 输出些制表数据看看？）中提到的参数输出函数被再次定义和使用。
额外的提示是， ggplot 的坐标概念，面对一张白纸，我们如果想要放东西上去，那首先得确定放在哪儿。这就是坐标，只要正确的映射坐标，再联合其形状大小颜色等其他参数，我们就能绘制绝大多数的 2D 图形。一张森林图自然不在话下。请参考下述代码：
# 我提供了 18 个模型的情况用于参考，其输出结果即上图
forest_fit_info <- data.frame( rbind(
forest_info(fit1),
forest_info(fit2),
forest_info(fit3),
forest_info(fit4),
forest_info(fit5),
forest_info(fit6),
forest_info(fit7),
forest_info(fit8),
forest_info(fit9),
forest_info(fit10),
forest_info(fit11),
forest_info(fit12),
forest_info(fit13),
forest_info(fit14),
forest_info(fit15),
forest_info(fit16),
forest_info(fit17),
forest_info(fit18)))
forest_fit_info$labels <- c('Model 1',
'Model 2',
'Model 3',
'Model 4',
'Model 5',
'Model 6',
'Model 7',
'Model 8',
'Model 9',
'Model 10',
'Model 11',
'Model 12',
'Model 13',
'Model 14',
'Model 15',
'Model 16',
'Model 17',
'Model 18')
# 搞搞图形映射的数据
forest_fit_info$HR1 <- as.numeric(forest_fit_info$HR)
forest_fit_info$lower1 <- as.numeric(forest_fit_info$lower)
forest_fit_info$upper1 <- as.numeric(forest_fit_info$upper)
# 根据HR从小到大排个序
forest_fit_info <- forest_fit_info[order(forest_fit_info$HR1), ]
forest_fit_info <- rbind( c( 'P.value', 'HR','','', '95% CI','MODELS',NA,NA,NA), forest_fit_info)
# 再转成数值型用于图形映射
forest_fit_info$HR1 <- as.numeric(forest_fit_info$HR1)
forest_fit_info$lower1 <- as.numeric(forest_fit_info$lower1)
forest_fit_info$upper1 <- as.numeric(forest_fit_info$upper1)
# 定序
forest_fit_info$order <- c(19:1)
# 画画图
forestplot_cox <-
ggplot(forest_fit_info, aes(y = order))+
# 画画三线图，画画 x = 1 那条虚线
geom_segment(aes(x =-0.15, y = 18.5, xend = 1.05, yend = 18.5),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =-0.15, y = 19.5, xend = 1.05, yend = 19.5),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =-0.15, y = 0, xend = 1.05, yend = 0),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =1, y = 0, xend = 1, yend = 18.4),color="#DA3B22",size=0.5,linetype=3)+
# 画画误差线
geom_point(aes(x= HR1), shape=15, size=2)+
geom_linerange(aes(xmin= lower1, xmax= upper1))+
# 搞搞误差线两边的小齿
geom_segment(aes(x =lower1, y =order-0.15 , xend = lower1, yend = order+0.15),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =upper1, y =order-0.15 , xend = upper1, yend = order+0.15),color="black",size=0.5,linetype=1)+
# 搞搞坐标轴的小齿
geom_segment(aes(x =0.7, y = 0 , xend = 0.7, yend = -0.25),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =0.8, y = 0 , xend = 0.8, yend = -0.25),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =0.9, y = 0 , xend = 0.9, yend = -0.25),color="black",size=0.5,linetype=1)+
geom_segment(aes(x =1, y = 0 , xend = 1, yend = -0.25),color="black",size=0.5,linetype=1)+
# 搞搞坐标轴的值
geom_text(aes(x = 0.68,y = -0.7, label = 0.7), hjust = 0, fontface = 'bold', family = 'Times New Roman')+
geom_text(aes(x = 0.78,y = -0.7, label = 0.8), hjust = 0, fontface = 'bold', family = 'Times New Roman')+
geom_text(aes(x = 0.88,y = -0.7, label = 0.9), hjust = 0, fontface = 'bold', family = 'Times New Roman')+
geom_text(aes(x = 0.9975,y = -0.7, label = 1), hjust = 0, fontface = 'bold', family = 'Times New Roman')+
# 搞搞表格的宽度
# 搞搞坐标的表格内容
coord_cartesian(xlim=c(-0.1, 1))+
geom_text(aes(x = -0.1, label = labels), hjust = 0, fontface = 'bold', family = 'Times New Roman')+
geom_text(aes(x = 0.1, label = HR),hjust = 0,
fontface = ifelse(forest_fit_info$HR == "HR", "bold", "plain"),family = 'Times New Roman')+
geom_text(aes(x =0.25, label = CI95),hjust = 0,
fontface = ifelse(forest_fit_info$CI95 == "95% CI", "bold", "plain"), family = 'Times New Roman')+
geom_text(aes(x = 0.5, label = P.value),hjust = 0,
fontface = ifelse(forest_fit_info$P.value == "P.value", "bold", "plain"), family = 'Times New Roman')+
# 调调宽窄
scale_x_continuous(expand = c(-0.5,1))+
# 搞搞主题
theme_classic()+
theme(axis.line.y = element_blank(),
axis.ticks.y= element_blank(),
axis.text.y= element_blank(),
axis.title.y= element_blank(),
axis.line.x = element_blank(),
axis.ticks.x= element_blank(),
axis.text.x= element_blank(),
axis.title.x= element_blank())
# 存成矢量图，想要多大有多大
ggsave('results/forestplot_cox.svg',forestplot_cox)
3. 关于代码
上述代码容易看出存在大量的重复片段，最近我还和一个作现实数据模拟的朋友讨论过，我的代码就像一块块砖块整整齐齐的，包含着大量重复片段。他的代码则是我们常常为之称道的“看不太懂的”代码风格。我们讨论了一下，觉得可能是因为我目前的分析任务所需的代码量其实并不大。但需要在分析过程中根据数据情况和结果调整之前的代码。而他做模拟是往往没有真实数据，因此可以直接用逻辑持续写下去。在逻辑上排 bug 即可，但这不适用于我目前的工作。
最开始写的时候，因为之前的学习经历，总想用蹩脚的循环判断、逻辑与或非、排序来解决，动不动就循环嵌套，恨不能一行代码搞定所有事情。但逐渐就发现这种写法在排 bug 和调整参数的时候，十分痛苦。后来转变成上述风格之后，效率反而有提升。这算是这是关于上述不太高明的代码风格的辩解（hhh），但我能肯定的是，它一定准确表达了我的想法（hhh）。
全文完，感谢浏览。
“这就是本次的全部内容啦。我期望我虽然是一个搬运工，但给了它们一点额外的修饰，以助我们大家稍微拓宽练习的边界。
作者介绍