R实战｜从文献入手谈谈logistic回归、Cox回归以及Lasso分析(一)

rm(list = ls())
data <- read.csv('data.csv')
head(data) 
summary(data)
# 将分类变量设置为factor，删除缺失值记录
## 因子化
VarsC<-c("sex","cp","fbs","restecg","exang","slope","ca","thal")#分类变量
for(i in VarsC){
  data[,i] <- as.factor(data[,i])
}#利用循环因子化
## 将hd(heart disease)中的数字（0~4）改成"Healthy" 和"Unhealthy"
data$hd <- ifelse(test=data$hd == 0, yes="Healthy", no="Unhealthy")
data$hd <- as.factor(data$hd) 
str(data)
# NA值处理：如果NA较少，可以直接删除；如果NA较多，需要进行缺失值处理
## 查看NA个数
table(is.na(data))
##展示含有NA的记录
data[!complete.cases(data),]
##去除含NA的样本
data <- na.omit(data)
#多因素logistic回归
## 构建Logistic模型
fitMul <- glm(formula = hd ~ ., data =data, family = 'binomial')
summary(fitMul) ##logistic回归结果

> head(ResultMul)
               Estimate Std. Error    z value     Pr(>|z|)           OR        2.5 %
(Intercept) -6.25397767 2.96039852 -2.1125459 0.0346396558  0.001922791 4.732245e-06
age         -0.02350753 0.02512158 -0.9357505 0.3494016427  0.976766619 9.290676e-01
sex1         1.67015202 0.55248644  3.0229738 0.0025030392  5.312975393 1.858511e+00
cp2          1.44839631 0.80913642  1.7900521 0.0734455437  4.256283277 8.875191e-01
cp3          0.39335276 0.70033787  0.5616614 0.5743467297  1.481941062 3.820686e-01
cp4          2.37328693 0.70909374  3.3469298 0.0008171191 10.732611757 2.823872e+00
                97.5 %
(Intercept)  0.5585898
age          1.0258173
sex1        16.4363241
cp2         21.7494748
cp3          6.0931562
cp4         46.5187371

#逐步回归
logit.step <- step(fitMul,direction="both")
summary(logit.step)

#重新建立回归模型
fitstep <- glm(formula = hd ~ sex + cp + trestbps + thalach + exang + oldpeak + slope + ca + thal, family = "binomial", data = data)
fitSum<-summary(fitstep)
ResultMul<-c()#准备空向量，用来储存结果
ResultMul<-rbind(ResultMul,fitSum$coef)
OR<-exp(fitSum$coef[,'Estimate'])
ResultMul<-cbind(ResultMul,cbind(OR,exp(confint(fitstep))))
# 导出结果
write.csv(ResultMul,file="Mul_log.csv")

# 批量单因素logistic
## 准备进行分析的自变量
varsU<-names(data[,1:13])#自变量
Result<-c()
for (i in 1:length(varsU)){
  fit<-glm(substitute(hd~x,list(x=as.name(varsU[i]))),data=data,family=binomial())
  fitSum<-summary(fit)
  result1<-c()
  result1<-rbind(result1,fitSum$coef)
  OR<-exp(fitSum$coef[,'Estimate'])
  result1<-data.frame(cbind(result1,cbind(OR,exp(confint(fit)))))
  result1$Characteristics<-varsU[i]   #添加变量名
  Result<-rbind(Result,result1[-1,])#[-1,],删除常数项
}
## 为R语言logistic回归默认将分类变量的第一个factor设置为参照
Uni_log<-data.frame(Result[,c(1,4:8)]) #提取"P","OR","CIlower","CIupper"和变量名
colnames(Uni_log)[2:5]<-c("P","OR","CIlower","CIupper")#变量重命名
head(Uni_log)
ExtractVar<-unique(Uni_log$Characteristics[Uni_log$"P"<0.05])#提取有意义的变量
write.csv(Uni_log,file="Uni_log.csv")#输出文档