The V Lab
适用对象公共卫生、流行病学及健康科学学习者
学习时长约 120–180 分钟
先修要求描述统计、置信区间与基础 R

关于本教程的数据 所有记录均由固定随机种子模拟，不包含真实个人健康信息。数据特意包含轻微非线性、异方差和部分缺失，以便练习模型诊断。模拟关系只是教学装置，不代表真实人群效应。

如何使用本教程

建议依次完成“问题定义 → 数据检查 → 拟合 → 诊断 → 解释 → 验证 → 报告”。先阅读每段解释，再运行代码，并尝试在展开答案前完成练习。代码默认显示，也可通过页面工具折叠。

学习目标

完成本教程后，你应能够：

  • 判断线性回归是否与研究问题、结局类型和数据结构相匹配；
  • 解释简单与多元线性回归中的截距、连续变量和分类变量系数；
  • 用普通最小二乘法的几何直觉理解模型拟合；
  • 报告系数、95% 置信区间、p 值、R2R^2 与调整 R2R^2；
  • 使用有意义的变量单位、交互项和非线性项回答更准确的问题；
  • 检查残差、正态 Q–Q 图、异方差、影响点和多重共线性；
  • 仅用 base R 计算 HC3 异方差稳健标准误；
  • 区分均值响应的置信区间与个体结果的预测区间；
  • 评价训练集与测试集表现，并识别过拟合和数据泄漏；
  • 说明缺失数据、研究设计与因果解释的限制；
  • 形成一段清晰、可复现、不过度解读的结果报告。

1 从研究问题到线性模型

1.1 什么时候使用线性回归？

线性回归主要描述连续型结局的条件均值如何随一个或多个预测变量变化。例如：

  • 平均收缩压是否随年龄变化？
  • 在年龄、吸烟状态和居住区域相同时，BMI 与平均收缩压有何关联？
  • 控制其他变量后，当前吸烟者与非当前吸烟者的平均收缩压相差多少？
  • 年龄与收缩压的关联是否因吸烟状态不同而改变？

一般模型写作：

Yi=β0+β1Xi1+⋯+βpXip+εi,E(εi∣Xi)=0. Y_i = \beta_0 + \beta_1X_{i1} + \cdots + \beta_pX_{ip} + \varepsilon_i, \qquad E(\varepsilon_i\mid X_i)=0.

这里的“线性”首先指模型对未知参数 β\beta 的线性组合。预测变量本身可以经过平方、分段或其他预先说明的变换，因此线性回归并不要求所有曲线都必须是直线。

先确认结局，而不是先选择函数 若结局是二元事件、计数、比例或生存时间，普通线性回归通常不是首选。方法应匹配结局分布、估计目标和抽样设计；不能仅因为 lm() 容易运行就使用它。

1.2 预测、描述与因果是不同任务

同一个回归公式可能服务于不同目标，但评价标准不同：

任务 主要问题 重点
描述关联 在已观测样本中，条件均值如何变化？ 系数、区间、函数形式与透明报告
预测 对尚未观察的新个体，能否准确预测？ 样本外误差、校准、适用人群
因果估计 若干预改变暴露，结局会如何改变？ 研究设计、时间顺序、混杂控制与可识别性假设

回归“调整了若干协变量”并不会自动把关联变成因果效应。后文所有解释默认是条件关联，除非研究设计和因果假设另有充分依据。

1.3 普通最小二乘法的直觉

对每位参与者，模型给出拟合值 Ŷi\hat Y_i；残差为 ei=Yi−Ŷie_i=Y_i-\hat Y_i。普通最小二乘法（ordinary least squares, OLS）选择系数，使残差平方和最小：

SSE=∑i=1n(Yi−Ŷi)2=∑i=1nei2. \text{SSE}=\sum_{i=1}^{n}(Y_i-\hat Y_i)^2=\sum_{i=1}^{n}e_i^2.

平方会让正负残差不互相抵消，也会让较大的偏差受到更大惩罚。OLS 的点估计并不要求结局本身服从正态分布；正态性更直接关系到小样本下经典 t 检验和置信区间的精确性。

1.3.1 矩阵视角（可选）

若设计矩阵 XX 满列秩，OLS 解为：

𝛃̂=(XTX)−1XT𝐲. \hat{\boldsymbol\beta}=(X^TX)^{-1}X^T\mathbf{y}.

这也解释了为什么完全共线会导致系数无法唯一估计：此时 XTXX^TX 不可逆。实际分析使用数值上更稳定的分解算法；无需手工求逆来拟合模型。

2 认识模拟数据

2.1 变量与分析单位

每行代表一名模拟参与者。主要结局为收缩压（mmHg），候选预测变量包括年龄、BMI、当前吸烟状态、生理性别、居住区域和每周体力活动时间。模拟完整数据在任何结局探索之前已随机保留 25% 作为测试集；本节起的 ph_data 仅包含开发数据，保留集要到“样本外预测表现”一节才会首次汇总。

str(ph_data)
## 'data.frame':    420 obs. of  11 variables:
##  $ participant_id            : chr  "P001" "P004" "P005" "P006" ...
##  $ age                       : num  23 31 38 67 44 28 33 30 48 53 ...
##  $ sex                       : Factor w/ 2 levels "Female","Male": 1 1 2 1 1 2 1 2 2 1 ...
##  $ neighborhood              : Factor w/ 4 levels "Central","North",..: 1 3 4 1 1 4 2 1 4 4 ...
##  $ smoking_status            : Factor w/ 2 levels "Not current",..: 1 1 1 1 1 2 1 1 1 1 ...
##  $ bmi                       : num  27.8 20.4 26.5 18.8 27.6 29 18.6 26.3 27.9 16.8 ...
##  $ physical_activity_min_week: num  113 69 116 105 32 361 NA 93 121 84 ...
##  $ systolic_bp               : num  117 109 120 151 111 ...
##  $ age_c10                   : num  -2.7 -1.9 -1.2 1.7 -0.6 -2.2 -1.7 -2 -0.2 0.3 ...
##  $ bmi_c5                    : num  0.56 -0.92 0.3 -1.24 0.52 0.8 -1.28 0.26 0.58 -1.64 ...
##  $ physical_activity_30      : num  3.77 2.3 3.87 3.5 1.07 ...

变量名、角色与单位应写进分析计划，而不是留到看结果后再决定。

变量 类型 本教程中的角色或单位
systolic_bp 连续 结局，mmHg
age 连续 岁
bmi 连续 kg/m²
smoking_status 二分类 当前吸烟 vs 当前不吸烟
sex 二分类 男性 vs 女性
neighborhood 多分类 北、南、西区 vs 中心区
physical_activity_min_week 连续 分钟/周；含缺失

2.2 建模前的完整性检查

