ガウス型線形回帰モデルによる中古マンション価格の分析

# はじめに 
本稿では, 中古マンションの売り出し価格に寄与する変数を, 線形回帰モデルを用いて確認します. 

断言調になっている箇所のほとんどは, 調べたことをまとめているに過ぎないので予めご了承ください. 


# 使用したデータ 

使用したデータは非公開で, 内容は以下のとおりです. 

- データ数: 20件　
- 関心の駅名で検索をかけるとヒットする中古マンションの1室が1つのレコードに対応
- 目的変数は売り出し中の販売価格　
- 不動産情報サイトを参考に筆者が手作業で作成

# 環境　
- Google Colaboratory 

## 関連ライブラリのインポート 
```python　
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 
``` 

## データの読み込み 
```python 
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)
"""
``` 

## データの加工 
そのままでは使いづらいかもしれないので, 以下のコードで簡単なクリーニングを行いました. 
```python 
# 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にそれぞれ格納します. 

```python 
# 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,))
"""
``` 

# 本編　
以下の線形モデルを仮定します. 

```math
\boldsymbol{y} = X\boldsymbol{\beta} + \epsilon
``` 
```math
\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)最小二乗法　
**最小二乗法**は誤差の二乗和を最小にするような回帰式を求める方法. 

誤差について以下の仮定を置くことにより, 良い推定量であることが保証されます. 

```math 
\begin{align}
E[\boldsymbol{\epsilon}] 
    &= \boldsymbol{0}^T \\
E[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^T] 
    &= \sigma^2 I_n \\
\end{align}
``` 
```math 
\begin{align}
\boldsymbol{0}
    &: 要素がすべて0のn次元ベクトル \\
I_n
    &: n \times n 単位行列\\
\sigma^2
    &: 誤差分散 \\ 
\end{align}
``` 


それぞれの観測点における誤差の期待値及び分散が全て0, $\sigma^2$で, 且つ互いに無相関ということになります. 

(今回の目的変数は中古マンションの売り出し価格で, その数値の増加に伴いばらつきも増加しそうですので, 直感として上記の仮定は誤りなるものと理解)


また, 誤差の二乗和は, 

```math 
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}$は, 

```math 
\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}}$は, 

```math 
\hat{\boldsymbol{\beta}} 
= (X^{T} X)^{-1} X^{T} \boldsymbol{y}
``` 
となります. 

**予測値ベクトルy$, $残差ベクトルe**をそれぞれ, 

```math 
\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}$)であることに注意すると, **残差平方和**は, 

```math 
\boldsymbol{e}^{T} \boldsymbol{e} 
    = \boldsymbol{y}^{T} [I_{n} - X (X^{T} X)^{-1} X^{T}] \boldsymbol{y}
``` 

のように表されます. 

最小二乗推定量$\boldsymbol{\hat{\beta}}$の期待値, 分散共分散行列は, 

```math
\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()`を使用して, モデルをフィットし, サマリーを表示することが可能. 

```python 
# 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)最尤法　
**最尤法**は今回(**ガウス型線形回帰モデル**)の場合, 目的変数が正規分布に従って観測されたと仮定し, 尤度関数を最大化するようなモデルを学習します. 

```math
\begin{align}

\boldsymbol{y} = X \boldsymbol{\beta} + \boldsymbol{\epsilon} 

\qquad 

\boldsymbol{\epsilon} \sim N(0, \sigma^2)

\end{align}
```

上式のように, 誤差項の従う分布を明示します. 

所定の対数尤度関数を最大化するための尤度方程式を解くことによって, $\beta$と$\sigma$の最尤推定値はそれぞれ, 

```math 
\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()`を使用して, モデルをフィットし, サマリーを表示することが可能. 

```python 
# 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のコードです. 

```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コードと実行結果です. 

```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()`という関数があるようですが, シンプルに以下のコードで計算しました. 

```python 
# 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})$を加えた関数

```math
\begin{align}
S_{\gamma}(\boldsymbol{\beta})
    &= (\boldsymbol{y} - X \boldsymbol{\beta})^{T} 
        (\boldsymbol{y} - X \boldsymbol{\beta})
        + \gamma R(\boldsymbol{\beta}) \\
\end{align}
```

```math
\begin{align}
\gamma(>0)
    &: 正則化パラメータ(或いは平滑化パラメータ)
\end{align}
```

の最小化を行います. 

ペナルティ項$R(\boldsymbol{\beta})$は2次形式 

$$
\boldsymbol{\beta}^{T} K \boldsymbol{\beta}
$$

で表現することができる場合が多く, 

```math
\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}
```
```math
\begin{align}
K 
    &: 非負値定符号行列
\end{align}
```

とすると, 

```math
\hat{\boldsymbol{\beta}} 
    = (X^T X + \gamma K)^{-1} 
        X^T \boldsymbol{y}
``` 

と, 解析的に陽に求めることが可能. 

#### (b) 正則化最尤法 

**ガウス型線形回帰モデル**の対数尤度関数は, 

```math 
\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}
``` 

上式に最尤推定量を代入すると, 

```math 
l(\boldsymbol{\hat{\theta}}) 
    =  - \frac{n}{2}log(2\pi)
        - \frac{n}{2}log\hat{\sigma}^2 
        - \frac{n}{2} 
``` 
分散の最尤推定値は残差平方和から計算されるため, 対数尤度の値はモデルの複雑さに伴い増加します. 

そこで, **正則化最小二乗法**と同様に, ペナルティ項を導入し, 

```math 
\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$の**正則化最尤推定量**は, 

```math
\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コードとその実行結果です. 

```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")
``` 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/469ceb54-2199-f6f5-3249-a593298a4fb8.png)

上図の横軸ラベルやコード内の'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コードとその実行結果です. 

```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()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/26cca1dc-6b6c-4ac0-558e-647c5546ae6e.png)


グレー色のプロットは, 横軸の$\alpha$の値に対応するそれぞれのフォールドの**MSR**の平方根です. 

さらに, 青色のプロットは全てのフォールドの平均です. 

赤色のプロットは, (`sklearn.linear_models.lassocv`を用いて)最適化された$\alpha$を表しています.  


最適化された$\alpha$では, 今回の場合全てのフォールドの平均値としての**MSR**が最小になっているように見えます. 

先ほどの**クロスバリデーション**の手法を用いて, 同様の結果になるか確認してみます. 

以下がPythonコードとその実行結果です. 

```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.lassolarscv`
- `sklearn.linear_models.lassolarsic`

今回はいずれの手法を用いても, ガウス型線形回帰モデルの総当たりによる特徴量選択の結果を上回る評価指標を得ることができませんでした. 

それぞれの手法により最適化された$\alpha$を取得するためのPythonコードとその実行結果は以下のとおりです. 

```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**により評価した.