「先生！BMIと死亡リスクの関係を見たくて多変量Cox比例ハザードモデルを組んだのですが、$P > 0.05$ で有意差が出ませんでした。でもデータをグラフで見ると、どう見てもU字型（低体重でも肥満でも死亡リスクが高い）なんです！直線モデルのまま論文に書いて良いのでしょうか？」
臨床研究において、連続変数（BMI、血圧、eGFR、CRP、年齢など）と疾患リスクの関係を解析する際、このような「手詰まり」に直面する研究者は非常に多い。
生体反応の多くは直線（線形）ではなく、「適正範囲」や「閾値（Threshold）」を持つ非線形な関係（U字型・J字型・L字型）を示す。これを無理やり直線モデルに当てはめると、リスクの上昇と低下が相殺されて「有意差なし」という誤った結論を導いてしまう。
かといって、連続変数を「4分位群（Q1〜Q4）」や「65歳以上/未満」などのカテゴリーに安易に分割して逃げると、情報損失が大きく査読者から「なぜそのカットオフなのか？ P-hacking（データ駆動的なカットオフ選択）ではないか？」と一発で叩かれる原因になる。
この問題を根本から解決し、現代医学統計のデファクトスタンダードとなっている手法が Restricted Cubic Splines（RCS / 制限付き3次スプライン） である。
本記事では、連続変数をカテゴリー化する「3つの罠」から、RCSが世界標準である理由、査読者を納得させる「非線形性の検定（$P$ for non-linearity）」の読み方、rms パッケージを用いたRでの実践手順、そして論文でそのまま使える英文テンプレートまでを徹底解説する。
1. 連続変数をカテゴリー化して解析する「3つの罠」
臨床データを解析する際、連続変数の扱い方を誤ると研究の質が著しく低下する。臨床医が陥りがちな「3つの罠」を整理する。
罠1：線形性の強制適用（Linearity Trap）
標準的なCox比例ハザードモデルやロジスティック回帰モデルに連続変数をそのまま投入すると、「1単位増加するごとのリスク変化（ハザード比：HR）が常に一定である」という直線（線形）モデルが自動的に適用される。
例えば、BMIと全死亡リスクのように「低体重（BMI < 18.5）でも高度肥満（BMI > 30）でも死亡率が上がる」U字型リスク構造を持つデータを直線モデルに入れると、低体重側のリスク上昇と肥満側のリスク上昇が互いに打ち消し合い、回帰係数の傾きが平坦（$\text{HR} \approx 1.0$）になって $P > 0.05$（有意差なし）という偽陰性の結論に至る。
罠2：安易なカテゴリー化（Arbitrary Categorization Trap）
直線モデルで有意差が出ないため、連続変数を「4分位群（Q1〜Q4）」や「臨床的カットオフ（例: 65歳以上）」でカテゴリー化して多変量解析に入れる手法がよく使われる。
しかし、カテゴリー化には以下の致命的なデメリットが存在する：
- 統計的検出力（Power）の著しい低下: 連続値を区分化することで、データが持つ豊富な情報量が大量に破棄される。
- 不自然な境界効果: カットオフの境界線（例: 64歳と65歳）でリスクが急激に跳ね上がるという、臨床的にあり得ない不自然なモデルになる。
- 査読者からの指摘: 「なぜそのカットオフを設定したのか？先行研究の根拠はあるのか？解析結果が良く見えるカットオフを探したのではないか？」という指摘を受ける最大の要因となる。
罠3：高次多項式の挿入（Polynomial Trap）
曲線を描くために、$X$ に加えて $X^2$（2次式）や $X^3$（3次式）をモデルに入れる手法もある。
しかし、高次多項式はデータの両端（最小値・最大値付近）で曲線が急激に大きく飛び跳ねる現象（ルンゲ現象）を起こしやすい。その結果、データが存在しない極端な値でのリスク予測が著しく不安定になり、臨床的な解釈に耐えられなくなる。
2. なぜ Restricted Cubic Splines (RCS) が現代医学統計のデファクトスタンダードなのか？
カテゴリー化や高次多項式の罠を全て回避し、連続変数の真のリスク構造を客観的に描き出す手法が Restricted Cubic Splines (RCS / 制限付き3次スプライン) である。
[多項式 (Polynomial)] [通常の Spline] [RCS (制限付き3次スプライン)]
端で急激に飛び跳ねる 区間ごとに滑らかに接続 端（両端のKnotsの外側）を
＼ ／ 〜〜〜〜 「直線（線形）」に制限！
＼______／ → 安定したカーブを描画
RCSの優れた仕組み
スプライン（Spline）とは、データの範囲を複数の区間に分け、各区間を3次関数（Cubic）で滑らかに繋ぎ合わせる技法である。
RCSの最大の特徴は 「Restricted（制限付き）」 という制約にある。データの両端（最も外側の結節 / Knotsの外側）において、関数を「3次式」ではなく「1次式（直線）」に強制制限する。これにより、データの端で曲線が異常に跳ね上がる多項式の弱点を完全に克服し、安定したリスクカーブを描画することができる。
Knots（結節）の決定指針
スプライン曲線を繋ぎ止めるポイントを Knots（結節） と呼ぶ。
結節の数は 3〜5個 が推奨されており、4個 がモデルの適合度と複雑さのバランスにおいてベストとされることが多く、標準的な第一選択となる。
重要であるのは、結節の位置は研究者が手動で選ぶのではなく、データのパーセンタイル（分位点：例 5%, 35%, 65%, 95%）に基づいて自動配置する点である。これにより、研究者の主観や P-hacking の介入を完全に排除できる。
基準点（Reference value）に基づく直感的解釈
RCSを用いることで、「特定の基準点（例: BMI = 22.0 $\text{kg/m}^2$）におけるハザード比を $\text{HR} = 1.0$ と設定したとき、BMIの変化に伴う連続的なハザード比の推移（95%信頼区間付き）」を1つの連続カーブとして提示できる。
これにより、「BMI 22 付近で最も死亡リスクが低く、20 未満および 28 以上で統計学的に有意にリスクが上昇する」といった臨床的にきわめて明快な解釈が可能となる。
※補足（関連記事）:
ロジスティック回帰（二値アウトカム）におけるRCS曲線の具体的な描き方については ロジスティック回帰におけるRCS曲線描画 を参照されたい。
3. 査読者を黙らせる「非線形性の検定（P for non-linearity）」の読み方
論文にRCSの綺麗なグラフを載せるだけでは、統計査読者を完全に納得させることはできない。重要なのは、「この曲線関係は、統計学的に意味のある非線形性と言えるのか？」を検定で証明することである。
全体効果（All effects） vs 非線形効果（Non-linear effect）の分離
Rの rms パッケージで anova() を実行すると、RCSモデルに含まれる各変数の効果を以下の2つに分離して $P$ 値を算出してくれる：
- Overall（全体としての関係性）: 変数がアウトカムに何らかの影響を与えているか（直線か曲線かを問わない全体効果の検定）。
- Non-linear（非線形効果）: 直線モデル（1次式）からのずれ（歪み）が統計的に有意かどうかの検定（$P$ for non-linearity）。
論文での主張シナリオ
- $P$ for non-linearity $< 0.05$ の場合:「連続変数とアウトカムの間には明確な非線形関係（U字型や閾値効果など）が存在する」と自信を持って宣言する。単一のハザード比（例: $\text{HR} = 1.02$）を記載することは不適切であり、RCSによる連続HR曲線（図）を提示することが主たる解析結果となる。
- $P$ for non-linearity $\ge 0.05$ の場合:「非線形性は統計学的に有意ではなかった（直線関係で近似可能）」と判断し、通常の直線モデル（1単位あたりの単一HR）を採用する根拠とする。
4. 【実践Rコード】rms パッケージによる Cox回帰 ＋ RCS モデル構築と連続HR曲線の可視化（Step by Step）
ここからは、医学統計における世界標準である Frank Harrell 教授の rms パッケージ を用いて、Cox回帰でのRCSモデル構築から非線形性検定、美麗な連続HR曲線の描画までを Step by Step で解説する。
Step 1: 疑似臨床データの作成
患者500人分、追跡時間（time）、死亡イベント（status）、BMI、年齢、性別のデータを生成する。BMI 22.0 近辺で最も死亡リスクが低く、低体重および肥満でリスクが上昇するU字型生存データを意図的に設計する。
library(rms)
library(tidyverse)
# 乱数シードの設定
set.seed(123)
n <- 500
# 疑似データの生成
bmi <- runif(n, min = 18, max = 35)
age <- round(rnorm(n, mean = 65, sd = 10))
gender <- factor(sample(c("Male", "Female"), n, replace = TRUE))
# U字型のリスク（BMI=22で最小ハザード）を表現する対数ハザードの生成
log_hazard <- 0.02 * (bmi - 22)^2 + 0.03 * (age - 65) + 0.3 * (gender == "Male") - 1.5
lambda <- exp(log_hazard)
# 指数分布に基づく生存時間と打切りの生成
time <- rexp(n, rate = lambda * 0.05)
status <- if_else(time > 36, 0, 1)
time <- pmin(time, 36)
df_rcs <- tibble(
id = 1:n,
time = round(time, 1),
status = status,
bmi = round(bmi, 1),
age = age,
gender = gender
)
# データの確認
head(df_rcs)
Step 2: rms::datadist によるデータ分布情報の設定と基準値指定
rms パッケージを使用する際は、モデルを構築する前に datadist() を用いてデータの分布情報（分位点や中央値）を保持させる必要がある。
ここで、ハザード比 $\text{HR} = 1.0$ の基準点（Reference value）となるBMI値を明確に指定する（今回は臨床的標準値である BMI = 22.0 にセットする）。
# データ分布情報の保持
dd <- datadist(df_rcs)
# BMIのリファレンス値（HR=1.0の基準点）を 22.0 に固定指定
dd$limits$bmi[2] <- 22.0
# オプションに登録
options(datadist = "dd")
Step 3: rms::cph + rcs() によるCox比例ハザードモデル構築
rms::cph() 関数（Cox Proportional Hazards）を用い、説明変数 bmi を rcs(bmi, 4) と指定して4つのKnotsを持つRCSモデルを適用する。年齢と性別も共変量として調整する。
※注意: 後で
Predict()やanova()を正しく動かすため、x = TRUE, y = TRUEオプションを必ず付与すること。
# RCSを組み込んだCox回帰モデルの構築
fit_cox_rcs <- cph(
Surv(time, status) ~ rcs(bmi, 4) + age + gender,
data = df_rcs,
x = TRUE,
y = TRUE
)
# モデル概要の出力
print(fit_cox_rcs)
Step 4: rms::anova による非線形性検定（$P$ for non-linearity）の抽出
モデルにおける非線形性の有意性を検証するため、anova() 関数を実行する。
# 非線形性の検定（Test for Non-linearity）
anova_res <- anova(fit_cox_rcs)
print(anova_res)
> print(anova_res)
Wald Statistics Response: Surv(time, status)
Factor Chi-Square d.f. P
bmi 273.40 3 <.0001
Nonlinear 41.43 2 <.0001
age 65.81 1 <.0001
gender 8.74 1 0.0031
TOTAL 308.91 5 <.0001
bmi(Factor全体): $\text{Chi-Square} = 273.40, P < 0.001$ $\rightarrow$ BMI全体として死亡リスクに強く関連している。Nonlinear(非線形効果): $\text{Chi-Square} = 41.43, P < 0.001$ $\rightarrow$ $P$ for non-linearity $< 0.001$ であり、明確な非線形関係（U字型）が存在することが統計的に証明された。
Step 5: Predict() ＋ ggplot2 による連続HR曲線（95% CI帯付き）の描画
rms::Predict() 関数を用いて、BMIの変化に伴う連続ハザード比とその95%信頼区間を計算する。
ref.zero = TRUE: 設定したリファレンス値（BMI = 22.0）での対数ハザードを 0 とする。fun = exp: 対数ハザード（Log-hazard）を指数変換してハザード比（Hazard Ratio）に変換する。
# 連続HRの予測値を計算
pred_hr <- Predict(fit_cox_rcs, bmi, ref.zero = TRUE, fun = exp)
# Predictオブジェクトをggplot2で美麗に描画
ggplot(pred_hr) +
geom_hline(yintercept = 1.0, linetype = "dashed", color = "gray50", size = 0.8) +
geom_vline(xintercept = 22.0, linetype = "dotted", color = "red", size = 0.8) +
theme_minimal(base_size = 14) +
labs(
x = expression(paste("Body Mass Index (kg/m"^2, ")")),
y = "Hazard Ratio for All-Cause Mortality (95% CI)",
title = "Continuous Association between BMI and Mortality Risk",
subtitle = "Restricted Cubic Spline Cox Regression (Reference: BMI = 22.0 kg/m²)"
) +
theme(
plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
描画されたグラフにより、「BMI = 22.0 を基準とした連続的なハザード比の推移」と「95%信頼区間の帯」がきれいに可視化され、低体重および肥満におけるリスク上昇が一目で理解できるようになる。
# Predictオブジェクトをデータフレームに変換
pred_df <- as.data.frame(pred_hr)
# 95%信頼区間の下限(lower)が 1.0 を超えている（＝統計学的に有意にHR>1.0である）BMI領域を抽出
pred_df %>%
filter(lower > 1.0)グラフでは、BMI 25 あたりから、95% 信頼区間の下限が 1 を超えているように見える。実際にいくつを超えると 95% 信頼区間の下限が 1 を超えるかは、上記のようにリストを出力して確認することができる。
5. 論文でそのまま使える英文テンプレート（Methods & Results）
論文の Methods および Results セクションにそのまま記述できる標準的な英文表現テンプレートを提示する。
Methods セクション英文テンプレート
“To evaluate potential non-linear associations between continuous variables (e.g., body mass index [BMI]) and all-cause mortality, restricted cubic splines (RCS) with 4 knots located at the 5th, 35th, 65th, and 95th percentiles were incorporated into Cox proportional hazards models using the
rmspackage in R.The reference value for calculating hazard ratios (HRs) was specified at BMI = 22.0 kg/m². Non-linearity was formally evaluated using the Wald test ($P$for non-linearity). All models were adjusted for age and gender.”
Results セクション英文テンプレート
“Multivariable Cox regression with restricted cubic splines revealed a significant non-linear, U-shaped association between BMI and all-cause mortality ($P$for overall association$< 0.001$,$P$for non-linearity$< 0.001$).
Using BMI = 22.0 kg/m² as the reference point ($\text{HR} = 1.0$), higher BMI levels ($> 24.6 \text{ kg/m}^2$) were significantly associated with an increased risk of all-cause mortality.”
6. まとめ＆一括コピペ用Rスクリプト
本記事のポイントおさらい
- 連続変数のカテゴリー化は悪手: 情報損失、不自然な境界リスク、査読者からの P-hacking 疑いをもたらす。
- RCSは非線形解析の正攻法: 区間を3次式で滑らかに繋ぎ、データの両端（エッジ）を直線に制限することで安定したリスクカーブを描画できる。
- 結節（Knots）は4個が基本: パーセンタイル自動配置により、研究者の主観を排除して客観的なモデルを構築する。
- $P$ for non-linearity の提示が必須:
rms::anova()で非線形性の $P$ 値を算出・提示し、グラフ（連続HRカーブ）とともに説得力を持って主張する。
一括実行用Rスクリプト
以下のコードをコピー＆ペーストすることで、疑似データ作成から datadist 設定、cph + rcs モデル構築、anova 非線形性検定、Predict + ggplot2 による95% CI帯付き連続HRカーブ出力までを一気通貫で走らせることができる。
# ==============================================================================
# Restricted Cubic Splines (RCS) による Cox回帰 連続HRカーブ描画スクリプト
# パッケージ: rms, tidyverse
# ==============================================================================
if (!requireNamespace("rms", quietly = TRUE)) install.packages("rms")
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
library(rms)
library(tidyverse)
# 1. 疑似臨床データの作成（U字型リスクデータ）
set.seed(123)
n <- 500
bmi <- runif(n, min = 18, max = 35)
age <- round(rnorm(n, mean = 65, sd = 10))
gender <- factor(sample(c("Male", "Female"), n, replace = TRUE))
log_hazard <- 0.02 * (bmi - 22)^2 + 0.03 * (age - 65) + 0.3 * (gender == "Male") - 1.5
lambda <- exp(log_hazard)
time <- rexp(n, rate = lambda * 0.05)
status <- if_else(time > 36, 0, 1)
time <- pmin(time, 36)
df_rcs <- tibble(
id = 1:n,
time = round(time, 1),
status = status,
bmi = round(bmi, 1),
age = age,
gender = gender
)
# 2. データ分布情報の設定とリファレンス値の固定（BMI = 22.0）
dd <- datadist(df_rcs)
dd$limits$bmi[2] <- 22.0
options(datadist = "dd")
# 3. rms::cph + rcs によるCoxモデル構築（Knots = 4）
fit_cox_rcs <- cph(
Surv(time, status) ~ rcs(bmi, 4) + age + gender,
data = df_rcs,
x = TRUE,
y = TRUE
)
# 4. 非線形性の検定（P for non-linearity）の抽出
cat("\n--- 非線形性の検定結果 (rms::anova) ---\n")
print(anova(fit_cox_rcs))
# 5. 連続HR予測値の算出と描画
pred_hr <- Predict(fit_cox_rcs, bmi, ref.zero = TRUE, fun = exp)
p_rcs <- ggplot(pred_hr) +
geom_hline(yintercept = 1.0, linetype = "dashed", color = "gray50", size = 0.8) +
geom_vline(xintercept = 22.0, linetype = "dotted", color = "red", size = 0.8) +
theme_minimal(base_size = 14) +
labs(
x = expression(paste("Body Mass Index (kg/m"^2, ")")),
y = "Hazard Ratio for All-Cause Mortality (95% CI)",
title = "Continuous Association between BMI and Mortality Risk",
subtitle = "Restricted Cubic Spline Cox Regression (Reference: BMI = 22.0 kg/m²)"
) +
theme(
plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
# グラフの表示
print(p_rcs)
# 統計学的に有意にHR>1.0であるBMI領域を抽出
# Predictオブジェクトをデータフレームに変換
pred_df <- as.data.frame(pred_hr)
# 95%信頼区間の下限(lower)が 1.0 を超えている（＝統計学的に有意にHR>1.0である）BMI領域を抽出
pred_df %>%
filter(lower > 1.0)
コメント