Pythonのscikit-learn に用意されている Diabetes dataset を用いて、
BMIから疾患進行度を予測する
というテーマで、単回帰モデルを作ることを通して単回帰分析を勉強する。単回帰分析では、1つの説明変数 $x$ から目的変数 $y$を予測するモデルを考える。今回の例では、
- 説明変数 $x$：BMI
- 目的変数 $y$：疾患進行度
とします。
※あくまで統計学の勉強のためであり、本結果から医学的内容を示唆しているわけではありません。
1. 単回帰分析の手順
単回帰モデルでは、次のような直線を考えます。
y = a+bx
-
データをグラフで確認する
-
$y$の平均値$\bar{y}$を求める
-
$x$の平均値$\bar{x}$を求める
-
yの個々のデータから$\bar{y}$を引いた値$(y-\bar{y})$を求める
-
xの個々のデータから$\bar{x}$を引いた値$(x-\bar{x})$を求める
-
4で求めた値の2乗の値$(y-\bar{y})^2$を求める
-
5で求めた値の2乗の値$(x-\bar{x})^2$を求める
-
4と5を乗算した値$(y-\bar{y})(x-\bar{x})$を求める
-
$y$の偏差平方和$S_{yy}$(6の合計値)を求める
-
$x$の偏差平方和$S_{xx}$(7の合計値)を求める
-
偏差積和$S_{xy}$(8の合計値)を求める
-
回帰係数$b$を求める
-
切片$a$を求める
-
回帰直線を求める
-
回帰係数$b$の検定を行う
-
回帰係数$b$の信頼区間を求める
-
決定係数$R^2$を求め、モデルの当てはまりを評価する
-
予測値の95%信頼区間を求める
-
予測値の95%予測区間を求める
回帰係数$b$は以下で求めることができる。
\boxed{b=\frac{S_{xy}}{S_{xx}}}
上記を用いて切片$a$は以下で求めることができる。
\boxed{a=\bar{y}-b\bar{x}}
求めた回帰直線$y = a+bx$は、「平均値($\bar{x}, \bar{y}$)を通る」性質がある。
回帰係数の検定
求めた回帰係数$b$が、偶然のばらつきによって得られたものではないかを検定する。次の仮説を設定する。
H_0:b=0
H_1:b\neq0
帰無仮説$H_0$は、今回の例では
「BMIと疾患進行度には線形な関係がない」
ことを意味する。
回帰係数$b$の標準誤差を$SE(b)$とすると、検定統計量$t$は自由度n-2のt分布に従い、
\boxed{
t=\frac{b}{SE(b)} \bigg( = \frac{求めた回帰係数の値}{標準誤差} \bigg)
}
で求めることができる。この$t$値から$p$値を求め、あらかじめ設定した有意水準$\alpha$と比較する。
標準誤差$SE(b)$は以下の過程で計算する。まず、残差平方和を$SSE$とすると、回帰直線$y = a+bx$による推定値$\hat{y_i}$を用いて以下のように表される。
SSE = \sum_{i=1}^n (y_i - \hat{y_i})^2
自由度$n-2$で割って誤差分散$s^2$とすると、
s^2 = \frac{1}{n-2}SSE
標準誤差$SE(b)$は、誤差分散$s^2$の平方根より、
SE(b) = \sqrt{s^2} = \sqrt{\frac{1}{n-2}\sum_{i=1}^n (y_i - \hat{y_i})^2}
回帰係数の信頼区間
回帰係数$b$の95%信頼区間は、自由度$n-2$のt分布を用いて以下で求めることができる。
\boxed{
b\pm t_{0.975,n-2}SE(b)
}
95%信頼区間に$b=0$が含まれていなければ、両側5%の検定においても統計的に有意と判断できる。
決定係数
直線との当てはまりを評価する決定係数$R^2$は以下で求めることができる。
\boxed{R^2=\frac{{S_{xy}}^2}{S_{xx}S_{yy}}}
これは相関係数$r$の2乗としても表すことができ、相関係数は以下式である。
r = \frac{S_{xy}}{\sqrt{S_{xx}S_{yy}}}
また、決定係数$R^2$は「回帰によって説明できた変動 ÷ 全体の変動」という意味を表しており、データがすべて回帰直線で説明できたら1となる。よって、回帰で説明できない変動から表現すると
\boxed{R^2=1-\frac{SSE}{S_{yy}}}
となる。
予測値の信頼区間・予測区間
回帰モデルの$x$に値$x_0$を入力すれば予測値$\hat{y_0}$が計算されるが、回帰モデルを扱う場合は、その予測の不確実性もセットで確認しなければならない。
\hat{y_0}=a+b x_0
予測値$\hat{y_0}$の周りには、
- 平均応答の95%信頼区間
- 個々の観測値に対する95%予測区間
平均応答の信頼区間は、例えば$x_0=20$とすると以下を推定するための区間となる。
「BMIが20の集団における平均的な疾患進行度はどの程度か」
一方、予測区間は以下を予測するための区間となる。
「BMIが20の個人の疾患進行度はどの程度になるか」
個々のデータにはばらつきがあるため、一般に予測区間の方が信頼区間よりも広くなる。
平均応答の標準誤差は、
\boxed{SE(\hat{y_0})=s\sqrt{\frac{1}{n}+\frac{(x_0-\bar{x})^2}{S_{xx}}}}
で計算することができ、信頼区間は以下で表すことができる。
\boxed{
\hat{y_0}\pm t_{0.975,n-2}SE(\hat{y_0})
}
個々の予測値$\hat{y_0}$に対する標準誤差は、
\boxed{SE_{prediction}=s\sqrt{1+\frac{1}{n}+\frac{(x_0-\bar{x})^2}{S_{xx}}}}
で計算することができ、予測区間は以下で表すことができる。
\boxed{
\hat{y_0}\pm t_{0.975,n-2}SE_{prediction}
}
2. Pythonプログラムによる単回帰分析(Part I)
この章では、以上の分析の流れをPythonプログラムにより実装し、順次内容を確認していきながら理解を進めていく。
必要なライブラリ読み込み
今回の分析では、Pandas、NumPy、SciPy、Matplotlib、Seabornを利用する。
import pandas as pd
import numpy as np
from scipy import stats
from matplotlib import pyplot as plt
import seaborn as sns
np.set_printoptions(precision=3)
pd.set_option('display.precision', 3)
sns.set()
データセットの準備
まず、scikit-learn の load_diabetes() を使ってデータを読み込む。デフォルトでは説明変数が標準化された状態となっているため、BMIの値を直感的に理解できるようscaled=False を指定して標準化される前の値としてデータを用意する。
from sklearn.datasets import load_diabetes
import pandas as pd
diabetes = load_diabetes(as_frame=True, scaled=False)
# 単回帰分析で使う BMI と疾患進行度だけを取り出す
data = diabetes.frame.rename(columns={
'bmi': 'BMI',
'target': 'disease_progression'
})[['BMI', 'disease_progression']]
data.head()
実行結果は、最初の5件が表示される。
BMI disease_progression
--- ------ ---------------------
0 32.1 151.0
1 21.6 75.0
2 30.5 141.0
3 25.3 206.0
4 23.0 135.0
なお、disease_progressionは糖尿病患者について測定された1年後の疾患進行度を表す定量的な指標です。医学的な重症度そのものを表す値ではありません。
データをグラフで確認する
BMIと疾患進行度の関係を散布図で確認する。
sns.scatterplot(
x='BMI',
y='disease_progression',
data=data,
color='black'
)
BMIが高くなるほど、疾患進行度も高くなる傾向が見られるが、かなりばらつきもみられる。そこで、
「BMIから疾患進行度をどの程度説明・予測できるのか？」
を、単回帰モデルを使って定量的に確認する。
x と y の平均値を求める
説明変数$x$と目的変数$y$を以下として変数を用意し求める。
-
x：BMI -
y：疾患進行度 disease_progression
x = data['BMI'].tolist()
y = data['disease_progression'].tolist()
n = len(y)
y_bar = sum(y) / n
x_bar = sum(x) / n
print('データ数:', len(x))
print('BMI:', x[:5])
print('疾患進行度:', y[:5])
print('xの平均値:', x_bar)
print('yの平均値:', y_bar)
データ数: 442
BMI: [32.1, 21.6, 30.5, 25.3, 23.0]
疾患進行度: [151.0, 75.0, 141.0, 206.0, 135.0]
xの平均値: 26.37579185520362
yの平均値: 152.13348416289594
x と y の偏差平方和Sxx, Syy と偏差積和Sxyを求める
まず、$x$と$y$をそれぞれの平均との差をとり偏差をx_deviationとy_deviationに格納する。そして、それらを2乗した偏差平方と、乗算した偏差積を求める。最後に各種それぞれ合計を計算することで偏差平方和$S_{xx}$, $S_{yy}$ と偏差積和$S_{xy}$を求める。
x_deviation = [xi - x_bar for xi in x] # x偏差
y_deviation = [yi - y_bar for yi in y] # y偏差
x_deviation_squared = [xv**2 for xv in x_deviation] # x偏差平方
y_deviation_squared = [yv**2 for yv in y_deviation] # y偏差平方
xy_deviation_product = [xv*yv for xv, yv in zip(x_deviation, y_deviation)] # xとyの偏差積
Sxx = sum(x_deviation_squared)
Syy = sum(y_deviation_squared)
Sxy = sum(xy_deviation_product)
print('xの偏差:', x_deviation[:5])
print('yの偏差:', y_deviation[:5])
print('xの偏差平方:', x_deviation_squared[:5])
print('yの偏差平方:', y_deviation_squared[:5])
print('xとyの偏差積:', xy_deviation_product[:5])
print('Sxx:', Sxx)
print('Syy:', Syy)
print('Sxy:', Sxy)
xの偏差: [5.724208144796382, -4.775791855203618, 4.1242081447963805, -1.0757918552036188, -3.3757918552036195]
yの偏差: [-1.1334841628959396, -77.13348416289594, -11.13348416289594, 53.86651583710406, -17.13348416289594]
xの偏差平方: [32.766558884953234, 22.808187844229217, 17.009092821604803, 1.1573281157224438, 11.395970649659095]
yの偏差平方: [1.284786347535909, 5949.574379107718, 123.9544696054547, 2901.6015284289824, 293.55627956020595]
xとyの偏差積: [-6.488299277246647, 368.3734654286357, -45.916806064576946, -57.94915900575329, 57.83907628836432]
Sxx: 8608.230972850679
Syy: 2621009.124434389
Sxy: 88089.12828054298
回帰係数b、切片a、回帰直線を求める
b = Sxy / Sxx
a = y_bar - b * x_bar
print('切片 a:', a)
print('回帰係数 b:', b)
print(f'回帰直線: y = {a:.3f} + {b:.3f}x')
切片 a: -117.77336656656522
回帰係数 b: 10.233127870100773
回帰直線: y = -117.773 + 10.233x
散布図に回帰直線を描画してみる。
# 散布図
sns.scatterplot(x='BMI',y='disease_progression',data=data,color='black')
# 回帰直線を描くためのBMIの範囲
x_line = np.linspace(min(x), max(x), 100)
# 推定した回帰式から予測値を計算
y_line = a + b * x_line
# 回帰直線
plt.plot(x_line,y_line,color='red')
plt.xlabel('BMI')
plt.ylabel('disease_progression')
黒い点がデータセットの値、赤い直線が回帰直線である。データ全体の傾向を最もよく表す直線となっている。
回帰係数bの検定を行う
「b の標準誤差 → 検定統計量t値 → p値 → 優位水準との判定」の流れで進める。
# 残差平方和 SSE
y_hat = [a + b * xi for xi in x]
residual = [yi - yhi for yi, yhi in zip(y, y_hat)]
residual_squared = [e**2 for e in residual]
SSE = sum(residual_squared)
# 誤差分散
df = n - 2
s2 = SSE / df
# 回帰係数 b の標準誤差
SE_b = (s2 / Sxx) ** 0.5
# 検定統計量t値
t_value = b / SE_b
# 両側検定のp値
p_value = 2 * (1 - stats.t.cdf(abs(t_value), df=df))
print('残差平方和 SSE:', SSE)
print('自由度:', df)
print('誤差分散 s²:', s2)
print('回帰係数 b:', b)
print('標準誤差 SE(b):', SE_b)
print('t値:', t_value)
print('p値:', p_value)
stats.t.cdfは$-\infty$から引数で指定した値までの面積を計算するため、1 - stats.t.cdfは上側の片側面積を計算していることになる。そのため2倍することで両側検定としてのp値を求めている。
実行結果は以下のようになり、
残差平方和 SSE: 1719581.8107738823
自由度: 440
誤差分散 s²: 3908.1404790315505
回帰係数 b: 10.233127870100773
標準誤差 SE(b): 0.6737955329480582
t値: 15.187289570365301
p値: 0.0
優位水準5% とすると、p<0.05が成立するため帰無仮説($H_0:b=0$)を棄却する。
回帰係数bの信頼区間を求める
# 有意水準
alpha = 0.05
# 95%信頼区間に使用するt値
t_critical = stats.t.ppf((1-alpha/2),df=df)
# 回帰係数 b の95%信頼区間
b_ci_lower = b - t_critical * SE_b
b_ci_upper = b + t_critical * SE_b
print('95%信頼区間の臨界値:', t_critical)
print('回帰係数 b:', b)
print(f'95%信頼区間: [{b_ci_lower:.3f}, {b_ci_upper:.3f}]')
stats.t.ppfも累積の確率値に対応するt値が求まる。両側検定の場合上側棄却域は$\alpha/2$となる。あるいは、引数に$\alpha/2$(下側の棄却域)を設定して求めたt値の絶対値で求めてもよい。
実行結果は以下となった。
95%信頼区間の臨界値: 1.9653701191315147
回帰係数 b: 10.233127870100773
95%信頼区間: [8.909, 11.557]
決定係数を求め、モデルの当てはまりを評価する
残差平方和から求めた場合と値が一致するかも一緒に確認する。
R2_1 = (Sxy ** 2) / (Sxx * Syy)
R2_2 = 1 - SSE/Syy
print('決定係数 R^2:', R2_1)
print('残差から計算した決定係数 R^2:', R2_2)
決定係数 R^2: 0.3439237602253802
残差から計算した決定係数 R^2: 0.3439237602253803
散布図を確認した際、データにばらつきがあることが分かっていたため、回帰直線との当てはまりはそこまで高くない数値である。
以上の結果から、BMIと疾患進行度には統計的に有意な関係がある一方で、決定係数$R^2$が0.344程度であるため、BMIだけですべての疾患進行度のばらつきを説明できるわけではないということが判断できる。つまり、
「統計的に有意」＝「予測性能が高い」ではない
これらのモデルを評価・理解したうえで次に予測を議論する。
予測値の95%信頼区間、予測区間を求める
x_new = 20 # 予測したいBMI
y_pred = a + b * x_new # 回帰モデルによる予測値
# BMI = 20 における平均応答の標準誤差
s = s2**0.5
SE_mean = s * ((1/n) + ((x_new-x_bar)**2/Sxx)) ** 0.5
# 95%信頼区間
mean_ci_lower = y_pred - t_critical * SE_mean
mean_ci_upper = y_pred + t_critical * SE_mean
# BMI = 20 における個々の観測値の予測標準誤差
SE_prediction = s * (1 + (1/n) + ((x_new-x_bar)**2/Sxx)) ** 0.5
# 95%予測区間
prediction_lower = y_pred - t_critical * SE_prediction
prediction_upper = y_pred + t_critical * SE_prediction
print(f'BMI = {x_new} の予測値: {y_pred:.3f}')
print(f'平均応答の標準誤差: {SE_mean:.3f}')
print(f'95%信頼区間: [{mean_ci_lower:.3f}, {mean_ci_upper:.3f}]')
print(f'95%予測区間: [{prediction_lower:.3f}, {prediction_upper:.3f}]')
BMI = 20 の予測値: 86.889
平均応答の標準誤差: 5.225
95%信頼区間: [76.621, 97.158]
95%予測区間: [-36.405, 210.183]
分散・共分散行列を計算する関数を使用する場合
# 分散共分散を求める。(1,1)=xの分散、(2,2)=yの分散、(1,2)=(2,1)=x,yの共分散
cov_mat = np.cov(x, y, ddof=0)
cov_mat
array([[ 19.476, 199.297],
[ 199.297, 5929.885]])
以下のように処理を簡略化できる。
x_bar = np.mean(x)
y_bar = np.mean(y)
Sxx = cov_mat[0,0]
Syy = cov_mat[1,1]
Sxy = cov_mat[0,1]
b = Sxy / Sxx
a = y_bar - b * x_bar
print('切片 a:', a)
print('回帰係数 b:', b)
print(f'回帰直線: y = {a:.3f} + {b:.3f}x')
切片 a: -117.7733665665651
回帰係数 b: 10.23312787010077
回帰直線: y = -117.773 + 10.233x
関数を使用しないで計算した結果と一致している。
3. Pythonプログラムによる単回帰分析(Part II)
statsmodelsを利用すると簡単に回帰モデルを作ることができる。
import statsmodels.formula.api as smf
# statsmodelsによる回帰分析
lm_model = smf.ols(formula='disease_progression ~ BMI',data=data).fit()
print(lm_model.summary())
以下の出力結果の画面中段より、それぞれ左から「推定された回帰係数」「標準誤差」「帰無仮説:回帰係数=0とした時のt値」「p値」「下限信頼限界と上限信頼限界」を示している。
t値は係数の値と標準誤差の除算で得られる。
決定係数 R-squared で表示され、このモデルはデータをどの程度説明できているのかを評価することができる。
OLS Regression Results
===============================================================================
Dep. Variable: disease_progression R-squared: 0.344
Model: OLS Adj. R-squared: 0.342
Method: Least Squares F-statistic: 230.7
Date: Fri, 14 Aug 2026 Prob (F-statistic): 3.47e-42
Time: 02:12:21 Log-Likelihood: -2454.0
No. Observations: 442 AIC: 4912.
Df Residuals: 440 BIC: 4920.
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -117.7734 18.019 -6.536 0.000 -153.187 -82.359
BMI 10.2331 0.674 15.187 0.000 8.909 11.557
==============================================================================
Omnibus: 11.674 Durbin-Watson: 1.848
Prob(Omnibus): 0.003 Jarque-Bera (JB): 7.310
Skew: 0.156 Prob(JB): 0.0259
Kurtosis: 2.453 Cond. No. 162.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Part I で定義にそって計算した値とstatsmodelsによる値を確認すると一致していることが分かり、定義の理解の正しさおよびstatsmodelsライブラリの出力結果を理解することができる。
print('--- 定義式による計算 ---')
print(f'b : {b:.6f}')
print(f'SE(b) : {SE_b:.6f}')
print(f't値 : {t_value:.6f}')
print(f'p値 : {p_value:.6e}')
print(f'95% CI : [{b_ci_lower:.6f}, {b_ci_upper:.6f}]')
print('\n--- statsmodels ---')
print(f'b : {lm_model.params["BMI"]:.6f}')
print(f'SE(b) : {lm_model.bse["BMI"]:.6f}')
print(f't値 : {lm_model.tvalues["BMI"]:.6f}')
print(f'p値 : {lm_model.pvalues["BMI"]:.6e}')
ci = lm_model.conf_int().loc['BMI']
print(f'95% CI : [{ci[0]:.6f}, {ci[1]:.6f}]')
--- 定義式による計算 ---
b : 10.233128
SE(b) : 0.673796
t値 : 15.187290
p値 : 0.000000e+00
95% CI : [8.908870, 11.557385]
--- statsmodels ---
b : 10.233128
SE(b) : 0.673796
t値 : 15.187290
p値 : 3.466006e-42
95% CI : [8.908870, 11.557385]
statsmodels を用いると、$x=20$ のときの信頼区間と予測区間も簡単に確認することができる。
# BMI 20 (=x_new)のとき
# 手計算による予測
prediction_manual = a+b*x_new
# statsmodelsによる予測
prediction_model = lm_model.predict(pd.DataFrame({'BMI': [x_new]}))[0]
print(f"手計算による予測値 : {prediction_manual:.2f}")
print(f"statsmodelsによる予測値: {prediction_model:.2f}")
手計算による予測値 : 86.89
statsmodelsによる予測値: 86.89
statsmodels の get_prediction() と summary_frame() を利用する。
# BMI 20 (=x_new)のときの信頼区間と予測区間
pred_interval = lm_model.get_prediction(pd.DataFrame({'BMI':[x_new]}))
pred_frame = pred_interval.summary_frame(alpha=0.05)
print(pred_frame)
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 86.889 5.225 76.621 97.158 -36.405 210.183
それぞれの意味は次の通りである。
-
mean：予測値 -
mean_se：予測値の標準誤差 -
mean_ci_lower/mean_ci_upper：平均値に対する95%信頼区間 -
obs_ci_lower/obs_ci_upper：個々の観測値に対する95%予測区間
今回の BMI = 20の場合、
- 予測値：86.889
- 95%信頼区間：約76.621 ～ 97.158
- 95%予測区間：約-36.405 ～ 210.183
次に、信頼区間と予測区間をグラフで確認する。まずは、seabornで信頼区間も表示したデータがプロットできる。
sns.lmplot(x='BMI', y='disease_progression', data=data, height=4, aspect=1.5)
statsmodels で得られるデータから信頼区間と予測区間のデータを作成する。
pred_all = lm_model.get_prediction() # 引数になにも指定しない場合は訓練データへの当てはめ値が出力
pred_frame_all = pred_all.summary_frame(alpha=0.05)
pred_graph = pd.concat([data.BMI, pred_frame_all], axis=1) # 説明変数を付ける
pred_graph = pred_graph.sort_values("BMI")
# 散布図
sns.scatterplot(x='BMI', y='disease_progression', data=data, color='steelblue', alpha=0.8)
# 回帰直線
sns.lineplot(x='BMI', y='mean', data=pred_graph, color='red', alpha=0.8)
# 信頼区間
sns.lineplot(x='BMI', y='mean_ci_lower', data=pred_graph, color='black', linestyle='dashed', alpha=0.6)
sns.lineplot(x='BMI', y='mean_ci_upper', data=pred_graph, color='black', linestyle='dashed', alpha=0.6)
# 予測区間
sns.lineplot(x='BMI', y='obs_ci_lower', data=pred_graph, color='black', linestyle='dotted', alpha=0.6)
sns.lineplot(x='BMI', y='obs_ci_upper', data=pred_graph, color='black', linestyle='dotted', alpha=0.6)
信頼区間は「BMIごとの平均的な疾患進行度」を推定する不確実性を表し、予測区間は「BMIがその値の個人の疾患進行度」を予測する不確実性を表すため、個々のデータのばらつきも含み、信頼区間より広くなる。
4. 各種計算式の導出
(編集予定)