一、背景
数据集展示了X市常住外来人口的基本情况，包括人口学变量和一些行为特征。假定这些变量的取值在观测期间内都保持不变，仔细查看和分析数据情况，试利用生存分析法完成下面的题目。

二、要求和代码

#*********************************前期数据处理***********************************
#1
#①利用R读取数据。注意：不要事先改动样本的数据内容。
data8 <- read.csv(file="F:/hxpRlanguage/homework8.csv",header=TRUE,sep=",",stringsAsFactors = F,na.strings=c("","NA"))
#②显示数据的前10条记录
data8[1:10,]
#③对变量重新命名，一律用英文字母命名变量
cnames <- c("number","searchtime","usetime","IP","source","detailsource",
"bornyear","bornmonth","gender","spouse","children","education","income",
"status","inyear","inmonth","outyear","outmonth","member")
colnames(data8)<-cnames
#④显示重命名数据集中变量的属性情况
str(data8)
#有6个需要分类的变量
#⑤对所需变量进行分类处理,包括性别、婚姻状况、孩子状况、教育程度、现居地和居住成员
#对性别进行分类,女0男1
data8$Gender[data8$gender=="B. 女" ] <- 0
data8$Gender[data8$gender=="A. 男" ] <- 1
#对婚姻状况进行分类,有配偶1没配偶0
data8$Spouse[substring(data8$spouse,0,1)=='A'] <- 0
data8$Spouse[substring(data8$spouse,0,1)=='B'] <- 1
#对孩子状况进行分类,无孩子0其他分别对应有的123
data8$Children[substring(data8$children,0,1)=='A'] <- 0
data8$Children[substring(data8$children,0,1)=='B'] <- 1
data8$Children[substring(data8$children,0,1)=='C'] <- 2
data8$Children[substring(data8$children,0,1)=='D'] <- 3
#对教育程度进行分类,小学中学高中中专=1,大专本科=2,硕士博士=3
data8$Education[data8$education=="A. 小学" | data8$education=="B. 初中" | data8$education=="C. 高中" | data8$education=="D. 中专"] <- "1"
data8$Education[data8$education=="E. 大专" | data8$education=="F. 大学本科"] <- "2"
data8$Education[data8$education=="G. 硕士" | data8$education=="H. 博士"] <- "3"
#对现居地/迁移状况进行分类,在北京=0,不在北京=1
data8$Status[data8$status=="A. 北京" ] <- 0
data8$Status[data8$status=="B. 北京以外的地区" ] <- 1
#对居住成员状况进行分类,ABCD分别对应1234
data8$Member[substring(data8$member,0,1)=='A'] <- 1
data8$Member[substring(data8$member,0,1)=='B'] <- 2
data8$Member[substring(data8$member,0,1)=='C'] <- 3
data8$Member[substring(data8$member,0,1)=='D'] <- 4
#⑥对分类变量转换为因子类型变量
data8$fGender <- as.factor(data8$Gender)
data8$fSpouse <- as.factor(data8$Spouse)
data8$fChildren <- as.factor(data8$Children)
data8$fEducation <- as.factor(data8$Education)
data8$fStatus <- as.factor(data8$Status)
data8$fMember <- as.factor(data8$Member)
data8$fIncomelevel <- factor(data8$income)
#⑦将收入转换成对数
data8$logIncome <- log(data8$income)
#⑧查看转成因子变量后的数据集结构
str(data8)

