如何基于含多个解释变量的gamlss模型计算z分数与百分位数？

由于centiles.pred仅支持单解释变量的gamlss模型，针对多变量模型，我们可以通过手动获取分布参数并计算的方式实现需求，具体步骤如下：

1. 拟合多变量gamlss模型并获取分布参数 #

首先拟合目标模型，然后利用predict()函数提取每个观测（或新数据）的分布参数（如BCPE家族的mu、sigma、nu、tau）：

library(gamlss)
# 拟合含多个解释变量的模型
FIT <- gamlss(mpg ~ disp + qsec, data = mtcars, family = BCPE)

# 准备新数据
NEWDATA <- data.frame(disp = 300, qsec = 18, mpg = 17)

# 获取新数据的分布参数
new_params <- predict(FIT, newdata = NEWDATA, type = "parameters")

2. 计算z分数（对应centiles.pred的type = "z-scores"） #

z分数的计算逻辑是：先通过模型分布的累积分布函数(CDF)得到观测值y对应的累积概率p，再将p转换为标准正态分布的分位数：

# 调用BCPE分布的CDF函数计算累积概率
p <- pBCPE(NEWDATA$mpg, 
           mu = new_params[,"mu"], 
           sigma = new_params[,"sigma"], 
           nu = new_params[,"nu"], 
           tau = new_params[,"tau"])

# 转换为z分数
z_score <- qnorm(p)
z_score

3. 计算百分位数 #

3.1 对应centiles.pred的family = "centiles" #

直接使用模型分布的分位数函数，计算指定百分位的数值：

# 计算第5、50、95百分位
centiles <- qBCPE(c(0.05, 0.5, 0.95), 
                  mu = new_params[,"mu"], 
                  sigma = new_params[,"sigma"],
                  nu = new_params[,"nu"], 
                  tau = new_params[,"tau"])
names(centiles) <- c("5th", "50th", "95th")
centiles

3.2 对应centiles.pred的family = "standard-centiles" #

这种类型的百分位是基于标准正态分布的转换，直接将累积概率p乘以100即可得到：

standard_centiles <- p * 100
standard_centiles

4. 封装为通用函数 #

为了方便重复使用，可以封装一个类似centiles.pred的函数，自动适配不同的gamlss分布族：

multi_centiles_pred <- function(model, newdata, yval = NULL, centiles = c(0.05, 0.5, 0.95), type = c("z-scores", "centiles", "standard-centiles")) {
  type <- match.arg(type)
  # 获取模型对应的分布族名称
  family_name <- family(model)$family
  # 匹配对应的CDF和分位数函数
  cdf_fun <- get(paste0("p", family_name))
  quant_fun <- get(paste0("q", family_name))
  
  # 获取新数据的分布参数
  params <- predict(model, newdata = newdata, type = "parameters")
  
  if (type == "z-scores" || type == "standard-centiles") {
    if (is.null(yval)) stop("需要提供yval参数来计算z分数或标准百分位")
    # 计算累积概率
    p <- do.call(cdf_fun, c(list(q = yval), as.list(params)))
    if (type == "z-scores") {
      res <- qnorm(p)
    } else {
      res <- p * 100
    }
  } else if (type == "centiles") {
    # 计算指定百分位
    res <- do.call(quant_fun, c(list(p = centiles), as.list(params)))
    names(res) <- paste0(centiles * 100, "th")
  }
  return(res)
}

# 函数使用示例
multi_centiles_pred(FIT, newdata = NEWDATA, yval = NEWDATA$mpg, type = "z-scores")
multi_centiles_pred(FIT, newdata = NEWDATA, type = "centiles")
multi_centiles_pred(FIT, newdata = NEWDATA, yval = NEWDATA$mpg, type = "standard-centiles")

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