このセクションで使う R パッケージ一覧

library(broom)
library(car)
library(datasets)
library(fBasics)
library(grid)
library(gridExtra)
library(jtools)
library(patchwork)
library(QuantPsyc)
library(stargazer)
library(summarytools)
library(tidyverse)

1. 回帰分析

回帰分析の目的：

1. 現象の「説明」

2. 未来の「予測」

3. x の値に人為的に「介入」したときの効果

4. 制御：y を特定の範囲に収めるためには x の値をどうするか

1.1 単回帰分析 (Simple Linear Regression)

一つの説明変数と一つの応答変数との「直線的な」関係を求め、説明変数から応答変数を推定する方法

回帰分析における 2 つの変数は様々な名称で呼ばれるので注意が必要

影響を与える変数 (X)	影響を与えられる変数 (Y)
独立変数 (`independent variable`)	従属変数 (`dependent variable`)
説明変数 (`explanatory variable`)	応答変数・目的変数 (`response variable`)
説明変数 (`explanatory variable`)	被説明変数 (`explained variable`)
リグレッサー (`regressor`)	リグレッサンド (`regressand`)

条件付き期待値を 1 次関数としてモデル化し、誤差項（ε）を付けたものが単回帰モデル
ある母集団における 2 つの変数 Y と X の間の関係が次のような母集団回帰関数 (PRF: population regression function) で表せるとする

Y: outcome (response variable, dependent variable)
X: predictor (explanatory variable, independent variable)
\(α\) と \(β\)：パラメータ（係数: coefficients)
\(α\)：X = 0 の時の Y の平均値
\(β\)：X が 1 単位増加した時の Y の増加分の平均値
\(ε\)：error term (disturbance) = 誤差項
　　観測値が回帰直線から逸脱している程度を示す
　　偶然によって生じる確率変数
母集団のパラメータ \(α\) と \(β\) の値は通常わからない
集められたサンプルデータから推定 (estimate) する

母集団と標本

母集団における 2 つの変数 Y と X の間の関係 \(Y = α + βX + ε\)
→　\(X\) を条件とした時の \(Y\) の期待値: \(E[Y|X] = α + βX\)
・母集団からサンプルをとった後は、母集団のパラメータ \(α\) と \(β\) の推定値は \(\hat{α}\) や \(\hat{β}\) と表し、それぞれ「アルファハット」「ベータハット」と呼ぶ
予測値 \(\hat{Y}\) は次のように表す

\[\hat{Y} = \hat{α} + \hat{β}x\]

この回帰直線を使って、x に対応する予測値\(\hat{Y}\) (predicted value) を計算できる
通常、予測値 \(\hat{Y}\)は観測値 \(Y\)と完全に一致しない
予測値 \(\hat{Y}\) は観測値 \(Y\) との差 \(\hat{ε}\)（エプシロンハット）で表す

\[\hat{ε} = Y - \hat{Y}\]

観測値 \(Y\) は、予測値 \(\hat{Y}\) に誤差項 \(\hat{ε}\) を加えて次のように表せる

\[Y = \hat{α} + \hat{β}X + \hat{ε}\]

\(\hat{ε}\) は「残差」(residual or predictin error) と呼ばれる
「残差」は母集団の error term の推定値 (estimate) なので、ハットを付ける
以上の関係を図で表すと次のようになる

最小二乗法 (`OLS: Ordinary Least Squares`)

残差平方和を最小化する切片 \(\hat{α}\) と傾き \(\hat{β}\) を使って線を引く方法

\[\hat{α} = \bar{Y} - \hat{β}\bar{X}\]

\[\hat{β} = \frac{\sum (X_i - \bar{X})(Y_i-\bar{Y})}{\sum(X_i - \bar{X})^2}\]

予測値 \(\hat{Y}\)は次の式で表すことができる

\[\hat{Y} = \hat{α} + \hat{β}x \]

\(x\) の値がその平均値 \(\bar{X}\) の時の \(Y\) の値を調べてみる
\(x\) → \(\bar{X}\)と置き換えると

\[\hat{Y} = \hat{α} + \hat{β}\bar{X}\]

\(\hat{α} = \bar{Y} - \hat{β}\bar{X}\) をこの式に代入すると

\[\hat{Y} = (\bar{Y} - \hat{β}\bar{X}) + \hat{β}\bar{X} = \bar{Y}\]

\[\bar{Y} - \hat{β}\bar{X} = \hat{α}\]

\(x\) の値がその平均値 \(\bar{X}\) の時、\(Y\) もその平均値 \(\bar{Y}\) になる
つまり、回帰直線は常にデータの平均値

\[(\bar{X}, \bar{Y})\]

を通る

係数を推定するために最小二乗法を使う際、回帰直線に基づく予測は、平均すると正確である
その根拠 → 残差 ( \(\hat{ε}\) ) の平均はゼロになるから

1.2 手計算で単回帰式を計算してみる

手計算で単回帰式を計算

ここでは次のデータを使って、手計算で単回帰式を計算してみる
動画による解説はこちら
単回帰式の求め方をわかりやすく解説しています
この動画で使っている csv.file はこちらからダウンロードできます

Source: https://youtu.be/GhrxgbQnEEU

1.3 回帰係数に関する仮説

帰無仮説:

\[H_0:β = 0\]

対抗仮説:

\[H_A: β ≠ 0\]

回帰式 \(Y_i = α + βX_i + ε_i\) において \(β=0\)　→　「傾きが」
→　「\(X_i\)から\(Y_i\)への線形の影響がない」という意味

2. 実証分析のプロセス

① パズルを見つける
② 理論の提示（パズルのメカニズムの説明）
③ 仮説（検証可能な仮説「もしこの理論が正しければ・・・のはずである」）
④ 仮説の検定

命題が「真」であることを直接的に証明することはできない
⇒ 否定されるべき仮説（＝帰無仮説）を立てる
　（パソコンの統計ソフトは帰無仮説を自動的に設定）
⇒ 帰無仮説を棄却する
⇒ 理論の正しさを間接的に示す

イタリアにおける地方政府のパフォーマンスを事例として考えてみる

2.1 パズルを見つける

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

データ：

region : イタリアの地方政府の略称
gov_p : 政府のパフォーマンス

出典：飯田健『計量政治分析』p.28の図を修正

2.2 理論の提示

パズルのメカニズムの説明
ここでは次の説明（＝理論）を使う

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

出典：Robert Putnam, (1994) Making Democracy Work: Civic Traditions in Modern Italy,
Princeton, NJ: Princeton University Press,
(河田潤ー訳『哲学する民主主義:伝統と改革の市民的構造』NTT 出版、2001 年)

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

2.4 仮説の提示

観察可能な応答変数と説明変数を示し、検証可能な仮説「もしこの理論が正しければ　・・・　のはずである」を明示する

【応答変数】　　

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

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

【説明変数】　　

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

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

仮説：

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

3. 回帰分析による仮説検証

3.1 データの読み取り (`putnam.csv`)

上記の政府のパフォーマンスと市民共同体指標の回帰分析を R を使って再現してみる
まず政府のパフォーマンスに関するデータ (putnam.csv) をダウンロードし、RProject フォルダー内の data フォルダに保存する
データのダウンロードが終わったら、データを読み込み putnam と名前を付ける

putnam <- read_csv("data/putnam.csv")

データフレーム putnam に含まれている変数を表示する

names(putnam)

[1] "region"   "gov_p"    "cc"       "econ"     "location"

データ：

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

データフレーム putnam の中を確認

DT::datatable(putnam)

3.2 記述統計量を示す

putnam に含まれている変数の記述統計量を示す:summary() を使う方法
手軽に記述統計量をチェックしたいとき

summary(putnam)

    region              gov_p             cc              econ      
 Length:20          Min.   : 1.50   Min.   : 1.000   Min.   : 2.50  
 Class :character   1st Qu.: 6.25   1st Qu.: 3.875   1st Qu.: 6.25  
 Mode  :character   Median :10.00   Median :15.000   Median :11.75  
                    Mean   : 9.15   Mean   :11.350   Mean   :10.43  
                    3rd Qu.:11.25   3rd Qu.:16.250   3rd Qu.:14.50  
                    Max.   :16.00   Max.   :18.000   Max.   :19.00  
   location        
 Length:20         
 Class :character  
 Mode  :character

putnam に含まれている変数の記述統計量を示す: stargazer() を使う方法
stargazer() を使うと、記述統計を見やすく表示できる

library("stargazer")

タイトルを付けることもできる

stargazer(as.data.frame(putnam), 
          type = "text", 
          title = "パットナムデータの記述統計", 
          digits = 2)


パットナムデータの記述統計
======================================
Statistic N  Mean  St. Dev. Min   Max 
--------------------------------------
gov_p     20 9.15    3.96   1.50 16.00
cc        20 11.35   6.26   1.00 18.00
econ      20 10.43   5.08   2.50 19.00
--------------------------------------

さらに covariate.labels を使って変数名を変更したり out を使って作成したテーブルを保存できる

stargazer(as.data.frame(putnam), 
          type = "text", 
          title = "パットナムデータの記述統計", 
          digits = 2,
          
covariate.labels = c("政府のパフォーマンス", # 変数名を変更  
                   "市民共同体指標",
                   "経済指標"))


パットナムデータの記述統計
=======================================
Statistic  N  Mean  St. Dev. Min   Max 
---------------------------------------
政府のパフォーマンス 20 9.15    3.96   1.50 16.00
市民共同体指標    20 11.35   6.26   1.00 18.00
経済指標       20 10.43   5.08   2.50 19.00
---------------------------------------

スタイルを指定し、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>

次のコマンドを コードチャンク内に打ち込む
R-Markdown で表示する際、type = "html" と指定し、チャンクオプションで{r, results = "asis"} と指定する

stargazer(as.data.frame(putnam), 
          type = "html", 
          title = "パットナムデータの記述統計", 
          digits = 2, 
          
covariate.labels = c("政府のパフォーマンス", # 変数名を変更  
                   "市民共同体指標",
                   "経済指標"))

**パットナムデータの記述統計**

Statistic	N	Mean	St. Dev.	Min	Max

政府のパフォーマンス	20	9.15	3.96	1.50	16.00
市民共同体指標	20	11.35	6.26	1.00	18.00
経済指標	20	10.43	5.08	2.50	19.00

3.3 散布図を描く

R 上で gov_p と cc の散布図を描いてみる

Macユーザーのみ

