1. Rossi数据集概述与背景解析
Rossi数据集是生存分析领域的一个经典数据集,最早出现在Rossi等人1980年发表的关于累犯研究的论文中。这个数据集记录了432名罪犯在释放后的生存时间(即从释放到再次被捕的时间间隔),以及各种可能影响累犯率的协变量信息。该数据集长期以来被用作生存分析教学的示范案例,特别是在Cox比例风险模型的应用方面。
数据集包含以下关键变量:
- week:释放后的周数(时间变量)
- arrest:是否再次被捕(事件指示变量,1=是,0=否)
- fin:是否获得经济援助(1=是,0=否)
- age:释放时的年龄(岁)
- race:种族(1=黑人,0=其他)
- wexp:释放前是否有工作经历(1=是,0=否)
- mar:释放时婚姻状况(1=已婚,0=未婚)
- paro:是否获得假释(1=是,0=否)
- prio:先前犯罪次数
重要提示:在使用Rossi数据集时需要注意,虽然变量编码中的1/0值在原始研究中具有特定含义,但在不同分析软件中可能需要重新编码以适应模型要求。例如在R中,因子变量的参考水平设置会影响结果解释。
2. 数据加载与预处理实战
2.1 在R环境中加载Rossi数据集
Rossi数据集被内置在多个R包中,最常用的是survival包。以下是加载数据的标准方法:
r复制# 安装并加载survival包(如果尚未安装)
if(!require(survival)) install.packages("survival")
library(survival)
# 加载Rossi数据集
data(rossi, package = "survival")
# 查看数据结构
str(rossi)
head(rossi)
对于Python用户,可以通过以下方式获取Rossi数据集:
python复制import statsmodels.api as sm
rossi = sm.datasets.get_rdataset("rossi", "survival").data
2.2 数据质量检查与清洗
在进行生存分析前,必须进行彻底的数据质量检查:
r复制# 检查缺失值
colSums(is.na(rossi))
# 检查变量分布
summary(rossi)
# 检查时间变量的合理性
quantile(rossi$week, probs = seq(0, 1, 0.1))
常见的数据问题及处理方法:
- 极端时间值:检查是否存在异常大或小的生存时间
- 事件指示符一致性:确保arrest变量与week变量逻辑一致
- 分类变量水平:检查各分类变量的水平数量是否符合预期
2.3 数据转换与特征工程
根据分析需求,可能需要对原始变量进行转换:
r复制# 将分类变量转为因子
rossi$fin <- factor(rossi$fin, levels = c(0,1), labels = c("No","Yes"))
rossi$race <- factor(rossi$race, levels = c(0,1), labels = c("Other","Black"))
# 创建年龄分组变量
rossi$age_group <- cut(rossi$age,
breaks = c(0, 20, 30, 40, Inf),
labels = c("<20", "20-30", "30-40", "40+"))
# 计算先前的犯罪率(每生存年的犯罪次数)
rossi$prio_rate <- rossi$prio / (rossi$age - 18) # 假设18岁开始犯罪
3. 生存分析核心方法实现
3.1 Kaplan-Meier生存曲线估计
Kaplan-Meier估计是非参数生存分析的基础方法:
r复制# 拟合整体生存曲线
km_fit <- survfit(Surv(week, arrest) ~ 1, data = rossi)
# 绘制生存曲线
plot(km_fit, xlab = "Weeks", ylab = "Survival Probability",
main = "Kaplan-Meier Survival Curve")
grid()
# 按经济援助分组绘制生存曲线
km_fin <- survfit(Surv(week, arrest) ~ fin, data = rossi)
plot(km_fin, col = c("red", "blue"), lty = 1:2,
xlab = "Weeks", ylab = "Survival Probability")
legend("topright", legend = levels(rossi$fin), col = c("red","blue"), lty = 1:2)
3.2 Cox比例风险模型构建
Cox模型是生存分析中最常用的半参数方法:
r复制# 构建基础Cox模型
cox_fit <- coxph(Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio,
data = rossi)
# 查看模型摘要
summary(cox_fit)
# 检查比例风险假设
cox.zph(cox_fit)
3.3 模型诊断与验证
模型诊断是确保分析可靠性的关键步骤:
r复制# 残差分析
residuals <- residuals(cox_fit, type = "martingale")
plot(predict(cox_fit), residuals, xlab = "Linear Predictor",
ylab = "Martingale Residuals")
abline(h = 0, col = "red")
# 影响点检测
dfbeta <- residuals(cox_fit, type = "dfbeta")
par(mfrow = c(2,4))
for(j in 1:ncol(dfbeta)){
plot(dfbeta[,j], ylab = names(coef(cox_fit))[j])
abline(h = 0, col = "red")
}
4. 高级分析与结果可视化
4.1 时间依赖变量分析
当比例风险假设不成立时,需要考虑时间依赖变量:
r复制# 创建时间依赖变量
rossi$fin_time <- rossi$fin * log(rossi$week + 1)
# 拟合含时间依赖项的模型
cox_td <- coxph(Surv(week, arrest) ~ fin + tt(fin) + age + race + prio,
data = rossi,
tt = function(x, t,...) x * log(t+1))
summary(cox_td)
4.2 竞争风险分析
当存在多种事件类型时,需要竞争风险分析:
r复制if(!require(cmprsk)) install.packages("cmprsk")
library(cmprsk)
# 假设我们有两种事件类型:1=因犯罪被捕,2=因其他原因被捕
# 这里简化演示,实际应根据研究问题定义事件类型
rossi$arrest_type <- ifelse(rossi$arrest == 1, 1, 0) # 简化处理
# 计算累积发生率函数
cif <- cuminc(rossi$week, rossi$arrest_type)
plot(cif, col = c("red", "blue"), lwd = 2)
4.3 结果可视化技巧
专业的结果可视化能显著提升分析报告的质量:
r复制# 风险比森林图
if(!require(survminer)) install.packages("survminer")
library(survminer)
ggforest(cox_fit, data = rossi)
# 生存曲线分面图
ggsurvplot(km_fin, data = rossi, risk.table = TRUE,
conf.int = TRUE, pval = TRUE,
palette = c("#E7B800", "#2E9FDF"),
break.time.by = 50)
5. 常见问题与解决方案
5.1 比例风险假设违反处理
当cox.zph检验显示某些变量违反比例风险假设时:
- 分层分析:对违反假设的变量进行分层
r复制cox_stratified <- coxph(Surv(week, arrest) ~ fin + age + race + strata(wexp),
data = rossi)
-
添加时间交互项:如前面所示的时间依赖变量方法
-
使用参数模型或转换时间尺度
5.2 模型选择策略
对于多变量选择,建议采用以下方法:
r复制# 逐步回归(基于AIC)
full_model <- coxph(Surv(week, arrest) ~ ., data = rossi)
step_model <- step(full_model, direction = "both")
# 或使用LASSO方法
if(!require(glmnet)) install.packages("glmnet")
library(glmnet)
x <- model.matrix(~ fin + age + race + wexp + mar + paro + prio, rossi)[,-1]
y <- with(rossi, Surv(week, arrest))
cv_fit <- cv.glmnet(x, y, family = "cox")
plot(cv_fit)
coef(cv_fit, s = "lambda.min")
5.3 小样本处理技巧
当样本量较小时(如某些亚组分析):
- 使用精确方法而非近似方法
r复制coxph(Surv(week, arrest) ~ fin, data = rossi, method = "exact")
- 考虑贝叶斯Cox模型
r复制if(!require(rstanarm)) install.packages("rstanarm")
library(rstanarm)
bayes_fit <- stan_surv(Surv(week, arrest) ~ fin + age, data = rossi)
- 使用bootstrap获取更稳健的标准误
6. 实际应用案例扩展
6.1 预测模型构建与验证
构建可用于新数据预测的生存模型:
r复制# 拆分训练集和测试集
set.seed(123)
train_idx <- sample(1:nrow(rossi), 300)
train <- rossi[train_idx,]
test <- rossi[-train_idx,]
# 训练模型
train_fit <- coxph(Surv(week, arrest) ~ fin + age + prio, data = train)
# 计算测试集的预测风险评分
test$pred_risk <- predict(train_fit, newdata = test, type = "risk")
# 评估预测性能(时间依赖ROC)
if(!require(timeROC)) install.packages("timeROC")
library(timeROC)
roc <- timeROC(T = test$week, delta = test$arrest,
marker = test$pred_risk, cause = 1,
times = c(52, 104)) # 1年和2年的ROC
plot(roc, time = 52)
6.2 敏感性分析
评估关键假设对结果的影响:
r复制# 检查删失模式的随机性
# 使用逻辑回归检验删失是否与协变量相关
censored <- 1 - rossi$arrest
summary(glm(censored ~ fin + age + race + prio,
data = rossi, family = binomial))
# 不同时间尺度的影响
cox_logt <- coxph(Surv(log(week+1), arrest) ~ fin + age + prio, data = rossi)
6.3 亚组分析与交互效应
识别不同亚组中的差异化效应:
r复制# 年龄与金融援助的交互效应
cox_interaction <- coxph(Surv(week, arrest) ~ fin * age + prio + race, data = rossi)
summary(cox_interaction)
# 可视化交互效应
library(ggplot2)
rossi$age_centered <- rossi$age - mean(rossi$age)
interaction_fit <- coxph(Surv(week, arrest) ~ fin * age_centered + prio, data = rossi)
newdata <- expand.grid(fin = factor(c(0,1), labels = c("No","Yes")),
age_centered = seq(-15, 15, by = 5),
prio = mean(rossi$prio))
newdata$hr <- predict(interaction_fit, newdata = newdata, type = "risk")
ggplot(newdata, aes(x = age_centered + mean(rossi$age), y = hr, color = fin)) +
geom_line(size = 1.5) +
labs(x = "Age", y = "Predicted Hazard Ratio", color = "Financial Aid") +
theme_minimal()
在实际应用中,Rossi数据集的分析不应仅限于技术层面的模型拟合,更重要的是理解每个分析步骤背后的犯罪学意义。例如,经济援助(fin)的显著保护效应可能提示政策制定者考虑扩大对释放囚犯的经济支持计划。而先前犯罪次数(prio)的强预测力则强调了早期干预的重要性。