data_check <- data.frame(
  变量 = names(ph_data),
  类型 = vapply(ph_data, function(x) class(x)[1], character(1)),
  缺失数 = vapply(ph_data, function(x) sum(is.na(x)), integer(1)),
  唯一值数 = vapply(ph_data, function(x) length(unique(x[!is.na(x)])), integer(1)),
  check.names = FALSE
)
knitr::kable(data_check, caption = "变量类型、缺失与唯一值检查")
变量类型、缺失与唯一值检查
变量 类型 缺失数 唯一值数
participant_id participant_id character 0 420
age age numeric 0 67
sex sex factor 0 2
neighborhood neighborhood factor 0 4
smoking_status smoking_status factor 0 2
bmi bmi numeric 0 152
physical_activity_min_week physical_activity_min_week numeric 29 201
systolic_bp systolic_bp numeric 0 288
age_c10 age_c10 numeric 0 67
bmi_c5 bmi_c5 numeric 0 152
physical_activity_30 physical_activity_30 numeric 29 201

还应检查不可能值、重复标识符、单位错误、编码变化和纳入标准。自动生成的 summary() 有用，但不能替代变量字典和领域知识。

summary(ph_data[c(
  "age", "bmi", "physical_activity_min_week", "systolic_bp",
  "smoking_status", "sex", "neighborhood"
)])
##       age            bmi       physical_activity_min_week  systolic_bp   
##  Min.   :18.0   Min.   :16.0   Min.   :  0                Min.   : 91.7  
##  1st Qu.:35.0   1st Qu.:22.2   1st Qu.: 47                1st Qu.:117.6  
##  Median :47.0   Median :24.9   Median : 94                Median :127.0  
##  Mean   :46.6   Mean   :25.0   Mean   :110                Mean   :128.0  
##  3rd Qu.:58.0   3rd Qu.:27.7   3rd Qu.:154                3rd Qu.:136.7  
##  Max.   :85.0   Max.   :38.5   Max.   :535                Max.   :170.7  
##                                NAs    :29                                
##      smoking_status     sex       neighborhood
##  Not current:353    Female:201   Central:136  
##  Current    : 67    Male  :219   North  : 85  
##                                  South  :110  
##                                  West   : 89  
##                                               
##                                               
##

2.3 探索性图形：先看形状，再拟合直线

plot(
  ph_data$age, ph_data$systolic_bp,
  pch = 16, cex = 0.65,
  col = rgb(0, 114 / 255, 178 / 255, 0.38),
  xlab = "年龄（岁）", ylab = "收缩压（mmHg）",
  main = "先检查线性关系是否合理"
)
abline(
  lm(systolic_bp ~ age, data = ph_data),
  col = palette_ph["vermillion"], lwd = 2.5
)
lines(
  lowess(ph_data$age, ph_data$systolic_bp, f = 2 / 3),
  col = palette_ph["teal"], lwd = 2.5, lty = 2
)
legend(
  "topleft",
  legend = c("OLS 直线", "LOWESS 平滑"),
  col = c(palette_ph["vermillion"], palette_ph["teal"]),
  lwd = 2.5, lty = c(1, 2), bty = "n"
)
年龄与收缩压散点图，叠加线性拟合线与平滑曲线。

收缩压与年龄的散点图。直线表示简单线性拟合，曲线表示 LOWESS 平滑。

平滑线不是最终模型，而是一种函数形式检查。如果平滑线与直线系统性分离，应考虑领域上合理的非线性项，而不是盲目加入高阶多项式。

boxplot(
  systolic_bp ~ smoking_status,
  data = ph_data,
  names = c("当前不吸烟", "当前吸烟"),
  col = c("#D9EAF3", "#F4C6A6"),
  border = palette_ph["navy"],
  ylab = "收缩压（mmHg）", xlab = "",
  main = "分类预测变量的组间分布"
)
stripchart(
  systolic_bp ~ smoking_status,
  data = ph_data,
  vertical = TRUE, method = "jitter", pch = 16,
  col = rgb(23 / 255, 50 / 255, 77 / 255, 0.22),
  add = TRUE
)
当前吸烟者与当前不吸烟者收缩压的箱线图和散点。

按吸烟状态比较收缩压分布；箱线图之外叠加了抖动后的个体观测。

原始组间差异可能同时反映年龄、BMI、区域构成等变量的差异。多元回归可以描述在模型所含协变量相同时的条件均值差，但其可信度仍依赖模型设定和研究设计。

3 简单线性回归

3.1 一个预测变量的模型

先用年龄解释平均收缩压：

E(SBP∣age)=β0+β1age. E(\text{SBP}\mid \text{age})=\beta_0+\beta_1\text{age}. 

simple_model <- lm(systolic_bp ~ age, data = ph_data)
summary(simple_model)
## 
## Call:
## lm(formula = systolic_bp ~ age, data = ph_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -35.45  -7.84   0.70   7.41  32.12 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 101.6771     1.7177    59.2   <2e-16 ***
## age           0.5660     0.0349    16.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.3 on 418 degrees of freedom
## Multiple R-squared:  0.386,  Adjusted R-squared:  0.384 
## F-statistic:  263 on 1 and 418 DF,  p-value: <2e-16
knitr::kable(
  coefficient_table(simple_model),
  digits = 3,
  caption = "年龄与收缩压的简单线性回归"
)
年龄与收缩压的简单线性回归
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 101.677 1.718 98.301 105.054 0
年龄（每 1 岁） 0.566 0.035 0.497 0.635 0

年龄系数为每增加 1 岁 0.57 mmHg（95% CI：0.5 至 0.63）。它描述样本中年龄相差 1 岁者的平均收缩压差异，并不证明年龄变化本身造成该差异。

截距是年龄为 0 岁时的外推平均值；该年龄不在本教程的成人样本中，因此截距主要用于定位直线，通常没有实质解释。将年龄中心化可让截距对应一个有意义的参照年龄。

3.2 拟合值、残差与观测值

simple_components <- data.frame(
  participant_id = ph_data$participant_id[1:8],
  observed = ph_data$systolic_bp[1:8],
  fitted = fitted(simple_model)[1:8],
  residual = residuals(simple_model)[1:8]
)
knitr::kable(
  simple_components,
  digits = 2,
  caption = "前 8 名参与者的观测值、拟合值与残差"
)
前 8 名参与者的观测值、拟合值与残差
participant_id observed fitted residual
1 P001 117 115 2.20
4 P004 109 119 -10.52
5 P005 120 123 -3.59
6 P006 151 140 11.70
7 P007 111 127 -15.18
8 P008 138 118 20.07
9 P009 111 120 -9.46
10 P010 122 119 3.34

含截距的 OLS 模型中，残差和在数值误差范围内为 0；这不等于模型对每个人都预测准确。

c(
  residual_sum = sum(residuals(simple_model)),
  residual_mean = mean(residuals(simple_model)),
  residual_rmse = sqrt(mean(residuals(simple_model)^2))
)
##  residual_sum residual_mean residual_rmse 
##      3.55e-13      1.27e-15      1.13e+01

4 多元线性回归

4.1 同时纳入连续与分类预测变量

multiple_model <- lm(
  systolic_bp ~ age + bmi + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(multiple_model),
  digits = 3,
  caption = "收缩压的多元线性回归"
)
收缩压的多元线性回归
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 84.046 3.500 77.165 90.926 0.000
年龄（每 1 岁） 0.557 0.033 0.491 0.622 0.000
BMI（每 1 单位） 0.640 0.126 0.393 0.886 0.000
当前吸烟（参照：当前不吸烟） 4.209 1.433 1.392 7.025 0.003
男性（参照：女性） 2.083 1.058 0.003 4.162 0.050
北区（参照：中心区） 2.840 1.482 -0.074 5.754 0.056
南区（参照：中心区） -2.728 1.369 -5.419 -0.037 0.047
西区（参照：中心区） 2.262 1.460 -0.609 5.132 0.122