チャンク内で下のコマンドを最初に実行する

# Macでggplot2内で日本語が使える様にする
theme_set(theme_bw(base_family = "HiraKakuProN-W3"))

Windows ユーザーは次のコマンドを入力する

putnam %>% 
  ggplot(aes(cc, gov_p)) +
  geom_point() +
  labs(x = "市民共同体指標 (cc)", y = "政府のパフォーマンス (gov_p)",
         title = "政府のパフォーマンスと市民共同体指標の散布図") + 
  stat_smooth(method = lm) +  # (method = lm, se = FALSE) → 95% 信頼区間が消える
  geom_text(aes(y = gov_p + 0.2, 
                label = region), 
            size = 4, 
            vjust = 0)

gov_p と cc の相関係数と p-value を求める

cor.test(putnam$gov_p, putnam$cc)


    Pearson's product-moment correlation

data:  putnam$gov_p and putnam$cc
t = 8.6583, df = 18, p-value = 7.806e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7558102 0.9593028
sample estimates:
      cor 
0.8979883

gov_p と cc の相関関係は 1% の有意水準 (p-value = 8.806e-08) で統計的に有意
強い正の相関 (相関係数 r = 0.898) がある
しかし、相関関係は2 つの変数間の「共変関係」を意味するだけ
→ 因果関係を探るためには、cc を説明変数、gov_p を応答変数とした回帰分析が必要

3.4 単回帰分析モデル

gov_p を応答変数、cc を説明変数とした単回帰分析モデルを Model_1 と名前をつけて実行する

Model_1 <- lm(gov_p ~ cc, data = putnam)

3.4.1 回帰分析の結果の表示方法

回帰分析の結果の表示方法は複数あるが、ここでは代表的な 3 つの方法を紹介する
目的に応じて使い分ける

回帰分析の結果の表示方法: summary()関数

summary(Model_1)


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

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

回帰分析の結果の表示方法: tidy()関数

tidy(Model_1)

# A tibble: 2 × 5
  term        estimate std.error statistic      p.value
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>
1 (Intercept)    2.71     0.844       3.21 0.00485     
2 cc             0.567    0.0655      8.66 0.0000000781

回帰分析の結果の表示方法: stargazer()関数
- 回帰分析の結果を見やすい表にまとめる時は stargazer() 関数を使う
- R-Markdown で html 表示する際にはチャンクオプションで {r, results = "asis"} と指定

stargazer(Model_1, type = "html")


	Dependent variable:

	gov_p

cc	0.567^***
	(0.066)

Constant	2.711^***
	(0.844)


Observations	20
R²	0.806
Adjusted R²	0.796
Residual Std. Error	1.789 (df = 18)
F Statistic	74.967^*** (df = 1; 18)

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

この結果から、次の標本回帰関数 (SRF) の回帰式を求めることができる

\[\hat{gov_p}\ = 2.7 + .567cc\]

参考:

分析結果で表示する項目の詳細を設定できる

stargazer(Model_1,     
          digits = 2,      　　　　　　　　　# 小数点以下 2 位まで示す
          keep.stat = c("n", "adj.rsq", "f"),# 表示する統計量を指定する
          dep.var.caption = "応答変数", 　  # 変数の種類を明示
          dep.var.labels = "政府のパフォーマンス", # 応答変数の名前
          title = "Model 1: 単回帰分析結果", # タイトル
          type ="html")

**Model 1: 単回帰分析結果**

	応答変数

	政府のパフォーマンス

cc	0.57^***
	(0.07)

Constant	2.71^***
	(0.84)


Observations	20
Adjusted R²	0.80
F Statistic	74.97^*** (df = 1; 18)

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

表示する統計量は keep.stat = c("...", "...", "...") を使って、次の一覧表から選ぶことが可能
有意確率を示すアステリスクを表示したくなければ、star.cutoffs = NA と omit.table.layout = “n” を付け加える

3.4.2 結果の解釈

tidy() を使うと簡素な分析結果を表示できる

tidy(Model_1)

# A tibble: 2 × 5
  term        estimate std.error statistic      p.value
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>
1 (Intercept)    2.71     0.844       3.21 0.00485     
2 cc             0.567    0.0655      8.66 0.0000000781

cc の係数 (.567) の解釈:
市民共同体指標 (cc) の値が 1 ポイントだけ異なる地方を比べると、cc の値が大きい地方が、平均して 0.567ポイントだけ政府のパフォーマンス (gov_p) が高い
しかし、この結果は「標本における結果」
→　知りたいのは、母集団におけるパラメタ \({β_1}\)
→ 次の帰無仮説と対立仮説を設定する

帰無仮説　\(H_0: β_1 = 0\)
対立仮説　\(H_a: β_1 ≠ 0\)

帰無仮説の意味
→ 母集団における cc の傾き（\(β_1\)) は 0
→ つまり「cc は gov_p に影響を与えない」
対立仮説の意味
→ 母集団における cc の傾き（\(β_1\)) は 0 ではない
→ つまり「cc が gov_p に与える影響を与えている」
次の式を使って t 値を求め仮説検定を行う

\[t = \frac{\hat{β_1} - β_1}{SE}\]

\({SE}\) は \(\hat{β_1}\) の標準誤差 (Std. Error) =　推定値（係数）の標準偏差のこと
→　推定値は平均値であり、散らばりがある
→　その散らばりの程度を表すのが標準誤差（= 推定値の標準偏差）
→　サンプルから得られた推定値（\(\hat{β_1}\)）は 0.5673
→　ここでの帰無仮説は「母集団における cc の傾き（\({β_1}\)）は 0 」

→ Std. Error = 0.06552 と傾き \(\hat{β_1}= 0.5673\)、 \(β_1= 0\) を代入すると t 値が得られる

\[t = \frac{\hat{β_1} - β_1}{SE} = \frac{0.5673 - 0}{0.06552}= 8.658\]

このサンプルから t = 8.658 という統計量が得られた

帰無仮説 \(H_0\) が正しい（つまり、cc の gov_p に対する影響はない）とき、
ここで求めた t 値（t = 8.658）より絶対値が大きい t 値を与えるようなデータが出現する確率
= \(p\) 値

