大家好！欢迎关注小号：医学统计数据分析，今天我们来介绍一下医学统计学中常用统计学方法的R语言基本操作。
R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个自由、免费、源代码开放的软件，它是一个用于统计计算和统计制图的优秀工具。
我们先下载R与Rstudio的安装包，先安装R，再安装Rstudio。RStudio是一款R语言的集成开发环境（IDE），R自带的环境操作起来可能不是方便，而Rstudio很好地解决了这个问题，RStudio只是辅助你使用R进行编辑的工具，因为它自身并不附带R程序。
R下载地址(官网)：
https://cran.r-project.org/bin/windows/base/
Rstudio(官网):
https://rstudio.com/
打开RStudio之后，会出现上图所示的窗口，其中有四个独立的面板。RStudio界面分为左上角的源码编辑、脚本显示，左下角的代码执行、控制台，右上角的代码历史记录、数据对象列表，右下角的代码组织管理、包安装、更新、绘图。
一、R导入Excel数据
我们先把需要分析的Excel表格另存为.csv格式，使用<-read.csv("")语句读取待分析的表格。
数据读取后，在脚本显示窗口可见数据预览：
使用summary()语句，对导入的数据大致分布做一下基本了解：
二、单因素（方差分析与卡方检验）
使用 anova<-aov(研究因素~分组因素)
summary(anova) 语句，可得方差分析结果：DF/SS/MS/F值/P值
使用plot(分组因素,研究因素)语句，可画出三组比较的箱形图：
使用chisq.test(分组因素,研究因素)做卡方检验，可见下图输出：
三、相关与回归
使用cor.test(testdata$1,testdata$2)语句做相关性分析，使用plot(因素1,因素2)语句画出散点图。
使用> lm(Y~X+1)语句及> summary(lm(Y~X+1))语句得到一元线性回归的B值、R2、残差分析等结果。
四、时间序列图等
根据整理好的时间序列资料，使用plot()语句，即可直接画出时间序列图（散点）：
安装包：
install.packages("zoo")
install.packages("xts")
载入包：
library(zoo)
library(xts)
计算与画图：
ts<- xts(data3$分析数据, as.Date(data3$日期, format='%Y/%m/%d'))
plot(ts)
xts可画出线图：
同样，我们也可以画出五日平均线等平滑曲线：
install.packages(pROC) # 下载 pROC 包
install.packages("ggplot2") # 下载 ggplot2 包
library(pROC) # 加载pROC包
library(ggplot2) # 调用ggplot2包以利用ggroc函数
roc1 <- roc(data1$结局,data1$`指标1`);roc1
roc2 <- roc(data1$结局,data1$`指标2`);roc2
roc3 <- roc(data1$结局,data1$`指标3`);roc3
roc4 <- roc(data1$结局,data1$`指标4`);roc4
roc5 <- roc(data1$结局,data1$`指标5`);roc5
roc6 <- roc(data1$结局,data1$`指标6`);roc6
roc7 <- roc(data1$结局,data1$`指标7`);roc7
plot(roc1,
max.auc.polygon=FALSE, # 填充整个图像
smooth=F, # 绘制不平滑曲线
main="Comparison of ROC curves", # 添加标题
col="red", # 曲线颜色为红色
legacy.axes=TRUE) # 使横轴从0到1，表示为1-特异度
plot.roc(roc2,
add=T, # 增加曲线
col="orange", # 曲线颜色为橙色
smooth = F) # 绘制不平滑曲线
plot.roc(roc3,
add=T, # 增加曲线
col="yellow", # 曲线颜色为黄色
smooth = F) # 绘制不平滑曲线
plot.roc(roc4,
add=T, # 增加曲线
col="green", # 曲线颜色为绿色
smooth = F) # 绘制不平滑曲线
plot.roc(roc5,
add=T, # 增加曲线
col="skyblue", # 曲线颜色为青色
smooth = F) # 绘制不平滑曲线
plot.roc(roc6,
add=T, # 增加曲线
col="blue", # 曲线颜色为兰色
smooth = F) # 绘制不平滑曲线
plot.roc(roc7,
add=T, # 增加曲线
col="purple", # 曲线颜色为紫色
smooth = F) # 绘制不平滑曲线
legend("bottomright", legend = c("指标1", "指标2","指标3","指标4","指标5","指标6","指标7"),col = c("red", "orange","yellow","green","skyblue","blue","purple"),
lwd = 1, cex = 0.5)
有时候我们做了一些单变量的ROC曲线，又做了一些联合指标的ROC曲线，判断联合指标是否在结局预测方面优于单变量，只看AUC的大小的话，也可以；同样我们也可以对这两条ROC曲线进行两两比较的统计学检验。那么怎么操作呢？下面我们用R语言演示一下。
Test for two correlated ROC curves，主要有三种方法，分别为："delong"、“bootstrap”或“venkatraman”。R语言实现方法分别为：
roc.test(roc1,roc2,method = "delong")
roc.test(roc1,roc2,method = "bootstrap")
roc.test(roc1,roc2,method = "venkatraman")
结果输出如下：
我们先画出roc1的图形：
plot(roc1,
max.auc.polygon=FALSE, # 填充整个图像
smooth=F, # 绘制不平滑曲线
main="Comparison of ROC curves", # 添加标题
col="red", # 曲线颜色
legacy.axes=TRUE) # 使横轴从0到1，表示为1-特异度
添加roc2的图形：
plot.roc(roc2,
add=T, # 增加曲线
col="orange", # 曲线颜色为红色
smooth = F) # 绘制不平滑曲线
比较两组曲线有无统计学差异：
testobj <- roc.test(roc1,roc2) # 检验两条曲线
text(0.8, 0.2, labels=paste("P value =", format.pval(testobj$p.value)), adj=c(0, .5)) # 在图上添加P值
legend(0.35,0.30, # 图例位置
bty = "n", # 图例样式，默认为 "o"
title="", # 引号内添加图例标题
legend=c("roc1","roc2"), # 添加分组
col=c("red","orange"), # 颜色跟前面一致
lwd=2) # 线条粗细
可以得到如下图形：
我们也可以通过更多的参数设置，增加曲线下面积及置信区间、最佳截断值、两两比较P值等。
详细绘制roc1:
plot(roc1,
print.auc=TRUE, print.auc.x=0.4, print.auc.y=0.5,
# 图像上输出AUC值,坐标为（x，y）
auc.polygon=TRUE, auc.polygon.col="#fff7f7", # 设置ROC曲线下填充色
max.auc.polygon=FALSE, # 填充整个图像
grid=c(0.5, 0.2), grid.col=c("black", "black"), # 设置间距为0.1，0.2，线条颜色
print.thres=TRUE, print.thres.cex=0.9, # 图像上输出最佳截断值，字体缩放倍数
smooth=F, # 绘制不平滑曲线
main="Comparison of two ROC curves", # 添加标题
col="red", # 曲线颜色
legacy.axes=TRUE) # 使横轴从0到1，表示为1-特异度
详细绘制roc2:
plot.roc(roc2,
add=T, # 增加曲线
col="orange", # 曲线颜色为红色
print.thres=TRUE, print.thres.cex=0.9, # 图像上输出最佳截断值，字体缩放倍数
print.auc=TRUE, print.auc.x=0.4,print.auc.y=0.4,
# 图像上输出AUC值,坐标为（x，y）
smooth = F) # 绘制不平滑曲线
详细绘制roc1VSroc2比较的P值:
testobj <- roc.test(roc1,roc2) # 检验两条曲线
text(0.8, 0.2, labels=paste("P value =", format.pval(testobj$p.value)), adj=c(0, .5)) # 在图上添加P值
legend(0.35,0.30, # 图例位置
bty = "n", # 图例样式，默认为 "o"
title="", # 引号内添加图例标题
legend=c("roc1","roc2"), # 添加分组
col=c("red","orange"), # 颜色跟前面一致
lwd=2) # 线条粗细
需要注意的是，有的ROC曲线根据坐标位置的约登指数大小，可以有两个或多个最佳截断值。
同样，我们也可以把roc3和roc6进行比较：
森林图是一种以统计指标和统计分析方法为基础，通过数值运算结果绘制出的图型。它主要用于综合展示每个被纳入研究的效应量以及汇总的合并效应量，帮助研究者直观地理解每个研究对总体效应的贡献。
森林图在医学研究中有着广泛的应用，包括但不限于：
Meta分析：森林图是Meta分析中最常用的结果综合表达形式，用于展示多个研究结果的汇总效应量及其可信区间。
亚组分析：除了在Meta分析中的应用，森林图还在临床试验的亚组分析中发挥作用，帮助研究者识别不同亚组之间的差异。
通过森林图，研究人员可以清晰地看到每个研究的效应量大小及其95%可信区间，以及所有研究合并后的效应量及可信区间，这对于评估干预措施的有效性、比较不同治疗方法的效果等方面具有重要意义。此外，森林图的使用不仅限于医学领域，还可以应用于其他学科的研究中，只要涉及到多个研究的汇总分析，都可以利用森林图来展示结果。
今天我们使用R语言的forestplot（）包演示一下森林图的画法和美化。
例如我们有如下待分析数据：
我们使用SPSS打开该EXCEL数据库
使用卡方检验对分类变量进行比较分析
对定量变量进行独立样本t检验的比较分析
将需要分析的变量纳入logistic回归，查看输出结果。
SPSS目前版本没有绘制森林图的功能，那么我们怎么用logistic回归结果画出森林图呢？
我们先把数据整理成如下格式，包括OR(95% CI）、OR_1、 OR_2、 OR_mean、P等。
然后打开Rstudio,
#安装加载森林图所需的包
install.packages(forestplot","grid","checkmate","abind)
library(forestplot)
library(grid)
library(checkmate)
library(abind)
#读取数据（Excel或CSV格式）
fdata <- read_excel("路径/森林图数据.xlsx")
fdata<-read.csv("路径/森林图数据.csv",header=T)
#运行得出森林图 Export 导出为图片或PDF格式
forestplot(labeltext=as.matrix(fdata[,1:3]),
mean=fdata$OR_mean,
lower=fdata$OR_1,
upper=fdata$OR_2,
zero=1,
boxsize=0.2,
graph.pos=2)
#对输出图形进行美化优化
forestplot(labeltext=as.matrix(fdata[,1:3]),
mean=fdata$OR_mean,
lower=fdata$OR_1,
upper=fdata$OR_2,
zero=1,
boxsize=0.2,
lineheight = unit(7,'mm'),
colgap=unit(2,'mm'),
lwd.zero=1,
lwd.ci=0.5,
col=fpColors(box='#458B00',
summary='#8B008B',
lines = 'black',
zero = '#7AC5CD'),
xlab="OR",
lwd.xaxis =1,
txt_gp = fpTxtGp(ticks = gpar(cex = 0.85),
xlab = gpar(cex = 0.8),
cex = 0.9),
lty.ci = "solid",
title = "Forestplot",
line.margin = 0.08,
graph.pos=2)
DCA曲线（Decision Curve Analysis）是一种评估临床预测模型实用性和临床应用价值的方法，它通过绘制决策曲线来比较不同预测模型在不同治疗决策阈值下的临床效益。DCA的核心思想是将患者或决策者的偏好整合到分析中，从而更全面地评估预测模型的实用性和临床应用价值。
一、准备数据：首先，你需要准备用于分析的数据。这些数据应该包括患者的相关信息，如诊断结果、治疗选择、以及可能影响诊断结果的各项指标。数据需要以数字形式表示，并且需要包含一个结果变量（例如，是否患病）和多个预测变量（例如，年龄、病症阶段、基因表达等）。
二、选择合适的软件和R包：对于使用R语言进行DCA分析，可以选择rmda包或其他相关包（如ggDCA）来帮助绘制DCA曲线。确保你已经安装了必要的R包，并且了解如何加载和使用这些包。
三、构建模型：使用选定的R包构建DCA模型。这通常涉及到指定模型公式、设置阈值范围以及进行必要的参数调整。例如，在使用rmda包时，你可能需要指定结果变量和预测变量，设置阈值范围，并进行一定数量的bootstrap重抽样来评估模型的稳定性。
四、绘制DCA曲线：一旦模型构建完成，就可以使用R包的绘图功能来生成DCA曲线。这些曲线将显示在不同治疗决策阈值下，使用预测模型进行治疗决策所获得的利益。利益可以定义为治疗正确的病例数减去治疗错误的病例数，其中治疗正确指的是根据预测模型的输出结果进行正确的治疗决策。
五、解释结果：最后，解释DCA曲线的结果。比较不同预测模型在各个阈值下的临床效益，确定哪个模型在不同治疗决策阈值下具有更高的临床效益。这有助于在不同预测模型之间做出最优选择，以达到更好的临床预测效果。
今天我们使用R语言的rmda包来演示一下DCA曲线的画法。
比如我们有如下图数据（同ROC曲线绘制数据）：
rm(list=ls()) #清空所有数据变量列表
install.packages() #安装包
library() #加载包
#载入数据集及安装包
install.packages(c("rmda","readxl")
library(readxl)
library(rmda)
ddata <- read_excel("C:/Users/hyy/Desktop/DCA数据.xlsx")
#将数据读取到R
#DCA参数运算
model1 <- decision_curve(结局~指标1,
family=binomial(link='logit'),
data = ddata,
thresholds = seq(0, .8, by = .05),
bootstraps = 10)
#绘制DCA曲线
plot_decision_curve(model1,
curve.names = "Model1",
col = c("red"),
confidence.intervals=T #不显示置信区间
)
#我们可以得到模型1的DCA曲线
#同样，我们可以对其他模型指标组合进行运算和绘图
model2 <- decision_curve(结局~指标1+指标2,
family=binomial(link='logit'),
data = ddata,
thresholds = seq(0, .8, by = .05),
bootstraps = 10)
plot_decision_curve(model2,
curve.names = "Model2",
col = c("blue"),
confidence.intervals=T #不显示置信区间
)
那么我们如何将这两个模型的曲线绘制到一起呢？
我们用list() 和c() 列表组合功能，将曲线绘制到一起。
#绘制两条DCA曲线
plot_decision_curve(list(model1, model2),
curve.names = c("model1", "model2"),
col = c("red", "blue"),#设置线的颜色
lty = c(1,2), # 设置线形
lwd = c(3,3, 2, 2),
confidence.intervals = F,
legend.position = "topright")
DCA曲线与ROC曲线都是用于评估诊断或预测模型的工具，它们共同关注于模型的性能评估。
DCA曲线（Decision Curve Analysis）和ROC曲线（Receiver Operating Characteristic）都是用来评估诊断或预测模型的性能。DCA曲线通过衡量模型在不同阈值下的净获益来评价模型的临床效用，将患者或决策者的偏好整合到分析中。ROC曲线则通过绘制真正例率（True Positive Rate）与假正例率（False Positive Rate）之间的关系，来评估模型的诊断准确性。
两者都旨在提供关于模型性能的直观理解，帮助决策者选择最佳模型。DCA曲线通过考虑不同阈值下的净获益，为临床决策提供更实用的指导；而ROC曲线则通过提供不同分类阈值下的性能指标，帮助评估模型的诊断能力。
DCA曲线与ROC曲线的主要区别在于它们的评估重点和应用场景不同。
评估重点不同：
ROC曲线（Receiver Operating Characteristic curve）主要关注于模型的区分度，通过绘制真正例率（True Positive Rate, TPR）与假正例率（False Positive Rate, FPR）之间的关系来评估模型的性能。它衡量的是模型在不同阈值下的表现，特别是在特异性和敏感性方面的表现，追求的是准确性。
DCA曲线（Decision Curve Analysis）则更侧重于临床决策的实际应用，将患者或决策者的偏好整合到分析中，考虑特定模型的临床效用。它不仅评估模型的诊断准确性，还考虑了模型的临床价值，即在特定情况下模型是否值得使用，以及如何选择最佳模型。
应用场景不同：
ROC曲线广泛应用于各种领域的模型评估，尤其是在需要评估分类模型的性能时，如疾病诊断、欺诈检测等，它提供了一个全面的视角来理解模型在不同阈值下的表现。
DCA曲线更适合于临床预测模型的评估，特别是在需要权衡假阳性和假阴性对临床决策影响的情况下。例如，在癌症筛查中，DCA可以帮助医生在不同筛查模型之间做出最优选择，以提高筛查的准确性和效率。
综上所述，ROC曲线和DCA曲线各有其适用场景和优势。ROC曲线提供了一个客观的、量化准确性的评估，而DCA曲线则更侧重于将临床决策者的偏好和需求纳入考虑，以评估模型的实际临床应用价值。
在变量间依存关系分析中，当因变量是连续型的定量变量且自变量与因变量间具有线性关系时我们常采用多重线性回归进行分析。
然而，在生物医学研究领域，特别是流行病学研究中，我们经常需要探索某疾病的发病原因或其影响因素，在医学研究中的碰到的因变量(结局指标，或称效应指标)通常是分类变量例如，“生存或死亡”“发病与不发病"等，并且影响因素(自变量)与其联系更多的是非线性关系，对于这类问题，应用线性回归显然缺乏合理性。
在二分类因变量数据的依存关系分析中，常用的方法包括logistic回归分析、Probit回归分析以及混合效应模型。这些方法适用于因变量只有两种可能结果（例如，患病与未患病、离职与未离职等）的情况。
今天我们分享一下用R语言各种程序包，进行二分类因变量数据相关分析的较全的过程，包括独立样本t检验、卡方检验、logistic回归、ROC曲线绘制、列线图绘制等。
首先我们打开上的测试数据库，因变量结局为二分类变量，自变量中有定量变量和分类变量，如下图：
#加载数据包并读取数据
install.packages(caret","lattice","glmnet",
"Matrix","pROC","Hmisc","rms)
library(ggplot2) #绘图用包
library(caret) #分层抽样拆分数据库
library(lattice) #绘图可视化包
library(readxl) #读取Excel包
library(gmodels) #卡方检验包
library(glmnet) #回归分析用包
library(Matrix) #矩阵运算用包
library(pROC) #ROC曲线用包
library(Hmisc) #列线图用包
library(rms) #列线图用包
data <- read_excel("C:/Users/hyy/Desktop/测试数据.xlsx")
可将数据在RStudio中打开：
#独立样本t检验
t.test(data$指标1 ~ 结局, data = data, var.equal = TRUE)
t.test(data$指标2 ~ 结局, data = data, var.equal = TRUE)
t.test(data$指标3 ~ 结局, data = data, var.equal = TRUE)
#卡方检验
CrossTable(data$结局,data$指标8,expected = T,chisq = T,fisher = T, mcnemar = T, format = "SPSS")
#logistic回归
data$Group<-as.factor(data$结局)
model1 <- glm(Group ~ 指标1, data = data, family = "binomial")
summary(model1)
model2 <- glm(Group ~ 指标1+指标2, data = data, family = "binomial")
summary(model2)
model3 <- glm(Group ~ 指标1+指标2+指标3, data = data, family = "binomial")
summary(model3)
#ROC曲线
roc1 <- roc(data$结局,data$指标1);roc1
roc2 <- roc(data$结局,data$指标2);roc2
roc3 <- roc(data$结局,data$指标3);roc3
plot(roc1,
max.auc.polygon=FALSE, # 填充整个图像
smooth=F, # 绘制不平滑曲线
main="Comparison of ROC curves", # 添加标题
col="red", # 曲线颜色
legacy.axes=TRUE) # 使横轴从0到1，表示为1-特异度
plot.roc(roc2,
add=T, # 增加曲线
col="orange", # 曲线颜色为红色
smooth = F) # 绘制不平滑曲线
plot.roc(roc3,
add=T, # 增加曲线
col="yellow", # 曲线颜色为红色
smooth = F) # 绘制不平滑曲线
#列线图
dd<-datadist(data)
options(datadist="dd")
data$Group<-as.factor(data$结局)
f_lrm<-lrm(Group~指标1+指标2+指标3+指标4+指标5+指标6+指标7+指标8,data=data)
summary(f_lrm)
par(mgp=c(1.6,0.6,0),mar=c(5,5,3,1))
nomogram <- nomogram(f_lrm,fun=function(x)1/(1+exp(-x)),
fun.at=c(0.01,0.05,0.2,0.5,0.9,0.99),
funlabel =data$结局,
conf.int = F,
abbrev = F )
plot(nomogram)
校准曲线是用于评价临床预测模型性能的一种图形工具，它反映了模型预测的绝对风险与实际观察到的风险之间的匹配程度。校准曲线通过对比不同分位数上的预测概率与观察概率，来评估预测概率的准确性。一个理想的校准曲线应该与斜率为1的理想线重合，表明模型预测的结果与实际结果高度一致。如果校准曲线越接近理想线，说明模型的预测结果越可靠，模型的校准性能越好。
校准曲线的绘制基于已知的实际观察数据和模型预测的数据。首先，收集一系列已知被测量值和对应的测量值的数据，然后将这些数据以散点图的形式表示，横坐标为已知的被测量值，纵坐标为相应的测量值。接着，使用适当的拟合方法（如线性回归、多项式回归等）将散点图上的点拟合为一条曲线，并通过最小二乘法等优化算法确定拟合曲线的参数。最后，在校准曲线上标注参数信息，如斜率、截距、相关系数等，同时注明被测量的范围和单位等信息。
在临床预测模型的评估中，校准曲线是一种重要的工具，用于确保模型的预测结果在实际应用中是可靠和有效的。通过校准曲线的分析，可以判断模型是否高估或低估了实际风险，从而帮助改进模型或提供更准确的预测。
Calibration curve，直译过来就是校准曲线或校准图。其实，校准曲线就是实际发生率和预测发生率的散点图。实质上，校准图曲线是Hosmer-Lemeshow拟合优度检验的结果可视化。目前校准曲线常用来评价logistic回归和cox回归模型。
#安装并加载包及数据
install.packages(readxl)#安装包
install.packages("rms")#安装包
library(readxl) #读取Excel包
library(rms) #校准曲线图用包
ceshi <- read_excel("C:/Users/hyy/Desktop/测试数据.xlsx")
#绘图所需变量转换
ceshi <-na.omit(ceshi)
dd<-datadist(ceshi)
options(datadist="dd")
model1 <- glm(结局 ~ 指标1+指标2+指标3+指标4+指标5+指标6+指标7+指标8, data = ceshi,family = "binomial")
summary(model1)
#对测试数据进行校准曲线绘制
ceshi$结局<-as.factor(ceshi$结局)
fit1<- lrm(结局~指标1+指标2+指标3+指标4+指标5,data=ceshi,x=TRUE,y=TRUE)
cal1<-calibrate(fit1,method="boot",B=1000)
plot(cal1)
#对测试数据进行校准曲线绘制
ceshi$结局<-as.factor(ceshi$结局)
fit2<- lrm(结局~指标1+指标2+指标3+指标4+指标5+指标6+指标7+指标8,data=ceshi,x=TRUE,y=TRUE)
cal2<-calibrate(fit2,method="boot",B=1000)
plot(cal2)
#对测试数据进行校准曲线绘制
ceshi$结局<-as.factor(ceshi$结局)
fit3<- lrm(结局~指标2+指标6+指标7+指标8,data=ceshi,x=TRUE,y=TRUE)
cal3<-calibrate(fit3,method="boot",B=1000)
plot(cal3)