连续变量系数是其他模型变量保持相同时，该变量增加一个单位对应的条件均值差。例如，年龄系数比较 BMI、吸烟状态、生理性别和区域相同而年龄相差 1 岁的参与者。

这种“保持相同”是模型中的条件比较，不意味着样本中一定存在完全匹配的真实个体，也不保证比较具有因果含义。

4.2 分类变量与参照水平

R 默认把因子的第一个水平作为参照。应在拟合前明确设置，而不是根据最显著的结果改变参照组。

list(
  smoking_status = levels(ph_data$smoking_status),
  sex = levels(ph_data$sex),
  neighborhood = levels(ph_data$neighborhood)
)
## $smoking_status
## [1] "Not current" "Current"    
## 
## $sex
## [1] "Female" "Male"  
## 
## $neighborhood
## [1] "Central" "North"   "South"   "West"

因此，smoking_statusCurrent 比较当前吸烟者与当前不吸烟者；三个区域系数分别比较北、南、西区与中心区。含 kk 个水平的无序因子通常产生 k−1k-1 个系数。

若研究问题要求以北区为参照，可显式重设水平：

ph_reference_demo <- ph_data
ph_reference_demo$neighborhood <- relevel(
  ph_reference_demo$neighborhood,
  ref = "North"
)
reference_demo_model <- lm(
  systolic_bp ~ age + bmi + smoking_status + sex + neighborhood,
  data = ph_reference_demo
)
coef(reference_demo_model)[grep("neighborhood", names(coef(reference_demo_model)))]
## neighborhoodCentral   neighborhoodSouth    neighborhoodWest 
##              -2.840              -5.568              -0.578

改变参照水平会改变部分系数的表达，但不会改变每人的拟合值、残差或整体模型拟合。

4.3 用有意义的单位重标度与中心化

“每 1 岁”与“每 1 个 BMI 单位”的效应可能太细碎。本教程定义：

  • age_c10 = (age - 50) / 10：每 10 岁，且 0 对应 50 岁；
  • bmi_c5 = (bmi - 25) / 5：每 5 个 BMI 单位，且 0 对应 BMI 25。
scaled_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(scaled_model),
  digits = 3,
  caption = "以有意义单位表达的多元模型"
)
以有意义单位表达的多元模型
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 127.86 1.085 125.730 129.997 0.000
年龄（每 10 岁；中心 50 岁） 5.57 0.335 4.907 6.224 0.000
BMI（每 5 单位；中心 25） 3.20 0.628 1.963 4.432 0.000
当前吸烟（参照：当前不吸烟） 4.21 1.433 1.392 7.025 0.003
男性（参照：女性） 2.08 1.058 0.003 4.162 0.050
北区（参照：中心区） 2.84 1.482 -0.074 5.754 0.056
南区（参照：中心区） -2.73 1.369 -5.419 -0.037 0.047
西区（参照：中心区） 2.26 1.460 -0.609 5.132 0.122

调整其他模型变量后，年龄相差 10 岁对应的平均收缩压差为 5.57 mmHg。BMI 相差 5 个单位对应的平均差为 3.2 mmHg。中心化没有改变拟合优度，只使截距对应 50 岁、BMI 25、当前不吸烟、女性、中心区这一参照组合。

4.4 系数的置信区间

点估计给出最符合样本的数值，置信区间表达重复抽样不确定性。区间宽度受样本量、残差变异、预测变量分布和共线性影响。

scaled_ci <- confint(scaled_model)
keep_coef <- rownames(scaled_ci) != "(Intercept)"
coef_values <- coef(scaled_model)[keep_coef]
ci_values <- scaled_ci[keep_coef, , drop = FALSE]
plot(
  coef_values, seq_along(coef_values),
  xlim = range(ci_values), yaxt = "n",
  pch = 19, col = palette_ph["blue"],
  xlab = "平均收缩压差（mmHg）", ylab = "",
  main = "系数与 95% 置信区间"
)
segments(
  ci_values[, 1], seq_along(coef_values),
  ci_values[, 2], seq_along(coef_values),
  col = palette_ph["navy"], lwd = 2
)
abline(v = 0, lty = 2, col = "grey45")
axis(2, at = seq_along(coef_values),
     labels = label_terms(names(coef_values)), las = 1, cex.axis = 0.72)
显示各回归系数点估计和水平置信区间的图。

多元线性回归系数及 95% 置信区间；为便于比较省略截距。

不要把“区间包含 0”等同于“没有关联”，也不要把“不包含 0”等同于“具有公共卫生重要性”。效应大小、区间范围、测量单位和研究背景应共同解释。

5 模型拟合程度：R² 与调整 R²

5.1 两个指标回答什么？

R2=1−∑i(Yi−Ŷi)2∑i(Yi−Y‾)2. R^2=1-\frac{\sum_i(Y_i-\hat Y_i)^2}{\sum_i(Y_i-\bar Y)^2}.

R2R^2 表示模型在当前样本中解释的结局总变异比例。加入任何预测变量都不会使训练样本 R2R^2 下降，即使新变量没有实际价值。调整 R2R^2 对模型复杂度作惩罚，可能下降，但它仍不是样本外表现的替代品。

fit_statistics <- function(model, model_name) {
  s <- summary(model)
  data.frame(
    模型 = model_name,
    样本量 = nobs(model),
    参数数 = length(coef(model)),
    R平方 = s$r.squared,
    调整R平方 = s$adj.r.squared,
    残差标准误 = s$sigma,
    check.names = FALSE
  )
}

fit_comparison <- rbind(
  fit_statistics(simple_model, "仅年龄"),
  fit_statistics(scaled_model, "年龄 + BMI + 吸烟 + 性别 + 区域")
)
knitr::kable(fit_comparison, digits = 3, caption = "样本内拟合指标")
样本内拟合指标
模型 样本量 参数数 R平方 调整R平方 残差标准误
仅年龄 420 2 0.386 0.384 11.3
年龄 + BMI + 吸烟 + 性别 + 区域 420 8 0.461 0.452 10.6

较低的 R2R^2 不一定使一个估计良好的关联失去价值；较高的 R2R^2 也不证明模型无偏、可推广或具有因果解释。

6 交互作用：关联是否因组别而异？

6.1 建立年龄 × 吸烟状态交互

无交互的模型假定两组年龄斜率相同。交互项允许当前吸烟者的年龄斜率不同：

interaction_model <- lm(
  systolic_bp ~ age_c10 * smoking_status +
    bmi_c5 + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(interaction_model),
  digits = 3,
  caption = "含年龄与吸烟状态交互的模型"
)
含年龄与吸烟状态交互的模型
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 127.89 1.082 125.760 130.012 0.000
年龄（每 10 岁；中心 50 岁） 5.30 0.359 4.599 6.010 0.000
当前吸烟（参照：当前不吸烟） 4.30 1.429 1.491 7.107 0.003
BMI（每 5 单位；中心 25） 3.21 0.626 1.983 4.444 0.000
男性（参照：女性） 2.07 1.054 0.000 4.145 0.050
北区（参照：中心区） 2.61 1.482 -0.301 5.524 0.079
南区（参照：中心区） -2.84 1.365 -5.525 -0.157 0.038
西区（参照：中心区） 2.05 1.459 -0.823 4.914 0.162
年龄 × 当前吸烟 1.93 0.976 0.010 3.847 0.049

