11 Chapter 11: 迴歸與相關方法 (regression and correlation methods)

從一起變動到調整共變項：別讓散佈圖孤單

11.1 本章學習目標

前面幾章我們學會比較一組、兩組或多組資料，也學會處理比例與類別資料。本章開始處理另一種常見問題：兩個或多個變項是否一起變動？如果 BMI 較高，HbA1c 是否也較高？如果調整年齡與身體活動後，BMI 與 HbA1c 的關係還在嗎？

相關 (correlation) 用來描述兩個變項一起變動的方向與強度。迴歸 (regression) 則進一步建立一個模型，用一個或多個解釋變項來描述或預測結果變項。相關像是說「兩個人常一起出現」，迴歸則會多問一句：「其中一個增加一單位，另一個平均會改變多少？」

讀完本章後，你應該能夠：

解釋 Pearson correlation 與 Spearman rank correlation 的差異。
使用散佈圖 (scatter plot) 檢查兩個連續變項的關係。
建立簡單線性迴歸 (simple linear regression) 並解讀斜率。
理解截距、斜率、殘差與 \(R^2\) 的意義。
建立多變項線性迴歸 (multiple linear regression) 並解讀調整後係數。
使用殘差圖檢查線性迴歸的基本假設。
區分解釋、預測與因果推論。

11.2 相關不等於迴歸

相關係數 (correlation coefficient) 通常介於 -1 到 1。正值代表一個變項增加時，另一個也傾向增加；負值代表一個增加時，另一個傾向下降；接近 0 則表示線性或單調關係不明顯。

Pearson correlation 衡量線性關係，適合兩個連續變項且關係大致呈直線。Spearman rank correlation 則先把資料轉成排序，再衡量單調關係，對離群值與非線性但單調的關係較穩健。

相關係數有兩個重要提醒。第一，相關不代表因果。第二，相關係數只是一個摘要，不能取代圖。醫學資料很會在散佈圖裡藏劇情，單一數字常常只講了預告片。

11.3 範例資料：BMI、年齡、活動量與 HbA1c

我們建立一組糖尿病門診資料，包含年齡、BMI、身體活動分數與 HbA1c。這是模擬資料，但設計成接近臨床教學會遇到的情境。

n <- 48
# 寫死隨機生成的數值以確保範例數值不變
age <- c(60.0, 61.0, 58.0, 60.0, 66.0, 44.0, 68.0, 57.0, 65.0, 53.0, 50.0, 45.0, 60.0, 50.0, 66.0, 62.0, 60.0, 62.0, 62.0, 60.0, 46.0, 60.0, 38.0, 62.0, 49.0, 63.0, 46.0, 46.0, 76.0, 46.0, 72.0, 46.0, 73.0, 34.0, 56.0, 51.0, 64.0, 51.0, 49.0, 61.0, 55.0, 62.0, 65.0, 75.0, 46.0, 66.0, 65.0, 65.0)
bmi <- c(18.1, 28.8, 25.0, 29.8, 29.9, 29.6, 29.3, 20.1, 25.3, 30.4, 27.4, 15.3, 32.7, 20.4, 29.6, 26.4, 22.8, 22.7, 25.3, 19.5, 22.2, 24.0, 22.1, 26.9, 23.5, 21.7, 25.7, 26.9, 32.5, 33.8, 27.1, 27.8, 26.4, 24.0, 17.5, 23.1, 32.0, 33.0, 27.2, 32.3, 24.1, 22.5, 31.1, 27.7, 25.5, 30.9, 21.4, 30.1)
activity <- c(5.2, 4.7, 2.3, 7.2, 1.0, 7.4, 9.1, 6.3, 3.4, 6.6, 7.6, 10.0, 9.0, 4.3, 3.7, 3.6, 6.5, 3.8, 3.6, 7.7, 4.5, 6.8, 7.0, 7.6, 5.6, 5.5, 4.2, 5.1, 5.6, 7.6, 7.2, 6.4, 4.1, 4.0, 9.7, 1.0, 6.8, 6.9, 2.7, 5.2, 5.7, 5.3, 3.0, 8.2, 4.4, 9.3, 2.2, 3.6)
hba1c <- c(7.85, 9.98, 9.32, 8.96, 10.92, 9.3, 8.87, 7.66, 8.74, 8.38, 8.23, 6.89, 9.28, 8.31, 10.28, 9.79, 8.06, 9.0, 9.17, 9.01, 8.71, 9.08, 7.6, 8.49, 8.91, 8.48, 9.13, 7.79, 11.04, 8.31, 10.17, 8.33, 9.83, 7.88, 7.83, 8.59, 10.4, 9.32, 9.28, 9.18, 8.37, 8.8, 11.01, 9.29, 8.69, 9.72, 9.29, 9.79)

