限制性立方樣條與劑量反應曲線

這一頁在補一個一直被推薦、卻沒有人示範過的東西

站上已經在三個地方叫你「改用限制性立方樣條」：線性迴歸 談非線性時、迴歸模型診斷看殘差圖時、 以及 Cox 比例風險模型檢查連續變項假設時。這一頁是那三個指標的落點。

先講清楚要放掉的是哪一條假設。把一個連續變項直接丟進迴歸式：

這裡的 $\beta \cdot x$ 隱含一句很強的話：$x$ 每增加一單位，對結果的影響都一樣大—— 從 0 到 1 顆淋巴結，跟從 15 到 16 顆，對 log hazard 的貢獻被強制設成同一個數字。 臨床上幾乎沒有人相信這件事，但它是預設值，而且不會有任何警告。

限制性立方樣條（restricted cubic spline, RCS；在 R 裡最常用的實作是 splines::ns()， natural spline）的做法是：把 $x$ 的範圍用幾個節點（knot）切成幾段，每段配一條三次曲線， 在 knot 處要求接得平滑，並且限制兩端為直線（這就是「restricted」與「natural」的意思—— 不加這個限制，資料最稀疏的兩端會翹得最厲害）。

這一頁的例子

survival::rotterdam 是荷蘭 Rotterdam 腫瘤中心的乳癌世代， 2982 位病人、1713 個 recurrence-free survival 事件 （復發或死亡，取先發生者）。暴露變項是 nodes，腋下陽性淋巴結數。

選這個變項有三個理由，而且三個都是資料本身的性質，不是為了舉例湊出來的：

1. 它與復發的關係真的是非線性的——前幾顆陡升、之後轉平，這正是直線畫不出來的形狀。
2. 2982 個人裡有 1436 個的 nodes 是零（48.2%）， 這個零膨脹會把教科書的 knot 放置規則撞壞，見下一節。
3. 零本身是臨床上認得的參考點（node-negative disease），曲線有一個讀者一看就懂的錨。

library(survival)
library(splines)
data(cancer, package = "survival")

rot <- rotterdam
rot$rfs     <- pmax(rot$recur, rot$death)
rot$rfstime <- ifelse(rot$recur == 1, rot$rtime, rot$dtime)

# 線性版本：一個參數，強制等距效果
fit_lin <- coxph(Surv(rfstime, rfs) ~ nodes + age + size + meno, data = rot)

# 樣條版本：knot 自己指定（為什麼不用預設，見下一節）
basis   <- ns(rot$nodes, knots = c(1, 4, 9), Boundary.knots = c(0, 20))
fit_spl <- coxph(Surv(rfstime, rfs) ~ basis + age + size + meno, data = rot)

# 完全不放 nodes，用來檢定「整體關聯」
fit_non <- coxph(Surv(rfstime, rfs) ~ age + size + meno, data = rot)

anova(fit_non, fit_spl)   # 整體關聯：有沒有關係
anova(fit_lin, fit_spl)   # 非線性：關係是不是直線

import numpy as np
import pandas as pd
import statsmodels.api as sm
from patsy import dmatrix

# statsmodels 的 get_rdataset("rotterdam", "survival") 會失敗：Rdatasets 的索引
# 把它登記在 survival 的 cancer bundle 底下，查不到獨立的 rotterdam 條目。
# 直接讀 CSV 最省事。
URL = "https://vincentarelbundock.github.io/Rdatasets/csv/survival/rotterdam.csv"
rot = pd.read_csv(URL)
rot["rfs"] = np.maximum(rot["recur"], rot["death"])
rot["rfstime"] = np.where(rot["recur"] == 1, rot["rtime"], rot["dtime"])

# constraints="center" 不能省。cr() 未加約束的基底每一列總和恰為 1，
# 而 Cox model 沒有截距，那個方向不可辨識 —— 少了它，這一頁唯一在談的
# nodes 那幾列 SE 全是 nan（係數照印，推論欄整排空白）。
# lower_bound / upper_bound 對齊 R 的 Boundary.knots = c(0, 20)。
SPL = 'cr(nodes, knots=[1, 4, 9], lower_bound=0, upper_bound=20, constraints="center")'

X = dmatrix(SPL + " + age + C(size) + meno",
            rot, return_type="dataframe").drop(columns="Intercept")
# ties="efron" 也對齊 R：coxph() 預設 Efron，PHReg 預設 Breslow。
fit = sm.PHReg(rot["rfstime"], X, status=rot["rfs"], ties="efron").fit()
print(fit.summary())

# 非線性檢定要自己組：兩個模型的 log-likelihood 相減再乘 2
X_lin = dmatrix("nodes + age + C(size) + meno",
                rot, return_type="dataframe").drop(columns="Intercept")