在含交互的模型中：

  • age_c10 是当前不吸烟者每 10 岁的斜率；
  • smoking_statusCurrent 是在 age_c10 = 0，即 50 岁时两组均值差；
  • 交互系数是当前吸烟者与当前不吸烟者年龄斜率之差；
  • 当前吸烟者的年龄斜率等于年龄主效应加交互效应。

6.2 用线性组合得到两组斜率

linear_combination_ci <- function(model, weights, level = 0.95) {
  beta <- coef(model)
  V <- vcov(model)
  L <- setNames(rep(0, length(beta)), names(beta))
  L[names(weights)] <- weights
  estimate <- sum(L * beta)
  se <- sqrt(drop(t(L) %*% V %*% L))
  critical <- qt(1 - (1 - level) / 2, df.residual(model))
  c(
    estimate = estimate,
    standard_error = se,
    lower = estimate - critical * se,
    upper = estimate + critical * se
  )
}

slope_not_current <- linear_combination_ci(
  interaction_model,
  c(age_c10 = 1)
)
slope_current <- linear_combination_ci(
  interaction_model,
  c(age_c10 = 1, "age_c10:smoking_statusCurrent" = 1)
)
interaction_slopes <- rbind(
  "当前不吸烟" = slope_not_current,
  "当前吸烟" = slope_current
)
knitr::kable(
  interaction_slopes,
  digits = 2,
  caption = "按吸烟状态计算的每 10 岁收缩压斜率（mmHg）"
)
按吸烟状态计算的每 10 岁收缩压斜率（mmHg）
estimate standard_error lower upper
当前不吸烟 5.30 0.36 4.60 6.01
当前吸烟 7.23 0.91 5.45 9.02 
age_grid <- seq(min(ph_data$age), max(ph_data$age), length.out = 120)
interaction_grid <- rbind(
  data.frame(
    age = age_grid,
    age_c10 = (age_grid - 50) / 10,
    bmi_c5 = 0,
    smoking_status = factor("Not current", levels = levels(ph_data$smoking_status)),
    sex = factor("Female", levels = levels(ph_data$sex)),
    neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
  ),
  data.frame(
    age = age_grid,
    age_c10 = (age_grid - 50) / 10,
    bmi_c5 = 0,
    smoking_status = factor("Current", levels = levels(ph_data$smoking_status)),
    sex = factor("Female", levels = levels(ph_data$sex)),
    neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
  )
)
interaction_grid$prediction <- predict(interaction_model, newdata = interaction_grid)

plot(
  systolic_bp ~ age, data = ph_data,
  pch = 16, cex = 0.45, col = "grey80",
  xlab = "年龄（岁）", ylab = "预测平均收缩压（mmHg）",
  main = "用图形解释交互作用"
)
for (group_index in seq_along(levels(ph_data$smoking_status))) {
  rows <- interaction_grid$smoking_status == levels(ph_data$smoking_status)[group_index]
  lines(
    interaction_grid$age[rows], interaction_grid$prediction[rows],
    col = c(palette_ph["blue"], palette_ph["vermillion"])[group_index],
    lwd = 3
  )
}
legend(
  "topleft", legend = c("当前不吸烟", "当前吸烟"),
  col = c(palette_ph["blue"], palette_ph["vermillion"]),
  lwd = 3, bty = "n"
)
两条按吸烟状态区分的年龄与预测收缩压回归线。

模型预测的年龄—收缩压关系，按吸烟状态分层；BMI 固定为 25、性别为女性、区域为中心区。

层级原则 若模型包含 X×ZX\times Z，通常保留 XX 与 ZZ 两个主效应，即使其中某个 p 值较大。交互应由预先定义的科学问题驱动，并在有意义的变量范围内用分层预测或边际效应解释。

7 非线性关系

7.1 用中心化二次项表示曲率

若年龄斜率并非常数，可加入中心化年龄的平方项：

linear_age_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)
quadratic_model <- lm(
  systolic_bp ~ age_c10 + I(age_c10^2) +
    bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(quadratic_model),
  digits = 3,
  caption = "含中心化年龄二次项的模型"
)
含中心化年龄二次项的模型
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 125.556 1.154 123.286 127.825 0.000
年龄（每 10 岁；中心 50 岁） 6.152 0.347 5.471 6.834 0.000
年龄二次项 0.889 0.180 0.536 1.243 0.000
BMI（每 5 单位；中心 25） 3.225 0.611 2.023 4.426 0.000
当前吸烟（参照：当前不吸烟） 4.566 1.396 1.823 7.310 0.001
男性（参照：女性） 2.630 1.035 0.595 4.665 0.011
北区（参照：中心区） 3.198 1.444 0.360 6.036 0.027
南区（参照：中心区） -3.307 1.337 -5.935 -0.679 0.014
西区（参照：中心区） 1.940 1.422 -0.855 4.735 0.173

当二次项存在时，age_c10 是 50 岁附近的局部斜率，而不是所有年龄共同的斜率。不要单独解释二次项；应结合两个年龄项绘制预测曲线。

knitr::kable(
  as.data.frame(anova(linear_age_model, quadratic_model)),
  digits = 3,
  caption = "嵌套线性年龄模型与二次年龄模型的比较"
)
嵌套线性年龄模型与二次年龄模型的比较
Res.Df RSS Df Sum of Sq F Pr(>F)
412 46714 NA NA NA NA
411 44094 1 2620 24.4 0

这个 F 检验只比较两个嵌套模型，不能证明二次曲线是真实机制。函数形式还应根据图形、预测目的、先验知识和外部验证来判断。

quadratic_grid <- data.frame(
  age_c10 = (age_grid - 50) / 10,
  bmi_c5 = 0,
  smoking_status = factor("Not current", levels = levels(ph_data$smoking_status)),
  sex = factor("Female", levels = levels(ph_data$sex)),
  neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
)
quadratic_prediction <- predict(
  quadratic_model,
  newdata = quadratic_grid,
  interval = "confidence"
)
plot(
  age_grid, quadratic_prediction[, "fit"], type = "n",
  ylim = range(quadratic_prediction),
  xlab = "年龄（岁）", ylab = "预测平均收缩压（mmHg）",
  main = "不要从单个二次项系数想象曲线"
)
polygon(
  c(age_grid, rev(age_grid)),
  c(quadratic_prediction[, "lwr"], rev(quadratic_prediction[, "upr"])),
  col = rgb(86 / 255, 180 / 255, 233 / 255, 0.28), border = NA
)
lines(age_grid, quadratic_prediction[, "fit"],
      col = palette_ph["blue"], lwd = 3)
年龄与预测平均收缩压的弯曲关系，带阴影置信区间。

调整其他变量后，年龄与预测平均收缩压的二次关系及 95% 置信带。