diabetes_df <- data.frame(
  age = age,
  bmi = bmi,
  activity = activity,
  hba1c = hba1c,
  stringsAsFactors = FALSE
)

head(diabetes_df)

  age  bmi activity hba1c
1  60 18.1      5.2  7.85
2  61 28.8      4.7  9.98
3  58 25.0      2.3  9.32
4  60 29.8      7.2  8.96
5  66 29.9      1.0 10.92
6  44 29.6      7.4  9.30

do.call(data.frame, 
        list(
          mean = sapply(diabetes_df, mean),
          sd = sapply(diabetes_df, sd),
          min = sapply(diabetes_df, min),
          q25 = sapply(diabetes_df, quantile, 0.25),
          median = sapply(diabetes_df, median),
          q75 = sapply(diabetes_df, quantile, 0.75),
          max = sapply(diabetes_df, max)
        )
) %>% round(2)

          mean   sd   min   q25 median   q75   max
age      57.44 9.59 34.00 49.75  60.00 64.25 76.00
bmi      26.07 4.46 15.30 22.78  26.40 29.65 33.80
activity  5.59 2.22  1.00  3.95   5.55  7.20 10.00
hba1c     8.99 0.91  6.89  8.36   8.98  9.32 11.04

11.4 Pearson 與 Spearman 相關

先看 BMI 與 HbA1c 的關係。散佈圖是第一步，相關係數是第二步。順序不要反過來；資料視覺化是臨床統計的聽診器。

p1 <- ggplot(diabetes_df, aes(x = bmi, y = hba1c)) +
  geom_point(size = 3, color = "#457b9d", alpha = 0.8) +
  geom_smooth(method = "lm", color = "#d62828", fill = "#d62828", alpha = 0.15, linewidth = 1.2) +
  labs(
    title = "BMI 與 HbA1c 的線性關係",
    x = "BMI",
    y = "HbA1c (%)"
  )

p1

Figure 11.1: BMI 與 HbA1c 的散佈圖及線性迴歸線。每個點代表一位門診病人。

pearson_result <- cor.test(diabetes_df$bmi, diabetes_df$hba1c, method = "pearson")
spearman_result <- cor.test(diabetes_df$bmi, diabetes_df$hba1c, method = "spearman", exact = FALSE)

cor_df <- data.frame(
  method = c("Pearson correlation", "Spearman rank correlation"),
  estimate = c(as.numeric(pearson_result$estimate), as.numeric(spearman_result$estimate)),
  p_value = c(pearson_result$p.value, spearman_result$p.value),
  stringsAsFactors = FALSE
)
cor_df$estimate <- round(cor_df$estimate, 4)
cor_df$p_value <- round(cor_df$p_value, 4)

cor_df

                     method estimate p_value
1       Pearson correlation   0.6222       0
2 Spearman rank correlation   0.5771       0

若 Pearson 與 Spearman 結果相近，通常表示關係大致線性且沒有被極端值嚴重支配。若兩者差很多，應回到散佈圖檢查是否有離群值、曲線關係或資料分群。

11.5 簡單線性迴歸

簡單線性迴歸用一個解釋變項 \(X\) 描述一個連續結果 \(Y\)：

\[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i \]

其中 \(\beta_0\) 是截距 (intercept)，\(\beta_1\) 是斜率 (slope)，\(\epsilon_i\) 是殘差 (residual)。在本例中，\(Y\) 是 HbA1c，\(X\) 是 BMI。斜率的解讀是：BMI 每增加 1 單位，HbA1c 平均改變多少百分點。

model_simple <- lm(hba1c ~ bmi, data = diabetes_df)
summary_simple <- summary(model_simple)