fit_lin = sm.PHReg(rot["rfstime"], X_lin, status=rot["rfs"], ties="efron").fit()
lrt = 2 * (fit.llf - fit_lin.llf)   # 99.686，df = 3

knot 放在哪裡，以及零膨脹的變項會出什麼事

最常被引用的規則來自 Harrell：三個 knot 放在第 10、50、90 百分位，四個放 5/35/65/95， 五個放 5/27.5/50/72.5/95。這條規則有它的道理——knot 放在資料密的地方，每一段才有足夠的人撐著。

但它是對「分佈連續」的變項說的。 這份資料的 nodes 分佈長這樣：

前兩個百分位都是零，因為將近一半的人（1436 位，48.2%） 沒有陽性淋巴結。照百分位放 knot，等於把好幾個 knot 疊在同一個點上，而那個點又剛好是下邊界。 R 的反應是自己把 knot 挪開，並丟一句警告：

所以這一頁改成明確指定：ns(nodes, knots = c(1, 4, 9), Boundary.knots = c(0, 20))。 內部 knot 選 1、4、9 顆 （臨床上會拿來分期的位置，也都落在資料密的地方）， 邊界訂在 0 與 20 顆—— 上邊界不取最大值，因為超過那裡的人太少，那一段畫出來的是外插不是估計。

參考點：曲線在講「跟誰比」

樣條的係數本身沒有辦法唸——basis1 到 basis4 各自不對應任何臨床量。 能唸的是曲線：固定其他共變項，把 $x$ 從參考值移到某個值，風險變成幾倍。

所以每一張劑量反應曲線都得先選一個參考點。本頁選 nodes = 0 （node-negative），理由是臨床上認得。常見的其他選法是中位數或第 10 百分位； 選哪一個不影響曲線的形狀，只是把整條曲線上下平移，也不影響任何檢定的 p 值。

兩個 p 值，兩個問題，一定要分開報

這是審稿最常抓的一點，也是 spec 裡點名的要求。樣條模型可以做兩個不同的概似比檢定：

兩個問題是獨立的，四種組合都可能發生：

• 整體顯著、非線性不顯著：關係存在，而且本研究沒有偵測到偏離直線。線性模型夠用。
• 整體顯著、非線性也顯著（本頁的情形）：關係存在，而且直線描述不了它。
• 整體不顯著、非線性顯著：少見但會發生——U 形關係的兩端往上、中間往下， 用一個線性係數去總結會互相抵消。只報線性模型的 p 會漏掉它。
• 兩個都不顯著：本研究未偵測到關聯。注意這不等於「沒有關聯」， 只代表這份資料沒有提供足以偵測的證據。

線性模型到底錯在哪裡

把線性模型隱含的那條直線疊上去，就看得到代價。

三件事值得停下來看：

1. 低估最嚴重的是中段，不是最前面。 直線與曲線的比值在約 7.6 顆處掉到最低， 那裡線性模型給的 HR 只有樣條的 0.57 倍；1 顆處的比值是 0.88，落差小得多。 圖上標出的 4 顆（樣條 2.09、線性 1.34） 取的是臨床上病人最多的位置，不是落差最大的位置。 而從 1 顆到 15 顆這一整段，直線都落在樣條的信賴區間之外。
2. 尾端反過來被高估。 20 顆時線性模型外推到 4.37， 樣條給 3.55（2.80–4.50）。 這裡直線還在區間內，但那是因為那一段的人太少、區間已經寬到分不出兩者—— 區間寬不是兩個模型一致的證據。
3. 斜率在前幾顆之後就明顯放緩。 從 0 到 4 顆，HR 從 1 升到 2.09； 從 4 到 9 顆再升到 3.39；而從 9 顆一路到 15 顆，只再升到 3.84，之後就不再上升。 臨床上這是說得通的：淋巴結轉移到某個程度之後，再多幾顆帶來的額外訊息有限。 線性模型結構上沒有能力表達這種轉平，它只能繼續往上。

這一頁與其他頁的關係

• 想知道「那乾脆切成四分位不就好了」為什麼通常更糟—— 見連續變項要不要切組。 簡短版：切組是拿檢定力去換掉線性假設，而樣條可以只換掉假設、不付那個代價。
• 殘差圖看出彎曲之後要怎麼辦——回到迴歸模型診斷。
• Cox 模型的連續變項假設——見 Cox 比例風險模型。 要注意這條假設（log hazard 對共變項為線性）與比例風險假設是兩件不同的事， 用樣條處理前者，不會修好後者。
• 樣條也會出現在校準那一頁，但角色不同： 那裡是用樣條去畫「預測風險 vs 實際風險」的彈性校準曲線，暴露變項是模型自己的線性預測值。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-05-splines.R

回報這一頁的問題