高阶多项式在数据边缘可能剧烈摆动，外推尤其危险。更复杂的关系可考虑分段线性、样条或领域指定的变换，但复杂度应与样本量和研究目的匹配。

8 建模假设与诊断

8.1 假设清单

经典线性回归推断通常关注：

  1. 条件均值设定合理：遗漏的非线性、交互或重要结构不会系统性留在残差中；
  2. 观测独立：重复测量、家庭、学校、社区或医院聚类需要专门处理；
  3. 同方差：给定预测变量后，误差变异近似恒定；
  4. 残差近似正态：主要影响小样本经典区间和检验，不要求每个预测变量正态；
  5. 无完全共线：任何预测变量都不是其他设计矩阵列的精确组合；
  6. 测量与抽样可信：严重测量误差、选择偏倚和不恰当权重不会由残差图自动暴露；
  7. 模型适用范围明确：预测不应无依据地外推到样本未覆盖的人群或数值范围。

诊断不是一次“通过/不通过”考试。它是发现模型与数据不一致、评估结论敏感性并改进报告的过程。

8.2 一个用于完整诊断的工作模型

工作模型同时允许年龄曲率与年龄—吸烟交互：

working_model <- lm(
  systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) +
    bmi_c5 + sex + neighborhood,
  data = ph_data
)
summary(working_model)
## 
## Call:
## lm(formula = systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) + 
##     bmi_c5 + sex + neighborhood, data = ph_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.272  -6.761   0.021   6.692  26.530 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    125.623      1.152  109.02  < 2e-16 ***
## age_c10                          5.916      0.372   15.92  < 2e-16 ***
## smoking_statusCurrent            4.636      1.393    3.33  0.00095 ***
## I(age_c10^2)                     0.871      0.180    4.84  1.8e-06 ***
## bmi_c5                           3.238      0.610    5.31  1.8e-07 ***
## sexMale                          2.610      1.033    2.53  0.01187 *  
## neighborhoodNorth                2.994      1.445    2.07  0.03889 *  
## neighborhoodSouth               -3.392      1.334   -2.54  0.01140 *  
## neighborhoodWest                 1.762      1.422    1.24  0.21618    
## age_c10:smoking_statusCurrent    1.657      0.952    1.74  0.08253 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.3 on 410 degrees of freedom
## Multiple R-squared:  0.495,  Adjusted R-squared:  0.484 
## F-statistic: 44.6 on 9 and 410 DF,  p-value: <2e-16

8.3 残差、Q–Q 图、尺度位置图与影响图

old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(working_model, which = 1:4, caption = rep("", 4))
线性模型的残差拟合值图、正态 Q-Q 图、尺度位置图和残差杠杆图。

工作模型的四个标准诊断图。应关注系统形状、尾部偏离、漏斗形和高影响观测。 

par(old_par)
  • Residuals vs Fitted：围绕 0 的无结构云团较理想；弯曲提示函数形式问题；
  • Normal Q–Q：重点看系统偏离和尾部；轻微偏离在大样本中未必严重；
  • Scale–Location：明显上升或漏斗形提示条件方差变化；
  • Cook’s distance：识别对整体拟合有较大影响的观测，不能作为自动删除规则。

8.4 异方差的定量筛查

先把拟合值分为五组，比较残差标准差：

fitted_groups <- cut(
  fitted(working_model),
  breaks = quantile(fitted(working_model), probs = seq(0, 1, 0.2)),
  include.lowest = TRUE
)
variance_by_fit <- data.frame(
  拟合值组 = levels(fitted_groups),
  样本量 = as.vector(table(fitted_groups)),
  残差标准差 = as.vector(tapply(residuals(working_model), fitted_groups, sd)),
  check.names = FALSE
)
knitr::kable(variance_by_fit, digits = 2, caption = "按拟合值五分组的残差变异")
按拟合值五分组的残差变异
拟合值组 样本量 残差标准差
[107,119] 84 10.65
(119,124] 84 9.55
(124,130] 84 9.55
(130,137] 84 9.63
(137,166] 84 11.73

下面给出只有一个辅助预测变量的 Breusch–Pagan 风格筛查。它用残差平方对拟合值回归，统计量近似服从自由度 1 的卡方分布。

hetero_auxiliary <- lm(I(residuals(working_model)^2) ~ fitted(working_model))
hetero_statistic <- nobs(working_model) * summary(hetero_auxiliary)$r.squared
hetero_p_value <- pchisq(hetero_statistic, df = 1, lower.tail = FALSE)
c(statistic = hetero_statistic, df = 1, p_value = hetero_p_value)
## statistic        df   p_value 
##     1.319     1.000     0.251

这只是简化筛查，并非最终裁决。显著结果可能来自错误函数形式；不显著也不证明完全同方差。应结合残差图、数据生成过程和稳健敏感性分析。

8.5 杠杆值、学生化残差与 Cook 距离

influence_summary <- data.frame(
  participant_id = ph_data$participant_id,
  leverage = hatvalues(working_model),
  studentized_residual = rstudent(working_model),
  cooks_distance = cooks.distance(working_model)
)
top_influence <- order(
  influence_summary$cooks_distance,
  decreasing = TRUE
)[1:8]
knitr::kable(
  influence_summary[top_influence, ],
  digits = 3,
  caption = "Cook 距离最大的 8 个观测"
)
Cook 距离最大的 8 个观测
participant_id leverage studentized_residual cooks_distance
486 P486 0.108 -1.91 0.044
399 P399 0.035 -3.01 0.032
235 P235 0.026 -2.91 0.022
543 P543 0.037 -2.35 0.021
71 P071 0.034 -2.40 0.020
375 P375 0.135 -1.10 0.019
367 P367 0.056 1.76 0.018
213 P213 0.035 -2.23 0.018

常见启发式阈值包括杠杆值 2p/n2p/n、Cook 距离 4/n4/n 和学生化残差绝对值 2 或 3，但它们不是普适的删除标准。

p_working <- ncol(model.matrix(working_model))
n_working <- nobs(working_model)
c(
  leverage_2p_over_n = 2 * p_working / n_working,
  cooks_4_over_n = 4 / n_working,
  count_abs_studentized_over_2 = sum(abs(rstudent(working_model)) > 2)
)
##           leverage_2p_over_n               cooks_4_over_n 
##                      0.04762                      0.00952 
## count_abs_studentized_over_2 
##                     18.00000

发现高影响点时，应核对原始记录、单位和纳入标准，并比较保留与排除后的结果。只有可辩护的数据质量或目标人群理由才能支持排除；“让 p 值变小”不是理由。

8.6 多重共线性

共线性不会必然造成预测偏差，但会放大个别系数的标准误，使系数对数据扰动敏感。下面以设计矩阵的每一列为单位计算 VIF：

vif_columns <- function(model) {
  X <- model.matrix(model)
  X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
  values <- vapply(seq_len(ncol(X)), function(j) {
    outcome_column <- X[, j]
    other_columns <- X[, -j, drop = FALSE]
    r_squared_j <- summary(lm(outcome_column ~ other_columns))$r.squared
    if (1 - r_squared_j < .Machine$double.eps^0.5) Inf else 1 / (1 - r_squared_j)
  }, numeric(1))
  data.frame(
    设计矩阵列 = label_terms(colnames(X)),
    VIF = values,
    row.names = NULL,
    check.names = FALSE
  )
}

