如何轻松地在 R（线图）中绘制针对混杂变量进行调整的数据？

Question

由于我现在了解了如何使用 ggpubr 包在 R 中用多个时间点绘制数据线，所以我想绘制已针对数据集中的一些已知混杂变量进行调整的数据。

我已经尝试过计算线性模型：结果〜混杂因素并绘制残差。然而，当我这样做时（由于残差的性质），值在 0 附近波动，并且我丢失了许多有关原始值的信息。 有什么办法，我可以绘制针对一个或多个混杂变量调整后的结果变量，这样我仍然可以掌握数据的原始维度（例如，我考虑从原始值中减去残差）。

我有一个关于代谢研究的数据集，其中包含葡萄糖耐量测试中血浆代谢物（葡萄糖、氨基酸等）的浓度。我想证明患有肝脂肪变性的人比没有肝脂肪变性的人具有更高的代谢物浓度。

这是我已经做的。请注意，我首先必须从数据集中排除缺失值，以便连接在调整后起作用。最后，我想绘制根据 BMI 调整的脂肪变性组（是/否）中的甘氨酸。

通过聊天 GPT 为您提供的最小数据集：D

# Creating a modified dataset
set.seed(123)

# Sample data with increasing Glycine values
main_study_long_T2D <- data.frame(
  BMI = rep(rnorm(10, mean = 25, sd = 5), 5),
  Glycine = rep(seq(50, 80, length.out = 10), 5) + ifelse(rep(c("yes", "no"), each = 5) == "yes", 10, 0) + rnorm(50, sd = 5),
  steatosis = rep(c("yes", "no"), each = 25),
  MMT_TIMEPOINT = rep(c(0, 30, 60, 120, 180), times = 10)
)

## Linear model to adjust Glycine for BMI
model_glycine <- lm(Glycine ~ BMI, data = main_study_long_T2D)
 
## Check
sum(!is.na(main_study_long_T2D$BMI))
sum(!is.na(main_study_long_T2D$Glycine)) ## Some patients / timepoints are missing in my original data
 
which(!is.na(main_study_long_T2D$Glycine))
 
## Calculate residuals
residuals_glycine <- residuals(model_glycine)
residuals_glycine
dim(residuals_glycine)
dim(as.data.frame(residuals_glycine))
 
## Patients / timepoints without missing glycine
identifiers <- main_study_long_T2D[!is.na(main_study_long_T2D$Glycine), c("PATNO", "MMT_TIMEPOINT")]
dim(identifiers)
length(residuals_glycine)
 
## Combine
residuals_data <- cbind(identifiers, residuals_glycine)
residuals_data
dim(residuals_data)
 
## Extend the dataset with residuals
main_study_long_T2D <- left_join(main_study_long_T2D, residuals_data,  by = c("PATNO", "MMT_TIMEPOINT"))
 
## ggline plot of residuals
main_study_long_T2D %>%
  filter(!(MMT_TIMEPOINT %in% c(-10, 90)) &
           DMTYP %in% c(0,2)) %>%
  mutate(DMTYP = as.factor(DMTYP)) %>%
  ggline(x="MMT_TIMEPOINT",
         y="residuals_glycine",  ## Use the calculated residuals
         color="steatosis",
         add=c("mean_se", "jitter"),
         add.params = list(alpha = 0.3),
         palette="jco") +
  stat_compare_means(aes(group = steatosis), label = "p.signif", method="t.test")

Answer 1

我们可以通过首先回归甘氨酸对 BMI 和脂肪变性的影响来消除 BMI 的影响。

library(tidyverse)
library(ggpubr)

mod <- lm(Glycine ~ BMI + steatosis, data = main_study_long_T2D)

这为我们提供了 BMI 对甘氨酸影响的系数。然后，我们选择要选择的 BMI 作为绘图目的的参考值。通常，这是研究中 BMI 值的平均值，但如果您愿意，可以将其设置为 25。BMI 与参考值之间的差值乘以系数就是我们添加到甘氨酸水平的量估计参考值的校正甘氨酸水平：

plot_df <- main_study_long_T2D %>%
  mutate(dif = coef(mod)[2] * (mean(BMI) - BMI),
         Glycine_corrected = Glycine + dif)

现在我们可以使用

ggline

查看校正后的甘氨酸水平图。我们将其与未校正的值进行比较。正如您所看到的，这些图相似但不相同：

p1 <- ggline(plot_df,
       x = "MMT_TIMEPOINT",
       y = "Glycine",
       color = "steatosis",
       add = c("mean_se", "jitter"),
       add.params = list(alpha = 0.3),
       palette = "jco") +
  stat_compare_means(aes(group = steatosis), 
                     label = "p.signif", method = "t.test") +
  ylim(50, 90) +
  ggtitle("Uncorrected")

p2 <- ggline(plot_df,
             x = "MMT_TIMEPOINT",
             y = "Glycine_corrected",
             color = "steatosis",
             add = c("mean_se", "jitter"),
             add.params = list(alpha = 0.3),
             palette = "jco") +
  stat_compare_means(aes(group = steatosis), 
                     label = "p.signif", method = "t.test") +
  ylim(50, 90) +
  ggtitle("Corrected for BMI")

library(patchwork)

p1 + p2