Lasso回帰で糖尿病データセットを解析
はじめに
PRML(Pattern Recognition and Machine Learning)を参考にして、実際のデータを使って学習します。主に3.1.4の内容です。
最尤推定による解の導出は、からっぽのしょこを参考にしました(この記事では導出には触れません)。
Lasso回帰の正則化項には主に過学習を防ぐ役割がありますが、その他にも正則化係数項の係数
Lasso回帰とは
- 説明変数
\left\{\mathbf{x_1}, \mathbf{x_2}, \cdots \mathbf{x_N}\right\}, \mathbf{x_n} = (x_{n, 1}, x_{n, 2}, \cdots, x_{n_D})^T - 目的変数
\left\{t_1, t_2, \cdots t_N\right\} - 基底関数
\mathbf{\phi(x_n)} = (\phi_0(\mathbf{x_n}), \phi_1(\mathbf{x_n}), \cdots \phi_{M-1}(\mathbf{x_n}))^T - 重み
\mathbf{w} = (w_0, w_1, \cdots, w_{M-1})^T
二乗和の誤差関数
また、Lasso回帰では以下の正則化項
次のような誤差関数
ここで、
トイデータで学習
import文
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
np.random.seed(28)
今回使用するデータセット
scikit-learnの「糖尿病データセット」というトイデータセットを使います。このデータセットでは、年齢、性別、BMI、平均血圧、6つの血清検査測定値という10の説明変数から、1年後の糖尿病の進行状況を予測します。データ数は442です。このうちの8割に当たる353個を訓練用データとし、残りをテスト用データとします。また、データは標準化しておきます。
input_variables = ['age', 'sex', 'bmi', 'map', 'tc', 'ldl', 'hdl', 'tch', 'ltg', 'glu']
data_url = 'https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt'
raw_df = pd.read_csv(data_url, sep="\s+")
diabetes_df = raw_df.set_axis(input_variables + ['y'], axis=1)
y_origin_df = diabetes_df.loc[:, 'y']
y_df = (y_origin_df - y_origin_df.mean()) / y_origin_df.std()
x_origin_df = diabetes_df.loc[:, input_variables]
x_df = (x_origin_df - x_origin_df.mean()) / x_origin_df.std()
train_x_df, test_x_df, train_y_df, test_y_df = train_test_split(x_df, y_df, test_size=0.2)
train_x = train_x_df.to_numpy(copy=True)
train_y = train_y_df.to_numpy(copy=True)
test_x = test_x_df.to_numpy(copy=True)
test_y = test_y_df.to_numpy(copy=True)
diabetes_df.head()
index | age | sex | bmi | map | tc | ldl | hdl | tch | ltg | glu | y |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 59 | 2 | 32.1 | 101.0 | 157 | 93.2 | 38.0 | 4.0 | 4.8598 | 87 | 151 |
1 | 48 | 1 | 21.6 | 87.0 | 183 | 103.2 | 70.0 | 3.0 | 3.8918 | 69 | 75 |
2 | 72 | 2 | 30.5 | 93.0 | 156 | 93.6 | 41.0 | 4.0 | 4.6728 | 85 | 141 |
3 | 24 | 1 | 25.3 | 84.0 | 198 | 131.4 | 40.0 | 5.0 | 4.8903 | 89 | 206 |
4 | 50 | 1 | 23.0 | 101.0 | 192 | 125.4 | 52.0 | 4.0 | 4.2905 | 80 | 135 |
基底関数
296個の基底関数を作ります。ガウス基底関数やシグモイド基底関数も使っています。今回は
\phi_0(\mathbf{x_n}) = 1 \phi_j(\mathbf{x_n}) = x_{n, k}x_{n, l} \quad(1 \leq j \leq 45) \phi_j(\mathbf{x_n}) = x_{n, k} \quad(46 \leq j \leq 55) \phi_j(\mathbf{x_n}) = x_{n, k}^2 \quad(56 \leq j \leq 65) \phi_j(\mathbf{x_n}) = x_{n, k}^3 \quad(66 \leq j \leq 75) \phi_j(\mathbf{x_n}) = exp(-\frac{(x_{n, k} - \mu_{l})^2}{2s^2}) \quad(76 \leq j \leq 185) \phi_j(\mathbf{x_n}) = \frac{1}{1 + exp(-\frac{x_{n, k} - \mu_{l}}{s})} \quad(186 \leq j \leq 295)
def init_base_func_list():
base_func_list = [0 for i in range(45)]
ind = 0
for i in range(10):
for j in range(i+1, 10):
base_func_list[ind] = (i, j)
ind += 1
return base_func_list
M = 296
base_func_list = init_base_func_list()
s = 0.5
def base_func(x_n, m):
if m == 0:
return np.ones(len(x_n))
elif m < 46:
return x_n[:, base_func_list[m-1][0]] * x_n[:, base_func_list[m-1][1]]
elif m < 56:
return x_n[:, m-46]
elif m < 66:
return x_n[:, m-56]**2
elif m < 76:
return x_n[:, m-66]**3
elif m < 186:
m -= 76
x_ind = m // 11
gauss_mean = -1.5 + 0.3 * (m % 11)
return np.exp(-(x_n[:, x_ind] - gauss_mean)**2 / (2 * s**2))
elif m < 296:
m -= 186
x_ind = m // 11
sigmoid_mean = -1.5 + 0.3 * (m % 11)
return 1 / (1 + np.exp((sigmoid_mean - x_n[:, x_ind]) / s))
raise ValueError()
その他必要な関数
まずは、基底関数に説明変数を入力したものです。
def Phi(x, M):
x_nm = np.zeros((len(x), M))
for m in range(M):
x_nm[:, m] = base_func(x, m)
return x_nm
phi_x_nm = Phi(train_x, M)
phi_x_test = Phi(test_x, M)
パラメータの最尤推定の際に必要となる場合分けをする関数です。
def ML_thresholding(S, lam, phi_x_n):
if S > lam:
return (S - lam) / np.sum(phi_x_n**2)
elif S < -lam:
return (S + lam) / np.sum(phi_x_n**2)
else:
return 0
二乗平均平方根誤差(RMSE)を求める関数です。目的変数の分散をかけることで、標準化する前のRMSEを求めています。
def RMSE(phi, y, w):
error_sum = sum((y - np.dot(phi, w)) ** 2) * y_origin_df.std() ** 2
return np.sqrt(error_sum / len(y))
学習のコード
Lasso回帰の最尤推定による解は直接求めることができないので、座標降下法によって重りパラメータを繰り返し更新することで求めます。
MAX_ITER = 200
LAMBDAS = [2**i for i in range(10)]
w_ms = np.empty((len(LAMBDAS), M))
w_m = np.random.normal(0.0, 0.05, M)
for lam in range(len(LAMBDAS)):
for i in range(MAX_ITER):
w_m[0] = 0
w_m[0] = sum(train_y - np.dot(phi_x_nm, w_m)) / len(train_y)
for m in range(1, M):
w_m[m] = 0
S = sum((train_y - np.dot(phi_x_nm, w_m)) * phi_x_nm[:, m])
w_m[m] = ML_thresholding(S, LAMBDAS[lam], phi_x_nm[:, m])
w_ms[lam] = w_m
train_RMSE = RMSE(phi_x_nm, train_y, w_m)
test_RMSE = RMSE(phi_x_test, test_y, w_m)
print("LAMBDA, train_RMSE, test_RMSE | {0:3d}, {1:.3f}, {2:.3f}".format(LAMBDAS[lam], train_RMSE, test_RMSE))
LAMBDA, train_RMSE, test_RMSE | 1, 45.287, 60.693
LAMBDA, train_RMSE, test_RMSE | 2, 46.375, 58.352
LAMBDA, train_RMSE, test_RMSE | 4, 47.683, 55.806
LAMBDA, train_RMSE, test_RMSE | 8, 49.728, 55.149
LAMBDA, train_RMSE, test_RMSE | 16, 51.904, 57.073
LAMBDA, train_RMSE, test_RMSE | 32, 54.529, 60.269
LAMBDA, train_RMSE, test_RMSE | 64, 58.275, 64.897
LAMBDA, train_RMSE, test_RMSE | 128, 61.292, 67.756
LAMBDA, train_RMSE, test_RMSE | 256, 65.199, 69.819
LAMBDA, train_RMSE, test_RMSE | 512, 74.262, 75.675
テストデータのRMSEが一番低いのは
結果と考察
\lambda に対する0でないパラメータの数
様々なfor i in range(10):
cnt = 0
for m in range(M):
if w_ms[i][m] != 0:
cnt += 1
print("LAMBDA, number of non-zero parameters | {0:3d}, {1:2d}".format(LAMBDAS[i], cnt))
LAMBDA, number of non-zero parameters | 1, 83
LAMBDA, number of non-zero parameters | 2, 65
LAMBDA, number of non-zero parameters | 4, 48
LAMBDA, number of non-zero parameters | 8, 33
LAMBDA, number of non-zero parameters | 16, 20
LAMBDA, number of non-zero parameters | 32, 12
LAMBDA, number of non-zero parameters | 64, 9
LAMBDA, number of non-zero parameters | 128, 7
LAMBDA, number of non-zero parameters | 256, 6
LAMBDA, number of non-zero parameters | 512, 2
ラッソ回帰では一般に
パラメータの値
for m in range(M):
if abs(w_ms[3][m]) > 0.1:
print("m, w_ms[3][m] | {0:3d}, {1:.2f}".format(m, w_ms[3][m], m))
m, w_ms[3][m] | 0, -0.16
m, w_ms[3][m] | 48, 0.20
m, w_ms[3][m] | 49, 0.11
m, w_ms[3][m] | 52, -0.15
m, w_ms[3][m] | 54, 0.35
m, w_ms[3][m] | 105, 0.13
上記は、
また、上で表示したパラメータに対応する基底関数は以下の通りです。善玉コレステロールの重みパラメータのみが負となっていて、そのほかの(バイアスを除いた)重みパラメータは正となっています。
- 0: 1(バイアス)
- 48: bmi(Body Mass Index)の一次関数
- 49: map(平均血圧)の一次関数
- 52: hdl(善玉コレステロール)の一次関数
- 54: ltg(「possibly log of serum triglycerides level」??)の一次関数
- 105: bmiの平均29.03のガウス基底関数
考察
Lasso回帰を行うと正則化項が大きいほど多くのパラメータが0となるため、目的変数の予測ができるだけでなく、0にならずに残ったパラメータからデータの大まかな特徴をとらえることもできます。例えば相関係数を見ることでもデータの特徴をつかめると思いますが、Lasso回帰をすることで、そのほかの重要な特徴が見つかるかもしれません。
Discussion