雑記

仮に愛知県民の 20% がきのこの山の方が好きで、残る 80% はたけのこの里の方が好きだとする。
このとき、「きのこの山が好きである」のオッズは 0.2 / 0.8 = 0.25 になる。
さらに、「きのこの山が好きである」のロジットは log(0.25) = -1.39 になる。
これらはとりもなおさず以下が成り立つことを意味する。

 〈愛知県のきのこ派人口〉= 0.2〈愛知県の人口〉
 〈愛知県のきのこ派人口〉= 0.25〈愛知県のたけのこ派人口〉
 〈愛知県のきのこ派人口〉= exp( -1.39 )〈愛知県のたけのこ派人口〉

仮に静岡県では 60% の県民がきのこの山の方が好きであったとすると、このオッズとロジットは以下になる。

 〈静岡県きのこ派人口〉= 0.6〈静岡県の人口〉
 〈静岡県きのこ派人口〉= 1.5〈静岡県たけのこ派人口〉
 〈静岡県きのこ派人口〉= exp( 0.405 )〈静岡県たけのこ派人口〉

このとき、「愛知県民であること」が「きのこの山を好きである」要因になっているかを知りたい
(ただし世界には愛知県民と静岡県民しかいないものとする)。
そこで、個々人がきのこ派であるかのロジットを予測する線形回帰モデルを考える(ロジスティック回帰)。
つまり以下のモデルを考える。

 〈その人がきのこ派であるロジット〉= α + β〈その人が愛知県民であるか否か〉

この β の最尤推定量は -1.7918 になる(本記事末尾のスクリプト)。スクリプトでは愛知県の人口と静岡県の人口をつかっているが、何人ずつでも構わない(無論 p 値には影響する)。なぜなら静岡県民のきのこ派率のみが切片 α を決定し、その上で愛知県民のきのこ派率が β を決定するからである。

 〈愛知県民がきのこ派であるオッズ〉= exp(α + β)
 〈静岡県民がきのこ派であるオッズ〉= exp(α)

さらにこれより β の値 -1.7918 とは、愛知県民のオッズと静岡県民のオッズの比が exp(-1.7918) という意味になる。
実際 exp(-1.7918) とオッズ比は等しく 0.166666 である。

print(np.exp(-1.7918))  # 0.166666
print(0.25 / 1.5)  # 0.166666

逆にいうとオッズ比さえ 0.166666 ならば β は -1.7918 になる。
愛知県のきのこ派率が 0.333 で静岡県きのこ派率が 0.75 でも β は -1.7918 になる。
愛知県のきのこ派率が 0.11111 で静岡県きのこ派率が 0.42857 でも β は -1.7918 になる。
愛知県のきのこ派率が 0.02439 で静岡県きのこ派率が 0.130438 でも β は -1.7918 になる。
どれも「愛知県民であることが、きのこ派であるオッズを 0.166666 倍にする」ということになる。

しかし、「オッズを exp(β) 倍にする」といわれても解釈しにくいものである。
この例では「愛知県民であること」によって「きのこ派である確率」が削られる幅が 0.6 - 0.2 = 0.4 だが、
同じ β=-1.7918 でも削られる幅が 0.75 - 0.333 = 0.417 であることもあるし、
0.130438 - 0.02439 = 0.106048 であることもある。
これは元々愛知県民でないときにきのこ派である確率の水準がどれくらいかに依存する。
なので異なる応答変数に対するロジスティック回帰の β の値で効果の大きさを比較できないのはもちろん、
同じデータでの同じ応答変数に対するロジスティック回帰であっても異なる説明変数を含めたケースとは比較できないはずである。


import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

n_aichi = 7500000
n_shizuoka = 3700000
p_aichi = 0.2
p_shizuoka = 0.6

df = pd.DataFrame({
    'n': [
        p_aichi * n_aichi, (1 - p_aichi) * n_aichi,
        p_shizuoka * n_shizuoka, (1 - p_shizuoka) * n_shizuoka
    ],
    'kinoko': [1, 0, 1, 0],
    'aichi': [1, 1, 0, 0],
})
family = sm.families.Binomial()
result = smf.glm(
    formula='kinoko ~ aichi',
    data=df,
    family=sm.families.Binomial(),
    freq_weights=np.asarray(df['n']),
).fit()
print(result.summary())
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.4055      0.001    382.085      0.000       0.403       0.408
aichi         -1.7918      0.001  -1280.005      0.000      -1.795      -1.789
==============================================================================