「入門はじめての多変量解析」をPythonで写経 Vol.8 ～ ２章「はじめての重回帰分析」⑦ダミー変数、変数変換

２章「はじめての重回帰分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」２章「はじめての重回帰分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、ダミー変数 と 変数変換 に取り組みます。
何と言うか…自由度が大きい記事になりました！

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

２章 はじめての重回帰分析

この記事は２章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import statsmodels.formula.api as smf

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

ダミー変数をはじめる

ダミー変数は変数の変換手法の１つです。
質的変数（カテゴリ変数）を０または１の数値に変換した変数です。
機械学習の文脈ではこの変換を One-Hot エンコーディングと呼んだりもします。

重回帰分析は「数値」の変数を扱います。
質的変数は「名称等の文字列」や「加減乗除の意味がない数値」だったりするので、

の計算に当てはめられるように数値に変換します。

🎬 ダミー変数の例

$$
\begin{array}{l:l|l:l}
変数 & 要素 & ダミー変数 & ダミー変数の要素(意味) \\
\hline
\\
性別 & 男性,\ 女性 & 性別 & 0\ (男性),\ 1\ (女性) \\
\\
色 & 赤,\ 青,\ 緑 & 色\_青 & 0\ (青以外),\ 1\ (青) \\
 & & 色\_緑 & 0\ (緑以外),\ 1\ (緑) \\
\\
回答 & 好き,\ 普通,\ 嫌い & 回答\_好き & 0\ (好き以外),\ 1\ (好き) \\
 & & 回答\_嫌い & 0\ (嫌い以外),\ 1\ (嫌い) \\
\end{array}
$$

🎬 変数「色」の場合のダミー変数の設定（例）

$$
\begin{array}{c|c:c}
変数「色」 & ダミー変数「色\_青」 & ダミー変数「色\_緑」 \\
\hline
赤 & 0 & 0 \\
青 & 1 & 0 \\
緑 & 0 & 1 \\
\end{array}
$$

🎬 ダミー変数のつくり方
要素が３つ以上の場合、多重共線性対策として 要素を１つ、ダミー変数化から外します。

• 変数に含まれる要素が２つの場合
  ⇒ $${0}$$、$${1}$$ に変換

• 変換に含まれる要素が３つの場合
  ⇒ 要素のうち２つを選んで、２つの変数を作成
  ⇒ 変数に対応した要素を $${1}$$、違う要素を $${0}$$ に変換

• 変換に含まれる要素が $${n}$$ 個の場合
  ⇒ 要素のうち $${n-1}$$ 個を選んで、$${n-1}$$ 個の変数を作成
  ⇒ 変数に対応した要素を $${1}$$、違う要素を $${0}$$ に変換

ダミー変数を重回帰分析で使う

テキストの「銀行員 265 人の給与データ」をお借りしたかったのですが、テキストに 265 行の全レコードの記載がなく、Webサイトでのデータ配布もなさそうでした。

🏠 引用データの概要

そこで、以下の家賃データをお借りします。
ありがとうございます！

【引用表示】
Webサイト：
千葉商科大学 商経学部 情報コース（永岡） 2024年度 応用情報処理
第9回 回帰分析(3)―重回帰分析 アパートの家賃(1)
URL：https://www3.cuc.ac.jp/~nagaoka/2011/ouyou/11/rent/index.html
データ：市川物件データ（rent.xlsx）

（参考）この物件データを用いた重回帰分析の解説サイト

■ データの確認

① データの読み込み
WebサイトのExcelファイルを直接、pandas データフレームに読み込みます。

### データの読み込み

# ExcelファイルのURLの設定
url = r'https://www3.cuc.ac.jp/~nagaoka/2011/ouyou/11/rent/rent.xlsx'

# WebサイトからExcelファイルを直接読み込み
data4 = pd.read_excel(url, index_col=0)

# データの概要
print('data4.shape: ', data4.shape)
data4.head()

