このセクションで使う 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\)）	予測値と観測値の差

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

Masahiko Asano

2025-01-30

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\) を計算してみる

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\) を計算してみる　　