#########數(shù)據(jù)處理
#1.載入R包
BiocManager::install("grid")
library(survival)
library(plyr)
#2.清理工作環(huán)境
rm(list = ls())
#3.讀入數(shù)據(jù)
setwd("C:\\Users\\xx\\xx")
aa<- read.csv('D-V1.csv') #文件名稱
#4.查看數(shù)據(jù)數(shù)據(jù)性質(zhì)
str(aa)
#5.查看結(jié)局，0=生存，1死亡
aa$status<-factor(aa$removal_state)
summary(aa$status)



#############批量單因素回歸
#1.構(gòu)建模型的y
y<- Surv(time=aa$removal_time,event=aa$removal_state==1)
#2.批量單因素回歸模型建立：Uni_cox_model
Uni_cox_model<- function(x){
  FML <- as.formula(paste0 ("y~",x))
  cox<- coxph(FML,data=aa)
  cox1<-summary(cox)
  HR <- round(cox1$coefficients[,2],2)
  PValue <- round(cox1$coefficients[,5],3)
  CI5 <-round(cox1$conf.int[,3],2)
  CI95 <-round(cox1$conf.int[,4],2)
  Uni_cox_model<- data.frame('Characteristics' = x,
                            'HR' = HR,
                            'CI5' = CI5,
                            'CI95' = CI95,
                            'p' = PValue)
  return(Uni_cox_model)} 

#3.將想要進(jìn)行的單因素回歸變量輸入模型
#3-(1)查看變量的名字和序號(hào)
names(aa)
#3-(2)輸入變量序
variable.names<- colnames(aa)[c(3:23)] #例這里選擇了3-23號(hào)變量
#4.輸出結(jié)果
Uni_cox <- lapply(variable.names, Uni_cox_model)
Uni_cox<- ldply(Uni_cox,data.frame)
#5.優(yōu)化表格，這里舉例HR+95% CI+P 風(fēng)格
Uni_cox$CI<-paste(Uni_cox$CI5,'-',Uni_cox$CI95)
Uni_cox<-Uni_cox[,-3:-4]
#查看單因素cox表格
View(Uni_cox)
write.csv(Uni_cox,file = "Uni_cox.csv",row.names = F)

若單因素分析中有顯著性差異因素過多，可進(jìn)行共線性診斷，剔除掉存在共線性的因素。
進(jìn)行共線性診斷分析的方法有許多，如特征值法；條件指數(shù)法和方差膨脹因子法。雖然方法不同，但是判定結(jié)果是一致的。下面以方差膨脹因子vif舉例。

###共線性診斷###
#6.共線性診斷-VIF值判斷法
#先安裝car包
options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
install.packages("BiocManager")
BiocManager::install("car")
library(car)
#進(jìn)行共線性診斷分析，~前觀察生存時(shí)間，~后為單因素中有顯著性差異的因素
fit<-lm(aa$removal_time ~ aa$IIRC+ aa$blood_qf+ aa$corneal_edema
        + aa$neovascularization_iris+ aa$high_pressure+ aa$IAC 
        + aa$local_therapy + aa$cal+ aa$intravitreal_hemorrhosis,data=aa)
vif(fit)
###VIF大于10表明與其他因素存在共線性。根據(jù)臨床專業(yè)判斷，先去掉一個(gè)不太重要的變量，再看下共線性，不行的話，就再去掉一個(gè)

將不存在共線性的變量納入多因素分析中

######多因素回歸
install.packages(c("survival", "survminer"))
library(survival)
library(survminer)
head(aa)#查看示例數(shù)據(jù)

#多因素分析
multi_cox <- coxph(Surv(removal_time, removal_state) ~ IIRC + blood_qf+
                  high_pressure+IAC+local_therapy+intravitreal_hemorrhosis+cal, 
                  data =  aa)
summary(multi_cox)
ggforest(multi_cox,
         data = aa,
         main = "Hazard ratio",
         cpositions = c(0.02, 0.22, 0.4),#前三列的位置
         fontsize = 0.7,#字體大小
         refLabel = "reference",#顯示變量
         noDigits = 3)#小數(shù)點(diǎn)后保留位數(shù)