在变量间依存关系分析中,当因变量是连续型的定量变量且自变量与因变量间具有线性关系时我们常采用多重线性回归进行分析。
然而,在生物医学研究领域,特别是流行病学研究中,我们经常需要探索某疾病的发病原因或其影响因素,在医学研究中的碰到的因变量(结局指标,或称效应指标)通常是分类变量例如,“生存或死亡”“发病与不发病"等,并且影响因素(自变量)与其联系更多的是非线性关系,对于这类问题,应用线性回归显然缺乏合理性。
在二分类因变量数据的依存关系分析中,常用的方法包括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)