#2
#遵循科学合理的原则，依据自己的需求对数据进行处理。
#注意：常住外来人口是指居住时间至少为6个月的人口。
#处理调查日期,即填调查问卷的日期,提取年月方便后面的计算
data8$searchtime <- as.character(data8$searchtime)
data8$survey.year <- rep(0,nrow(data8))
data8$survey.month <- rep(0,nrow(data8))
for (i in 1:nrow(data8))
{
a <- data8$searchtime[i] #将调查日期赋给a
data8$survey.year[i] <- substr(a,1,4) #取调查日期取值的前4个字符，赋给survey.year
data8$survey.month[i] <- substr(a,6,6) #取调查日期取值的第6个字符，赋给survey.month
}
data8$survey.year <- as.integer(data8$survey.year)
data8$survey.month <- as.integer(data8$survey.month)
#View(data8)
#利用R编写程序计算个体的年龄,以月为单位,仍在北京居住的个体用调查年月减出生年月,迁离北京的个体用调查年月减出生年月
#仍在北京（=0）的个体年龄
a0=(data8$survey.year*12+data8$survey.month)-(data8$bornyear*12+data8$bornmonth)
#迁离北京（=1）的个体年龄
a1=(data8$outyear*12+data8$outmonth)-(data8$bornyear*12+data8$bornmonth)
data8$Age <- ifelse(data8$outyear==1,a0,a1)
data8$fAgelevel <- factor(data8$Age)
#删除年龄小于16岁（192个月）或大于65岁（780个月）的记录
data8 <- data8[-c(which(data8$Age<192|data8$Age>780)),]
#删除收入小于500元的记录
data8 <- data8[-c(which(data8$income<500)),]
#计算北京的居住时间,d1表示还在北京的,d2表示迁离北京了的
d1=(data8$survey.year*12+data8$survey.month)-(data8$inyear*12+data8$inmonth)
d2=(data8$outyear*12+data8$outmonth)-(data8$inyear*12+data8$inmonth)
data8$Duration <- ifelse(data8$outyear==1,d1,d2)
data8 <- data8[data8$Duration>0,]

#3
#需要写明变量的定义和操作化方法，必要时展示变量的编码情况。
#设置完剩余242条数据

#4
#分别用无参数法、半参数法和参数法分析相关因素如何影响外来人口的迁留行为，并对结果进行解释。
#载入程序包survival
library(survival)
#根据Duration取值进行排序
data8 <- data8[order(data8$Duration),]
#创建生存数据对象y并查看
y <- Surv(data8$Duration,data8$fStatus==1)
#显示右删减数据情况
y 
#生成存在事件发生的时间点的集合E
y.curves <- survfit(y~1) #创建survfit对象
y.curves #显示中位数生存时间
summary(y.curves) #显示集合E
#显示生存曲线y.curves的图形
plot(y.curves,col=c(2,1,1),xlab="居住时间(单位：月)",ylab="留居概率",xlim=c(0,300))
#分析各因素不同水平的生存时间
#******************************非参数分析log-rank*****************************
#分析性别不同水平的中位生存时间
survfit(y~data8$fGender)
#画出性别的生存曲线
plot(survfit(y~data8$fGender),col=c(1,2),xlab="居住时间（单位，月）",ylab="留居概率",xlim=c(0,300))
legend(40,1,c("男性","女性"),col=c(1,2),lty=c(1,2))
#利用log-rank法检验检验因素组间差异的显著性
survdiff(y~fGender,data=data8,rho=0)

#分析教育程度不同水平的中位生存时间
survfit(y~data8$fEducation)
#画出教育程度的生存曲线
plot(survfit(y~data8$fEducation),col=c(1,2,3),lty=c(1,2,3),xlab="居住时间（单位，月）",ylab="留居概率",
xlim=c(0,300))
legend(40,1,c("小学, 初中, 高中, 中专","大专,大学本科"," 硕士, 博士"),
col=c(1,2,3),lty=c(1,2,3))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~fEducation,data = data8,rho=0)

#分析子女个数不同的中位生存时间
survfit(y~data8$fChildren)
#画出子女个数情况的生存曲线
plot(survfit(y~data8$fChildren),col=c(1,2,3,4),lty=c(1,2,3,4),xlab="居住时间（单位：月）",ylab="留居概率"
,xlim=c(0,300))
legend(40,1,c("无子女","有1个子女","有2个子女","有3个及以上子女"),col=c(1,2,3,4),lty=c(1,2,3,4))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~fChildren,data = data8,rho=0)

#分析配偶情况不同的生存时间
survfit(y~data8$fSpouse)
#画出生存曲线
plot(survfit(y~data8$fSpouse),col=c(1,2),lty=c(1,2),xlab="居住时间（单位：月）",ylab="留居概率"
,xlim=c(0,300))
legend(40,1,c("无配偶","有配偶"),col=c(1,2),lty=c(1,2))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~fSpouse,data = data8,rho=0)

#分析居住情况不同的生存时间
survfit(y~data8$fMember)
#画出生存曲线
plot(survfit(y~data8$fMember),col=c(1,2,3,4),lty=c(1,2,3,4),xlab="居住时间（单位：月）",ylab="留居概率"
,xlim=c(0,300))
legend(40,1,c("独自居住生活","住在一起生活的有主要家庭成员","住在一起生活的没有主要家庭成员但有其他亲属","其他"),col=c(1,2,3,4),lty=c(1,2,3,4))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~fMember,data = data8,rho=0)

