経済情報処理 講義ノート
第
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 順序選択モデル . . . 171 データの読み込み 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であることが示されている。
2
回帰分析
2.1
散布図と回帰直線
直線: curve, lines, ablinex-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.sp は formula とよばれるオブジェクトで,age.sf = α + β age.sp
という線形回帰直線を表現するための記法である。lmコマンドで計算した切片や傾きの情報を
regという名前で保存しておいて,これをablineに渡すことで散布図に回帰直線を描画してい
る。abline(lm(age.sf age.sp))と記述しても同じ結果が得られる。lmコマンドで計算した内
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
き(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 推定値の分散・共分散
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)
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時
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
例題 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歳代
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) = pi,V (Yi) = pi(1− pi) である。潜在変数による回帰 仕事を「するか」,「しないか」という意思決定をモデル化する場合,(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回講義ノートを参照。
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をα, βについて最大化すると最尤推定量α, ˆˆ βが得られる。
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を乗じた値で定義する。LRやAIC などの統計量は,この表現の方
がすっきりする(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"であることを指定する。推定 結果は次のようになる。
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)式より,対数尤度を利
表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 擬似決定係数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のようにまとめることができる。
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では次のように計算する。
表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 擬似決定係数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 を選択することを示している。切片は一つではなく区切り点に応じて設定さ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は説明変数の数,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コマンドを使用する。
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),
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%だけ上昇する。
表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)でパッケージを読み
込んでおくこと。)
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