分析に必要なパッケージをロードする

library(estimatr)     # for lm_robust()  
library(modelsummary) # for modelsummary()  
library(tidyverse)
library(stargazer)

1. ダミー変数でわかること

1.1 ダミー変数とは

ある属性があるかどうかを示す 2 値変数
特定のカテゴリーに属している場合・・・1
そのカテゴリーに属さない場合・・・0
例）
性別（女性 = 1、男性 = 0）
選挙結果（当選 = 1、落選 = 0）
社会の状態（戦時 = 1、平和時 = 0）
州の位置（北部 = 1、南部 = 0）などなど・・・

1.2 ダミー変数でわかること

「9. 回帰分析 1（単回帰と重回帰）」のロバート・パットナムの事例
経済状況（横軸）と政府のパフォーマンスの間には正の相関・・・左側の図
「イタリア州政府の位置（南北）」という northダミーをモデルに入れる
→　経済状況（横軸）と政府のパフォーマンスの間には負の相関（無相関）・・・右側の図

2. 南北格差ダミーと経済の近代化

2.1 理論と仮説

リサーチクエスチョン：
「イタリア地方政府のパフォーマンスに著しい違いがあるのはなぜか？」

理論：
Social Capital（社会関係資本）が政府のパフォーマンスを高める

地方政府のパフォーマンスの違いは、その地域の社会関係資本の蓄積の度合いによって説明できる
「社会関係資本」・・・個人の結びつき → 互恵性の社会ネットワークや規範
（＝見知らぬ相手と協力関係を構築する一助となるもの）
社会関係資本の蓄積の高い地域
→　互いに信頼し協力しあう
→　政府のパフォーマンスを高める

【応答変数】

「政府のパフォーマンス」を作業化 → gov_p : 12の指標から構成

地方政府の内閣の安定性
予算通過の早さ
統計- 情報サービスの提供

【説明変数】

「社会関係資本の蓄積の度合」を作業化 → cc: Civic Community Index（市民共同体指標）

比例代表での個人名記入投票の割合 = Clientelism（政治的恩顧主義）の度合
住民投票での投票率 = 地域社会への関心の度合
新聞購読者の割合 = 市民的な熟慮能力の度合
スポーツ- 文化団体の割合 = 市民の社交的生活の度合

仮説
もしこの理論が正しいなら、
cc（市民共同体指標）が大きくなるほど gov_p（政府のパフォーマンス）も大きくなるはず

応答変数：gov_p（政府のパフォーマンス）
説明変数：cc（Civic Community Index 市民共同体指標）
コントロール変数：econ（地方政府の経済指標）

2.2 oldberg による Putnam 批判とその検証

【イタリアにおける南北格差を考慮すべき】　　

oldberg (1996) による Putnam (1994)への批判
イタリアにおいて北部と南部とは全く異なる歴史、伝統、文化をもつ
政治、経済、社会の状態の違いはすべて南北地域の違いで説明できる
社会関係資本の蓄積の度合は、南北の地域の違いを反映している
分析において、南北の地域差を考慮する必要がある
北部　→　社会関係資本が多い　　→　政府のパフォーマンスが高い
南部　→　社会関係資本が少ない　→　政府のパフォーマンスが低い

検証 1: 南北によって政府のパフォーマンスは違うのか？

Question 1:
社会関係資本の蓄積の度合 (cc) は南北の地域の違い (location) を反映しているだけであって、政府のパフォーマンス (gov_p) とは無関係なのでは？

データ (putnam.csv)：

変数の種類	変数名	詳細
	`region`	イタリア州政府の略称
応答変数	`gov_p`	政府のパフォーマンス
説明変数	`cc`	Civic Community Index （市民共同体指標）
制御変数	`econ`	地方政府の経済指標（大きい程、経済が良好）
制御変数	`location`	イタリア北部地域ダミー（北部なら north、南部なら south）

✔ ここでは「制御変数」として新たに location というダミー変数をを付け加えている

このセクションの分析で使うデータ putnam.csvをダウンロードし、RProject Folder に保存する

df1 <- read_csv("data/putnam.csv")  # 欠損処理をしないので read_csv を使う
df1

# A tibble: 20 × 5
   region gov_p    cc  econ location
   <chr>  <dbl> <dbl> <dbl> <chr>   
 1 Ab       7.5   8     7   south   
 2 Ba       7.5   4     3   south   
 3 Cl       1.5   1     3   south   
 4 Cm       2.5   2     6.5 south   
 5 Em      16    18    13   north   
 6 Fr      12    17    14.5 north   
 7 La      10    13    12.5 north   
 8 Li      11    16    15.5 north   
 9 Lo      11    17    19   north   
10 Ma       9    15.5  10.5 north   
11 Mo       6.5   3.5   2.5 south   
12 Pi      13    15.5  17   north   
13 Pu       5.5   3.5   4   south   
14 Sa       5.5   8.5   8.5 south   
15 Si       4.5   3.5   5.5 south   
16 To      13    17.5  14.5 north   
17 Tr      11    18    12.5 north   
18 Um      15    15.5  11   north   
19 Va      10    15    15   north   
20 Ve      11    15    13.5 north

- 政府のパフォーマンスが南北地域間で異なるかどうか t 検定する

箱ひげ図を描く

df1 %>% 
  ggplot(aes(x = location, y = gov_p, fill = location)) +
    geom_boxplot() +
  labs(x = "イタリア北部- 南部(location)", y = "政府のパフォーマンス (gov_p)",
         title = "イタリアの地域別　政府のパフォーマンス") + 
  stat_smooth(method = lm, se = FALSE) +  # se = FALSE → 95% 信頼区間を消す
  theme_bw(base_family = "HiraKakuProN-W3")

視覚的に南北間で政府のパフォーマンスに明らかに差がある
政府のパフォーマンスのデータは unpaired だから default で t 検定する

t.test(df1$gov_p[df1$location == "north"],
       df1$gov_p[df1$location == "south"])


    Welch Two Sample t-test

data:  df1$gov_p[df1$location == "north"] and df1$gov_p[df1$location == "south"]
t = 6.8253, df = 14.552, p-value = 6.737e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 4.607777 8.808890
sample estimates:
mean of x mean of y 
 11.83333   5.12500

Question 1 結果・北部 gov_pの平均値　・・・11.833
・南部の gov_p の平均値・・・5.125
・その差 -6.708 は 1% 水準で統計的に有意 (p-value = 6.737e-06)
→ Goldberg (1996) が主張するように、政府のパフォーマンスには南北地域差あり

検証 2: 南北「内」でも「経済」は「政府のパフォーマンス」を説明できるのか？

Question 2:

経済の近代化の度合 (econ) が高い地域ほど政府のパフォーマンス (gov_p) が高い。これは南北地域「内」でもみられるか？
econ と gov_p の散布図を描く

df1 %>% 
  ggplot(aes(econ, gov_p)) +
  geom_point() +
  theme_bw() +
  labs(x = "経済状況 (econ)", y = "政府のパフォーマンス (gov_p)",
         title = "政府のパフォーマンスと経済状況") + 
  stat_smooth(method = lm, se = FALSE) +  # se = FALSE → 95% 信頼区間を消す
  theme_bw(base_family = "HiraKakuProN-W3")

経済の近代化の度合いと政府のパフォーマンス (gov_p) には正の相関がある
econ と gov_p には正の相関がある
- 地方政府の経済指標 (econ) が高い地域ほど、政府のパフォーマンス (gov_p) が高い
回帰式を求める

