• 検索結果がありません。

1 15 R Part : website:

N/A
N/A
Protected

Academic year: 2021

シェア "1 15 R Part : website:"

Copied!
22
0
0

読み込み中.... (全文を見る)

全文

(1)

経済情報処理 講義ノート

15

回 

R

による統計分析

Part 4

2017

7

24

日(月)

4

担当教員: 唐渡 広志 website: http://www3.u-toyama.ac.jp/kkarato/ email: kkarato@eco.u-toyama.ac.jp

目次

1 データの読み込み 2 2 回帰分析 3 2.1 散布図と回帰直線. . . 3 2.2 最小2乗推定 . . . 4 2.3 多重回帰モデル . . . 6 3 一般化線形モデル 10 3.1 ロジット・リンク. . . 10 3.2 最尤法 . . . 12 3.3 一般化線形モデルの評価 . . . 13 3.4 ロジット・モデルの推定 . . . 13 3.5 プロビット・モデル . . . 16 3.6 順序選択モデル . . . 17

(2)

1 データの読み込み 2

1

データの読み込み

データやスクリプト(プログラム)を読み書きする場合は,あらかじめ作業フォルダ(ディレク トリ)を決めておくとよい。Rを起動したら,ファイル・メニューのディレクトリの変更から,自 分の作業ディレクトリを指定することができる。指定したフォルダにデータやスクリプト(プログ ラム)を置いておくと,ファイルの読み込みや書き出しを行うときに便利である。 まず,以下の準備を行う。 1. ファイル・メニューのディレクトリの変更から,自分の作業ディレクトリを指定する 2. ファイル・メニューの新しいスクリプトから,Rに内蔵されているエディタを起動する(「無 題 Rエディタ」と表示される)。エディタにはソースコードを記述していく。 3. ファイル・メニューの保存で作業ディレクトリに名前をつけて保存しておく。拡張子はrと する。 4. data 11.csvをダウンロードして,自分の作業ディレクトリにファイルを移動しておく。 準備ができたら,次のソースコードをエディタに記述する。   data<- read.csv("data_11.csv")  

ここで,read.csvはCSV(Comma Separated Value,カンマ区切り)データを読み込むためのコ マンドであり,括弧の中に読み込むファイル名を(拡張子もつけて)ダブル・クォーテーションで 括っている。読み込んだものをdataという名前を付けて記憶している。 書き終わったら,右クリックをして「カーソル行または選択中のRコードを実行」を選択してコ マンドを実行する。または,ショートカットのctrl + Rでも同じく実行できる。エラーが出てい なければ,続いて以下のソースコードを書く。   attach(data) obs<- length(data[,1]) str(data) # データの集約 summary(data)   上記の4行を書いたら,4行分をまとめて選択して実行する。ここで, • attach はデータの変数へのアクセスを簡単にするコマンド, • length(data[,1])は(1列目の)データの長さ(観測値の数)を測るコマンド, • strはdataの中にあるオブジェクト(つまり変数)の内容を簡易表示するコマンド, • # の行はRでは無視されるので,コメント文やメモを書いておきたい場合に利用する • summaryはdataの簡易記述統計(四分位数や平均など)を示すコマンドである。 str(data)の結果から,観測値の数は2162,変数の数は44であることが示されている。

(3)

2

回帰分析

2.1

散布図と回帰直線

直線: curve, lines, abline

x-y座標平面に直線を描くにはいくつかの方法がある。curveコマンドを利用して直線 y = 3− 2x を描くには,   curve(3-2*x)   と記述する。linesコマンドは2点の座標値を結ぶ直線を描く。   lines(c(0.2,0.8),c(2,2.5))   これは,第1 の座標を (0.2, 2),第 2 の座標を (0.8, 2.5)としたときの 2 点を結んだ直線にな る。ablineコマンドは直線の切片と傾きを指定して描く方法である。例えば,切片 1,傾き2 (y = 1 + 2x)の直線の場合   abline(1,2)  

と な る 。ま た ,abline(v=0.5) と す る と 横 軸 x = 0.5 の 位 置 に 垂 線 (vertical line) を ,

abline(h=1.5)とすると縦軸y = 1.5の位置に水平線(horizontal line)を描くことができる。

回帰直線の描画 2つのデータの散布図を作成して,最小2乗法で求めた回帰直線を当てはめる場合にもabline を利用する。回答者本人の年齢[18]age.sfと配偶者の年齢[25]age.spの場合,次のように書く。   plot(age.sf, age.sp) reg<- lm(age.sf~age.sp) abline(reg)  

ここで,lmは直線モデル当てはめのコマンド(fitting linear models)であり,最小2乗推定値を 計算する。 age.sf~age.spformula とよばれるオブジェクトで,age.sf = α + β age.sp

という線形回帰直線を表現するための記法である。lmコマンドで計算した切片や傾きの情報を

regという名前で保存しておいて,これをablineに渡すことで散布図に回帰直線を描画してい

る。abline(lm(age.sf age.sp))と記述しても同じ結果が得られる。lmコマンドで計算した内

(4)

2 回帰分析 4

 

reg

Call:

lm(formula = age.sf ~ age.sp)

Coefficients: (Intercept) age.sp 2.131 0.959  

2.2

最小

2

乗推定

summaryによる分析結果の表示 lmコマンドで計算した内容の詳細を表示するにはsummaryコマンドを使う。summary(reg)ま たは summary(lm(age.sf~age.sp)) と入力する。   summary(reg) Call:

lm(formula = age.sf ~ age.sp)

Residuals:

Min 1Q Median 3Q Max -20.5875 -2.6887 -0.3003 2.7407 26.4690

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 2.130918 0.393006 5.422 6.55e-08 *** age.sp 0.958977 0.006999 137.011 < 2e-16 ***

---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.581 on 2160 degrees of freedom Multiple R-squared: 0.8968, Adjusted R-squared: 0.8968 F-statistic: 1.877e+04 on 1 and 2160 DF, p-value: < 2.2e-16

 