simple_coef <- data.frame(
  term = c("截距", "BMI"),
  estimate = summary_simple$coefficients[, "Estimate"],
  se = summary_simple$coefficients[, "Std. Error"],
  t = summary_simple$coefficients[, "t value"],
  p_value = summary_simple$coefficients[, "Pr(>|t|)"],
  stringsAsFactors = FALSE
)
simple_coef$estimate <- round(simple_coef$estimate, 4)
simple_coef$se <- round(simple_coef$se, 4)
simple_coef$t <- round(simple_coef$t, 4)
simple_coef$p_value <- round(simple_coef$p_value, 4)

simple_coef

            term estimate     se      t p_value
(Intercept) 截距   5.6607 0.6255 9.0494       0
bmi          BMI   0.1275 0.0237 5.3908       0

simple_r2_df <- data.frame(
  quantity = c("R-squared"),
  value = c(summary_simple$r.squared),
  stringsAsFactors = FALSE
)
simple_r2_df$value <- round(simple_r2_df$value, 4)

simple_r2_df

   quantity  value
1 R-squared 0.3872

\(R^2\) 表示模型解釋了結果變異的比例。例如 \(R^2 = 0.30\) 可解讀為此模型約解釋 HbA1c 變異的 30%。但 \(R^2\) 高不代表因果正確，\(R^2\) 低也不代表模型沒用。臨床研究裡，很多重要因素本來就只能解釋一小部分個體差異。

11.6 殘差與模型檢查

殘差是觀察值與模型預測值的差。線性迴歸常見假設包括線性、殘差變異大致相同、殘差近似常態，以及觀察值彼此獨立。

殘差圖是檢查模型的第一站。如果殘差圖出現彎曲形狀，表示線性模型可能漏掉非線性關係；如果殘差隨預測值變大而擴散，表示變異可能不等。

diabetes_df$fitted_simple <- model_simple$fitted.values
diabetes_df$resid_simple <- model_simple$residuals

p2 <- ggplot(diabetes_df, aes(x = fitted_simple, y = resid_simple)) +
  geom_point(size = 3, color = "#457b9d", alpha = 0.8) +
  geom_hline(yintercept = 0, color = "#d62828", linetype = "dashed", linewidth = 1) +
  labs(
    title = "殘差圖：檢查線性迴歸假設",
    x = "預測 HbA1c",
    y = "殘差"
  )

p2

Figure 11.2: 簡單線性迴歸的殘差圖。理想狀況下，點應大致隨機分布在 0 附近。

11.7 多變項線性迴歸

簡單迴歸只看 BMI 與 HbA1c，但臨床上 BMI 不是孤島。年齡、身體活動、用藥、病程長短都可能影響 HbA1c。多變項線性迴歸可以同時放入多個解釋變項：

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{3i} + \epsilon_i \]

在多變項模型中，BMI 係數的意思變成：在年齡與身體活動相同的情況下，BMI 每增加 1 單位，HbA1c 平均改變多少。這句「在其他變項相同的情況下」很重要，講太起就會滑倒。

model_multi <- lm(hba1c ~ bmi + age + activity, data = diabetes_df)
summary_multi <- summary(model_multi)

coef_df <- data.frame(
  term = c("截距", "BMI", "年齡", "身體活動"),
  estimate = summary_multi$coefficients[, "Estimate"],
  se = summary_multi$coefficients[, "Std. Error"],
  t = summary_multi$coefficients[, "t value"],
  p_value = summary_multi$coefficients[, "Pr(>|t|)"],
  stringsAsFactors = FALSE
)
coef_df$lower <- coef_df$estimate - 1.96 * coef_df$se
coef_df$upper <- coef_df$estimate + 1.96 * coef_df$se

coef_df %>%
  mutate(across(c(estimate, se, t, p_value, lower, upper), ~ round(., 4)))

                term estimate     se       t p_value   lower   upper
(Intercept)     截距   4.3380 0.5589  7.7611   0e+00  3.2425  5.4335
bmi              BMI   0.1037 0.0162  6.4151   0e+00  0.0720  0.1353
age             年齡   0.0473 0.0075  6.2854   0e+00  0.0325  0.0620
activity    身體活動  -0.1377 0.0315 -4.3673   1e-04 -0.1995 -0.0759

plot_coef <- coef_df %>% filter(term != "截距")
plot_coef$term <- factor(plot_coef$term, levels = c("身體活動", "年齡", "BMI"))