単回帰分析 Model_1 の結果を summary(Model_1) で確認
tidy(Model_1）でも（簡素版の）同じ結果が得られる

分析結果画面の cc と t valueの交差する数字 (8.658) が \(t\) 値
分析結果画面の cc と \(Pr(> | t |)\)の交差する数字 (7.81e-08) が\(p\)値
\(p\)値が 7.81e-08 = 0.000000781 （ゼロが 8 つ）で、0.001より小さい
→ 帰無仮説 \(H_0: β_1 = 0\) を 0.1% の有意水準で棄却できる

【解説】

帰無仮説 \(H_0\) が正しい（つまり、cc の gov_p に対する影響はない）とき、単純無作為抽出で同じ大きさの標本を選ぶと、0.000000781 の確率で t 値の絶対値が 8.656 より大きいデータが得られる
帰無仮説 \(H_0\) が正しい（つまり、cc の gov_p に対する影響はない）とき、私たちが扱っているようなデータが出現する確率は 0.000000781（＝ほとんどない）
→ 帰無仮説 \(H_0\) が間違っている
→ 帰無仮説 \(H_0\) を棄却する

結論 「市民共同体指標 (cc) の値が大きい地方ほど、政府のパフォーマンス (gov_p) の値が高い」

R-squared の値が 0.8064
→ 地方政府のパフォーマンス (gov_p) の分散の約 81% が市民共同体指標 (cc) の分散で説明できた
→　だからといって、市民共同体指標が地方政府のパフォーマンスを引き起こしているとは限らない
R-squared = 0.8064 が意味すること
→このモデルがどれだけ上手に応答変数の変化を「説明」しているか
＝R-squared =「そのモデルがどれだけ正確か」についての指標
単回帰分析の限界・・・一つの説明変数以外の要因が含まれていない → これは問題
ほとんどの現象は一つ以上の要因から影響をうけているはず
セレクションバイアスがある可能性が高い

解決策 → 重回帰分析
重回帰分析は他の要因の影響をコントロールした分析を可能にしてくれる
→　ある程度、セレクションバイアスの影響を回避できる
→　重回帰分析は最も基本的なセレクションバイアスの削減方法

手計算で単回帰式を計算してみる（確認）

手計算で単回帰式を計算

ここでは次のデータを使って、手計算で単回帰式を計算してみる
手計算で使う xlsx.file と模範解答はこちらからダウンロードできます

最小二乗法による予測値 \(\hat{Y}\) の式

\[\hat{Y} = \hat{α} + \hat{β}x\]

傾き \(β\) は次の式で表せる

\[\hat{β} = r \frac{s_y}{s_x}\]

\(r\) →　ピアソンの相関係数
\(s_y\) →　y の標準偏差
\(s_x\) →　x の標準偏差

切片 \(α\) は次の式で表せる

\[\hat{α} = \bar{Y} - β\bar{x}\]

ピアソンの相関係数 \(r\) は次の式で求められる

\[r = \frac{\sum((x-\bar{x})(y-\bar{y}))}{\sqrt{Σ(x-\bar{x})^2Σ(y-\bar{y})^2}}= \frac{422.95}{\sqrt{(745.55)(297.55)}}= \frac{422.95}{{471}}= 0.898\]

\(s_y\) →　 y の標準偏差

\[s_y= \sqrt{\frac{Σ(y-\bar{y})^2}{{n-1}}}= \sqrt{\frac{297.55}{{19}}}=3.957\]
\(s_x\) →　 x の標準偏差

\[s_x= \sqrt{\frac{Σ(x-\bar{x})^2}{{n-1}}}= \sqrt{\frac{745.55}{{19}}}=6.264\]
\[\hat{β} = r \frac{s_y}{s_x}= 0.898*\frac{3.956}{6.264}=0.5671\]

\[\hat{α} = \bar{Y} - β\bar{x}= 9.15 - 0.5671*11.35 =2.71\]

\[\hat{Y}=2.71+0.567x\]

4. 重回帰モデル

単回帰分析の結果「市民共同体指標 (cc) の値が 1 ポイント高いと、政府のパフォーマンス (gov_p) の値が 0.567 ポイント高く」この結果は 1% の有意水準で統計的に有意であることがわかった
しかし、ある結果がたった一つの原因で引き起こされているとは限らない
むしろ、複数の原因が結果に作用していることが多い
結果に影響を与えると思われる他の要因をコントロールする必要がある
それを可能にさせるのが重回帰分析
複数の説明変数 (cc と econ) と一つの応答変数 (gov_p) との「直線的な」関係を考える

説明変数は「主要な説明変数 (cc)」と「コントロール変数 (econ)」から構成される
ここでは最も大事な変数を cc と想定し、econ は「主要な説明変数 (cc) 以外に gov_p に影響を与えると思われる変数」（=コントロール変数）と位置づけている
もし cc も econ も同レベルで重要な変数と想定するなら、cc と econ の両方が「主要な説明変数」になる
分析で使う変数は次のとおり

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

重回帰分析の流れ

重回帰式を求める意義を調べるため説明変数と応答変数の散布図を描く
重回帰式 \(\hat{Y} = \hat{α} + \hat{β_1}\ cc + \hat{β_2}\ econ\) を求める（最小二乗法）
偏回帰係数（\(α, β_1, β_2\)）の検定
モデル全体の有意性検定（F 検定）
重回帰式の精度 = 決定係数（\(R^2\)）を確認する
修正済み決定係数 (Adjusted R-squared)を確認
標準偏回帰係数（β ：ベータ）を確認

この手続きに従って、イタリアにおける地方政府のパフォーマンスに関するデータを使って重回帰分析を行う　　

4.1 散布図

分析で使う変数は次の三つ

変数の種類	変数名	詳細
応答変数	`gov_p`	政府のパフォーマンス
説明変数	`cc`	Civic Community Index （市民共同体指標）
コントロール変数	`econ`	地方政府の経済指標: 大きい程経済が良好

重回帰式を求める意義を調べるため次の2 つの散布図を描く

gov_p と cc
gov_p と econ

散布図を描く前に分析で使う変数の記述統計を確認

データの記述統計は次のとおり

stargazer(as.data.frame(putnam), 
          type = "text", 
          title = "パットナムデータの記述統計", 
          digits = 2)


パットナムデータの記述統計
======================================
Statistic N  Mean  St. Dev. Min   Max 
--------------------------------------
gov_p     20 9.15    3.96   1.50 16.00
cc        20 11.35   6.26   1.00 18.00
econ      20 10.43   5.08   2.50 19.00
--------------------------------------

R 上で gov_p と cc の散布図を描く

plt_1 <- putnam |> 
  ggplot(aes(cc, gov_p)) +
  geom_point() +
  labs(x = "市民共同体指標 (cc)", 
       y = "政府のパフォーマンス (gov_p)",
       title = "政府のパフォーマンスと市民共同体指標の散布図") + 
  stat_smooth(method = lm) +  # (method = lm, se = FALSE) → 95% 信頼区間が消える
  geom_text(aes(y = gov_p + 0.2, 
                label = region), 
            size = 4, 
            vjust = 0) +
  theme_set(theme_bw(base_family = "HiraKakuProN-W3"))

plt_2 <- putnam |> 
  ggplot(aes(econ, gov_p)) +
  geom_point() +
  labs(x = "経済指標 (econ)", 
       y = "政府のパフォーマンス (gov_p)",
       title = "政府のパフォーマンスと経済指標の散布図") + 
  stat_smooth(method = lm) +  # (method = lm, se = FALSE) → 95% 信頼区間が消える
  geom_text(aes(y = gov_p + 0.2, 
                label = region), 
            size = 4, 
            vjust = 0)

plt_1 + plt_2

どちらも正の関係があることがわかる

4.2 二つの単回帰分析

gov_p を応答変数、cc を応答変数とした単回帰分析モデルを Model_1 と名前をつける

Model_1 <- lm(gov_p ~ cc, data = putnam)

gov_p を応答変数、econ を応答変数とした単回帰分析モデルを Model_2 と名前をつける

Model_2 <- lm(gov_p ~ econ, data = putnam)

分析結果の表示法: stargazer()を使う

回帰分析の結果をシンプルな表にまとめる時は stargazer() 関数を使う
分析結果で表示する項目の詳細を設定できる
stargazer(Model_1, type = "html") という簡単なコマンドでも結果を表示できる
R-Markdown で表示する際にはチャンクオプションで {r, results = "asis"} と指定する

stargazer(Model_1, Model_2,     
          digits = 3,         # 小数点以下 3 位まで示す
          style = "ajps", # AJPS仕様  
          dep.var.caption = "応答変数", # 変数の種類を明示
          dep.var.labels = "政府のパフォーマンス", # 応答変数の名前
          title = "分析結果", # タイトル
          type ="html")

**分析結果**

	政府のパフォーマンス
	Model 1	Model 2

cc	0.567^***
	(0.066)
econ		0.589^***
		(0.120)
Constant	2.711^***	3.011^**
	(0.844)	(1.385)
N	20	20
R-squared	0.806	0.572
Adj. R-squared	0.796	0.549
Residual Std. Error (df = 18)	1.789	2.659
F Statistic (df = 1; 18)	74.967^***	24.097^***

p < .01; p < .05; p < .1

統計的有意水準に関して Stargazer() では自動的にアステリスクが付される
p < .01（1% 水準で有意）- - - 「\(***\)」
p < .05（5% 水準で有意）- - - 「\(**\)」
p < .1 （10% 水準で有意）- - - 「\(*\)」

分析結果 (Model_1)

「cc が 1 ポイント高いと gov_p が 0.567 ポイント高い」
cc の係数にアステリスクが三つ (***) 付いている
→ この結果は 1% の有意水準で統計的に有意

分析結果 (Model_2)

「econ が 1 ポイント高いと gov_p が 0.589 ポイント高い」
econ の係数にアステリスクが三つ (***) 付いている
→ この結果は 1% の有意水準で統計的に有意

ポイント：cc も econ も単回帰分析では統計的にそれぞれ有意

参考：

t 値や p 値を確認したければ、summary()関数や tidy() 関数を使う

tidy(Model_1, conf.int = TRUE) # 95%信頼区間を表示

# A tibble: 2 × 7
  term        estimate std.error statistic      p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>    <dbl>     <dbl>
1 (Intercept)    2.71     0.844       3.21 0.00485         0.937     4.49 
2 cc             0.567    0.0655      8.66 0.0000000781    0.430     0.705

tidy(Model_2, conf.int = TRUE) # 95%信頼区間を表示

# A tibble: 2 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)    3.01      1.38       2.17 0.0433      0.102     5.92 
2 econ           0.589     0.120      4.91 0.000113    0.337     0.841

4.3 単回帰＋単回帰＝重回帰？

複数の単回帰分析を同時に実行する = 重回帰分析
単回帰分析では x 軸と y 軸の二次元分析だが、ここでは gov_p軸、cc軸、econ 軸の三次元分析になる
三次元分析することで、cc と econ が gov_p に与える影響を「同時に」確認することができ「変数をコントロール」できる
ここで、gov_p を応答変数、cc を説明変数、econ をコントロール変数とした重回帰分析モデルを Model_3 と名前をつけて、結果を表示する

3つのモデルの分析結果を表示

Model_3 <- lm(gov_p ~ cc + econ, data = putnam)
summary(Model_3)


Call:
lm(formula = gov_p ~ cc + econ, data = putnam)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2593 -0.6329 -0.0574  0.7944  2.8674 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.2358     0.9122   3.547  0.00248 ** 
cc            0.7538     0.1520   4.958  0.00012 ***
econ         -0.2534     0.1873  -1.353  0.19388    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.749 on 17 degrees of freedom
Multiple R-squared:  0.8252,    Adjusted R-squared:  0.8046 
F-statistic: 40.13 on 2 and 17 DF,  p-value: 3.645e-07

この重回帰分析の結果をシンプルな表にまとめる時は stargazer() 関数を使う
R-Markdown で表示する際にはチャンクオプションで ```{r, results = "asis"} と指定する

stargazer(Model_1, Model_2, Model_3,     
          digits = 3,  　# 小数点以下 3 位まで示す
          style = "ajps", # AJPS仕様  
          covariate.labels = c("市民共同体指標 (cc)", "経済指標 (econ)"),
          dep.var.caption = "応答変数", # 変数の種類を明示
          dep.var.labels = "政府のパフォーマンス", # 応答変数の名前
          title = "分析結果（イタリア州政府のパフォーマンス）", # タイトル
          type ="html")

**分析結果（イタリア州政府のパフォーマンス）**

	政府のパフォーマンス
	Model 1	Model 2	Model 3

市民共同体指標 (cc)	0.567^***		0.754^***
	(0.066)		(0.152)
経済指標 (econ)		0.589^***	-0.253
		(0.120)	(0.187)
Constant	2.711^***	3.011^**	3.236^***
	(0.844)	(1.385)	(0.912)
N	20	20	20
R-squared	0.806	0.572	0.825
Adj. R-squared	0.796	0.549	0.805
Residual Std. Error	1.789 (df = 18)	2.659 (df = 18)	1.749 (df = 17)
F Statistic	74.967^*** (df = 1; 18)	24.097^*** (df = 1; 18)	40.126^*** (df = 2; 17)

p < .01; p < .05; p < .1

有意確率を示すアステリスクを表示したくなければ、star.cutoffs = NA と omit.table.layout = “n” を付け加える

Model_3 の重回帰式を求める（最小二乗法）

上記 R による重回帰分析での結果から、Model_3 の回帰式は次の様になる

\[\hat{gov_p}\ = 3.236 + 0.754cc - 0.253econ\]

4.4 重回帰分析のメカニズム

上の例では、重回帰分析は複数の単回帰分析を「同時に行っている」ようにみえるが、そうではない
つまり、応答変数 \(\hat{y}\) と 2 つの説明変数 \(x_1, x_2\) の重回帰

\[\hat{gov_p}\ = 3.236 + 0.754cc - 0.253econ\] と 2 つの単回帰

\[\hat{gov_p} = 2.711 + 0.567cc\]
\[\hat{gov_p} = 3.011 + 0.589econ\]
を考えたとき、重回帰分析で得られる説明変数の係数と、2 つの単回帰分析で得られる係数は一致していない

	`cc の係数`	`econ の係数`
単回帰	0.567	0.589
重回帰	0.754	0.253

重回帰分析で使った複数の説明変数のそれぞれを使って単回帰分析を行っても、重回帰分析と同じ結果(係数)を得ることはできない

重回帰で行っていること = 二段階推定

説明変数同士で回帰して推定する
→　残差を計算する
応答変数を残差に回帰する
→　交絡因子を取り除いた影響力を推定できる

重回帰分析における `cc の係数 (0.754)` について考える

\[\hat{gov_p}\ = 3.236 + 0.754cc - 0.253econ\]

Step 1:

2 つの共変量を使って、cc を econ に回帰する

step1 <- lm(cc ~ econ, data = putnam)
coef(step1)

(Intercept)        econ 
 -0.2984271   1.1173551

下図左側の ESS が「econ によって説明された cc の部分 (= 1.117)」
単回帰における econ の係数 (1.117)が得られた
これより次の回帰式が回帰直線が得られる

\[\hat{cc} = -0.298 + 1.117econ\]

「econ によって説明されれなかった cc の残り部分」が残差: e_cc

Step 2:

e_cc を使って gov_p を説明する
resid()関数を使って step1 の残差 e_cc を計算し、中身を確認する

e_cc <- resid(step1)
e_cc

         1          2          3          4          5          6          7 
 0.4769413  0.9463618 -2.0536382 -4.9643812  3.7728106  1.0967779 -0.6685119 
         8          9         10         11         12         13         14 
-1.0205772 -3.9313201  4.0661984  1.0050393 -3.1966099 -0.6709934 -0.6990914 
        15         16         17         18         19         20 
-2.3470260  1.5967779  4.3314881  3.5075208 -1.4618997  0.2141330

gov_p を e_cc に回帰する

step2 <- lm(gov_p ~ e_cc, data = putnam)
coef(step2)

(Intercept)        e_cc 
   9.150000    0.753799

ここで得られた e_cc の係数 (=0.754) が、重回帰分析における cc の係数 (= 0.754)
→　econ をコントロールしたときの cc の gov_p に対する影響力

通常、このような二段階推定はせずに重回帰分析を実行する

model <- lm(gov_p ~ cc + econ,
            data = putnam)
coef(model)

(Intercept)          cc        econ 
  3.2357962   0.7537990  -0.2533731

confint(model, level = 0.95)

                 2.5 %    97.5 %
(Intercept)  1.3112706 5.1603219
cc           0.4330438 1.0745542
econ        -0.6485673 0.1418211

二段階推定（もしくは重回帰分析）を実行することで、交絡因子を取り除いて cc の影響力を推定できた

4.5 モデル全体の有意性検定（F 検定）

Model_3 の結果最後にある F 値と p 値を確認

 F-statistic: 40.13 on 2 and 17 DF,  p-value: 3.645e-07

F-statistic は cc と econ が同時に 0 であるかどうか検定できる
ここでの帰無仮説:「切片以外の全てのパラメータ (cc と econ) が同時に 0」
F-statistic は「40.13」、p-value は「3.645e-07」
→この結果は 1% の有意水準で統計的に有意
→市民共同体指標 (cc) の影響と経済指標 (econ) の影響が同時にないとは言えない」
→母集団を表すモデルとして適切である

F 検定の解説・F 検定では複数の帰無仮説を同時に検定する
・複数の帰無仮説を同時に検定する理由:
→ 1 つずつ検定すると有意でないが、全体で検定する有意なことがありうるから

検定の方法:

→ 2 つのモデルを作る
(1) 無制限モデル (unrestricted model) — cc と econ を含むモデル
(2) 制限モデル (restricted model) — cc と econ を含まないモデル
→　それぞれの \(R^2\) を計算
→　「cc と econ を含むことによって、 \(R^2\) が十分に増えたか？」を検定
・この場合に使われる F 値 (F-statistic) は 2 つの自由度によって形が決まる F 分布に従う

→　F 検定は「統計的に有意でない」変数をモデルから外すべきかどうかの判断に有益　　

4.6 偏回帰係数の検定（`Model_3`）

「cc が 1 ポイント高いと gov_p が 0.754 ポイント高い」
cc の p-value は 0.00012 なので、この結果は 1% の有意水準で統計的に有意である
「econ が 1 ポイント高いと gov_p が 0.2534 ポイント低い」
econ の p-value は 0.19388 なので、この結果は統計的に有意ではない

統計的に有意ではないことの意味：
→ 点推定量は 0.19388と推定されるが、これは誤差の範囲と区別がつかない

大事な点 - econ と gov_p を単回帰分析したときには econ は統計的に有意
- cc と gov_p を単回帰分析したときには cc は統計的に有意
- cc を含めた重回帰分析では、gov_p に対する econ の影響は消えた
→ econ と gov_p は spurious correlation（疑似相関）
→ gov_p に影響を与えているのは cc

4.7 修正済み決定係数を確認

\(R^2\) は「決定係数」「寄与率」とも呼ばれる
重回帰分析では \(R^2\) ではなく Adjusted R-squared（修正済み \(R^2\)）が使われる

 Adjusted R-squared:  0.8046

応答変数 (gov_p) の分散の80.46% は cc と econ の分散によって説明される
→修正済み \(R^2\) は、そのモデルがどれだけ上手に応答変数の変化を「説明」しているかの指標
=「そのモデルがどれだけ正確か」についての指標
この場合、このモデルで gov_p の分散の約20%が説明されないまま (unexplained) 残るという意味

4.8 標準偏回帰係数（`β`：ベータ）を確認　　

通常、説明変数の単位は異なるため（例えば、一つは cm でもう一つは kg）、変数の係数値を比較するだけでは、どちらの変数が応答変数により大きな影響を与えているか比較できない
β値 (Standardized Beta Coefficient)とは、説明変数の値を１標準偏差増やした時に、応答変数の標準偏差がどれだけ変化するかを示す
β値を比較することで、複数の説明変数のうち、応答変数に対してどちらがより大きな影響を与えているのかわかる

library("QuantPsyc")  # β値を計算するパッケージ

Model_3 <- lm(gov_p ~ cc + econ, data = putnam)
lm.beta(Model_3)

        cc       econ 
 1.1932019 -0.3255233

「cc が 1 標準偏差高いと gov_p が 1.19標準偏差増える」
「econ が 1 標準偏差高いと gov_p が 0.33標準偏差減る」

4.9 重回帰分析結果の表示方法

従来の表示方法

**分析結果（イタリア州政府のパフォーマンス）**

	政府のパフォーマンス
	Model 1	Model 2	Model 3

市民共同体指標 (cc)	0.567^***		0.754^***
	(0.066)		(0.152)
経済指標 (econ)		0.589^***	-0.253
		(0.120)	(0.187)
Constant	2.711^***	3.011^**	3.236^***
	(0.844)	(1.385)	(0.912)
N	20	20	20
R-squared	0.806	0.572	0.825
Adj. R-squared	0.796	0.549	0.805
Residual Std. Error	1.789 (df = 18)	2.659 (df = 18)	1.749 (df = 17)
F Statistic	74.967^*** (df = 1; 18)	24.097^*** (df = 1; 18)	40.126^*** (df = 2; 17)

p < .01; p < .05; p < .1

注意：

ここでは 3 つのモデルの分析結果を表示している
その理由は、回帰分析のメカニズムを説明する教育的配慮のため
しかし、実際に学術論文に分析結果を示す場合には、表記方法に工夫が必要
例えば、最も大切な変数が cc だとする
→　私たちが知りたいのは「他のコントロール変数をモデルに含めても cc が応答変数と関連があるかどうか」ということ
→　全てのモデルに最も大事な変数 (cc) が入っている必要がある
→　この場合、結果の表示は次のようにする

**分析結果（イタリア州政府のパフォーマンス）**

	政府のパフォーマンス
	Model 1	Model 2

市民共同体指標 (cc)	0.567^***	0.754^***
	(0.066)	(0.152)
経済指標 (econ)		-0.253
		(0.187)
Constant	2.711^***	3.236^***
	(0.844)	(0.912)
N	20	20
R-squared	0.806	0.825
Adj. R-squared	0.796	0.805
Residual Std. Error	1.789 (df = 18)	1.749 (df = 17)
F Statistic	74.967^*** (df = 1; 18)	40.126^*** (df = 2; 17)

p < .01; p < .05; p < .1

このように表記すると「変数 cc は単回帰 (Model_1) でも統計的に有意だし、econ というコントロール変数を含めたモデル (Model_2) でも有意」という情報を伝えることができる
econ 以外のコントロール変数が増えても同様に考える
分析結果 summary() で得られた出力を単にコピペするのではなく、stargazer() を使って分析のタイトルや変数名を明示して見やすい表で示すこと

近年の表示方法

近年、社会科学の学術論文では重回帰分析結果を可視化する傾向が著しい
上記のような分析結果を示す従来の表は Appendix に入れ、可視化した表記法が好まれる　　
ここでは R パッケージを使った2 つの方法 (sjPlot, jtools) を紹介する

(1) sjPlot を使った分析結果の可視化

coef_sjplot <- sjPlot::plot_model(Model_3, 
                          show.values = T,       #推定値の係数を表示する
                          show.p = T,            #統計的有意性を*で表示する  
                          vline.color = "black", #推定値=0の縦線の色を黒に指定  
                          order.terms = c(1, 2),　
                          # 説明変数を表示する順番        
                          # model_3で cc は 1 番目の説明変数、econは 2 番目の説明変数と指定
                          axis.labels = c("地方政府の経済指標 (econ)", "市民共同体指標 (cc)"),
                          # order.termsと表示順番が逆なので注意!
                          axis.title = "推定値", title = "イタリアにおける地方政府のパフォーマンス",
                          digits = 3) # 推定値の表示を小数点3桁に指定

coef_sjplot

confinf関数を使って信頼区間の数値を表示できる

confint(Model_3, levels=0.95)

                 2.5 %    97.5 %
(Intercept)  1.3112706 5.1603219
cc           0.4330438 1.0745542
econ        -0.6485673 0.1418211

【解釈の仕方】：

ドット「●」は各変数の推定値 (= 係数：coefficient)
横線は95%信頼区間（推定値から ±1.96 x 2 標準誤算範囲）
有意水準 α = 0.05 (= 5%) (defaultの有意水準は 5%)
青の横線・・・縦の黒線に触れない →　統計的に有意（有意水準 5%）
赤の横線・・・縦の黒線に触れる →　統計的に有意でない（有意水準 5%）

【結果】：

市民共同体指標 (cc)・・・有意水準 5% で統計学的に有意
地方政府の経済指標 (econ)・・・有意水準 5% で統計学的に有意でない

標準偏回帰係数（β：ベータ）の可視化

type = "std" と指定することで β 値の結果を可視化できる

beta <- plot_model(Model_3, 
                          show.values = T,       #推定値の係数を表示する
                          show.p = T,            #統計的有意性を*で表示する
                          vline.color = "black", #推定値=0の縦線の色を黒に指定  
                          order.terms = c(1, 2), # 説明変数を表示する順番        
                                                 # model_3で cc は 1 番目の説明変数、econは 2 番目の説明変数と指定
                          axis.labels = c("地方政府の経済指標 (econ)", "市民共同体指標 (cc)"),
                          # order.termsと表示順番が逆なので注意!
                          axis.title = "推定値", title = "",
                          digits = 3,   # 推定値の表示を小数点3桁に指定
                          type = "std" ) # 標準偏回帰係数（β：ベータ）を確認

beta

上記で数値で示した β値が可視化されていることがわかる

・「cc が 1 標準偏差高いと gov_p が 1.193標準偏差増える」

・「econ が 1 標準偏差高いと gov_p が 0.326標準偏差減る」

(2) 近年の表示方法 (jtools)

jtools::plot_summs(Model_3)

回帰分析結果の表示法まとめ・最も重要な説明変数を全てのモデルに含める
・論文では sfPlot や jtools を使って可視化した「近年の表示方法」で結果を示す
・「従来の表示方法」は論文の Appendix に入れる

5.「決定係数」と「修正済み決定係数」

5.1 決定係数 (\(R^2\))

モデルの「当てはまりのよさ」を表す指標
回帰分析では、最小二乗法によって観測値と予測値のズレ (= 残差 residual) ができるだけ小さくなるよう直線を引く

その試みがうまくいっている程度を示すのが決定係数 \(R^2\)
決定係数 \(R^2\)は次の式で表すことができる：

\[R^2\ = \frac{モデルで説明された部分}{総変動（応答変数の変化）}\] \[ = 1 - \frac{残差（モデルで説明されない部分）}{総変動（応答変数の変化）}\]

つまり決定係数は観測された「総変動（応答変数のばらつき）」のうち何パーセントがモデルの「説明変数」によって説明できるかを示す数値
決定係数が\(R^2\)と呼ばれる理由：
回帰モデルによって予測された予測値 \(\hat{y_i}\) と現実の観測値 \(y_i\) との間の相関 (correlation) 係数を二乗すると \(R^2\) の値になるから
決定係数は 0 以上 1 以下の範囲の値をとる
値が 1 に近いほど回帰分析の当てはまりがよい
例）決定係数\(R^2 = 0.85\)
→観測された応答変数 \(y_i\) のばらつきの 85％が、予測値（\(\hat{y_i}\)）のばらつきによって説明される
決定係数は回帰分析の大まかな当てはまりのよさを示してくれる
決定係数がいくつ以上であればよい分析であるといえるような基準はない
社会科学の研究では決定係数が 0.1 程度のものもたくさんある

\(R^2\)が大きなモデルほど「いいモデル」とはいえない

決定係数には、説明変数の数が増えるほど大きな値をとるという性質がある（決して減少しない）
実際には当てはまりがそれほどよくなっていなくても、説明変数の数を増やすと決定係数が 1 に近づく

その理由

モデルに新たな説明変数を追加する
→ 今までのモデルで説明できなかった部分を、新たな変数が説明する
→ 新たなモデルによって説明できる部分が増加する
→ 決定係数 \(R^2\) 値が大きくなる
→ 理論的には、全ての目的変数をモデルに含めると- - -
→ 決定係数 \(R^2 = 1\) （応答変数の 100% が説明された!!）が可能
しかし、実際にはそのようなことはしない

その理由

モデルを複雑にするほど、現実をより良く説明できるが、扱いにくく、説明する上での「切れ味」が悪い
「結局、何がいいたいの？」という研究になりがち
→ 問題関心に関連した必要な範囲内で、できるだけ単純な (parsimony)モデルほど、理解しやすく活用しやすい

5.2 修正済み決定係数 (\(Adjusted-R^2\))

モデルに説明変数を追加するメリット

説明変数を追加する　→　モデルの説明力が増加する

モデルに説明変数を追加するデメリット

説明変数を追加する　→　モデルがより複雑になる
説明変数を追加する　→　\(R^2\) の値が自動的に大きくなる

対処法 → 修正済み決定係数 (\(Adjusted-R^2\))

説明変数の数を考慮に入れ、決定係数の値を割り引いて考えた指数
重回帰分析では修正済み決定係数 がよく使われる
修正済み決定係数 \(Adjusted-R^2\)は次の式で表すことができる：
(n：データ数、k: 説明変数の個数)

\[Adjusted-R^2\ = 1 - \frac{残差（モデルで説明されない部分）/(n-k-1)}{総変動（応答変数の変化）/(n-1)}\]

k (説明変数の個数) が分子にある
→ 説明変数の数が増える
→ \(Adjusted-R^2\)の値は小さくなる
\(R^2\) の場合と異なり、複数のモデルを比較して\(Adjusted-R^2\)が最大のモデルを選ぶ
前節に取り上げた「イタリア州政府のパフォーマンス」の分析結果を確認

**分析結果（イタリア州政府のパフォーマンス）**

	政府のパフォーマンス
	Model 1	Model 2	Model 3

cc	0.567^***		0.754^***
	(0.066)		(0.152)
econ		0.589^***	-0.253
		(0.120)	(0.187)
Constant	2.711^***	3.011^**	3.236^***
	(0.844)	(1.385)	(0.912)
N	20	20	20
R-squared	0.806	0.572	0.825
Adj. R-squared	0.796	0.549	0.805
Residual Std. Error	1.789 (df = 18)	2.659 (df = 18)	1.749 (df = 17)
F Statistic	74.967^*** (df = 1; 18)	24.097^*** (df = 1; 18)	40.126^*** (df = 2; 17)

p < .01; p < .05; p < .1

Model_1, Model_2, Model_3 いずれにおいても \(R^2\) よりも \(Adjusted-R^2\) の方が小さい
三つのモデルの\(Adjusted-R^2\) を比較すると、Model_3 (0.805) が最も大きく、モデルの当てはまりが良い

注意：説明変数を増やすと\(Adjusted-R^2\) の値は減ることもあり得る

5.3 適切なモデルの選び方

高い \(Adjusted-R^2\) が得られるようにモデルをいじるべきではない
論文のテーマを決め、モデルを決める際に重要なのは「理論に導かれて変数を選ぶ」こと
理論的に考えて「この変数は応答変数に影響を与えているに違いない」と思えば（例えその変数をモデルに加えることで \(Adjusted-R^2\) が低下しても）モデルに加えるべき
その変数を加えると \(Adjusted-R^2\) が上昇するとしても、理論的に考えて不適切だと思うならその変数をモデルに含めてはいけない
\(Adjusted-R^2\) がだいたいどの程度かという判断は分析対象によって異なる

★応答変数が多くの要因によって左右される場合:

例：応答変数が「企業の業績」- - - \(Adjusted-R^2\) は 0.1〜0.3 程度
→ 「企業の業績」に影響する要因は非常に多く、それらのデータを入手するのは困難
→　応答変数の 10% 〜 30%しか説明できないのは仕方がない

★説明変数に前年のデータ（ラグ変数）などが含まれる場合:

例：応答変数が「大学への出願者数」など - - - \(Adjusted-R^2\) は 0.9 程度
→ ここで「ラグ変数」とは「前年度の出願者数」のこと
→ 出願者数は毎年それほど大きな変動はない
→「今年度の出願者数」は「前年度の出願者数」でほとんど説明されてしまう
→　応答変数の 90%が説明されるのは当然

5.4 Model 1 の \(R^2\) を計算してみる　　

Model 1 は市民共同体指標 (cc)　がイタリアにおける地方政府のパフォーマンス (gov_p) に与える効果を表す

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

\(TSS\)	: `Total Sum of Squares`（総変動）	応答変数 \(gov_{p}\) の総変動
\(var[gov_{pi}]\)	: 応答変数 \(gov_p\) の不偏分散
\(ESS\)	: `Explained Sum of Squares`（説明された部分）	予測値と平均値の差
\(USS\)	: `Unexplained Sum of Squares`（残差 = \(e_i\)）	予測値と観測値の差

データを読み取る

putnam <- read.csv("data/putnam.csv")

必要は変数だけを選ぶ

df0 <- putnam  %>%
  dplyr::select(gov_p, cc)

head(df0)

TSS (Total Sum of Squares)

変数 \(Y\) の総変動

\[TSS = \sum_{i=1}^n(Y_i - \bar{Y})^2\]

変数 \(Y\) の個々の値 \(Y_i\) と変数 \(Y\) の平均値 \(\bar{Y}\) との距離（= 偏差）: \(Y_i - \bar{Y}\)
その偏差を二乗してから足して偏差平方和を作る: \(\sum_{i=1}^n (Y_i - \bar{Y})^2\)
ここで求めたいのは変数 \(gov_p\) の総変動だから

\[TSS = \sum_{i=1}^{20}(gov_{pi} - \bar{gov_p})^2\]

hensa <- with(df0, gov_p - mean(gov_p))  # gov_pの個々の値とgov_pの平均値との距離（偏差）を計算
hensa2 <- hensa^2 # 偏差を二乗する  
tss <- sum(hensa2) # 二乗した偏差を足す（偏差平方和）
tss

[1] 297.55

変数 \(Y_i\) の不偏分散

\[var[Y_i] = \frac{1}{n-1}\sum_{i=1}^n (Y_i - \bar{Y})^2 \\ = \frac{1}{n-1}TSS\]

ここで求めたいのは変数 \(gov_p\) の不偏分散だから

\[var[gov_{pi}] = \frac{1}{n-1}TSS \\ = \frac{1}{20-1}・297.55 = 15.66053 \]

hensa <- with(df0, gov_p - mean(gov_p))  # gov_pの個々の値とgov_pの平均値との距離（偏差）を計算
hensa2 <- hensa^2 # 偏差を二乗する  
tss <- sum(hensa2) # 二乗した偏差を足す（偏差平方和）
tss/(20-1) # 偏差平方和を 20-1 で割る

[1] 15.66053

不変分散は var() 関数を使っても計算できる

var(df0$gov_p)

[1] 15.66053

ESS (Explained Sum of Squares)

\(X\) を説明変数とする単回帰モデルから説明できる \(Y\) の変動部分: \((\hat{Y_i} - \bar{Y})\)

\[ESS= \sum_{i=1}^n (\hat{Y_i} - \bar{Y})^2\]

応答変数 \(Y\) の予測値 \(\hat{Y_i}\) と平均値 \(\bar{Y_i}\) の差を計算する: \(\hat{Y_i} - \bar{Y}\)
この差を二乗してから足して平方和を作る:　\(\sum_{i=1}^n (\hat{Y_i} - \bar{Y_i})^2\)
ここで求めたいのは応答変数 \(gov_p\) の ESS なので

\[ESS = \sum_{i=1}^{20} (\hat{gov_{pi}} - \bar{gov_p})^2\]

\(gov_p\) の予測値は \(\hat{gov_p} = 2.711 + 0.567cc_i\) なので、これを yhat1 と定義する

yhat1 <- 2.711 + 0.567 * df0$cc

応答変数 \(Y\) の予測値 \(\hat{Y_i}\) と平均値 \(\bar{Y_i}\) の差を二乗した値を yhat2 と定義する

yhat2 <- (yhat1 - mean(df0$gov_p))^2

yhat2 を足して平方和を作る　→　これが ESS

ess <- sum(yhat2)
ess

[1] 239.6864

USS (Unexplained Sum of Squares)

\(X\) を説明変数とする単回帰モデルから説明できない \(Y\) の変動部分: \((Y_i - \hat{Y})\)

\[USS= \sum_{i=1}^n (Y_i - \hat{Y})^2 = \sum_{i=1}^n\epsilon_i^2\]

応答変数 \(Y\) の観測値 \(Y_i\) と予測値 \(\hat{Y_i}\) の差を計算する: \(Y_i - \hat{Y}\)
この差を二乗してから足して平方和を作る: \(\sum_{i=1}^n(Y_i - \hat{Y_i})^2\)
ここで求めたいのは応答変数 \(gov_p\) の USS なので

\[USS = \sum_{i=1}^{20} (gov_{pi} - \hat{gov_{pi}})^2\]

\(gov_p\) の予測値は \(\hat{gov_{pi}} = 2.711 + 0.567cc_i\) なので、これを yhat1 と定義する

yhat1 <- 2.711 + 0.567 * df0$cc

応答変数 \(Y\) の観測値 \(Y_i\) の差を二乗した値を e1 と定義する

e1 <- (df0$gov_p - yhat1)^2

e1 を足して平方和を作る　→　これが USS

uss <- sum(e1)
uss

[1] 57.61108

決定係数 \(R^2\) を計算する

決定係数 \(R^2\)は次の式で表すことができる：

\[R^2\ = \frac{モデルで説明された部分: ESS}{総変動:TSS}\]　　

\[ = 1 - \frac{残差（モデルで説明されない部分）: US S}{総変動: TSS} \\ = 1 - \frac{57.61108}{297.55} = 0.80638\]

\(R\) を使ってダブルチェックしてみる

R_sq <- 1 - (uss/tss)
R_sq

[1] 0.8063819

Model_1 の回帰分析結果を表示してみる

summary(Model_1)


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

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

Multiple R-squared: 0.8064 が確認できる

5.5 ベン図を使った決定係数の解説

上では応答変数が gov_p、説明変数が cc の単回帰モデルにおける \(R^2\) を計算した
応答変数が gov_p、説明変数が cc と econ の重回帰モデルおける TSS, ESS, USS の関係をベン図を使って考えてみる

5.5.1 説明変数の間に相関がない場合

説明変数 cc は共変量だが、交絡因子ではない場合
説明変数の間に相関がない場合、新たな説明変数として econ を加えると、ESS が緑色の部分だけ増える
→　モデル全体で説明できる部分 gov_p の総変動 (TSS) が増える →　gov_p を説明する力が向上する
→　新たな説明変数 econ をモデルに加えるべき

応答変数が gov_p、説明変数が econ の単回帰モデル
\[\hat{gov_t} = \hat{β_0} + \hat{β_1}econ_i・・・(1)\]
応答変数が gov_p、説明変数が econ と cc の重回帰モデル

\[\hat{gov_t} = \hat{β_0} + \hat{β_1}econ_i + \hat{β_2}cc_i・・・(2)\]　　

cc と econ が無相関なら　→　式(1) と (2) の \(\hat{β_1}\) の値は一致する
式(1)で推定しても式(2) で推定しても、\(\hat{β_1}\) は偏りなく推定できる
→ econ をモデルに入れなくても、\(\hat{β_1}\) は \({β_1}\) の不偏推定量
説明変数の間に相関がない場合
→　式(2)の標準誤差の方が式(1)の標準誤差よりも小さくなる
→　より推定の精度が向上する

5.5.2 説明変数の間に相関がある場合

説明変数 cc は共変量であり、交絡因子である場合
市民共同体指標 (cc) が gov_p と econ の両方に影響を与える交絡因子
この場合の ESS を可視化して表すと次のようになる

この場合、図の赤い部分が「cc からの間接的な効果」（＝交絡）　
この部分を取り除いて分析しなければ、gov_p に対する econ の影響が正しく推定されない
cc をモデルに含めることで、この部分を取り除いた分析が可能になる
応答変数が gov_p、説明変数が econ の単回帰モデル
\[\hat{gov_t} = \hat{β_0} + \hat{β_1}econ_i・・・(1)\]
応答変数が gov_p、説明変数が econ と cc の重回帰モデル

\[\hat{gov_t} = \hat{β_0} + \hat{β_1}econ_i + \hat{β_2}cc_i・・・(2)\]　　

cc と econ が相関しているので　→　式(1) と (2) の \(\hat{β_1}\) の値は一致しない
交絡変数 (cc) を含めた式(2)のモデルを使って \(\hat{β_1}\) を推定する必要がある

5.5.3 「回帰モデルの標準誤差」と「回帰係数の標準誤差」

標準誤差とは、母数を推定する標本統計量の標準偏差のこと
母集団から標本を抽出するたびに異なる標本が得られる
標本統計量は、標本ごとに異なる
標準誤差はこの「ばらつき」の程度を示す数値
「得られた標本統計量が、母数から典型的にどれだけ離れているか」という統計的推測における不確実性を示す
標準誤差が小さいほど、推定の精度が高い

\[t = \frac{\hat{β_1} - β_1}{SE} = \frac{変数の偏回帰係数}{標準誤差}\]

→　標準誤差が小さいほど、\(t\) 値が大きい　→　\(p\) が小さい

回帰モデル全体の標準誤差: \(s_{reg}\)

まず回帰モデル全体の標準誤差 \(s_{reg}\) (standard error of the regression) を計算する

\[s_{reg}= \sqrt\frac{\sum(Y_i-\hat{Y_i})^2}{n-k}= \sqrt{\frac{USS}{n-k}}\]

\(n-k\) は自由度：\(k\) は推定すべき母数（パラメータ）の数：ここでは \(n-k = 20-2=18\)

\[s_{reg}= \sqrt{\frac{USS}{n-k}} = \sqrt{\frac{57.61108}{20-2}}= 1.789\]

summary(Model_1) の出力を確認

回帰係数 \(\beta_1\) の標準誤差: \(s.e.(\hat{\beta_1})\)

\[s.e.(\hat{\beta_1}) = \frac{s_{reg}}{\sqrt{\sum(X_i-\bar{X})^2}}\]

\(\sum(X_i-\bar{X})^2\) を計算する

st1 <- with(df0, cc - mean(cc))  # ccの個々の値とccの平均値との距離（偏差）を計算
st2 <- st1^2 # 偏差を二乗する  
st3 <- sum(st2) # 二乗した偏差を足す（偏差平方和）
st3

[1] 745.55

\[s.e.(\hat{\beta}) = \frac{s_{reg}}{\sqrt{\sum(X_i-\bar{X})^2}} = \frac{1.789}{\sqrt{745.55}} = \frac{1.789}{27.3}=0.0655\]

summary(Model_1) の出力を確認

次の式に Std. Error = 0.06552 と傾き \(\hat{β_1}= 0.56730\)、 \(β_1= 0\) を代入すると t 値が得られる

\[t 値 = \frac{\hat{β_1} - β_1}{SE} =\frac{偏回帰係数}{標準誤差}\\ = \frac{0.56730}{0.06552}= 8.658\]

6. 二乗項を含めたモデル

6.1 データの準備 (`hr96-21.csv`)

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

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

download.file(url = "http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/Data/hr96-21.csv",
              destfile = "data/hr96-21.csv")

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

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

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

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

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

読み取り方法 1

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

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

読み取り方法 2

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

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

読み取り方法 3

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

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

hr96_17.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"

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

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

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

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

hr96-21.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 に含まれる変数を確認する

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"

変数名	詳細
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	世襲元の政治家の氏名と関係

このデータフレーム hr の中から2005年の衆議院選挙データだけ取り出し、次の 3 つの変数 (voteshare, exppv, ldp) を使って重回帰分析を行う

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

hr12 <- hr %>% 
  dplyr::filter(year == 2012) %>% 
  mutate(exp = exp/1000000) %>% 
  mutate(exp_sq = exp^2) %>% 
  mutate(previous_sq = previous^2)

ggplot(hr12, aes(exp, voteshare)) +
  geom_point() +
  labs(x = "選挙費用", y = "得票率",
         title = "選挙費用と得票率の散布図")

ggplot(hr12, aes(exp_sq, voteshare)) +
  geom_point() +
  stat_smooth(method = lm) +
  labs(x = "選挙費用（二乗）", y = "得票率",
         title = "選挙費用と得票率の散布図")

model_1 <- lm(voteshare ~ exp, data = hr12)
model_2 <- lm(voteshare ~ exp + exp_sq, data = hr12)

stargazer(model_1, model_2, type = "html")


	Dependent variable:

	voteshare
	(1)	(2)

exp	2.687^***	4.725^***
	(0.089)	(0.226)

exp_sq		-0.131^***
		(0.013)

Constant	7.735^***	2.527^***
	(0.629)	(0.809)


Observations	1,280	1,280
R²	0.417	0.458
Adjusted R²	0.417	0.457
Residual Std. Error	13.046 (df = 1278)	12.592 (df = 1277)
F Statistic	915.539^*** (df = 1; 1278)	538.891^*** (df = 2; 1277)

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

ggplot(hr12, aes(previous_sq, voteshare)) +
  geom_point() +
  stat_smooth(method = lm)

model_3 <- lm(voteshare ~ previous, data = hr12)
model_4 <- lm(voteshare ~ previous + previous_sq, data = hr12)

stargazer(model_3, model_4, type = "html")


	Dependent variable:

	voteshare
	(1)	(2)

previous	4.714^***	8.653^***
	(0.176)	(0.404)

previous_sq		-0.534^***
		(0.050)

Constant	17.912^***	16.676^***
	(0.430)	(0.428)


Observations	1,294	1,294
R²	0.356	0.409
Adjusted R²	0.355	0.408
Residual Std. Error	13.735 (df = 1292)	13.165 (df = 1291)
F Statistic	713.873^*** (df = 1; 1292)	446.223^*** (df = 2; 1291)

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

7. 条件付き期待値としての回帰分析

回帰モデルは、説明変数 \(x\) がある値をとる時の応答変数 \(y\) の条件付き平均値（＝期待値）を計算している
観測値 \(Y_i\) は次の回帰式で表すことができる

\[E[Y_i|X_i] = \beta_0 + \beta_1X_i + \epsilon_i\]

\(Y_i\)	応答変数
\(X_i\)	説明変数
\(ε_i\)	誤差項 (`error term)`

X と Y 全ての値が一直線上に収まるとは限らない
→ 誤差的な変動を表すために誤差項 (\(\epsilon_i\))を想定する
\(\epsilon_i\) の期待値はゼロと想定
ここでは次の 3 つの条件を仮定している

\(E[ε]=0\)

誤差項 (\(\epsilon_i\))の期待値がゼロである
誤差であればゼロを中心にデータがばらつくはずと想定

\(V[ε_i]=\sigma^2\)

\(\sigma_i\) ではなくて \(\sigma\) であることに注意
全ての個体に対して、ばらつきの大きさは同じという意味

\(cov[x,ε]=0\)

説明変数 \(X\) と誤差項 \(ε\) の共分散はゼロ（＝独立である）
説明変数 X は誤差項 \(ε\) と相関がない
i 番目の個体と j 番目の個体の攪乱項に関しては関係が無い

3 つの仮定を置いた上で、\(\beta_0\) と \(\beta_1\) という未知のパラメータをデータから推定する

事例を使った解説

y <- c(1,2,3,4,5,6,7,8,9)
x <- c(1,1,1,2,2,2,3,3,3)
xy <- tibble(x,y)
xy

# A tibble: 9 × 2
      x     y
  <dbl> <dbl>
1     1     1
2     1     2
3     1     3
4     2     4
5     2     5
6     2     6
7     3     7
8     3     8
9     3     9

散布図を描いてみる

xy %>% 
  ggplot(aes(x, y)) + 
  geom_point()

回帰分析してみる

lm(y ~ x, data = xy)


Call:
lm(formula = y ~ x, data = xy)

Coefficients:
(Intercept)            x  
         -1            3

\[\hat{y}= -1 + 3x_i\]

\(x\) の値が 1 の時の \(y\) の予測値 \(\hat{y}\) は

-1 + 3 * 1

[1] 2

これは \(x=1\) という条件時の \(y\) の期待値が 2 という意味: \(E[Y_i|X_i = 1]=2\)

\(x\) の値が 2 の時の \(y\) の予測値 \(\hat{y}\) は

-1 + 3 * 2

[1] 5

これは \(x=2\) という条件時の \(y\) の期待値が 5 という意味: \(E[Y_i|X_i = 2]=5\)

\(x\) の値が 3 の時の \(y\) の予測値 \(\hat{y}\) は

-1 + 3 * 3

[1] 8

これは \(x=3\) という条件時の \(y\) の期待値が 8 という意味: \(E[Y_i|X_i = 3]=8\)
回帰直線を加えた散布図を描いてみる

xy %>% 
  ggplot(aes(x, y)) + 
  geom_point() +
  stat_smooth(method = lm, se = FALSE)

\(x=1, 2, 3\) それぞれに対応する \(y\) の予測値がそれぞれ \(y = 2, 5, 8\)

ポイント 回帰モデルは、説明変数 \(x\) がある値をとる時の応答変数 \(y\) の条件付き平均値（＝期待値）を計算している

8. アンスコムの例 (Anscombe’s qurtet)

統計学者のフランク・アンスコムが1973年に紹介した回帰分析において注意すべき例
回帰分析において、散布図はそれぞれ異なるのに回帰直線やその他の統計量が同じになってしまう現象のこと

教訓：

回帰分析を実行する前に散布図を確認し傾向を把握することは大切
←　その理由：外れ値が統計量に大きな影響を与えるから

ここで必要な R package

library(fBasics)
library(ggplot2)
library(grid)
library(gridExtra)
library(datasets)

datasets パッケージに内包されているデータを読み込み df_qurtet と名前を付ける

df_qurtet <- datasets::anscombe

df_qurtet の中身を表示させる

DT::datatable(df_qurtet)

df_qurtet の記述統計を表示させる

stargazer::stargazer(df_qurtet,
  type = "html")


Statistic	N	Mean	St. Dev.	Min	Max

x1	11	9.000	3.317	4	14
x2	11	9.000	3.317	4	14
x3	11	9.000	3.317	4	14
x4	11	9.000	3.317	8	19
y1	11	7.501	2.032	4.260	10.840
y2	11	7.501	2.032	3.100	9.260
y3	11	7.500	2.030	5.390	12.740
y4	11	7.501	2.031	5.250	12.500

x1, x2, x3, x4 の統計量（観測値、平均値、標準偏差、最小値、最大値）はほとんど同じ
y1, y2, y3, y4 の統計量（観測値、平均値、標準偏差、最小値、最大値）はほとんど同じ

回帰分析

model_1 <- lm(y1 ~ x1, data = df_qurtet)
model_2 <- lm(y2 ~ x2, data = df_qurtet)
model_3 <- lm(y3 ~ x3, data = df_qurtet)
model_4 <- lm(y4 ~ x4, data = df_qurtet)

stargazer::stargazer(model_1, model_2, model_3, model_4,
  type = "html")


	Dependent variable:

	y1	y2	y3	y4
	(1)	(2)	(3)	(4)

x1	0.500^***
	(0.118)

x2		0.500^***
		(0.118)

x3			0.500^***
			(0.118)

x4				0.500^***
				(0.118)

Constant	3.000^**	3.001^**	3.002^**	3.002^**
	(1.125)	(1.125)	(1.124)	(1.124)


Observations	11	11	11	11
R²	0.667	0.666	0.666	0.667
Adjusted R²	0.629	0.629	0.629	0.630
Residual Std. Error (df = 9)	1.237	1.237	1.236	1.236
F Statistic (df = 1; 9)	17.990^***	17.966^***	17.972^***	18.003^***

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

Model_1 から Model_4 それぞれの回帰式は次のとおり

\[y_1 = 0.5x_1 + 3\]

\[y_2 = 0.5x_2 + 3\]

\[y_3 = 0.5x_3 + 3\]

\[y_4 = 0.5x_4 + 3\]

x1とy1, x2とy2, x3とy3, x4とy4それぞれの散布図を描いてみる

M1 <- df_qurtet |> 
  ggplot() +
  geom_point(aes(x1, y1), color = "darkorange", size = 1.5) +
  scale_x_continuous(breaks = seq(0,20,2)) +
  scale_y_continuous(breaks = seq(0,12,2)) +
  expand_limits(x = 0, y = 0) +
  labs(x = "x1", y = "y1",
       title = "Model_1" ) +
  theme_bw() +
  geom_abline(intercept = 3, slope = 0.5, color = "blue")

M2 <- df_qurtet |> 
  ggplot() +
  geom_point(aes(x2, y2), color = "darkorange", size = 1.5) +
  scale_x_continuous(breaks = seq(0,20,2)) +
  scale_y_continuous(breaks = seq(0,12,2)) +
  expand_limits(x = 0, y = 0) +
  labs(x = "x2", y = "y2",
       title = "Model_2" ) +
  theme_bw() +
  geom_abline(intercept = 3, slope = 0.5, color = "blue")

M3 <- df_qurtet |> 
  ggplot() +
  geom_point(aes(x3, y3), color = "darkorange", size = 1.5) +
  scale_x_continuous(breaks = seq(0,20,2)) +
  scale_y_continuous(breaks = seq(0,12,2)) +
  expand_limits(x = 0, y = 0) +
  labs(x = "x3", y = "y3",
       title = "Model_3" ) +
  theme_bw() +
  geom_abline(intercept = 3, slope = 0.5, color = "blue")

M4 <- df_qurtet |> 
  ggplot() +
  geom_point(aes(x4, y4), color = "darkorange", size = 1.5) +
  scale_x_continuous(breaks = seq(0,20,2)) +
  scale_y_continuous(breaks = seq(0,12,2)) +
  expand_limits(x = 0, y = 0) +
  labs(x = "x4", y = "y4",
       title = "Model_4" ) +
  theme_bw() +
  geom_abline(intercept = 3, slope = 0.5, color = "blue")

library(patchwork)

(M1 + M2 + M3 + M4)

Model_1

\(x_1\) と \(y_1\) は線形関係　→　回帰分析ができる

Model_2

\(x_2\) と \(y_2\) は線形関係ではない　→　回帰分析は不適切

Model_3

\(x_3\) と \(y_3\) には線形関係だが外れ値がある
→　ロバスト回帰分析 (robust regression) をする

Model_4

\(x_4\) と \(y_4\) には線形関係ではないかもしれない
→　右上の一つの外れ値をチェックすべき

9. Excercise

Excercise 8.1

smile_nagoya.csv は2015年に実施された政令指定都市の市議会選挙（名古屋市）の選挙結果データ
smile は選挙公報に掲載された候補者の顔写真（ポスターと同じ写真）を（株）オムロンの okao vision を使って測定した「笑顔度」
候補者の写真と「笑顔度」の例

このデータには次の16の変数が含まれている

変数名	詳細
1. district	選挙区名
2. eligible	有権者数
3. turnout	投票率
4. dm	選挙区定数
5. nocand	立候補者数
6. rank	ランク
7. wl	当選: win、落選: lose
8. status	現職: incumbent、新人: challenger、元職: former_inc
9. votes	得票数
10. voteshare	得票率 (%)
11. male	男性: male、女性: female
12. age	年齢
13. smile	笑顔度 (0-100)
14. inc	現職: inc、それ以外: not_inc
15. previous	当選回数（2015年選挙結果を含む）
16. party	所属政党

次の問題にこたえなさい

Q1: party 変数を使って ldp という名称の「自民党ダミー（自民党候補者なら 1、それ以外は 0）」を作成し、データフレームに加えなさい（「ダミー変数」の作り方は「3．データの可視化1（基礎）」を参照）

Q2: stargazer パッケージを使って smile_nagoya.csv の記述統計を示しなさい

Q3: voteshare を応答変数、smile を説明変数とした散布図を描きなさい（95％信頼区間を示す回帰直線も含む）その際、geom_jitter() 関数を使ってデータを適度に散らすこと

Q4: この散布図から smile と voteshare の間にはどのような関係があると言えるか述べなさい

Q5: 立候補者の笑顔度 (smile) から得票率 (voteshare) を予測したい

次の 13 の変数の中でコントロール変数として含めるべき変数を全て選び、「含めるべき理由」と「含めるべきでない理由」をそれぞれ述べなさい

変数名	詳細	統制変数に含めるべき？その理由
1. district	選挙区名
2. eligible	有権者数
3. turnout	投票率
4. dm	選挙区定数
5. nocand	立候補者数
6. rank	ランク
7. wl	当選: win、落選: lose
8. votes	得票数
9. male	男性: male、女性: female
10. age	年齢
11. inc	現職: inc、それ以外: not_inc
12, previous	当選回数（2015年選挙結果を含む）
13. ldp	自民党ダミー（自民党候補者なら 1、それ以外は 0）

Q6: Q5で皆さんが選んだコントロール変数を一つだけ含めて、voteshare を応答変数、smile を説明変数とした重回帰分析を実行し model_2 と名前を付けなさい

model_2 に皆さんが選んだコントロール変数を一つずつ含めて model_3, model_4,… と複数の回帰分析を実行しなさい
stargazer()関数を使って、実行した回帰分析結果を表示し、解釈しなさい
解釈する際、必要に応じてsummary()関数、もしくは tidy()関数を使うこと

Q7: sjPlot もしくは jtools パッケージを使って、最も適切と思われるモデルを一つ選び、その分析結果を可視化しなさい

Q8: 残差プロットによる診断によって、誤差の独立性の診断しなさい

Q9: 正規QQプロットによる診断によって、誤差が正規分布に従うことを確認しなさい

Q10: VIF による診断によって、多重共線性の有無を確認しなさい

Excercise 8.2

・hr96-21.csv は1996年から2021年に実施された総選挙結果である　　

このデータセットには 23 個の変数が入っている

変数名	詳細
1. year	選挙年 (1996-2021)
2. pref	都道府県名
3. ku	小選挙区名
4. kun	小選挙区
5. rank	当選順位
6. nocand	立候補者数
7. seito	候補者の所属政党
8. j_name	候補者の氏名（日本語）
9. name	候補者の氏名（ローマ字）
10. previous	当選回数
11. gender	立候補者の性別: “male”, “female”
12. age	立候補者の年齢
13. wl	選挙の当落: 1 = 小選挙区当選、2 = 復活当選、0 = 落選
14. wlsmd	選挙の当落: 1 = 当選（小選挙区）、0 = 落選（小選挙区）
15. exp	立候補者が使った選挙費用（総務省届け出）
16. status	候補者のステータス: 0 = 非現職、1 現職、2 = 元職
17. vote	得票数
18. voteshare	得票率 (%)
19. eligible	小選挙区の有権者数
20. turnout	小選挙区の投票率 (%)
21. seshu_dummy	世襲候補者ダミー: 1 = 世襲、0 = 非世襲（地盤世襲 or 非世襲）
22. jiban_seshu	地盤の受け継ぎ元の政治家の氏名と関係
23. nojiban_seshu	世襲元の政治家の氏名と関係

hr96-21.csv から 2009年のデータだけを抜き出し、次の問題にこたえなさい

Q1: seito 変数を使って dpj という名称の「民主党ダミー（民主党候補者なら 1、それ以外は 0）」を作成し、データフレームに加えなさい

Q2: exp の単位は「円」だが、この単位を「百万円」に変換し、変数名を expm と変更しデータフレームに加えなさい　

Q3: status 変数を使って inc という名称の「現職ダミー（現職なら 1、それ以外は 0）」を作成し、データフレームに加えなさい

Q4: stargazer パッケージを使ってデータフレームの記述統計を示しなさい

Q5: voteshare を応答変数、expm を説明変数とした散布図を描きなさい（95％信頼区間を示す回帰直線も含む）

Q6: この散布図から 2 つの変数の間にはどのような関係があると言えるか述べなさい　　

Q7: 立候補者が使う選挙費用 (expm) から得票率 (voteshare) を予測したい

次の 14 の変数の中でコントロール変数として含めるべき変数を全て選びなさい
変数として「含めるべき理由」と「含めるべきでない理由」をそれぞれ述べなさい

変数名	詳細	統制変数に含めるべき？その理由
1. rank	当選順位
2. nocand	立候補者数
3. j_name	候補者の氏名（日本語）
4. previous	当選回数
5. gender	立候補者の性別: “male”, “female”
6. age	立候補者の年齢
7. wl	選挙の当落: 1 = 小選挙区当選、2 = 復活当選、0 = 落選
8. wlsmd	選挙の当落: 1 = 当選（小選挙区）、0 = 落選（小選挙区）
9. vote	得票数
10. eligible	小選挙区の有権者数
11. turnout	小選挙区の投票率 (%)
12. seshu_dummy	世襲候補者ダミー: 1 = 世襲、0 = 非世襲（地盤世襲 or 非世襲）
13. dpj	民主党ダミー（民主党候補者なら 1、それ以外は 0）
14. inc	現職ダミー: 0 = 非現職、1 現職（小選挙区当選者）

Q8: Q7で皆さんが選んだコントロール変数を一つだけ含めて、voteshare を応答変数、expm を説明変数とした重回帰分析を実行し model_2 と名前を付けなさい

model_2 に皆さんが選んだコントロール変数を一つずつ含めて model_3, model_4,… と複数の回帰分析を実行しなさい
stargazer()関数を使って、実行した回帰分析結果を表示し、解釈しなさい
解釈する際、必要に応じてsummary()関数、もしくは tidy()関数を使うこと

Q9: sjPlot もしくは jtools パッケージを使って、最も適切と思われるモデルを一つ選び、その分析結果を可視化しなさい
Q10: 残差プロットによる診断によって、誤差の独立性の診断しなさい
Q11: 正規QQプロットによる診断によって、誤差が正規分布に従うことを確認しなさい
Q12: VIF による診断によって、多重共線性の有無を確認しなさい

参考文献

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

ロバート・パットナム（河田潤ー訳）『哲学する民主主義』NTT 出版、2001年.

森田果『実証分析入門』日本評論社、2014年.

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

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

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

高橋将宣『統計的因果推論の理論と実装』、共立出版、2022年

浅野正彦, 矢内勇生.『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, 2021

14. 重回帰分析 1 (単回帰と重回帰）

Masahiko Asano

2024-08-22

1. 回帰分析

回帰分析の目的：

1. 現象の「説明」

2. 未来の「予測」

3. x の値に人為的に「介入」 したときの効果

4. 制御 ：y を特定の範囲に収めるためには x の値をどうするか

1.1 単回帰分析 (Simple Linear Regression)

母集団と標本

最小二乗法 (OLS: Ordinary Least Squares)

1.2 手計算で単回帰式を計算してみる

1.3 回帰係数に関する仮説

2. 実証分析のプロセス

2.1 パズルを見つける

2.2 理論の提示

2.4 仮説の提示

3. 回帰分析による仮説検証

3.1 データの読み取り (putnam.csv)

3.2 記述統計量を示す

3.3 散布図を描く

Macユーザーのみ

Windows ユーザーは次のコマンドを入力する

3.4 単回帰分析モデル

3.4.1 回帰分析の結果の表示方法

3.4.2 結果の解釈

手計算で単回帰式を計算してみる（確認）

最小二乗法による予測値 \(\hat{Y}\) の式

傾き \(β\) は次の式で表せる

切片 \(α\) は次の式で表せる

\[\hat{Y}=2.71+0.567x\]

4. 重回帰モデル

重回帰分析の流れ

4.1 散布図

4.2 二つの単回帰分析

4.3 単回帰＋単回帰＝重回帰？

4.4 重回帰分析のメカニズム

重回帰で行っていること = 二段階推定

重回帰分析における cc の係数 (0.754) について考える

Step 1:

Step 2:

通常、このような二段階推定はせずに重回帰分析を実行する

4.5 モデル全体の有意性検定（F 検定）

4.6 偏回帰係数の検定（Model_3）

4.7 修正済み決定係数を確認

4.8 標準偏回帰係数（β：ベータ）を確認

4.9 重回帰分析結果の表示方法

5.「決定係数」と「修正済み決定係数」

5.1 決定係数 (\(R^2\))

5.2 修正済み決定係数 (\(Adjusted-R^2\))

5.3 適切なモデルの選び方

5.4 Model 1 の \(R^2\) を計算してみる

TSS (Total Sum of Squares)

変数 \(Y_i\) の不偏分散

ESS (Explained Sum of Squares)

USS (Unexplained Sum of Squares)

決定係数 \(R^2\) を計算する

5.5 ベン図を使った決定係数の解説

5.5.1 説明変数の間に相関がない場合

5.5.2 説明変数の間に相関がある場合

5.5.3 「回帰モデルの標準誤差」と「回帰係数の標準誤差」

回帰モデル全体の標準誤差: \(s_{reg}\)

回帰係数 \(\beta_1\) の標準誤差: \(s.e.(\hat{\beta_1})\)

6. 二乗項を含めたモデル

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

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

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

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

読み取り方法 1

読み取り方法 2

読み取り方法 3

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

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

7. 条件付き期待値としての回帰分析

\(E[ε]=0\)

\(V[ε_i]=\sigma^2\)

\(cov[x,ε]=0\)

3 つの仮定を置いた上で、\(\beta_0\) と \(\beta_1\) という未知のパラメータをデータから推定する

事例を使った解説

3. x の値に人為的に「介入」したときの効果

4. 制御：y を特定の範囲に収めるためには x の値をどうするか

最小二乗法 (`OLS: Ordinary Least Squares`)

3.1 データの読み取り (`putnam.csv`)

重回帰分析における `cc の係数 (0.754)` について考える

4.6 偏回帰係数の検定（`Model_3`）

4.8 標準偏回帰係数（`β`：ベータ）を確認　　

5.4 Model 1 の \(R^2\) を計算してみる　　

6.1 データの準備 (`hr96-21.csv`)