【実行結果】
297 件の千葉県市川市の家賃データです。
目的変数は「家賃」です。家賃を予測する回帰モデルを作ります。

② データの概観の確認
データの型と欠損値の有無を確認します。

### データの概観の確認

# データの型の確認
data4.info()

【実行結果】
量的変数は int64 型、float64 型 の変数です。
質的変数は object 型の変数です。

# 欠損値の確認
data4.isnull().sum().rename('欠損値の数').to_frame().T

【実行結果】
欠損値はありません。

■ データの分析

③ 基礎的な統計量の確認
量的変数の要約統計量を表示します。

### データの基礎統計量の確認

# データの要約統計量の確認（量的変数のみ）
data4.describe().T.round(3)

【実行結果】
家賃の平均（および中央値）はおよそ６万円です。
面積の平均（および中央値）はおよそ 20 ㎡。
築年数の平均は 16年程度、駅徒歩の平均は８分程度。

量的変数間の相関係数を表示します。

# データの相関係数の確認（量的変数のみ）
data4.corr(numeric_only=True).round(3)

【実行結果】
家賃との相関係数の絶対値が高いのは、築年数と面積です。
説明変数間の相関係数は高くないので多重共線性は無さそうです。

質的変数の要約統計量を表示します。

# データの要約統計量の確認（質的変数のみ）
data4.describe(include='object').T

【実行結果】
unique は変数に含まれる要素種類の数です。
所在地は 10 箇所あり、新田が最も多く 100 個あります。

質的変数の要素名とその個数をカウントしましょう。

# 質的変数のユニーク要素の確認
for col in data4.columns[data4.dtypes=='object']:
    display(data4[col].value_counts().to_frame())

【実行結果】
新田、1K、南向き、木造が圧倒的に多い印象です。

④ 目的変数と説明変数の関係の可視化
説明変数が量的変数の場合は散布図、質的変数の場合はヴァイオリンプロットで目的変数との関係を確認します。
なお、量的変数のうち「カテゴリ別」のほうが見やすい変数は、質的変数の方に含めます。

変数を量的変数と質的変数に層別して、目的変数 vs 量的変数の散布図を描きます。

### 目的変数との関係性の可視化

# 変数の種類別にリスト化
TARGET = '家賃'
num_cols = ['面積', '築年数', '駅徒歩分']
cat_cols = data4.columns[~data4.columns.isin([TARGET]+num_cols)]

# 目的変数 vs 量的変数：散布図
fig, axes = plt.subplots(1, 3, figsize=(10, 3), sharey=True, tight_layout=True)
for col, ax in zip(num_cols, axes.flat):
    sns.regplot(data=data4, x=col, y=TARGET, ax=ax,
                scatter_kws={'s': 30, 'ec': 'white'})
    ax.set_title(col)
plt.show()

【実行結果】
面積と家賃は正の相関、築年数と家賃は負の相関がありそうです。
駅徒歩分と家賃も薄っすらと負の相関があるかもです。

目的変数と質的変数のヴァイオリンプロットを描きます。

# 目的変数 vs 質的変数：ヴァイオリンプロット
fig, axes = plt.subplots(2, 3, figsize=(10, 6), sharey=True, tight_layout=True)
for col, ax in zip(cat_cols, axes.flat):
    if col == '所在地':
        ax.xaxis.set_tick_params(rotation=90)
    sns.violinplot(data=data4, x=col, y=TARGET, fill=False, alpha=0.7, ax=ax)
    ax.set_title(col)
plt.show()

【実行結果】
なんだか微妙な、相関関係が分かりにくい感じです。。。
ひとまず「間取」と「方位」を説明変数に残してみましょう。

■ 重回帰モデルの構築

⑤ データの前処理
説明変数を「面積」「築年数」「駅徒歩分」「間取」「方位」に絞り込みます。
質的変数の「間取」「方位」をダミー変数にします。
pandas の get_dummies 関数で、質的変数をダミー変数に変換します。
処理結果を先頭５行表示します。

