R語言學(xué)習(xí)筆記
參考內(nèi)容：https://www.bioinfo-scrounger.com/archives/647/?#注意該文章中，部分代碼前面多了>符號
? ? ? ? ? ? ? ? ? http://www.itdecent.cn/p/2da6645e0a86?#兩篇文章采取的函數(shù)略有不同
library("survival")
library("survminer")
data("lung") #載入lung數(shù)據(jù)庫
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??#此處為數(shù)據(jù)的注釋，無需理會
* 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結(jié)果
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")
? ? ? ) # 將以上的內(nèi)容進(jìn)行可視化，繪K-M圖
以下是自己的數(shù)據(jù)（SHC）采用的代碼
library("survival")
library("survminer")
#導(dǎo)入excel數(shù)據(jù)
fit <- survfit(Surv(time1, status1) ~ sex, data = training)? #單因素分析k-m
print (fit) # 輸出結(jié)果
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")
) # 將以上的內(nèi)容進(jìn)行可視化，繪K-M圖
#當(dāng)數(shù)據(jù)量不夠時，缺少繪圖需要的必要元素，如95%置信區(qū)間，可能會報錯，提示信息不全，如若對china_stage進(jìn)行分析，則會報錯
surv_diff <- survdiff(Surv(time1, status1) ~ sex, data = training) #進(jìn)行l(wèi)og-rank檢驗
surv_diff #展示結(jié)果
fit2 <- coxph(Surv(time1, status1) ~ age + sex, data = training) #多因素分析cox
summary (fit2) #完整報告
#注意，K-M分析和cox會有細(xì)微差別，此次也可采用coxph函數(shù)做單因素分析
#對多個變量批量處理
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)
#提取結(jié)果，此段代碼無需修改
