生活習慣と収縮期血圧の関連を線形回帰で調べたところ、年齢、BMI、食塩摂取量の係数はいずれも想定した向きに推定されたのに、決定係数は0.15だった、ということがあります。この0.15を見て説明変数を追加する判断につながることがありますが、そこで変数を足せばよいとは限りません。
決定係数は、結果変数の変動のうち当てはめた回帰式で説明される割合を表しますが、その値は説明変数と結果変数の関係だけでなく、説明変数の数や対象集団のばらつきによっても動きます。そのため、0.15という数字だけからモデルの良し悪しを判断することはできません。
ここでは、最小二乗法で当てはめた切片を含む線形回帰を対象にします。乱数で生成した血圧データを使い、決定係数が変数追加、対象集団の範囲、標本内評価と交差検証でどのように動くかを確かめ、予測性能を評価したい場合に決定係数と何を一緒に見るべきかまで考えます。
決定係数の0.15は何を表しているか
回帰式そのものの意味と、単回帰と重回帰で同じ係数が異なる比較条件を表すことについては、回帰分析の目的と単回帰分析と重回帰分析の違いで扱っています。
対象者数を$n$、$i$番目の対象者で観測された結果変数を$y_i$、回帰式から得られた予測値を$\hat{y}_i$、結果変数の平均を$\bar{y}$とします。切片を除く説明変数の数は$p$とします。
切片を含む通常の最小二乗回帰では、決定係数は次のように書けます。
$$R^2 = 1 – \frac{\sum_{i=1}^{n}(y_i – \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i – \bar{y})^2} \tag{1}$$
式(1)の分母は観測値が平均からどれだけ散らばっているかを表す全平方和、分子は回帰式を当てはめた後にも残った残差平方和です。したがって、$R^2$は全平方和に対して残差平方和がどれだけ減ったかを割合で表しています。
乱数で生成した240例について、収縮期血圧を年齢だけで説明すると$R^2$は0.0809でしたが、BMIと食塩摂取量を加えると0.1512となり、調整済み決定係数は0.1404、残差標準誤差は11.113mmHgでした。データ生成時に誤差の標準偏差を11mmHgと置いているため、残差のばらつきは設定した誤差の大きさに近い値へ戻っています。
そして、ここで得られた11.113mmHgは、この標本に回帰式を当てはめた後に残った誤差の標準偏差を推定した値で、新しい患者に対する予測誤差そのものではありません。ただし、結果変数と同じmmHgで表されるため、モデルを当てはめた後にどの程度のばらつきが残っているのかは、$R^2=0.1512$より直感的に読み取れます。
Kvålseth（1985）は、$R^2$には複数の定義があり、特に切片を含まないモデルや通常の線形回帰以外では値の性質が変わることを指摘しています。今回のような切片を含む最小二乗回帰では式(1)で問題ありませんが、「1に近いほど常に良いモデル」と読み替えると意味が変わります。
変数を足すと決定係数は下がらない
年齢、BMI、食塩摂取量の3変数に、血圧とは無関係に生成した乱数変数を追加し、60例だけを取り出すと3変数モデルの$R^2$は0.1567でした。乱数変数を5個追加すると0.2193、10個で0.2682、20個で0.4411、30個では0.5550まで上がりました。
この上昇は最小二乗法の仕組みから起こり、既存の説明変数を残したまま新しい変数を加えると元の係数をそのまま使う選択肢も残るため、残差平方和は増えません。そのため$R^2$は下がらず、有限標本では結果と無関係な変数を追加した場合でも通常は少しずつ上がります。
たとえば、真のデータ生成式に入れた説明変数は3個ですが、$R^2$だけを基準にすると33個の説明変数を入れたモデルの方が高く評価されます。
この変数数に対する増加を補正する指標が、調整済み決定係数です。
$$R_{\mathrm{adj}}^2 = 1 – (1 – R^2)\,\frac{n – 1}{n – p – 1} \tag{2}$$
式(2)では、説明変数が増えるほど残差自由度$n-p-1$が小さくなります。そのため、新しい変数を加えて得られた残差平方和の減少が変数1個分の自由度に見合わなければ、調整済み決定係数は下がります。
前記の60例の3変数モデルでは調整済み決定係数が0.1115でしたが、乱数変数を30個追加すると通常の$R^2$は0.5550まで上がった一方、調整済み決定係数は-0.0097でした。この負の値は式(2)から生じるもので、$R^2<p/(n-1)$なら調整済み決定係数は0を下回ります。
ただし、調整済み決定係数を使えば過剰適合まで補正できると読むと、その役割を広く捉えすぎます。式(2)が反映するのは標本サイズと説明変数の数であり、同じデータを見ながら候補変数を選んだことによる楽観性や、モデル選択を何度も試した影響まで評価しているわけではありません。
Babyak（2004）が問題にしているのも、このような標本への過剰な適合で、候補変数が多い状況やデータ駆動型の変数選択では、標本内で得られた性能が新しい標本でも再現されるとは限りません。私は、調整済み決定係数を「過剰適合を解決した$R^2$」とは扱わず、説明変数数に対する自由度補正として読みます。
同じ傾きでも集団が変われば値は動く
今度はモデルを変えず、対象集団だけを変えます。収縮期血圧に対する年齢の傾きを0.45mmHg/歳、誤差の標準偏差を11mmHgに固定し、年齢のばらつきだけが異なる300例の集団を2つ生成しました。
ここで、年齢の標準偏差が12.891歳だった集団では$R^2=0.2744$、標準偏差2.994歳の集団では$R^2=0.0325$となり、8倍を超える差が生じました。一方、残差標準誤差は11.293mmHgと10.693mmHgで、どちらも生成時に設定した11mmHg付近でした。
そして、単回帰モデル$Y=\alpha+\beta X+\varepsilon$で、$X$と誤差が独立し、誤差分散を$\sigma^2$とすると、母集団での決定係数に対応する量は次のように書けます。
$$\rho^2 = \frac{\beta^2\,\mathrm{Var}(X)}{\beta^2\,\mathrm{Var}(X) + \sigma^2} \tag{3}$$
式(3)を見ると、傾き$\beta$と誤差分散$\sigma^2$が同じでも、$\mathrm{Var}(X)$が変われば$\rho^2$も変わることが分かります。傾き0.45、誤差標準偏差11のまま、年齢の標準偏差を12とすると理論値は0.1942、標準偏差3なら0.0148です。
この例からも、年齢分布、疾患の重症度、施設の選択基準などが変われば、同じ回帰関係でも結果変数と説明変数の分散が変わるため、別々の研究に報告された$R^2$を、そのままモデルの優劣として並べる方法が推奨されないことが分かります。
絶対的な誤差の大きさを比較したいのであれば、同じアウトカム定義と測定尺度を使っていることを確認したうえで、残差標準誤差や平均絶対誤差など、結果変数と同じ単位を持つ指標も一緒に見る方が解釈しやすくなります。
見かけの0.44が交差検証では負になった
60例へ乱数変数を20個追加した23変数モデルでは、当てはめに使った同じ60例で計算した$R^2$が0.4411、調整済み決定係数が0.0840でした。ここまでは全て標本内の評価です。
同じモデル構築手順が、当てはめに使っていない対象者でもどの程度の予測性能を持つかを見るため、一つ抜き交差検証を行い、1例を除いた59例で回帰式を当てはめて除いた1例を予測する操作を60回繰り返します。その予測誤差平方和をPRESSとすると、ここでは次の量を計算しています。
$$R_{\mathrm{CV}}^2 = 1 – \frac{\mathrm{PRESS}}{\sum_{i=1}^{n}(y_i-\bar{y})^2} \tag{4}$$
式(4)は、交差検証で得られた予測残差平方和を、元の標本における平均からの全平方和と比較したPRESS型の交差検証指標です。式(1)の標本内$R^2$と同じ計算対象ではないため、両者は区別して表記した方が安全です。
また、23変数モデルでは、見かけの$R^2$が0.4411だったのに対し、$R_{\mathrm{CV}}^2$は-0.5431でした。この負の値は、交差検証で得られたPRESSが、この標本の平均$\bar{y}$を基準にした全平方和より大きかったことを表します。
したがって、「44%を説明できたモデルが、新しい対象者でも44%を説明できる」とは読めません。今回の23変数モデルでは、標本内評価と交差検証による予測評価の差が非常に大きくなりました。
そして、真の生成式に入れた3変数だけのモデルでも、見かけの$R^2$は0.1567、$R_{\mathrm{CV}}^2$は0.0134でした。有限標本ではモデルを当てはめたデータでの評価が楽観的になりやすく、今回のデータでは不要な変数を多く加えた23変数モデルで標本内評価との差が特に大きくなっています。
Harrell、Lee、Mark（1996）は、回帰モデルの予測性能をモデル開発に使った同じ標本だけで評価すると楽観的になりうることを論じ、ブートストラップや交差検証による内部妥当化を扱っています。予測モデルの記事で$R^2$を見る場合、私は数値の大きさより先に、開発標本上の見かけの値なのか、内部妥当化を経た値なのかを確認します。
論文では決定係数だけで性能を語らない
$R^2$は、結果変数の変動に対して回帰式がどの程度の説明力を持っていたかを1つの割合で示せるため、報告する意味があります。ただし、研究目的が予測なら、$R^2$だけでモデル性能を説明するのは足りません。
予測モデルの報告を扱うTRIPODの解説論文で、Moonsら（2015）は予測性能の主要な側面として較正と識別を挙げ、較正は予測した値と実際に観測された値がどの程度一致しているかを表すため、可能なら図で確認します。また、モデル全体の性能を表す指標として$R^2$などの説明された変動も扱われており、$R^2$は性能評価の一部であって単独で完結する指標ではありません。
線形回帰以外で「決定係数」と呼ばれる値を見るときは定義の確認が必要で、ロジスティック回帰では式(1)の平方和分解がそのまま成立しません。Nagelkerke（1991）は尤度に基づく$R^2$の一般化を扱い、Cox-Snell型の指標が離散アウトカムでは1に到達しないことから上限が1になるよう再尺度化した指標を提案し、CameronとWindmeijer（1997）も非線形回帰モデルにおける適合度指標を論じています。
そのため、名前が同じでも、何を分子と分母に置いているかが異なるため、線形回帰の$R^2=0.15$とロジスティック回帰のNagelkerke $R^2=0.15$を、同じ尺度の15%として直接比較することはできません。
もう一つ分けておきたいのが因果解釈で、高い$R^2$はモデルが観測された結果変数の変動をよく説明していることを表しますが、説明変数の係数に因果効果としての意味を与える根拠にはなりません。因果効果として解釈するには時間的な順序、交絡、選択、測定などの仮定を別に考える必要があり、$R^2$はその判断をする指標ではありません。
実際の解析で私が$R^2$を見るときは「高いか低いか」から入らず、研究目的が関連の記述なのか予測なのかを確認し、予測なら絶対誤差や内部妥当化後の性能まで見ます。説明変数を追加する場合も、$R^2$を上げるためではなく、研究目的とモデルに必要な変数だから入れるという順序にします。
対象集団の範囲についても、年齢を狭い範囲に限定すれば式(3)から予想されるように年齢が説明する変動の割合は小さくなります。正しく指定された線形モデルで、選択が説明変数の範囲だけに基づき誤差との関係を変えない条件なら、範囲を狭めること自体が傾きの期待値を変えるわけではありませんが、年齢の情報量が減るため係数の標準誤差は大きくなり、今回の生成データでも0.0507から0.2066へ増えました。
Rで3つの動きを確かめる
ここまでの数値は、base Rとstatsだけで再現できます。最初に、年齢、BMI、食塩摂取量から収縮期血圧を生成した240例を用意します。
set.seed(42)
n_main <- 240
age <- rnorm(n_main, mean = 62, sd = 9)
bmi <- rnorm(n_main, mean = 24.5, sd = 3.2)
salt <- rnorm(n_main, mean = 9.5, sd = 2.0)
sbp <- 92 + 0.45 * age + 0.90 * bmi + 1.60 * salt +
rnorm(n_main, mean = 0, sd = 11)
main_data <- data.frame(sbp, age, bmi, salt)
fit_age <- lm(sbp ~ age, data = main_data); fit_full <- lm(sbp ~ age + bmi + salt, data = main_data)
round(c(summary(fit_age)[["r.squared"]], summary(fit_full)[["r.squared"]],
summary(fit_full)[["adj.r.squared"]], summary(fit_full)[["sigma"]]), 4)#> [1] 0.0809 0.1512 0.1404 11.1130左から、年齢だけの$R^2$、3変数モデルの$R^2$、調整済み決定係数、残差標準誤差です。0.1512と11.113は同じ「当てはまり」を別の尺度で言い換えた値ではなく、前者は全変動に対する割合、後者は回帰後に残ったばらつきをmmHgで表した値です。
60例へ不要な変数を追加した図1は、同じ乱数変数30個を使い、モデルへ入れる個数だけを変えて作れます。
small_data <- main_data[1:60, ]; noise <- matrix(rnorm(60 * 30), nrow = 60)
colnames(noise) <- paste0("noise", 1:30); small_data <- cbind(small_data, as.data.frame(noise))
k <- c(0, 5, 10, 20, 30)
r2_table <- t(sapply(k, function(j) {
vars <- c("age", "bmi", "salt", if (j > 0) paste0("noise", 1:j))
fit <- lm(reformulate(vars, response = "sbp"), data = small_data)
c(R2 = summary(fit)[["r.squared"]], Adjusted = summary(fit)[["adj.r.squared"]])
}))
matplot(k, r2_table, type = "b", xlab = "Number of noise variables", ylab = "R-squared")
round(r2_table[c(1, 5), ], 4)#> R2 Adjusted
#> [1,] 0.1567 0.1115
#> [2,] 0.5550 -0.0097同じ候補モデルに対して、通常の$R^2$と調整済み決定係数が異なる方向へ動いています。ここから言えるのは、調整済み決定係数の方が常に正しいということではなく、通常の$R^2$には変数追加への罰則がないということです。
続いて、傾きと誤差の大きさを固定し、年齢の範囲だけを変えます。図2の2つのパネルもこのデータから作れます。
make_cohort <- function(age_sd) {
age_c <- rnorm(300, mean = 62, sd = age_sd)
data.frame(age = age_c, sbp = 92 + 0.45 * age_c + rnorm(300, mean = 0, sd = 11))
}
wide <- make_cohort(12); narrow <- make_cohort(3)
fit_wide <- lm(sbp ~ age, data = wide); fit_narrow <- lm(sbp ~ age, data = narrow)
par(mfrow = c(1, 2)); plot(sbp ~ age, data = wide, xlab = "Age (years)", ylab = "SBP (mmHg)")
abline(fit_wide); plot(sbp ~ age, data = narrow, xlab = "Age (years)", ylab = "SBP (mmHg)"); abline(fit_narrow)
round(c(summary(fit_wide)[["r.squared"]], summary(fit_narrow)[["r.squared"]]), 4)#> [1] 0.2744 0.0325同じ傾き0.45、同じ誤差標準偏差11を置いていても、年齢の広がりを変えるだけで$R^2$は0.2744と0.0325に分かれました。そのため、$R^2$の違いをすぐにモデルの質の違いへ読み替えることは適切ではありません。
最後に、23変数モデルについて、標本内$R^2$とPRESS型の$R_{\mathrm{CV}}^2$を比べます。
noisy_data <- small_data[, 1:24]; fit_noisy <- lm(sbp ~ ., data = noisy_data)
apparent_pred <- fitted(fit_noisy)
loo_pred <- noisy_data[["sbp"]] - residuals(fit_noisy) / (1 - hatvalues(fit_noisy))
press <- sum((noisy_data[["sbp"]] - loo_pred)^2)
tss_small <- sum((noisy_data[["sbp"]] - mean(noisy_data[["sbp"]]))^2)
cv_r2 <- 1 - press / tss_small
par(mfrow = c(1, 2)); plot(noisy_data[["sbp"]], apparent_pred, xlab = "Observed", ylab = "Predicted")
abline(0, 1); plot(noisy_data[["sbp"]], loo_pred, xlab = "Observed", ylab = "LOOCV prediction"); abline(0, 1)
round(c(summary(fit_noisy)[["r.squared"]], cv_r2), 4)#> [1] 0.4411 -0.54310.4411と-0.5431は、同じデータから出た数字でも評価している場面が異なり、前者はモデルを当てはめた標本上の適合度、後者は1例ずつ除外して得た予測誤差を使った内部妥当化の指標です。予測を目的とする研究なら、この区別を方法と結果の両方にはっきり書くことで、モデル性能の読み違いを避けやすくなります。
$R^2$、調整済み決定係数、予測誤差のどれを報告するかは、研究目的が関連の推定なのか予測なのか、説明変数を事前に決めたのかデータから選んだのか、内部妥当化を行ったのかによって異なります。また、係数の標準誤差については、ロバスト標準誤差で扱った分散の仮定も別に確認する必要があり、Dr.データサイエンスでは、こうした解析方針と報告内容のご相談を承っています。
参考文献
- Kvålseth TO. Cautionary Note about R². The American Statistician. 1985;39(4):279-285.
- Nagelkerke NJD. A note on a general definition of the coefficient of determination. Biometrika. 1991;78(3):691-692.
- Cameron AC, Windmeijer FAG. An R-squared measure of goodness of fit for some common nonlinear regression models. Journal of Econometrics. 1997;77(2):329-342.
- Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine. 1996;15(4):361-387.
- Babyak MA. What you see may not be what you get: a brief, nontechnical introduction to overfitting in regression-type models. Psychosomatic Medicine. 2004;66(3):411-421.
- Moons KGM, Altman DG, Reitsma JB, Ioannidis JPA, Macaskill P, Steyerberg EW, Vickers AJ, Ransohoff DF, Collins GS. Transparent Reporting of a multivariable prediction model for Individual Prognosis or Diagnosis (TRIPOD): explanation and elaboration. Annals of Internal Medicine. 2015;162(1):W1-W73.
