Rを使った分析の準備

ここで使うRのパッケージは次のとおり。

library(corrplot)
library(jtools)
library(margins)
library(ROCR)
library(patchwork)
library(prediction)
library(stargazer)
library(tidyverse)

1. 線形回帰と非線形回帰の違い

選挙での「候補者の成績」と「選挙で使うお金の額」の関係を調べたいとする
→　つまり、お金を使う候補者ほど選挙結果がいいのか？という疑問
説明したい対象（＝応答変数）としては 2 つ想定できる：
①得票率（もしくは得票数）　　　　　・・・連続変数 => 回帰分析
②当落・・・当選したら1 、落選したら 0 ・・・2値変数 => ロジスティック回帰分析
回帰分析とロジスティック回帰分析の違いを視覚的に示してみる

	回帰分析	ロジスティック回帰分析
応答変数	得票率(%)	当落(0 or 1)
説明変数	選挙費用（百万円）	選挙費用（百万円）

主要な違い

・回帰分析の応答変数は得票率＝「連続変数」
・ロジスティック回帰分析の応答変数は選挙での当落＝ 2 値変数 (1 or 0)
・それぞれの図に回帰直線を図に書き入れてみる

・回帰分析では、選挙費用を使うほど得票率が大きくなっている => 問題なし
・ロジスティック回帰分析でも、選挙費用を使うほど得票率が大きくなっている
→　しかし、3つ問題点がある
1. 当選確率（0〜1）が 1 (100%) を超えてしまうこと
2. 当選確率（0〜1）がマイナスの値をとってしまうこと
3. 当てはまり感がない（むりやり当てはめている感が強い）

2. ロジスティック回帰分析とは

「ロジスティック回帰分析」とは、2値変数をロジット変換した回帰分析

ロジット変換とは

2値を表す散布図を「率」に変換すると S 字になる（左の図 → 中央の図）
率をロジット変換すると直線になる（中央の図 → 右の図）
ロジスティック回帰分析ではこの数学の特徴を使う

選挙での「当落」（左側の図）の代わりに「当選率」を考えてみる（中央の図）
選挙費用ごとにどれくらいの割合が当選しているかを棒グラフで表してみる
中央の「当選率」の棒グラフを見ると・・・

・選挙費用が少ない200万円以下の人はほとんど当選していない
・選挙費用が1000万円の候補者は50%程当選している
・選挙費用が1600万円の候補者は75%以上が当選している
・選挙費用ごとの「当選率」の対数をとってロジット変換する
→　Ｓ字が直線になる
→　しかし、実際のデータを使うと完璧な直線ではないものの、かなり直線ぽくなる
→　この曲線を使ってロジスティック回帰分析を推定する

2値変数をロジットに変換してから回帰分析を行う→ロジスティック回帰分析と呼ばれる

・ロジスティック回帰分析では「割合」ではなく「オッズ」で結果が示される

Logit Model の基礎知識：オッズ

オッズ (`Odds`) とは

例えば、選挙費用を800万円使う候補者は100人いる
そのうち当選したのは60人で落選したのは40人
この場合の当選者の「割合」と「オッズ」は次のように計算できる

「割合」は 1 を超えることはない
しかし、「オッズ」は 1 を超えることはある
「割合」も「オッズ」もどちらも起こりやすさを示す指数

オッズのまとめ

ある事象が起こる確率 \(\mathrm{p}\) と起こらない確率 \(\mathrm{1 - p}\) との比

p: ある事象が起こる確率
1-p: ある事象が起こらない確率

\[Odds = \frac{p}{1-p}\]

例えば、p = 0.01（起こる確率が1%）の場合の Odds を計算してみる
p = 0.01, 1-p = 1-0.01= 0.99 を Odds を求める式に代入する

\[Odds = \frac{p}{1-p} = \frac{0.01}{0.99} = \frac{1}{99} = 0.01\]

Odds の特徴・オッズが大きいほど、その事象 (p) が起こりやすい
・オッズの最小値は 0
・オッズ = 1 は「その事象 (p) が起こる確率が 50%」という意味
・オッズの最大値は無限大 (∞)

