## はじめに
連続変数（年齢、BMI、検査値など）とアウトカムの関連を調べるとき、線形の仮定を置くと実態を見逃すことがあります。スプライン回帰を使えば、変数のどの範囲でリスクが高いのかを柔軟にモデル化し、ハザード比の曲線として可視化できます。
この記事では `mkspline` で線形スプライン変数を作成し、Cox回帰で推定したあと、`xblc` コマンドでハザード比を図示する手順を紹介します。
## 必要なパッケージ
`xblc` は外部コマンドです。初回のみインストールが必要です。福井大学内からインストールする場合はプロキシ設定が必要です（[[Stata - 福井大学内でのパッケージインストール]] を参照）。
```stata
net install st0215_1, from(http://www.stata-journal.com/software/sj11-3)
```
## データの準備
心臓移植データ `stan3` を使用します。
```stata
clear
webuse stan3
stset stime, failure(died)
```
```
Survival-time data settings
Failure event: died!=0 & died<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
172 total observations
0 exclusions
--------------------------------------------------------------------------
172 observations remaining, representing
75 failures in single-record/single-failure data
60,591 total analysis time at risk and under observation
```
## 線形スプライン変数の作成
`mkspline` は指定したノット（結び目）で変数を区間ごとに分割し、区分線形のスプライン変数を生成します。ここでは年齢を 30、40、50歳の3点で分割し、4つのスプライン変数を作ります。
```stata
mkspline age_s1 30 age_s2 40 age_s3 50 age_s4 = age
```
これにより以下の変数が作成されます。
| 変数 | 区間 |
|------|------|
| `age_s1` | age < 30 の部分の傾き |
| `age_s2` | 30 ≤ age < 40 の部分の傾き |
| `age_s3` | 40 ≤ age < 50 の部分の傾き |
| `age_s4` | 50 ≤ age の部分の傾き |
各区間で異なる傾きを持つことで、年齢とハザードの関連が直線に制約されなくなります。
## 共変量の中心化
調整変数は平均値で中心化しておくと、`xblc` での基準値の解釈がしやすくなります。
```stata
sum year
gen year_c = year - r(mean)
```
複数の変数を中心化したい場合はループを使います。
```stata
foreach var of varlist age height weight {
quietly sum `var'
gen `var'_c = `var' - r(mean)
}
```
## Cox回帰の実行
スプライン変数と中心化した共変量を投入します。
```stata
stcox age_s* year_c
```
```
Cox regression with Breslow method for ties
No. of subjects = 172 Number of obs = 172
No. of failures = 75
Time at risk = 60,591
LR chi2(5) = 15.52
Log likelihood = -336.18518 Prob > chi2 = 0.0084
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age_s1 | .912946 .042566 -1.95 0.051 .8332165 1.000305
age_s2 | 1.114656 .0776674 1.56 0.119 .972368 1.277766
age_s3 | .9969981 .0430169 -0.07 0.944 .9161531 1.084977
age_s4 | 1.064414 .0474151 1.40 0.161 .9754236 1.161524
year_c | .852065 .0598065 -2.28 0.023 .742552 .9777293
------------------------------------------------------------------------------
```
各スプライン変数の係数は「その区間における年齢1歳あたりのハザード比」を表しています。個別の係数よりも、全体の曲線パターンを見ることが重要です。
## xblc によるハザード比の算出
`xblc` は回帰モデルの結果からスプライン変数を合成し、特定の基準値に対するハザード比と信頼区間を算出します。
```stata
levelsof age
quietly: xblc age_s*, covname(age) at(`r(levels)') reference(40) eform ///
line yscale(log) generate(pa hr lb ub)
```
主なオプションの意味は以下の通りです。
| オプション | 説明 |
|------------|------|
| `covname(age)` | 元の変数名を指定する |
| `at(...)` | ハザード比を計算する年齢の値を指定する |
| `reference(40)` | 基準値を40歳に設定する（HR = 1 となる点） |
| `eform` | 指数変換してハザード比として表示する |
| `generate(pa hr lb ub)` | 予測値の変数を生成する（pa: 年齢値、hr: ハザード比、lb/ub: 信頼区間下限/上限） |
## グラフの作成
`xblc` が生成した変数を使って、信頼区間付きのハザード比曲線を描きます。
```stata
* ノット位置と両端にドットを表示する
quietly sum pa
gen show_dot = 1 if pa==`r(min)' | pa==`r(max)' | pa==30 | pa==40 | pa==50
twoway rarea lb ub pa, col(gs13) || ///
scatter hr pa if show_dot==1, ///
msymbol(circle) mcol(black) msize(medium) connect(l) ||, ///
yscale(log) ///
ylab(0.1 0.2 0.5 1 2 5 10, nogrid angle(horizontal) format(%5.1f)) ///
legend(off) scheme(s2mono) ///
ytitle("Hazard ratio") xtitle("Age") ///
name(f1, replace)
```
![[spline_hr_plot.png|500]]
灰色の帯が95%信頼区間、黒い線と点がハザード比を表しています。基準値（40歳、HR = 1）から離れるほど信頼区間が広がっていることがわかります。
## グラフの読み方
- **HR = 1 の線**が基準です。40歳を基準にすると、30歳付近では HR ≈ 0.3 で死亡リスクが低く、60歳以上では HR > 2 でリスクが高いことが読み取れます。
- **信頼区間の幅**はその年齢帯のデータ量を反映しています。端の年齢ではデータが少ないため帯が広がります。
- **ノットの位置**（30、40、50歳）で傾きが変化しています。ノットの選び方で曲線の形が変わるため、臨床的知識やデータの分布に基づいて設定します。
## ノットの選び方
ノットの数と位置はモデルの柔軟性を左右します。一般的な方針は以下の通りです。
- データの分布（パーセンタイル）に基づいて配置する（例: 25th、50th、75th）
- 臨床的に意味のあるカットポイント（例: 30、40、50歳）を使う
- ノットが多すぎると過適合、少なすぎると非線形性を捉えられないため、3-5個程度が一般的です
## 多変量調整での応用
実際の研究では複数の交絡因子を調整します。調整変数を中心化してモデルに追加すれば、同じ手順でスプライン曲線を描けます。
```stata
* 複数変数の中心化
foreach var of varlist age height weight {
quietly sum `var'
gen `var'_c = `var' - r(mean)
}
* 調整済みモデル
stcox age_s* year_c height_c weight_c
```
`xblc` は最後に推定したモデルの結果を使うため、追加の変数が調整変数として反映されます。
## 制限付き3次スプライン（RCS）との違い
`mkspline` で作成するのは線形スプラインです。より滑らかな曲線が必要な場合は、制限付き3次スプライン（restricted cubic spline）を使う方法もあります。
```stata
mkspline age_rcs = age, cubic nknots(4)
```
`cubic` オプションを加えると3次スプライン変数が作成されます。`xblc` はどちらのスプラインにも対応しています。
## まとめ
| 手順 | コマンド |
|------|----------|
| スプライン変数の作成 | `mkspline` |
| 生存時間回帰 | `stcox` |
| ハザード比の算出と可視化 | `xblc` |
| グラフの描画 | `twoway rarea` + `scatter` |
スプライン回帰は、連続変数を安易にカテゴリ化することなく、非線形の用量反応関係を探索できる有用な手法です。