R语言可视化及作图13--nomogram和生存曲线

R语言绘图系列：

1. nomogram

library(rms)
n <- 1000 #define sample size
d <- data.frame(age=rnorm(n,50,10),
                blood.pressure=rnorm(n,120,15),
                cholesterol=rnorm(n,200,25),
                sex=factor(sample(c('female','male'),n,TRUE)))
#Specify population model for log odds that Y=1
#Simulate binary y to have prob(y=1)=1/[1+exp(-L)]
d <- upData(d,
            L=.4*(sex=='male')+.045*(age-50)+
              (log(cholesterol-10)-5.2)*(-2*(sex=='female')+2*(sex=='male')),
            y=ifelse(runif(n) < plogis(L),1,0))
ddist <- datadist(d);options(datadist='ddist')
#这一步是必不可少的
#The datadist object that was in effect when the model was fit is used to specify the limits of the axis
#for continuous predictors when the user does not specify tick mark locations in the nomogram call.
f <- lrm(y~lsp(age,50)+sex*rcs(cholesterol,4)+blood.pressure,
         data=d) #logistic regression
nom <- nomogram(f,fun = function(x)1/(1+exp(-x)), #or fun=plogis,logistic distribution
                fun.at = c(.001,.01,.05,seq(.1,.9,by=.1),.95),
                #fun.at:function values to label on axis
                funlabel = 'Risk of Death')
#lp:是否显示liner predictor
#fun:对linear predictor进行转化，并声称一条新轴
#Instead of fun.at, could have specified fun.lp.at=logit of
#sequence above-faster and slightly more accurate
plot(nom,xfrac=.35)
print(nom)

nom <- nomogram(f,age=seq(10,90,by=10))#对age进行重新定义
plot(nom,xfrac = .45)

nom <- nomogram(f,age=seq(10,90,by=10))#对age进行重新定义
plot(nom,xfrac = .45)

g <- lrm(y~sex+rcs(age,3)*rcs(cholesterol,3),data=d)
#对胆固醇变量进行了rcs转换
nom <- nomogram(g,interact = list(age=c(20,40,60)),#列出年龄在20，40，60的情况
                conf.int = F)
plot(nom,col.conf = c(1,.5,.2),naxes = 7)
#naxes:至多在一张图上展示几条轴

w <- upData(d,
            cens=15*runif(n),
            h=.02*exp(.04*(age-50)+.8*(sex=='Female')),
            d.time=-log(runif(n))/h,
            death=ifelse(d.time <= cens,1,0),
            d.time=pmin(d.time,cens))
f <- psm(Surv(d.time,death)~sex*age,data=w,dist='lognormal')
med <- Quantile(f)
surv <- Survival(f) #This would also work if f was from cph
plot(nomogram(f,fun=function(x) med(lp=x),funlabel = 'Median Survival Time'))
nom <- nomogram(f,fun=list(function(x) surv(3,x),
                           function(x) surv(6,x)),
                funlabel = c('3-Month Survival Probability',
                             '6-Month Survival Probability'))
plot(nom,xfrac = .7)

2. 生存曲线

2.1 survival包

library(survival)
data('GBSG2',package='TH.data')
plot(survfit(Surv(time,cens)~horTh,data = GBSG2),lty=c(2,1),col=c(2,1),mark.time = F)
legend('bottomleft',legend=c('yes','no'),title = 'Hormonal Therapy',lty = c(2,1),col = c(2,1),bty = 'n')

2.2 survminer包 (基于ggplot2)

library(survminer)
library(ggeffects)
data('lung')
fit <- survfit(Surv(time,status)~sex,data = lung)
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 typr by  group
           surv.median.line='hv', #specify median survival
           ggtheme=theme_bw(), #change ggplot2 theme
           palette=c('#E7B800','#2E9FDF'))

ggsurvplot(
  fit, #survfit object with calculated statistics
  pval = TRUE, #show p-value of log-rank test
  conf.int = TRUE, #show confidence intervals for point estimates of survival curves
  conf.int.style='ribbon', #customize style of conficence intervals
  xlab='Time in days', #customize X axis label
  break.time.by=200, #break X axis in time intervals by 200
  ggtheme = theme_light(), #customize plot and risk table with a theme
  risk.table = 'abs_pct', #absolute number and percentage at risk
  risk.table.y.text.col=T, #colour risk table text annotations
  risk.table.y.text=FALSE, #show bars instead of names in text annotations in legend of risk table
  ncensor.plot=TRUE, #plot the number of censored subjects at time t
  surv.median.line = 'hv', #add the median survival pointer
  legend.labs=c('Male','Female'), #change legend labels
  palette = c('#E7B800','#2E9FDF') #custom color palettes
)

ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col='strata',
           ggtheme = theme_bw(),
           palette = c('#E7B800','#2E9FDF'),
           fun = 'event')

ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col='strata',
           ggtheme = theme_bw(),
           palette = c('#E7B800','#2E9FDF'),
           fun = 'cumhaz')

不同情况下的生存情况
fit2 <- survfit(Surv(time,status)~sex+rx+adhere,
                data=colon)
ggsurv <- ggsurvplot(fit2,fun = 'event',conf.int = TRUE,ggtheme = theme_bw())
ggsurv $plot+theme_bw()+theme(legend.position = 'right')+facet_grid(rx~adhere)

R语言可视化及作图13--nomogram和生存曲线

1. nomogram

2. 生存曲线

推荐阅读更多精彩内容