Up

度数・頻度データの分析

ポアッソン指数型回帰モデル

 

度数表(表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.

 

 

Up