Rでロジスティック回帰をしてみる

2023年11月9日2023年12月14日

ロジスティック回帰は説明変数が連続変数で、目的変数が2値（陽性・陰性、ある・なし、など）の場合の分類方法です。

今回のお題は次のようなものにします。

最近は政府系のデータはかなりアクセスがよくなっていまして、さまざまな統計のデータを簡単にとることができます。

今回は、令和3年度（2022年）学校保健統計調査から各学年の身長とその標準偏差のデータを持ってきます。この「統計の窓口」についてもどこかで詳しく触れたいですね。

データの生成

さて、学校保健統計調査から作った表の一部です。

この黄色いところを用いて次のようにデータをランダムに生成します。

#データの作成
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) #小学６年女子

これで各50人ずつのデータがランダムに生成されます。

データを視覚化する

次にヒストグラムを描いてみます。これは100cmから170cmを5cm刻みで作りました。

#ヒストグラムの作成
#個々のヒストグラムを表示する場合
h3b <- hist(e3b, breaks=seq(100,170,5))
h3g <- hist(e3g, breaks=seq(100,170,5))
h6b <- hist(e6b, breaks=seq(100,170,5))
h6g <- hist(e6g, breaks=seq(100,170,5))

これの結果を示します。

これをみるとやはり身長にそれなりの差があるように見えます。

これらを重ねてみます。

plot(0, 0, type="n", xlim=c(100,170), ylim=c(0, 30), xlab="", ylab="")
lines(h3b$mids, h3b$counts, col="blue")
lines(h3g$mids, h3g$counts, col="red")
lines(h6b$mids, h6b$counts, col="cyan")
lines(h6g$mids, h6g$counts, col="magenta")

3年生の男子が青、女子が赤、6年生の男子が水色、女子が紫色になっています。

こうしてみると、135cmの辺りに境目がありそうにみえます。

データをまとめる

次に、データを一つのデータフレームにまとめます。変数としてGrade(0=3年生、1=6年生)、Jender(0=男子,1=女子)を付け加えます。

#データの結合
df1 <- data.frame(Grade=0, Jender=0, Height=e3b)
df2 <- data.frame(Grade=0, Jender=1, Height=e3g)
df3 <- data.frame(Grade=1, Jender=0, Height=e6b)
df4 <- data.frame(Grade=1, Jender=1, Height=e6g)

df_all <- rbind(df1,df2,df3,df4)

これをRStudioで表示させると次のようなものができていると思います。Heightはランダムなので実際は違う値になっているはずです。

分析の実行

さて、ロジスティック回帰はGLM（一般化線形モデル）と呼ばれる分析方法の一つです。詳しいことは別記事で扱う予定ですが、オッズの対数を線形回帰することに相当します。

さて、実際の実行は次のようなコードになります。glm関数の中身はおまじないと思ってもらっていいですが、

Grade~Height　で、Gradeが目的変数、Heightが説明変数であることを指定しています。

#分析の実行
result <- glm(Grade~Height, family=binomial(link="logit"), df_all)
summary(result)

a <- result$coefficients[2]
b <- result$coefficients[1]

curve(1/(1+exp(-(a*x+b))),100,170)

今回の実行では、$a = 0.665, b = -90.8$の結果が得られました。出力されるグラフは下のようなものになります。これは身長$x$に対して、その子どもが6年生と判断される確率を表しています。

得られた係数を解釈します。そのために、ロジスティック回帰の式(\ref{logi})がちょうど確率$p=\frac{1}{2}$になるところを探します。そうすると、$ax+b=0$ならばよいので、$x_0 = -\frac{b}{a}$のところが確率$p=\frac{1}{2}$になります。

一方、aについては最初のオッズの式(\ref{odds})から、

\begin{equation}
\frac{1-p}{p} = e^{ax + b}
\end{equation}

ですから、身長$x$が1cm大きいとオッズは$e^a$倍になることが分かります。つまりこれはオッズ比になります。というわけで$a$が大きいほどオッズ比が大きい、すなわち分類が明確にできるということになります。

これらを計算しているところが次の部分です。

print(exp(coef(result)))
print(exp(confint(result)))

border <- -(b/a)
or <- exp(a)

borderが$x_0$でちょうど境目になっている身長、orがオッズ比(odds ratio)です。

今回のデータではborder=137、or=1.94ということになりました。境目が137cmというのはさきほどのヒストグラムを目で見た感覚とだいたい合っていると思います。

上の2つのprint文で出力される結果はそれぞれの推定値と信頼区間を表示するもので、次のようなものになっていると思います。

(Intercept)       Height 
3.538239e-40 1.944215e+00 

                   2.5 %       97.5 %
(Intercept) 2.547481e-58 4.350324e-28
Height      1.585396e+00 2.641038e+00

ここではHeightの方だけが重要で、Heightの推定値1.944、95％信頼区間(confidential Interval)が[1.585, 2.641]ということが分かります。オッズ比が1にならないので分類はきちんとできていそうだということが分かります。

答え合わせ

さて、ではこのロジスティック回帰でどのくらい正しく予測できるのでしょうか？

そこで回帰の結果$p>0.5$となったときに6年生、そうでないときに3年生と判断することにして、その結果を元のデータフレームに付け加えてみます。元のデータフレームをいじらないように新しいデータフレームに代入しています。

df_r <- transform(df_all, Predict=ifelse(result$fitted.values>0.5,1,0))
df_r <- transform(df_r, Correct=ifelse(Grade==Predict,1,0))

print(mean(df_r$Correct))

transformは変数を付け加える関数です。result＄fitted.valuesは各データにロジスティック回帰の結果を当てはめた結果の$p(x)$が計算されて入っています。それが0.5より大きいかどうかで1と0を分けています。今までの議論から分かる通り、これはboderより高ければ1、低ければ0と分けていることと変わりません。

2個めのtransformは、実際のGradeと予測されたPredictが同じならば正解の1、そうでなければ不正解の0を代入します。

この２つの変数を付け加えたデータフレームはこんな感じになります。

137cmぐらいを境目にしてPredictが0と1に分かれていることが分かります。

最後にCorrectの平均値を計算することで「正解率」がでます。今回の実行では0.935と表示されました。まあまあの結果ではないでしょうか。