[R] 東京都の所得階級分布から元の分布を推定する方法

概要

自分が以前書いていたものを見返していたら,
ill-identified.hatenablog.com
において, 都内の市区町村別の平均世帯収入の見積もり方が少し気になったので, もう少し厳密な方法で推定してみた.

階級の境界値とカウント数のみがわかっている, 集計加工済みのデータに対して, 仮定した元の分布のパラメータを推定する方法である.

記事を書いた当時は天下り的に, というか面倒臭かったのでよそ*1の計算手順をそのまま流用したのだが, 今回はもう少しまじめに考えてみたい.
収入の階級値がわかっていて境界値が分かってないなどの場合はまた別の方法になるが, 自分は詳しくないしそこまでカバーするのは時間がかかるので今回は書かない.
本文は約12 KB

階級分布の特徴

元のデータは, 2010 (平成20) 年の「土地統計調査確報集計都道府県編 13東京都市区町村第29表」であり, 世帯収入の階級の境界値と, それに含まれる世帯数のみがわかっている. 階級で表す場合は階級値, つまり階級ごとの平均値なども掲載することが多いが, 今回はこれがわかっておらず, 「100万円以上200万円未満」のように階級を区分する境界値だけが記載されている.

以前の記事では, 境界値が「X万円以上Y万円未満」のように上下に境界が設けられている場合は, 2つの中央値*2を階級値とし, 「1500万円以上」のように上限がない最上位の階級では「2000万円」と決めたうえで, 各階級の世帯数を掛けて平均世帯年収を計算していた. 今回はこのやり方をもう少し洗練させるというのが, 狙いになる.

とりあえず, 従来のやり方で便宜的に収入階級値ごとの世帯人数をグラフにしてみる. どの点がどの市区に対応しているかは書いてもわかりづらいので表記してない.

f:id:ill-identified:20160923145824p:plain

元のデータは構成人数別など, 世帯の特徴ごとにも集計を分けていたりするが, 今回は市区町村の全世帯の値だけを見ている. また, 土地統計調査は 2015 年版もあるのだが, こちらは世帯年収の区分が6階級に削減されている. それが以下の図.

f:id:ill-identified:20160923145850p:plain

2015年版では多くの市区では300~400万円台がモードとなるのにもかかわらず, 階級区分を削減したことで分布が潰れてしまい, 形状がわかりにくくなっている. このような場合でも, 後で紹介する方法でできないことは無いが, 計算が収束しづらくなる可能性がある.

推定に使うモデル

分布の形状を推定する方法といえば, カーネル密度推定法などがあるが, 今回のデータの特徴として, 実際の数値が階級の境界値ごとにわけられており, いわば非常に強い丸め誤差が存在するようなもので, そのまま適用しようとしてもうまくいかない*3. さらに, 境界値にも上下の値があるため, 階級に属する世帯数にどう対処すればよいか, という問題もある.

f:id:ill-identified:20160923145923p:plain

先に書いたように, 階級値がなく, かと言って「境界値の間をとった値で代用する」というのも明確な根拠に欠ける.

そこで, 今回の方法では, 「階級への分類」として捉えることで問題を解決する. つまり, 世帯収入がある分布 $F(\theta)$ に従っているとして, ある世帯の収入を $x$ とする. しかし, 今はこの $x$ は不明で,代わりに階級 $y$ に数え上げられる. $K$ 個の階級は $\mathcal{C}_{k}$ で表し, 今回は12階級あるので $K=12$ となる. つまり, 以下のようなルールになる.

