R語言繪圖系列:
- R語言可視化及作圖1--基礎(chǔ)繪圖(par函數(shù)户矢,散點(diǎn)圖,盒形圖,條形圖避消,直方圖)
- R語言可視化及作圖2--低級繪圖函數(shù)
- R語言可視化及作圖3--圖形顏色選取
- R語言可視化及作圖4--qplot和ggplot2美學(xué)函數(shù)
- R語言可視化及作圖5--ggplot2基本要素和幾何對象匯總
- R語言可視化及作圖6--ggplot2之點(diǎn)圖碳锈、條形圖顽冶、盒形圖、直方圖售碳、線圖
- R語言可視化及作圖7--ggplot2之標(biāo)簽渗稍、圖例和標(biāo)題繪制
- R語言可視化及作圖8--坐標(biāo)軸自定義和坐標(biāo)系轉(zhuǎn)換
- R語言可視化及作圖9--主題函數(shù)
- R語言可視化及作圖10--ggplot2的theme函數(shù)
- R語言可視化及作圖11--圖片分面函數(shù)和一頁多圖
- R語言可視化及作圖12--venn圖、火山圖和熱圖
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進(jìn)行轉(zhuǎn)化团滥,并聲稱一條新軸
#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進(jìn)行重新定義
plot(nom,xfrac = .45)
nom <- nomogram(f,age=seq(10,90,by=10))#對age進(jìn)行重新定義
plot(nom,xfrac = .45)
g <- lrm(y~sex+rcs(age,3)*rcs(cholesterol,3),data=d)
#對膽固醇變量進(jìn)行了rcs轉(zhuǎn)換
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)