Rで学ぶ計量経済学と機械学習 4
回帰分析
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2023-10-07
1 本講義の目的
- Rを用いた回帰分析の方法を学ぶ。
- 要約統計量や回帰分析の結果の実践的な出力方法を学ぶ
2 パッケージの読み込み
まず、{tidyverse}や{AER}や{stargazer}など、今回用いるパッケージのインストールと読み込みを行う。新しい関数については、適宜説明していく。
3 データの用意
以下ではRで回帰分析を行っていくが、今回は{AER}パッケージ（Applied Econometrics with Rパッケージ）に収録されているCPS1985というデータセットを例に使用する。
今回使用するCPS1985データセットは、1985年アメリカの時給や教育年数などの変数が入ったものである。
data()関数でデータセットの読み込みを行う。
## wage education experience age ethnicity region gender occupation
## 1 5.10 8 21 35 hispanic other female worker
## 1100 4.95 9 42 57 cauc other female worker
## 2 6.67 12 1 19 cauc other male worker
## 3 4.00 12 4 22 cauc other male worker
## 4 7.50 12 17 35 cauc other male worker
## 5 13.07 13 9 28 cauc other male worker
## sector union married
## 1 manufacturing no yes
## 1100 manufacturing no yes
## 2 manufacturing no no
## 3 other no no
## 4 other no yes
## 5 other yes no4 記述統計の表示
回帰分析を行う前に、ほとんどの場合、分析に用いる統計の要約統計量（summary statistics)を載せる。
要約統計量をチェックするだけならば、summary()やskimr::skim(df)を使えばよい。
## wage education experience age
## Min. : 1.000 Min. : 2.00 Min. : 0.00 Min. :18.00
## 1st Qu.: 5.250 1st Qu.:12.00 1st Qu.: 8.00 1st Qu.:28.00
## Median : 7.780 Median :12.00 Median :15.00 Median :35.00
## Mean : 9.024 Mean :13.02 Mean :17.82 Mean :36.83
## 3rd Qu.:11.250 3rd Qu.:15.00 3rd Qu.:26.00 3rd Qu.:44.00
## Max. :44.500 Max. :18.00 Max. :55.00 Max. :64.00
## ethnicity region gender occupation sector
## cauc :440 south:156 male :289 worker :156 manufacturing: 99
## hispanic: 27 other:378 female:245 technical :105 construction : 24
## other : 67 services : 83 other :411
## office : 97
## sales : 38
## management: 55
## union married
## no :438 no :184
## yes: 96 yes:350
##
##
##
## | Name | CPS1985 |
| Number of rows | 534 |
| Number of columns | 11 |
| _______________________ | |
| Column type frequency: | |
| factor | 7 |
| numeric | 4 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| ethnicity | 0 | 1 | FALSE | 3 | cau: 440, oth: 67, his: 27 |
| region | 0 | 1 | FALSE | 2 | oth: 378, sou: 156 |
| gender | 0 | 1 | FALSE | 2 | mal: 289, fem: 245 |
| occupation | 0 | 1 | FALSE | 6 | wor: 156, tec: 105, off: 97, ser: 83 |
| sector | 0 | 1 | FALSE | 3 | oth: 411, man: 99, con: 24 |
| union | 0 | 1 | FALSE | 2 | no: 438, yes: 96 |
| married | 0 | 1 | FALSE | 2 | yes: 350, no: 184 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| wage | 0 | 1 | 9.02 | 5.14 | 1 | 5.25 | 7.78 | 11.25 | 44.5 | ▇▃▁▁▁ |
| education | 0 | 1 | 13.02 | 2.62 | 2 | 12.00 | 12.00 | 15.00 | 18.0 | ▁▁▂▇▃ |
| experience | 0 | 1 | 17.82 | 12.38 | 0 | 8.00 | 15.00 | 26.00 | 55.0 | ▇▇▃▂▁ |
| age | 0 | 1 | 36.83 | 11.73 | 18 | 28.00 | 35.00 | 44.00 | 64.0 | ▆▇▅▃▃ |
だが、これらはそのまま論文などに掲載することはできない。そこで、stargazerパッケージを使い、要約統計量の表を作成する。
まず、text形式で直接表示する。
##
## ===========================================
## Statistic N Mean St. Dev. Min Max
## -------------------------------------------
## wage 534 9.024 5.139 1.000 44.500
## education 534 13.019 2.615 2 18
## experience 534 17.822 12.380 0 55
## age 534 36.833 11.727 18 64
## -------------------------------------------次に、html形式で出力し、フォルダに保存する。
#htmlで表示し、summary_CPS1985.htmlという名前での保存
stargazer(CPS1985,
type = "html",
out = "summary_CPS1985.html")| Statistic | N | Mean | St. Dev. | Min | Max |
| wage | 534 | 9.024 | 5.139 | 1.000 | 44.500 |
| education | 534 | 13.019 | 2.615 | 2 | 18 |
| experience | 534 | 17.822 | 12.380 | 0 | 55 |
| age | 534 | 36.833 | 11.727 | 18 | 64 |
Wordで作成した論文に張り付ける場合はhtml形式のほうがよい。またstargazerの表はもっと細かく調整できる。ヘルプを参照したりググったりして、自分用にカスタマイズした表を作ることができる。また、latexを使った出力も可能である。
（※stargazer ver.5.2.3現在、{tibble}パッケージのデータフレームを投入すると結果が表示されないバグがあるので注意）
##
## =================================
## Statistic N Mean St. Dev. Min Max
## =================================5 回帰分析の基礎
回帰分析（regression analysis）は、ある変数\(Y\)（被説明変数）を、別の変数\(X\)（説明変数）で説明することで両変数の関係を分析する手法である。
最も基本的な回帰モデルは次のような線形で説明変数が１つの回帰モデルである。
\[ Y_i = \beta_0 + \beta_1 X_i + u_i \]
こうした線形の回帰モデルは線形回帰（linear regression）と呼ばれ、また説明変数\(X_i\)が１つの回帰モデルは単回帰モデル（simple regression model）と呼ばれる。
5.1 最小二乗法
線形回帰のパラメータ\(\beta_0, \beta_1\)はデータ\(X_i,Y_i\)から推定する。最も有名で強力な推定方法は最小二乗法（ordinary least squares method, OLS method）である。この方法では、実測値\(Y_i\)と予測値\(\hat{Y}_i\)との残差二乗和\(\sum_{i=1}^n (Y_i - \hat{Y}_i)^2\)を最小化するように推定量\(\hat{\beta}_0, \hat{\beta}_1\)を求める。
最小二乗法による線形回帰はlm()関数で行うことができる。
例えば、被説明変数を賃金\(\text{Wage}\)、説明変数を教育年数\(\text{Education}\)とする
\[
\text{Wage}_i = \beta_0 + \beta_1 \text{Education}_i + u_i
\]
のような回帰モデルを作って推定したい場合は、lm()のformula引数にwage ~ educationと指定し、使用するデータフレームをdata引数に指定すればよい。
##
## Call:
## lm(formula = wage ~ education, data = CPS1985)
##
## Coefficients:
## (Intercept) education
## -0.7460 0.7505ここでCoefficients:の部分に表示されている(Intercept)は切片の係数\(\beta_0\)で、educationは傾き係数\(\beta_1\)である。
推定したlmオブジェクトは、summary()関数でより詳細な情報を見ることができる。
##
## Call:
## lm(formula = wage ~ education, data = CPS1985)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.911 -3.260 -0.760 2.240 34.740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.74598 1.04545 -0.714 0.476
## education 0.75046 0.07873 9.532 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.754 on 532 degrees of freedom
## Multiple R-squared: 0.1459, Adjusted R-squared: 0.1443
## F-statistic: 90.85 on 1 and 532 DF, p-value: < 2.2e-165.2 回帰分析の結果表示
summary()による結果表示は論文に使用できるような整った形ではない。そこで再び{stargazer}パッケージを使って整形する。
stargazer()にlm()関数のオブジェクトを入れ、type引数に"text"を指定することで以下のように表示される
##
## ===============================================
## Dependent variable:
## ---------------------------
## wage
## -----------------------------------------------
## education 0.750***
## (0.079)
##
## Constant -0.746
## (1.045)
##
## -----------------------------------------------
## Observations 534
## R2 0.146
## Adjusted R2 0.144
## Residual Std. Error 4.754 (df = 532)
## F Statistic 90.852*** (df = 1; 532)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01type = "html", out = "reg1.html"とおけばhtml形式で保存できる
6 対数変換
データ分析の過程では、データを対数変換したい場合もある。
例えば、賃金のように「低い人（低所得者層）が多く、高い人（高所得者層）が少ない」ような分布はL字型のような右に歪んだ分布になる。こうした分布は対数変換を行うことで正規分布に近づく。
g1 <- ggplot(CPS1985, aes(x = wage)) +
geom_histogram(bins = 10)
g2 <- ggplot(CPS1985, aes(x = log(wage))) +
geom_histogram(bins = 10)
grid.arrange(g1, g2) # gridExtraパッケージによるグラフの結合。他にもpatchworkパッケージなどがある。
また、対数変換を行うことによって係数の解釈も次のように変わる。
| モデル | 係数の解釈 |
|---|---|
| \(Y = \beta_0 + \beta_1 X\) | 「\(X\)が1単位増加すると、\(Y\)が\(\beta_1\)単位増加する」 |
| \(Y = \beta_0 + \beta_1 \ln(X)\) | 「\(X\)が1%増加すると、\(Y\)が\(\beta_1 / 100\)単位増加する」 |
| \(\ln(Y) = \beta_0 + \beta_1 X\) | 「\(X\)が1単位増加すると、\(Y\)が\((\beta_1 \times 100)\)%増加する」 |
| \(\ln(Y) = \beta_0 + \beta_1 \ln(X)\) | 「\(X\)が1%増加すると、\(Y\)が\(\beta_1\)%増加する」 |
回帰分析での変数の対数変換はlm()関数内で変数をlog()で囲うだけでよい（なお、log()はデフォルトでは自然対数\(\ln(\cdot)\)となる）。
先程の単回帰モデルの被説明変数\(Wage\)を対数変換した回帰モデル
\[ \ln(\text{Wage}_i) = \beta_0 + \beta_1 \text{Education}_i + u_i \]
は以下のコードで実行できる。
対数変換をしなかった線形回帰の結果regと今回のreg_logwageをstargazer()に入れることで、まとめて結果表示できる。
##
## ===========================================================
## Dependent variable:
## ----------------------------
## wage log(wage)
## (1) (2)
## -----------------------------------------------------------
## education 0.750*** 0.077***
## (0.079) (0.008)
##
## Constant -0.746 1.060***
## (1.045) (0.107)
##
## -----------------------------------------------------------
## Observations 534 534
## R2 0.146 0.145
## Adjusted R2 0.144 0.143
## Residual Std. Error (df = 532) 4.754 0.489
## F Statistic (df = 1; 532) 90.852*** 90.006***
## ===========================================================
## Note: *p<0.1; **p<0.05; ***p<0.017 重回帰モデル
7.1 重回帰モデルとは
複数の説明変数を使用した回帰モデルを重回帰モデル（multiple regression model）という。
\[ Y_i = \beta_0 + \beta_1X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + u_i \]
先程の回帰モデルに新たな変数\(\text{Experience}\)（潜在経験年数）を追加した回帰モデル
\[ \ln(\text{Wage}_i) = \beta_0 + \beta_1 \text{Education}_i + \beta_2 \text{Experience}_i + u_i \]
をRで推定する場合、lm()関数の引数formulaに+ experienceを追加し、
と書けばよい。
##
## ===================================================================
## Dependent variable:
## -----------------------------------------------
## log(wage)
## (1) (2)
## -------------------------------------------------------------------
## education 0.077*** 0.096***
## (0.008) (0.008)
##
## experience 0.012***
## (0.002)
##
## Constant 1.060*** 0.594***
## (0.107) (0.124)
##
## -------------------------------------------------------------------
## Observations 534 534
## R2 0.145 0.211
## Adjusted R2 0.143 0.209
## Residual Std. Error 0.489 (df = 532) 0.470 (df = 531)
## F Statistic 90.006*** (df = 1; 532) 71.214*** (df = 2; 531)
## ===================================================================
## Note: *p<0.1; **p<0.05; ***p<0.017.2 ミンサー型賃金関数と２乗項
ミンサー型賃金関数
賃金\(\text{Wage}\)を教育年数\(\text{Education}\)、潜在経験年数\(\text{Experience}\)、潜在経験年数の２乗\(\text{Experience}^2\)、観察不能な賃金決定要因の関数\(u\)として示す式
\[ \ln(\text{Wage}_i) = \beta_0 + \beta_1 \text{Education}_i + \beta_2 \text{Experience}_i + \beta_3 \text{Experience}_i^2 + u_i \]
はミンサー型賃金関数と呼ばれ、世界各国の様々な時点の賃金分布をよく説明することが知られている。
２乗項を含む重回帰モデル
ミンサー型賃金関数を実際のデータにあてはめて分析する場合は、
- 潜在経験年数の２乗\(\text{Experience}^2\)の変数を作成し
- 重回帰モデルとして分析する
という手順で分析すればよい。
まず、２乗の変数を作成する。
mutate(新しい変数の名前 = 既存の変数の名前^2)のように書くことで任意の変数の2乗項を作成できる。
次に、このデータフレームdfと新たな変数名experience2を使ってlm()関数を書く。
text形式で分析結果表を作成すると以下のようになる。
##
## ===========================================================================================
## Dependent variable:
## -----------------------------------------------------------------------
## log(wage)
## (1) (2) (3)
## -------------------------------------------------------------------------------------------
## education 0.077*** 0.096*** 0.090***
## (0.008) (0.008) (0.008)
##
## experience 0.012*** 0.035***
## (0.002) (0.006)
##
## experience2 -0.001***
## (0.0001)
##
## Constant 1.060*** 0.594*** 0.520***
## (0.107) (0.124) (0.124)
##
## -------------------------------------------------------------------------------------------
## Observations 534 534 534
## R2 0.145 0.211 0.238
## Adjusted R2 0.143 0.209 0.234
## Residual Std. Error 0.489 (df = 532) 0.470 (df = 531) 0.462 (df = 530)
## F Statistic 90.006*** (df = 1; 532) 71.214*** (df = 2; 531) 55.229*** (df = 3; 530)
## ===========================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01html形式で分析結果表を作成すると以下のようになる。なおRMarkdown（このウェブ資料）ではNote以下の部分が正しく表示されていないので、出力・保存されたhtmlファイルを見ること。
| Dependent variable: | |||
| log(wage) | |||
| (1) | (2) | (3) | |
| education | 0.077*** | 0.096*** | 0.090*** |
| (0.008) | (0.008) | (0.008) | |
| experience | 0.012*** | 0.035*** | |
| (0.002) | (0.006) | ||
| experience2 | -0.001*** | ||
| (0.0001) | |||
| Constant | 1.060*** | 0.594*** | 0.520*** |
| (0.107) | (0.124) | (0.124) | |
| Observations | 534 | 534 | 534 |
| R2 | 0.145 | 0.211 | 0.238 |
| Adjusted R2 | 0.143 | 0.209 | 0.234 |
| Residual Std. Error | 0.489 (df = 532) | 0.470 (df = 531) | 0.462 (df = 530) |
| F Statistic | 90.006*** (df = 1; 532) | 71.214*** (df = 2; 531) | 55.229*** (df = 3; 530) |
| Note: | p<0.1; p<0.05; p<0.01 | ||
8 不均一分散への対応
8.1 不均一分散に対して頑健な標準誤差
lmによるOLS推定における標準誤差は、誤差項の分散が均一（均一分散）であると仮定している。しかし、実際の計量経済分析において均一分散の仮定が満たされていると想定できるケースはまれであり、通常は、不均一分散に対して頑健な標準誤差(ロバスト標準誤差,
heteroskedasticity robust standard error)を用いる。
8.2 lm_robust
ロバスト標準誤差は、lm()の代わりに{estimatr}パッケージのlm_robust()を用いて計算し、分析結果表に組み込むことができる。
パッケージを読み込み、lmと同じような形でlm_robustを用いる、ロバスト標準誤差にはいくつかの種類があるが、ここでは
HC1と指定する。
# ロバスト標準誤差を用いたOLS推定
reg_logwage_r <-
lm_robust(formula = log(wage) ~ education,
data = CPS1985, se_type = "HC1")
reg_multiple_r <-
lm_robust(formula = log(wage) ~ education + experience,
data = CPS1985,
se_type = "HC1")
reg_mincer_r <-
lm_robust(formula = log(wage) ~ education + experience + experience2,
data = df,
se_type = "HC1")なお、HC1は計量経済学でもっともよく使われるロバスト標準誤差であり、経済学者によく使用される計量分析ソフトのStataのロバスト標準誤差と同一のものであり、lm_robustでは、HC1のかわりにstataと表記してもよい。
推定結果は、lmと同じくsummary()関数で見ることができる。
##
## Call:
## lm_robust(formula = log(wage) ~ education, data = CPS1985, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 1.05989 0.105624 10.035 8.107e-22 0.8524 1.26738 532
## education 0.07676 0.008072 9.509 6.607e-20 0.0609 0.09262 532
##
## Multiple R-squared: 0.1447 , Adjusted R-squared: 0.1431
## F-statistic: 90.42 on 1 and 532 DF, p-value: < 2.2e-16##
## Call:
## lm_robust(formula = log(wage) ~ education + experience, data = CPS1985,
## se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 0.59417 0.117546 5.055 5.940e-07 0.36326 0.8251 531
## education 0.09641 0.008030 12.007 1.500e-29 0.08064 0.1122 531
## experience 0.01177 0.001794 6.563 1.255e-10 0.00825 0.0153 531
##
## Multiple R-squared: 0.2115 , Adjusted R-squared: 0.2085
## F-statistic: 77.27 on 2 and 531 DF, p-value: < 2.2e-16##
## Call:
## lm_robust(formula = log(wage) ~ education + experience + experience2,
## data = df, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 0.5203218 0.1218633 4.270 2.320e-05 0.2809274 0.7597161 530
## education 0.0897561 0.0081694 10.987 1.937e-25 0.0737077 0.1058045 530
## experience 0.0349403 0.0061116 5.717 1.812e-08 0.0229343 0.0469464 530
## experience2 -0.0005362 0.0001357 -3.952 8.783e-05 -0.0008028 -0.0002697 530
##
## Multiple R-squared: 0.2382 , Adjusted R-squared: 0.2339
## F-statistic: 53.46 on 3 and 530 DF, p-value: < 2.2e-16論文などで使用するための分析結果表については、lm_robustとstargazerの相性はよくないため、そのままstargazerを使うことができない。それでも両者を使う方法は参考文献のリンク先に記されているが、ここでは分析結果表のための別のパッケージである{texreg}を使う。
# 結果の出力
screenreg(list(reg_logwage_r,
reg_multiple_r,
reg_mincer_r),
include.ci = FALSE,
digits = 3)##
## ==================================================
## Model 1 Model 2 Model 3
## --------------------------------------------------
## (Intercept) 1.060 *** 0.594 *** 0.520 ***
## (0.106) (0.118) (0.122)
## education 0.077 *** 0.096 *** 0.090 ***
## (0.008) (0.008) (0.008)
## experience 0.012 *** 0.035 ***
## (0.002) (0.006)
## experience2 -0.001 ***
## (0.000)
## --------------------------------------------------
## R^2 0.145 0.211 0.238
## Adj. R^2 0.143 0.209 0.234
## Num. obs. 534 534 534
## RMSE 0.489 0.470 0.462
## ==================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05htmlもしくはWordで出力する場合は、次のようにする。
9 （参考）不均一分散への対応 (lm_robustを使わない方法)
lm_robust()登場前のRでのロバスト標準誤差の計算および分析結果表への組み込みの仕方を紹介する。
9.1 vcovHC
{sandwich}パッケージのvcovHC関数を用いる。なおsandwichという名称は、ロバスト標準誤差の推定量の式の形状から来ている。
まず{sandwich}パッケージをインストールして読み込む。
上記の3種類の重回帰モデルの推定結果の標準誤差を、ロバスト標準誤差に置き換える。
まず、vcovHC関数を用いてロバスト標準誤差を得る。
# (1) reg_lowage
vcov_reg_lowage <- vcovHC(reg_logwage,
type = "HC1") # 分散共分散行列
robustSE_reg_logwage <- sqrt(diag(vcov_reg_lowage)) # ロバスト標準誤差（分散共分散行列の対角成分の平方根）
# (2) reg_multiple
vcov_reg_multiple <- vcovHC(reg_multiple,
type = "HC1") # 分散共分散行列
robustSE_reg_multiple <- sqrt(diag(vcov_reg_multiple)) # ロバスト標準誤差（分散共分散行列の対角成分の平方根）
# (3) reg_mincer
vcov_reg_mincer <- vcovHC(reg_mincer,
type = "HC1") # 分散共分散行列
robustSE_reg_mincer <- sqrt(diag(vcov_reg_mincer)) # ロバスト標準誤差（分散共分散行列の対角成分の平方根）ここで、type = "HC1"とは、ロバスト標準誤差の種類の指定であり、HC1は計量経済学でもっともよく使われるものである。
{stargazer}で表を作成する際には、以下のように標準誤差をロバスト標準誤差に置き換える。
stargazer(reg_logwage,
reg_multiple,
reg_mincer,
se = list(robustSE_reg_logwage,
robustSE_reg_multiple,
robustSE_reg_mincer), # robust SEに置き換え
type = "text")##
## ===========================================================================================
## Dependent variable:
## -----------------------------------------------------------------------
## log(wage)
## (1) (2) (3)
## -------------------------------------------------------------------------------------------
## education 0.077*** 0.096*** 0.090***
## (0.008) (0.008) (0.008)
##
## experience 0.012*** 0.035***
## (0.002) (0.006)
##
## experience2 -0.001***
## (0.0001)
##
## Constant 1.060*** 0.594*** 0.520***
## (0.106) (0.118) (0.122)
##
## -------------------------------------------------------------------------------------------
## Observations 534 534 534
## R2 0.145 0.211 0.238
## Adjusted R2 0.143 0.209 0.234
## Residual Std. Error 0.489 (df = 532) 0.470 (df = 531) 0.462 (df = 530)
## F Statistic 90.006*** (df = 1; 532) 71.214*** (df = 2; 531) 55.229*** (df = 3; 530)
## ===========================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01通常は、均一分散を仮定した標準誤差よりもロバスト標準誤差のほうが大きくなる傾向があるが、今回はほとんど変わらない。なお回帰係数の推定値は全く同じとなる。
9.2 不均一分散に対して頑健なF検定
上の表において、F統計量は均一分散を仮定したもののままである。下記では、lmでの回帰分析結果を利用して、不均一分散に対して頑健なF検定を行う方法だけ記す。
次に、F検定に使用する「制約なし」の回帰モデル（被説明変数を定数項のみに回帰するモデル）を実行し、オブジェクトに保存する。
均一分散を仮定したF検定は、{AER}パッケージのwaldtest関数を使って、以下のような形で行うことができる。上記の（均一分散下での）回帰分析でのF検定の結果と一致することが確認できる。
## Wald test
##
## Model 1: log(wage) ~ education
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 532
## 2 533 -1 90.006 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Wald test
##
## Model 1: log(wage) ~ education + experience
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 531
## 2 533 -2 71.214 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Wald test
##
## Model 1: log(wage) ~ education + experience + experience2
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 530
## 2 533 -3 55.229 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1一方、不均一分散に対して頑健なF検定は、以下のように実施できる。
# (1) reg_lowage
waldtest(reg_logwage,
reg_no_X,
vcov = vcovHC(reg_logwage,
type = "HC1")) # HC1を用いたF検定## Wald test
##
## Model 1: log(wage) ~ education
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 532
## 2 533 -1 90.416 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# (2) reg_multiple
waldtest(reg_multiple,
reg_no_X,
vcov = vcovHC(reg_multiple,
type = "HC1")) # HC1を用いたF検定## Wald test
##
## Model 1: log(wage) ~ education + experience
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 531
## 2 533 -2 77.265 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# (3) reg_mincer
waldtest(reg_mincer,
reg_no_X,
vcov = vcovHC(reg_mincer,
type = "HC1")) # HC1を用いたF検定## Wald test
##
## Model 1: log(wage) ~ education + experience + experience2
## Model 2: log(wage) ~ 1
## Res.Df Df F Pr(>F)
## 1 530
## 2 533 -3 53.462 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1F値が均一分散を仮定したものとは異なっていることが確認できる。
なお、上記の方法についてのさらに詳しい解説は、このページなどを参照すること。ただし上記では、このサイトとは若干異なる方法でF検定を行っている。
10 参考文献
回帰分析や均一分散/不均一分散については、「計量経済学１」の講義資料などを参照すること。
lm_robstを利用した分析については、以下のサイトを参考にした。
Getting started using estimatr