p3 <- ggplot(plot_coef, aes(x = estimate, y = term)) +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0.1, color = "#2a9d8f", linewidth = 1) +
  geom_point(color = "#2a9d8f", size = 3) +
  geom_vline(xintercept = 0, color = "#6c757d", linetype = "dashed", linewidth = 1) +
  labs(
    title = "多變項線性迴歸係數",
    x = "迴歸係數估計值與 95% 信賴區間",
    y = "變項"
  )

p3

Figure 11.3: 多變項線性迴歸係數與 95% 信賴區間。係數代表在其他變項固定下的平均差異。

multi_r2_df <- data.frame(
  quantity = c("Multiple regression R-squared"),
  value = c(summary_multi$r.squared),
  stringsAsFactors = FALSE
)
multi_r2_df$value <- round(multi_r2_df$value, 4)

multi_r2_df

                       quantity  value
1 Multiple regression R-squared 0.7432

多變項迴歸常用於調整混雜因子 (confounder)。但要小心，放進模型的變項不是越多越好。樣本數太小、變項太多，模型會變得不穩；如果放入中介變項或碰撞變項，也可能讓解釋偏掉。模型不是抽屜，不是什麼都塞進去就會變整齊。

11.8 非線性關係

線性迴歸的「線性」指的是模型參數的線性，不表示所有醫學關係都真的是直線。劑量反應可能有平台期、U 型關係或閾值效應。若散佈圖顯示曲線，硬套直線會讓模型看起來很有紀律，但其實在裝鎮定。

med_n <- 44
# 寫死隨機生成的數值以確保範例數值不變
dose <- c(5.0, 6.7441860465116275, 8.488372093023255, 10.232558139534884, 11.976744186046512, 13.720930232558139, 15.465116279069768, 17.209302325581397, 18.953488372093023, 20.69767441860465, 22.441860465116278, 24.186046511627907, 25.930232558139537, 27.674418604651162, 29.41860465116279, 31.162790697674417, 32.90697674418604, 34.651162790697676, 36.3953488372093, 38.13953488372093, 39.883720930232556, 41.627906976744185, 43.372093023255815, 45.116279069767444, 46.86046511627907, 48.604651162790695, 50.348837209302324, 52.093023255813954, 53.83720930232558, 55.58139534883721, 57.325581395348834, 59.06976744186046, 60.81395348837209, 62.55813953488372, 64.30232558139535, 66.04651162790697, 67.7906976744186, 69.53488372093024, 71.27906976744185, 73.02325581395348, 74.76744186046511, 76.51162790697674, 78.25581395348837, 80.0)
sbp_drop <- c(5.837844557811708, 5.907471259225027, 1.7839795746373102, 5.535283648843863, 7.565777273386345, 5.21077099591696, 8.70496105652988, 4.569595683043113, 6.916866759657563, 7.355477304672896, 7.777696678530289, 9.913619715403884, 8.986719805168795, 8.363678041818417, 9.252394005193914, 9.535810997682997, 7.368083895208159, 7.236488747764479, 7.559719019445076, 7.933943563085618, 7.02708483127353, 12.271749253180861, 6.023330468242213, 7.803720458891527, 12.166963285261163, 12.254311618559953, 9.834318622927281, 6.075174515924076, 11.788422780238829, 12.807920481663698, 9.584397770726484, 13.673986063401891, 8.875409198192381, 10.297755277770236, 9.255464794666091, 7.832140037645975, 8.537762337975726, 11.060525137961276, 12.042861969902395, 10.474019248070409, 10.427191345717237, 11.225738485144223, 6.84545858658862, 11.405038043490304)

dose_df <- data.frame(dose = dose, sbp_drop = sbp_drop, stringsAsFactors = FALSE)

head(dose_df)

       dose sbp_drop
1  5.000000 5.837845
2  6.744186 5.907471
3  8.488372 1.783980
4 10.232558 5.535284
5 11.976744 7.565777
6 13.720930 5.210771

