廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

StanとRで最低賃金と失業率の関係を調べる

今日の川柳

下記の内容について「先行研究も調べずにがさつな分析で結論を出すのはよくない」(引用は不正確)というようなコメントを頂戴し、そりゃそうだとおもったので最低賃金と雇用の関係について勉強になりそうな文献へのリンクをいくつか貼ります。

以下の文章はあまり真に受けず他山の石としてしていただければ幸いです。

記述統計

ちょっと古い話題になるのかもしれないけど、最低賃金引き上げは失業率を上昇させるか? という議論がある。

最低賃金引き上げは失業率を上昇させるか? - himaginaryの日記

最低賃金を上がると人を雇いづらくなるから、失業率が上がる」と言われたら、そういうこともあるかもしれないなあと思う。

むずかしい理屈はよくわからないから、データからわかる範囲で調べよう。

都道府県別の最低賃金のデータは
地域別最低賃金の全国一覧 |厚生労働省
にある。

PDFファイルになってて扱いにくかったので手作業でCSVにした。
最低賃金.csv · GitHub
手作業なのでどこかミスってるかもしれない。ミスってたらすみません。

労働力調査の結果は
統計局ホームページ/<参考>労働力調査(基本集計)都道府県別結果
にある。こっちはエクセルファイルになってる。

データを成形して

library(tidyverse)
library(readxl)

saitei <-read_csv("最低賃金.csv")
rodo <-read_xls("lt02y.xls",skip = 7)
shitsugyo <- read_xls("lt04y.xls",skip=7)

shitsugyo2 <-shitsugyo %>% 
  mutate(year=as.integer(substr(X__2,2,5))) %>% 
  select(-X__1,-X__2) %>% 
  gather(prefecture,shitsugyo,-year) %>% 
  mutate(shitsugyo=as.integer(shitsugyo)) %>% 
  dplyr::filter(!is.na(shitsugyo)) %>% 
  mutate(prefecture=gsub("府|県","",prefecture))%>% 
  mutate(prefecture=gsub("東京都","東京",prefecture)) 

rodo2 <- rodo %>% 
  mutate(year=as.integer(substr(X__2,2,5))) %>% 
  select(-X__1,-X__2) %>% 
  gather(prefecture,rodo,-year) %>% 
  mutate(rodo=as.integer(rodo)) %>% 
  dplyr::filter(!is.na(rodo)) %>% 
  mutate(prefecture=gsub("府|県","",prefecture)) %>% 
  mutate(prefecture=gsub("東京都","東京",prefecture)) 

saitei2<-saitei %>% 
  gather(prefecture,chingin,-year)

joint_df <-left_join(saitei2,shitsugyo2,by=c("year","prefecture")) %>% 
  left_join(rodo2,by=c("year","prefecture")) %>% 
  mutate(shitsugyo_ritsu=shitsugyo/rodo) %>% 
  dplyr::filter(!is.na(shitsugyo_ritsu)) %>% 
  group_by(year) %>% 
  mutate(med_chingin=median(chingin),prefecture=factor(prefecture))

いくつかグラフを書いてみる。

最初は散布図。

theme_set(theme_bw(12,"HiraKakuProN-W3"))
p_point <-ggplot(joint_df,aes(y=shitsugyo_ritsu,x=chingin,colour=prefecture))+
  geom_point()+
  theme(legend.position = "none")+
  xlab("最低賃金")+ylab("完全失業者数/労働力人口")
print(p_point)

f:id:abrahamcow:20180224023436p:plain

都道府県で色分けしている。
これを見た感じ、最低賃金が高いときはむしろ失業率が低いんじゃないかと思える。

相関係数の検定(無相関検定)をやったところ、p値はめっちゃ小さくなった。

> cor.test(joint_df$chingin,joint_df$shitsugyo_ritsu)

	Pearson's product-moment correlation

data:  joint_df$chingin and joint_df$shitsugyo_ritsu
t = -9.3594, df = 703, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3969498 -0.2655545
sample estimates:
       cor 
-0.3328669 

これでおしまい、最低賃金引き上げは失業率を上昇させません、ということにはしない。

このやり方は都道府県の情報や時系列の情報をつぶしている。

対応のある標本の分析をわざわざStanでやる - 廿TT で触れたように、「個体差や場所差が識別できてしまうようなデータのとりかたをしている」ときは、個体差や場所差を考慮して分析するとしないとで結果が大きく変わることがよくある。

時系列の情報を見るために箱ひげ図を書いてみよう。

p_box <-ggplot(joint_df,aes(x=year,y=shitsugyo_ritsu,group=year,fill=med_chingin))+
  geom_boxplot()+
  xlab("年")+ylab("完全失業者数/労働力人口")+
  labs(fill="最低賃金(中央値)")
print(p_box)

f:id:abrahamcow:20180224024353p:plain

色はその年の全国の最低賃金の中央値に対応している。2002年から2006年までは最低賃金が低いままで失業率が下がっている。2009年以降は失業率の減少と最低賃金の引き上げが連動しているように見える。2016年は最低賃金が高いのに失業率は低く抑えられていて良い感じだ。

折れ線グラフも書こう。

