もう1つ回帰分析の例として，教育の収益率の推定をやってみよう。
ここで用いるのは，Wooldridgeのテキストで使われている”wage1”というデータ。元データは，1976年のアメリカのCurrent Population Survey。古いデータではあるが，練習用には便利。
散布図のマトリクスを作成するためにGGallyパッケージを利用するので，インストールしておこう。
library(tidyverse)
library(GGally)
library(stargazer)
Please cite as: Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables. R package version 5.2.3. https://CRAN.R-project.org/package=stargazer library(wooldridge)
data("wage1")データの詳細は，ヘルプで確認できる。
?wage1glimpse(wage1)Rows: 526
Columns: 24
$ wage <dbl> 3.10, 3.24, 3.00, 6.00, 5.30, 8.75, 11.25, 5.00, 3.60, 18.18,…
$ educ <int> 11, 12, 11, 8, 12, 16, 18, 12, 12, 17, 16, 13, 12, 12, 12, 16…
$ exper <int> 2, 22, 2, 44, 7, 9, 15, 5, 26, 22, 8, 3, 15, 18, 31, 14, 10, …
$ tenure <int> 0, 2, 0, 28, 2, 8, 7, 3, 4, 21, 2, 0, 0, 3, 15, 0, 0, 10, 0, …
$ nonwhite <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ female <int> 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1…
$ married <int> 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0…
$ numdep <int> 2, 3, 2, 0, 1, 0, 0, 0, 2, 0, 0, 0, 2, 0, 1, 1, 0, 0, 3, 0, 0…
$ smsa <int> 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ northcen <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ south <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ west <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ construc <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ ndurman <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ trcommpu <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ trade <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ services <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ profserv <int> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1…
$ profocc <int> 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1…
$ clerocc <int> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0…
$ servocc <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0…
$ lwage <dbl> 1.1314021, 1.1755733, 1.0986123, 1.7917595, 1.6677068, 2.1690…
$ expersq <int> 4, 484, 4, 1936, 49, 81, 225, 25, 676, 484, 64, 9, 225, 324, …
$ tenursq <int> 0, 4, 0, 784, 4, 64, 49, 9, 16, 441, 4, 0, 0, 9, 225, 0, 0, 1…526人の労働者について24の変数が格納されている。変数の中には，基本的な変数の対数変換や，二乗なども含まれている (通常これらは自分で作成する)。
ここでは，賃金と修学年数との関係に注目する。賃金と修学年数の分布 (ヒストグラム)と散布図を描いてみよう。
wage1 %>%
ggplot(aes(educ)) +
geom_histogram(
bins = 20,
fill = "skyblue",
color = "black"
) +
labs(
title = "修学年数の分布",
x = "修学年数",
y = "人数"
)
wage1 %>%
ggplot(aes(wage)) +
geom_histogram(
bins = 20,
fill = "skyblue",
color = "black"
) +
labs(
title = "時間あたり賃金の分布",
x = "時間あたり賃金",
y = "人数"
)
wage1 %>%
ggplot(
aes(x = educ, y = wage)
) +
geom_point(color = "skyblue") +
labs(
title = "修学年数と時間あたり賃金の関係",
x = "修学年数",
y = "時間あたり賃金"
)
散布図を見ると，教育水準が高いほど時間あたり賃金も高いことがわかる。ここで，われわれが知りたいのは，修学年数が1年増加すると，どの程度賃金が上昇するかということである。
そのため，散布図で示された修学年数と時間あたり賃金との関係に，最も良く当てはまる直線 (回帰直線)を引いてみよう。
wage1 %>%
ggplot(
aes(x = educ, y = wage)
) +
geom_point(color = "skyblue") +
geom_smooth(
method = "lm",
formula = y ~ x
) +
labs(
title = "回帰直線",
x = "修学年数",
y = "時間あたり賃金"
)
回帰直線は右上がりであり，修学年数が増加すると時間あたり賃金が上昇するという関係を表している。回帰直線の式は，最小二乗法 (OLS)で求める。
wage1 %>%
lm(wage ~ educ, data = .) %>%
summary()
Call:
lm(formula = wage ~ educ, data = .)
Residuals:
Min 1Q Median 3Q Max
-5.3396 -2.1501 -0.9674 1.1921 16.6085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.90485 0.68497 -1.321 0.187
educ 0.54136 0.05325 10.167 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.378 on 524 degrees of freedom
Multiple R-squared: 0.1648, Adjusted R-squared: 0.1632
F-statistic: 103.4 on 1 and 524 DF, p-value: < 2.2e-16すなわち，回帰直線の式は以下で表される。 \[ \widehat{wage} = -0.905 + 0.541 educ \] したがって，修学年数が1年増加すると，時間あたり賃金は0.54ドル増加することがわかる。
次に，時間あたり賃金の対数値を被説明変数として回帰分析を行ってみよう。
wage1 %>%
ggplot(
aes(x = educ, y = lwage)
) +
geom_point(color = "skyblue") +
geom_smooth(
method = "lm",
formula = y ~ x
) +
labs(title = "修学年数と対数賃金の関係",
x = "修学年数",
y = "対数賃金")
wage1 %>%
lm(lwage ~ educ, data = .) %>%
summary()
Call:
lm(formula = lwage ~ educ, data = .)
Residuals:
Min 1Q Median 3Q Max
-2.21158 -0.36393 -0.07263 0.29712 1.52339
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.583773 0.097336 5.998 3.74e-09 ***
educ 0.082744 0.007567 10.935 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4801 on 524 degrees of freedom
Multiple R-squared: 0.1858, Adjusted R-squared: 0.1843
F-statistic: 119.6 on 1 and 524 DF, p-value: < 2.2e-16すなわち，回帰直線の式は， \[ \widehat{\log (wage)} = 0.584 + 0.083 educ \] 修学年数が1年増加すれば，時間あたり賃金の対数値が0.08増加する。これは，修学年数が1年増加すれば賃金が約8%増加するとことを意味する。このように，被説明変数として賃金の対数値を用いることで，修学年数の係数を1年あたりの教育の収益率と見なすことができる。
従属変数と独立変数にレベル (もとの値)を使うか，log (対数)を使うかによって，回帰係数の解釈は以下のようになる。
| 従属変数 | 独立変数 | 係数\(\beta\)の解釈 |
|---|---|---|
| \(y\) | \(x\) | 限界効果 (\(x\)が1単位変化 → \(y\)は\(\beta\)単位変化) |
| \(\log y\) | \(x\) | \(x\)が1単位変化 → \(y\)は\(\beta \times 100\)パーセント変化 |
| \(y\) | \(\log x\) | \(x\)が1%変化 → \(y\)は\(\beta / 100\)単位変化 |
| \(\log y\) | \(\log x\) | 弾力性 (\(x\)が1%変化 → \(y\)は\(\beta\)%変化) |
時間あたり賃金に影響を与えるのは，修学年数だけではない。たとえば，労働市場における経験年数は時間あたり賃金を決める重要な要素だと考えられる。そこで，以下のような回帰式を考える。
\[ \widehat{\log (wage)} = \hat{\beta_0} + \hat{\beta_1} educ + \hat{\beta_2} exper\]
この式も最小二乗法で推定することができるが，その前に変数間の関係を確認しておこう。GGallyパッケージのggpairs関数を使うと，変数間の関係を視覚的に確認することができる。
ggpairs(
wage1[ c("lwage", "educ", "exper")],
diag = list(continuous = wrap("barDiag", bins = 20)),
upper = list(continuous = wrap("points"), combo = wrap("facethist")),
lower = list(continuous = wrap("cor"), combo = wrap("box"))
)
経験年数と対数賃金の相関係数は0.11であるが，経験年数と対数賃金の関係は線形というより，逆U字型である。すなわち，経験年数が増加すると賃金は上昇するが，その上昇率はだんだん小さくなっていき，一定の年齢を超えると賃金は低下する。また，経験年数が増加するほど修学年数が減少していることがわかる。これは，教育水準が上昇しており，経験年数の短い若い世代の労働者ほど修学年数が長いためである。
それでは最小二乗法による重回帰分析をやってみよう。
wage1 %>%
lm(lwage ~ educ + exper, data = .) %>%
stargazer(type = "text")
===============================================
Dependent variable:
---------------------------
lwage
-----------------------------------------------
educ 0.098***
(0.008)
exper 0.010***
(0.002)
Constant 0.217**
(0.109)
-----------------------------------------------
Observations 526
R2 0.249
Adjusted R2 0.246
Residual Std. Error 0.461 (df = 523)
F Statistic 86.862*** (df = 2; 523)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01推定結果を式で表すと， \[ \widehat{\log (wage)} = 0.53 + 0.098 educ + 0.010 exper \]
この式が意味するのは，
ということである。
単回帰の場合と重回帰の場合とで，修学年数の係数が異なることに注意しよう。これは，修学年数と経験年数が相関しているためである。若い労働者ほど修学年数が長いので，単回帰では修学年数が増加することの影響に，経験年数が短いことの影響が含まれてしまい，それによって修学年数の効果が小さく推定される。一方で，重回帰では，経験年数を一定とした修学年数の影響を計測しているため，係数が単回帰の場合に比べて大きくなっている。
経験年数と対数賃金の非線形 (逆U字型)の関係を捉えるため，経験年数の二乗項を加えた回帰式を考えよう。
\[ \widehat{\log (wage)} = \hat{\beta_0} + \hat{\beta_1} educ + \hat{\beta_2} exper + \hat{\beta_3}exper^2\]
この回帰式も最小二乗法によって推定することができる。そのためには，\(exper\)の二乗の変数を作成して，独立変数に加えれば良い。ここでは，データセットにもともと\(exper\)の二乗を計算した\(expersq\)が含まれているので，それを使う。
wage1 %>%
lm(lwage ~ educ + exper + expersq, data = .) %>%
stargazer(type = "text")
===============================================
Dependent variable:
---------------------------
lwage
-----------------------------------------------
educ 0.090***
(0.007)
exper 0.041***
(0.005)
expersq -0.001***
(0.0001)
Constant 0.128
(0.106)
-----------------------------------------------
Observations 526
R2 0.300
Adjusted R2 0.296
Residual Std. Error 0.446 (df = 522)
F Statistic 74.668*** (df = 3; 522)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01推定結果から，\(exper\)の係数は0.041，\(exper^2\)の係数は-0.00071である。対数賃金と経験年数の関係は二次関数であり，経験年数がゼロのときには，経験年数が1年増加すれば4.1%賃金が上昇する。しかし，その上昇率は1年につき0.14%ずつ逓減していく。経験年数が約28年になると上昇率はマイナスとなり，賃金は経験年数とともに下がり始める。
対数賃金を独立変数を修学年数と経験年数，経験年数の二乗で説明しようとする回帰式をミンサー型賃金関数という。ミンサー型賃金関数は，どのような国のデータにも当てはまりが良く，一定の条件の下で修学年数の係数\(\hat{\beta_1}\)は，教育の内部収益率として解釈できる。
ミンサー型賃金関数は修学年数と年齢が記録されたデータがあれば推定できるため，多くの研究で用いられてきた。しかし，推定された\(\hat{\beta_1}\)を教育の賃金に対する因果効果と見なすためには，いくつかの条件が満たされる必要がある。
いま，真のモデルが以下のように表されるとしよう。
\[\log (wage) = \beta_0 + \beta_1 educ + \beta_2 exper + \beta_3 exper^2 + u\] これは真のモデルなので，\(\beta_0\), \(\beta_1\), \(\beta_2\), \(\beta_3\)は未知のパラメータである。最小二乗法は，これらの未知のパラメータを推定する方法の1つである。また，\(u\)は誤差項と呼ばれ，修学年数や経験年数以外に賃金に影響を与えるすべての要因が含まれている。
最小二乗法によってこれらのパラメータをバイアスなく推定できるのは，独立変数と誤差項\(u\)とが無相関であるときである。しかし，この条件は多くの場合満たされない。
たとえば，経験年数が観測不可能だとすれば，経験年数の影響はすべて誤差項\(u\)に含まれることになる。さきほど見たように修学年数と経験年数には相関があるため，観測データを使って最小二乗法によって推定された\(\hat{\beta_1}\)は，修学年数の効果だけではなく，経験年数の効果も含んでしまう。
実際に経験年数はデータに記録されているため，独立変数に経験年数を含めることで，経験年数を一定として修学年数が1年増加したときに賃金がどう変化するかを推定することができる。また，さまざまな個人属性を独立変数に追加することで，観察可能な要因をコントロールすることができる (データセットの内容をよく見て，独立変数を追加してみよう)。
しかし，賃金に影響を与える要因の多くは観察不可能である。たとえば，生まれながらの能力はデータには記録されていないが，賃金に影響を与えていると考えられる。加えて，能力の高い人は修学年数も長くなる傾向があるだろう。すると，推定された\(\hat{\beta_1}\)は，修学年数の効果だけではなく，能力が高いことの効果を含んでいることになる。このように，観察不可能な要因が独立変数と誤差項との間の相関をもたらすことを欠落変数バイアス (Omitted variable bias)という。とくに，ここでは能力が欠落変数となっているので，それによりもたらされるバイアスを能力バイアスという。
大学卒の労働者は高校卒の労働者よりも賃金が高いが，その理由は大学で教育を受けることにより知識や技能といった人的資本が蓄積されるからなのか，それとも大学卒の労働者は高校卒の労働者よりも生まれながらの能力が高いからなのかは，修学年数と賃金のデータだけを見てもわからない。それは，回帰分析を使っても同じことである。
説明変数と誤差項に相関があるために，パラメータの推定にバイアスが生じることを内生性の問題という。理想的な実験データを除けば，およそどのようなデータを用いて回帰分析を行っても，少なからず内生性の問題は存在するが，それによって具体的にどのような問題が生じるかを考えることが重要である。