【统计模型】生存分析实例——加速失效模型和Cox等比例风险模型的建立

最新推荐文章于 2024-04-30 19:24:04 发布

数据人的自我救赎

最新推荐文章于 2024-04-30 19:24:04 发布

阅读量718

点赞数

分类专栏：统计模型文章标签： python 开发语言

本文链接：https://blog.csdn.net/m0_72318954/article/details/131292102

版权

统计模型专栏收录该内容

11 篇文章 0 订阅

订阅专栏

一、代码

a=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/生存分析/生存分析数据集.csv",header=1)
nrow(a)
#变量重命名
names(a)=c("hk","gender","age","Y","C")
#将连续型变量变为分类型变量
a$age2=floor(a$age/10)#每隔十年划分一个年龄段
a[c(1:5),]
c(dim(a),mean(a$C))

library(survival)
a$YS=Surv(a$Y,a$C)#YS是生存数据，合并Y和C数据
a$YS
View(a)
a[1:5,]
head(a$YS,100)

#做生存曲线
km.fit=survfit(a$YS~1)#计算S(t)，生存函数值
summary(km.fit)
plot(km.fit,xlab="生存时间",ylab="生存概率")#绘制生存曲线
#分户口绘制生存曲线
plot(survfit(YS~hk,data=a),col=c(1,2),lty=c(1,2),xlab="生存时间",ylab="生存概率")
legend(24,0.9,c("本地户口","异地户口"),col=c(1,2),lty=c(1,2),lwd=c(1,2),cex=0.9)
#分性别
plot(survfit(YS~gender,data=a),col=c(1,2),lty=c(1,2),xlab="生存时间",ylab="生存概率")
legend(24,0.9,c("男性","女性"),col=c(1,2),lty=c(1,2),lwd=c(1,2),cex=1)
#分年龄段
plot(survfit(YS~as.factor(age2),data=a),col=c(1,2,3),lty=c(1,2,3),xlab="生存时间",ylab="生存概率")
legend(24,0.9,c("20-30","30-40","40-50"),col=c(1,2,3),lty=c(1,2,3),lwd=c(1,2),cex=1)#图例

#建立加速失效模型,YS是因变量(模型自动对因变量取对数)，任职时长
#模型结果和线性回归模型很像，牛顿迭代算法计算
model.aft=survreg(YS~hk+gender+as.factor(age2),data=a)	
summary(model.aft)										

#建立Cox等比例风险模型，模型系数正负和加速失效模型正好相反
model.cox=coxph(YS~hk+gender+as.factor(age2),data=a)
summary(model.cox)

#模型应用
pct = 1:99/100
ptime = predict(model.aft,newdata=data.frame(hk="异地",gender="男",age2=2),
                type="quantile",p=pct)#加速失效模型预测生存概率 
matplot(ptime,1-pct,xlab="在职时长",ylab="在职概率",type="l",col="darkorchid1",lwd=3)
quantile(ptime)

coefCPH = coef(model.cox) 
meanhk  = sum(a$hk == "异地")/length(a$hk)    #异地户口占比
meangender  = sum(a$gender == "女")/length(a$gender)                    #女性占比
age3 = sum(a$age2 == 3)/length(a$age2)   #30-40岁年龄段占比
age4 = sum(a$age2 == 4)/length(a$age2) #40-50岁年龄段占比
rMean = exp(coefCPH[1]*meanhk + coefCPH[2]*meangender+coefCPH[3]*age3+coefCPH[4]*age4)
r12 = exp(coefCPH[1]*1 + coefCPH[2]*0+coefCPH[3]*0 + coefCPH[4]*0)
names(r12) = names(rMean) = NULL
r12/rMean
predict(model.cox,newdata=data.frame(hk="异地",gender="男",age2=2),type="risk")
summary(predict(model.cox,newdata=a,type="risk") )