如何基于尿肌酐与年龄校正尿钠变量并生成校正后列？

问题背景 #

我正在研究通过尿肌酐（PRECR24mmol）和年龄（PREALD）校正尿钠（PRENA24），以便将校正后的变量用于后续分析。

请问如何生成校正后的新变量？是将NA24除以肌酐和年龄，还是相乘？恳请指导。

我已运行如下线性模型，但不清楚如何利用其结果：

Call:
lm(formula = PRENA24 ~ PRECR24mmol * PREALD, data = c1.3)

Residuals:
     Min       1Q   Median       3Q      Max 
-228.439  -43.024   -5.215   37.790  274.414 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        66.84482   29.60684   2.258  0.02423 *  
PRECR24mmol         7.00565    2.10989   3.320  0.00094 ***
PREALD             -0.66555    0.60912  -1.093  0.27488    
PRECR24mmol:PREALD  0.06335    0.04392   1.442  0.14962    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 65.94 on 798 degrees of freedom
Multiple R-squared:  0.2963,    Adjusted R-squared:  0.2937 
F-statistic:   112 on 3 and 798 DF,  p-value: < 2.2e-16

我需要校正PRENA24值，并生成名为PRENA24.ADJ的新列。我知道以下代码不正确，但不清楚如何利用线性模型的信息。此外，实验后数据还按处理类型分组。

c1 <- c1.3 %>%
  mutate(PRENA24.ADJ = (PRENA24-66.84482+(7.00565*PRECR24mmol)+(-0.66555*PREALD)))
c2 <- c1 %>%
  mutate(NA24.ADJ = (NA24-24.59443+(10.54905*CR24mmol)+(0.58894*ALD)))

解决方案 #

1. 校正变量的核心逻辑：用模型残差代表校正后的值 #

当你用含交互项的线性模型校正变量时，校正后的变量本质是模型的残差——即原始尿钠值减去模型预测的尿钠值。这个残差代表去除PRECR24mmol、PREALD以及两者交互作用影响后的尿钠水平，适合用于后续分析。

简单的“除以肌酐/年龄”或“相乘”属于比例校正，仅适用于肌酐是单纯浓度稀释因子的场景；但你的模型加入了年龄和交互项，说明变量间存在更复杂的线性关系，残差法能更准确地控制这些混杂因素。

2. 正确生成校正变量的代码实现 #

方法一：用broom包自动提取残差（推荐）

broom包的augment()函数可直接从拟合好的模型中提取残差，无需手动计算系数：

library(dplyr)
library(broom)

# 先拟合模型
model <- lm(PRENA24 ~ PRECR24mmol * PREALD, data = c1.3)

# 生成包含残差的新数据框，将残重命名为目标列名
c1.3_adj <- augment(model, data = c1.3) %>%
  rename(PRENA24.ADJ = .resid)

方法二：手动计算残差（无需额外包）

若不想用broom，可根据模型系数手动计算预测值，再用原始值减去预测值得到残差：

c1.3_adj <- c1.3 %>%
  mutate(
    # 计算模型预测值：包含截距、主效应和交互项
    pred_PRENA24 = 66.84482 + (7.00565 * PRECR24mmol) + (-0.66555 * PREALD) + (0.06335 * PRECR24mmol * PREALD),
    # 校正后的值 = 原始值 - 预测值
    PRENA24.ADJ = PRENA24 - pred_PRENA24
  )

注：你之前的错误代码漏掉了交互项的计算，这是核心问题。

3. 按处理类型分组校正的方法 #

若需按处理类型分组分别校正（每组单独拟合模型），可结合group_by()和do()实现：

# 假设处理类型列名为treatment
c1.3_group_adj <- c1.3 %>%
  group_by(treatment) %>%
  do(
    augment(lm(PRENA24 ~ PRECR24mmol * PREALD, data = .), data = .) %>%
      rename(PRENA24.ADJ = .resid)
  ) %>%
  ungroup()

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