R语言学习笔记
参考内容：https://www.bioinfo-scrounger.com/archives/647/ #注意该文章中，部分代码前面多了>符号
//www.greatytc.com/p/2da6645e0a86 #两篇文章采取的函数略有不同
library("survival")
library("survminer")
data("lung") #载入lung数据库
head(lung)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
* inst: Institution code #此处为数据的注释，无需理会
* time: Survival time in days
* status: censoring status 1=censored, 2=dead
* age: Age in years
* sex: Male=1 Female=2
* ph.ecog: ECOG performance score (0=good 5=dead)
* ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
* pat.karno: Karnofsky performance score as rated by patient
* meal.cal: Calories consumed at meals
* wt.loss: Weight loss in last six months
fit <- survfit(Surv(time, status) ~ sex, data = lung) #单因素分析
print (fit) #输出fit结果
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF")
) # 将以上的内容进行可视化，绘K-M图
以下是自己的数据（SHC）采用的代码
library("survival")
library("survminer")
#导入excel数据
fit <- survfit(Surv(time1, status1) ~ sex, data = training) #单因素分析k-m
print (fit) # 输出结果
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF")
) # 将以上的内容进行可视化，绘K-M图
#当数据量不够时，缺少绘图需要的必要元素，如95%置信区间，可能会报错，提示信息不全，如若对china_stage进行分析，则会报错
surv_diff <- survdiff(Surv(time1, status1) ~ sex, data = training) #进行log-rank检验
surv_diff #展示结果
fit2 <- coxph(Surv(time1, status1) ~ age + sex, data = training) #多因素分析cox
summary (fit2) #完整报告
#注意，K-M分析和cox会有细微差别，此次也可采用coxph函数做单因素分析
#对多个变量批量处理
covariates<-c("age","sex")
univ_formulas <- sapply(covariates, function(x) as.formula(paste('Surv(time1, status1)~', x)))
univ_formulas
univ_models<-lapply(univ_formulas,function(x){coxph(x,data=training)})
univ_models
univ_results<-lapply(univ_models,function(x){x<-summary(x)p.value<-signif(x$wald["pvalue"],digits=2)wald.test<-signif(x$wald["test"],digits=2)beta<-signif(x$coef[1],digits=2);#coeficient beta HR<-signif(x$coef[2],digits=2);#exp(beta)HR.confint.lower<-signif(x$conf.int[,"lower .95"],2)HR.confint.upper<-signif(x$conf.int[,"upper .95"],2)HR<-paste0(HR," (",HR.confint.lower,"-",HR.confint.upper,")")res<-c(beta,HR,wald.test,p.value)names(res)<-c("beta","HR (95% CI for HR)","wald.test","p.value")return(res)#return(exp(cbind(coef(x),confint(x))))})class(univ_results)##[1]"list"str(univ_results)
res<-t(as.data.frame(univ_results,check.names=FALSE))as.data.frame(res)
#提取结果，此段代码无需修改
