準備

library("tidyverse")
## Macユーザは次の行も実行する
#theme_set(theme_gray(base_size = 10, base_family = "HiraginoSans-W3"))
library("broom")
#install.packages("coefplot", dependencies = TRUE)  # coefplotをインストール済みでない者は実行
library("coefplot")

INC <- read_csv("data/6_1_income.csv")
glimpse(INC)

## Observations: 4,299
## Variables: 5
## $ exper   <dbl> 7, 8, 8, 10, 10, 11, 12, 13, 14, 14, 15, 15, 16, 16, 16,…
## $ exper2  <dbl> 49, 64, 64, 100, 100, 121, 144, 169, 196, 196, 225, 225,…
## $ yeduc   <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,…
## $ income  <dbl> 100.00, 150.00, 150.00, 200.00, 300.00, 150.00, 200.00, …
## $ lincome <dbl> 4.605170, 5.010635, 5.010635, 5.298317, 5.703783, 5.0106…

myd <- INC %>%
  dplyr::select(exper, yeduc, income)

glimpse(myd)

## Observations: 4,299
## Variables: 3
## $ exper  <dbl> 7, 8, 8, 10, 10, 11, 12, 13, 14, 14, 15, 15, 16, 16, 16, …
## $ yeduc  <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, …
## $ income <dbl> 100.00, 150.00, 150.00, 200.00, 300.00, 150.00, 200.00, 1…

rm(INC)

回帰分析

ミンサー方程式を推定してみよう。この回帰分析における結果変数は所得（万円）の自然対数、説明変数は修学年数、就業可能年数と就業可能年数の二乗である。

まず、lm() を使って推定を行う。推定式の中で自然対数を計算するために log() を、二乗の計算を行うために、二乗した式を I() で囲む。

fit_mincer <- lm(log(income) ~ yeduc + exper + I(exper^2), data = myd)

結果を確認してみよう。

tidy(fit_mincer, detail = TRUE)

## # A tibble: 4 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  2.49     0.111         22.4 1.64e-105
## 2 yeduc        0.118    0.00706       16.6 2.31e- 60
## 3 exper        0.196    0.00749       26.2 2.75e-140
## 4 I(exper^2)  -0.00638  0.000316     -20.2 1.32e- 86

教科書 p.135 と同じ結果が得られた。結果変数が対数になっていて、修学年数 (yeduc) の係数の推定値 (coef.est) が 0.12なので、教育の収益率は12%ほどであることがわかる。また、\(p\)値が0.05より小さいので、5パーセントの有意水準で、この効果は統計的に有意である。

問： では、この効果は実質的に重要だろうか？

信頼区間も求めよう。 95%信頼区間は以下のコマンドで求めらる。

confint(fit_mincer)

##                    2.5 %       97.5 %
## (Intercept)  2.268311730  2.702692653
## yeduc        0.103704953  0.131388465
## exper        0.181482422  0.210864855
## I(exper^2)  -0.007001041 -0.005761256

点推定値と信頼区間を図示するには、coefplot::coefplot() を使う。

coefplot(fit_mincer)

私たちがあまり興味を持たない切片の値だけ他の値とかけ離れていて図が見にくいので、切片を除外する。また、lwdOuter = 1 として 95%信頼区間を示す。

coefplot(fit_mincer, intercept = FALSE, lwdOuter = 1)

これで、係数の推定値とその95%信頼区間を示すことができた。

さらに微調整を加えよう。

coefplot(fit_mincer, intercept = FALSE, lwdOuter = 1,
         xlab = "係数の推定値", ylab = "説明変数", 
         title = "回帰分析の結果：結果変数は対数年収",
         newNames = c(yeduc = "修学年数", 
                      exper = "就業可能年数", 
                      `I(exper^2)` = "就業可能年数の二乗"))

修学年数の影響とは異なり、就業可能年数の効果は表からはわかりにくい。なぜなら、就業可能年数を動かすと、就業可能年数の二乗も一緒に動いてしまうからだ。そこで、修学年数を固定し、就業可能年数と就業可能年数を動かすと、結果変数である年収の対数がどのように動くか図示してみよう。

まず、データ内の修学年数の平均値を求める。

(mean_yeduc <- mean(myd$yeduc))

## [1] 13.85997

次に、就業可能年数の最小値と最大値を求め、その範囲の値をとる長さ1,000のベクトルを作る。

exper_vec <- seq(from = min(myd$exper), to = max(myd$exper), length = 1000)

