オッズ(odds)
- オッズは、競馬では「倍率」と同義で、馬券が的中すると、オッズに応じて払い戻しが行われます。かりに10倍のオッズがついていたら、賭け金が10倍になって返ってくるわけですが･･･
- 統計においては、オッズとは、ある事象が起きる確率(P)と起きない確率(1-P)との比のことで、求め方は$\displaystyle \frac{P}{1-P}$です。
- たとえば、薬品ＡとＢがあって、それぞれの治験の結果を「効果あり」「効果なし」という２つの事象として、各オッズを示します。
| 効果あり P | 効果なし 1-P | オッズ P/(1-P) | |
|---|---|---|---|
| 薬品Ａ | 0.2 | 0.8 | 0.250 |
| 薬品Ｂ | 0.05 | 0.95 | 0.053 |
オッズ比(odds ratio)
- オッズ比とは、２つのオッズの比をとったものです。薬品Ａと薬品Ｂのオッズ比は、Ａのオッズ÷Ｂのオッズで、$0.250/0.053=4.75$となります。
- オッズ比が大きいほど（または小さいほど）事象の間に関連があることを意味し、オッズ比が1.0だとすると、事象間の関係が独立していて互いに影響していないことになります。さらに･･･
- オッズ比の対数をとったものが対数オッズ比で、回帰係数は、説明変数を1単位変化させた時の対数オッズ比となります。
ロジスティック回帰におけるオッズ比と回帰係数との関係
を考察するために、勉強時間数からテストの合否（合格ならば1, 不合格ならば0）を予測するモデルを仮定します。
⑴ ライブラリのインポート
# 数値計算に使うライブラリ
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
# グラフを描画するライブラリ
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
# 統計モデルを推定するライブラリ
import statsmodels.formula.api as smf
import statsmodels.api as sm
# 表示桁数の指定
%precision 3
- scipyは、高度な科学計算を行うためのライブラリで、statsという統計関数をまとめたモジュールをもっています。
- https://docs.scipy.org/doc/scipy/reference/stats.html
- seabornは、matplotlibをベースとする可視化ライブラリで、matplotlibのグラフをより表現力のあるものにしてくれます。
- https://seaborn.pydata.org/introduction.html
- statsmodelsは、ずばりその名の通り、いろいろな統計モデルを推定するためのクラスおよび関数を提供するモジュールです。
- https://www.statsmodels.org/devel/
⑵ データの読み込みと確認
# データの取得
url = 'https://raw.githubusercontent.com/yumi-ito/sample_data/master/6-3-1-logistic-regression.csv'
# データの読み込み
df = pd.read_csv(url)
# # データの先頭5行を出力
df.head()
# データの基本統計量を出力
df.describe().apply(lambda s: s.apply(lambda x: format(x, 'g')))
- 説明変数は、勉強時間(hours)のみで「0, 1, 2, ･･･, 8, 9時間」の全10項目あります。
- 目的変数はテストの合否(result)で、平均が0.46ですから100人中46人が合格しています。
# 棒グラフを描画
sns.set()
sns.barplot(x = "hours", y = "result", data = df, palette='summer_r')
-
sns.set()はグラフ領域のスタイルを指定、デフォルトの状態です。あるいはsns.set(style="whitegrid")とすると白地に灰色の目盛線になります。 -
seaborn.barplot(x列名, y列名, DataFrame名, カラーパレット名)という引数になっています。 - デフォルトで95％信頼区間を示すエラーバー（誤差範囲）が描画されます。引数に
ci=Noneを加えれば非表示になります。 - https://seaborn.pydata.org/generated/seaborn.barplot.html
- 勉強時間が長いほど合格率が高くなる傾向があり、また、5時間以下か6時間以上かで合格率の水準が大きく変動しています。
- 勉強時間ごとの合格率を計算してみます。
pass_rate = df.groupby("hours").mean()
- pandasの関数
groupby()をつかって引数に指定した列で同一要素をグループ化し、それぞれの平均を.mean()で求めます。 - 勉強時間が6時間以上になると合格率が飛躍的に上がっています。
⑶ モデルの推定と結果の確認
# モデルを推定
mod_glm = smf.glm(formula = "result ~ hours",
data = df,
family = sm.families.Binomial()).fit()
- statsmodelsで一般化線形モデルを推定する場合は、
smf.glm()関数を使います。glmとはGeneralized Linear Modelsの略。 - １つ目の引数
formulaでモデルの構造を、目的変数がresult、説明変数がhoursと指定します。 - ３つ目の引数
familyは確率分布の指定で、この例は二項分布(Binomial Distribution)なのでsm.families.Binomial()を適用します。
# 推定結果のサマリーを出力
mod_glm.summary()
# 回帰曲線を描く
sns.lmplot(x = "hours", y = "result",
data = df,
logistic = True,
scatter_kws = {"color": "green"},
line_kws = {"color": "black"},
x_jitter = 0.1, y_jitter = 0.02)
- 引数
logisticはTrueの場合、yをバイナリ変数と仮定し、statsmodelsを使用してロジスティック回帰モデルを推定します。バイナリ変数とは0と1の2値のみをとり得る変数のことです。 - 引数
scatter_kwsとline_kwsは、plt.scatterおよびplt.plotに渡す追加のキーワード引数で、散布図のドットや回帰曲線の色指定をしています。 - 引数
x_jitterとy_jitterは、ドットを少し上下に散らばせるという指定で、単に見た目だけのことです。合否は1か0しかとらないため、ドットが重なってしまうのを制御しています。
- 勉強時間ごとの合格率を数値で取得してみます。
# 列名をhoursとする等差数列(0～9)のDataFrameを作成
predicted_value = pd.DataFrame({"hours": np.arange(0, 10, 1)})
# 合格率の予測値を算出
pred = mod_glm.predict(predicted_value)
- モデル
mod_glmに対し、関数predict()で、作成したデータフレームpredicted_valueに合わせて予測値を算出させます。
⑷ 対数オッズ比を求め、係数と比較してみる
# １時間と２時間の合格率を取得
pred_1 = pred[1]
pred_2 = pred[2]
# それぞれのオッズを計算
odds_1 = pred_1 / (1 - pred_1)
odds_2 = pred_2 / (1 - pred_2)
# 対数オッズ比を算出
print("対数オッズ比：", round(sp.log(odds_2 / odds_1), 3))
# モデルの係数を算出
value = mod_glm.params["hours"]
print("モデルの係数：", round(value, 3))
- 対数オッズ比は、回帰係数と一致します。
- また、回帰係数のexp(指数関数)をとると、オッズ比に一致します。
-
numpy.exp(x)は自然対数$e$の$x$乗を返します。
# 回帰係数のexpをとる
exp = sp.exp(mod_glm.params["hours"])
print("係数のexp：", round(exp, 3))
# オッズ比を算出
odds = odds_2 / odds_1
print("オッズ比：", round(odds, 3))
- つまり、説明変数が1だけ変化するとき、目的変数はオッズ比倍だけ変化することになります。
- 回帰係数よりも、オッズ比に変換したほうが、納得しやすかったり、説明しやすかったりする場合もあるかと思います。