在R中优化分位数回归结合分段回归分析及解决AAPC计算问题

问题背景 #

基于重复横断面数据（单条观测/参与者），结合分段回归与分位数回归估计体能指标（600米跑百分位数ri_t600）随年份leto_meritev的变化，控制体重状态weight_status2，核心目标包括计算年度百分比变化（APC）、平均年度百分比变化（AAPC），并选择最优断点。

可复现数据 #

boys_one <- structure(list(
  id = c(6127454L, 3980838L, 2774166L, 1073113L, 6773578L, 5519253L, 6838723L, 6016016L, 3215266L, 1963399L, 182440L, 4669807L, 4145886L, 7097019L, 4444052L, 4997667L, 935319L, 4996953L, 1277313L, 4071109L, 5277830L, 1128991L, 2985130L, 5289288L, 3601256L, 3588296L, 749642L, 3364739L, 6415968L, 6418673L, 1051895L, 884990L, 4820890L, 4005036L, 1275617L, 5458071L, 229093L, 6250384L, 5501375L, 2355478L, 5346473L, 6913525L, 2940600L, 2258047L, 2384073L, 881742L, 5758330L, 4198975L, 650372L, 1923040L),
  ri_t600 = c(0.179709154, 0.978963377, 0.390405133, 0.02564258, 0.238324573, 0.994369221, 0.538462935, 0.128138805, 0.036528587, 0.340610927, 0.963138367, 0.797643121, 0.180022262, 0.803258459, 0.11975408, 0.001910298, 0.987400068, 0.44888971, 0.239764074, 0.160539785, 0.991878782, 0.32160318, 0.4296986, 0.108017596, 0.781181377, 0.861045956, 0.397665975, 0.100899654, 0.088997989, 0.096091582, 0.271999709, 0.195473708, 0.645156168, 0.797750286, 0.100194241, 0.51141675, 0.13893534, 0.618567227, 0.053032507, 0.317074434, 0.194654348, 0.237211849, 0.606357737, 0.521430525, 0.550012191, 0.268146643, 0.056041775, 0.705752766, 0.322134155, 0.877469326),
  leto_meritev = c(2016L, 2001L, 2007L, 2004L, 1997L, 1993L, 1997L, 1994L, 1991L, 2013L, 2017L, 1994L, 2001L, 2019L, 1995L, 2016L, 1993L, 1989L, 2016L, 2007L, 1996L, 2018L, 2009L, 1997L, 2009L, 2008L, 1993L, 1990L, 2018L, 1990L, 1990L, 2015L, 2005L, 2004L, 2018L, 2011L, 2018L, 2018L, 1995L, 2007L, 2016L, 2015L, 1993L, 2015L, 2007L, 1990L, 1996L, 2013L, 1999L, 1990L),
  weight_status2 = c("normal", "overweight", "overweight", "normal", "normal", "overweight", "normal", "normal", "normal", "normal", "overweight", "overweight", "normal", "normal", "normal", "normal", "normal", "overweight", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "overweight", "overweight", "overweight", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal", "normal")
), class = "data.frame", row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50"))

代码优化与运行效率提升 #

1. 修复segmented对象类型问题 #

当前结果非segmented类，核心原因是segmented.default对rq对象的兼容性不足，需显式指定模型类型与自定义协方差函数：

library(segmented)
library(quantreg)

# 自定义分位数回归协方差函数（避免覆盖全局方法）
vcov_rq <- function(object, ...) {
  sm <- summary(object, se = "ker", covariance = TRUE)$cov
  rownames(sm) <- colnames(sm) <- names(coef(object))
  sm
}

# 先拟合基础分位数回归
base_rq <- rq(ri_t600 ~ leto_meritev + weight_status2, tau = 0.50, data = boys_one)

# 使用segmented()主函数，显式传递vcov函数
seg_fit <- segmented(base_rq, seg.Z = ~ leto_meritev, npsi = 1, 
                     control = seg.control(display = FALSE, vcov = vcov_rq))

2. 大数据集效率优化 #

针对74万条观测的数据集，可通过以下方式加速运算：

• 使用更快的分位数回归方法：

  base_rq <- rq(ri_t600 ~ leto_meritev + weight_status2, tau = 0.50, 
                data = boys_one, method = "fn") # 内点法，比默认"br"更快
  

• 中心化年份变量，减少数值计算误差与收敛时间：

  boys_one$leto_centered <- boys_one$leto_meritev - mean(boys_one$leto_meritev)
  base_rq <- rq(ri_t600 ~ leto_centered + weight_status2, tau = 0.50, 
                data = boys_one, method = "fn")
  seg_fit <- segmented(base_rq, seg.Z = ~ leto_centered, npsi = 1, 
                       control = seg.control(display = FALSE, vcov = vcov_rq))
  

• 并行计算（需加载doParallel）：

  library(doParallel)
  cl <- makeCluster(4) # 根据CPU核心数调整
  registerDoParallel(cl)
  base_rq <- rq(ri_t600 ~ leto_centered + weight_status2, tau = 0.50, 
                data = boys_one, method = "fn", parallel = TRUE)
  stopCluster(cl)
  

替代实现方法 #

1. 混合分位数回归（flexmix包） #

通过混合模型模拟分段效果，适合多断点场景：

library(flexmix)
# 拟合2段混合分位数回归
mix_fit <- flexmix(ri_t600 ~ leto_meritev + weight_status2, 
                   data = boys_one, 
                   model = FLXMRquantreg(tau = 0.5),
                   k = 2)
# 提取断点（基于分组年份中位数）
cluster_years <- split(boys_one$leto_meritev, clusters(mix_fit))
breakpoint <- median(c(max(cluster_years[[1]]), min(cluster_years[[2]])))

2. 手动分段分位数回归 #

若已知或初步估计断点，可直接构建分段变量拟合模型：

# 假设断点为2005
boys_one$leto_post_2005 <- ifelse(boys_one$leto_meritev >= 2005, 
                                  boys_one$leto_meritev - 2005, 0)
# 拟合分段模型
seg_rq_manual <- rq(ri_t600 ~ leto_meritev + leto_post_2005 + weight_status2, 
                    tau = 0.50, data = boys_one, method = "fn")

AAPC计算与最优断点选择 #

1. 手动实现AAPC计算（含标准误） #

aapc()函数对分位数回归兼容性差，可基于Clegg等人的公式手动实现，同时考虑断点不确定性：

calc_aapc <- function(seg_model, data) {
  # 提取斜率与断点
  slopes <- coef(seg_model)[grep("leto_centered", names(coef(seg_model)))]
  psi_centered <- seg_model$psi[, "Est."]
  psi <- psi_centered + mean(data$leto_meritev) # 还原中心化断点
  year_min <- min(data$leto_meritev)
  year_max <- max(data$leto_meritev)
  
  # 计算时间段权重
  weights <- c(psi - year_min, year_max - psi)
  weights <- weights / sum(weights)
  
  # 计算APC与AAPC
  apcs <- (exp(slopes) - 1) * 100
  aapc <- sum(apcs * weights)
  
  # 使用delta方法计算AAPC标准误
  cov_mat <- vcov_rq(seg_model$model)
  grad <- weights * exp(slopes) * 100
  aapc_se <- sqrt(t(grad) %*% cov_mat %*% grad)
  
  return(list(AAPC = round(aapc, 2), SE = round(aapc_se, 2), 
              APCs = round(apcs, 2), Breakpoint = round(psi, 0)))
}

# 计算结果
aapc_result <- calc_aapc(seg_fit, boys_one)
print(aapc_result)

2. 基于BIC选择最优断点 #

segmented包不支持分位数回归的BIC，可遍历断点数量手动计算：

# 定义BIC计算函数（分位数回归似然近似）
bic_rq_segmented <- function(npsi, base_model, data) {
  fit <- try(segmented(base_model, seg.Z = ~ leto_centered, npsi = npsi, 
                       control = seg.control(display = FALSE, vcov = vcov_rq)),
             silent = TRUE)
  if(inherits(fit, "try-error")) return(Inf)
  # 分位数回归似然近似
  ll <- sum(log(pmax(1e-10, fit$residuals * (0.5 - (fit$residuals < 0)))))
  k <- length(coef(fit)) + npsi # 参数数：系数+断点
  bic <- -2*ll + k*log(nrow(data))
  return(bic)
}

# 遍历0-3个断点
npsi_candidates <- 0:3
bic_values <- sapply(npsi_candidates, bic_rq_segmented, 
                     base_model = base_rq, data = boys_one)
# 选择BIC最小的断点数量
best_npsi <- npsi_candidates[which.min(bic_values)]

常见问题解决 #

• 结果非segmented类：使用segmented()主函数而非segmented.default()，并显式传递自定义vcov函数。
• aapc()函数报错：放弃依赖该函数，使用上述自定义AAPC计算逻辑直接处理分位数回归结果。
• 大数据收敛慢：采用中心化变量、method="fn"分位数回归、并行计算等优化手段。

内容的提问来源于stack exchange，提问作者Antonio Martinko