推定パラメータ(係数)をα, β,誤差項をui としたときの線形回帰モデル

age.sf = α + β age.sp + ui

の推定結果について,Residualsは最小2乗残差の最小,最大,四分位数の記述統計を示している。

Coefficientsは係数 α, β の推定値(Estimate),標準誤差(Std. Error),t値(t value),p

(5)

き(age.spの係数推定値)βˆである。アスタリスク * は有意性の判断基準を示す記号Signif. codesであり, • *** ならば有意水準0.1%で有意, • ** ならば有意水準1%で有意, • * ならば有意水準5%で有意, • . ならば有意水準10%で有意 となる。その他の情報は次のようになっている。

• Residual standard error 回帰の標準誤差σ =ˆ ∑uˆ2i/自由度 とその自由度degrees of freedom

• Multiple R-squared (重)決定係数R2

• Adjusted R-squared 自由度調整済み決定係数 adj.R2

• F-statistic 説明変数の係数がすべて0であるという帰無仮説(この場合はH0 : β = 0) に対するF検定統計量と分子・分母の自由度DF,およびp値 summaryの成分 summary コ マ ン ド で 示 さ れ る 計 算 結 果 は ,い く つ か の 成 分 を ま と め て 表 示 し て い る 。 names(lm(formula))と入力すると全ての成分を表示できる。   names(reg)

[1] "coefficients" "residuals" "effects" "rank"

[5] "fitted.values" "assign" "qr" "df.residual" [9] "xlevels" "call" "terms" "model"

  個々の成分を見る場合には,summaryコマンドの後ろに成分の番号やValueを付加する。例えば, 回帰の標準誤差(残差分散の平方根)を見るにはsummary(reg)[6]あるいはsummary(reg)$sig と記述する(summary(reg)$sigmaと記述してもよい)。代表的なものを表1に示した。 表1 summary(lm(formula))の代表的な成分 成分の番号 Value 出力内容 [3] $res 残差 [4] $coef 係数推定値,標準誤差,t値,p値 [6] $sig 回帰の標準誤差 [8] $r.sq 決定係数 [9] $adj.r.sq 自由度調整済み決定係数 [10] $fstat F検定統計量,分子の自由度,分母の自由度 [11] $cov 推定値の分散・共分散

(6)

2 回帰分析 6

2.3

多重回帰モデル

ダミー変数

就労時間[21] hour.sfを説明するために,年齢[18] age.sfと性別[5] sexを要因とする多重回 帰モデルを考える。性別のような名義尺度は次のようにしてダミー変数に変換する。   D<- ifelse(sex==2,1,0)   sex==2は女性であることを示しているので,Dは女性ダミー変数になる。 データの条件指定 なお,就労時間は仕事をしていない場合観察できない(非該当)データであるので,hour.sf!=888 を条件にする(非該当データを除去する)。回帰分析に利用するデータセットを次のように定義 する。   y<- hour.sf[hour.sf!=888] x<- age.sf[hour.sf!=888] D<- D[hour.sf!=888]   推定モデルは y = β1+ β2x + β3D + u (1) なので,次のようにして推定する。   reg<- lm(y~x + D) summary(reg)   被説明変数をy,説明変数をx2, x3, x4とするときの多重回帰モデルは,」 y~x2+x3+x4 のようにformulaを記述していく(回帰係数を書く必要がない)。 subsetの使用 あらかじめ,条件指定したデータセットで分析をする場合には,subsetコマンドでデータを抽 出しておく。例えば,都道府県が富山県 pref==16 の就労時間[21] hour.sf データを抽出する には   subset(hour.sf,pref==16)  

と書く。また,富山県pref==16 の男性sex==1就労時間[21] hour.sfデータを抽出するには

 

subset(hour.sf,pref==16 & sex==1)

 

(7)

data全体を抽出するには,

 

subset(data,pref==16 & sex==1)

  と書く。就労時間が非該当のデータhour.sf==888を除いて,data全体を抽出するには,   subset(data,hour.sf!=888)   と書く。 data.frameの使用とfactor によるダミー変数

csvファイルを data<- read.csv("data_11.csv")のように読み込むとき, data のこと

を data.frame (データ・フレーム)とよぶ。元のデータ・フレームから別のデータ・フレーム

を作成するには,subsetコマンドでデータを条件抽出しておき,data.frame関数で別の名前

data2を付けておく。

 

data2<- data.frame(subset(data, hour.sf!=888))

 

(1)を推定するには

 

reg<- lm(hour.sf~age.sf + factor(sex), data=data2) summary(reg)   のように書く。lmコマンドで最小2乗推定をする際に,どのデータ・フレームを使うのかというこ とをdata=data2で指定している。factorは量的な尺度を名義尺度に変換するときに使う命令で ある。元々のデータsexは量的な尺度としての1や2という数値であるが,factorによって1や 2という数値を分類上の番号として定義し直している。lmコマンドのformulaにおいて,factor を使用すると自動的に名義尺度に対するダミー変数として扱われる。 推定結果(係数推定値)は次のようになっている。   Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 63.65005 1.76697 36.02 <2e-16 *** age.sf -0.35703 0.03288 -10.86 <2e-16 *** factor(sex)2 -15.95804 0.80737 -19.77 <2e-16 ***

 

推定結果は lm(y~ x + D) と全く同じものになる。factor(sex)2は名義尺度のうち sex==2に

対する女性ダミー変数を示している。仮にこのモデルが正しいとすれば,被説明変数は(週当たり) 就労時間hour.sfであるから,(他の条件が等しければ)年齢が1歳上がると就労時間は0.357時

(8)

2 回帰分析 8

2乗項

年齢と就労時間の関係を次の2次式で表現する場合がある。

hour.sf = β1age.sf + β2age.sf2+ β3factor(sex)2 + u (2)

年齢の2乗項を formula に含めるには I関数(Inhibit Interpretation/Conversion of Objects)

を利用してI(age.sf^2) を追加する。

 

