cox<-read.csv("경로/cox1.csv")
cox<-na.omit(cox)
head(cox)

#Cox 비례위험 회귀분석
install.packages("survival")
library(survival)

#치료방법별 사망률
r1<-coxph(Surv(time,event==1)~factor(clinic),data=cox)
summary(r1)

#성별 사망률
r2<-coxph(Surv(time,event==1)~factor(sex),data=cox)
summary(r2)

#연령별 사망률
r3<-coxph(Surv(time,event==1)~age,data=cox)
summary(r3)

#흡연여부별 사망률
r4<-coxph(Surv(time,event==1)~factor(smoking),data=cox)
summary(r4)

#폐활량 결과별 사망률
r5<-coxph(Surv(time,event==1)~lungtest,data=cox)
summary(r5)

#유의성 검증
rtotal<-coxph(Surv(time,event==1)~factor(clinic)+factor(sex)+age+factor(smoking)+lungtest,data=cox)
summary(rtotal)

#위험비 확인
install.packages("moonBook")
library(moonBook)

#시기별 위험비
cox$TS=Surv(cox$time,cox$event==1)
out=mycph(TS~factor(clinic)+factor(sex)+age+factor(smoking)+lungtest, data=cox)
out

#위험비 그래프
HRplot(rtotal, type=2, show.CI=TRUE, sig.level=0.05, main="Hazard ratios of all individual variables")

#논문에 자주 쓰이는 양식
install.packages("stargazer")
library(stargazer)

stargazer(r1,r2,r3,r4,r5,type="text",no.space = T)
stargazer(rtotal,type="text",no.space = T)

#생존곡선/위험곡선 시각화
install.packages("jskm")
library(jskm)

r11 <- survfit(Surv(cox$time,cox$event) ~ factor(clinic), data = cox)
jskm(r11,  cumhaz = T,  mark = T, ylab = "Cumulative hazard (%)", surv.scale = "percent", pval =T,table = T)