$\begin{align} x\sim & F(\theta),\\ y= & \begin{cases} \mathcal{C}_{1} & \text{if }x\in]b_{0},\,b_{1}[\\ \mathcal{C}_{2} & \text{if }x\in[b_{1},\,b_{2}[\\ & \vdots\\ \mathcal{C}_{K-1} & \text{if }x\in[b_{K-2},\,b_{K-1}[\\ \mathcal{C}_{K} & \text{if }x\in[b_{K-1},b_{K}[ \end{cases}\end{align}$

ただし, $b_{i}$ は階級の境界値で, $b_{0}=-\infty$ , $b_{K}=\infty$ である. $x$ がどの階級に属するかは, 分布 $F$ と階級の境界値で決まる.

$\begin{align} p_{k}= & P\{x\in\mathcal{C}_{k}\}\\ = & \int_{b_{k-1}}^{b_{k}}f(x;\,\theta)dx\end{align}$

よって, あとは適切な目的関数を決めて評価すれば, 母数 $\theta$ が決まり分布の形が特定できる. では, どういう目的関数を使うかというと, たとえば松浦 (1993, 世帯主年齢階級別所得分布の推定) では $\chi^{2}$ 値

$\begin{align} \chi^{2}= & n\sum_{k=1}^{K}\frac{(n_{k}/(n-p_{k}))^{2}}{p_{k}}\end{align}$

を最小化することでパラメータを求めていた. $n_{k}$ は $\mathcal{C}_{k}$ に属する世帯数, $n$ は全世帯数になる.

今ひとつの方法は, McDonald & Ransom (1979, Functional Forms, Estimation Techniques and the Distribution of Income) にあるように, $K$ 個の階級のいずれに属するかを多項分布で表し, 最尤推定するというものになる. 各階級には $n_{k}$ 世帯が含まれるならば, 1世帯がその階級に含まれる確率を $p_{k}$ とすれば, 各階級ごとに何世帯が含まれるのかというのを多項分布を用いて尤度関数を以下のように表現できる:

$\begin{align} \mathcal{L}= & n!\prod_{k=1}^{K}\frac{p_{k}^{n_{k}}}{n_{k}}.\end{align}$

この尤度関数を最大化することで分布のパラメータを求める. (どうでもいいが, 自分が最初に考えた方法はこちらに近い.). R は最小カイ二乗推定をするパッケージもあるらしいが, 今回は後者でやってみる.

仮定する分布は, 対数正規分布で当てはめてみるが, ガンマ分布や対数ロジスティック分布や Singh-Maddala 分布*4などを使うこともできるようだ.

対数尤度は,

$\begin{align} \ln\mathcal{L}= & \sum_{k=1}^{K}\ln\frac{p_{k}^{n_{k}}}{n_{k}}+\ln n!\\ = & \sum_{k=1}^{K}n_{k}\ln p_{k}-\sum_{k=1}^{K}\ln n_{k}+\sum_{i=1}^{n}\ln i\end{align}$

となるので, これを元に stat4 パッケージの mle 関数で計算する*5 *6.

なお, In plnorm(b.l, meanlog = lmu, sdlog = lsigma) : 計算結果が NaN になりました という警告が出ることがあるが, これは問題ないはず. plnorm() は確率点がゼロ以下でもエラーを返さないので, 分散パラメータの最適値を探索する際にマイナスの値を指定することがあるのが原因と思われる.

2016/10/15: 解の探索中に分散パラメータがマイナスになることによる警告は, 変数を指数変換することで解決できる。この場合, 返ってくる分散パラメータの解は対数なので, 結果をプロットする際は再び指数変換しなければならない. これは数値計算では基本的なテクニックだったらしい

なおガンマ分布はスケールパラメータ $\beta$ の初期値をいい加減にするとすぐ計算エラーがでるので桁くらいは合わせたほうが良い.

対数正規分布の場合のコードは以下の通り. mle() は対数尤度を**マイナスにした**値を返す関数を与えないとダメという罠があるので使ったことがない人は注意。

gistf1826ad6a7e552823c5664416e9a3914

このコードは入力データを以下のような形式にして用いる. b.u が境界の上限値, b.l が下限値, n が世帯数である.市区町村ごとの計算なので，市区町村ごとにレコードを抽出して計算する

       region  b.u  b.l    n
1  あきる野市  100 -Inf  730
2  あきる野市  200  100 2030
3  あきる野市  300  200 3500
4  あきる野市  400  300 4140
5  あきる野市  500  400 3600
6  あきる野市  600  500 2870
7  あきる野市  700  600 2020
8  あきる野市  800  700 2010
9  あきる野市  900  800 1620
10 あきる野市 1000  900 1330
11 あきる野市 1500 1000 1600
12 あきる野市  Inf 1500  440
13     稲城市  100 -Inf 1200
14     稲城市  200  100 2310
15     稲城市  300  200 4030
16     稲城市  400  300 3840
17     稲城市  500  400 3620
18     稲城市  600  500 3260
19     稲城市  700  600 2460
20     稲城市  800  700 2420

結果と結論

推定された分布を, 新宿区と港区のみグラフにする.わかりにくいが, 港区の平均は約803万, 新宿区は約574万である.

f:id:ill-identified:20160923150206p:plain
f:id:ill-identified:20161015002800p:plain

最後に, 推定された各市区町村の平均世帯収入を計算する. 対数正規分布なので, 推定されたパラメータで, $\exp(\mu+\sigma^2/2)$ をとれば平均収入がわかる.2016/10/15: 結果表が一部おかしな数値を示していたので差し替え。

f:id:ill-identified:20161015005733p:plain

従来の方法と比較したところ, 総じて今回の方法での推定値のほうが大きくなり, 千代田区や港区などは数十万の差がある. 一方で平均値の低い区では, 差が小さいように見える. おそらく, 1500万以上世帯の階級値を恣意的に一律に決めてしまったことで, 高所得世帯の多い区ほど差が開いたのだろう. 今回試した方法では, 統計的な誤差はもちろん存在するものの, 主観の入り込まない階級の上下限値を用いて推定しているという違いがある. ここからわかることは, (1)階級値を主観的に一律に決める方法は, 高所得世帯の多い区ほど平均の誤差が大きくなるという問題がある. (2)対数正規分布のような歪んだ分布の平均を目分量で判断するのは難しい,ということである. ただし, 従来の方法を否定したからといって, 即座に対数正規分布のあてはめが正しいこと意味するわけではない. 例えば今回はやらなかったガンマ分布や対数ロジスティック分布なども検討し, 厳密に見ていく必要があるだろう (ただしもっと細かい話は現状関心がないので私はやらない. 方針は示したのでやる気のある方がんばってください.).

参考文献

McDonald, J. B., & Ransom, M. R. (1979). Functional Forms, Estimation Techniques and the Distribution of Income. Econometrica, 47(6), 1513. doi:https://doi.org/10.2307/1914015

松浦和幸 (1993). 世帯主年齢階級別所得分布の推定. 産業経営, 19, 193–218. Retrieved from https://dspace.wul.waseda.ac.jp/dspace/handle/2065/37120

*1:http://tmaita77.blogspot.jp/2012/12/blog-post.htmlhttp://tmaita77.blogspot.jp/2012/12/blog-post.html

*2:実質, 平均値と同じ

*3:カーネル密度推定をうまく応用する方法を思いついてないだけかもしれない.

*4:初めて知ったのだが, Type-XII Burr 分布とも呼ばれ, 一般化対数ロジスティック分布の特殊形らしい.

*5:ただし, 人数nが多いと尤度が桁溢れすることがあるため, パラメータに依存しない後ろの2つの項は足さずに第1項だけを最大化してもよい.

*6:カイ二乗推定もこのパッケージで計算できそうではある.