knitr::kable(
  vif_columns(working_model),
  digits = 2,
  caption = "工作模型设计矩阵列层面的 VIF"
)
工作模型设计矩阵列层面的 VIF
设计矩阵列 VIF
年龄（每 10 岁；中心 50 岁） 1.35
当前吸烟（参照：当前不吸烟） 1.02
年龄二次项 1.17
BMI（每 5 单位；中心 25） 1.03
男性（参照：女性） 1.05
北区（参照：中心区） 1.33
南区（参照：中心区） 1.35
西区（参照：中心区） 1.33
年龄 × 当前吸烟 1.17

VIF 取 5 或 10 作为警戒线只是经验规则。含多水平因子、交互或多项式时，单列 VIF 尤其需要谨慎解释。中心化可降低非必要的数值共线，但不能解决两个不同变量几乎测量同一概念的问题。

9 HC3 异方差稳健标准误

9.1 为什么需要稳健标准误？

若条件均值模型合理但误差方差不恒定，OLS 系数仍可描述条件均值关联，但经典标准误可能不可靠。异方差一致协方差矩阵用每个观测的残差信息替代共同方差假设。

HC3 对高杠杆观测作较强修正，其“肉”矩阵使用：

Ω̂ii=ei2(1−hii)2,V̂HC3=(XTX)−1XTΩ̂X(XTX)−1. \widehat\Omega_{ii}=\frac{e_i^2}{(1-h_{ii})^2}, \qquad \widehat V_{HC3}=(X^TX)^{-1}X^T\widehat\Omega X(X^TX)^{-1}.

9.2 仅用 base R 实现 HC3

hc3_table <- function(model, level = 0.95) {
  if (anyNA(coef(model))) {
    stop("模型含不可估计系数；请先处理完全共线。")
  }
  X <- model.matrix(model)
  residual <- residuals(model)
  leverage <- hatvalues(model)
  omega_hc3 <- residual^2 / (1 - leverage)^2

  bread <- solve(crossprod(X))
  meat <- crossprod(X, X * omega_hc3)
  covariance_hc3 <- bread %*% meat %*% bread
  standard_error <- sqrt(diag(covariance_hc3))

  estimate <- coef(model)
  statistic <- estimate / standard_error
  degrees_freedom <- df.residual(model)
  critical <- qt(1 - (1 - level) / 2, degrees_freedom)

  data.frame(
    term = names(estimate),
    estimate = unname(estimate),
    robust_se_hc3 = standard_error,
    lower = estimate - critical * standard_error,
    upper = estimate + critical * standard_error,
    p_value = 2 * pt(abs(statistic), df = degrees_freedom, lower.tail = FALSE),
    row.names = NULL
  )
}
working_hc3 <- hc3_table(working_model)
classic_se <- summary(working_model)$coefficients[, "Std. Error"]
hc3_comparison <- data.frame(
  变量或对比 = label_terms(working_hc3$term),
  估计值 = working_hc3$estimate,
  经典标准误 = classic_se[working_hc3$term],
  HC3标准误 = working_hc3$robust_se_hc3,
  HC3下限 = working_hc3$lower,
  HC3上限 = working_hc3$upper,
  HC3_p值 = working_hc3$p_value,
  row.names = NULL,
  check.names = FALSE
)
knitr::kable(
  hc3_comparison,
  digits = 3,
  caption = "经典标准误与 HC3 稳健标准误"
)
经典标准误与 HC3 稳健标准误
变量或对比 估计值 经典标准误 HC3标准误 HC3下限 HC3上限 HC3_p值
截距 125.623 1.152 1.150 123.363 127.883 0.000
年龄（每 10 岁；中心 50 岁） 5.916 0.372 0.384 5.160 6.672 0.000
当前吸烟（参照：当前不吸烟） 4.636 1.393 1.549 1.592 7.680 0.003
年龄二次项 0.871 0.180 0.188 0.500 1.241 0.000
BMI（每 5 单位；中心 25） 3.238 0.610 0.602 2.054 4.422 0.000
男性（参照：女性） 2.610 1.033 1.051 0.544 4.676 0.013
北区（参照：中心区） 2.994 1.445 1.528 -0.010 5.998 0.051
南区（参照：中心区） -3.392 1.334 1.373 -6.091 -0.693 0.014
西区（参照：中心区） 1.762 1.422 1.377 -0.946 4.469 0.202
年龄 × 当前吸烟 1.657 0.952 0.990 -0.290 3.603 0.095

稳健标准误不能修复所有问题 HC3 只改变不确定性估计，不改变 OLS 系数。它不能修复错误的条件均值、严重测量误差、未控制混杂、选择偏倚、聚类依赖或无依据外推。聚类数据需要与聚类结构匹配的标准误或模型。

10 置信区间与预测区间

10.1 两类区间回答不同问题

  • 置信区间：给定协变量组合的人群平均收缩压有多不确定？
  • 预测区间：具有该协变量组合的一个新个体，其收缩压可能落在哪里？

预测区间还包含个体围绕条件均值的残差变异，所以通常更宽。

new_people <- data.frame(
  age = c(35, 50, 70),
  bmi = c(23, 25, 30),
  smoking_status = factor(
    c("Not current", "Current", "Not current"),
    levels = levels(ph_data$smoking_status)
  ),
  sex = factor(c("Female", "Male", "Female"), levels = levels(ph_data$sex)),
  neighborhood = factor(
    c("Central", "North", "West"),
    levels = levels(ph_data$neighborhood)
  )
)
new_people$age_c10 <- (new_people$age - 50) / 10
new_people$bmi_c5 <- (new_people$bmi - 25) / 5

mean_intervals <- predict(
  working_model, newdata = new_people,
  interval = "confidence", level = 0.95
)
individual_intervals <- predict(
  working_model, newdata = new_people,
  interval = "prediction", level = 0.95
)
interval_table <- cbind(
  new_people[c("age", "bmi", "smoking_status", "sex", "neighborhood")],
  平均拟合值 = mean_intervals[, "fit"],
  均值CI下限 = mean_intervals[, "lwr"],
  均值CI上限 = mean_intervals[, "upr"],
  个体PI下限 = individual_intervals[, "lwr"],
  个体PI上限 = individual_intervals[, "upr"]
)
knitr::kable(interval_table, digits = 1, caption = "均值置信区间与个体预测区间")
均值置信区间与个体预测区间
age bmi smoking_status sex neighborhood 平均拟合值 均值CI下限 均值CI上限 个体PI下限 个体PI上限
35 23 Not current Female Central 117 115 120 97 138
50 25 Current Male North 136 132 139 115 156
70 30 Not current Female West 146 142 149 125 166

predict.lm(..., interval = "prediction") 使用经典同方差模型。若异方差明显，不同预测位置的个体变异可能不同，经典预测区间未必校准；稳健系数标准误本身也不能自动产生可靠的个体预测区间。

