pacman::p_load(
tidyverse,
skimr,
estimatr,
AER,
modelsummary,
stargazer,
gridExtra
)Rで学ぶ計量経済学と機械学習 4
回帰分析
This version: 2026-10-14
1 本講義の目的
- Rを用いた回帰分析の方法を学ぶ。
- 要約統計量や回帰分析の結果の実践的な出力方法を学ぶ
2 パッケージの読み込み
まず、{tidyverse}や{AER}や{modelsummary}など、今回用いるパッケージのインストールと読み込みを行う。新しい関数については、適宜説明していく。
3 データの用意
以下ではRで回帰分析を行っていくが、今回は{AER}パッケージ（Applied Econometrics with Rパッケージ）に収録されているCPS1985というデータセットを例に使用する。
今回使用するCPS1985データセットは、1985年アメリカの時給や教育年数などの変数が入ったものである。
data()関数でデータセットの読み込みを行う。
# データの読み込み
data("CPS1985")
head(CPS1985) 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)を使えばよい。
summary(CPS1985) 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
skimr::skim(CPS1985)| 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 | ▆▇▅▃▃ |
だが、これらはそのまま論文などに掲載することはできない。そこで、{modelsummary}パッケージのdatasummary_skim()関数を使い、要約統計量の表を作成する。
datasummary_skim(CPS1985, type = "numeric")| Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|
| wage | 238 | 0 | 9.0 | 5.1 | 1.0 | 7.8 | 44.5 |  |
| education | 17 | 0 | 13.0 | 2.6 | 2.0 | 12.0 | 18.0 |  |
| experience | 52 | 0 | 17.8 | 12.4 | 0.0 | 15.0 | 55.0 |  |
| age | 47 | 0 | 36.8 | 11.7 | 18.0 | 35.0 | 64.0 |  |
type = "numeric"は数値の変数だけを表にするという指定である（type = "categorical"とすると、genderやoccupationなどのカテゴリ変数の度数の表になる）。
output引数にファイル名を指定すると、表をファイルとして保存できる。拡張子を.docxにすればWord形式、.htmlにすればhtml形式で保存される。
# Word形式で保存（Wordで作成した論文に貼り付ける場合はこちらが便利）
datasummary_skim(CPS1985, type = "numeric", output = "summary_CPS1985.docx")
# html形式で保存
datasummary_skim(CPS1985, type = "numeric", output = "summary_CPS1985.html"){modelsummary}の表はもっと細かく調整できる。ヘルプや公式サイトを参照して、自分用にカスタマイズした表を作ることができる。また、LaTeX形式（拡張子.tex）での出力も可能である。
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引数に指定すればよい。
# Regression
reg <- lm(formula = wage ~ education,
data = CPS1985) # 教育年数が賃金を決めるモデル
reg
Call:
lm(formula = wage ~ education, data = CPS1985)
Coefficients:
(Intercept) education
-0.7460 0.7505 ここでCoefficients:の部分に表示されている(Intercept)は切片の係数\(\beta_0\)で、educationは傾き係数\(\beta_1\)である。
推定したlmオブジェクトは、summary()関数でより詳細な情報を見ることができる。
summary(reg)
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()による結果表示は論文に使用できるような整った形ではない。そこで再び{modelsummary}パッケージを使って整形する。
modelsummary()にlm()関数の結果を入れると、以下のような表が作成される。stars = TRUEを指定すると、係数の有意水準を示す*がつく。括弧内の数値は標準誤差である。
modelsummary(reg, stars = TRUE)| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| (Intercept) | -0.746 |
| (1.045) | |
| education | 0.750*** |
| (0.079) | |
| Num.Obs. | 534 |
| R2 | 0.146 |
| R2 Adj. | 0.144 |
| AIC | 3184.4 |
| BIC | 3197.3 |
| Log.Lik. | -1589.209 |
| F | 90.852 |
| RMSE | 4.75 |
表の下部にはサンプルサイズ（Num.Obs.）や決定係数（R2）などの統計量が表示される。gof_omit引数を使うと、不要な統計量を省略できる。以下では、AIC・BIC・対数尤度（Log.Lik.）・RMSEを省略して表示する。
output = "reg1.docx"（Word形式）やoutput = "reg1.html"（html形式）とおけばファイルとして保存できる。
modelsummary(reg, # 回帰分析の結果
stars = TRUE, # 有意水準の*を表示
output = "reg1.docx") # 保存するファイル名（拡張子でファイル形式が決まる）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_logwage <- lm(formula = log(wage) ~ education,
data = CPS1985)対数変換をしなかった線形回帰の結果regと今回のreg_logwageをlist()でまとめてmodelsummary()に入れることで、まとめて結果表示できる。
modelsummary(list(reg, reg_logwage),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | -0.746 | 1.060*** |
| (1.045) | (0.107) | |
| education | 0.750*** | 0.077*** |
| (0.079) | (0.008) | |
| Num.Obs. | 534 | 534 |
| R2 | 0.146 | 0.145 |
| R2 Adj. | 0.144 | 0.143 |
| F | 90.852 | 90.006 |
7 重回帰モデル
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を追加し、
reg_multiple <- lm(formula = log(wage) ~ education + experience,
data = CPS1985)と書けばよい。
modelsummary(list(reg_logwage, reg_multiple),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 1.060*** | 0.594*** |
| (0.107) | (0.124) | |
| education | 0.077*** | 0.096*** |
| (0.008) | (0.008) | |
| experience | 0.012*** | |
| (0.002) | ||
| Num.Obs. | 534 | 534 |
| R2 | 0.145 | 0.211 |
| R2 Adj. | 0.143 | 0.209 |
| F | 90.006 | 71.214 |
7.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乗項を作成できる。
# mutate関数で新しい変数を追加したデータをdfオブジェクトに代入
df <- CPS1985 |>
dplyr::mutate(experience2 = experience^2) 次に、このデータフレームdfと新たな変数名experience2を使ってlm()関数を書く。
reg_mincer <- lm(
formula = log(wage) ~ education + experience + experience2,
data = df
)3つのモデルの推定結果を並べると以下のようになる。
modelsummary(list(reg_logwage, reg_multiple, reg_mincer),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | (3) | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| (Intercept) | 1.060*** | 0.594*** | 0.520*** |
| (0.107) | (0.124) | (0.124) | |
| 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) | |||
| Num.Obs. | 534 | 534 | 534 |
| R2 | 0.145 | 0.211 | 0.238 |
| R2 Adj. | 0.143 | 0.209 | 0.234 |
| F | 90.006 | 71.214 | 55.229 |
Word形式で保存する場合は、以下のようにする。
modelsummary(list(reg_logwage, reg_multiple, reg_mincer),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE",
output = "04reg_table.docx")8 不均一分散への対応
8.1 不均一分散に対して頑健な標準誤差
lmによるOLS推定における標準誤差は、誤差項の分散が均一（均一分散）であると仮定している。しかし、実際の計量経済分析において均一分散の仮定が満たされていると想定できるケースはまれであり、通常は、不均一分散に対して頑健な標準誤差(ロバスト標準誤差, heteroskedasticity robust standard error)を用いる。
8.2 lm_robust
ロバスト標準誤差は、lm()の代わりに{estimatr}パッケージのlm_robust()を用いて計算し、分析結果表に組み込むことができる。
パッケージを読み込み、lmと同じような形でlm_robustを用いる、ロバスト標準誤差にはいくつかの種類があるが、ここでは HC1と指定する。
# ロバスト標準誤差を用いたOLS推定
reg_logwage_robust <-
lm_robust(formula = log(wage) ~ education,
data = CPS1985, se_type = "HC1")
reg_multiple_robust <-
lm_robust(formula = log(wage) ~ education + experience,
data = CPS1985,
se_type = "HC1")
reg_mincer_robust <-
lm_robust(formula = log(wage) ~ education + experience + experience2,
data = df,
se_type = "HC1")なお、HC1は計量経済学でもっともよく使われるロバスト標準誤差であり、経済学者によく使用される計量分析ソフトのStataのロバスト標準誤差と同一のものであり、lm_robustでは、HC1のかわりにstataと表記してもよい。
推定結果は、lmと同じくsummary()関数で見ることができる。
summary(reg_logwage_robust)
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-16summary(reg_multiple_robust)
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-16summary(reg_mincer_robust)
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-16lm_robust()の結果も、lm()の結果と同じようにmodelsummary()に入れて表にすることができる。
modelsummary(list(reg_logwage_robust, reg_multiple_robust, reg_mincer_robust),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | (3) | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| (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) | |||
| Num.Obs. | 534 | 534 | 534 |
| R2 | 0.145 | 0.211 | 0.238 |
| R2 Adj. | 0.143 | 0.209 | 0.234 |
8.3 modelsummaryでロバスト標準誤差を指定する
modelsummary()では、lm()の結果に対してvcov = "HC1"と指定するだけでも、ロバスト標準誤差を用いた表を作ることができる。係数の推定値は同じで、括弧内の標準誤差がロバスト標準誤差に置き換わる。表の最下部のStd.Errorsの行に、どの種類の標準誤差を用いたかが表示される。
modelsummary(list(reg_logwage, reg_multiple, reg_mincer),
vcov = "HC1",
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | (3) | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| (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) | |||
| Num.Obs. | 534 | 534 | 534 |
| R2 | 0.145 | 0.211 | 0.238 |
| R2 Adj. | 0.143 | 0.209 | 0.234 |
| F | 90.416 | 77.265 | 53.462 |
| Std.Errors | HC1 | HC1 | HC1 |
上のlm_robust()の表と、係数・標準誤差が一致していることが確認できる。また、この場合は表のF統計量（Fの行）も、不均一分散に対して頑健なものに置き換わっている（vcovを指定していない前節の表のF統計量とは値が異なる）。
Word形式で保存する場合は、以下のようにする。
modelsummary(list(reg_logwage, reg_multiple, reg_mincer),
vcov = "HC1",
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE",
output = "04reg_table_r.docx")9 （参考）ロバスト標準誤差を自分で計算する方法
ロバスト標準誤差が実際にどのように計算されているかを確認するため、{sandwich}パッケージを使って自分で計算する方法を紹介する（lm_robust()が登場する前は、この方法が一般的だった）。
9.1 vcovHC
{sandwich}パッケージのvcovHC関数を用いる。なおsandwichという名称は、ロバスト標準誤差の推定量の式の形状から来ている。
まず{sandwich}パッケージをインストールして読み込む。
pacman::p_load(sandwich)上記の3種類の重回帰モデルの推定結果の標準誤差を、ロバスト標準誤差に置き換える。
まず、vcovHC関数を用いてロバスト標準誤差を得る。
# (1) reg_logwage
vcov_reg_logwage <- vcovHC(reg_logwage,
type = "HC1") # 分散共分散行列
robustSE_reg_logwage <- sqrt(diag(vcov_reg_logwage)) # ロバスト標準誤差（分散共分散行列の対角成分の平方根）
# (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は計量経済学でもっともよく使われるものである。
modelsummary()のvcov引数に、計算した分散共分散行列をリストで渡すと、表の標準誤差がそれに置き換わる（Std.Errorsの行はCustomと表示される）。
modelsummary(list(reg_logwage, reg_multiple, reg_mincer),
vcov = list(vcov_reg_logwage,
vcov_reg_multiple,
vcov_reg_mincer), # 自分で計算した分散共分散行列
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | (3) | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| (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) | |||
| Num.Obs. | 534 | 534 | 534 |
| R2 | 0.145 | 0.211 | 0.238 |
| R2 Adj. | 0.143 | 0.209 | 0.234 |
| Std.Errors | Custom | Custom | Custom |
vcov = "HC1"と指定した前節の表と同じ標準誤差になっていることが確認できる。通常は、均一分散を仮定した標準誤差よりもロバスト標準誤差のほうが大きくなる傾向があるが、今回はほとんど変わらない。なお回帰係数の推定値は全く同じとなる。
9.2 不均一分散に対して頑健なF検定
前節でvcov = "HC1"と指定した表では、F統計量も不均一分散に対して頑健なものに置き換わっていた（なお、上のように分散共分散行列を自分で渡した場合は、F統計量は表示されない）。以下では、lmでの回帰分析結果を利用して、均一分散を仮定したF検定と不均一分散に対して頑健なF検定を自分で行う方法を記す。
まず、F検定に使用する「制約あり」の回帰モデル（被説明変数を定数項のみに回帰するモデル）を実行し、オブジェクトに保存する。
reg_no_X <- lm(log(wage) ~ 1,
data = CPS1985) # no_X:説明変数なし均一分散を仮定したF検定は、{AER}パッケージのwaldtest関数を使って、以下のような形で行うことができる。vcovを指定せずに作成した表（重回帰モデルの節の表）のF統計量と一致することが確認できる。
# (1) reg_logwage
waldtest(reg_logwage,
reg_no_X)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# (2) reg_multiple
waldtest(reg_multiple,
reg_no_X)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# (3) reg_mincer
waldtest(reg_mincer,
reg_no_X)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_logwage
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値が均一分散を仮定したものとは異なり、vcov = "HC1"と指定した表のF統計量と一致していることが確認できる。
なお、ロバスト標準誤差や頑健な検定についてのさらに詳しい解説は、Introduction to Econometrics with Rの5.4節などを参照すること。ただしそこでは、上記とは若干異なる方法（{lmtest}パッケージのcoeftest()や{car}パッケージのlinearHypothesis()）で検定を行っている。
10 （参考）stargazerを使う場合
{stargazer}パッケージも、回帰分析の結果表を作るパッケージとして長く使われてきた。現在は更新が止まっている（2022年のver.5.2.3が最新）が、古い解説記事などで今もよく見かけるので、基本的な使い方を紹介しておく。
type = "text"とすると、表をテキスト形式で表示する。
stargazer(reg_logwage,
reg_multiple,
reg_mincer,
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.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.01type = "html", out = "ファイル名.html"とすれば、html形式で保存できる。
ロバスト標準誤差を表示したい場合は、se引数に、上で計算したロバスト標準誤差をリストで渡す。
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ただし、この場合も表の下部のF統計量は均一分散を仮定したもののままであることに注意。また、lm_robust()の結果はそのままではstargazer()に入れることができない。
要約統計量の表もstargazer()にデータフレームを入れれば作成できるが、{tibble}パッケージのデータフレームを入れると結果が表示されないバグがあるので注意（ver.5.2.3時点）。
stargazer(CPS1985, type = "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
-------------------------------------------stargazer(tibble::tibble(CPS1985), type = "text") # tibble：表示されない
=================================
Statistic N Mean St. Dev. Min Max
=================================11 参考文献
回帰分析や均一分散/不均一分散については、「計量経済学１」の講義資料などを参照すること。
- Christoph Hanck, Martin Arnold, Alexander Gerber, Martin Schmelzer Introduction to Econometrics with R
- Stock & Watsonの計量経済学の教科書に沿って、Rでの実践方法を解説した無料のウェブ教科書（英語）。本講義と同じ
{AER}パッケージを使っており、回帰分析やロバスト標準誤差（5.4節）なども扱っている。
- Stock & Watsonの計量経済学の教科書に沿って、Rでの実践方法を解説した無料のウェブ教科書（英語）。本講義と同じ
- 矢内勇生 KUT 計量経済学応用
- 第1回でも紹介した、Rを使った計量経済学（統計的因果推論）の講義資料（日本語）。
- modelsummary 公式サイト
- 本講義で用いた
modelsummary()・datasummary_skim()の使い方や、表のカスタマイズ方法の公式解説（英語）。
- 本講義で用いた
- Getting started using estimatr
lm_robust()を含む{estimatr}パッケージの公式解説（英語）。
12 （おまけ）Pythonで同様の作業を行う場合
本講義の演習はRで行うが、参考として、ここまで学んだ内容に対応するPythonのコードを簡単に紹介する。
Pythonで計量経済学的な回帰分析を行う場合は、主にstatsmodelsというライブラリを使う。smf.ols()を使うと、Rのlm()とほぼ同じ書き方（"y ~ x1 + x2"という式）で回帰分析ができる。
# ターミナル上で実行（Rのinstall.packages()に相当）
# pip install pandas statsmodels matplotlib seabornimport numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.iolib.summary2 import summary_col
import matplotlib.pyplot as plt
import seaborn as sns12.1 データの用意と記述統計
Rのパッケージに収録されているデータセットは、sm.datasets.get_rdataset()でインターネット経由で読み込める（Rのdata("CPS1985")に相当）。
cps = sm.datasets.get_rdataset("CPS1985", "AER").data
cps.head()
cps.describe() # 要約統計量（Rの summary() に相当）（※Macで公式サイト版のPythonを使っている場合、ここでSSL: CERTIFICATE_VERIFY_FAILEDというエラーが出ることがある。その場合は、PythonのインストールフォルダにあるInstall Certificates.commandを一度実行すればよい。）
12.2 回帰分析
.fit()で推定し、.summary()で詳細な結果を表示する（Rのlm()とsummary()に相当）。
reg = smf.ols("wage ~ education", data=cps).fit()
reg.summary()12.3 対数変換・重回帰モデル・２乗項
式の中でnp.log()を使えば対数変換ができる。２乗項の変数は、あらかじめ列として追加しておく（Rのmutate()に相当）。
reg_logwage = smf.ols("np.log(wage) ~ education", data=cps).fit()
reg_multiple = smf.ols("np.log(wage) ~ education + experience", data=cps).fit()
cps["experience2"] = cps["experience"] ** 2
reg_mincer = smf.ols("np.log(wage) ~ education + experience + experience2", data=cps).fit()賃金の分布の対数変換前後の比較は、以下のように描ける（Rのgrid.arrange()に相当）。
fig, axes = plt.subplots(2, 1)
sns.histplot(cps["wage"], bins=10, ax=axes[0])
sns.histplot(np.log(cps["wage"]), bins=10, ax=axes[1])
plt.show()12.4 回帰分析の結果表示
複数のモデルを並べた表はsummary_col()で作成できる（Rのmodelsummary()に相当）。
summary_col([reg_logwage, reg_multiple, reg_mincer],
stars=True,
info_dict={"N": lambda x: f"{int(x.nobs)}",
"F": lambda x: f"{x.fvalue:.3f}"})なお、summary_col()の*の基準（* p<.1, ** p<.05, *** p<.01）は、modelsummary()の基準（+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001）とは異なるので注意。
12.5 ロバスト標準誤差
.fit()の中でcov_type="HC1"と指定すると、ロバスト標準誤差（HC1）を用いた推定結果になる（Rのlm_robust(..., se_type = "HC1")やmodelsummary(..., vcov = "HC1")に相当）。
reg_logwage_robust = smf.ols("np.log(wage) ~ education", data=cps).fit(cov_type="HC1")
reg_multiple_robust = smf.ols("np.log(wage) ~ education + experience", data=cps).fit(cov_type="HC1")
reg_mincer_robust = smf.ols("np.log(wage) ~ education + experience + experience2",
data=cps).fit(cov_type="HC1")
summary_col([reg_logwage_robust, reg_multiple_robust, reg_mincer_robust],
stars=True,
info_dict={"N": lambda x: f"{int(x.nobs)}",
"F": lambda x: f"{x.fvalue:.3f}"})この場合、.fvalue（F統計量）も不均一分散に対して頑健なものになり、Rのwaldtest(..., vcov = vcovHC(..., type = "HC1"))やmodelsummary(..., vcov = "HC1")のF統計量と一致する。
reg_logwage_robust.fvalue # 頑健なF統計量
reg_logwage.fvalue # 均一分散を仮定したF統計量