如何用gamlss::centiles.pred估计分位数的标准误？

我用以下代码拟合了LMS模型：

m0 <- lms(weight,age, data=df)

之后，用该模型对指定年龄预测分位数：

newages <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
centiles.pred(m0, xname="age", xvalues= newages, cent= c(5,25,50,75,95))

能否对预测得到的分位数估计误差？具体方法是什么？

以下是gamlss包内置数据的可复现代码示例：

library("gamlss")
library("gamlss.data")
data("dbbmi")
m0 <- lms(bmi,age, data=dbbmi, trans.x=TRUE, k=2)
newages <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
percentiles <- centiles.pred(m0, xname="age", xvalues= newages, cent= c(5,25,50,75,95))

分位数预测误差的估计方法 #

方法1：利用predict()函数直接获取标准误

gamlss的predict()函数支持直接预测分位数并返回标准误，基于模型参数的渐近正态性计算：

cent <- c(5,25,50,75,95)
# 预测分位数及标准误
pred_result <- predict(m0, newdata = data.frame(age = newages), 
                       type = "quantile", cent = cent, se.fit = TRUE)

# 提取分位数预测值、标准误，计算95%置信区间
centile_vals <- pred_result$fit
centile_se <- pred_result$se.fit
centile_ci <- cbind(
  lower = centile_vals - 1.96 * centile_se,
  upper = centile_vals + 1.96 * centile_se
)

输出的centile_ci就是每个分位数预测值的95%置信区间，标准误centile_se可直接作为误差估计。

方法2：Bootstrap抽样估计误差

通过对原始数据重复抽样、重新拟合模型并预测分位数，用分位数的bootstrap分布来估计误差，更适合小样本或模型假设不严格的场景：

library(boot)
cent <- c(5,25,50,75,95)

# 定义bootstrap统计量函数
boot_func <- function(data, idx) {
  boot_data <- data[idx, ]
  # 重新拟合LMS模型
  boot_model <- lms(bmi, age, data = boot_data, trans.x = TRUE, k = 2)
  # 预测分位数并提取数值
  pred <- centiles.pred(boot_model, xname = "age", xvalues = newages, cent = cent)
  return(as.vector(pred$centiles))
}

# 执行bootstrap抽样（示例用100次，实际可增加到500-1000次）
set.seed(123)
boot_output <- boot(data = dbbmi, statistic = boot_func, R = 100)

# 计算单个分位数的置信区间（以age=1的5分位数为例）
boot.ci(boot_output, index = 1, type = "perc")

# 批量处理所有分位数的置信区间
ci_list <- lapply(1:length(boot_output$t0), function(i) {
  boot.ci(boot_output, index = i, type = "perc")$percent[4:5]
})
centile_boot_ci <- do.call(rbind, ci_list)

bootstrap方法得到的置信区间更稳健，不受渐近正态性假设限制。

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