どのぐらい受験勉強したら大学に合格するか， １日の平均勉強時間ごとの合格確率を知りたい。

データ

１日の平均勉強時間［hour］とある大学の合否 （1：合格，0：不合格） の仮想的データ

d <- read.csv(file = 'https://stats.dip.jp/01_ds/data/univ_exam_data.csv')

library(DT)
datatable(d, caption = '勉強時間と大学合否')

ロジスティック回帰モデル

\[ \begin{align} \hat{p}&=\frac{1}{1+\exp\{-(\beta_0+\beta_1 x_{i}+\cdots+\beta_k x_{k})\}}\\ &ここで，\\ &\hat{p}：確率の推定値\\ &x_{k}： 説明変数\\ &\beta_k： 偏回帰係数\\ \end{align} \]

ロジスティック回帰分析

フィッティング

説明変数で目的変数（0か1）に完全に分類できるときは次のような警告がでる。
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred このようなときは，分類木による分類がより分かり易いので，そちらを使用することを推奨。

fit <- glm(合格 ~ 勉強時間, data = d, family = 'binomial')
summary(fit)
## 
## Call:
## glm(formula = 合格 ~ 勉強時間, family = "binomial", data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -2.7611     1.2649  -2.183   0.0291 *
## 勉強時間      0.6635     0.2736   2.425   0.0153 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.726  on 19  degrees of freedom
## Residual deviance: 17.662  on 18  degrees of freedom
## AIC: 21.662
## 
## Number of Fisher Scoring iterations: 5
#library(sjPlot)
#tab_model(fit, show.aic = T)

グラフ

xp <- seq(0, 10, 0.1) # 勉強時間（初項0，公差0.1の等差数列）
yp <- predict(fit, type = 'response',
              newdata = data.frame(勉強時間 = xp))

matplot(x = d$勉強時間, y = d$合格, pch = 1,
        ylim = c(0, 1.1), xlim = c(0, 10),
        main = 'ロジスティック回帰分析',
        xlab = '平均勉強時間',
        ylab = '合格確率')
grid()

matlines(x = xp, y = yp, lty = 1, col = 2, lwd = 2)

library(latex2exp)
text(x = 3.8, y = 0.4, adj = 0, 
     labels = TeX('$\\leftarrow \\hat{p}=\\frac{1}{1+\\exp\\{-(b_0 + b_1 x)\\}}$'))
legend('topleft', pch = c(1, NA), lty = c(NA, 1), col = c(1, 2),
       legend = c('合格(1)，不合格(0)', 'ロジスティック回帰モデル'))

インタラクティブグラフ

library(plotly)
plot_ly(type = 'scatter') |>
  add_trace(x = d$勉強時間, y = d$合格,
            mode = 'markers', name = '合格(1)，不合格(0)') |>
  add_trace(x = xp, y = yp, 
            mode = 'lines',   name = 'ロジスティック回帰モデル') |> 
  layout(title = 'ロジスティック回帰分析', 
         xaxis = list(title = '平均勉強時間'),
         yaxis = list(title = '合格確率')) |>
  add_annotations(x = 3.6, y = 0.4, ax = 150, 
                  text = '$\\hat{p}=\\frac{1}{1+\\exp\\{-(b_0 + b_1 x)\\}}$') |> 
  config(mathjax = 'cdn')

予測

平均勉強時間が5.5時間の受験生の大学合格確率を予測せよ。

yp <- predict(fit, type = 'response', newdata = data.frame(勉強時間 = 5.5))
sprintf('合格確率：%2.1f％', yp * 100)
## [1] "合格確率：70.9％"

合否など二項に分類するときはp=0.5より 大きいか小さいかで分ける場合が多い。
この場合はyp > 0.5以上なので合格予想となる。

Python

## データ

import pandas as pd
import numpy as np
d = pd.read_csv('https://stats.dip.jp/01_ds/data/univ_exam_data.csv')
x = d['勉強時間'].to_numpy().reshape(-1, 1)
y = d['合格']
#y = pd.get_dummies(y) # 文字列変数の場合はダミー変数に変換する

## ロジスティック回帰分析

from sklearn.linear_model import LogisticRegression
model = LogisticRegression(random_state = 1)
model.fit(x, y)
LogisticRegression(random_state=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
## グラフ

xp = np.arange(0, 10, 0.1).reshape(-1,1)
yp = model.predict_proba(xp)

import matplotlib.pyplot as plt
#plt.switch_backend('agg')

# ラベル
plt.title('ロジスティック回帰分析')
plt.xlabel('平均勉強時間［Hour］')
plt.ylabel('合格（１），不合格（０）')

# 格子線（grid lines）
plt.grid(linestyle = '--', color = (0.9, 0.9, 0.9, 0.25))
plt.plot(x, y, 'o', label = 'Data points')
plt.plot(xp, yp[:, 1], linestyle = 'solid', 
         label = '$\hat{y}_i=\\frac{1}{1-\\exp(b_0 + b_1 x_i)}$')

# 凡例（はんれい）
plt.legend(loc = 'lower right')

plt.show()

xp1 = [[5.5]]
yp1 = model.predict_proba(xp1)
print('合格確率：', yp1[0, 1].round(1) * 100, '%')
## 合格確率： 70.0 %

インタラクティブグラフ

【注意】 add_traceのたびに描画されるのはrmarkdownのBugか？（2024-09-11現在）

xp = np.arange(0, 10, 0.1)

import plotly.graph_objects as go
fig = go.Figure()
fig.add_trace(go.Scatter(x=d['勉強時間'], y=d['合格'], 
                         mode = 'markers', name = '合格(1)，不合格(0)'))
fig.add_trace(go.Scatter(x=xp, y=yp[:, 1],
                         mode = 'lines', name = 'ロジスティック回帰モデル'))
fig.update_layout(title='ロジスティック回帰分析',
                   xaxis_title='平均勉強時間［Hour］',
                   yaxis_title='合格（１），不合格（０）')  
fig.show()

以下はPythonデータをRでdatatableを使って表示させている。 （Rチャンク）

library(reticulate)
library(DT)
datatable(py$d, caption = '勉強時間と大学合否')

演習課題

次のデータは米国男女の属性データである。
【単位】身長：cm，体重：kg，年収：US$

  1. 身長ごとに男性である確率を推定し図示せよ。
  2. 身長170cmの人の性別を推定せよ。

データ

d <- read.csv(file = 'https://stats.dip.jp/01_ds/data/gender.csv')

library(DT)
datatable(d, caption = '性別データ')

目的変数を2値（バイナリ）化： 女性：0，男性：1

d$性別 <- ifelse(d$性別 == '男性', 1, 0)