如何在R中从GAMLSS对象提取分位数制作年龄组参考表

1. 定义目标参数 #

先明确要计算的年龄节点和百分位数：

# 目标百分位数（单位：%）
target_centiles <- c(2.5, 3, 5, 10, 25, 50, 75, 90, 95, 97, 97.5)

# 生成各年龄区间的代表性年龄点（可按需调整，比如按年取值或取区间端点/中点）
age_points <- c(seq(14, 18, 1), seq(19, 30, 1), seq(31, 65, 2))

2. 计算百分位参考值 #

直接用gamlss包的centiles()函数，基于已拟合的模型计算指定年龄和百分位的参考值：

# 计算参考值
ref_table <- centiles(m1_IAT_m, xname = "age", xvalues = age_points, cent = target_centiles)

# 转换为数据框并添加年龄分组标签
ref_df <- as.data.frame(ref_table)
ref_df$age_group <- cut(ref_df$age, 
                        breaks = c(13.9, 18, 30, 65),
                        labels = c("14-18岁", "18-30岁", "30-65岁"))

3. 整理并输出表格 #

可以用基础R或dplyr整理表格，方便查看或导出：

# 用dplyr按年龄分组整理（需先安装加载dplyr）
library(dplyr)
final_table <- ref_df %>%
  group_by(age_group) %>%
  select(age, all_of(paste0("C", target_centiles))) %>%
  arrange(age)

# 查看完整结果
print(final_table, n = Inf)

# 导出为CSV文件（可选）
write.csv(final_table, "IAT_reference_values.csv", row.names = FALSE)

备选方案（手动计算参数） #

如果centiles()函数使用有问题，可手动预测分布参数后计算分位数：

# 构造新数据框
new_data <- data.frame(age = age_points)

# 预测BCPE分布的四个参数
mu <- predict(m1_IAT_m, newdata = new_data, parameter = "mu")
sigma <- predict(m1_IAT_m, newdata = new_data, parameter = "sigma")
nu <- predict(m1_IAT_m, newdata = new_data, parameter = "nu")
tau <- predict(m1_IAT_m, newdata = new_data, parameter = "tau")

# 计算各百分位对应的参考值
ref_values <- sapply(target_centiles/100, function(p) {
  qBCPE(p, mu = mu, sigma = sigma, nu = nu, tau = tau)
})

# 组合成最终数据框
ref_df_manual <- cbind(new_data, ref_values)
colnames(ref_df_manual)[-1] <- paste0("C", target_centiles)

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