はじめに
本稿では, 中古マンションの売り出し価格に寄与する変数を, 線形回帰モデルを用いて確認します.
断言調になっている箇所のほとんどは, 調べたことをまとめているに過ぎないので予めご了承ください.
使用したデータ
使用したデータは非公開で, 内容は以下のとおりです.
- データ数: 20件
- 関心の駅名で検索をかけるとヒットする中古マンションの1室が1つのレコードに対応
- 目的変数は売り出し中の販売価格
- 不動産情報サイトを参考に筆者が手作業で作成
環境
- Google Colaboratory
関連ライブラリのインポート
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.preprocessing import (
StandardScaler
)
from sklearn.linear_model import (
LassoLarsIC,
lasso_path,
LassoCV,
Lasso,
LassoLarsCV,
LassoLarsIC
)
import itertools
from itertools import cycle
データの読み込み
fn = 'data.csv'
raw_df = pd.read_csv(fn)
raw_df.shape, raw_df.dtypes
"""
((20, 13),
id int64
name object
renovated int64
area_square_meter float64
walking_minute int64
floor int64
has_service_room int64
x_ldk int64
service_charge_yen int64
selling_price_yen int64
has_elevator int64
year_construction int64
repair_reserve_fund int64
dtype: object)
"""
データの加工
そのままでは使いづらいかもしれないので, 以下のコードで簡単なクリーニングを行いました.
# make a copy of "raw_df" as "df":
df = raw_df.copy()
# define a function for data cleaning:
def modify_dtypes_to_category(df, cols):
# change dtypes of "cols" to "category":
df[cols] = df[cols].astype("category")
return df
# define another function for data cleaning:
def drop_cols(df, cols):
# remove the variables:
df = df.drop(cols, axis=1)
return df
# implement the functions defined above:
## prepare a list of column names for modifying dtypes
cols_to_category = [
'renovated',
'has_service_room',
'has_elevator'
]
## change dtypes of intended columns to "category"
df = modify_dtypes_to_category(df=df, cols=cols_to_category)
## prepare a list of column names to be dropped
cols_to_drop = [
'name',
'id'
]
## drop the columns
df = drop_cols(df=df, cols=cols_to_drop)
df.shape, df.dtypes
"""
((20, 11),
renovated category
area_square_meter float64
walking_minute int64
floor int64
has_service_room category
x_ldk int64
service_charge_yen int64
selling_price_yen int64
has_elevator category
year_construction int64
repair_reserve_fund int64
dtype: object)
"""
更に, データを変数Xとyにそれぞれ格納します.
# split the data to X and y:
def split_data(df, target):
X = df.drop(labels=[target], axis=1)
y = df[target]
return X, y
target = 'selling_price_yen'
X, y = split_data(df=df, target=target)
X.shape, y.shape
"""
((20, 10), (20,))
"""
本編
以下の線形モデルを仮定します.
\boldsymbol{y} = X\boldsymbol{\beta} + \epsilon
\begin{align}
\boldsymbol{y}&: 目的変数に関するn個のデータからなるn次元観測値ベクトル \\
X &: n \times (p+1)次元計画行列 \\
\boldsymbol{\beta}&: (p+1)次元回帰係数ベクトル \\
\boldsymbol{\epsilon}&: n次元誤差ベクトル \\
\end{align}
今回は$n, p = 20, 10$.
モデルの推定
(a)最小二乗法
最小二乗法は誤差の二乗和を最小にするような回帰式を求める方法.
誤差について以下の仮定を置くことにより, 良い推定量であることが保証されます.
\begin{align}
E[\boldsymbol{\epsilon}]
&= \boldsymbol{0}^T \\
E[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T]
&= \sigma^2 I_n \\
\end{align}
\begin{align}
\boldsymbol{0}
&: 要素がすべて0のn次元ベクトル \\
I_n
&: n \times n 単位行列\\
\sigma^2
&: 誤差分散 \\
\end{align}
それぞれの観測点における誤差の期待値及び分散が全て0, $\sigma^2$で, 且つ互いに無相関ということになります.
(今回の目的変数は中古マンションの売り出し価格で, その数値の増加に伴いばらつきも増加しそうですので, 直感として上記の仮定は誤りなるものと理解)
また, 誤差の二乗和は,
S(\boldsymbol{\beta})
= \sum^n_{i=1} \epsilon_i^2
= \boldsymbol{\epsilon}^T \boldsymbol{\epsilon}
= (\boldsymbol{y} - X \boldsymbol{\beta})^T
(\boldsymbol{y} - X \boldsymbol{\beta})
のように表されます.
上式を最小化するための$\boldsymbol{\beta}$は,
\frac{\partial S(\boldsymbol{\beta})}{\partial \beta}
= -2 X^T \boldsymbol{y} + 2 X^T X \boldsymbol{\beta}
= \boldsymbol{0}
で与えられる正規方程式の解$\boldsymbol{\hat{\beta}}$となります.
もしも$X^{T} X$の逆行列が存在すれば, $\boldsymbol{\hat{\beta}}$は,
\hat{\boldsymbol{\beta}}
= (X^{T} X)^{-1} X^{T} \boldsymbol{y}
となります.
予測値ベクトルy$, $残差ベクトルeをそれぞれ,
\begin{align}
\hat{\boldsymbol{y}}
&= X \boldsymbol{\hat{\beta}}
= X (X^{T} X)^{-1} X^{T} \boldsymbol{y} \\
\boldsymbol{e}
&= \boldsymbol{y} - \hat{\boldsymbol{y}}
= [I_{n} - X (X^{T} X)^{-1} X^{T}] \boldsymbol{y}
\end{align}
のようにし,
$$
P = I_{n} - X (X^{T} X)^{-1} X^{T}
$$
が, $X$の直交補空間への射影行列で、冪等行列($P = P^{2}$)であることに注意すると, 残差平方和は,
\boldsymbol{e}^{T} \boldsymbol{e}
= \boldsymbol{y}^{T} [I_{n} - X (X^{T} X)^{-1} X^{T}] \boldsymbol{y}
のように表されます.
最小二乗推定量$\boldsymbol{\hat{\beta}}$の期待値, 分散共分散行列は,
\begin{align}
&E[
\boldsymbol{\hat{\beta}}
]
= \boldsymbol{\beta} \\
&E[
(\boldsymbol{\hat{\beta}} - \boldsymbol{\beta})
(\boldsymbol{\hat{\beta}} - \boldsymbol{\beta})^T
]
= \sigma^2(X^{T} X)^{-1}
\end{align}
となります.
Pythonの場合, statsmodel.formula.api.ols()を使用して, モデルをフィットし, サマリーを表示することが可能.
# Construct the formula for multiple linear regression:
# Use the variable names to build the formula.
dependent_variable = target
independent_variables = " + ".join(X.columns.values)
formula = f"{dependent_variable} ~ {independent_variables}"
# Create and fit the multiple linear regression model:
# Use the constructed formula and the dataset to create the model.
model_multiple_lr = smf.ols(formula=formula, data=df)
results = model_multiple_lr.fit()
# Display the summary of the regression results:
# Print the summary of the fitted model.
print(results.summary())
コードを実行すると, テキストベースのサマリーが出力されます.
縦方向に長いので, 3つに分けて確認したいと思います.
下図は1つ目のセクションです.
OLS Regression Results
==============================================================================
Dep. Variable: selling_price_yen R-squared: 0.940
Model: OLS Adj. R-squared: 0.874
Method: Least Squares F-statistic: 14.21
Date: Thu, 23 Nov 2023 Prob (F-statistic): 0.000239
Time: 15:57:54 Log-Likelihood: -322.58
No. Observations: 20 AIC: 667.2
Df Residuals: 9 BIC: 678.1
Df Model: 10
Covariance Type: nonrobust
上図を見ると, 自由度調整後の決定係数が0.874と, 一見すると悪くない数値に見えます.
続いて, 2つ目のセクションです.
=========================================================================================
coef std err t P>|t| [0.025 0.975]
-----------------------------------------------------------------------------------------
Intercept -6.977e+08 2.8e+08 -2.492 0.034 -1.33e+09 -6.43e+07
renovated[T.1] 8.841e+06 2.16e+06 4.101 0.003 3.96e+06 1.37e+07
has_service_room[T.1] -1.678e+06 3.45e+06 -0.486 0.639 -9.49e+06 6.14e+06
has_elevator[T.1] 5.738e+06 3.78e+06 1.516 0.164 -2.82e+06 1.43e+07
area_square_meter 2.931e+05 1.53e+05 1.910 0.088 -5.41e+04 6.4e+05
walking_minute -1.087e+06 2.99e+05 -3.628 0.005 -1.76e+06 -4.09e+05
floor 2.833e+05 2.64e+05 1.072 0.312 -3.15e+05 8.81e+05
x_ldk 3.377e+06 2.93e+06 1.151 0.279 -3.26e+06 1e+07
service_charge_yen 449.0882 348.329 1.289 0.229 -338.887 1237.063
year_construction 3.496e+05 1.41e+05 2.473 0.035 2.98e+04 6.69e+05
repair_reserve_fund 2.0551 62.176 0.033 0.974 -138.596 142.706
==============================================================================
Omnibus: 0.182 Durbin-Watson: 1.936
Prob(Omnibus): 0.913 Jarque-Bera (JB): 0.280
Skew: -0.190 Prob(JB): 0.869
Kurtosis: 2.563 Cond. No. 9.83e+06
==============================================================================
信頼区間に0を含まない変数は,
- Intercept
- renovated
- walking_minute
- year_construction
となりました.
一方で,
- area_square_meter
- has_elevator
など, 売り出し中の販売価格に影響を与えそうな変数のうちいくつかは, 偏回帰係数が0と有意に異なると計算されなかったことが分かります.
今回は, 1つの駅近傍の売り出し中の中古マンションのデータを使用しているため、データ数が20と少ないです.
そのため, 回帰係数の確認はこのあたりにしておきたいと思います.
次に3つ目のセクションです.
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.83e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
condition numberなる数値が高く, 多重共線性やその他の数値的な問題があるようです.
以下は公式ドキュメントの説明です.
Calculated as ratio of largest to smallest singular value of the exogenous variables. This value is the same as the square root of the ratio of the largest to smallest eigenvalue of the inner-product of the exogenous variables.
inner-product of the exogenous variablesとは, $X^T X$のことですね.
ググってみると, condition numberが大きいと, $X^T X$の逆行列の数値が, $X^T X$の変化に対して不安定になるという記載がありますが, そのことが線形回帰モデルの評価とどのように関係しているのかは理解が及びませんでした.
ひとまずは, 大きなcondition numberは多重共線性が発生している可能性を示唆している, とだけ理解しておきます.
実際, 例えばarea_square_meterとrepair_reserve_fundには強い相関があるなど, 多重共線性が発生しているのだと思います.
(b)最尤法
最尤法は今回(ガウス型線形回帰モデル)の場合, 目的変数が正規分布に従って観測されたと仮定し, 尤度関数を最大化するようなモデルを学習します.
\begin{align}
\boldsymbol{y} = X \boldsymbol{\beta} + \boldsymbol{\epsilon}
\qquad
\boldsymbol{\epsilon} \sim N(0, \sigma^2)
\end{align}
上式のように, 誤差項の従う分布を明示します.
所定の対数尤度関数を最大化するための尤度方程式を解くことによって, $\beta$と$\sigma$の最尤推定値はそれぞれ,
\begin{align}
\boldsymbol{\hat{\beta}}
&= X (X^{T} X)^{-1} X^{T} \boldsymbol{y} \\
\hat{\sigma}^2
&= \frac{1}{n} (\boldsymbol{y} - \hat{\boldsymbol{y}})^{T} (\boldsymbol{y} - \hat{\boldsymbol{y}})
\end{align}
と表されます.
これは, 理論上最小二乗法と同じ回帰係数を得ることになるので, 確認してみたいと思います.
Pythonの場合, statsmodel.formula.api.glm()を使用して, モデルをフィットし, サマリーを表示することが可能.
# Construct the formula for multiple linear regression with generalized linear models:
# Use the variable names to build the formula.
dependent_variable = target
independent_variables = " + ".join(X.columns.values)
formula = f"{dependent_variable} ~ {independent_variables}"
# Create and fit the multiple linear regression model using generalized linear models:
# Use the constructed formula, dataset, and Gaussian family with Identity link.
model_multiple_lr = smf.glm(
formula=formula, data=df,
family=sm.families.Gaussian(link=sm.genmod.families.links.Identity())
)
results = model_multiple_lr.fit()
# Display the summary of the regression results:
# Print the summary of the fitted model.
print(results.summary())
最尤法による回帰モデルは一般化線形モデル(generalized linear models)というフレームワークに整理されており, 実際statsmodels.apiでも, 最尤法はglm()という一般化線形モデルに関するクラスから実施します.
一般化線形モデルでは誤差が従う分布について, 指数分布族と呼ばれる確率分布のファミリーの中から指定することが可能で, 今回はガウス型線形回帰モデルを想定しているため, 引数familyにGaussian()を指定しています.
同様に, 一般化線形モデルではリンク関数という概念を用いて, 特徴量の線形結合と, 目的変数(y: endogenous variable)の橋渡しを行います.
今回のガウス型線形回帰モデルの場合は, 特徴量の線形結合をそのまま目的変数とするため, sm.families.Gaussian()の引数linkにsm.genmod.families.links.Identity()を指定しました.
以下は上記コードの実行結果です.
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: selling_price_yen No. Observations: 20
Model: GLM Df Residuals: 9
Model Family: Gaussian Df Model: 10
Link Function: Identity Scale: 1.3297e+13
Method: IRLS Log-Likelihood: -322.58
Date: Thu, 23 Nov 2023 Deviance: 1.1967e+14
Time: 21:56:54 Pearson chi2: 1.20e+14
No. Iterations: 3 Pseudo R-squ. (CS): 0.9994
Covariance Type: nonrobust
=========================================================================================
coef std err z P>|z| [0.025 0.975]
-----------------------------------------------------------------------------------------
Intercept -6.977e+08 2.8e+08 -2.492 0.013 -1.25e+09 -1.49e+08
renovated[T.1] 8.841e+06 2.16e+06 4.101 0.000 4.62e+06 1.31e+07
has_service_room[T.1] -1.678e+06 3.45e+06 -0.486 0.627 -8.45e+06 5.09e+06
has_elevator[T.1] 5.738e+06 3.78e+06 1.516 0.130 -1.68e+06 1.32e+07
area_square_meter 2.931e+05 1.53e+05 1.910 0.056 -7679.393 5.94e+05
walking_minute -1.087e+06 2.99e+05 -3.628 0.000 -1.67e+06 -5e+05
floor 2.833e+05 2.64e+05 1.072 0.284 -2.35e+05 8.01e+05
x_ldk 3.377e+06 2.93e+06 1.151 0.250 -2.37e+06 9.13e+06
service_charge_yen 449.0882 348.329 1.289 0.197 -233.624 1131.801
year_construction 3.496e+05 1.41e+05 2.473 0.013 7.25e+04 6.27e+05
repair_reserve_fund 2.0551 62.176 0.033 0.974 -119.807 123.917
=========================================================================================
確かにOLS(ordinary least squares: 最小二乗法)の結果と一致しました.
変数選択
情報量規準AIC(Akaike's information criterion)による変数選択を試みたいと思います.
AICは一般的に,
$$
AIC = -2(モデルの最大対数尤度) + 2(モデルの自由パラメータ数)
$$
の式で与えられます.
10個の説明変数すべての組み合わせに対して, ガウス型線形回帰モデルを計算し, 評価します.
以下がPythonのコードです.
# List of descriptors (independent variables)
desc_list = X.columns
def generate_formulas(descriptors, num):
"""
Generate combinations of formulas with a given number of descriptors.
Args:
- descriptors: List of descriptors (independent variables).
- num: Number of descriptors in each formula.
Returns:
- A generator yielding formulas.
"""
if num > len(descriptors):
return None
else:
for f in itertools.combinations(descriptors, num):
formula = f'{target} ~ ' + '+'.join(f)
yield formula
def modeling_linear(formula, data):
"""
Fit a generalized linear model using the given formula and data.
Args:
- formula: Formula for the model.
- data: DataFrame containing the data.
Returns:
- Tuple containing formula and fitted model.
"""
model = smf.glm(
formula=formula,
data=data,
family=sm.families.Gaussian(link=sm.genmod.families.links.Identity())
).fit()
return formula, model
# Initialize lists to store results
n_var = []
formula = []
model = []
AIC = []
pseudo_r2 = []
sigma_hat_squared = []
llf = []
# Null model
null_model = smf.glm(
formula=f'{target} ~ 1',
data=data,
family=sm.families.Gaussian(link=sm.genmod.families.links.Identity())
).fit()
n_var.append(0)
formula.append(f'{target} ~ 1')
model.append(null_model)
AIC.append(null_model.aic)
pseudo_r2.append(null_model.pseudo_rsquared())
sigma_hat_squared.append(((null_model.predict(df) - y)**2).sum() / len(df))
llf.append(null_model.llf)
# Model building loop
for i in range(1, len(desc_list)):
for f in generate_formulas(desc_list, i):
eq, m = modeling_linear(formula=f, data=data)
n_var.append(i)
formula.append(f)
model.append(m)
AIC.append(m.aic)
pseudo_r2.append(m.pseudo_rsquared())
sigma_hat_squared.append(((m.predict(df) - y)**2).sum() / len(df))
llf.append(m.llf)
# Create a DataFrame to store the results and sort by AIC
tmp_df = pd.DataFrame({
"n_var": n_var,
"sigma_hat_squared": sigma_hat_squared,
"log_likelihood": llf,
'AIC': AIC,
"formula": formula,
"model": model
}).sort_values(by="AIC", ascending=True)
# Display the top 10 models based on AIC
top_10_models = tmp_df[["n_var", "sigma_hat_squared", "log_likelihood", "AIC"]].head(10)
print(top_10_models)
説明変数の数が10だけでしたので, 全組み合わせを総当たりしました.
以下が実行結果で, AICを昇順(良い評価の順)に10個分表示させたものです.
n_var sigma_hat_squared log_likelihood AIC
404 5 7.656279e+12 -325.044243 662.088486
669 6 7.421884e+12 -324.733311 663.466623
666 6 7.428721e+12 -324.742519 663.485039
660 6 7.452729e+12 -324.774785 663.549569
878 7 6.775834e+12 -323.822606 663.645213
978 8 6.199225e+12 -322.933225 663.866450
672 6 7.601233e+12 -324.972087 663.944174
650 6 7.634222e+12 -325.015393 664.030786
181 4 9.394501e+12 -327.090227 664.180454
402 5 8.643370e+12 -326.256908 664.513815
インデックス404のモデルが最小のAICを与えるモデルで, 使用されている説明変数は5個でした.
summary()で中身を確認します.
以下はPythonコードと実行結果です.
tmp_df.loc[404]["model"].summary()
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: selling_price_yen No. Observations: 20
Model: GLM Df Residuals: 14
Model Family: Gaussian Df Model: 5
Link Function: Identity Scale: 1.0938e+13
Method: IRLS Log-Likelihood: -325.04
Date: Thu, 23 Nov 2023 Deviance: 1.5313e+14
Time: 23:50:23 Pearson chi2: 1.53e+14
No. Iterations: 3 Pseudo R-squ. (CS): 0.9998
Covariance Type: nonrobust
=====================================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------------
Intercept -8.485e+08 2.33e+08 -3.648 0.000 -1.3e+09 -3.93e+08
renovated[T.1] 7.706e+06 1.67e+06 4.618 0.000 4.44e+06 1.1e+07
has_elevator[T.1] 5.559e+06 3.12e+06 1.783 0.075 -5.52e+05 1.17e+07
area_square_meter 4.45e+05 8.32e+04 5.350 0.000 2.82e+05 6.08e+05
walking_minute -9.893e+05 1.88e+05 -5.260 0.000 -1.36e+06 -6.21e+05
year_construction 4.274e+05 1.17e+05 3.666 0.000 1.99e+05 6.56e+05
=====================================================================================
P-値が有意となる変数が多いですね.
回帰係の解釈としては,
- リノベーションの有の場合, +770万円
- エレベータ有の場合, +560万円
- 1平米あたり+45万円
- 駅からの徒歩1分あたり-100万円
- 築造年(西暦ベース)1年あたり43万円
- 切片-8億5千万円
となり, 切片を除きそれなりに常識的な係数が推定されていることが分かります.
切片については, 全ての特徴量を使用したモデルと同様, 築造年(西暦ベース)で露骨な調整が起こっているようです.
築造年が有益な特徴量なので切片がこのようになったと思われますが, 例えば2,000など, きりの良い年で中心化した方が, 直感的な理解が得やすそうです.
モデルの評価
先ほどは変数選択の基準としてAICを利用し, 所与の特徴量の中から, 最小のAICを持つ特徴量の組み合わせ及び対応するモデルを特定しました.
ここでは, クロスバリデーションを用いて, 特定されたモデルの残差を評価します.
今回は$n$個の観測データの中か$i$番目のデータを除く残り$(n-1)$個のデータに基づき重みベクトル$\boldsymbol{\beta}$を推定し, これを$\hat{\boldsymbol{\beta}}^{(-i)}$とします.
次に, 推定した回帰モデル$y = u(\boldsymbol{x}_i ; \hat{\boldsymbol{\beta}}^{(-i)})$に対して残差を計算し, クロスバリデーション
$$
CV
= \frac{1}{n} \sum_{i=1}^{n}
{y_i - u(\boldsymbol{x}_i ; \hat{\boldsymbol{\beta}}^{(-i)}) }^2
$$
を得ます.
Pythonでは, sklearnにLeaveOneOut()という関数があるようですが, シンプルに以下のコードで計算しました.
# Initialize a list to store residuals for Leave-One-Out Cross Validation
loocv_residuals = []
# Perform Leave-One-Out Cross Validation
for i in range(len(df)):
# Split the dataset into test and training data
test_data = df.iloc[[i]]
train_data = df.drop(index=i)
# Train the model on the training data
model = smf.ols(
formula='selling_price_yen ~ renovated + area_square_meter + walking_minute + has_elevator + year_construction',
data=train_data
).fit()
# Predict using the test data
y_pred = model.predict(test_data)
# Calculate residuals
residual = y_pred - test_data[target]
loocv_residuals.append(residual)
# Calculate mean squared residuals for Leave-One-Out Cross Validation
msr = np.sum(np.square(loocv_residuals)) / len(df)
sqrt_msr = np.sqrt(msr)
print(f"Square root of mean squared residuals: {sqrt_msr:,.0f}")
# Calculate mean absolute residuals for Leave-One-Out Cross Validation
mar = np.sum(np.abs(loocv_residuals)) / len(df)
print(f"Mean absolute residuals: {mar:,.0f}")
"""
Square root of mean squared residuals: 3,975,210
Mean absolute residuals: 3,478,479
"""
実行結果に表示されている3,975,210は, 上述のクロスバリデーションによる残差二乗和の平均の正の平方根です.
同様の手順で, 全ての特徴量を使用したモデルについても計算したところ,
Square root of mean squared residuals: 17,433,699
Mean absolute residuals: 8,441,324
を得ました.
確かに, AICを最小にするモデルは, AICだけでなくクロスバリデーションによるMSRも小さくなるようです.
なお, AICを最小にするモデルが, 上述の（Leave-One-Out法の）クロスバリデーションによるMSRを最小にするモデルなのかどうかについては更なる確認を要します.
正則化
(a) 正則化最小二乗法
先ほどは, 10種類の特徴量をモデルの説明変数として使用するのかどうかについて, AICを評価指標とした上で, 全ての組み合わせで総当たり計算し
, 最小のAICを与える特徴量の組み合わせを得ました.
一般的にこのような方法は計算資源の問題に直面することが多いため, 今回は正則化を利用した変数選択の結果も試みます.
正則化(或いは罰則付き)最小二乗法(regularized(penalized) least squares method)では, 使用する特徴量の増加などといったモデルの複雑さが増すにつれて減少する誤差の二乗和の項に, 逆にモデルの複雑さに伴って大きくなる項$R(\boldsymbol{\beta})$を加えた関数
\begin{align}
S_{\gamma}(\boldsymbol{\beta})
&= (\boldsymbol{y} - X \boldsymbol{\beta})^{T}
(\boldsymbol{y} - X \boldsymbol{\beta})
+ \gamma R(\boldsymbol{\beta}) \\
\end{align}
\begin{align}
\gamma(>0)
&: 正則化パラメータ(或いは平滑化パラメータ)
\end{align}
の最小化を行います.
ペナルティ項$R(\boldsymbol{\beta})$は2次形式
$$
\boldsymbol{\beta}^{T} K \boldsymbol{\beta}
$$
で表現することができる場合が多く,
\begin{align}
S_{\gamma}(\boldsymbol{\beta})
&= (\boldsymbol{y} - X \boldsymbol{\beta})^{T}
(\boldsymbol{y} - X \boldsymbol{\beta})
+ \gamma \boldsymbol{\beta}^{T} K \boldsymbol{\beta} \\
\end{align}
\begin{align}
K
&: 非負値定符号行列
\end{align}
とすると,
\hat{\boldsymbol{\beta}}
= (X^T X + \gamma K)^{-1}
X^T \boldsymbol{y}
と, 解析的に陽に求めることが可能.
(b) 正則化最尤法
ガウス型線形回帰モデルの対数尤度関数は,
\begin{align}
l(\boldsymbol{\theta})
&= \sum_{i=1}^{n} log f(y_i | \boldsymbol{x}_i ; \boldsymbol{\beta}, \sigma^2) \\
&= - \frac{n}{2}log(2\pi\sigma^2)
- \frac{1}{2\sigma^2}(\boldsymbol{y} - X \boldsymbol{\beta})^T(\boldsymbol{y} - X \boldsymbol{\beta}) \\
\boldsymbol{\theta} &= (\boldsymbol{\beta}^T, \sigma^2)^T
\end{align}
上式に最尤推定量を代入すると,
l(\boldsymbol{\hat{\theta}})
= - \frac{n}{2}log(2\pi)
- \frac{n}{2}log\hat{\sigma}^2
- \frac{n}{2}
分散の最尤推定値は残差平方和から計算されるため, 対数尤度の値はモデルの複雑さに伴い増加します.
そこで, 正則化最小二乗法と同様に, ペナルティ項を導入し,
\begin{align}
l_{\lambda}(\boldsymbol{\theta})
&= \sum_{i=1}^{n} \log f(y_i | \boldsymbol{x}_i ; \boldsymbol{\beta}, \sigma^2)
- \frac{\lambda}{2}R(\boldsymbol{\beta}) \\
\\
\lambda(>0)
&: 正則化パラメータ(或いは平滑化パラメータ)
\end{align}
としたものは正則化対数尤度関数, 或いは罰則付き対数尤度関数と呼ばれます.
$R(\boldsymbol{\beta})$が2次形式で与えられる場合, $\boldsymbol{\beta}$と$\sigma ^2$の正則化最尤推定量は,
\begin{align}
\hat{\boldsymbol{\beta}}
&= (X^T X + \lambda \hat{\sigma}^2 K)^{-1} X^T \boldsymbol{y} \\
\hat{\sigma}^2
&= \frac{1}{n}(\boldsymbol{y} - \hat{\boldsymbol{y}})^T
(\boldsymbol{y} - \hat{\boldsymbol{y}})
\end{align}
で与えられます.
今回の場合, 正則化最尤法と, 正則化最小2乗法は,
$$
\lambda \sigma^2 = \gamma
$$
とすれば, 同じ$\hat{\boldsymbol{\beta}}$の推定量を得ることが分かります.
(c) リッジ回帰とlasso推定
$$
R(\boldsymbol{\beta}) = \boldsymbol{\beta}^T I_m \boldsymbol{\beta} = \boldsymbol{\beta}^T \boldsymbol{\beta} = \sum_{i=1}^{p+1} |\beta_i|^2
$$
は, リッジ型の推定量と呼ばれ, この考え方はパラメータベクトル$\boldsymbol{\beta}$の長さを一般化した$L_q$ノルムを用いた正則化最小二乗法という概念に拡張されます.
$q=2$の場合はリッジ回帰, $q=1$の場合はlasso(least absolute shrinkage and selection operator)推定と呼ばれます.
lasso推定はその制約の性質からパラメータの一部を完全に0とすることが可能であるため, 変数選択の手法として利用されます.
そのため今回はlasso推定による変数選択を行いたいと思います.
回帰係数の絶対値からなる$L_1$ノルムは微分不可能なため, アルゴリズムを用いて計算することになり, sklearnを用いて実施します.
はじめにsklearn.linear_models.lasso_path
を用いて, 正則化パラメータを増加させたときに, 偏回帰係数が0になる順番を確認したいと思います.
以下がPythonコードとその実行結果です.
# Standardize data using StandardScaler
X_scaled = StandardScaler().fit_transform(X)
# Compute optimized alpha with LassoCV
lasso_cv_model = LassoCV(
cv=20,
random_state=42
).fit(X_scaled, y)
alpha = lasso_cv_model.alpha_
# Compute paths for Lasso
eps = 5e-3
alphas_lasso, coefs_lasso, _ = lasso_path(
X_scaled, y, eps=eps
)
# Find the last nonzero alphas for each feature
last_nonzero_alphas = np.max(np.where(coefs_lasso != 0, alphas_lasso, 0), axis=1)
sorted_indices = np.argsort(last_nonzero_alphas)[::-1]
# Use seaborn color palette for better visualization
colors = sns.color_palette("husl", n_colors=X.shape[1])
log_alphas_lasso = np.log10(alphas_lasso)
# Plot Lasso coefficients for each feature
for coef_l, c, label in zip(
coefs_lasso[sorted_indices],
colors,
X.columns.values[sorted_indices]
):
plt.plot(log_alphas_lasso, coef_l, c=c, label=label)
# Plot vertical line at the estimated alpha with LassoCV
plt.axvline(
np.log10(alpha),
c="black",
linestyle="--",
label="estimated_alpha with LassoCV"
)
# Set plot labels and title
plt.xlabel("Log10(alpha)")
plt.ylabel("Coefficients")
plt.title("Lasso for Standardized X")
# Move the legend outside the plot for better visibility
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.axis("tight")
上図の横軸ラベルやコード内の'alpha'は正則化パラメータに相当し, プロットは横軸の$\alpha$の値に対応するそれぞれの特徴量の回帰係数の推定値です.
x軸と垂直なプロットは, sklearn.linear_models.lassocvを用いて最適化された$\alpha$を表しています.
グラフ右側の凡例の表示順序は, 回帰係数が0に収束するのが遅いものから降順としました.
また, グラフタイトルをご覧のとおり, 特徴量は標準化した上で計算しています.
上図から、$\alpha$を増加させても, 係数が0になりにくいものを上から5つ挙げると,
- year_construction
- walking_minute
- service_charge_yen
- area_square_meter
- renovated
であることが分かります.
特徴量の標準化の有無が異なるので単純比較して良いのか定かではありませんが, 先ほど確認した最小のAICを与えるモデルとは, 選択された特徴量の種類が少し異なる結果になりました.
sklearn.linear_models.lassocvクラスのalphas_, mse_path_という属性を使用して, 正則化パラメータ$\alpha$選定の過程を可視化することが可能です.
以下がPythonコードとその実行結果です.
# Number of cross-validation folds
num_folds = 20
# Fit LassoCV with specified number of folds
lasso_cv_model = LassoCV(
cv=num_folds,
random_state=42
).fit(X_scaled, y)
# Get the alphas used in cross-validation
alphas = lasso_cv_model.alphas_
# Plot the root mean squared error (RMSE) for each fold
for i in range(num_folds):
plt.plot(
np.log10(alphas),
np.sqrt(lasso_cv_model.mse_path_[:, i]),
c="grey",
label=f"fold {i+1}"
)
# Plot the average RMSE across all folds
plt.plot(
np.log10(alphas),
np.sqrt(lasso_cv_model.mse_path_.mean(axis=1)),
label='average RMSE of 20 folds'
)
# Plot vertical line at the estimated alpha with LassoCV
plt.axvline(
np.log10(lasso_cv_model.alpha_),
c="r",
linestyle="--",
label="estimated_alpha with LassoCV"
)
# Set plot labels and title
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.xlabel('Log10(alpha)')
plt.ylabel('Square Root of Mean Squared Residuals (RMSE)')
plt.title('Alphas and RMSRs')
# Set y-axis limit for better visualization
plt.ylim(0, 2e7)
# Show the plot
plt.show()
グレー色のプロットは, 横軸の$\alpha$の値に対応するそれぞれのフォールドのMSRの平方根です.
さらに, 青色のプロットは全てのフォールドの平均です.
赤色のプロットは, (sklearn.linear_models.lassocvを用いて)最適化された$\alpha$を表しています.
最適化された$\alpha$では, 今回の場合全てのフォールドの平均値としてのMSRが最小になっているように見えます.
先ほどのクロスバリデーションの手法を用いて, 同様の結果になるか確認してみます.
以下がPythonコードとその実行結果です.
# Initialize a list to store residuals for Leave-One-Out Cross Validation
loocv_residuals = []
# Fit LassoCV to determine optimal alpha
lasso_cv_model = LassoCV(
cv=20,
random_state=42
).fit(X_scaled, y)
optimal_alpha = lasso_cv_model.alpha_
# Initialize StandardScaler and fit to the entire dataset
scaler = StandardScaler()
scaler.fit(X)
# Perform Leave-One-Out Cross Validation
for i in range(len(df)):
# Split the dataset into test and training data
test_X = X.iloc[[i]]
train_X = X.drop(index=i)
test_y = y.iloc[[i]]
train_y = y.drop(index=i)
# Scale the training data
train_X_scaled = scaler.transform(train_X)
# Fit Lasso with the optimal alpha to the training data
lasso_model = Lasso(alpha=optimal_alpha).fit(train_X_scaled, train_y)
# Predict using the scaled test data
y_pred = lasso_model.predict(scaler.transform(test_X))
# Calculate residuals
residual = y_pred - test_y
loocv_residuals.append(residual)
# Calculate mean squared residuals for Leave-One-Out Cross Validation
msr = np.sum(np.square(loocv_residuals)) / len(df)
sqrt_msr = np.sqrt(msr)
print(f"Square root of mean squared residuals: {sqrt_msr:,.0f}")
# Calculate mean absolute residuals for Leave-One-Out Cross Validation
mar = np.sum(np.abs(loocv_residuals)) / len(df)
print(f"Mean absolute residuals: {mar:,.0f}")
"""
Square root of mean squared residuals: 4,590,545
Mean absolute residuals: 3,997,761
"""
上記のとおり, RMSRが4,590,545と, 同じような結果になりました.
注意点として, 上記のコードと, sklearn.linear_models.lassocvのいずれについても, 分割前の特徴量行列に対して標準化を実施しています.
試しに, 上記のコードで, 分割ごとに特徴量の標準化を行ったところクロスバリデーションによるMSRの数値が大幅に増加しました.
どのタイミングで標準化を行うかについては, とりわけモデルを比較する際には注意を要しそうです.
クロスバリデーションによるMSRを評価指標とした比較では, ガウス型線形回帰モデルの総当たりによる特徴量選択の結果の方が, sklearn.linear_models.lassocvを用いて最適化された正則化パラメータによるlasso推定モデルの結果よりも優れているという結果になりました.
なお, sklearn.linear_models.lassocvの他に, 正則化パラメータ$\alpha$の最適化のために, sklearnには以下のようないくつかのクラスが用意されています.
sklearn.linear_models.lassolarscvsklearn.linear_models.lassolarsic
今回はいずれの手法を用いても, ガウス型線形回帰モデルの総当たりによる特徴量選択の結果を上回る評価指標を得ることができませんでした.
それぞれの手法により最適化された$\alpha$を取得するためのPythonコードとその実行結果は以下のとおりです.
reg = LassoCV(
cv=20,
random_state=42
).fit(X_scaled, y)
print(f"Optimized alpha with LassoCV: {np.log10(reg.alpha_)}")
reg = LassoLarsCV(
cv=20
).fit(X_scaled, y)
print(f"Optimized alpha with LassoLarsCV: {np.log10(reg.alpha_)}")
reg = LassoLarsIC(
criterion='aic'
).fit(X_scaled, y)
print(f"Optimized alpha with LassoLarsIC using AIC: {np.log10(reg.alpha_)}")
reg = LassoLarsIC(
criterion='bic'
).fit(X_scaled, y)
print(f"Optimized alpha with LassoLarsIC using BIC: {np.log10(reg.alpha_)}")
print("Values obtained by taking the common logarithm are displayed as alphas.")
"""
Optimized alpha with LassoCV: 5.628925368029047
Optimized alpha with LassoLarsCV: 5.629822044608447
Optimized alpha with LassoLarsIC using AIC: 5.572488495533479
Optimized alpha with LassoLarsIC using BIC: 5.572488495533479
Values obtained by taking the common logarithm are displayed as alphas.
"""
まとめ
-
ガウス型線形回帰モデルと対応する最小二乗法による線形回帰モデルの結果が同等であることを
sklearnを用いて確認した. - AICを用いた総当たりによる線形回帰モデルに使用する特徴量の探索を行ない, 探索したモデルを所定のクロスバリデーションによるMSRにより評価した.
- 所定の方法で最適化された正則化パラメータを用いたlassoモデルを上述のMSRにより評価した.