題目

使用R自帶數(shù)據(jù)集veteran進(jìn)行操作，是被診斷有肺癌的老兵生存數(shù)據(jù)。

載入數(shù)據(jù)集

library（survival） library（survminer) data()#列出已載入的包的所有數(shù)據(jù)集 data(package=.packages(all.available=TRUE))#列出已安裝的包的所有數(shù)據(jù)集

可以從彈出的窗口中看到veteran數(shù)據(jù)集已經(jīng)在里面
注：可以使用?veteran， 查看數(shù)據(jù)集解釋（包括變量解釋）

str(veteran)#顯示veteran的數(shù)據(jù)基本信息 'data.frame': 137 obs. of 8 variables:$ trt : num 1 1 1 1 1 1 1 1 1 1 ...$ celltype: Factor w/ 4 levels "squamous","smallcell",..: 1 1 1 1 1 1 1 1 1 1 ...$ time : num 72 411 228 126 118 10 82 110 314 100 ...$ status : num 1 1 1 1 1 1 1 1 1 0 ...$ karno : num 60 70 60 60 70 20 40 80 50 70 ...$ diagtime: num 7 5 3 9 11 5 10 29 18 6 ...$ age : num 69 64 38 63 65 49 69 68 43 70 ...$ prior : num 0 10 0 10 10 0 10 0 0 0 ...

trt是treatment type；
cell type是癌細(xì)胞類型，包含squamous，adeno，large，small；
time是survival或censoring time（days）；
status是indicator（1指death，0指censoring）；
karno為performance status（scale of 0-100 which is karnofsky performance score）；
diagtime為disease duration（months）；
age為病人年齡；
prior指是否給予了prior therapy。

Cox PHM實(shí)現(xiàn)

擬合Cox PHM

treatment為唯一協(xié)變量（binary）

phm.trt=coxph(Surv(time,status)~trt)#擬合Cox PHM d.phm.trt=coxph.detail(phm.trt)#得到模型細(xì)節(jié) summary(phm.trt)#顯示概要 Call: coxph(formula = Surv(time, status) ~ trt)n= 137, number of events= 128 coef exp(coef) se(coef) z Pr(>|z|) trt 0.01774 1.01790 0.18066 0.098 0.922exp(coef) exp(-coef) lower .95 upper .95 trt 1.018 0.9824 0.7144 1.45Concordance= 0.525 (se = 0.026 ) Likelihood ratio test= 0.01 on 1 df, p=0.9 Wald test = 0.01 on 1 df, p=0.9 Score (logrank) test = 0.01 on 1 df, p=0.9

得到擬合模型為：
$$h(t,x)=h_0(t,\alpha )exp(0.01774x)$$
注：continuous變量（如age）同理。

celltype為唯一協(xié)變量（categorical）

此處celltype在R中為因子型變量，也是分類變量，一共四個(gè)種類，擬合所得Cox PHM將僅包含其中三個(gè)變量，因?yàn)樗膫€(gè)變量之和為1（可以理解為自由度是3）。得到如下結(jié)果：

Call: coxph(formula = Surv(time, status) ~ celltype)n= 137, number of events= 128 coef exp(coef) se(coef) z Pr(>|z|) celltypesmallcell 1.0013 2.7217 0.2535 3.950 7.83e-05 *** celltypeadeno 1.1477 3.1510 0.2929 3.919 8.90e-05 *** celltypelarge 0.2301 1.2588 0.2773 0.830 0.407 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1exp(coef) exp(-coef) lower .95 upper .95 celltypesmallcell 2.722 0.3674 1.656 4.473 celltypeadeno 3.151 0.3174 1.775 5.594 celltypelarge 1.259 0.7944 0.731 2.168Concordance= 0.608 (se = 0.028 ) Likelihood ratio test= 24.85 on 3 df, p=2e-05 Wald test = 24.09 on 3 df, p=2e-05 Score (logrank) test = 25.51 on 3 df, p=1e-05

$$h(t,x)=h_0(t,\alpha )exp(1.0013x_1+1.1477x_2+0.2301x_3)$$

Tips(刪除環(huán)境中數(shù)據(jù)）

有時(shí)候模型建立的很多，因?yàn)樾枰茨男┳兞渴秋@著有效的，所以環(huán)境中的數(shù)據(jù)就很多很雜，這時(shí)候可以刪除一些不必要的數(shù)據(jù)對(duì)象。

rm(phm.age)#刪除phm.age rm(list=ls())#刪除所有對(duì)象

利用PHM畫(huà)圖

此處以celltype為covariate擬合得到PHM后畫(huà)圖
factor型的變量無(wú)法計(jì)算均值，三種方法處理

方案一

手動(dòng)計(jì)算均值，然后填入

d.phm.celltype=coxph.detail(phm.celltype) times=c(0,d.phm.celltype$time) h0=c(0,d.phm.celltype$hazard) #hazard baseline function S0=exp(-cumsum(h0)) #survival baseline function beta=phm.celltype$coefficients #系數(shù)

coxph.detail返回the estimate for an average individual（即協(xié)變量是基于樣本的平均估計(jì)出來(lái)的）而非真正的baseline hazard，故要把協(xié)變量的均值減去。

meanx=c(0.197080292,0.197080292,0.350364964) #協(xié)變量均值 x=c(0)-meanx Sx=S0^exp(t(beta)%*%x) #生存函數(shù) xlb='t'; ylb=expression(hat(S)(t)) plot(times,Sx,type='s',xlab=xlb,ylab=ylb)

方案二

也可以選擇利用as.numeric將covariate變成continuous type，但是由于系統(tǒng)自帶數(shù)據(jù)集中的變量是無(wú)法更改的，所以要將轉(zhuǎn)化后的數(shù)據(jù)賦到一個(gè)新變量身上，對(duì)新變量建立cox phm，如：

cell=as.numeric(celltype) meanx=mean(cell) x1=1-meanx x2=2-meanx x3=3-meanx #分別對(duì)應(yīng)三類的取值1,2，3，之所以只有三個(gè)是因?yàn)槿艚r(shí)四個(gè)變量都用會(huì)導(dǎo)致多重共線性 Sx1=S0^exp(t(beta[1])%*%x1) Sx2=S0^exp(t(beta[2])%*%x2) Sx3=S0^exp(t(beta[3])%*%x3) #三個(gè)生存函數(shù)，t()是轉(zhuǎn)置函數(shù) par(mar=c(3,2,2,1)) #調(diào)整圖片邊距（若有報(bào)錯(cuò)邊距過(guò)大則要調(diào)整，否則不需要） plot(times,Sx1,pch=20) #Sx1的散點(diǎn)圖 lines(times,Sx1,col=1,type='s') lines(times,Sx2,col=2,type='s') lines(times,Sx3,col=3,type='s') #lines函數(shù)加上其它幾個(gè)函數(shù)的線圖 legend("topright",1,c('squamous','smallcell','adeno'),lty=1,col=c(1,2,3),bty='n') #調(diào)整圖例，線條，顏色等

得到結(jié)果如下：

方案三

將因子變量轉(zhuǎn)化為啞變量
利用class.ind函數(shù)

library("nnet")#class.ind()函數(shù)需要包 dummycell=class.ind(celltype) #注意class.ind得到的是四個(gè)啞變量，只能使用其中三個(gè)，否則會(huì)有共線性 ...(步驟一樣省略) meanx=colMeans(dummycell)#對(duì)每列求均值 meanx=meanx[1:length(meanx)-1]#去掉最后一個(gè)均值 ...（省略）

總結(jié)

以上是生活随笔為你收集整理的生存分析 R语言（二）——Cox PHM(1)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。