### データの前処理

## 説明変数の絞り込み
# 削除する説明変数の設定
drop_cols = ['敷金', '礼金', '所在地', '構造']

# 説明変数の絞り込み
data4_select = data4.drop(columns=drop_cols)

## カテゴリ変数をダミー変数化
data4_dummy = pd.get_dummies(data4_select, drop_first=True, dtype=int)

## 結果の確認
print('data4_dummy.shape: ', data4_dummy.shape)
data4_dummy.head()

【実行結果】

ダミー変数から除外した要素は以下のとおりです。
　・間取：1DK
　・方位：北 

⑥ 重回帰分析の実行

## 重回帰分析の実行
result4 = smf.ols(formula='家賃' + ' ~ ' + ' + '.join(data4_dummy.columns[1:]),
                  data=data4_dummy).fit()
result4.summary2()

【実行結果】

【モデル全体の状態】
決定係数 $${0.670}$$、自由度調整済み決定係数 $${0.657}$$ は、まずまずの当てはまりの良さだと思います。
$${F}$$ 値の $${p}$$ 値は $${< 0.05}$$ であり、有意水準 $${5\%}$$ で回帰モデルは有意と考えられます。

【偏回帰係数の有意性の検定】
有意水準 $${5\%}$$ で有意な偏回帰係数（切片を除く）は次のとおりです。
　・面積
　・築年数
　・駅徒歩分
　・間取_1K
　・間取_1R

ダミー変数では「間取」が有意となり、「方位」は有意ではありません。

【重回帰式の解釈】
有意な偏回帰係数に関して、説明変数が目的変数に与える影響（傾向）をまとめます。

【メモ：残差診断】

• 残差の自己相関に関しては、ダービン・ワトソン比（Durbin-Watson） $${2.134}$$ が２に近似しているので、残差の自己相関は無さそうです。

• 残差の正規性に関しては、正規性の検定（Omnibus、Jarque-Bera）がいずれも有意水準 $${5\%}$$ で有意であり、残差が正規分布に従っていない、と判断されます。

🔚 モデルの追求、今日はここまでにしましょう🍀

対数変換 p.82

次のテーマに進みます。

テキストの Section 2.11 で「量的変数の変数変換」が紹介されています。
仮想データを作成して、変数変換を行います。

■ 仮想データの作成
２つの変数はいずれも「峰が左側に偏り、右側に裾が長い分布」です。
実態は、２変量正規分布乱数を指数変換したものです。

### 仮想データの作成

## 設定と準備 
# ２変量正規分布のσ1, σ2, ρの設定
sigma1, sigma2, rho = 1, 0.8, 0.8
# ２変量正規分布の共分散の算出
cov = sigma1 * sigma2 * rho
# 乱数生成器
rng = np.random.default_rng(seed=3)

## ２変量正規分布のデータの生成
# ２変量正規分布の設定
dist = stats.multivariate_normal(
    mean=[0, 0], cov=[[sigma1**2, cov], [cov, sigma2**2]], seed=rng)
# ２変量正規分布乱数の生成 
x_orgn, y_orgn = dist.rvs(size=1000).T

## 乱数を指数化してデータフレーム化
data5 = pd.DataFrame({'x': np.exp(x_orgn), 'y': np.exp(y_orgn)})

## データの可視化
sns.jointplot(data=data5, x='x', y='y');

【実行結果】
０付近に偏っています（意図的に偏らせています）。

■ 対数変換
目的変数と説明変数の両方を対数変換します。
numpy の log() で自然対数を計算できます。

## 対数変換

# x,yを対数変換
data5['log_x'] = np.log(data5['x'])
data5['log_y'] = np.log(data5['y'])

# 可視化
sns.jointplot(data=data5, x='log_x', y='log_y');

【実行結果】
対数変換の結果、
周辺の（１変数の）ヒストグラムは正規分布に近づきました。
散布図も右上がりになっていて、強めの正の相関を示しています。