reg<- lm(hour.sf~age.sf + I(age.sf^2) + factor(sex), data=data2) summary(reg)

 

(2)式の推定結果は次のようになる

 

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 27.605957 5.633785 4.900 1.07e-06 *** age.sf 1.169697 0.229275 5.102 3.84e-07 *** I(age.sf^2) -0.015166 0.002255 -6.726 2.56e-11 *** factor(sex)2 -16.279732 0.795979 -20.452 < 2e-16 ***   2乗項を含むときの回帰式の意味を考えるために,就業年数と年齢の関係性だけに注目してみよ う。summary(reg)$coefより推定値や標準誤差などを表示できるので,β2の推定値をb2,β3の 推定値をb3 とすると,   b2<- summary(reg)$coef[2,1] b3<- summary(reg)$coef[3,1]   より,b2 = 1.169697とb3 = -0.015166を保存できる。就業年数の理論値 Yˆ は概ね次の形を しているといえる。 ˆ Y = 1.170x− 0.015x2+ extras (3) ここで,xは年齢,extrasは固定されたその他の項を示している。これを図示するには次のように 書く。   curve(b2*x+b3*x^2, 20,80, axes=F) axis(1)   ただし,axes=Fは軸を表示させないオプション,axis(1)は横軸だけを表示させるコマンドであ る。この結果を見ると,30歳代後半が就労時間のピークになっていることがわかる。実際に,就 労時間が最大となる年齢は,(3)をxについて微分して0とおくことで,x = 38.56歳であること がわかる。Rでは次のように計算したらよい。   -b2/(2*b3) [1] 38.5621  

(9)

  例題 1 データ・フレーム data2 を利用して,次の説明変数で就労時間hour.sfを説明す るモデルを推定しなさい。 1. 年齢,年齢の2乗,性別ダミー,就労形態ダミー(job.sf) 2. 10歳階級年齢ダミー(age.sf10),性別ダミー,就労形態ダミー(job.sf) ただし,factorを使用して名義尺度に対するダミー変数を定義すること。   例題1の 1の推定モデルとその推定は次のように書く。  

reg<- lm(hour.sf~age.sf + I(age.sf^2) + factor(sex) + factor(job.sf), data=data2)

summary(reg)

 

 

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 36.195526 5.322200 6.801 1.56e-11 *** age.sf 1.002056 0.208864 4.798 1.78e-06 *** I(age.sf^2) -0.012484 0.002071 -6.029 2.13e-09 *** factor(sex)2 -9.741887 0.859124 -11.339 < 2e-16 *** factor(job.sf)2 -5.885321 1.442993 -4.079 4.80e-05 *** factor(job.sf)3 -21.253979 1.584140 -13.417 < 2e-16 *** factor(job.sf)4 -12.676375 3.859618 -3.284 0.00105 ** factor(job.sf)5 -10.263944 1.667614 -6.155 9.88e-10 *** factor(job.sf)6 -5.737439 2.217513 -2.587 0.00978 **  

factor(job.sf)2からfactor(job.sf)6は定義より,経営者・役員 job.sf==1 を基準とする

ダミー変数である。つまり,「経営者・役員」と比べた就労時間は(他の条件が等しければ)「常時

雇用の一般従業者」では5.8時間だけ短い,「臨時雇用」では21.2時間だけ短い,などの結果が得

られる。

例題1の 2 の推定モデルとその推定は次のように書く。

 

reg<- lm(hour.sf~factor(age.sf10) + factor(sex) + factor(job.sf), data=data2)

summary(reg)

 

age.sf10は10歳階級ごとの年齢なので,factor(age.sf10)30からfactor(age.sf10)80ま

では,20歳代age.sf10==20を基準とするダミー変数である。20歳代と比較すると,30, 40, 50

歳代の就労時間は有意な差がないが,60歳代は20歳代に比べて8.0時間短く,70歳代は20歳代

(10)

3 一般化線形モデル 10

 

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 57.3306 2.2986 24.942 < 2e-16 *** factor(age.sf10)30 -2.0756 1.9640 -1.057 0.290799 factor(age.sf10)40 -1.6515 1.9093 -0.865 0.387197 factor(age.sf10)50 -3.1518 1.9292 -1.634 0.102551 factor(age.sf10)60 -8.0257 2.0113 -3.990 6.95e-05 *** factor(age.sf10)70 -14.6730 2.6060 -5.630 2.19e-08 *** factor(age.sf10)80 -23.0555 5.7030 -4.043 5.58e-05 *** factor(sex)2 -9.7094 0.8668 -11.202 < 2e-16 *** factor(job.sf)2 -5.9239 1.4554 -4.070 4.97e-05 *** factor(job.sf)3 -21.3883 1.5956 -13.405 < 2e-16 *** factor(job.sf)4 -13.3221 3.8877 -3.427 0.000629 *** factor(job.sf)5 -10.7186 1.6755 -6.397 2.18e-10 *** factor(job.sf)6 -5.8525 2.2318 -2.622 0.008831 **  

3

一般化線形モデル

3.1

ロジット・リンク

