生存分析 R语言（四）——Cox PHM(3) 检查PHA和scale

Assess the appropriateness of PHA

数据说明

使用R自带数据集ovarian。检查age是否满足proportional hazards assumption。

graphical methods

expected vs. observed

即上一篇中提到的画图方法，若两条线很接近则认为满足PHA，否则可以考虑不满足。
很明显不满足

log-log plot

age是连续性变量，所以如果要检查需要对其分组，一般分两组，以均值为界限

attach(ovarian)
str(ovarian)
?ovarian
k=survfit(Surv(futime,fustat)~(age<56))
k$label=c(rep(1,k$strata[1]),rep(2,k$strata[2]))
#此处用rep()函数做循环，将生存对象的label按分组变成了1和2
t1=c(0,subset(k$time,k$label==1))
t2=c(0,subset(k$time,k$label==2))
St1=c(1,subset(k$surv,k$label==1))
St2=c(1,subset(k$surv,k$label==2))
#分别选取了对应分组的子集
plot(t1,-log(-log(St1)),col=1,type='s',xlab='t',ylab=expression(hat(S)(t)),ylim=c(-3,5))
lines(t2,-log(-log(St2)),col=2,type='s')
#plot画图，lines函数在原图基础上追加线条
legend("topright",c('age>=56','age<56'),lty=1,col=c(1,2),bty='n')

两个曲线看上去平行，似乎可以认为PHA成立
综上可见矛盾，这种方法太主观了，假设检验的方法会更加客观

Test

attach(ovarian)
S=Surv(futime,fustat)
phm.trt=coxph(S~age)
cox.zph(phm.trt,transform='rank',global=F)
#cox.zph()用来检验PHA的函数

  chisq df    p
age 0.739  1 0.39

由此可见，应该接受原假设，age满足PHA。

Check scale——AIC

数据说明

The Worcester Heart Attack Study data(notes 4)

检查哪类变量

一般选择continuous type进行scale check

如何检查

使用AIC进行检查，对不同的模型计算其AIC，小则优。

rm(list=ls())
WHAS=read.csv(choose.files())
#数据的变量名不易识别的时候可以用col指定每列变量名
attach(WHAS)
S=Surv(lenfol,fstat)
x1=age
x3=hr
x5=diasbp
x6=bmi
x2=gender
x9=sho
x10=chf
x13=mitype
phm.1.1=coxph(S~x1+x2+x3+x5+x6+x9+x10+x13)
phm.1.2=coxph(S~x1+I(x1^2)+x2+x3+x5+x6+x9+x10+x13)
phm.1.3=coxph(S~sqrt(x1)+x2+x3+x5+x6+x9+x10+x13)
phm.1.4=coxph(S~x1+I(x1^2)+I(x1^3)+x2+x3+x5+x6+x9+x10+x13)
phm.1.5=coxph(S~log(x1)+x2+x3+x5+x6+x9+x10+x13)
aic1.1=extractAIC(phm.1.1)[2]
aic1.2=extractAIC(phm.1.2)[2]
aic1.3=extractAIC(phm.1.3)[2]
aic1.4=extractAIC(phm.1.4)[2]
aic1.5=extractAIC(phm.1.5)[2]
#选sqrt(x1)
phm.2.1=coxph(S~sqrt(x1)+x2+x3+x5+x6+x9+x10+x13)
phm.2.2=coxph(S~sqrt(x1)+x3+I(x3^2)+x2+x5+x6+x9+x10+x13)
phm.2.3=coxph(S~sqrt(x1)+sqrt(x3)+x2+x5+x6+x9+x10+x13)
phm.2.4=coxph(S~sqrt(x1)+x3+I(x3^2)+I(x3^3)+x2+x5+x6+x9+x10+x13)
phm.2.5=coxph(S~sqrt(x1)+log(x3)+x2+x5+x6+x9+x10+x13)
aic2.1=extractAIC(phm.2.1)[2]
aic2.2=extractAIC(phm.2.2)[2]
aic2.3=extractAIC(phm.2.3)[2]
aic2.4=extractAIC(phm.2.4)[2]
aic2.5=extractAIC(phm.2.5)[2]
#选log(x3)
phm.3.1=coxph(S~sqrt(x1)+log(x3)+x2+x5+x6+x9+x10+x13)
phm.3.2=coxph(S~sqrt(x1)+log(x3)+x2+x5+I(x5^2)+x6+x9+x10+x13)
phm.3.3=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+x6+x9+x10+x13)
phm.3.4=coxph(S~sqrt(x1)+log(x3)+x2+x5+I(x5^2)+I(x5^3)+x6+x9+x10+x13)
phm.3.5=coxph(S~sqrt(x1)+log(x3)+x2+log(x5)+x6+x9+x10+x13)
aic3.1=extractAIC(phm.3.1)[2]
aic3.2=extractAIC(phm.3.2)[2]
aic3.3=extractAIC(phm.3.3)[2]
aic3.4=extractAIC(phm.3.4)[2]
aic3.5=extractAIC(phm.3.5)[2]
#选sqrt(x5)
phm.4.1=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+x6+x9+x10+x13)
phm.4.2=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+x6+I(x6^2)+x9+x10+x13)
phm.4.3=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+sqrt(x6)+x9+x10+x13)
phm.4.4=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+x6+I(x6^2)+I(x6^3)+x9+x10+x13)
phm.4.5=coxph(S~sqrt(x1)+log(x3)+x2+sqrt(x5)+log(x6)+x9+x10+x13)
aic4.1=extractAIC(phm.4.1)[2]
aic4.2=extractAIC(phm.4.2)[2]
aic4.3=extractAIC(phm.4.3)[2]
aic4.4=extractAIC(phm.4.4)[2]
aic4.5=extractAIC(phm.4.5)[2]
summary(phm.4.2)
#选第二个模型

每个阶段的结果是：

最后选择的模型是：
$\sqrt{x_1}+log{x}_3+x_2+\sqrt{x_5}+x_6+I(x_6^2)+x_9+x_{1}0+x_{1}3$