修学年数を平均値に固定し、就業可能年数と就業可能年数を動かして予測値を求める。

pred_mean <- coef(fit_mincer)[1] + 
  coef(fit_mincer)[2] * mean_yeduc +
  coef(fit_mincer)[3] * exper_vec +
  coef(fit_mincer)[4] * exper_vec^2

就業可能年数の変化に応じて年収の自然対数がどのように変化するか図示してみよう。

dd <- tibble(exper = exper_vec,
             pred_mean  = pred_mean)   # ggplot2 で使うためにデータフレームを作る
plt_1 <- ggplot(dd, aes(x = exper, y = pred_mean)) +
  geom_line() +
  labs(x = "就業可能年数", y = "対数年収の予測値", 
       title = "修学年数が平均値（13.9年）のとき")
plot(plt_1)

これで、就業年数と対数年収の間にある非線形の関係が図示できた。 対数年収はわかりにくいので、exp() で元の単位に戻そう。

plt_2 <- ggplot(dd, aes(x = exper, y = exp(pred_mean))) +
  geom_line() +
  labs(x = "就業可能年数", y = "年収の予測値（万円）", 
       title = "修学年数が平均値（13.9年）のとき")
plot(plt_2)

これで、修学年数が平均値のとき、就業可能年数と年収の間にある関係が図示できた。

修学年数が最小値のときと、最大値のときはどうだろうか。

修学年数が最小値のときは、

(min_yeduc <- min(myd$yeduc))

## [1] 9

pred_min <- coef(fit_mincer)[1] + 
  coef(fit_mincer)[2] * min_yeduc +
  coef(fit_mincer)[3] * exper_vec +
  coef(fit_mincer)[4] * exper_vec^2
dd$pred_min <- pred_min
plt_3 <- ggplot(dd, aes(x = exper, y = exp(pred_min))) +
  geom_line() +
  labs(x = "就業可能年数", y = "年収の予測値（万円）", 
       title = "修学年数が最小値（9年）のとき")
plot(plt_3)

修学年数が最大値のときは、

(max_yeduc <- max(myd$yeduc))

## [1] 18

pred_max <- coef(fit_mincer)[1] + 
  coef(fit_mincer)[2] * max_yeduc +
  coef(fit_mincer)[3] * exper_vec +
  coef(fit_mincer)[4] * exper_vec^2
dd$pred_max <- pred_max
plt_4 <- ggplot(dd, aes(x = exper, y = exp(pred_max))) +
  geom_line() +
  labs(x = "就業可能年数", y = "年収の予測値（万円）", 
       title = "修学年数が最大値（18年）のとき")
plot(plt_4)

すべてを同じ図にしたいときは、以下のようにする（やや難しい）。 まず、データを変換する。 現在のデータは、

head(dd)

## # A tibble: 6 x 4
##    exper pred_mean pred_min pred_max
##    <dbl>     <dbl>    <dbl>    <dbl>
## 1 0           4.11     3.54     4.60
## 2 0.0260      4.12     3.55     4.61
## 3 0.0521      4.12     3.55     4.61
## 4 0.0781      4.13     3.56     4.62
## 5 0.104       4.14     3.56     4.62
## 6 0.130       4.14     3.57     4.63

dim(dd)

## [1] 1000    4

のように横長 (wide) になっている。これを縦長 (long) に変換する。そのために、tidyr::gather() を使う。

dd_long <- dd %>%
  gather(key = "yeduc", value = "predicted", -exper)

どのように変換されたか確認してみよう。

head(dd_long)

## # A tibble: 6 x 3
##    exper yeduc     predicted
##    <dbl> <chr>         <dbl>
## 1 0      pred_mean      4.11
## 2 0.0260 pred_mean      4.12
## 3 0.0521 pred_mean      4.12
## 4 0.0781 pred_mean      4.13
## 5 0.104  pred_mean      4.14
## 6 0.130  pred_mean      4.14

dim(dd_long)

## [1] 3000    3

1000行4列のデータが、3000行3列になったことがわかる。これで、3つの曲線を1つの図に示せる。 以下のようにする。

plt5 <- ggplot(dd_long, aes(x = exper, y = exp(predicted), color = yeduc)) +
  geom_line() +
  labs(x = "就業可能年数", y = "年収の予測値（万円）") +
  scale_color_discrete(name = "修学年数",
                       labels = c("最大値 (18年）", "平均値（13.9年）", "最小値（9年）"))
plot(plt5)