線形確率モデル work.sfは • work.sf==1「仕事をした」 • work.sf==2「休んだ(病気,休暇)」 • work.sf==3「仕事をしていない」 というデータであり,これをダミー変数 Yi= { 1 仕事をした 0 それ以外 i = 1, 2,· · · , n (4) のように2項変数に変換して, Yi= α + βXi+ ui (5) を推定する場合について考える(これを線形確率モデルとよぶ)。(5)式の推定は,左辺に2項変数 (ダミー変数)があるという点で特徴があるが,最小2乗法による推定は可能である(ただし,誤 差項は分散不均一である)。しかしながら,その理論値が必ずしも0と1の間に収まるとは限らず, 「就労確率」を0と1の間で定義するときには推定結果を理解するのは困難である。 Yiが0か1のデータということは,誤差項Yiはベルヌーイ分布(試行回数1の二項分布)にし たがっていると考えることができる。つまり,真の成功確率(仕事をする確率)pi が存在すると き,Yi = 0や Yi = 1は1回の試行による実現値である。また,E(Yi) = piV (Yi) = pi(1− pi) である。

(11)

潜在変数による回帰 仕事を「するか」,「しないか」という意思決定をモデル化する場合,(5)式に替えて,次の2項 選択モデルを定義する。 Yi∗= α + βXi+ ui (6) Yi∗< 0ならばYi= 0 Yi∗≧ 0ならばYi= 1 ここで,Yi∗ は観察できないが,何かを選択するときに間接的に観察される潜在変数 (latent variable)である。潜在変数Yi∗は2項変数 Yiとは異なって,−∞からの範囲で連続的な値 をとる量的変数である。α + βXi+ ui によってYi∗ < 0ならば,Yi = 0という意思決定をし, Yi≧ 0ならば,Yi= 1という意思決定をする(「仕事をする」を選択する)。つまり,潜在変数が 0という閾値を超えるか否かを計測することによって,選択確率を推定する。 リンク関数 通常,回帰モデルでは誤差項uiが正規分布であることを想定しているので,(5)を最小2乗法で 計算することはできるものの,最小2乗推定量の良い性質は失われてしまう。誤差項が正規分布以 外のケースに回帰モデルを適用できるように一般化したものを「一般化線形モデル」(generalized linear model; GLM)とよんでいる。一般化線形モデルは「系統的部分」,「リンク関数」,「誤差項」 の3点で記述することができる。「系統的部分」とは被説明変数Yiを説明変数Xiでどのように説 明するのかを示すα + βXiのことである。「リンク関数」とは,Yiの期待値(E(Yi) = pi)と系統 的部分との関係を繋ぐ関数のことであり,関数gを用いて g(E(Yi)) = α + βXi (7) と表現する。つまり,被説明変数の期待値をリンク関数gで変換したときの形が系統的部分に等し くなる。例えば,通常の線形回帰モデルでui∼ N(0, σ2)の場合E(Yi) = α + βXi となる。これ を恒等リンクという。 ロジット変換 (4)式のような0 か1 の2 項変数の場合は,リンク関数をロジット変換で定義できる。確率 Pr(Yi= 1) = pi をロジスティック分布の累積確率で表現するとき,pi= e α+βXi 1+eα+βXi であり,これ をロジット変換するとln pi 1−pi になる *1。ベルヌーイ分布より,E(Y i) = piなので,リンク関数 g(pi)を ln pi 1− pi = α + βXi (8) と書ける。これをロジット・リンクとよんでいる。 *1第10回講義ノートを参照。

(12)

3 一般化線形モデル 12

3.2

最尤法

ロジット・モデルを推定するために最尤法 (maximum likelihood method; ML)とよばれる推

定手法を用いる。最尤法とはモデルの尤度(likelihood)を最大化するパラメータを求める手法であ

り,尤度とは与えられた標本のもとでの確率の積で表現できる。*2

L = Pr(Y1= y1)× Pr(Y2= y2)× · · · × Pr(Yn = yn) (9)

ここで,yiは0または1の実現値である。 ロジット・モデルの尤度は次のように考える。まず,試行回数n,成功確率pの2項分布B(n, p) にしたがう確率変数をY とおくとき,y回成功する確率はPr(Y = y) =nCypy(1− p)n−y である ので,試行回数n = 1のベルヌーイ分布の場合は,成功回数Y = 0, 1について Pr(Y = 0) = p0(1− p)1−0= 1− p, Pr(Y = 1) = p1(1− p)1−1= p となる。つまり,i番目のデータが成功するか(Yi= 1)しないか(Yi= 0)を表現する確率は Pr(Yi= yi) = pyii(1− pi)1−yi (10) と書ける。ここで,(7)より,与えられたXiのもとでYi= 1となる確率は Pr(Yi= 1) = Pr(Yi∗≧ 0) (11) = Pr(α + βXi+ ui≧ 0) = Pr(ui≧ −α − βXi) = pi であり,Yi= 0となる確率は Pr(Yi= 0) = Pr(Yi∗< 0) (12) = Pr(α + βXi+ ui< 0) = Pr(ui<−α − βXi) = 1− pi である。したがって,(12),(13)を(10)に代入して尤度は次のように書ける。

L = Pr(Y1= y1)× Pr(Y2= y2)× · · · × Pr(Yn= yn)

= Pr(u1≧ −α − βX1)y1Pr(u1<−α − βX1)1−y1

× Pr(u2≧ −α − βX2)y2Pr(u2<−α − βX2)1−y2

× · · ·

× Pr(un≧ −α − βXn)ynPr(un<−α − βXn)1−yn

X1,· · · , Xn を与えると,尤度L はパラメータα, β の値によって決まる関数であることがわか

る。最尤法は尤度Lを最大化するα, βを見つける方法である。このとき,尤度の対数(対数尤度)

LL = ln Lα, βについて最大化すると最尤推定量α, ˆˆ βが得られる。

(13)

3.3

一般化線形モデルの評価

最尤推定量の仮説検定β = 0 は標準正規分布で行う。すべての説明変数の係数が0であるとい

う帰無仮説を検定するには尤度比検定(likelihood ratio test; LR)を利用する。実際には尤度の比

率ではなく,尤度の対数差分で検定を行う。

LR =−2(LL0 − LL) ∼ χ2(k) (13)

ここで,LL0はすべての説明変数の係数を0と仮定した場合の最大対数尤度,kは説明変数の数で

ある。自由度kのカイ2乗分布で有意性を検定する。これは,線形回帰モデルでのF 検定に該当

する。

一般化線形モデルを評価する場合,赤池情報基準量 (Akaike’s Information Criterion; AIC)が 用いられる。 AIC =−2LL + 2(k + 1) (14) いくつかのモデルを比較する場合,AICが最小のモデルを選択することが望ましい。 一般化線形モデルでは,線形回帰モデルで用いられる残差2乗和に替えて,残差逸脱度(residual deviance) Devという指標を用いる。 Dev =−2LL (15) 逸脱度は最大対数尤度に−2を乗じた値で定義する。LRAIC などの統計量は,この表現の方

がすっきりする(LR = Dev0− Dev, AIC = Dev + 2k)。

対数尤度や逸脱度は他のモデルとの比較では重要であるが,その値自体でモデルの適合度 を判断できない。一般化線形モデルではマクファーデンの擬似決定係数(McFadden’s pseudo R-squared)という指標が利用される。 pR2= 1 LL LL0 (16) 計算は単純でありLR検定を行うときに利用したものと同じ2つの対数尤度を使う。

3.4

ロジット・モデルの推定

glmコマンド Rでロジット・モデル(ロジスティック回帰モデル)を推定するにはglmコマンドを利用する。   y<- ifelse(work.sf==1 1,0)

y.glm<- glm(y~ age.sf + factor(sex), family=binomial(link="logit")) summary(y.glm)   2 項 変 数 (4) を work.sf デ ー タ か ら 作 成 し て y と す る 。y を 年 齢 age.sf と 女 性 ダ ミ ー factor(sex)で説明するモデルをglmコマンド内で記述して,さらにfamilyで目的変数が2項 変数binomialであること,利用するリンク関数がlink="logit"であることを指定する。推定 結果は次のようになる。

(14)

3 一般化線形モデル 14

 

Coefficients:

Estimate Std. Error z value Pr(>|z|) (Intercept) 5.527970 0.273707 20.20 <2e-16 *** age.sf -0.077703 0.004228 -18.38 <2e-16 *** factor(sex)2 -1.381965 0.109007 -12.68 <2e-16 ***

---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2877.6 on 2161 degrees of freedom Residual deviance: 2360.4 on 2159 degrees of freedom AIC: 2366.4

Number of Fisher Scoring iterations: 4

 

ここで,ロジット・モデルのパラメータは標準正規分布で仮説検定を行うので,z valueを表示し

ている。推定モデルは2つの説明変数を利用しているのでk = 2である。Null devianceは2つ

の説明変数がどちらも0であるという帰無仮説のもとでの逸脱度 (Dev0 =−2LL0) を示してい

る。サンプル・サイズn = 2162に対して,Null devianceの自由度はn− 1である。Residual devianceは最大対数尤度LLに対してDev =−2LLであり,自由度はn− (k + 1)である。赤 池情報基準量はAIC = Dev + 2(k + 1) = 2360.4 + 2× 3となっている。 モデルの評価 sumamryコマンドでは,対数尤度,尤度比検定,擬似決定係数などの情報は表示されないが, devianceなどの情報を利用して簡単に計算することができる。   LL<- y.glm$dev/(-2) LL0<- y.glm$null.dev/(-2) LR<- -2*(LL0-LL) bhat<- y.glm$coef k<- length(bhat)-1 LRp<- 1-pchisq(LR,df=k) AIC<- y.glm$aic pR2<- 1 - LL/LL0   逸脱度は$devで表示できるので,(15)よりy.glm$dev/(-2)が最大対数尤度LLになる。推定 値だけを取り出すには,$coef を利用する。lengthで推定パラメータの数を数えて,1を引くと k = 2が得られる。尤度比は(13)式より計算し,そのp値は自由度kのカイ2乗分布を利用して 計算する。赤池情報基準量は$aicで保存されている。擬似決定係数は(16)式より,対数尤度を利

(15)

表2 ロジット・モデルの推定結果(観測値の数2162) 変数 推定値 標準誤差 限界効果M E 定数項 5.528 0.274 *** -年齢 -0.078 0.004 *** -0.014 女性ダミー -1.382 0.109 *** -0.254 最大対数尤度LL -1180.209 尤度比 LR 517.220 尤度比LRp値 0.000 赤池情報基準量AIC 2366.417 擬似決定係数pR2 0.180 用して計算する。 限界効果 ロジット・モデルでは,説明変数と確率の関係は非線形になっているので,係数推定値だけでは 具体的な効果を知ることができない。そこで,説明変数Xが1単位変化したときのPr(Y = 1)と なる確率の変化である限界効果(marginal effects)を求める(第10回講義ノート)。ロジット・モ デルの限界効果は M E = ∂ Pr(Yi= 1) ∂Xi = ˆβp∗(1− p∗) (17) より計算する。ここで,p∗(1− p∗) はロジスティック分布の確率密度を示している。確率密度は Xiの値によって変動するので,通常は代表的な値を利用して評価する。ここでは,確率の理論値 ˆ pi= Pr(ui≧ −ˆα − ˆβXi)より p∗(1− p∗) = ∑ ipˆi(1− ˆpi) n を利用して平均的な限界効果を求める。 確率の理論値はy.glm$fit に格納されているので,次のように記述する。   phat<-y.glm$fit ME<- bhat*mean(phat*(1-phat)) ME (Intercept) age.sf D 1.01781992 -0.01430682 -0.25444983  

ここで,bhat は係数推定値,phat は確率の理論値,mean(phat*(1-phat))は確率密度の平均

値を示している。つまり,(平均的には)就労確率は,年齢age.sfが1歳上昇すると,0.0143 (1.43%)だけ低下し,女性Dは男性に比べて 0.2544 (25.44%)だけ低いことが示される。 以上より,ロジット・モデル ln pi 1− pi = β1+ β2age.sfi+ β3Di+ ui の推定結果一覧表を表2のようにまとめることができる。

(16)

3 一般化線形モデル 16

 

例題 2 ownhouse==1 ならば y = 1,それ以外はy = 0 となる 2 項変数を作り,年齢

age.sf,世帯年収factor(income.hh),都道府県 factor(pref)を説明変数とする持ち家

率のロジット・モデルを推定しなさい。また,年齢や年収をコントロールした場合,どの都道 府県の持ち家率が高いのかを限界効果で考えなさい。  

3.5

プロビット・モデル

プロビット・リンク ロジット・モデルは,(12)式のpiをロジスティック分布における確率と考えてリンク関数を定 義している。これとは別に,(12)式のpiを正規分布N (0, σ2)における確率と考えてリンク関数を 定義するのがプロビット・モデル(probit model)である。 pi= Pr(ui≦ α + βXi) = ∫ α+βXi −∞ f (ui)dui (18) ここで,f (ui)は正規分布の確率密度関数である。ロジット・モデルと同様に,プロビット・モデ ルも最尤法でパラメータ推定を行う。(18)におけるf (ui) の累積確率をF (α + βXi)と書くと, α + βXiについて解いた逆関数はF−1(pi)なので,プロビット・リンクは F−1(pi) = α + βXi (19) と書ける。 プロビット・モデルの推定 プロビット・モデルの推定もglmコマンドを利用する。yを仕事をするかしないかの2項変数と して,これを年齢と女性ダミーで説明するモデルを次のように記述する。  

y.glm<- glm(y~ age.sf + factor(sex), family=binomial(link="probit")) summary(y.glm)

 

推定結果を抜粋するとは次のようになる

 

Coefficients:

Estimate Std. Error z value Pr(>|z|) (Intercept) 3.262107 0.150845 21.63 <2e-16 *** age.sf -0.045510 0.002371 -19.19 <2e-16 *** factor(sex)2 -0.854726 0.063304 -13.50 <2e-16 ***   対数尤度やAIC はロジット・モデルのときと同じように計算できる。ただし,限界効果を計算す るには,標準正規分布の確率密度関数を利用する。 M E = ∂ Pr(Yi= 1) ∂Xi = ˆβif (α + βXi) n (20) Rでは次のように計算する。

(17)

表3 プロビット・モデルの推定結果(観測値の数2162) 変数 推定値 標準誤差 限界効果 定数項 3.262 0.151 *** -年齢 -0.046 0.002 *** -0.014 女性ダミー -0.855 0.063 *** -0.265 最大対数尤度LL -1180.950 尤度比 LR 517.220 尤度比LRp値 0.000 赤池情報基準量 AIC 2367.900 擬似決定係数pR2 0.180   phat<-y.glm$fit ME<- bhat*mean(dnorm(qnorm(phat))) ME

(Intercept) age.sf factor(sex)2 1.0103654 -0.0140957 -0.2647326   以上のプロビット推定の結果は表3にまとめられる。   例題 3 例題2をプロビット・モデルで推定しなさい。  

3.6

順序選択モデル

区切り点選択 順序選択モデル(順序ロジットあるいはプロビット・モデル)は,目的変数が順序尺度のデータ の場合に用いられる分析手法である。例えば,op56: 「子供の教育」は個人や家族の責任でしょうか,国や地方自治体の責任でしょうか? という質問について,回答の選択肢 Yi ={1, 2, 3, 4, 5} は値が大きくなるにつれて,国や自治体 (公)の責任をより重視し,値が小さくなるにつれて,個人や家族(個)の責任を重視するという順 序になっている。回答番号Yiを被説明変数にして最小2乗法で推定することもできるが,ロジッ ト・モデルで説明したように結果の解釈は難しくなる。 目的変数が順序尺度のデータの場合も(7)と同様に潜在変数を利用した回帰モデルを考える。 Yi∗= βXi+ ui (21) ζj−1< Yi∗< ζj ならばYi= j (j = 1, 2, 3, 4, 5) このモデルは,区切り点 (breakpoints) ζ0, ζ1, ζ2, ζ3, ζ4, ζ5 で区切られた区間に潜在変数が入ると き,選択肢番号j を選択することを示している。切片は一つではなく区切り点に応じて設定さ

(18)

3 一般化線形モデル 18 れる。 Yi= 1 ↔ ζ0< Yi∗< ζ1 ↔ ζ0− βXi< ui< ζ1− βXi Yi= 2 ↔ ζ1< Yi∗< ζ2 ↔ ζ1− βXi< ui< ζ2− βXi Yi= 3 ↔ ζ2< Yi∗< ζ3 ↔ ζ2− βXi< ui< ζ3− βXi Yi= 4 ↔ ζ3< Yi∗< ζ4 ↔ ζ3− βXi< ui< ζ4− βXi Yi= 5 ↔ ζ4< Yi∗< ζ5 ↔ ζ4− βXi< ui< ζ5− βXi ここで,ζ0=−∞, ζ5=であるので,推定モデルで決定されるのは ζ1, ζ2, ζ3, ζ4 となる。 Λをロジスティック分布(累積確率)とすると,順序ロジット・モデルのYi= jである(与えら れたXiのもとでの)確率は, pij = Pr(Yi= j) = Λ(ζj− βXi)− Λ(ζj−1− βXi) (22) となる。また,F を正規分布(累積確率)とすると,順序プロビット・モデルの確率は pij = Pr(Yi= j) = F (ζj− βXi)− F (ζj−1− βXi) (23) 2項選択ではなく,多項選択であるので確率分布も多項分布(pyi1 i1 )· (p yi2 i2 )· (p yi3 i3 )· (p yi4 i4 )· (p yi5 i5 )に なる。 モデルの評価 順序ロジット・モデルや順序プロビット・モデルの推定は,(22)や(23)を利用して尤度関数を 定義し,最尤法でパラメータ推定を行う。赤池情報基準量は AIC =−2LL + 2(k + kζ) (24) により計算する。ただし,kは説明変数の数, は推定すべき区切り点の数(選択肢が5段階の場 合ζ1, ζ2, ζ3, ζ4)である。 限界効果は密度関数の差分を利用して計算する。 順序ロジット・モデル: ∂ Pr(Yi= j) ∂Xi = [λ( ˆζj−1− ˆβXi)− λ(ˆζj− ˆβXi)] ˆβ (25) 順序プロビット・モデル: ∂ Pr(Yi= j) ∂Xi = [f ( ˆζj−1− ˆβXi)− f(ˆζj − ˆβXi)] ˆβ (26) ここで,λはロジスティック分布の密度関数,f は正規分布の密度関数である。 順序選択モデルの推定 目的変数となる順序尺度データを「子供の教育」op56として,説明変数を年齢age.sf,性別 ダミーfacto(sex),6歳未満の子供の数c06,6歳以上12歳未満の子供の数c0612,12歳以上 18歳未満の子供の数c1218,学歴ダミーfactor(edu.sf)とする。順序選択モデルを推定するに は,パッケージpkgを読み込む必要がある。パッケージには,あらかじめインストールされている ものと,外部からダウンロードして新しくインストールする必要があるものがある。順序選択モデ ルで利用するパッケージはMASSという名前で,Rには初期状態でインストールされている。MASS パッケージを読み込むにはlibraryコマンドを使用する。

(19)

  library(MASS)

 

順序選択モデルを推定するためにpolrコマンドを次のように使う。

 

ologit<- polr(as.factor(op56)~ age.sf + factor(sex)

+ c06 + c0612 + c1218 + factor(edu.sf), method = "logistic") summary(ologit)

 

ここで,op56 を順序尺度として認識させるために as.ordered というコマンドを使う。順序ロ

ジット・モデルのときはmethod = "logistic"とし,順序プロビット・モデルのときはmethod = "probit"とする。

推定モデルを評価するために,対数尤度やAIC などを以下のように計算する。

 

LL<- ologit$deviance/(-2)

ologit0<- polr(as.ordered(op56)~ 1, method = "logistic") LL0<- ologit0$deviance/(-2) LR<- -2*(LL0-LL); # 尤度比 LRp<- 1-pchisq(LR,df=k) # H0:「係数がすべて0」の尤度比検定のp値 bhat<- matrix(ologit$coef) # 推定値 k<- length(bhat) # 説明変数の数 kz<- length(ologit$zeta) # 順序尺度のカテゴリー数 - 1 AIC<- -2*LL + 2*(k+kz) # 赤池情報基準量 pR2<- 1-(LL/LL0) # 擬似決定係数   polrコマンドでは,説明変数の係数をすべて0と仮定した切片だけのモデルの対数尤度や逸脱度

(null.deviance)は出力されないので,自前で polr(as.ordered(op56)~ 1, method = "logistic")

を推定してLL0を計算する。 限界効果を求めるには誤差項uiの確率密度の差分を計算する必要がある。ologit$fitには (22)式で定義された確率の理論値が保存されているので,これを利用して確率密度を求める。   phat<- ologit$fit cphat<- array(0,c(ologit$n,4)) cphat[,1]<- phat[,1]

cphat[,2]<- phat[,1] + phat[,2]

cphat[,3]<- phat[,1] + phat[,2] + phat[,3]

cphat[,4]<- phat[,1] + phat[,2] + phat[,3] + phat[,4] lambda<- dlogis(qlogis(cphat))

 

ここで,lambdaは4列のデータで(25)式における確率密度λ(ζ1−βXi),λ(ζ2−βXi),λ(ζ3−βXi),

(20)

3 一般化線形モデル 20   dl1<- mean(-lambda[,1]) dl2<- mean(lambda[,1]-lambda[,2]) dl3<- mean(lambda[,2]-lambda[,3]) dl4<- mean(lambda[,3]-lambda[,4]) dl5<- mean(lambda[,4])   最後に(25)式にしたがって,各係数の限界効果を計算する。  

ME<- cbind(bhat*dl1, bhat*dl2, bhat*dl3, bhat*dl4, bhat*dl5) ME   限界効果M Eの出力結果は次のようになる。   [,1] [,2] [,3] [,4] [,5] [1,] -0.0013880099 -0.0010179285 0.0002881974 0.001078294 0.001039447 [2,] -0.0037981183 -0.0027854363 0.0007886166 0.002950620 0.002844318 [3,] -0.0203954382 -0.0149574575 0.0042347764 0.015844473 0.015273646 [4,] -0.0300370033 -0.0220283180 0.0062366885 0.023334654 0.022493979 [5,] -0.0336688021 -0.0246917801 0.0069907716 0.026156066 0.025213744 [6,] 0.0197321449 0.0144710162 -0.0040970545 -0.015329185 -0.014776922 [7,] 0.0009674851 0.0007095272 -0.0002008823 -0.000751604 -0.000724526 [8,] 0.0173066125 0.0126921970 -0.0035934327 -0.013444877 -0.012960500 [9,] -0.0097591267 -0.0071570770 0.0020263217 0.007581510 0.007308372   以上の推定結果は表4にまとめることができる。推定結果より,「子供の教育」に対する「公」と 「個」の選択には年齢,6歳以上12差未満の子供の数,12歳以上18差未満の子供の数が有意に正 の関連性を持っていることがわかる。その一方で,性別,6歳未満の子供の数,学歴の違いは選択 に影響しない。このことから,本人の学歴に関わらず,年齢が高いほど,そして就学している子供 の数が多いほど,「子供の教育」を国や自治体の責任と考える傾向が強いことがわかる。 限界効果で具体的な確率をみてみよう。年齢が1歳上昇すると, 「個人や家族の責任」と考える確率が0.14%だけ下落する。 「どちらかといえば個人や家族の責任」と考える確率が0.10%だけ下落する。 「どちらかといえば国や自治体の責任」と考える確率が0.11%だけ上昇する。 「国や自治体の責任」と考える確率が0.10%だけ上昇する。 また,6歳以上12差未満の子供の数が一人増えると, 「個人や家族の責任」と考える確率が3.00%だけ下落する。 「どちらかといえば個人や家族の責任」と考える確率が2.20%だけ下落する。 「どちらかといえば国や自治体の責任」と考える確率が2.33%上昇する。 「国や自治体の責任」と考える確率が2.25%だけ上昇する。

(21)

表4 順序ロジット・モデルの推定結果(観測値の数2162) 限界効果 推定値 標準誤差 Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 age.sf 0.010 0.004 ** -0.0014 -0.0010 0.0003 0.0011 0.0010 factor(sex)2 0.028 0.082 -0.0038 -0.0028 0.0008 0.0030 0.0028 c06 0.150 0.101 -0.0204 -0.0150 0.0042 0.0158 0.0153 c0612 0.220 0.080 *** -0.0300 -0.0220 0.0062 0.0233 0.0225 c1218 0.247 0.074 *** -0.0337 -0.0247 0.0070 0.0262 0.0252 factor(edu.sf)2 -0.145 0.282 0.0197 0.0145 -0.0041 -0.0153 -0.0148 factor(edu.sf)3 -0.007 0.275 0.0010 0.0007 -0.0002 -0.0008 -0.0007 factor(edu.sf)4 -0.127 0.291 0.0173 0.0127 -0.0036 -0.0134 -0.0130 factor(edu.sf)5 0.072 0.280 -0.0098 -0.0072 0.0020 0.0076 0.0073 区切り点: 1|2 ζ1 -0.972 0.409 2|3 ζ2 0.222 0.408 3|4 ζ3 1.529 0.409 4|5 ζ4 2.714 0.413 LL -3353.736 LR 22.522 LR p値 0.007 AIC 6733.473 pR2 0.003 12歳以上18差未満の子供の数が一人増えると, 「個人や家族の責任」と考える確率が3.37%だけ下落する。 「どちらかといえば個人や家族の責任」と考える確率が2.47%だけ下落する。 「どちらかといえば国や自治体の責任」と考える確率が2.62%だけ上昇する。 「国や自治体の責任」と考える確率が2.52%だけ上昇する。 なお,女性ダミーは有意でないが,男性に比べて女性は 「個人や家族の責任」と考える確率が0.38%だけ小さい。 「どちらかといえば個人や家族の責任」と考える確率が0.28%だけ小さい。 「どちらかといえば国や自治体の責任」と考える確率が0.30%だけ大きい。 「国や自治体の責任」と考える確率が0.28%だけ大きい。   例題 4 現在の仕事の満足度job.satiを順序ロジット・モデルで推定する。ただし,仕事 をしている人を対象にするために,次の新しいデータ・フレームを用意しておく。

data3<- data.frame(subset(data, work.sf == 1))

また,満足度の順序を「満足している」を5,「不満である」を1に変更するために,目的変数

をas.ordered(6-job.sati)とし,説明変数を年齢age.sf, 性別factor(sex),本人年

factor(income.sf),就労時間hour.sfとする。(library(MASS)でパッケージを読み

込んでおくこと。)

(22)

3 一般化線形モデル 22

推定結果は次のようになる。

 

> data3<- data.frame(subset(data, work.sf == 1))

> ologit<- polr(as.ordered(6-job.sati)~ age.sf + factor(sex) + + factor(income.sf) + hour.sf

+ , method = "logistic", data=data3) > summary(ologit)

Re-fitting to get Hessian

Call:

polr(formula = as.ordered(6 - job.sati) ~ age.sf + factor(sex) + factor(income.sf) + hour.sf, data = data3, method = "logistic")

Coefficients:

Value Std. Error t value age.sf 0.006856 0.004474 1.532 factor(sex)2 0.436100 0.131630 3.313 factor(income.sf)2 0.533906 0.144522 3.694 factor(income.sf)3 0.647301 0.158829 4.075 factor(income.sf)4 1.536530 0.236342 6.501 hour.sf -0.010312 0.003584 -2.877 Intercepts:

Value Std. Error t value 1|2 -3.2136 0.3643 -8.8205 2|3 -1.8616 0.3355 -5.5490 3|4 -0.3207 0.3277 -0.9786 4|5 1.4821 0.3301 4.4903 Residual Deviance: 3500.94 AIC: 3520.94  

表 2 ロジット・モデルの推定結果(観測値の数 2162 ) 変数 推定値 標準誤差 限界効果 M E 定数項 5.528 0.274 ***  -年齢 -0.078 0.004 *** -0.014 女性ダミー -1.382 0.109 *** -0.254 最大対数尤度 LL -1180.209 尤度比 LR 517.220 尤度比 LR の p 値 0.000 赤池情報基準量 AIC 2366.417 擬似決定係数 pR 2 0.180 用して計算する。 限界効果 ロジット・モデルでは,説明変数と確率
表 3 プロビット・モデルの推定結果(観測値の数 2162 ) 変数 推定値 標準誤差 限界効果 定数項 3.262 0.151 ***  -年齢 -0.046 0.002 *** -0.014 女性ダミー -0.855 0.063 *** -0.265 最大対数尤度 LL -1180.950 尤度比 LR 517.220 尤度比 LR の p 値 0.000 赤池情報基準量 AIC 2367.900 擬似決定係数 pR 2 0.180   phat&lt;-y.glm$fit ME&lt;- bhat
表 4 順序ロジット・モデルの推定結果(観測値の数 2162 ) 限界効果 推定値 標準誤差 Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 age.sf 0.010 0.004 ** -0.0014 -0.0010 0.0003 0.0011 0.0010 factor(sex)2 0.028 0.082 -0.0038 -0.0028 0.0008 0.0030 0.0028 c06 0.150 0.101 -0.0204 -0.0150 0.0042 0.0158 0.0153 c0612

参照

関連したドキュメント

Actually a similar property shall be first de- rived for a general class of first order systems including the transport equation and Schr¨odinger equations.. Then we shall consider

(The Elliott-Halberstam conjecture does allow one to take B = 2 in (1.39), and therefore leads to small improve- ments in Huxley’s results, which for r ≥ 2 are weaker than the result

[r]

“Breuil-M´ezard conjecture and modularity lifting for potentially semistable deformations after

lines. Notice that Theorem 4 can be reformulated so as to give the mean harmonic stability of the configuration rather than that of the separate foliations. To this end it is

S., Oxford Advanced Learner's Dictionary of Current English, Oxford University Press, Oxford

At the end of the section, we will be in the position to present the main result of this work: a representation of the inverse of T under certain conditions on the H¨older

支払方法 支払日 ※② 緊急時連絡先等 ※③.