p4 <- ggplot(dose_df, aes(x = dose, y = sbp_drop)) +
  geom_point(size = 2.5, color = "#457b9d", alpha = 0.8) +
  geom_smooth(aes(color = "線性"), method = "lm", formula = y ~ x, linewidth = 1.2) +
  geom_smooth(aes(color = "二次曲線"), method = "lm", formula = y ~ poly(x, 2), se = FALSE, linewidth = 1.2) +
  scale_color_manual(values = c("線性" = "#d62828", "二次曲線" = "#2a9d8f")) +
  labs(
    title = "非線性關係：直線不一定夠用",
    x = "藥物劑量 (mg)",
    y = "收縮壓下降量 (mmHg)",
    color = "模型"
  ) +
  theme(legend.position = "right")

p4

Figure 11.4: 藥物劑量與收縮壓下降量的非線性關係。曲線模型比單一直線更貼近資料形狀。

遇到非線性時，可考慮轉換變項、加入多項式項、使用 splines，或改用更適合的模型。重點不是追求複雜，而是讓模型尊重資料形狀。

11.9 解釋、預測與因果

迴歸可以服務不同目的。若目標是解釋，我們關心係數大小、方向與臨床意義。若目標是預測，我們關心模型對新病人的預測準確度。若目標是因果推論，就需要研究設計、時間順序、混雜控制與因果假設。

同一個迴歸表在不同目的下，解讀方式會不同。把預測模型的係數拿來當因果效果，或把因果模型拿來吹噓預測能力，都是常見的分析誤會。

11.10 常見錯誤

第一個錯誤是看到相關就說因果。BMI 與 HbA1c 相關，不代表單靠觀察資料就能證明 BMI 改變會導致 HbA1c 改變。

第二個錯誤是不看散佈圖。相關係數可能被離群值、非線性或分群資料誤導。

第三個錯誤是只報告 p 值，不報告係數與信賴區間。迴歸最有價值的是「每增加一單位，平均改變多少」。

第四個錯誤是把調整後係數解讀成「完全控制所有差異」。多變項迴歸只能調整模型中有放入、且測量品質足夠的變項。

第五個錯誤是忽略樣本數與模型複雜度。資料少、變項多，模型就容易不穩，信賴區間也會變寬。

11.11 本章重點整理

相關描述兩個變項一起變動的方向與強度；迴歸描述一個或多個解釋變項與結果變項之間的平均關係。Pearson correlation 適合線性關係，Spearman correlation 適合單調關係。簡單線性迴歸給出截距、斜率、殘差與 \(R^2\)；多變項迴歸可以在模型內調整其他變項。

不論模型多漂亮，都要回到研究問題、資料圖形與臨床意義。迴歸不是把資料送進機器後自動吐出真理；它比較像一份很有用但需要人類監督的病歷摘要。

11.12 小練習

使用一組包含年齡與收縮壓的資料，畫散佈圖並計算 Pearson correlation。
若散佈圖呈現單調但非線性關係，你會選 Pearson 還是 Spearman？為什麼？
建立「LDL-C = 截距 + BMI」的簡單線性迴歸，解釋 BMI 係數。
在上述模型加入年齡與性別後，BMI 係數變小。請提出兩個可能原因。
殘差圖出現漏斗形狀時，可能代表什麼問題？
請用自己的話說明「預測模型表現好」與「模型係數有因果意義」為何不同。

11.13 Glossary

中文術語	English term	說明
相關	correlation	描述兩個變項一起變動方向與強度的統計概念。
Pearson 相關	Pearson correlation	衡量兩個連續變項線性關係的相關係數。
Spearman 等級相關	Spearman rank correlation	以排序衡量兩個變項單調關係的相關係數。
迴歸	regression	用一個或多個解釋變項描述或預測結果變項的方法。
簡單線性迴歸	simple linear regression	只有一個解釋變項的線性迴歸。
多變項線性迴歸	multiple linear regression	包含多個解釋變項的線性迴歸。
截距	intercept	解釋變項為 0 時，模型預測的結果值。
斜率	slope	解釋變項增加一單位時，結果變項平均改變的量。
殘差	residual	觀察值與模型預測值的差。
\(R^2\)	\(R^2\)	模型解釋結果變異比例的摘要指標。
混雜因子	confounder	同時與暴露及結果相關、可能扭曲關係估計的變項。
預測	prediction	使用模型估計新觀察值的結果。
因果推論	causal inference	評估某暴露或介入是否造成結果改變的推論框架。