p_line <-ggplot(joint_df,aes(x=year,y=shitsugyo_ritsu,group=prefecture,colour=chingin))+
  geom_line()+
  geom_point()+
  xlab("年")+ylab("完全失業者数/労働力人口")+
  labs(colour="最低賃金")
print(p_line)

f:id:abrahamcow:20180226082529p:plain

一本一本の折れ線が都道府県に対応している。最低賃金高い県もあれば低い県もありますねという感じでよくわからない。

モデルとそのパラメータの推定結果

目的変数 y_i を完全失業者数として、説明変数 x_i最低賃金にする。都道府県  pref_i や時間  year_i の影響も入れる。


y_i \sim {\rm Binomial}({\rm logit}^{-1} (\alpha_{year_i}+\beta x_i +w_{pref_i}),n_i)
w_j \sim {\rm Normal}(0,\sigma^2)
\alpha_j \sim {\rm Normal}(\alpha_{j-1},\tau^2)
ここで n_i労働力人口。特に記述がないパラメータに関してはフラットプライヤを仮定している。

このモデルをStanで書くとこうなった。

//mod_binom.stan
data{
  int<lower=1> m;
  int<lower=0> year[m];
  int<lower=1> pref[m];
  int<lower=0> x[m];
  int<lower=0> n[m];
  int<lower=0> y[m];
}
transformed data{
  int<lower=1> l;
  l = max(year);
}
parameters{
  real w[max(pref)];
  real alpha[l];
  real beta;
  real<lower=0> sigma;
  real<lower=0> tau;
}
model{
  w ~ normal(0,sigma);
  alpha[2:l] ~ normal(alpha[1:(l-1)],tau);
  for(i in 1:m){
   y[i] ~ binomial_logit(n[i],alpha[year[i]]+beta*x[i]+w[pref[i]]); 
  }
}
library(rstan)
rstan_options(auto_write = TRUE)
mod_binom <-stan_model("mod_binom.stan")
datlist <- list(m=nrow(joint_df),year=joint_df$year-min(joint_df$year)+1,
                pref=as.integer(joint_df$prefecture),
                n=joint_df$rodo,y=joint_df$shitsugyo,x=joint_df$chingin)
fit_binom <- sampling(mod_binom,datlist,cores=parallel::detectCores(),
                      control=list(max_treedepth=12))

無事収束した模様です。

f:id:abrahamcow:20180224030425p:plain

> all(summary(fit_binom)$summary[,"Rhat"]<1.1)
[1] TRUE

さっそく最低賃金にかかる係数 \beta の事後分布を見てみよう。青いエラーバーは95%信用区間を示している。

ex_binom <- rstan::extract(fit_binom)
p_beta <-ggplot(,aes(x=ex_binom$beta))+
  geom_histogram(bins = nclass.FD(ex_binom$beta),fill="white",colour="black")+
  geom_errorbarh(aes(y=10,
                     xmin=quantile(ex_binom$beta,0.025),
                     xmax=quantile(ex_binom$beta,0.975)),height=10,
                 colour="royalblue",size=1)+
  xlab(expression(beta))
print(p_beta)

f:id:abrahamcow:20180226233421p:plain

ほぼ0だ。関係が「ない」ことをいうのは結構むずかしいんだけど、これだけ0に近いと最低賃金の引き上げが失業率を上げることはまずない、と言っちゃっていいと思う。

ただし、マクロには影響なくても統計データにのぼってこないような個々人への影響はあり得るのかもしれない。

せっかくなので他のパラメータの推定結果もざっと見る。

alpha_df <-t(apply(ex_binom$alpha,2,quantile,prob=c(0.025,0.975))) %>% 
  as_data_frame() %>% 
  set_names(c("lower","upper")) %>% 
  mutate(mean=apply(ex_binom$alpha,2,mean),year=2002:2016)

p_alpha <-ggplot(alpha_df,aes(x=year,y=mean,ymin=lower,ymax=upper))+
  geom_point()+
  geom_line()+
  geom_ribbon(alpha=0.3)+
  ylab(expression(alpha))
print(p_alpha)

f:id:abrahamcow:20180224031426p:plain

これは \alpha_j の事後期待値と95%信用区間\alpha_j がぐっと持ち上がってる2009年とかは、なにがあったのか調べておいたほうが良さそう。

都道府県の効果 w_k については、オッズ比とその95%信用区間を見よう。

w_df <-t(apply(exp(ex_binom$w),2,quantile,prob=c(0.025,0.975))) %>% 
  as_data_frame() %>% 
  set_names(c("lower","upper")) %>% 
  mutate(mean=apply(exp(ex_binom$w),2,mean),prefecture=sort(unique(joint_df$prefecture)))

p_w <-ggplot(w_df,aes(x=reorder(prefecture,mean),y=mean,ymin=lower,ymax=upper))+
  geom_pointrange()+
  coord_flip()+
  xlab("")+ylab("w")
print(p_w)

f:id:abrahamcow:20180224031822p:plain

失業率が高くなりやすい県のランキング。その原因も考えたほうが良さそうだけど、むずかしい。

参考にしたものなど

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3