simple_grid <- data.frame(age = seq(min(ph_data$age), max(ph_data$age), length.out = 160))
simple_confidence <- predict(simple_model, simple_grid, interval = "confidence")
simple_prediction <- predict(simple_model, simple_grid, interval = "prediction")
plot(
  ph_data$age, ph_data$systolic_bp,
  pch = 16, cex = 0.42, col = "grey75",
  xlab = "年龄（岁）", ylab = "收缩压（mmHg）",
  main = "置信带不等于预测带"
)
polygon(
  c(simple_grid$age, rev(simple_grid$age)),
  c(simple_prediction[, "lwr"], rev(simple_prediction[, "upr"])),
  col = rgb(230 / 255, 159 / 255, 0, 0.16), border = NA
)
polygon(
  c(simple_grid$age, rev(simple_grid$age)),
  c(simple_confidence[, "lwr"], rev(simple_confidence[, "upr"])),
  col = rgb(0, 114 / 255, 178 / 255, 0.28), border = NA
)
lines(simple_grid$age, simple_confidence[, "fit"],
      col = palette_ph["navy"], lwd = 2.5)
legend(
  "topleft",
  legend = c("拟合均值", "均值 95% CI", "个体 95% PI"),
  col = c(palette_ph["navy"], palette_ph["blue"], palette_ph["orange"]),
  lwd = c(2.5, 8, 8), bty = "n"
)
年龄与收缩压散点图，叠加较窄的均值置信带和较宽的个体预测带。

简单回归中，平均响应的 95% 置信带与新个体的 95% 预测带。

11 样本外预测表现

11.1 为什么要划分训练集与测试集？

在训练数据上加入变量总能让残差平方和不增，但新数据表现可能变差。若目标包含预测，应在建模前保留测试集，所有变量选择、变换和调参都只使用训练集。

train_data <- ph_data

prediction_formula <- systolic_bp ~
  age_c10 * smoking_status + I(age_c10^2) +
  bmi_c5 + sex + neighborhood

train_model <- lm(prediction_formula, data = train_data)
test_prediction <- predict(train_model, newdata = test_data)
split_composition <- data.frame(
  数据集 = c("开发/训练集", "保留测试集"),
  样本量 = c(nrow(train_data), nrow(test_data)),
  平均收缩压 = c(
    mean(train_data$systolic_bp),
    mean(test_data$systolic_bp)
  )
)
knitr::kable(
  split_composition,
  digits = 1,
  caption = "预先划分的训练集与首次打开的保留测试集"
)
预先划分的训练集与首次打开的保留测试集
数据集 样本量 平均收缩压
开发/训练集 420 128
保留测试集 140 131 
prediction_metrics <- function(observed, predicted) {
  c(
    RMSE = sqrt(mean((observed - predicted)^2)),
    MAE = mean(abs(observed - predicted)),
    R2_test = 1 - sum((observed - predicted)^2) /
      sum((observed - mean(observed))^2)
  )
}

model_metrics <- prediction_metrics(test_data$systolic_bp, test_prediction)
baseline_prediction <- rep(mean(train_data$systolic_bp), nrow(test_data))
baseline_metrics <- prediction_metrics(test_data$systolic_bp, baseline_prediction)
performance_table <- rbind(
  "回归模型" = model_metrics,
  "训练集均值基线" = baseline_metrics
)
knitr::kable(
  performance_table,
  digits = 2,
  caption = "保留测试集上的预测表现"
)
保留测试集上的预测表现
RMSE MAE R2_test
回归模型 10.7 8.73 0.51
训练集均值基线 15.7 11.87 -0.04

RMSE 对大误差更敏感；MAE 是绝对误差的平均；测试集 R2R^2 可以为负，表示模型比以测试集均值为基准的简单预测更差。一次随机划分本身有不确定性；正式预测研究通常需要交叉验证和外部验证。

避免数据泄漏 若先用全部数据选择变量、处理异常值或确定变换，再划分测试集，测试集已经间接参与训练。对缺失值插补、标准化和特征选择也应在每个训练折中估计规则，再应用于验证数据。

12 缺失数据与因果解释

12.1 完整案例分析改变了谁被分析

lm() 默认删除公式中任一变量缺失的行。加入体力活动后，样本量会下降：

activity_model <- lm(
  systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) +
    bmi_c5 + sex + neighborhood + physical_activity_30,
  data = ph_data,
  na.action = na.exclude
)

missing_model_comparison <- rbind(
  fit_statistics(working_model, "不含体力活动"),
  fit_statistics(activity_model, "加入体力活动；完整案例")
)
knitr::kable(
  missing_model_comparison,
  digits = 3,
  caption = "加入含缺失变量后的分析样本变化"
)
加入含缺失变量后的分析样本变化
模型 样本量 参数数 R平方 调整R平方 残差标准误
不含体力活动 420 10 0.495 0.484 10.3
加入体力活动；完整案例 391 11 0.497 0.483 10.4

直接比较两个模型系数时要小心：差异可能来自变量调整，也可能来自分析样本改变。应报告每个变量缺失量、进入模型的样本数、缺失机制假设与处理方法。

完整案例分析在缺失完全随机等较强条件下最容易解释。若缺失与已观测信息相关，可考虑多重插补；若与未观测值相关，还需要敏感性分析。本教程不以均值填补，因为单次均值填补会低估变异并扭曲变量关系。

12.2 调整变量不是越多越好

变量选择应由目标估计量和因果结构驱动：

  • 混杂变量可能需要控制；
  • 暴露之后发生的中介变量若被控制，会改变所估计的效应；
  • 同时受暴露和结局原因影响的碰撞变量，调整后可能引入偏倚；
  • 自动逐步选择会产生不稳定模型、夸大显著性并忽略领域知识；
  • 若数据来自复杂抽样、聚类或重复测量，普通 lm() 的独立误差推断不适用。

因果语言需要额外依据 “控制了年龄和 BMI 后仍显著”不足以证明因果。需要明确时间顺序、干预定义、混杂假设、可比性、缺失和选择机制，以及与研究设计相适应的估计方法。

13 如何报告线性回归

13.1 最低报告清单

一份可复核的报告至少应说明：

  1. 目标人群、分析样本、结局、主要暴露和估计目标；
  2. 连续变量单位、中心值、变换方式和分类变量参照组；
  3. 协变量选择依据，以及交互或非线性项是否预先指定；
  4. 分析样本量与每个关键变量的缺失量及处理方式；
  5. 系数、95% CI，必要时报告 p 值，但不以 p 值替代效应解释；
  6. R2R^2、调整 R2R^2 或与目标匹配的样本外性能；
  7. 残差、异方差、影响点、共线性和独立性检查；
  8. 是否使用经典或稳健标准误，以及选择理由；
  9. 研究设计、测量、选择、残余混杂和推广性的限制；
  10. 足以复现结果的软件、代码、随机种子和数据来源说明。

13.2 一个四句式结果示例

在 420 名模拟参与者中，我们用 OLS 描述收缩压与年龄、BMI、吸烟状态、生理性别及区域的条件关联，并预先允许年龄二次项和年龄—吸烟交互。对当前不吸烟、其他模型变量相同的参与者，50 岁附近年龄每增加 10 岁，平均收缩压相差 5.9 mmHg（HC3 95% CI：5.2 至 6.7）。模型的样本内调整 R2R^2 为 0.48；残差检查提示方差并非完全恒定，因此报告 HC3 稳健标准误。结果来自模拟横断面式数据，不能解释为年龄变化的因果效应，也不应外推到样本年龄和协变量范围之外。

