在醫學研究中，我們經常構建迴歸模型來分析自變量和因變量之間的關係。事實上，大多數的迴歸模型有一個重要的假設就是自變量和因變量呈線性關聯，這個條件實際很難滿足。常見的解決方法是將連續變量分類，但類別數目和節點位置的選擇往往帶有主觀性，並且分類往往會損失信息。因此，一個更好的解決方法是擬合自變量與因變量之間的非線性關係，

近年來在Lancet、BMJ等雜誌經常見到利用限制性立方樣條來擬合非線性關係。在之前醫咖會R語言入門課程第18課《展示非線性關係》中也簡單介紹過限制性立方樣條的應用（R課程第18期：多項式迴歸、分段迴歸、限制性立方樣條...）。

本文將以Lancet Diabetes Endocrinol雜誌2018年的文章《Association of BMI with overall and cause-specific mortality: a population-based cohort study of 3.6 million adults in the UK》和作者模擬的數據為例，帶領大家學習限制性立方樣條的使用和結果解讀。

在這篇文章中，作者使用限制性立方樣條繪製BMI與死因之間的關係，第一個圖形如下：

從圖中我們可以看到，不管是在全人群還是非吸菸人群中，BMI與all-cause, communicable, and non-communicable disease mortality之間呈J型關係，死亡率最低點的BMI大概是25Kg/m²。另外，對於injuries and external causes mortality而言，隨著BMI的增加，死亡率在降低，但當BMI超過27Kg/m²繼續增加時，死亡率開始增加，但是增加幅度很小。

通過以上的例子，我們可以看到，使用限制性立方樣條，可以很清晰的描述自變量與因變量之間的關係。事實上，限制性立方樣條的應用範圍非常廣，凡是想描述自變量和因變量的關係都可以在迴歸模型中加入限制性立方樣條，除了以上的例子中的COX迴歸外，還可以應用到線性迴歸、Logistic迴歸、以及Meta分析中劑量反應關係的Meta迴歸等。

什麼是立方樣條？

迴歸樣條(regression spline)本質上是一個分段多項式, 但它一般要求每個分段點上連續並且二階可導，這樣可以保證曲線的平滑性。而限制性立方樣條是在迴歸樣條的基礎上附加要求：樣條函數在自變量數據範圍兩端的兩個區間內為線性函數。

在利用限制性立方樣條繪製曲線關係時，通常需要設置樣條函數節點的個數(k)和位置（ti）。絕大多數情況下, 節點的位置對限制性立方樣條的擬合影響不大, 而節點的個數則決定曲線的形狀, 或者說平滑程度。當節點的個數為2時, 得到的擬合曲線就是一條直線，

在《Regression Modeling Strategies》這本書中，Harrell建議節點數為4時，模型的擬合較好，同時可以兼顧曲線的平滑程度和避免過擬合造成的精度降低。而當樣本量較大時，例如因變量為未刪失的連續變量並且大於100時，5個節點是更好的選擇。小樣本（如n<30）可以選擇3個節點。以下是Harrell推薦的節點數和相應的節點位置，大家可以參考。

案例說明（模擬數據）

目前SAS、STATA、R等軟件都可以進行限制性立方樣條分析。基於畫圖的方便，我們以R語言為例進行說明。首先參照rms包，生成一個模擬數據集，包括性別（sex），年齡(age)以及生存時間(time)和結局變量(death)。

若想分析年齡和生存率之間關係，傳統的方法可以在Cox迴歸中將年齡作為連續變量處理，也可以對年齡進行分組，這樣的做法都無法更直觀的呈現年齡與死亡風險之間的關聯。以下我們用限制性立方樣條來分析年齡與死亡風險之間的關係：

#####加載所需要的包

library(ggplot2)

library(rms) #####立方樣條所需要的包

###參照rms包，先生成一個模擬數據集，包括性別（sex），年齡(age)以及生存時間(time)和結局變量(death);

n

set.seed(731)

age

label(age)

sex

cens

h

time

label(time)

death

time

units(time)

data

#######開始正式畫圖

dd

options(datadist='dd') #為後續程序設定數據環境

####擬合cox迴歸模型，注意這裡的R命令是“cph”，而不是常見的生存分析中用到的“coxph"命令

fit

dd$limits$age[2]

fit=update(fit)

HR

P1

P1

######自己畫圖

P2

geom_ribbon(data=HR, aes(age,ymin = lower, ymax = upper),alpha = 0.1,fill="red")

######進一步設置圖形

P2

labs(title = "RCS", x="age", y="HR (95%CI)")

P2

####如果想看是否存在非線性關係，可以使用anova

anova(fit)

可以看到age整體是有意義的（包括線性或者非線性關聯），然後看P-Nonlinear =0.0168<0.05，這裡我們可以說年齡與死亡風險之間存在非線性關聯。

如果自變量與關注的結局變量存在非線性關係，如何在文章中對結果更詳細的描述呢，建議大家可以參照上文中提到的Lancet的文章。

個人還發現BMJ的一篇文章《Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study》對於非線性關係描述的非常好，摘抄一部分放在這裡供大家參考：

In figure 1, we used restricted cubic splines to flexibly model and visualize the relation of predicted fat mass and lean body mass with all cause mortality in men. The risk of all cause mortality was relatively flat until around 21 kg of predicted fat mass and then started to increase rapidly afterwards (P for non-linearity <0.001). The average BMI for men with 21 kg of predicted fat mass was 25. Above 21 kg, the hazard ratio per standard deviation higher predicted fat mass was 1.22 (1.18 to 1.26). Regarding the strong U shaped relation between predicted lean body mass and all cause mortality, the plot showed a substantial reduction of the risk within the lower range of predicted lean body mass, which reached the lowest risk around 56 kg and then increased thereafter (P for non-linearity <0.001). Below 56 kg, the hazard ratio per standard deviation higher predicted lean body mass was 0.87 (0.82 to 0.92).

那麼在方法部分如何描述使用了限制性立方樣條，還是可以參照這篇文章：We also used restricted cubic splines with five knots at the 5th, 35th, 50th, 65th, and 95th centiles to flexibly model the association of lean body mass, fat mass, and BMI with mortality.

以上就是對限制性立方樣條的簡單介紹，原理和操作都比較簡單。但加入到分析中，能更直觀的描述感興趣的自變量和因變量之間的關係，發現更有趣的點，可以為文章增色不少。對R不熟悉也不要緊，在sas和stata中也可以實現，感興趣的同學可以去嘗試。

[1] Bhaskaran K, Dos-Santos-Silva I, Leon DA, Douglas IJ, Smeeth L. Association of BMI with overall and cause-specific mortality: a population-based cohort study of 3·6 million adults in the UK. Lancet Diabetes Endocrinol. 2018;6(12):944–953. doi:10.1016/S2213-8587(18)30288-2.

[2] Lee DH, Keum N, Hu FB, et al. Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study. BMJ. 2018;362:k2575. Published 2018 Jul 3. doi:10.1136/bmj.k2575

[3] 羅劍鋒, et al. "限制性立方樣條在非線性迴歸中的應用研究%The Application of Restricted Cubic Spline in Nonlinear Regression." 中國衛生統計 027.003(2010):229-232.

1. R課程第18期：多項式迴歸、分段迴歸、限制性立方樣條...

2. 建立非線性迴歸預測模型，來看R教程！

3. 如何理解線性迴歸的方差齊性檢驗

關注醫咖會，及時獲取最新統計教程！

閱讀更多 醫咖會 的文章

關鍵字: 立方 樣條 R語言