正規分布に近づけるための変換 p.83

仮想データのように右側に裾が長い分布の場合、変数を変換して正規分布に近づけることで、重回帰分布などのモデルで扱いやすくなるようです。

ここでは、Box-Cox 変換と Yeo-Johnson 変換を試してみます。

📊 Box-Cox 変換の公式
データ $${x}$$ をパラメータ $${\lambda}$$ のもとで $${x^{(\lambda)}}$$ に変換する Box-Cox 変換の公式です。
ただし、$${x}$$ が正（$${x > 0}$$）の場合のみ計算可能です。

仮想データに Box-Cox 変換を適用して可視化します。
scipy.stats の boxcox() を利用します。

## Box-Cox変換

# x,yをBox-Cox変換
data5['boxcox_x'] = stats.boxcox(data5['x'])[0]
data5['boxcox_y'] = stats.boxcox(data5['y'])[0]

# 可視化
sns.jointplot(data=data5, x='boxcox_x', y='boxcox_y');

【実行結果】
変換の結果、
周辺の（１変数の）ヒストグラムは正規分布に近づきました。
散布図も右上がりになっていて、強めの正の相関を示しています。

📊 Yeo-Johnson 変換
変換元データが負の場合も扱える変換です。

scipy の公式サイトより、Yeo-Johnson 変換の公式を引用いたします。

y = ((x + 1)**lmbda - 1) / lmbda,                for x >= 0, lmbda != 0
    log(x + 1),                                  for x >= 0, lmbda = 0
    -((-x + 1)**(2 - lmbda) - 1) / (2 - lmbda),  for x < 0, lmbda != 2
    -log(-x + 1),                                for x < 0, lmbda = 2

負の値が含まれる仮想データを作成します。
最初の仮想データから１を差し引くだけの操作です。

## 負のデータの作成 ※BoxCox変換は負のデータが含まれていると適用できない。

# x,y から1を引いたデータを作成
data5['x_minus'] = data5['x'] - 1
data5['y_minus'] = data5['y'] - 1

# 可視化
sns.jointplot(data=data5, x='x_minus', y='y_minus');

【実行結果】
x軸、y軸ともに負のデータ点が含まれました。
分布は相変わらず、右側に裾が長い状態を継続しています。

負の値を含む仮想データに Yeo-Johnson 変換を適用して可視化します。
scipy.stats の yeojohnson() を利用します。

## Yeo-Johnson変換

# x,yをYeo-Johnson変換
data5['yeojohnson_x'] = stats.yeojohnson(data5['x_minus'])[0]
data5['yeojohnson_y'] = stats.yeojohnson(data5['y_minus'])[0]

# 可視化
sns.jointplot(data=data5, x='yeojohnson_x', y='yeojohnson_y');

【実行結果】
変換の結果、
周辺の（１変数の）ヒストグラムは正規分布に近づきました。
散布図も右上がりになっていて、やや強めの正の相関を示しています。

その他の話題のメモ

話題を変えます。

テキストが「いろいろな回帰分析」の見出しで、以下の回帰分析の名称を紹介しています。
・ロジスティック回帰分析
・カテゴリカル回帰分析
・順序回帰分析
・名義回帰分析
・非線形回帰分析
・時系列データの回帰分析（自己回帰モデルを利用します）

ただし、各回帰分析の内容には踏み込んでいません。

【ご参考】
上に列挙した「ロジスティック回帰分析やカテゴリカル回帰分析（多項ロジスティック回帰）」は目的変数が質的変数のケースを扱います。
機械学習の文脈では「分類」タスクと呼んだりします。

以下のブログで分類タスクの機械学習を取り扱っています。
ロジスティック回帰を含む７つのモデルの Python コード例を掲載していますので、よかったらお読みください！

記事の最後はChatGPTが締めくくります。
今回は新しい息吹に心を向けて。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。