实际报告还应解释交互，使年龄系数明确对应当前不吸烟参照组，并用预测图展示不同吸烟状态下的关系。

14 练习与答案

14.1 练习 1：识别合适的结局

以下哪个问题最适合普通线性回归？

  1. 未来 30 天是否住院；
  2. 未来 30 天住院次数；
  3. 出院时收缩压（mmHg）；
  4. 从入组到首次住院的时间。
查看答案 答案：3。 收缩压是连续结局，线性回归可描述其条件均值。二元住院结局、计数和事件时间通常分别需要与其分布及删失结构匹配的方法。

14.2 练习 2：解释连续变量系数

在 scaled_model 中，如何解释 bmi_c5 系数？它能否称为 BMI 的因果效应？

查看答案 该系数表示在年龄、吸烟状态、生理性别和区域相同的模型比较中，BMI 相差 5 个单位对应的平均收缩压差。仅凭这个观察性回归不能称为因果效应，因为混杂、选择、测量和时间顺序等条件尚未得到保证。

14.3 练习 3：解释分类变量参照组

multiple_model 中 neighborhoodSouth 的比较对象是谁？如果把 South 设为参照，模型预测会改变吗？

查看答案 它比较南区与因子的第一水平中心区，同时保持其他模型变量相同。把南区设为参照会重新表达截距与区域系数，但在同一模型空间下，拟合值、残差、R2R^2 和整体区域信息不会改变。

14.4 练习 4：自己拟合一个模型

拟合收缩压关于年龄、BMI 和吸烟状态的模型，并输出 95% CI。下面的答案代码仅使用 base R 和本教程函数。

查看答案与代码 
exercise_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status,
  data = ph_data
)
knitr::kable(
  coefficient_table(exercise_model),
  digits = 3,
  caption = "练习模型的系数与 95% 置信区间"
)
练习模型的系数与 95% 置信区间
变量或对比 估计值 标准误 CI下限 CI上限 p值
截距 129.24 0.595 128.07 130.41 0.000
年龄（每 10 岁；中心 50 岁） 5.39 0.340 4.72 6.06 0.000
BMI（每 5 单位；中心 25） 3.13 0.637 1.87 4.38 0.000
当前吸烟（参照：当前不吸烟） 4.17 1.458 1.30 7.03 0.004
年龄与 BMI 已分别按 10 岁和 5 个单位表达；吸烟系数以当前不吸烟者为参照。

14.5 练习 5：读懂交互

假设年龄主效应为 5.0 mmHg/10 岁，年龄 × 当前吸烟交互系数为 1.8 mmHg/10 岁。当前吸烟者的年龄斜率是多少？smoking_statusCurrent 主效应对应哪个年龄？

查看答案 当前吸烟者的年龄斜率为 (5.0+1.8=6.8) mmHg/10 岁。因为年龄以 50 岁为中心，吸烟主效应比较的是 50 岁、其他协变量处于同一指定值时的两组平均收缩压。

14.6 练习 6：诊断漏斗形残差

若 Residuals vs Fitted 图显示拟合值越大，残差散布越宽，你会采取哪些步骤？

查看答案 先检查结局单位、异常记录和函数形式，再查看分组残差方差以及误差变化是否有科学原因。若条件均值模型合理，可报告 HC3 等稳健标准误；若目标是个体预测，还需建立或验证随预测位置变化的方差。仅使用稳健标准误不能修复遗漏的非线性、聚类或偏倚。

14.7 练习 7：置信区间还是预测区间？

卫生部门想知道 60 岁、BMI 25、当前不吸烟女性的平均收缩压不确定性；临床人员想估计下一位具有这些特征者的收缩压范围。两者分别使用什么区间？

查看答案 前者使用条件均值的置信区间，后者使用新个体的预测区间。预测区间还包含个体残差变异，因此通常明显更宽。

14.8 练习 8：训练集 R² 很高意味着什么？

一个包含大量变量和交互的模型训练集 R2=0.90R^2=0.90，是否足以说明预测能力优秀？

查看答案 不足。训练集 R2R^2 会随模型复杂度增加而不下降，可能反映过拟合。应在未参与变量选择和调参的测试数据上评价 RMSE、MAE、校准和适用人群，并尽可能进行外部验证。

15 快速参考

15.1 常用 base R 代码

目的 代码模式
拟合线性模型 lm(y ~ x1 + x2, data = d)
查看完整摘要 summary(model)
系数 coef(model)
经典协方差矩阵 vcov(model)
95% 系数区间 confint(model)
拟合值与残差 fitted(model)；residuals(model)
标准诊断图 plot(model)
杠杆与 Cook 距离 hatvalues(model)；cooks.distance(model)
学生化残差 rstudent(model)
均值置信区间 predict(model, newdata, interval = "confidence")
个体预测区间 predict(model, newdata, interval = "prediction")
比较嵌套模型 anova(smaller, larger)
查看设计矩阵 model.matrix(model)

15.2 解释顺序

每次阅读模型结果，可按以下顺序：

问题与估计目标 → 分析样本 → 单位与参照组 → 函数形式 → 系数及区间 → 交互或非线性 → 诊断 → 样本外表现 → 局限与适用范围

15.3 核心公式

概念 公式或含义
条件均值 E(Y∣X)=XβE(Y\mid X)=X\beta
残差 ei=Yi−Ŷie_i=Y_i-\hat Y_i
OLS 目标 最小化 ∑iei2\sum_i e_i^2
OLS 解 β̂=(XTX)−1XTY\hat\beta=(X^TX)^{-1}X^TY
R2R^2 1−SSE/SST1-\mathrm{SSE}/\mathrm{SST}
RMSE n−1∑i(Yi−Ŷi)2\sqrt{n^{-1}\sum_i(Y_i-\hat Y_i)^2}
HC3 权重 ei2/(1−hii)2e_i^2/(1-h_{ii})^2
当前吸烟者年龄斜率 年龄主效应 + 年龄 × 吸烟交互效应

最后的模型检查

在发布结果前，请确认：

  • 结局是适合条件均值线性建模的连续变量；
  • 单位、编码、纳入标准和参照组已核对；
  • 函数形式与交互由问题、图形和领域知识支持；
  • 每个报告系数都能用完整句子解释；
  • 区间类型与问题一致；
  • 残差、异方差、影响点、共线性与依赖结构已检查；
  • 缺失数据和分析样本变化已报告；
  • 样本内拟合没有被误称为样本外预测能力；
  • “调整后关联”没有在缺乏依据时写成“因果效应”；
  • 代码、随机种子和软件环境足以支持复现。
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60
##  [5] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31
##  [9] lifecycle_1.0.5 cli_3.6.6       sass_0.4.10     jquerylib_0.1.4
## [13] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [17] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

本教程中的数据与结果全部为模拟内容，仅用于统计学教学。任何真实公共卫生决策都应结合研究设计、数据质量、领域证据与适当的专业审查。