度数・頻度データの分析
ポアッソン指数型回帰モデル
度数表(表1)における行r、列cにおける度数(頻度)データの値がであるとする。
表1 行r、列cにおける度数データ |
||||
|
・・・ |
列c |
・・・ |
|
|
|
|
|
|
行r |
・・・ |
|
・・・ |
|
|
|
|
|
|
|
・・・ |
|
・・・ |
|
このとき、行r、列cの期待度数をとおくと、行と列が独立であることは次式で表せる(Kruschke, 2011)。
これより、次式を得る。
期待度数をで、をで、をで、をで表すと
となり、式(1)の表す独立性は、式(2)における列と行の主効果の和で表される。これに対して、交互作用がある場合は
となる。
ポアッソン指数型回帰モデル(Poisson exponential regression model)では、観測度数は、パラメータのポワッソン分布から得られたとする。すなわち、
である。ここで、がパラメータのポワッソン分布に従うとき
である。ポアッソン分布の描画プログラムは<ポアッソン分布の描画>に用意している。
表1のデータに対する尤度関数は
事後分布は
で与えられる。ここで、は事前分布である。
分析例
表2は平成7年国勢調査による横浜市のデータである(2013年にダウンロードしたhttp://www.city.yokohama.lg.jp/ex/stat/jinko/sangyo-betu/h07.pdfから抜粋)。
表2 横浜市平成7年国勢調査 |
|||
|
第1次産業 |
第2次産業 |
第3次産業 |
男 |
6,547 |
390,973 |
688,777 |
女 |
3,877 |
105,220 |
483,473 |
表2のデータを式(2)あるいは式(3)のようなポアッソン指数型回帰モデルに従って分析するために表3のようにパラメータを設定する。
表3 のモデル |
|||
|
第1次産業 |
第2次産業 |
第3次産業 |
男 |
|
|
|
女 |
|
|
|
〜 が主効果を表し、と が交互作用を表す。表3のモデルをダミー変数を用いて表すと
となる。このとき表2のデータは表4のように表せる。
表4 データとダミー変数 |
||||||
|
|
|
|
|
|
|
6547 |
1 |
0 |
0 |
0 |
0 |
0 |
390973 |
1 |
1 |
0 |
0 |
0 |
0 |
688777 |
1 |
0 |
1 |
0 |
0 |
0 |
3877 |
1 |
0 |
0 |
1 |
0 |
0 |
105220 |
1 |
1 |
0 |
1 |
1 |
0 |
483473 |
1 |
o |
1 |
1 |
0 |
1 |
ダミー変数については、<岡本安晴「データ分析のための統計学入門」おうふう、2009、第4.3節「重回帰分析と分散分析」>などを参照されたい。
プログラムPoissonExpRegは、表4の形式のデータに対してパラメータおよびその線形結合の事後分布を求めるものである。
ソースコードファイルなどは圧縮ファイルPoissonExpReg.zipとしてまとめたので、ファイル名PoissonExpReg.zipをクリックしてダウロード解凍すれば、Visual Studio 2012で開くことができる。
実行形式ファイルPoissonExpReg.exeはファイル名PoissonExpReg.exeをクリックしてダウンロード実行することができる。実行時に「何々.dllがインストールされていません」というようなエラーがでるときは、
http://www.microsoft.com/visualstudio/jpn/downloads
に無料の開発ツールとして用意されているVisual Studio Express 2012 for Windows Desktopをダウンロード・インストールすれば実行に必要なライブラリもインストールされ、実行形式プログラムPoissonExpReg.exeが実行できるようになる。
表4のデータをプログラムの入力ファイルとして用意するときは、次のようにする。
===================== 次の行から ====================
Sample Data
サンプル データ
/ 独立変数の数
6
X0 = 1: ベース値
X1: 第2次産業
X2: 第3次産業
X3: 女性
X4: X1 * X3
X5: X2 * X3
/
6547 1 0 0 0 0 0
390973 1 1 0 0 0 0
688777 1 0 1 0 0 0
3877 1 0 0 1 0 0
105220 1 1 0 1 1 0
483473 1 0 1 1 0 1
/
/
0 -1 1 0 0 0
/
0 -1 1 0 -1 1
/
0 0 0 1 1 0
/
0 0 0 1 0 1
/E
================= 上の行まで =====================
データは、/ では始まる行の次の行から、ひとまとまりごとに書かれている。
1番目のスラッシュ / で始まる行に続いて、独立変数Xjの個数がまず書かれている(上のデータの場合は 6 が書かれている)。続いて、各独立変数の説明が書かれている。
2番目のスラッシュで / で始まる行と3番目のスラッシュ / で始まる行に挟まれて表4のデータが書かれている。
4番目以降のスラッシュ / で始まる行に続いて書かれているのはウェイトであり、各ウェイトの組に対して
が算出される。
/E で始まる行のところで、プログラムの処理は終了する。
プログラムは、実行開始時に図1のようなフォームが現れる。
図1
「Go」ボタンのクリックで入力データファイルの指定を求めるダイアログウィンドウが表示される(図2)。
図2
入力データファイル名(図2の場合は、DataYokohama.txt)を指定して「開く」ボタンをクリックすると、次に出力ファイル名の設定を求めるダイアログウィンドウが表示される(図3)。
図3
図3の場合はtemp.txtが設定されている。既存のファイル名を設定すると出力ファイルが上書きされる。
上のデータファイルを入力ファイルとしたときの出力ファイルは次のようになる。
=======================================
Data
File...D:\yasuharu\www\LaCoocan\booksetc\pm2010\PoissonRegression\PoissonExpReg\PoissonExpReg\DataYokohama.txt
Sample
Data
サンプル データ
/ 独立変数の数
X[1]
--> X0
= 1: ベース値
X[2]
--> X1:
第2次産業
X[3]
--> X2:
第3次産業
X[4]
--> X3:
女性
X[5]
--> X4:
X1 * X3
X[6]
--> X5:
X2 * X3
/
Sample
Data
サンプル データ
/ 独立変数の数
6
X0 =
1: ベース値
X1: 第2次産業
X2: 第3次産業
X3: 女性
X4:
X1 * X3
X5:
X2 * X3
/
Data...
0: 6547 1 0 0 0 0 0
1: 390973 1 1 0 0 0 0
2: 688777 1 0 1 0 0 0
3: 3877 1 0 0 1 0 0
4: 105220 1 1 0 1 1 0
5: 483473 1 0 1 1 0 1
tempB[0] =
8.78676220844108
tempB[1] =
4.08963157443366
tempB[2] =
4.65591063169172
tempB[3] =
-0.523945270770165
tempB[4] =
-0.788639836819782
tempB[5] =
0.170023180140749
tempSgmB =
3.45843292678339
Parameter
Mean
LB
Med
UB
B[0]:
8.78393 8.76216
8.78337
8.80750
B[1]:
4.09244
4.06875
4.09300
4.11471
B[2]:
4.65875
4.63514
4.65945
4.68055
B[3]:
-0.51761 -0.55220 -0.51698 -0.48046
B[4]:
-0.79493 -0.83244
-0.79549 -0.75973
B[5]:
0.16368
0.12653
0.16311
0.19894
Weights...
0 -1 1 0 0 0
Mean =
0.56631 LB =
0.56244 Median =
0.56631 UB =
0.57020
Weights...
0 -1 1 0 -1 1
Mean =
1.52491 LB =
1.51827 Median =
1.52488 UB =
1.53159
Weights...
0 0 0 1 1 0
Mean =
-1.31254 LB =
-1.31933 Median =
-1.31252 UB =
-1.30578
Weights...
0 0 0 1 0 1
Mean =
-0.35393 LB =
-0.35766 Median =
-0.35391 UB =
-0.35025
=======================================
「Parameter」で始まる行に続く、B[j]の行の値は、パラメータの事後分布の平均(Mean)、95%確信区間(Credible Interval)の左端の値(LB)、中央値(Med)、95%確信区間の右端の値(UB)である。いずれのパラメータの95%確信区間にも0が含まれていないので、パラメータはすべて5%水準で有意に0ではない。確信区間については<岡本安晴「計量心理学」培風館、2006、付録H「ベイズ的アプローチ」>などを参照されたい。
ウェイトは、例えば最初の
0 −1 1 0 0 0
は、表3のモデルにおいて「男、第3次産業」と「男、第2次産業」との差
を表すもので、この95%確信区間が0を含まないことから、「男、第3次産業」と「男、第2次産業」の差は5%水準で有意であることがわかる。
上のデータにおいて全体のデータを1000分の1の比率で小さくしたデータを次のように作成してみた。
=========================================
Sample Data
サンプル データ
/ 独立変数の数
6
X0 = 1: ベース値
X1: 第2次産業
X2: 第3次産業
X3: 女性
X4: X1 * X3
X5: X2 * X3
/
7 1 0 0 0 0 0
391 1 1 0 0 0 0
689 1 0 1 0 0 0
4 1 0 0 1 0 0
105 1 1 0 1 1 0
483 1 0 1 1 0 1
/
/
0 -1 1 0 0 0
/
0 -1 1 0 -1 1
/
0 0 0 1 1 0
/
0 0 0 1 0 1
/E
=========================================
上のデータに対する出力ファイルは以下のようになる。
=========================================
Data
File...D:\yasuharu\www\LaCoocan\booksetc\pm2010\PoissonRegression\PoissonExpReg\PoissonExpReg\DataYokohama_small.txt
Sample
Data
サンプル データ
/ 独立変数の数
X[1]
--> X0
= 1: ベース値
X[2]
--> X1:
第2次産業
X[3]
--> X2:
第3次産業
X[4]
--> X3:
女性
X[5]
--> X4:
X1 * X3
X[6]
--> X5:
X2 * X3
/
Sample
Data
サンプル データ
/ 独立変数の数
6
X0 =
1: ベース値
X1: 第2次産業
X2: 第3次産業
X3: 女性
X4:
X1 * X3
X5:
X2 * X3
/
Data...
0:
7 1 0 0 0 0 0
1:
391 1 1 0 0 0 0
2:
689 1 0 1 0 0 0
3:
4 1 0 0 1 0 0
4:
105 1 1 0 1 1 0
5:
483 1 0 1 1 0 1
tempB[0] =
1.94591014905531
tempB[1] =
4.02279741093005
tempB[2] =
4.58933112195835
tempB[3] =
-0.559615787935424
tempB[4] =
-0.75513142189242
tempB[5] =
0.204391170574336
tempSgmB =
2.12446633652249
Parameter
Mean
LB
Med
UB
B[0]:
1.72959
0.93683
1.73703
2.43706
B[1]:
4.23761
3.51942
4.22943
5.04491
B[2]:
4.80549
4.09365
4.79504
5.59769
B[3]:
-0.37691 -1.66294
-0.33019
0.69332
B[4]:
-0.94079 -2.03140 -0.98796
0.36313
B[5]:
0.02082 -1.06404
-0.02437
1.31282
Weights...
0 -1 1 0 0 0
Mean =
0.56788 LB =
0.44349 Median =
0.56892 UB =
0.69069
Weights...
0 -1 1 0 -1 1
Mean =
1.52949 LB =
1.32379 Median =
1.52723 UB =
1.74702
Weights...
0 0 0 1 1 0
Mean =
-1.31770 LB =
-1.53906 Median =
-1.31568 UB =
-1.10579
Weights...
0 0 0 1 0 1
Mean =
-0.35609 LB =
-0.47318 Median =
-0.35643 UB =
-0.23655
=========================================
元のデータの場合では、95%確信区間が0を含まなくても、度数を1000分の1にしたデータでは95%確信区間が0を含むもの(5%水準で有意に0でないといえないもの)も出てくることがわかる。すなわち、データが少ないと統計的分析結果において有意差が出難くなる。
上の結果の場合、産業の主効果との95%確信区間は0を含まないという意味で有意だが、性の主効果は95%確信区間が0を含むという意味で有意でない。男7人と女4人の差は非有意ということである。
ウェイトによる分析結果は、以下のように有意差が認められる。
男・第3次―男・第2次: ウェイト(0、 −1、 1、 0、 0、 0)
[LB, UB]=[0.443, 0.691]で有意。
女・第3次―女・第2次: ウェイト(0、 −1、 1、 0、 −1、 1)
[LB, UB]=[1.324, 1.747]で有意。
女・第2次―男・第2次: ウェイト(0、 0、 0、 1、 1、 0)
[LB, UB]=[-1.539, -1.106]で有意。
女・第3次―男・第3次: ウェイト(0、 0、 0、 1、 0、 1)
[LB, UB]=[-0.473, -0.237]で有意。
ここで、[LB, UB]は95%確信区間を表す。
ベイズ的データ分析の場合、データに含まれる情報は事後確率分布で表される。この事後確率分布はパラメータの同時確率分布であり、データに対応して決まるが、同時確率分布の表す情報は、上の場合のようにパラメータの関数としての重み付和の確率分布として簡単に取り出すことができる。
参考文献
Kruschke, J. K. (2011). Doing Bayesian Data Analysis: A tutorial with R and BUGS. Elsevier.