model_1 <- lm(gov_p ~ econ, data = df1)

summary(model_1)


Call:
lm(formula = gov_p ~ econ, data = df1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3386 -1.7733  0.0086  0.8336  5.5114 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.0108     1.3847   2.174 0.043264 *  
econ          0.5889     0.1200   4.909 0.000113 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.659 on 18 degrees of freedom
Multiple R-squared:  0.5724,    Adjusted R-squared:  0.5487 
F-statistic:  24.1 on 1 and 18 DF,  p-value: 0.0001131

\[\widehat{gov_p}\ = 3.01 + 0.589econ\] - df1 に含まれる変数の型をチェックする

str(df1)

spc_tbl_ [20 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ region  : chr [1:20] "Ab" "Ba" "Cl" "Cm" ...
 $ gov_p   : num [1:20] 7.5 7.5 1.5 2.5 16 12 10 11 11 9 ...
 $ cc      : num [1:20] 8 4 1 2 18 17 13 16 17 15.5 ...
 $ econ    : num [1:20] 7 3 3 6.5 13 14.5 12.5 15.5 19 10.5 ...
 $ location: chr [1:20] "south" "south" "south" "south" ...
 - attr(*, "spec")=
  .. cols(
  ..   region = col_character(),
  ..   gov_p = col_double(),
  ..   cc = col_double(),
  ..   econ = col_double(),
  ..   location = col_character()
  .. )
 - attr(*, "problems")=<externalptr>

df1 に含まれる変数 location は charactor なので、0, 1 のダミー変数 north を付け加える
変換後のデータフレームは df1 のまま

df1 <- df1 |> 
  mutate(north = if_else(location == "north", 1, 0)) 
                               # north = 1, south = 0 に変換 
head(df1) #変換したデータをの一部を表示

# A tibble: 6 × 6
  region gov_p    cc  econ location north
  <chr>  <dbl> <dbl> <dbl> <chr>    <dbl>
1 Ab       7.5     8   7   south        0
2 Ba       7.5     4   3   south        0
3 Cl       1.5     1   3   south        0
4 Cm       2.5     2   6.5 south        0
5 Em      16      18  13   north        1
6 Fr      12      17  14.5 north        1

- 経済の近代化の度合 (econ) が高い地域ほど、政府のパフォーマンス (gov_p) が高い
- これが南北地域「内」でもみられるかどうかを確認するためには、

→モデルに econ と north ダミーを同時に入れて重回帰分析を行う必要がある

model_2 <- lm(gov_p ~ econ + north, data = df1)

分析結果を表示する

summary(model_2)


Call:
lm(formula = gov_p ~ econ + north, data = df1)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6638 -1.1011 -0.2199  1.2497  4.1464 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  5.22207    1.34660   3.878  0.00121 **
econ        -0.01941    0.22037  -0.088  0.93083   
north        6.88386    2.22907   3.088  0.00667 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.19 on 17 degrees of freedom
Multiple R-squared:  0.7261,    Adjusted R-squared:  0.6939 
F-statistic: 22.53 on 2 and 17 DF,  p-value: 1.659e-05

この結果から、次の標本回帰関数 (SRF) の回帰式が得られる

\[\widehat{gov_p}\ = 5.22 - 0.019econ + 6.88north\]

経済の近代化の度合 (econ) は gov_p に有意な影響を与えていない
北部地域だと gov_p を 6.88 ポイント上げる 1% 水準で有意
north ダミーによって統制すると、次の 2 つの（傾きの等しい）回帰式を得る

南部 (north = 0) なら、

\[\widehat{gov_p}\ = 5.22 - 0.019econ\]

北部 (north = 1) なら、

\[\widehat{gov_p}\ = 12.11 - 0.019econ\]

以上の結果を散布図で表してみる

全体で見ると（上図左）
→ 経済の近代化の度合 (econ)が高くなるにつれて政府のパフォーマンス (gov_p) が高い
南北地域別に見ると（上図右）
→ 経済の近代化の度合 (econ) は政府のパフォーマンス (gov_p) に影響を与えていない
location は 1% 水準で統計的に有意 (p-value = 0.00667)
→　北部地域の方が南部地域より gov_p が 6.88 ポイント高い
econ は 5% 水準でも統計的に有意ではない
→ 南北地域ダミー (north) を考慮すると
→ 経済の近代化の度合 (econ) が政府のパフォーマンス (gov_p) に対する影響が消える

→ 経済の近代化の度合 (econ) と政府のパフォーマンス (gov_p)は偽の関係 (spurious correlation)

Question 2 結果・南北地域ダミー (north) を考慮すると、南北「内」で「経済の近代化の度合」は「政府のパフォーマンス」を説明できない

検証 3: 政府のパフォーマンスを説明するのは「経済」「南北」「社会関係資本」？

Question 3:
社会関係資本 (cc) が政府のパフォーマンス (gov_p) に影響を与えるという Putnam (1994) の主張は、南北地域差 (north) を考慮しても成立するのか？

cc と gov_p の散布図を描く
ggplot で日本語を表示したい場合には、マックユーザーは以下の行を入力

theme_bw(base_family = "HiraKakuProN-W3")

df1 %>% 
  ggplot(aes(cc, gov_p)) +
  geom_point() +
  labs(x = "市民共同体指標", y = "政府のパフォーマンス",
         title = "政府のパフォーマンスと市民共同体指標") + 
  stat_smooth(method = lm, se = FALSE)  + # se = FALSE → 95% 信頼区間を消す
  theme_bw(base_family = "HiraKakuProN-W3")

cc と gov_p には正の相関がある
- 社会関係資本 (cc) が高い地域ほど、政府のパフォーマンス (gov_p) が高い
回帰式を求めると

model_3 <- lm(gov_p ~ cc, data = df1)

summary(model_3)


Call:
lm(formula = gov_p ~ cc, data = df1)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5043 -1.3481 -0.2087  0.9764  3.4957 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.71115    0.84443   3.211  0.00485 ** 
cc           0.56730    0.06552   8.658 7.81e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.789 on 18 degrees of freedom
Multiple R-squared:  0.8064,    Adjusted R-squared:  0.7956 
F-statistic: 74.97 on 1 and 18 DF,  p-value: 7.806e-08

\[\widehat{gov_p}\ = 2.711 + 0.567cc\]

社会関係資本 (cc) が 1 単位増えると、政府のパフォーマンス (gov_p) が0.567ポイント上昇する
これは偶然とは考えられない（← p-value が 7.81e-08）
- これが南北地域「内」でもみられるかどうかを確認するためには、
→モデルに cc と north ダミーを同時に入れて重回帰分析を行う必要がある

model_4 <- lm(gov_p ~ cc + north, data = df1)

分析結果を表示する

summary(model_4)


Call:
lm(formula = gov_p ~ cc + north, data = df1)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5003 -1.3445 -0.2058  0.9773  3.4997 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  2.69850    1.12131   2.407   0.0278 *
cc           0.57094    0.21485   2.657   0.0166 *
north       -0.04781    2.67759  -0.018   0.9860  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.841 on 17 degrees of freedom
Multiple R-squared:  0.8064,    Adjusted R-squared:  0.7836 
F-statistic:  35.4 on 2 and 17 DF,  p-value: 8.689e-07

この結果から、次の標本回帰関数 (SRF) の回帰式が得られる
\[\widehat{gov_p}\ = 2.7 + 0.57cc - 0.048north\]
社会関係資本 (cc) は gov_p に有意な影響を与えている 5%水準で有意
north ダミーは gov_p に影響を与えていない
north ダミーによって統制すると、次の 2 つの（傾きの等しい）回帰式を得る

南部 (loation = 0) なら、

\[\widehat{gov_p}\ = 2.7 + 0.57cc\]

北部 (location = 1) なら、

\[\widehat{gov_p}\ = 2.65 + 0.57cc\]

全体で見ると（上図左）
→ 社会関係資本 (cc) が高くなるにつれて政府のパフォーマンス (gov_p) が高い
南北地域別に見ても（上図右）
→ 社会関係資本 (cc)が高くなるにつれて政府のパフォーマンス (gov_p) が高い
north は 5% 水準で統計的に有意ではない (p-value = 0.9860)
cc は 5% 水準でも統計的に有意 (p-value = 0.0166)
南北地域ダミー (north) を考慮しても、
→ 社会関係資本 (cc) が政府のパフォーマンス (gov_p) に対して影響を与えている

Question 3 結果・同レベルの社会関係資本 (cc) をもつ自治体では、北部であろうが南部であろうが、政府のパフォーマンス (gov_p) の高さに違いはない
・政府のパフォーマンス (gov_p) は「南北地域の違い」(north) によってではなく「社会関係資本」(cc) によって説明される

3. ダミー変数でわかること（総選挙）

日本の衆議院選挙結果データを使った次のモデルを考えてみよう

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

3.1.1 データのダウンロード方法

予めダウンロード先を指定する方法

RProject フォルダ内に data という名称のフォルダを作成する
下のコマンドを実行すると、csv ファイルがパソコンにダウンロードされ、data 内に自動的に保存される

download.file(url = "https://asanoucla.github.io/hr96-24.csv",
              destfile = "data/hr96-24.csv")

注意：一度ダウンロードを完了すれば、このコマンドを実行する必要はありません

ダウンロード先を指定しない方法

hr96-24.csv をクリックしてデータをパソコンにダウンロード　　
RProject フォルダ内に data という名称のフォルダを作成する
ダウンロードした hr96-24.csv を手動でRProject フォルダ内にある data フォルダに入れる

3.1.2 選挙データの読み取り方法

次のいずれかの方法で hr96-24.csv を読み取る

読み取り方法 1

na = "."というコマンドは「欠損値をドットで置き換える」という意味
欠損値を空欄のまま残すと、本来「数値 (numeric)」型のデータが「」文字型 (character)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処する

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

読み取り方法 2

読み取った値の日本語が文字化けする場合
locale()関数を使って日本語エンコーディング (cp932) を指定する

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

読み取り方法 3

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

3.1.3 読み取った選挙データを確認

hr96-24.csv は1996年に衆院選挙に小選挙区が導入されて以来実施された 10 回の衆議院選挙（1996, 2000, 2003, 2005, 2009, 2012, 2014, 2017, 2021, 2024）の結果のデータ
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"

df1 には 22 個の変数が入っている

変数名	詳細
year	選挙年 (1996-2024)
pref	都道府県名
ku	小選挙区名
kun	小選挙区
rank	当選順位
wl	選挙の当落: 1 = 小選挙区当選、2 = 復活当選、0 = 落選
nocand	立候補者数
seito	候補者の所属政党
j_name	候補者の氏名（日本語）
name	候補者の氏名（ローマ字）
previous	これまでの当選回数（当該総選挙結果は含まない）
gender	立候補者の性別: “male”, “female”
age	立候補者の年齢
exp	立候補者が使った選挙費用（総務省届け出）
status	候補者のステータス: 0 = 非現職、1 現職、2 = 元職
vote	得票数
voteshare	得票率 (%)
eligible	小選挙区の有権者数
turnout	小選挙区の投票率 (%)
seshu_dummy	世襲候補者ダミー: 1 = 世襲、0 = 非世襲（地盤世襲 or 非世襲）
jiban_seshu	地盤の受け継ぎ元の政治家の氏名と関係
nojiban_seshu	世襲元の政治家の氏名と関係

データの型をチェック

str(hr)

spc_tbl_ [10,773 × 22] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ year         : num [1:10773] 1996 1996 1996 1996 1996 ...
 $ pref         : chr [1:10773] "愛知" "愛知" "愛知" "愛知" ...
 $ ku           : chr [1:10773] "aichi" "aichi" "aichi" "aichi" ...
 $ kun          : num [1:10773] 1 1 1 1 1 1 1 2 2 2 ...
 $ wl           : num [1:10773] 1 0 0 0 0 0 0 1 0 2 ...
 $ rank         : num [1:10773] 1 2 3 4 5 6 7 1 2 3 ...
 $ nocand       : num [1:10773] 7 7 7 7 7 7 7 8 8 8 ...
 $ seito        : chr [1:10773] "新進" "自民" "民主" "共産" ...
 $ j_name       : chr [1:10773] "河村たかし" "今枝敬雄" "佐藤泰介" "岩中美保子" ...
 $ gender       : chr [1:10773] "male" "male" "male" "female" ...
 $ name         : chr [1:10773] "KAWAMURA, TAKASHI" "IMAEDA, NORIO" "SATO, TAISUKE" "IWANAKA, MIHOKO" ...
 $ previous     : num [1:10773] 2 2 2 0 0 0 0 2 0 0 ...
 $ age          : num [1:10773] 47 72 53 43 51 51 45 51 71 30 ...
 $ exp          : num [1:10773] 9828097 9311555 9231284 2177203 NA ...
 $ status       : num [1:10773] 1 2 1 0 0 0 0 1 2 0 ...
 $ vote         : num [1:10773] 66876 42969 33503 22209 616 ...
 $ voteshare    : num [1:10773] 40 25.7 20.1 13.3 0.4 0.3 0.2 32.9 26.4 25.7 ...
 $ eligible     : num [1:10773] 346774 346774 346774 346774 346774 ...
 $ turnout      : num [1:10773] 49.2 49.2 49.2 49.2 49.2 49.2 49.2 51.8 51.8 51.8 ...
 $ seshu_dummy  : num [1:10773] 0 0 0 0 0 0 0 0 1 0 ...
 $ jiban_seshu  : chr [1:10773] NA NA NA NA ...
 $ nojiban_seshu: chr [1:10773] NA NA NA NA ...
 - attr(*, "spec")=
  .. cols(
  ..   year = col_double(),
  ..   pref = col_character(),
  ..   ku = col_character(),
  ..   kun = col_double(),
  ..   wl = col_double(),
  ..   rank = col_double(),
  ..   nocand = col_double(),
  ..   seito = col_character(),
  ..   j_name = col_character(),
  ..   gender = col_character(),
  ..   name = col_character(),
  ..   previous = col_double(),
  ..   age = col_double(),
  ..   exp = col_double(),
  ..   status = col_double(),
  ..   vote = col_double(),
  ..   voteshare = col_double(),
  ..   eligible = col_number(),
  ..   turnout = col_double(),
  ..   seshu_dummy = col_double(),
  ..   jiban_seshu = col_character(),
  ..   nojiban_seshu = col_character()
  .. )
 - attr(*, "problems")=<externalptr>

数値は numeric 文字は character として認識されていることがわかる

3.2 データの作成と記述統計

hr96-24.csv は1996年に衆院選挙に小選挙区が導入されて以来実施された 9 回の衆議院選挙（1996, 2000, 2003, 2005, 2009, 2012, 2014, 2017, 2021）の結果データ
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"

hr には 22 個の変数が入っていることがわかる
これらの変数を使って、分析に必要な新たな変数 (exppv, ldp) をつくる

exppv （有権者一人あたり選挙費用）の作成

exp と eligible のデータ型を確認する

str(hr$exp)

 num [1:10773] 9828097 9311555 9231284 2177203 NA ...

str(hr$eligible)

 num [1:10773] 346774 346774 346774 346774 346774 ...

どちらも数値型 (numeric) で問題ない
exp と eligible を使って、有権者 1 人あたりに使う選挙費用 (exppv) を作る

hr <- hr %>% 
  dplyr::mutate(exppv = exp/eligible) # eligible は小選挙区ごとの有権者数

summary(hr$exppv)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
  0.0013   7.7494  17.4784  21.4130  30.2146 119.2479     1294

NA が 1974 もある！

ldp（自民党ダミー）の作成

seito に含まれる値を確かめる

unique(hr$seito)

 [1] "新進"                   "自民"                   "民主"                  
 [4] "共産"                   "文化フォーラム"         "国民党"                
 [7] "無所"                   "自由連合"               "政事公団太平会"        
[10] "社民"                   "新社会"                 "日本新進"              
[13] "新党さきがけ"           "青年自由"               "さわやか神戸・市民の会"
[16] "民主改革連合"           "市民新党にいがた"       "沖縄社会大衆党"        
[19] "緑の党"                 "公明"                   "諸派"                  
[22] "保守"                   "無所属の会"             "自由"                  
[25] "改革クラブ"             "保守新"                 "新党尊命"              
[28] "ニューディールの会"     "世界経済共同体党"       "新党日本"              
[31] "国民新党"               "新党大地"               "幸福"                  
[34] "みんな"                 "改革"                   "日本未来"              
[37] "日本維新の会"           "当たり前"               "アイヌ民族党"          
[40] "政治団体代表"           "安楽死党"               "次世"                  
[43] "維新"                   "生活"                   "立憲"                  
[46] "希望"                   "緒派"                   "Ｎ党"                  
[49] "国民"                   "れい"                   "参政"

分析には自民党ダミーが必要なので、if_else()関数を使って LDP = 1, それ以外の政党を 0 にしたダミー変数を作り hr に付け加える

hr <- hr %>% 
  mutate(ldp = if_else(seito == "自民", 1, 0))

小選挙区を表す変数 district の作成
- ku + kun => district
- 例）aichi + 1 => aichi_1

hr <- hr %>%
  mutate(district = str_c(ku, kun, sep = "_"))

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" "exppv"         "ldp"          
[25] "district"

最後に ldp が作られていることを確認
22個あった変数に 2 個加えたので、合計で 24個になっている
hr には次の 24 個の変数が入っている

変数名	詳細
year	選挙年 (1996-2021)
pref	都道府県名
ku	小選挙区名
kun	小選挙区
rank	当選順位
wl	選挙の当落: 1 = 小選挙区当選、2 = 復活当選、0 = 落選
nocand	立候補者数
seito	候補者の所属政党
j_name	候補者の氏名（日本語）
name	候補者の氏名（ローマ字）
previous	当選回数
gender	立候補者の性別: “male”, “female”
age	立候補者の年齢
exp	立候補者が使った選挙費用（総務省届け出）
status	候補者のステータス: 0 = 非現職、1 現職、2 = 元職
vote	得票数
voteshare	得票率 (%)
eligible	小選挙区の有権者数
turnout	小選挙区の投票率 (%)
seshu_dummy	世襲候補者ダミー: 1 = 世襲、0 = 非世襲（地盤世襲 or 非世襲）
jiban_seshu	地盤の受け継ぎ元の政治家の氏名と関係
nojiban_seshu	世襲元の政治家の氏名と関係
exppv	有権者一人あたりに費やす選挙費用（円）
ldp	自民党ダミー: 1 自民党候補者、0 = それ以外

このデータフレーム hr の中から2005年の衆議院選挙データだけ取り出し、次の 3 つの変数 (voteshare, exppv, ldp, district) を使って重回帰分析を行う
district は ku と kun を結合させて作成する

変数の種類	変数名	詳細
応答変数	`voteshare`	得票率 (%)
説明変数	`exppv`	有権者一人当たりに候補者が費やした選挙費用（円）
制御変数	`ldp`	自民党ダミー（自民党候補者 = 1、それ以外 = 0）
	`district`	小選挙区を表す変数：例）`aichi_1`

hr05 <- hr %>%　
  dplyr::filter(year == 2005) |> 　  　# 2005年のデータだけを選ぶ
  dplyr::mutate(district = str_c(ku, kun, sep = "_")) |> 
  dplyr::select(voteshare, exppv, ldp, district) # 3 つの変数だけを選ぶ

変数の中身を確認する

DT::datatable(hr05)

データのサマリーを表示

summary(hr05)

   voteshare         exppv             ldp           district        
 Min.   : 0.60   Min.   : 0.148   Min.   :0.0000   Length:989        
 1st Qu.: 8.80   1st Qu.: 8.352   1st Qu.:0.0000   Class :character  
 Median :34.80   Median :22.837   Median :0.0000   Mode  :character  
 Mean   :30.33   Mean   :24.627   Mean   :0.2932                     
 3rd Qu.:46.60   3rd Qu.:35.269   3rd Qu.:1.0000                     
 Max.   :73.60   Max.   :89.332   Max.   :1.0000                     
                 NA's   :4

text 表示する場合

stargazer(as.data.frame(hr05), type = "text")


==========================================
Statistic  N   Mean  St. Dev.  Min   Max  
------------------------------------------
voteshare 989 30.333  19.230  0.600 73.600
exppv     985 24.627  17.907  0.148 89.332
ldp       989 0.293   0.455     0     1   
------------------------------------------

html 表示する場合
R-Markdown で表示する際
→　type = "html" と指定、チャンクオプションで ```{r, results = "asis"} と指定

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


Statistic	N	Mean	St. Dev.	Min	Max

voteshare	989	30.333	19.230	0.600	73.600
exppv	985	24.627	17.907	0.148	89.332
ldp	989	0.293	0.455	0	1

html を綺麗に表示したいときには、下のコマンドを「チャンクの外に」貼り付けると、均整がとれた出力を得られる

<style>
table, td, th {
  border: none;
  padding-left: 1em;
  padding-right: 1em;
  min-width: 50%;
  margin-left: auto;
  margin-right: auto;
  margin-top: 1em;
  margin-bottom: 1em;
}
</style>

3.3 散布図の表示

exppv と voteshare の散布図を表示
ggplot で日本語を表示したい場合には、マックユーザーは以下の行を入力

theme_bw(base_family = "HiraKakuProN-W3")

hr05 %>% 
  ggplot(aes(exppv, voteshare)) +
  geom_point() +
  labs(x = "一人あたり選挙費用（円）", y = "得票率(%)",
         title = "候補者の得票率と選挙費用") +
  theme_bw(base_family = "HiraKakuProN-W3") + 
  geom_smooth(method = lm, se = FALSE)  # se = FALSE → 95% 信頼区間を消す

ldp と voteshare の散布図を表示

hr05 %>% 
  ggplot(aes(ldp, voteshare)) +
  geom_point() + 
  labs(x = "0 = 非自民党, 1 = 自民党", y = "得票率 (%)",
         title = "候補者の得票率と自民党所属") + 
  geom_jitter(width = 0.02) + # データを散らして表示  
  theme_bw(base_family = "HiraKakuProN-W3") +
  geom_smooth(method = lm, se = FALSE)  # se = FALSE → 95% 信頼区間を消す

3.4 ふたつの回帰分析 (Model_5 & Model_6)

voteshare と exppv の単回帰

model_5 <- lm(voteshare ~ exppv, 
              data = hr05)

voteshare と exppv + ldp の重回帰

model_6 <- lm(voteshare ~ exppv + ldp, 
              data = hr05)

3 つのモデルの結果を表示する
R-Markdown で表示する際、type = "html" を指定するときにはチャンクオプションで ```{r, results = "asis"} と指定する

stargazer(model_5, model_6, type = "html")


	Dependent variable:

	voteshare
	(1)	(2)

exppv	0.767^***	0.565^***
	(0.024)	(0.024)

ldp		15.852^***
		(0.961)

Constant	11.453^***	11.779^***
	(0.727)	(0.644)


Observations	985	985
R²	0.512	0.618
Adjusted R²	0.512	0.617
Residual Std. Error	13.422 (df = 983)	11.883 (df = 982)
F Statistic	1,031.634^*** (df = 1; 983)	794.203^*** (df = 2; 982)

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

model_6（右の列）の結果から、次の標本回帰関数 (SRF) の回帰式が得られる

\[\widehat{voteshare}\ = 11.78 + 0.57\cdot exppv + 15.85\cdot ldp\]

exppv（候補者が有権者一人当たりに使う選挙費用）が 1 円増えると、候補者の得票率が 0.57 % points 増える 1% で統計的に有意
自民党候補者は、そうでない候補者と比べると、平均して 15.85 % points 増える 1% で統計的に有意

自民党候補者でない (ldp = 0) なら、

\[\widehat{voteshare}\ = 11.78 + 0.57\cdot exppv\]

自民党候補者 (ldp = 1) なら、

\[\widehat{voteshare}\ = 27.63 + 0.57 \cdot exppv\]

この結果を図で表すと、傾きが 0.57 の二本の平行線になる

## 予測値のデータフレームを作り pred と名前をつける
pred <- with(hr05, expand.grid(
  exppv = seq(min(exppv, na.rm=TRUE), max(exppv, na.rm=TRUE), length = 100),
  ldp = c(0,1)
))

## mutate を使って、新たな変数である予測値 (pred) を作り計算する
pred <- mutate(pred, voteshare = predict(model_6, newdata = pred))

## 散布図のドット (point) は観測値 (hr05) で描き、回帰直線は予測値 (pred) で描く
p3 <- ggplot(hr05, aes(x = exppv, y = voteshare, color = as.factor(ldp))) +
  geom_point(size = 1) + geom_line(data = pred)
p3 <- p3 + labs(x = "一人あたり選挙費用", y = "得票率 (%)")
p3 <- p3 + scale_color_discrete(name = "所属政党", 
                                labels = c("非自民党","自民党")) + guides(color = guide_legend(reverse = TRUE)) +
  theme_bw(base_family = "HiraKakuProN-W3")

print(p3 + ggtitle("得票率と選挙費用の関係（候補者の所属政党ごと)"))

このモデルは 2 つの直線が平行になる（傾きが同じになる）ように設定されている
この図と分析結果から、候補者が有権者一人当たりに使う選挙費用が得票率を高めることと、自民党候補者なら16ポイント近く得票率が高いことがわかる
しかし、選挙費用が得票率に与える影響の大きさが、自民党候補者とそうでない候補者の間でどれだけ違うのかはわからないこのモデルでは、影響の大きさ（つまり傾き）は同じと想定しているため　
交差項をモデルに加えることで、自民党候補者とそうでない候補者それぞれの影響の大きさ（つまり傾き）を知ることができる

3.5 ダミー変数を含む回帰分析 (Model_6) のメカニズム

ダミー変数を含めた上記重回帰分析 (model_6) のメカニズムを理解するため、ここでは次の3 段階の回帰分析を試してみる

3 段階の回帰分析
Step1: 連続変数とダミー変数を使った単回帰分析 (voteshare, ldp)
Step2: 連続変数と連続変数の単回帰分析 (voteshare, exppv)
Step3: ダミー変数を含めた重回帰分析 (voteshare, exppv, ldp)

Step1: 連続変数とダミー変数を使った単回帰分析

候補者の得票率 (voteshare) を応答変数、候補者が自民党公認か否か (ldp) を説明変数とした回帰分析を考える
voteshare は連続変数、ldp はカテゴリカル変数
voteshare を縦軸、ldp を横軸にして両変数の相関関係を散布図で表す

ggplot(hr05, aes(ldp, voteshare)) +
  geom_point() +
  labs(x = "0 = 非自民党, 1 = 自民党", y = "得票率(%)",
         title = "候補者の得票率と候補者の自民党所属") + 
  geom_jitter(width = 0.02) + # データを散らして表示  
  geom_smooth(method = lm, se = FALSE) + # se = FALSE → 95% 信頼区間を消す
  theme_bw(base_family = "HiraKakuProN-W3")

このモデルの回帰式は次のようになる

\[voteshare_i \sim \mathrm{N}(\alpha_0 + \alpha_1 \cdot ldp_i, \sigma^2)\]

lm() を使って回帰式を求め、model_7と名前をつける

model_7<- lm(voteshare ~ ldp, data = hr05)

回帰分析結果は summary(　) を使って表示できる

summary(model_7)


Call:
lm(formula = voteshare ~ ldp, data = hr05)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.223 -14.313  -0.423  12.387  46.187 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.4132     0.5593   40.07   <2e-16 ***
ldp          27.0096     1.0329   26.15   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.79 on 987 degrees of freedom
Multiple R-squared:  0.4093,    Adjusted R-squared:  0.4087 
F-statistic: 683.8 on 1 and 987 DF,  p-value: < 2.2e-16

以上の結果から、$\widehat{\alpha}_0= 22.41$, $\widehat{\alpha}_1 = 27$, $\widehat{\sigma} = 14.79$ が得られる
この結果を踏まえると、model_7 の回帰式は次のようになる

\[\widehat{voteshare}\ = 22.41 + 27 \cdot ldp\]

model_7 の分析結果の解釈

F 検定の結果は最下行に示されており、 p 値は 2.2e-16 なので F 検定は 1% で有意
帰無仮説は「ldp の係数は 0 」
$Pr ( > | t | )$ 下の「2e-16」が p 値
→ p 値が 0.01 より小さいので、有意水準 $α = 0.01 (1 %) $で帰無仮説が棄却
→ 候補者が自民党候補者だと (ldp = 1) 得票率が 27 %ポイント上がる
Adjusted R-squared の値が 0.4087 → 応答変数 (voteshare) の分散の 40.87% が ldp によって説明できた
R では Signif. codes: 以下に有意水準の表示基準が示され、p 値右上に付されたアステリスク数によって、帰無仮説が棄却される有意水準が示されている
アステリスクク数が 3 つ・・・有意水準 $α = 0.001 (0.1\%)$ で有意
アステリスクク数が 2 つ・・・有意水準 $α = 0.01 (1\%)$ で有意
アステリスクク数が 1 つ・・・有意水準 $α = 0.05 (5\%)$ で有意
stargazer()関数を使うと結果を見やすく表示できる
R-Markdown で表示する際、type = "html" を指定するときにはチャンクオプションで ```{r, results = "asis"} と指定する

stargazer(model_7,
          type = "html",
          dep.var.labels = "Voteshare (%)", # 応答変数の名前
          title = "Table 1: LDP Affiliation and Voteshare in the 2005HR Election") # タイトル

**Table 1: LDP Affiliation and Voteshare in the 2005HR Election**

	Dependent variable:

	Voteshare (%)

ldp	27.010^***
	(1.033)

Constant	22.413^***
	(0.559)


Observations	989
R²	0.409
Adjusted R²	0.409
Residual Std. Error	14.788 (df = 987)
F Statistic	683.759^*** (df = 1; 987)

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

ダミー変数を使った単回帰分析と t 検定

この回帰直線の切片である13.98は、非自民党候補者の平均得票率（予測得票率）で、ldp = 0 を代入して計算できる:

$22.41 + 27.01 \cdot 0 = 22.41$

他方、予測値の式に ldp = 1 を代入すると、自民党候補者の平均得票率（予測得票率）が得られる:

\[22.41 + 27.01 \cdot 1 = 49.42\]

R で候補者ごとに平均得票率を求め、上で求めた2 つの予測値と一致するか確かめよう
非自民党候補者 (ldp = 0) の平均得票率（予測得票率）

hr05 %>%
  filter(ldp == 0) %>%  # filter を使って ldp = 0 だけのデータに限定
  with(mean(voteshare)) %>%    # with を使って votehshare の平均値を計算する
  round(2)                     # 小数点 2 位まで表示

[1] 22.41

自民党候補者 (ldp = 1) の平均得票率（予測得票率）

hr05 %>%
  filter(ldp == 1) %>%         # filter を使って ldp = 1 だけのデータに限定
  with(mean(voteshare)) %>%    # with を使って votehshare の平均値を計算する
  round(2)                     # 小数点 2 位まで表示

[1] 49.42

以上から、予測値は説明変数の値を与えられたときの、応答変数の平均値であることがわかる
実際に、自民党候補者と非自民党候補者の voteshare の平均を t 検定すると回帰分析と同様の結果が得られる
voteshare は unpaired だから default で t 検定する

t.test(hr05$voteshare[hr05$ldp == 1], 
       hr05$voteshare[hr05$ldp == 0])


    Welch Two Sample t-test

data:  hr05$voteshare[hr05$ldp == 1] and hr05$voteshare[hr05$ldp == 0]
t = 32.446, df = 897.39, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 25.37584 28.64335
sample estimates:
mean of x mean of y 
 49.42276  22.41316

Step2: 連続変数と連続変数を使った単回帰分析

データフレーム hr05 を使って、候補者の得票率 (voteshare) を応答変数、候補者が使う選挙費用 (exppv: 単位は円) を説明変数とした回帰分析を考える
voteshare も exppv もどちらも連続変数
voteshare を縦軸、exppv を横軸にして両変数の相関関係を散布図で表す

ggplot(hr05, aes(exppv, voteshare)) +
  geom_point() +
  labs(x = "一人あたり選挙費用（円）", y = "得票率(%)",
         title = "候補者の得票率と選挙費用") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  geom_smooth(method = lm, se = FALSE)  # se = FALSE → 95% 信頼区間を消す

lm() を使って回帰式を求め、model_8 と名前を付ける

model_8 <- lm(voteshare ~ exppv, data = hr05)

回帰分析結果を表示する

summary(model_8)


Call:
lm(formula = voteshare ~ exppv, data = hr05)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.596  -9.661  -3.493   9.909  42.851 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.45334    0.72744   15.74   <2e-16 ***
exppv        0.76745    0.02389   32.12   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.42 on 983 degrees of freedom
  (4 observations deleted due to missingness)
Multiple R-squared:  0.5121,    Adjusted R-squared:  0.5116 
F-statistic:  1032 on 1 and 983 DF,  p-value: < 2.2e-16

この結果を踏まえると、model_8 の回帰式は次のようになる
\[\widehat{voteshare}\ = 11.45 + 0.77 \cdot exppv\]

model_8 の分析結果の解釈

$F$ 検定結果は最下行に示されており $p$ 値は 2.2e-16 なので $F$ 検定は 1% で有意
帰無仮説は「exppv の係数は 0 」
$Pr ( > | t | )$ 下の「2e-16」が $p$ 値
→ $p$ 検定値が 0.01 より小さいので、有意水準 $α = 0.01 (1 %)$ で帰無仮説が棄却
→ 候補者が選挙費用を 1 円使うと、得票率が 0.767 %ポイント増える
Adjusted R-squared の値が 0.512 → 応答変数 (voteshare) の分散の 51.2% が exppv によって説明できた
stargazer()関数を使うと結果を見やすく表示できる
R-Markdown で表示する際、type = "html" を指定するときにはチャンクオプションで ```{r, results = "asis"} と指定する

stargazer(model_8,
          type = "html",　
          dep.var.labels = "Voteshare (%)", # 応答変数の名前
          title = "Table 2: Money and Voteshare in the 2005HR Election") # タイトル

**Table 2: Money and Voteshare in the 2005HR Election**

	Dependent variable:

	Voteshare (%)

exppv	0.767^***
	(0.024)

Constant	11.453^***
	(0.727)


Observations	985
R²	0.512
Adjusted R²	0.512
Residual Std. Error	13.422 (df = 983)
F Statistic	1,031.634^*** (df = 1; 983)

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

Step3: ダミー変数を含めた重回帰分析

【ポイント：自民党候補者であることは候補者の得票率に影響したのか？】

応答変数を voteshare（得票率）とし、ldp（自民党候補者）と exppv（有権者一人当たりの選挙費用）という2 つの説明変数で説明する重回帰分析を試みる
モデルには model_9 という名前をつける
- model_9 の前提:「自民党候補者であるかどうかによって得票率が異なるにしても、選挙費用が得票率に与える影響（傾きの大きさ）は同じ」

model_9 <- lm(voteshare ~ ldp + exppv, data = hr05)

# 分析結果を表示する
summary(model_9)


Call:
lm(formula = voteshare ~ ldp + exppv, data = hr05)

Residuals:
    Min      1Q  Median      3Q     Max 
-42.990  -8.828  -2.427   9.284  46.290 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.77885    0.64431   18.28   <2e-16 ***
ldp         15.85224    0.96087   16.50   <2e-16 ***
exppv        0.56538    0.02444   23.13   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.88 on 982 degrees of freedom
  (4 observations deleted due to missingness)
Multiple R-squared:  0.618, Adjusted R-squared:  0.6172 
F-statistic: 794.2 on 2 and 982 DF,  p-value: < 2.2e-16

この結果を踏まえると、model_9 の回帰式は次のようになる

\[\widehat{voteshare}\ = 11.78 + 15.85 \cdot ldp\ + 0.565 \cdot exppv`\]

model_9 の分析結果の解釈
$F$ 検定結果は最下行に示されており、$p$ 値は 2.2e-16 なので$F$ 検定は 1% で有意

帰無仮説は次の 2 つ：

$H_0$:「候補者が自民党公認か否かという変数 (ldp) の係数は 0」
$H_0$:「候補者が使う有権者一人あたりの選挙費用 (exppv) の係数は 0」

$Pr ( > | t | )$ 下の「2e-16」が $p$ 値
→ 両変数の $p$値が 0.01 より小さいので、有意水準 $α$ = 0.01 (1 %) で帰無仮説が棄却
→ exppv の値を平均に固定した時、自民党公認候補者だと (ldp = 1) 得票率が 15.85 %ポイント上がる
→ ldp の値を平均に固定した時、候補者が有権者一人当たりに選挙費用を 1 円使うと、得票率が 0.565 %ポイント上がる
切片の係数が 11.78 だから、有権者一人当たりに選挙費用 expp と自民党ダミー ldp の両者が 0 のとき、得票率の予測値は 11.78 %ポイント
Adjusted R-squared の値が 0.6172 → 応答変数 (voteshare) の分散の 61.72% が exppv と ldp によって説明できた
ldp の係数 15.85 とはどういう意味か？
→ 「ldp ダミーが 1 単位増えると、得票率が 15.85%ポイント上昇する」
＝「自民党公認候補者は、そうでない候補者と比べて得票率が 15.85%ポイント高い」
model_9 の結果を図示すると、次のようになる

## 予測値のデータフレームを作り pred と名前をつける
pred <- with(hr05, expand.grid(
  exppv = seq(min(exppv, na.rm=TRUE), max(exppv, na.rm=TRUE), length = 100),
  ldp = c(0,1)
))

## mutate を使って、新たな変数である予測値 (pred) を作り計算する
pred <- mutate(pred, voteshare = predict(model_9, newdata = pred))

## 散布図のドット (point) は観測値 (hr05) で描き、回帰直線は予測値 (pred) で描く
p9 <- ggplot(hr05, aes(x = exppv, y = voteshare, color = as.factor(ldp))) +
  geom_point(size = 1) + geom_line(data = pred)
p9 <- p9 + labs(x = "一人あたり選挙費用（円）", y = "得票率(%)")
p9 <- p9 + scale_color_discrete(name = "ldp candidate", 
                                labels = c("no","yes")) + guides(color = guide_legend(reverse = TRUE)) +
theme_bw(base_family = "HiraKakuProN-W3") 
  
print(p9 + ggtitle("選挙費用と得票率（自民党所属ごと）"))

このモデルは 2 つの直線が平行になる（傾きが同じになる）ように設定されている

model_9 の回帰

\[\widehat{voteshare} = 11.78 + 15.85 \cdot ldp\ + 0.565 \cdot exppv\]

model_9 の回帰式に、非自民党候補者 (ldp = 0) を代入すると、トマト色の回帰式が得られる

\[\widehat{voteshare}\ = 11.78 + 0.565 \cdot exppv\]

model_9 の回帰式に、自民党候補者 (ldp = 1) を代入すると、青色の回帰式が得られる

\[\widehat{voteshare}\ = 27.63 + 0.565 \cdot exppv\]

ldp ダミーの値が異なる2 つの式を比較すると、違うのは切片だけで、傾きは同じ
2 つの回帰直線は平行だが、その切片は異なる(11.78 と 27.63)
自民党候補者であるという要因が、回帰直線を上方に移動させている
（= 自民党候補者は、そうでない候補者と比べると、平均して得票率が 15.85ポイント増える)
全体的に自民党公認候補者の得票率（青色）が、そうでない候補者の得票率（トマト色）より高い
2 本の回帰直線は、それぞれ自民党公認候補者（青色）とそうでない候補者（トマト色）の、選挙費用と得票率の関係を捉えている
ダミー変数を説明変数に加えると、ダミー変数の値によって回帰直線が平行移動するような変化（つまり、自民党候補者であるか否かという要因は候補者の得票率に影響したのかということ）を捉えることができる
stagazer() を使うと、3 つの統計モデル (model_7, model_8, model_9)の分析結果を一つの表で表すこともできる
stargazer(model_7, model_8, model_9, type = "html")という簡単なコマンドでも結果を表示できる
R-Markdown で表示する際、type = "html" を指定するときにはチャンクオプションで ```{r, results = "asis"} と指定する

stargazer(model_7, model_8, model_9, 
          type = "html", 　　　                    
          dep.var.labels = "Voteshare (%)", # 応答変数の名前
          title = "Table 4: 選挙費用と得票率（自民党所属ごと）") # タイトル

**Table 4: 選挙費用と得票率（自民党所属ごと）**

	Dependent variable:

	Voteshare (%)
	(1)	(2)	(3)

ldp	27.010^***		15.852^***
	(1.033)		(0.961)

exppv		0.767^***	0.565^***
		(0.024)	(0.024)

Constant	22.413^***	11.453^***	11.779^***
	(0.559)	(0.727)	(0.644)


Observations	989	985	985
R²	0.409	0.512	0.618
Adjusted R²	0.409	0.512	0.617
Residual Std. Error	14.788 (df = 987)	13.422 (df = 983)	11.883 (df = 982)
F Statistic	683.759^*** (df = 1; 987)	1,031.634^*** (df = 1; 983)	794.203^*** (df = 2; 982)

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

`modelsummary()` を使った結果の表示

・stargazer の代わりに、modelsummary() を使って、分析結果を表示させる
・modelsummary パッケージを読み込む

library(modelsummary)

msummary(list(model_7, model_8, model_9),
  gof_omit = "IC|F|RMSE|Log",
    stars = TRUE)

	(1)	(2)	(3)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	22.413***	11.453***	11.779***
	(0.559)	(0.727)	(0.644)
ldp	27.010***		15.852***
	(1.033)		(0.961)
exppv		0.767***	0.565***
		(0.024)	(0.024)
Num.Obs.	989	985	985
R2	0.409	0.512	0.618
R2 Adj.	0.409	0.512	0.617

3.6 95 パーセント信頼区間を加えた散布図

上記の散布図に 95 パーセント信頼区間を加えてみる
信頼区間を求めるために標準誤差を利用する
標準誤差は t 分布に従うので、qt( ) で分布の 95% が収まる範囲を求める

## 信頼区間の下限と上限を計算する
## まず、予測値 (pred) と標準誤差 (err) を求める
err <- predict(model_9, newdata = pred, se.fit = TRUE)

## 予測値と標準誤差を使って信頼区間 (pred$lower, pred$upper) を求める
pred$lower <- err$fit + qt(0.025, df = err$df) * err$se.fit
pred$upper <- err$fit + qt(0.975, df = err$df) * err$se.fit

p9.ci95 <- p9 +
  geom_smooth(data = pred, aes(ymin = lower, ymax = upper), stat ="identity")

# 作成した図を表示する
print(p9.ci95 + ggtitle("選挙費用と得票率（自民党所属ごと）"))

3.7 頑強な標準誤差を使った回帰分析

R上の回帰分析で使う lm() では、不均一分散に対応できない
結果変数の分散が一定でないと
→　lm()が計算する標準誤差は回帰係数のばらつきを過小評価してしまう
不均一分散と理論的な対処法に関する詳細に関しては森氏幸之助「回帰分析ではlm()ではなく estimtr::lm_robust()を使おう」参照

解決策 → `estimatr::lm_robust()` を使う

頑強な標準誤差 (Robust Standard Error) で不均一分散に対処する
Table 4 で示した選挙費用と得票率の分析結果をさらに精緻化してみる
285の小選挙区の事情を考慮した分析
=> 小選挙区で cluster 化する
model_9 を小選挙区で cluster 化して回帰分析する
→ model_10 と名前を付けて保存
lm_robust() を使うために必要なパッケージを読み込む

library(estimatr)

DT::datatable(hr05)

lm() を使った普通の重回帰分析

model_9 <- lm_robust(voteshare ~ ldp + exppv, 
  data = hr05)

lm_robust() を使って、289の小選挙区ごとにクラスター化した重回帰分析

model_10 <- lm_robust(voteshare ~ ldp + exppv, 
  data = hr05,
  clusters = district)

stargazerを使って結果を表記してみる

stargazer::stargazer(model_9, model_10,
  type = "text")


% Error: Unrecognized object type.
% Error: Unrecognized object type.

エラーのため、結果が表示されない
stargazer は lm_robust() に対応していない

解決策 => `lm_robust()` に対応している `modelsummary` を使う

library(modelsummary)

modelsummary(list(model_9, model_10), 
  stars = TRUE,               # アステリスクを表示させる  
  statistic = 'std.error',    # 標準誤差はデフォルト  
  gof_omit = "IC|F|RMSE|Log") # 非表示の設定)

	(1)	(2)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	11.779***	11.779***
	(0.632)	(0.500)
ldp	15.852***	15.852***
	(0.933)	(0.852)
exppv	0.565***	0.565***
	(0.026)	(0.025)
Num.Obs.	985	985
R2	0.618	0.618
R2 Adj.	0.617	0.617
Std.Errors		by: district

lm()が計算する標準誤差は回帰係数のばらつきを過小評価している
lm で推定した (1) のモデルの標準誤差 (括弧内の数字) が lm_robust で推定した (2) のモデルの標準誤差よりも大きい
例）ldp の標準誤差 (0.933) > (0.852)
=> (2) のモデルの方が不均一分散に対処したよりよい推定

変数名の変更とタイトル

変数名の変更は coef_map を使う
あらかじめ、変数名の対応表をベクトルに指定しておく
変数を表示する順番も指定できる

var_nam = c("ldp"       = "自民党ダミー",
            "exppv"     = "一人あたり選挙費（円）",
            "(Intercept)" = "Constant")

modelsummary(list(model_9, model_10), 
  stars = TRUE,               # アステリスクを表示させる  
  statistic = 'std.error',    # 標準誤差はデフォルト  
  gof_omit = "IC|F|RMSE|Log", # 非表示の設定
  title = "2005年総選挙における得票率の分析",
  coef_map = var_nam)

2005年総選挙における得票率の分析
	(1)	(2)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
自民党ダミー	15.852***	15.852***
	(0.933)	(0.852)
一人あたり選挙費（円）	0.565***	0.565***
	(0.026)	(0.025)
Constant	11.779***	11.779***
	(0.632)	(0.500)
Num.Obs.	985	985
R2	0.618	0.618
R2 Adj.	0.617	0.617
Std.Errors		by: district

出力したい統計量の設定・statistic = 'p.value'・・・p値を出力したい時
・statistic = 'statistic'・・・t値を出力したい時
・statistic = 'conf.int'・・・95%信頼区間を出力したい時
・statistic = 'conf.int', conf_level = .99・・・99%信頼区間を出力したい時

4. Exercise

hr96-24.csv は1996年に衆院選挙に小選挙区が導入されて以来実施された 9 回の衆議院選挙（1996, 2000, 2003, 2005, 2009, 2012, 2014, 2017, 2021, 2024）の結果のデータ
このデータの中から2005年の衆議院選挙データだけ取り出し、次の 3 つの変数 (voteshare, exppv, jcp) を使って重回帰分析を行いなさい

変数の種類	変数名	詳細
応答変数	`voteshare`	得票率 (%)
説明変数	`exppv`	有権者一人当たりに候補者が費やした選挙費用（円）
制御変数	`jcp`	共産党ダミー（共産党候補者 = 1、それ以外 = 0）

注意：jcp ダミーは自分で作成すること
Q1: stargazer()を使って、2005年の衆議院選挙の 3 つの変数 (voteshare, exppv, jcp) の記述統計を表示させなさい
Q2: exppv と voteshare の散布図を描きなさいその際、回帰直線も表示すること
Q3: jcp と voteshare の散布図を描きなさいその際、回帰直線も表示すること
Q4: voteshare を応答変数、exppv と jcp を説明変数とした重回帰分析を実行し、その重回帰式を示し、重回帰分析結果を解釈しなさい
Q5: voteshare を応答変数、exppv と jcp を説明変数とした散布図を描きなさいその際、観測値と回帰直線はダミー変数 (jcp) で色分けすること
（白黒で描く場合には、散布図のドットの形を変えたり、回帰直線の種類を 2 種類使うこと）

参考文献

飯田健『計量政治分析』共立出版、2013年.

Ellis Goldberg (1996), Thinking about How Democracy Works, Politics & Society, Vol. 24, pp.7-18.

宋財泫 (Jaehyun Song)- 矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」

矢内勇生（高知工科大学）授業一覧

浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年

浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年

Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.

Kieran Healy, DATA VISUALIZATION, Princeton, 2019

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017

15. 回帰分析 2（ダミー変数）

Masahiko Asano

2024-12-17

1. ダミー変数でわかること

1.1 ダミー変数とは

1.2 ダミー変数でわかること

2. 南北格差ダミーと経済の近代化

2.1 理論と仮説

2.2 oldberg による Putnam 批判とその検証

検証 1: 南北によって政府のパフォーマンスは違うのか？

検証 2: 南北「内」でも「経済」は「政府のパフォーマンス」を説明できるのか？

検証 3: 政府のパフォーマンスを説明するのは「経済」「南北」「社会関係資本」？

3. ダミー変数でわかること（総選挙）

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

3.1.1 データのダウンロード方法

予めダウンロード先を指定する方法

ダウンロード先を指定しない方法

3.1.2 選挙データの読み取り方法

読み取り方法 1

読み取り方法 2

読み取り方法 3

3.1.3 読み取った選挙データを確認

3.2 データの作成と記述統計

3.3 散布図の表示

3.4 ふたつの回帰分析 (Model_5 & Model_6)

3.5 ダミー変数を含む回帰分析 (Model_6) のメカニズム

`modelsummary()` を使った結果の表示

3.6 95 パーセント信頼区間を加えた散布図

3.7 頑強な標準誤差を使った回帰分析

解決策 → `estimatr::lm_robust()` を使う

解決策 => `lm_robust()` に対応している `modelsummary` を使う

変数名の変更とタイトル

4. Exercise

15. 回帰分析 2（ダミー変数）

Masahiko Asano

2024-12-17

1. ダミー変数でわかること

1.1 ダミー変数とは

1.2 ダミー変数でわかること

2. 南北格差ダミーと経済の近代化

2.1 理論と仮説

2.2 oldberg による Putnam 批判とその検証

検証 1: 南北によって政府のパフォーマンスは違うのか？

検証 2: 南北「内」でも「経済」は「政府のパフォーマンス」を説明できるのか？

検証 3: 政府のパフォーマンスを説明するのは「経済」「南北」「社会関係資本」？

3. ダミー変数でわかること（総選挙）

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

3.1.1 データのダウンロード方法

予めダウンロード先を指定する方法

ダウンロード先を指定しない方法

3.1.2 選挙データの読み取り方法

読み取り方法 1

読み取り方法 2

読み取り方法 3

3.1.3 読み取った選挙データを確認

3.2 データの作成と記述統計

3.3 散布図の表示

3.4 ふたつの回帰分析 (Model_5 & Model_6)

3.5 ダミー変数を含む回帰分析 (Model_6) のメカニズム

modelsummary() を使った結果の表示

3.6 95 パーセント信頼区間を加えた散布図

3.7 頑強な標準誤差を使った回帰分析

解決策 → estimatr::lm_robust() を使う

解決策 => lm_robust() に対応している modelsummary を使う

変数名の変更とタイトル

4. Exercise

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

`modelsummary()` を使った結果の表示

解決策 → `estimatr::lm_robust()` を使う

解決策 => `lm_robust()` に対応している `modelsummary` を使う