📝

Lasso回帰で糖尿病データセットを解析

2022/03/17に公開

Python

機械学習

prml

tech

はじめに

PRML(Pattern Recognition and Machine Learning)を参考にして、実際のデータを使って学習します。主に3.1.4の内容です。
最尤推定による解の導出は、からっぽのしょこを参考にしました（この記事では導出には触れません）。
Lasso回帰の正則化項には主に過学習を防ぐ役割がありますが、その他にも正則化係数項の係数 $\lambda$ を大きくするほど、0になるパラメータが増加するという面白い性質があります( $\lambda$ をある程度大きくとればスパース推定となる）。この性質は、最尤推定解を見れば納得できると思います。

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$

二乗和の誤差関数 $E_D(\mathbf{w})$ を考えます。

\begin{aligned} E_D(\mathbf{w}) = \frac{1}{2}\sum^{N}_{n=1}\left\{t_n - \mathbf{w}^T\mathbf{\phi(x_n)}\right\} \end{aligned}

また、Lasso回帰では以下の正則化項 $E_W(\mathbf{w})$ を用います。正則化項は過学習を抑制するためにあります。

\begin{aligned} E_W(\mathbf{w}) = \frac{1}{2}\sum^{M-1}_{j=1}|w_j| \end{aligned}

次のような誤差関数 $E(\mathbf{w})$ を最小化することでLasso回帰の学習を行います。

\begin{aligned} E(\mathbf{w}) = E_D(\mathbf{w}) + \lambda E_W(\mathbf{w}) \end{aligned}

ここで、 $\lambda$ は二乗和の誤差関数 $E_D(\mathbf{w})$ と正則化項 $E_W(\mathbf{w})$ の比重を調整するハイパーパラメータです。

トイデータで学習

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個の基底関数を作ります。ガウス基底関数やシグモイド基底関数も使っています。今回は $s=0.5$ としています。

$\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=8$ の時だということが分かります。

結果と考察

様々な $\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

ラッソ回帰では一般に $\lambda$ が増加するにつれて、0でないパラメータの数が減少していきますが、この学習においてもそうなっていることが分かります。

パラメータの値

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

上記は、 $\lambda=8$ の時のパラメータの中で絶対値が $0.1$ より大きいものです。実際は0でないパラメータは33個ありますが、絶対値が0.1より大きいものだけを表示しました。この他にも2つの説明変数の積の関数、二次関数、三次関数のパラメータには0でないものがありましたが、シグモイド基底関数はこのデータに合ってないのか、全てのパラメータが0になっていました。

また、上で表示したパラメータに対応する基底関数は以下の通りです。善玉コレステロールの重みパラメータのみが負となっていて、そのほかの（バイアスを除いた）重みパラメータは正となっています。

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回帰をすることで、そのほかの重要な特徴が見つかるかもしれません。