plot()中不同变量着色咨询及rethinking包残差图绘制技术问题

问题1：如何在plot()函数中为不同变量着色？ #

在基础R的plot()函数里给不同分组/变量着色其实很灵活，核心就是把颜色参数col和你的分组变量绑定就行，给你几个实用的方法：

• 直接映射分组变量
  如果数据里有现成的分组列（比如group），直接把它传给col参数，R会自动给不同分组分配默认颜色：

  # 示例数据
  df <- data.frame(
    x = rnorm(50),
    y = rnorm(50),
    group = rep(c("组1", "组2"), each = 25)
  )
  # 用分组变量直接着色
  plot(df$x, df$y, col = df$group, pch = 16, xlab = "X变量", ylab = "Y变量")
  

• 自定义颜色映射
  要是你想指定更贴合需求的颜色，用ifelse()或者颜色向量手动匹配就行：

  # 用ifelse指定颜色
  plot(df$x, df$y, 
       col = ifelse(df$group == "组1", "#E63946", "#457B9D"), 
       pch = 16, xlab = "X变量", ylab = "Y变量")
  
  # 或者先定义颜色字典，再通过分组索引
  group_cols <- c("组1" = "#E63946", "组2" = "#457B9D")
  plot(df$x, df$y, 
       col = group_cols[df$group], 
       pch = 16, xlab = "X变量", ylab = "Y变量")
  

• 多分组自动配色
  如果分组超过2个，用调色板函数（比如rainbow()、viridis::viridis()）会更省心：

  # 生成3个分组的示例数据
  df$group <- rep(c("组1", "组2", "组3"), length.out = 50)
  # 生成对应数量的颜色
  cols <- rainbow(length(unique(df$group)))
  plot(df$x, df$y, 
       col = cols[as.factor(df$group)], 
       pch = 16, xlab = "X变量", ylab = "Y变量")
  

问题2：基于rethinking包的map模型绘制ASD和TD儿童的残差图 #

针对你已经构建好的m3.2模型，咱们一步步来计算残差并绘制分组残差图：

第一步：计算期望残差 #

你提到的残差是“预测ADOS值减去观测ADOS值”，先提取模型系数，计算每个儿童的预测值，再得到残差：

# 提取模型的截距a和斜率bV
model_coefs <- coef(m3.2)

# 计算每个样本的预测ADOS值（mu）
ASDTD$mu_pred <- model_coefs["a"] + model_coefs["bV"] * ASDTD$VerbalIQ.s

# 计算残差（预测值 - 观测值）
ASDTD$residual <- ASDTD$mu_pred - ASDTD$ADOS

# 要是需要常规的“观测值 - 预测值”残差，改成下面这行：
# ASDTD$residual <- ASDTD$ADOS - ASDTD$mu_pred

第二步：绘制分组残差与言语IQ的图 #

假设你的数据框ASDTD里有分组列（比如叫diagnosis，取值为"ASD"和"TD"），用基础plot()就能实现分组绘图：

# 定义分组配色
diagnosis_cols <- c("ASD" = "#1D3557", "TD" = "#A8DADC")

# 初始化空图（只画坐标轴，不画点）
plot(ASDTD$VerbalIQ.s, ASDTD$residual, 
     type = "n",  # 不绘制数据点
     xlab = "言语IQ (VerbalIQ.s)", 
     ylab = "残差值",
     main = "ASD与TD儿童的残差图")

# 添加ASD组的点
points(ASDTD$VerbalIQ.s[ASDTD$diagnosis == "ASD"], 
       ASDTD$residual[ASDTD$diagnosis == "ASD"],
       col = diagnosis_cols["ASD"], 
       pch = 16,  # 实心圆点
       cex = 1.2)

# 添加TD组的点
points(ASDTD$VerbalIQ.s[ASDTD$diagnosis == "TD"], 
       ASDTD$residual[ASDTD$diagnosis == "TD"],
       col = diagnosis_cols["TD"], 
       pch = 17,  # 实心三角
       cex = 1.2)

# 添加残差为0的参考虚线
abline(h = 0, lty = 2, col = "gray60")

# 添加图例区分两组
legend("topright", 
       legend = names(diagnosis_cols),
       col = diagnosis_cols,
       pch = c(16, 17),
       bty = "n")  # 去掉图例边框

进阶：用rethinking的link()函数更规范地获取预测值 #

如果想考虑模型参数的不确定性，可以用link()函数生成多个预测样本，再取均值作为期望预测值：

# 生成mu的后验预测样本
mu_samples <- link(m3.2)
# 计算每个样本的平均预测值
ASDTD$mu_pred <- apply(mu_samples, 2, mean)
# 后续残差计算和绘图和之前完全一致

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