R语言生存分析实战：用coxph包30分钟内构建风险比例模型

 # 1. R语言生存分析入门 生存分析是统计学的一个分支，专门处理与生存时间（time-to-event）相关的问题。在R语言中，生存分析是通过专门的包来实现的，如`survival`包。本章将带领读者对R语言在生存分析中的应用有一个初步的了解。 ## 1.1 生存分析简介 生存分析是分析生存时间数据的一组统计方法，它关注的是事件发生的时间以及是否发生。例如，在医学研究中，这可以是病人病情复发或死亡的时间；在工业研究中，可能是产品失效的时间。 ## 1.2 R语言的生存分析包 R语言中的`survival`包是进行生存分析的主要工具。该包提供了进行生存分析所需的各种函数，包括数据处理、模型拟合、图形绘制等。此外，还有`survminer`包等其他辅助工具，可以帮助我们更好地解释和可视化生存分析的结果。 ## 1.3 生存分析的应用场景 生存分析在多个领域都有着广泛的应用。在医学研究中，生存分析用于分析患者预后、研究不同治疗方法的效果等。在社会科学、金融领域等，生存分析也可以用于预测特定事件的发生，如失业、贷款违约等。 我们将从下一章开始，深入探讨生存数据的基础概念和处理方法。 # 2. 生存数据的基础概念和处理方法 ### 2.1 生存数据的特点与结构 生存数据是医学统计和生物统计中常见的一种数据类型，主要来源于临床试验和生存研究。它描述了从试验开始到某个事件（如疾病复发、死亡等）发生的时间长度。理解生存数据的特点和结构是进行生存分析的基础。 #### 2.1.1 生存时间与事件指标 生存时间是生存分析的核心变量，通常指从研究开始到研究结束或某个感兴趣事件发生的时间长度。事件指标是与生存时间相关联的变量，通常是一个二元变量，表示事件是否发生。例如，在医学研究中，生存时间可能表示从病人被诊断患有某种疾病到病情复发的时间，而事件指标则是病情复发的有无。 代码块展示如何在R中创建生存时间和事件指标： ```r # 安装和加载survival包 install.packages("survival") library(survival) # 假设我们有以下数据集，其中time是生存时间，status是事件指标 # status为1表示事件发生，为0表示右删失（即到研究结束事件未发生） # 创建一个生存对象 my_surv_obj <- Surv(time = c(11, 13, 17, 18, 19, 22, 24), event = c(1, 1, 1, 0, 1, 0, 1)) # 查看生存对象 print(my_surv_obj) ``` 该代码块创建了一个生存对象，其中包含了生存时间和对应的事件指标，这在后续的生存分析中将会使用。 #### 2.1.2 数据清洗与预处理 在开始生存分析之前，我们需要对数据进行清洗和预处理，以确保分析的准确性和可靠性。这包括处理缺失值、异常值、以及数据转换等。一个重要的数据预处理步骤是识别和处理删失数据。删失数据是指在研究结束时，感兴趣的事件尚未发生的数据点。 数据预处理的一个示例代码： ```r # 继续使用上述数据集 # 假设原始数据集中有缺失值和异常值，我们需要对其进行清洗 data <- data.frame(time = c(11, 13, 17, NA, 19, 22, 24), status = c(1, 1, 1, 0, 1, 0, 1)) # 删除缺失值 clean_data <- na.omit(data) # 处理异常值，例如将生存时间大于20的值设为NA clean_data$time[clean_data$time > 20] <- NA # 再次删除含有NA的数据行 final_data <- na.omit(clean_data) # 查看清洗后的数据集 print(final_data) ``` 通过这段代码，我们首先删除了含有缺失值的行，然后处理了生存时间大于20的异常值，并再次清理了含有缺失值的数据行。 ### 2.2 生存数据的可视化分析 在生存分析中，数据可视化是理解数据分布和分析结果的重要手段。通过图形我们可以直观地了解数据的分布特征，以及生存时间和事件发生之间的关系。 #### 2.2.1 Kaplan-Meier曲线 Kaplan-Meier曲线是生存分析中使用最广泛的非参数估计方法，用于估计生存函数。生存函数是在某一特定时间点，研究对象存活的概率。 使用R中的survival包，我们可以轻松绘制Kaplan-Meier曲线： ```r # 使用survival包中的survfit函数拟合生存曲线 fit <- survfit(my_surv_obj ~ 1) # 绘制Kaplan-Meier曲线 plot(fit, xlab = "Survival Time", ylab = "Survival Probability", main = "Kaplan-Meier Survival Curve") # 添加95%置信区间 conf <- confint(fit) polygon(c(conf[,1], rev(conf[,2])), c(0, rep(1, length(conf[,1])), 0), col = "lightblue", border = NA) ``` 在这段代码中，我们首先使用survfit函数对生存对象进行了拟合，并绘制了Kaplan-Meier曲线。然后添加了95%的置信区间，以提供对生存概率估计的不确定性认识。 #### 2.2.2 生存函数与风险函数 除了生存函数外，风险函数也是生存分析中重要的概念。风险函数描述了在特定时间点，生存个体在下一个时间单位内发生感兴趣事件的概率。 在R中计算和绘制风险函数： ```r # 计算风险函数 hazard_fit <- survfit(Surv(time = final_data$time, event = final_data$status) ~ 1, type = "fleming-harrington") # 绘制风险函数图 plot(hazard_fit, fun = "cumhaz", main = "Cumulative Hazard Function", xlab = "Time", ylab = "Cumulative Hazard") ``` 通过上述代码，我们使用了survfit函数并指定了类型为`fleming-harrington`来计算和绘制累积风险函数。图形显示了随着时间的变化累积风险函数的变化趋势。 ### 2.3 生存数据的探索性分析 在进行正式的生存分析之前，进行探索性分析有助于了解数据的基本特征，并可能揭示潜在的数据问题。 #### 2.3.1 描述性统计分析 描述性统计分析是分析生存数据的第一步，可以提供生存时间的中心趋势和离散程度的概括。 ```r # 使用summary函数对生存对象进行描述性统计 summary(my_surv_obj) ``` 该代码块使用summary函数提供生存对象的中位数生存时间、平均生存时间、以及生存概率等描述性统计信息。 #### 2.3.2 相关性分析与分组比较 在生存数据中，我们经常需要评估不同协变量（如年龄、性别等）和生存时间之间的相关性。此外，我们还可能对不同组别的生存时间进行比较。 ```r # 创建一个分组变量 group <- ifelse(final_data$status == 1, "Event", "Censored") # 使用生存包的survdiff函数进行分组比较 group_comparison <- survdiff(Surv(final_data$time, final_data$status) ~ group) # 输出分组比较结果 print(group_comparison) ``` 在这个代码块中，我们创建了一个表示分组的变量，并使用survdiff函数对两个分组（事