#分析年龄不同的生存分析，若年龄是连续型变量，画出年龄的生存曲线并进行显著性检验
data8$fAge<-1*(data8$Age>median(data8$Age)) #依据中位数将年龄分为两组
plot(survfit(y~data8$fAge),col=c(1,2),xlab="居住时间（单位，月）",ylab="留居概率",xlim=c(0,300))
legend(40,1,c("年龄小的组","年龄大的组"),col=c(1,2),lty=c(1,2))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~fAge,data=data8,rho=0)

#分析收入不同的生存分析
data8$flogIncome<-1*(data8$logIncome>median(data8$logIncome)) #依据中位数将年龄分为两组
plot(survfit(y~data8$flogIncome),col=c(1,2),xlab="居住时间（单位，月）",ylab="留居概率",xlim=c(0,300))
legend(40,1,c("收入低的组","收入高的组"),col=c(1,2),lty=c(1,2))
#利用log-rank法检验因素组间差异的显著性
survdiff(y~flogIncome,data = data8,rho=0)

#********************************参数法***********************************
#利用加速死亡模型拟合数据，残差默认假设为Weibull分布
fit1 <- survreg(y~fGender+fEducation+fChildren+fSpouse+fStatus+logIncome+Age,
dist="weibull",data=data8)
summary(fit1)#scale是比例参数，对数生存时间
#多重共线性检验
library(car)
vif(fit1)
#利用AIC标准选择最优模型
fit1.aic <- step(fit1,trace = F)
summary(fit1.aic)
#利用BIC标准选择最优模型
fit1.bic <- step(fit1,trace=F,k=log(nrow(data8)))
summary(fit1.bic)
#利用加速死亡模型拟合数据，残差假设为指数分布
fit2 <- survreg(y~fGender+fEducation+fChildren+fSpouse+fStatus+logIncome+Age,
dist="exponential",data=data8)
summary(fit2)#scale是比例参数，对数生存时间
#多重共线性检验
library(car)
vif(fit2)
#利用AIC标准选择最优模型
fit2.aic <- step(fit2,trace = F)
summary(fit2.aic)
#利用BIC标准选择最优模型
fit2.bic <- step(fit2,trace=F,k=log(nrow(data8)))
summary(fit2.bic)

#*******************************半参数法************************************
#Cox比例风险模型拟合数据，对收入和年龄分组
fit3 <- coxph(y~fGender+fEducation+fChildren+fSpouse+fStatus+fIncomelevel+fAgelevel,
data=data8)
step(fit3) #逐步回归选择协变量
#根据挑选的协变量，利用Cox比例风险模型再次拟合数据，对收入和年龄分组
fit4 <- coxph(y~fEducation+fSpouse+fStatus,data=data8)
summary(fit4) #参数估计
vif(fit4) #多重共线性检验
#使用Schoenfeld残差检验fit3的PH假设
library(survminer) #载入Schoenfeld程序包
#install.packages("survminer")
library(ggplot2) #载入程序包ggplot2
#install.packages("ggplot2")
test1 <- cox.zph(fit4,transform=rank)
test1 #P值大于0.05说明符合PH假定
ggcoxzph(test1) #通过Schoenfeld残差图来判断是否符合PH假定
#利用deviance residuals图检验模型fit4中的异常值
ggcoxdiagnostics(fit4,type="deviance",liner.predictions = FALSE, ggtheme = theme_bw())
#利用dfbeta values图检验模型fit3中各协变量的异常值
ggcoxdiagnostics(fit4,type="dfbeta",liner.predictions = FALSE, ggtheme = theme_bw())
#利用AIC标准选择最优模型
fit4.aic=step(fit4,trace=F)
summary(fit4.aic)
#利用BIC标准选择最优模型
fit4.bic=step(fit4,trace=F,k=log(nrow(data8)))
summary(fit4.bic)
#Cox比例风险模型拟合数据，收入和年龄都是连续型变量
fit5 <- coxph(y~fGender+fEducation+fChildren+fSpouse+fStatus+logIncome+Age,data=data8)
summary(fit5)
#利用Martingale残差检验对收入和年龄进行非线性诊断
ggcoxfunctional(y~logIncome,data=data8)
ggcoxfunctional(y~Age,data=data8) #纵坐标是生存时间求对数然后加一个负号
#****************************************************************************