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進(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)
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
禁止轉(zhuǎn)載灸姊,如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者拱燃。
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市力惯,隨后出現(xiàn)的幾起案子碗誉,更是在濱河造成了極大的恐慌召嘶,老刑警劉巖,帶你破解...
    沈念sama閱讀 221,576評論 6 515
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件哮缺,死亡現(xiàn)場離奇詭異弄跌,居然都是意外死亡,警方通過查閱死者的電腦和手機(jī)尝苇,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 94,515評論 3 399
  • 文/潘曉璐 我一進(jìn)店門铛只,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人糠溜,你說我怎么就攤上這事淳玩。” “怎么了非竿?”我有些...
    開封第一講書人閱讀 168,017評論 0 360
  • 文/不壞的土叔 我叫張陵蜕着,是天一觀的道長。 經(jīng)常有香客問我红柱,道長承匣,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 59,626評論 1 296
  • 正文 為了忘掉前任锤悄,我火速辦了婚禮悄雅,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘铁蹈。我一直安慰自己宽闲,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 68,625評論 6 397
  • 文/花漫 我一把揭開白布握牧。 她就那樣靜靜地躺著容诬,像睡著了一般。 火紅的嫁衣襯著肌膚如雪沿腰。 梳的紋絲不亂的頭發(fā)上览徒,一...
    開封第一講書人閱讀 52,255評論 1 308
  • 那天,我揣著相機(jī)與錄音颂龙,去河邊找鬼习蓬。 笑死,一個胖子當(dāng)著我的面吹牛措嵌,可吹牛的內(nèi)容都是我干的躲叼。 我是一名探鬼主播,決...
    沈念sama閱讀 40,825評論 3 421
  • 文/蒼蘭香墨 我猛地睜開眼企巢,長吁一口氣:“原來是場噩夢啊……” “哼枫慷!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 39,729評論 0 276
  • 序言:老撾萬榮一對情侶失蹤或听,失蹤者是張志新(化名)和其女友劉穎探孝,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體誉裆,經(jīng)...
    沈念sama閱讀 46,271評論 1 320
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡顿颅,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,363評論 3 340
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了足丢。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片粱腻。...
    茶點(diǎn)故事閱讀 40,498評論 1 352
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖霎桅,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情讨永,我是刑警寧澤滔驶,帶...
    沈念sama閱讀 36,183評論 5 350
  • 正文 年R本政府宣布,位于F島的核電站卿闹,受9級特大地震影響揭糕,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜锻霎,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,867評論 3 333
  • 文/蒙蒙 一著角、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧旋恼,春花似錦吏口、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,338評論 0 24
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至蜀细,卻和暖如春舟铜,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背奠衔。 一陣腳步聲響...
    開封第一講書人閱讀 33,458評論 1 272
  • 我被黑心中介騙來泰國打工谆刨, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留,地道東北人归斤。 一個月前我還...
    沈念sama閱讀 48,906評論 3 376
  • 正文 我出身青樓痊夭,卻偏偏與公主長得像,于是被迫代替她去往敵國和親脏里。 傳聞我的和親對象是個殘疾皇子生兆,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,507評論 2 359

推薦閱讀更多精彩內(nèi)容