Rで多重ロジスティック＋交互作用をみる

2023年12月14日2023年12月18日

今回の話はRでロジスティック回帰の続きになります。

前回は身長データから学年が小学３年生か小学６年生かを判別できるか、ということでロジスティック回帰をやってみました。

もともとのデータは男女別々に作ったものを合算していました。今回はその性別の情報を用いて精度を上げられないかということを考えます。

今回はいくつかの試行錯誤が入りますので、データがころころ変わってしまうと評価が難しくなってしまいます。そこで乱数から生成されるデータを固定します。

set.seed(12345)

#データの作成
e3b = rnorm(50, 128.3, 5.48) #小学３年男子
e3g = rnorm(50, 127.6, 5.68) #小学３年女子
e6b = rnorm(50, 145.9, 7.27) #小学６年男子
e6g = rnorm(50, 147.3, 6.47) #小学６年女子

set.seedの中の数字を乱数の「種」といいます。これを変えると生成されるデータが変化します。今回はあまりロジスティック回帰での正解率が高くないものを選びました。

次にデータの代表値を計算しておきます。これは回帰の結果を解釈する際に使います。

#代表値（平均・中央値）の出力
m_height = mean(df_all$Height)
m_height_b = mean(df_all[df_all$Gender==0,]$Height)
m_height_g = mean(df_all[df_all$Gender==1,]$Height)

s <- sprintf("mean of height=%.1f (boys=%.1f, girls=%.1f",m_height, m_height_b, m_height_g)
print(s)

med_height = median(df_all$Height)
med_height_b = median(df_all[df_all$Gender==0,]$Height)
med_height_g = median(df_all[df_all$Gender==1,]$Height)

s <- sprintf("median of height=%.1f (boys=%.1f, girls=%.1f",med_height, med_height_b, med_height_g)
print(s)

df_all[条件式,]で条件式に合う要素だけを取り出します。条件式のあとのコンマ(,)を忘れがちなので注意です。（たまにやります笑）
結果は下の表のようになります。

代表値	全体	男子	女子
平均値	138.1	137.5	138.7
中央値	137.3	136.4	139.5

乱数で作ったデータですので、もとの統計の値とは若干ずれます。

分析１：学年～身長

まずは前回と同じ身長だけでのロジスティック回帰を行います。ソースコードは前回とほぼ変わらないので割愛します。出力の重要な部分は次のようになります。

[1] "border=137.578/OR=1.539"
[1] "correctness=0.865"

正解率86.5%です。200人ですから外れているのは27人ということになります。
この境界値は代表値と見比べてみると、全体の中央値に近いことが分かります。境界値は「6年生と判別する確率が50％の身長」でしたから、これが中央値付近になることも納得がいきます。

分析２：学年～身長＋性別

次に性別を加えた回帰を行います。もともとは連続データに対する回帰を行うのがロジスティック回帰ですが、0と1で表現されたカテゴリー変数（ダミー変数といいます）を入れても回帰をしてくれます。

result <- glm(Grade~Height+Gender, family=binomial(link="logit"), df_all)

回帰式を与える部分にGender変数を加えるだけです。まずは結果の係数(coefficients)そのものの値とそれを指数関数に入れたものを出します。

#print((coef(result)))の結果＝係数
(Intercept)      Height      Gender 
-65.1912136   0.4781025  -1.3927909
# print(exp(coef(result)))＝e^係数を計算したもの
 (Intercept)       Height       Gender 
4.873216e-29 1.613011e+00 2.483811e-01 

Heightは指数関数に入れてOR（オッズ比）にするので下をみますが、ほかの２つはどのように解釈すればよいでしょうか？Heightを$x_1$、Genderを$x_2$とすると、ロジスティック回帰の中の式は、

\begin{equation}
f(x_1,x_2) = a x_1 + b + c x_2
\end{equation}

という形になります。ここでの係数$a,b$は前回の表記に合わせているので、ちょっと変な順番にみえますが、$x_2$の項が付け加わったわけです。これを$e^c$とすれば下の段の数値から「女子であるということで学年のオッズが0.248倍になる」ということになりますが、これは解釈がわかりにくいと思います。ここではGender変数$x_2$が0と1だけしか取らないことを思い出し、女子の場合には回帰式が

\begin{equation}
f(x_1,x_2) = a x_1 + (b + c)
\end{equation}

のように定数項部分を変化させるという見かたができます。そして定数項は「境界」に関係したのですから、「女子であることで境界が$-\frac{c}{a}$だけずれる」という解釈をした方がわかりやすいのではないでしょうか？

この観点から「境界のズレ」を計算した結果が次のような出力になります。

[1] "border=136.354/OR=1.613, gender difference=2.913"
[1] "correctness=0.905"

つまり、男子の境界値は136.354cmであり、女子は2.913cm足して139.267cmを境界にする、ということです。今回のデータの中央値は、男子が136.4cm、女子が139.5cmでしたからなかなかイイ線をいっているのではないでしょうか。正解率も90.5%まで上昇しました。

分析３：交互作用を考慮する

さて、境界が男女で変わるならば、ORは変わらないのでしょうか？

ある変数の値（この場合は性別）が、別の変数の係数（この場合、身長によるOR）を変化させるとき、これを交互作用(interaction)といいます。交互作用を考慮したモデルの記述は以下のようになります。

result <- glm(Grade~Height+Gender+Height:Gender, family=binomial(link="logit"), df_all)

Height:Genderと：で結んだ部分が交互作用になります。これにはまとめて

result <- glm(Grade~Height*Gender, family=binomial(link="logit"), df_all)

という表記方法もあります。＊で結ぶとその他の主効果も含めて自動的に入れてくれる、ちょっとした省略表記になります。