Logit Model の基礎知識：ロジット変換

ロジット変換とは

・ロジット変換とは 2 値変数 (0 or 1) のオッズの対数をとること
・Odds の下限は 0 なので、説明変数としては扱いにくい
→ Odds [= p / (1 - p)] を対数変換し「対数オッズ」(Log Odds)を計算する
→　これが logit

\[LogOdds =log\frac{p}{1-p} = logit\]

ロジット (logit) は対数オッズ (Log Odds) として解釈できる

Log Odds (= logit) の特徴・最小値も最大値もどちらも無限大 (∞)
・Log Odds = 0 が意味するのは「その事象 (p) が起こる確率が 50%」という意味

3. ロジスティック回帰分析の手順

ロジスティック回帰分析の手順は次のとおりである

対立仮説と帰無仮説を設定する
説明変数と応答変数の散布図を表示する
ロジスティック回帰式を求める
回帰係数の有意性を検定する
推定結果の意味を解釈する
ロジスティック回帰モデルの当てはまり具合を評価する

3.1 対立仮説と帰無仮説を設定する

ここで検証したい仮説（＝対立仮説）は次のとおり

仮説1・選挙費を使えば使うほど、小選挙区での当選確率は大きい

従って、否定されるために設定する帰無仮説は次のようになる

帰無仮説・選挙費の額は、小選挙区での当選確率とは関係がない

ロジスティック回帰分析における検定では、重回帰分析における検定と同様、得られた p 値 が有意水準よりも小さいときに帰無仮説を棄却し、対立仮説を受け容れる
仮説1を検証するために次のモデルを考える

Model 1

ダウンロードした 2021年衆院選挙データを使い、候補者の小選挙区における当落 (wlsmd) を縦軸に、選挙費用 (expm) を横軸にした散布図を表示する

データの準備 (hr96-24.csv)

ggplot2のテーマとフォントの設定を行う

if (.Platform$OS.type == "windows") { 
  if (require(fontregisterer)) {
    my_font <- "Yu Gothic"
  } else {
    my_font <- "Japan1"
  }
} else if (capabilities("aqua")) {
  my_font <- "HiraginoSans-W3"
} else {
  my_font <- "IPAexGothic"
}

theme_set(theme_gray(base_size = 9,
                     base_family = my_font))

衆議院議員総選挙の得票データhr96-24.csv をダウンロード

データのダウンロードが終わったら、データを読み込み hr と名前を付ける。

hr <- read_csv("data/hr96-24.csv", na = ".")

データフレーム hr の中身を表示する

names(hr)

 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"

select() 関数を使って year, wl, previous, exp という 4 つの変数だけを選ぶ
filter() 関数を使って 2021年衆院選だけのデータを残す
wl の中身を確認

unique(hr$wl)

[1] 1 0 2

0・・・小選挙区での落選
1・・・小選挙区での当選
2・・・復活当選
wl を使って、選挙結果を示す変数 wlsmd （1 = 小選挙区当選、0 = その他）を作る
exp を使って、選挙費用を示す変数 expm（単位は「百万円」）を作る

hr21 <- hr %>%
  select(year, wl, previous, exp) %>%
  filter(year == 2021) %>% 
  mutate(wlsmd = ifelse(wl == 1, 1, 0)) |> 
  mutate(expm = exp / 1000000) |> 
  
  select(year, wlsmd, expm, previous)

データフレーム hr21 の中身を表示する

hr21

# A tibble: 857 × 4
    year wlsmd  expm previous
   <dbl> <dbl> <dbl>    <dbl>
 1  2021     1 13.4         3
 2  2021     0  9.62        2
 3  2021     0 NA           0
 4  2021     1  8.82        8
 5  2021     0 11.4         0
 6  2021     1  7.42        8
 7  2021     0 11.9         3
 8  2021     1 13.0         3
 9  2021     0  7.48        6
10  2021     0  4.61        0
# ℹ 847 more rows

データ：

変数名	詳細
`year`	衆院選が行われた年
`wlsmd`	小選挙区での当落ダミー(当選 = 1, 落選 = 0)
`previous`	当選回数
`expm`	候補者の選挙費用（百万円）

