如何在R中基于系数与协方差矩阵构建预测模型并使用predict()

1. 手动定义模型参数 #

先把论文给出的系数和协方差矩阵转换成R可识别的对象：

# 定义模型系数
coefs <- c(
  intercept = 0.7889,
  age = 0.1107,
  Ln_TKV = 0.8207,
  Female = -0.0486,
  Age_Ln_TKV = -0.0160
)

# 定义协方差矩阵
cov_mat <- matrix(
  c(
    1.279758,  -0.031790, -0.175654, -0.001306,  0.004362,
   -0.031790,   0.008230,  0.004361, -0.000016, -0.000113,
   -0.175651,   0.004361,  0.024207, -0.000155, -0.000601,
   -0.001306,  -0.000016,  0.000155,  0.000708,  0.000002,
    0.004362,  -0.000113, -0.000601,  0.000002,  0.000016
  ),
  nrow = 5,
  dimnames = list(names(coefs), names(coefs))
)

2. 构建支持predict()的自定义模型对象 #

创建包含参数和公式的模型对象，并编写对应的predict方法，同时支持点估计和考虑模型变异的置信区间：

# 创建自定义模型对象
tkv_model <- list(
  coefficients = coefs,
  vcov = cov_mat,
  formula = ~ 1 + age + Ln_TKV + Female + age:Ln_TKV
)
class(tkv_model) <- "tkv_custom"

# 自定义predict方法
predict.tkv_custom <- function(object, newdata, interval = c("none", "confidence"), n_sim = 1000) {
  interval <- match.arg(interval)
  
  # 计算线性预测项
  X <- model.matrix(object$formula, newdata = newdata)
  linear_pred <- X %*% object$coefficients
  
  # 计算点预测值
  delta_tkv_pred <- exp(linear_pred) - 500
  
  if (interval == "none") {
    return(as.vector(delta_tkv_pred))
  } else if (interval == "confidence") {
    # 从多元正态分布抽样系数，模拟模型变异
    coef_samples <- MASS::mvrnorm(n_sim, mu = object$coefficients, Sigma = object$vcov)
    pred_samples <- exp(X %*% t(coef_samples)) - 500
    
    # 计算95%置信区间
    pred_ci <- t(apply(pred_samples, 1, quantile, probs = c(0.025, 0.975)))
    colnames(pred_ci) <- c("lower", "upper")
    
    return(data.frame(
      fit = as.vector(delta_tkv_pred),
      lower = pred_ci[, "lower"],
      upper = pred_ci[, "upper"]
    ))
  }
}

3. 测试模型预测功能 #

用模拟的患者数据验证模型：

# 模拟新患者数据集
new_patients <- data.frame(
  age = c(40, 55, 30),
  Ln_TKV = c(8.5, 9.2, 7.8),
  Female = c(1, 0, 1) # 1代表女性，0代表男性
)

# 获取点预测结果
predict(tkv_model, newdata = new_patients)

# 获取带95%置信区间的预测结果（包含模型变异）
predict(tkv_model, newdata = new_patients, interval = "confidence")

关键提示 #

• 自定义方法会自动处理公式中的交互项age:Ln_TKV，无需手动计算
• 生成置信区间依赖MASS包的mvrnorm函数，未安装时先运行install.packages("MASS")
• 可通过调整n_sim参数改变抽样次数，平衡精度和计算速度

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