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)
影響を与える変数 (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
)estimate
) するY
と X
の間の関係 \(Y = α + βX + ε\)\[\hat{Y} = \hat{α} + \hat{β}x\]
x
に対応する予測値\(\hat{Y}\) (predicted value
)
を計算できる\[\hat{ε} = Y - \hat{Y}\]
\[Y = \hat{α} + \hat{β}X + \hat{ε}\]
residual
or predictin error
)
と呼ばれるerror term
の推定値
(estimate
) なので、ハットを付けるOLS: Ordinary Least Squares
)\[\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{α} + \hat{β}x \]
\[\hat{Y} = \hat{α} + \hat{β}\bar{X}\]
\[\hat{Y} = (\bar{Y} - \hat{β}\bar{X}) + \hat{β}\bar{X} = \bar{Y}\]
\[\bar{Y} - \hat{β}\bar{X} = \hat{α}\]
\[(\bar{X}, \bar{Y})\]
を通る
手計算で単回帰式を計算
ここでは次のデータを使って、手計算で単回帰式を計算してみる
動画による解説はこちら
単回帰式の求め方をわかりやすく解説しています
この動画で使っている csv.file
はこちらからダウンロードできます
\[H_0:β = 0\]
\[H_A: β ≠ 0\]
① パズルを見つける
② 理論の提示(パズルのメカニズムの説明)
③
仮説(検証可能な仮説「もしこの理論が正しければ・・・のはずである」)
④ 仮説の検定
命題が「真」であることを直接的に証明することはできない
⇒ 否定されるべき仮説(=帰無仮説)を立てる
(パソコンの統計ソフトは帰無仮説を自動的に設定)
⇒ 帰無仮説を棄却する
⇒ 理論の正しさを間接的に示す
イタリアにおける地方政府のパフォーマンスを事例として考えてみる
リサーチクエスチョン:
「イタリア地方政府のパフォーマンスに著しい違いがあるのはなぜか?」
データ:
region
: イタリアの地方政府の略称gov_p
: 政府のパフォーマンス理論:社会関係資本が政府のパフォーマンスを高める
出典:Robert Putnam, (1994) Making Democracy Work: Civic Traditions
in Modern Italy,
Princeton, NJ: Princeton University Press,
(河田潤ー訳『哲学する民主主義:伝統と改革の市民的構造』NTT 出版、2001
年)
理論の概略
地方政府のパフォーマンスの違いは、その地域の社会関係資本の蓄積の度合いによって説明できる
「社会関係資本」(Social Capital
)
とは個人の結びつきのことで、互恵性の社会ネットワークや規範
「社会関係資本」は見知らぬ相手と協力関係を構築する一助となる
→
社会関係資本の蓄積が高い地域では、互いに信頼し協力しあうので、政府のパフォーマンスを高める
観察可能な応答変数と説明変数を示し、検証可能な仮説「もしこの理論が正しければ ・・・ のはずである」を明示する
【応答変数】
【説明変数】
Civic Community Index
(市民共同体指標)Clientelism
(政治的恩顧主義)の度合仮説:
putnam.csv
)RProject
フォルダー内の
data
フォルダに保存するputnam
と名前を付ける[1] "region" "gov_p" "cc" "econ" "location"
データ:
変数の種類 | 変数名 | 詳細 |
---|---|---|
region |
イタリア州政府の略称 | |
location |
イタリア北部・南部(north &
south ) |
|
応答変数 | gov_p |
政府のパフォーマンス |
説明変数 | cc |
Civic Community Index (市民共同体指標) |
説明変数 | econ |
地方政府の経済指標(大きい程、経済が良好) |
putnam
に含まれている変数の記述統計量を示す:summary()
を使う方法 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()
を使うと、記述統計を見やすく表示できる
パットナムデータの記述統計
======================================
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
---------------------------------------
<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 |
gov_p
と cc
の散布図を描いてみる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
を求める
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
) があるcc
を説明変数、gov_p
を応答変数とした回帰分析が必要gov_p
を応答変数、cc
を説明変数とした単回帰分析モデルを Model_1
と名前をつけて実行する回帰分析の結果の表示方法:
summary()関数
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()関数
# 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"}
と指定
Dependent variable: | |
gov_p | |
cc | 0.567*** |
(0.066) | |
Constant | 2.711*** |
(0.844) | |
Observations | 20 |
R2 | 0.806 |
Adjusted R2 | 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 |
\[\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")
応答変数 | |
政府のパフォーマンス | |
cc | 0.57*** |
(0.07) | |
Constant | 2.71*** |
(0.84) | |
Observations | 20 |
Adjusted R2 | 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”
を付け加えるtidy()
を使うと簡素な分析結果を表示できる# 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
) の値が 1
ポイントだけ異なる地方を比べると、cc
の値が大きい地方が、平均して
0.567ポイントだけ政府のパフォーマンス (gov_p
) が高い帰無仮説 \(H_0: β_1 = 0\)
対立仮説 \(H_a: β_1 ≠ 0\)
帰無仮説の意味
→ 母集団における cc
の傾き(\(β_1\)) は 0
→ つまり「cc
は gov_p
に影響を与えない」
対立仮説の意味
→ 母集団における c
c の傾き(\(β_1\)) は 0 ではない
→ つまり「cc
が gov_p
に与える影響を与えている」
次の式を使って t 値
を求め仮説検定を行う
\[t = \frac{\hat{β_1} - β_1}{SE}\]
Std. Error
)
= 推定値(係数)の標準偏差のこと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\)値7.81e-08 = 0.000000781
(ゼロが 8
つ)で、0.001
より小さい【解説】
cc
の gov_p
に対する影響はない)とき、単純無作為抽出で同じ大きさの標本を選ぶと、0.000000781
の確率で t 値の絶対値が 8.656 より大きいデータが得られるcc
の gov_p
に対する影響はない)とき、私たちが扱っているようなデータが出現する確率は
0.000000781
(=ほとんどない)R-squared
の値が 0.8064
→ 地方政府のパフォーマンス (gov_p
) の分散の約 81%
が市民共同体指標 (cc
) の分散で説明できた
→ だからといって、市民共同体指標が地方政府のパフォーマンスを引き起こしているとは限らない
R-squared = 0.8064
が意味すること
→このモデルがどれだけ上手に応答変数の変化を「説明」しているか
=R-squared
=「そのモデルがどれだけ正確か」についての指標
単回帰分析の限界・・・一つの説明変数以外の要因が含まれていない →
これは問題
ほとんどの現象は一つ以上の要因から影響をうけているはず
セレクションバイアスがある可能性が高い
解決策 →
重回帰分析
重回帰分析は他の要因の影響をコントロールした分析を可能にしてくれる
→ ある程度、セレクションバイアスの影響を回避できる
→ 重回帰分析は最も基本的なセレクションバイアスの削減方法
手計算で単回帰式を計算
xlsx.file
と模範解答はこちらからダウンロードできます\[\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 = \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\]
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 |
地方政府の経済指標 |
最小二乗法
)F 検定
)Adjusted R-squared
)を確認β
:ベータ)を確認この手続きに従って、イタリアにおける地方政府のパフォーマンスに関するデータを使って重回帰分析を行う
変数の種類 | 変数名 | 詳細 |
---|---|---|
応答変数 | gov_p |
政府のパフォーマンス |
説明変数 | cc |
Civic Community Index (市民共同体指標) |
コントロール変数 | econ |
地方政府の経済指標: 大きい程経済が良好 |
gov_p
と cc
gov_p
と econ
パットナムデータの記述統計
======================================
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
--------------------------------------
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)
gov_p
を応答変数、cc
を応答変数とした単回帰分析モデルを Model_1
と名前をつけるgov_p
を応答変数、econ
を応答変数とした単回帰分析モデルを Model_2
と名前をつける分析結果の表示法: 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
の係数にアステリスクが三つ (***
)
付いている分析結果 (Model_2
)
econ
が 1 ポイント高いと gov_p
が 0.589
ポイント高い」econ
の係数にアステリスクが三つ (***
)
付いている
ポイント:cc
も econ
も単回帰分析では統計的にそれぞれ有意
参考:
t 値
や
p 値
を確認したければ、summary()
関数や
tidy()
関数を使う# 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
# 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
gov_p
軸、cc
軸、econ
軸の三次元分析になるcc
と econ
が
gov_p
に与える影響を「同時に」確認することができ「変数をコントロール」できるgov_p
を応答変数、cc
を説明変数、econ
をコントロール変数とした重回帰分析モデルを
Model_3 と名前をつけて、結果を表示する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 の重回帰式を求める(最小二乗法)
Model_3
の回帰式は次の様になる\[\hat{gov_p}\ = 3.236 + 0.754cc - 0.253econ\]
\[\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\]
cc
を econ
に回帰する(Intercept) econ
-0.2984271 1.1173551
ESS
が「econ
によって説明された
cc
の部分 (= 1.117)」econ の係数 (1.117)
が得られた\[\hat{cc} = -0.298 + 1.117econ\]
econ
によって説明されれなかった cc
の残り部分」が残差:
e_cc
e_cc
を使って gov_p
を説明する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
に回帰する(Intercept) e_cc
9.150000 0.753799
e_cc
の係数 (=0.754)
が、重回帰分析における cc の係数 (= 0.754)
econ
をコントロールしたときの cc
の
gov_p
に対する影響力(Intercept) cc econ
3.2357962 0.7537990 -0.2533731
2.5 % 97.5 %
(Intercept) 1.3112706 5.1603219
cc 0.4330438 1.0745542
econ -0.6485673 0.1418211
cc
の影響力を推定できた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
」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 検定は「統計的に有意でない」変数をモデルから外すべきかどうかの判断に有益
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
Adjusted R-squared
(修正済み \(R^2\))が使われるAdjusted R-squared: 0.8046
gov_p
) の分散の80.46%
は
cc
と econ
の分散によって説明されるgov_p
の分散の約20%が説明されないまま (unexplained
)
残るという意味β
:ベータ)を確認 通常、説明変数の単位は異なるため(例えば、一つは cm
でもう一つは
kg
)、変数の係数値を比較するだけでは、どちらの変数が応答変数により大きな影響を与えているか比較できない
β値
(Standardized Beta Coefficient
)とは、説明変数の値を1標準偏差増やした時に、応答変数の標準偏差がどれだけ変化するかを示す
β値
を比較することで、複数の説明変数のうち、応答変数に対してどちらがより大きな影響を与えているのかわかる
cc econ
1.1932019 -0.3255233
cc
が 1 標準偏差高いと gov_p
が
1.19標準偏差増える」econ
が 1 標準偏差高いと gov_p
が
0.33標準偏差減る」従来の表示方法
政府のパフォーマンス | |||
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 |
注意:
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
に入れ、可視化した表記法が好まれる 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関数
を使って信頼区間の数値を表示できる 2.5 % 97.5 %
(Intercept) 1.3112706 5.1603219
cc 0.4330438 1.0745542
econ -0.6485673 0.1418211
【解釈の仕方】:
coefficient
)±1.96 x 2
標準誤算範囲)α = 0.05 (= 5%)
(defaultの有意水準は 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
)
sfPlot
や jtools
を使って可視化した「近年の表示方法」で結果を示すAppendix
に入れる
residual
)
ができるだけ小さくなるよう直線を引く\[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\)が大きなモデルほど「いいモデル」とはいえない
その理由
その理由
parsimony
)モデルほど、理解しやすく活用しやすいモデルに説明変数を追加するメリット
モデルに説明変数を追加するデメリット
対処法 → 修正済み決定係数 (\(Adjusted-R^2\))
n
:データ数、k
: 説明変数の個数)\[Adjusted-R^2\ = 1 - \frac{残差(モデルで説明されない部分)/(n-k-1)}{総変動(応答変数の変化)/(n-1)}\]
k
(説明変数の個数) が分子にある政府のパフォーマンス | |||
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\) の方が小さいModel_3
(0.805
)
が最も大きく、モデルの当てはまりが良い★応答変数が多くの要因によって左右される場合:
★説明変数に前年のデータ(ラグ変数)などが含まれる場合:
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\)) |
予測値と観測値の差 |