データフレーム hr21 の記述統計を表示させる

summary(hr21)

      year          wlsmd             expm              previous     
 Min.   :2021   Min.   :0.0000   Min.   : 0.009319   Min.   : 0.000  
 1st Qu.:2021   1st Qu.:0.0000   1st Qu.: 3.435080   1st Qu.: 0.000  
 Median :2021   Median :0.0000   Median : 5.899882   Median : 1.000  
 Mean   :2021   Mean   :0.3372   Mean   : 6.434585   Mean   : 2.162  
 3rd Qu.:2021   3rd Qu.:1.0000   3rd Qu.: 8.692808   3rd Qu.: 3.000  
 Max.   :2021   Max.   :1.0000   Max.   :27.443685   Max.   :17.000  
                                 NA's   :22

データフレーム hr21 の expm に欠損値 (missing data = NA's) が22個あることがわかる
na.omit() を使って欠測のない観測だけを残す

hr21 <- na.omit(hr21)

データフレーム hr21 の記述統計を表示して確認する

summary(hr21)

      year          wlsmd             expm              previous     
 Min.   :2021   Min.   :0.0000   Min.   : 0.009319   Min.   : 0.000  
 1st Qu.:2021   1st Qu.:0.0000   1st Qu.: 3.435080   1st Qu.: 0.000  
 Median :2021   Median :0.0000   Median : 5.899882   Median : 1.000  
 Mean   :2021   Mean   :0.3461   Mean   : 6.434585   Mean   : 2.211  
 3rd Qu.:2021   3rd Qu.:1.0000   3rd Qu.: 8.692808   3rd Qu.: 3.000  
 Max.   :2021   Max.   :1.0000   Max.   :27.443685   Max.   :17.000

22個の NA's が消えていることがわかる

3.2 説明変数と応答変数の散布図を表示する

分析で使うデータ (hr21) の記述統計を表示

stargazer::stargazer(as.data.frame(hr21), 
  type = "html")


Statistic	N	Mean	St. Dev.	Min	Max

year	835	2,021.000	0.000	2,021	2,021
wlsmd	835	0.346	0.476	0	1
expm	835	6.435	4.287	0.009	27.444
previous	835	2.211	2.921	0	17

当選回数と当落の散布図を描く.
2つの変数の関係を図示してみる
文字化けを避けるため、Macユーザは次の2行を実行する

theme_set(theme_gray(base_size = 10, 
                     base_family = "HiraginoSans-W3"))

wlsmd が 2 値変数 (0 or 1) なので、jitter()関数を使ってデータを散らして表示させる

p <- ggplot(hr21, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,        # データを散らして表示させる指定
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当落")

plot(p)

この図にむりやり通常の回帰直線を当てはめてみるとこのようになる

p + geom_smooth(method = "lm", se = FALSE) +
    annotate("label", 
           label = "当落 (0 or 1) = a + b選挙費用", 
           x = 6, y = 0.75,
           size = 5, 
           colour = "blue", 
           family = "HiraginoSans-W3")

応答変数が2値 (o or 1) なのに、むりやり回帰直線を当てはめた時の問題点
→ 応答変数（当選確率）がマイナスの値になったり、100%を超える場合がある

解決策

当選確率を「ロジット変換」する
当選確率 \(P\) は 0 と 1 の間の値をとる 2 値変数
\(\frac{P}{1-P}\) はオッズ（コラムLogit Model の基礎知識：オッズを参照）
オッズ＝当選する確率 \(\mathrm{p}\) と落選する確率 \(\mathrm{1 - p}\) との比
「当選確率\(P\)」と「選挙費用」の関係は次の式で表現できる

\[\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]

回帰式の左辺のオッズ \(\frac{p}{1-P}\)は曲線を表現している
回帰式の右辺 \(a + 選挙費用b\) は直線を表現している
→　表現が一致していない
回帰式の左辺をする＝オッズ \(\frac{p}{1-P}\)の対数をとる
当選確率 \(P\) は 0 から 1 の範囲しか取れないが、オッズの対数をとると全ての実数をとれる
→　曲線の関係を直線の関係に変更でき、推定が可能になる＝「ロジット変換」

\[log\frac{P}{1-p} = a + 選挙費用b + 当選回数c\]

オッズの対数をとった曲線

p3 <- ggplot(hr21, aes(x = expm, y = wlsmd)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  geom_smooth(method = "glm", 
    color = "red",
    method.args = list(family = binomial(link = "logit"))) +
  labs(x = "選挙費用（100万円）",
       y = "小選挙区での当選確率")
print(p3)

3.3 ロジスティック回帰式を求める

ロジスティック回帰式 (`Model 1`)

ここでは2021年の衆院選データを用い、選挙費用（expm: 単位100万円）で衆院選小選挙区での当落 (wlsmd) を説明するモデルを考える

Model 1

このモデルを式で表すとこのようになる

\[Pr(当選)=logit^{−1}(𝛼+\beta_1選挙費用 + \beta_2当選回数)\]

Rでロジスティック回帰分析を行うには、一般化線形モデル (Generalized Linear Models: GLM) を当てはめるための関数 glm() を使う
glm()はロジスティック回帰以外でも頻繁に使う関数
この関数をロジスティック回帰分析に使うには、引数 family = binomial(link = "logit") を指定
logit とは「オッズの対数」= Log Odds のこと

model_1 <- glm(wlsmd ~ expm + previous, 
            data = hr21, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定

tidy() を使って、推定結果を確認する

tidy(model_1, conf.int = TRUE)

# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   -2.65     0.194     -13.7  9.90e-43   -3.04     -2.28 
2 expm           0.170    0.0235      7.23 4.86e-13    0.125     0.217
3 previous       0.342    0.0355      9.62 6.38e-22    0.274     0.413

→　表示される回帰式の係数 estimate は「オッズの対数」 = Log Odds

Log Odds は解釈しにくい → 確率に変換する
当選確率を予測するためのロジスティック回帰式は次のとおり

\[Pr(当選) = logit^{-1}(\alpha + \beta_1選挙費用 + \beta_2当選回数)\]

\[= \frac{1}{1+exp(-[\alpha + \beta選挙費用])}\] \[= \frac{1}{1+exp(-[-2.65 + 0.170[expm] + 0.342[previous]])}\]

3.4 回帰係数の有意性を検定する

tidy() を使って、推定結果を再度確認してみる

tidy(model_1, conf.int = TRUE)

# A tibble: 3 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   -2.65     0.194     -13.7  9.90e-43   -3.04     -2.28 
2 expm           0.170    0.0235      7.23 4.86e-13    0.125     0.217
3 previous       0.342    0.0355      9.62 6.38e-22    0.274     0.413

expm の p.value は 1.46e-24
→　有意水準を 0.05 とすれば、選挙費用 (expm) の p値は 1.61e- 5 なので、「影響がない」という帰無仮説を棄却する
→ 選挙費用が当選確率に影響すると判断する
previous の p.value は 4.29e-14
→　有意水準を 0.05 とすれば、過去の当選回数 (previous) の p値は 4.29e-14 なので、「影響がない」という帰無仮説を棄却する
→ 過去の当選回数が当選確率に影響すると判断する

ロジスティック回帰分析の係数

→　表示される回帰式の係数 estimate は「オッズの対数」 = Log Odds

Log Odds は解釈しにくい → 確率に変換する
当選確率を予測するためのロジスティック回帰式は次のとおり

\[Pr(当選) = logit^{-1}(\alpha + \beta_1選挙費用 = \beta_2当選回数)\]

\[= \frac{1}{1+exp(-(\alpha + \beta選挙費用)}\] \[= \frac{1}{1+exp(-(-2.65 + 0.170expm + 0.342previous))}\]

例えば、選挙費用 (expm) の係数の推定値は約 0.170
これが重回帰分析の結果であれば、「選挙費用が 1 単位（100 万円）増えるごとに、当選確率が 0.170 ポイント大きくなる」という意味
しかし、ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない

予測当選確率と限界効果

予測当選確率

選挙費用の値ごとに候補者が当選する確率の予測値が当選予測確率: 0 〜 1

限界効果

ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない
ロジスティック曲線は直線ではない → 「曲線」
→ 説明変数 1単位の変化が応答変数に与える影響（＝傾き）は説明変数の値によって異なる
説明変数が応答変数に与える影響を考えるためには
→　説明変数の値を特定し、「特定された値における影響」を計算する必要あり
→　限界効果 (Marginal Effects)
上図の例では、xの値 1, 5, 10, 15, 20 それぞれの時点での直線の傾きが限界効果

予測当選確率

当選回数 (previous) が3回の候補者に関して次の二つのケースを考えてみる
例1: 選挙費用を0 円から 100 万円に変化させた場合の予測当選確率
例2: 選挙費用を100 万円から 200 万円に変化させた場合の予測当選確率

例1: 選挙費用 (0 円 → 100 万円)

ここでは当選回数が 3 回の候補者が、選挙費用を0 円から 100 万円に変化させると、応答変数にどのような変化があるか確かめる
→それぞれの予測当選確率を計算して差を取る
当選回数が 3 回で選挙費用が 0 円の候補者の予測当選確率は、次の回帰式に、previous = 3、expm = 0 を代入すれば求められる

\[Pr（当選） = \frac{1}{1+exp(-(-1.98 + 0.0735expm + 0.285previous))}\]

この値を R で求めてみる
当選回数が 3 回で選挙費用が 0 円の候補者の予測当選確率は

p_0 <- predict(model_1, type = "response", #予測当選確率を表示したい時の指定  
               newdata = data_frame(previous = 3, expm = 0))
p_0

        1 
0.1642294

0.1642294 (= 約16.4%) である
当選回数が 3 回で選挙費用が 100 万円の候補者の予測当選確率は

p_1 <- predict(model_1, type = "response",
               newdata = data_frame(previous = 3, expm = 1))
p_1

        1 
0.1889165

0.1889165 (= 約18.9%)だとわかる
当選回数が 3 回の候補者が選挙費用を 0 円から 100 万円に変化させたときの予測当選確率の変化は、上で求めた二つの確率の差 (p_1 - p_0)

p_1 - p_0

         1 
0.02468706

当選回数が 3 回の候補者が選挙費用を 0 円から 100 万円に増やすと、予測当選確率は約 2.47 %ポイント増える

例2: 選挙費用 (100万円 → 200万円)

ここでは当選回数が 3 回の候補者が、選挙費用を100万円から 200万円に変化させると、応答変数にどのような変化があるか確かめる
当選回数が 3 回で選挙費用が 200万円の候補者の予測当選確率は 0.2163539（約21.6%ポイント）

p_2 <- predict(model_1, type = "response",
               newdata = data_frame(previous = 3, expm = 2))
p_2

        1 
0.2163539

当選回数が 3 回の候補者が選挙費用を 100 万円から 200 万円に変化させたときの予測当選確率の変化は、上で求めた二つの確率の差 (p_2 - p_1) なので

p_2 - p_1

         1 
0.02743738

当選回数が 3 回の候補者が選挙費用を 100 万円から 200万円に増やすと、予測当選確率は 0.02743738 (= 約2.7% ポイント) 増える

まとめ・選挙費用を 0 円から 100 万円に増やすと、予測当選確率は 2.47%ポイント上昇する
・選挙費用を 100 万円から 200 万円に増やすと、予測当選確率は 2.7%ポイント上昇する
・選挙費用の0 円から 100 万円への変化も、100 万円から 200 万円への変化も、選挙費用の増分としては同じ
・しかし、それらの変化が応答変数（すなわち当選確率）に与える影響は異なる
・0 円から 100 万円への 100 万円の変化より、100 万円から 200 万円への 100 万円の変化のほうが、応答変数に与える影響が大きい
・実際に何%ポイント予測当選確率が上昇するかを、具体的な予測当選確率を計算せずに係数の推定値(expm = 0.170) から読み取るのは難しい


	Dependent variable:

	wlsmd

expm	0.170^***
	(0.024)

previous	0.342^***
	(0.036)

Constant	-2.652^***
	(0.194)


Observations	835
Log Likelihood	-397.969
Akaike Inf. Crit.	801.939

Note:	p<0.1; p<0.05; p<0.01

・この推定値は「オッズの対数」= Log Odds = logit
→ いくつか特定した expm の値で予測確率を計算し、その差を比較する作業が必要

限界効果

`margins` を使った限界効果の計算

上の例のように、予測当選確率を複数計算し、その差をとるという作業を繰り返すのは面倒
→ 特定の値における限界効果を直接計算するための marginsパッケージを使うのが便利
margins()関数を使うと、特定の選挙費用額ごとの限界効果を求めることができる

選挙費用の限界効果 (previous = 3, expm = 400万円)

例えば、当選回数が 3 回で選挙費用が 400 万円のとき、選挙費用 1 単位（= 100 万円）の増加が当選確率に与える影響の大きさ（＝限界効果）は

margins(model_1, variables = "expm",
        at = list(previous = 3, expm = 4))

 at(previous) at(expm)    expm
            3        4 0.03424

選挙費用が当選確率に与える影響（限界効果）は 0.03424 (約3.4%ポイント) だとわかる

選挙費用の限界効果（previous = 3, expm: 0 〜 2,800万円)

どの値での限界効果を求めるかは、引数 atに変数のリストを渡すことで指定する
複数の値、例えば、選挙費用が 0円 (previous = 0) から 2,800万円 (previous = 3) までの限界効果を同時に表示してみる

margins(model_1, variables = "expm",
at = list(previous = 3, expm = c(0:28))) #選挙費用を0円から2,800万円まで指定

 at(previous) at(expm)     expm
            3        0 0.023337
            3        1 0.026052
            3        2 0.028827
            3        3 0.031585
            3        4 0.034239
            3        5 0.036685
            3        6 0.038820
            3        7 0.040540
            3        8 0.041758
            3        9 0.042406
            3       10 0.042449
            3       11 0.041884
            3       12 0.040742
            3       13 0.039087
            3       14 0.037004
            3       15 0.034595
            3       16 0.031964
            3       17 0.029215
            3       18 0.026437
            3       19 0.023710
            3       20 0.021093
            3       21 0.018631
            3       22 0.016352
            3       23 0.014273
            3       24 0.012397
            3       25 0.010723
            3       26 0.009241
            3       27 0.007938
            3       28 0.006801

ロジスティック回帰分析の係数の限界効果は、説明変数の値によって変化する
→ 各説明変数の係数を見ただけで結果を理解するのは困難
→ 説明変数が応答変数に影響を与える様子をいくつかの条件について図示する
→ 説明変数の値に応じた限界効果 (ME:Marginal Effects) を図示する必要がある

3.5 推定結果の意味を解釈する

仮説 1 の結果の表示方法 (`Model_1`)

「予測当選確率」と「限界効果」は異なる

「予測当選確率」は特定の選挙費用を使った候補者が当選する確率
「限界効果」は選挙費用が予測当選確率に与える影響（＝傾き）
予測当選確率と限界効果は margins::cplot() を利用して図示する
・Y軸に予測当選確率を表示したい時は what = "prediction" と指定
・Y軸に限界効果を表示したい時は what = "effect" と指定

「選挙費用」と「予測当選確率」(`hr21_me`)

hr21_me <- cplot(model_1, 
                  x = "expm", 
                  what = "prediction") %>%
  as_data_frame() %>%
  ggplot(aes(x = xvals, y = yvals, ymin = lower, ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
labs(x = "選挙費用（単位：100万円）", 
     y = "予測当選確率の予測値", 
     title = "予測当選確率の予測値:2021総選挙")

「選挙費用」と「選挙費用の限界効果」」(`hr21_pred`)

hr21_pred <- cplot(model_1,  
                  x = "expm",    # x軸に据える変数
                  dx = "expm",   # 説明変数
                  what = "effect") %>%
  as_data_frame() %>%
  ggplot(aes(x = xvals, y = yvals, ymin = lower, ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
  geom_abline(intercept = 0, slope = 0, linetype = "dashed", color = "red") + 
  ylim(-0.0001, 0.06) +
labs(x = "選挙費用", 
     y = "選挙費用の限界効果", 
     title = "選挙費用の限界効果:2021総選挙")

仮説検証のまとめ（仮説 1）:

hr21_me + hr21_pred

仮説検証の結論（仮説 1）

左側の図からわかること

実質的有意性:
・選挙費用が大きくなるにつれて、候補者の当選確率は大きくなる
・選挙費用が小さいとき、予測当選確率は緩やかに上昇
・選挙費用が中程度になると曲線の傾きが少しずつ急になる
・選挙費用が大きくなると再び傾きが小さくなる

右側の図からわかること

実質的有意性:
・選挙費用の限界効果は、約1000万円の時が最大
・選挙費用がそれより多くても少なくても効果が小さい
→　選挙費用が約1300万円の時、当選確率に対する影響力が最大
統計的有意性: 全ての選挙費用の範囲で統計的に有意

選挙費用の限界効果の最小値と最大値を確認したければ、次のように入力する
この限界効果は、モデルに含まれているコントロール変数 (previous) の値を平均値に固定した値である

margins(model_1, variables = "expm",
  at = list(expm = c(0:28))) #選挙費用を0円から2,800万円まで指定

 at(expm)    expm
        0 0.01911
        1 0.02082
        2 0.02257
        3 0.02435
        4 0.02614
        5 0.02789
        6 0.02959
        7 0.03119
        8 0.03264
        9 0.03391
       10 0.03494
       11 0.03569
       12 0.03611
       13 0.03618
       14 0.03587
       15 0.03518
       16 0.03413
       17 0.03275
       18 0.03107
       19 0.02916
       20 0.02708
       21 0.02489
       22 0.02266
       23 0.02045
       24 0.01831
       25 0.01626
       26 0.01435
       27 0.01259
       28 0.01099

expm = 13 の時に限界効果が最大値 (0.03618) であることが確認できる

当選回数ごとの予測当選確率

予測当選確率

選挙費用の額に関わりなく、選挙費用が予測当選確率に正の影響を与えていることはわかった
それは、候補者の当選回数が異なってもいえるのか？
候補者の当選回数ごとに選挙費用が予測当選確率に与える影響をチェックする
候補者の当選回数 (previous) の記述統計を確認する

table(hr21$previous)


  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  16  17 
390  79  55 106  51  31  33  30  24  16   7   3   5   2   1   1   1

hist(hr21$previous)

候補者の大多数の当選回数が 0 だとわかる
17回の当選回数を0から2回ずつ区切り、0, 2, 4, …, 16まで分割する
→　previous = seq(0, 16, by = 2)
候補者の選挙費用 (expm) の記述統計を確認する

summary(hr21$expm)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
 0.009319  3.435080  5.899882  6.434585  8.692808 27.443685

hist(hr21$expm)

- 選挙費用は0円から2700万円強まで分布
→ 500万円ずつ分割する →　expm = seq(0, 20, by = 5))
- 当選回数ごとの予測当選確率を可視化してみる

df_pre <- expand.grid(previous = seq(0, 17, by = 2), 
  expm = seq(0, 20, by = 5)) %>% 
  as_data_frame() 
pred <- predict(model_1, 
  type = "response", 
  newdata = df_pre, 
  se.fit = TRUE) 

df_pre$fit <- pred$fit 
df_pre$lower <- with(pred, fit - 2 * se.fit) 
df_pre$upper <- with(pred, fit + 2 * se.fit) 
df_pre <- df_pre %>%

mutate(lower = ifelse(lower < 0, 0, lower), 
  upper = ifelse(upper > 1, 1, upper)) 
plt_prob <- ggplot(df_pre, aes(x = expm, y = fit)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = "gray") +
geom_line() +

facet_wrap(. ~ previous) + 
  labs(x = "選挙費用（100万円）", y = "当選確率の予測値") 
print(plt_prob)

当選回数が少ない候補者は、選挙費用が当選確率の予測値に強い正の影響を与えている
当選回数が増えるにつれて、選挙費用が当選確率の予測値に与える影響が小さくなる
→　当選回数16回の候補者は、選挙費用を使わなくても使ってもほぼ当選している

限界効果

次に、当選回数が変わっても、選挙費用が当選確率に予測値に与える影響が統計的に有意なのかどうかをチェック
当選回数 (previous) が増えるにつれて、選挙費用 (expm) の影響がどう変化するかを指定
→ x = "previous", dx = "expm"

plt_previous <- cplot(model_1, 
                  x = "previous", # x軸に据える変数
                  dx = "expm", 　 # 説明変数
                  what = "effect") %>%
  as_data_frame() %>%
  ggplot(aes(x = xvals, y = yvals, ymin = lower, ymax = upper)) +
  geom_ribbon(fill = "gray") +
  geom_line() +
  geom_abline(intercept = 0, slope = 0, linetype = "dashed", color = "red") + 
labs(x = "当選回数", 
     y = "選挙費用の限界効果", 
     title = "当選回数ごとの選挙費用の限界効果:2021総選挙")

print(plt_previous)

previous = 15までは、95%信頼区間が y = 0 の赤い点線に触れていない
→ 当選回数が変わっても、previous = 15 までは、選挙費用が当選確率に予測値に与える影響が統計的に有意

当選回数ごとの予測当選確率：まとめ

plt_prob + plt_previous

33. ロジスティック回帰分析

Masahiko Asano

2025-03-10

1. 線形回帰と非線形回帰の違い

主要な違い

2. ロジスティック回帰分析とは

ロジット変換とは

オッズ (`Odds`) とは

オッズのまとめ

ロジット変換とは

3. ロジスティック回帰分析の手順

3.1 対立仮説と帰無仮説を設定する

Model 1

3.2 説明変数と応答変数の散布図を表示する

解決策

オッズの対数をとった曲線

3.3 ロジスティック回帰式を求める

ロジスティック回帰式 (`Model 1`)

3.4 回帰係数の有意性を検定する

ロジスティック回帰分析の係数

予測当選確率

限界効果

予測当選確率

例1: 選挙費用 (0 円 → 100 万円)

例2: 選挙費用 (100万円 → 200万円)

限界効果

`margins` を使った限界効果の計算

3.5 推定結果の意味を解釈する

仮説 1 の結果の表示方法 (`Model_1`)

「予測当選確率」と「限界効果」は異なる

「選挙費用」と「予測当選確率」(`hr21_me`)

「選挙費用」と「選挙費用の限界効果」」(`hr21_pred`)

左側の図からわかること

右側の図からわかること

当選回数ごとの予測当選確率

予測当選確率

限界効果

33. ロジスティック回帰分析

Masahiko Asano

2025-03-10

1. 線形回帰と非線形回帰の違い

主要な違い

2. ロジスティック回帰分析とは

ロジット変換とは

オッズ (Odds) とは

オッズのまとめ

ロジット変換とは

3. ロジスティック回帰分析の手順

3.1 対立仮説と帰無仮説を設定する

Model 1

3.2 説明変数と応答変数の散布図を表示する

解決策

オッズの対数をとった曲線

3.3 ロジスティック回帰式を求める

ロジスティック回帰式 (Model 1)

3.4 回帰係数の有意性を検定する

ロジスティック回帰分析の係数

予測当選確率

限界効果

予測当選確率

例1: 選挙費用 (0 円 → 100 万円)

例2: 選挙費用 (100万円 → 200万円)

限界効果

margins を使った限界効果の計算

3.5 推定結果の意味を解釈する

仮説 1 の結果の表示方法 (Model_1)

「予測当選確率」と「限界効果」は異なる

「選挙費用」と「予測当選確率」(hr21_me)

「選挙費用」と「選挙費用の限界効果」」(hr21_pred)

左側の図からわかること

右側の図からわかること

当選回数ごとの予測当選確率

予測当選確率

限界効果

オッズ (`Odds`) とは

ロジスティック回帰式 (`Model 1`)

`margins` を使った限界効果の計算

仮説 1 の結果の表示方法 (`Model_1`)

「選挙費用」と「予測当選確率」(`hr21_me`)

「選挙費用」と「選挙費用の限界効果」」(`hr21_pred`)