実行結果は以下のようになります。

  (Intercept)        Height        Gender Height:Gender 
  -58.5449266     0.4291986   -16.8706193     0.1121687 
#指数関数
  (Intercept)        Height        Gender Height:Gender 
 3.751988e-26  1.536026e+00  4.711760e-08  1.118702e+00 
#信頼区間
(Intercept)   -89.8243431 -37.1219060
Height          0.2716414   0.6591436
Gender        -67.8776866  25.8010027
Height:Gender  -0.1980097   0.4792922
#指数関数にいれた信頼区間
Waiting for profiling to be done...
                     2.5 %       97.5 %
(Intercept)   9.767501e-40 7.553723e-17
Height        1.312116e+00 1.933136e+00
Gender        3.319673e-30 1.604106e+11
Height:Gender 8.203619e-01 1.614931e+00
[1] "border=136.405/OR=1.536, gender difference=39.307, interaction OR=1.119"
[1] "correctness=0.900/AIC=88.596"

性別による境界の変化の値が39.3とおかしなことになっています。それをみるために今回は信頼区間も入れてありますが、Genderの項の信頼区間が-68.9から25.8の非常に大きな幅を持ったものになっており、この値に意味があるようには思えません。これは多重共線性というもので、交互作用というはすでに入っている変数どおしの掛け算を新たな変数として付け加えるので、もとの変数と相関が強いものになってしまいます。こういう変数が入っていると解が安定しなくなります。

中心化による処方

詳しい理屈はここでは書きませんが、中心化という方法でこの問題を回避することができます。

これは簡単なことで身長のデータを平均値からの差にすることに相当します。

#中心化
df_all$Height = df_all$Height - m_height

result <- glm(Grade~Height+Gender+Height:Gender, family=binomial(link="logit"), df_all)

こうして身長の項目を平均身長からの差に変換して分析します。

a <- result$coefficients["Height"]
b <- result$coefficients["(Intercept)"]
c <- result$coefficients["Gender"]
i <- result$coefficients["Height:Gender"]

border <- -(b/a) + m_height
or <- exp(a)
d <- -c/a
b_g <- -(b+c)/(a+i) + m_height
i_or = exp(i)

ちょっといくつか説明を入れます。交互作用項というのは２つの変数の積です。上ではその係数をiとしていますので、次のような回帰式となっています。

\begin{equation}
f(x_1,x_2) = a x_1 + b + c x_2 + i x_1 x_2
\end{equation}

これの意味を見るために少し変形すると

\begin{equation}
f(x_1,x_2) = (a +i x_2) x_1 + (b + c x_2)
\end{equation}

とできます。つまり、女子の場合には回帰式が

\begin{equation}
f(x_1,x_2) = (a +i) x_1 + (b + c)
\end{equation}

となって身長に対する係数も$(a+i)$に変化するということを意味します。よって、女子の場合の境界は

\begin{equation}
border = -\frac{b+c}{a+i} + mean(Height)
\end{equation}

という形になっているわけです。

以上の結果より

  (Intercept)        Height        Gender Height:Gender 
    0.7300291     0.4291986    -1.3794301     0.1121687 
#指数関数
  (Intercept)        Height        Gender Height:Gender 
     2.075141      1.536026      0.251722      1.118702 
#信頼区間
Waiting for profiling to be done...
                    2.5 %     97.5 %
(Intercept)   -0.02342652  1.6192326
Height         0.27164135  0.6591436
Gender        -2.78209635 -0.1469858
Height:Gender -0.19800972  0.4792922
#指数関数に入れた信頼区間
Waiting for profiling to be done...
                   2.5 %    97.5 %
(Intercept)   0.97684575 5.0492143
Height        1.31211633 1.9331360
Gender        0.06190859 0.8633062
Height:Gender 0.82036188 1.6149310
[1] "border=(boys 136.405 girls 139.306)/OR=1.536,  interaction OR=1.119"
[1] "correctness=0.900/AIC=88.596"

となって、男子の境界が136.4cm、女子の境界が139.3cmとなりました。Genderの係数の信頼区間も幅はありますが0をまたがず、有意といってよい結果となりました。

どのモデルがよいのか？

ところで、最後の正解率のところですが、分析2が90.5%に対して、交互作用を入れたものでは90%とわずかに下がっています（一人分ハズレが増えた）。正解率だけがモデルの良さではないのですが、正解率のあとにAIC（赤池情報量規準）という量を入れてあります。これはモデルの良さを比較するための量で少ないほどよいことになります。ただし絶対的な数字に意味はなく、今回の４つのモデルを比較するために使います。

そこで以下の表に今回の結果をまとめてみます。

分析	1	2	3	4
境界	137.578	136.354	136.405	136.405
OR	1.539	1.613	1.536	1.536
性別	–	2.913	39.307	3.214
境界（女子）	–	139.267	175.712	139.306
交互作用	–	–	1.119	1.119
正解率	0.865	0.905	0.900	0.900
AIC	90.374	87.093	88.596	88.596

AICについてざっくり説明すると、モデルの当てはまりがいいと小さい値になりますが、変数（交互作用も含む）を増やすとペナルティがついて値が大きくなってしまいます。AICの観点からも分析２の交互作用を含まないモデルが最もよい、ということになります。このほか交互作用の検定などもありますが、これはまたどこかで取り上げられたらと思います。

モデル選択はいろいろな規準がありますので必ずこれがよいとも言えませんが、今回はあまり交互作用を考慮する必要はなさそうです。そもそも、今回の予測は境界と比べて大きいか小さいかのみでの判断ですからORは判別に関係していません。交互作用を入れることで境界が若干変わったことが結果の僅かな差になったわけですが、今回はあまり意味のある差ではないということが言えると思います。