担当 : 長倉大輔 ( ながくらだいすけ ) この資料は私の講義において使用するために作成した資料です WEB ページ上で公開しており 自由に参照して頂いて構いません ただし 内容について 一応検証してありますが もし間違いがあった場合でもそれによって生じるいかなる損害 不利益について責任を負いかねますのでご了承ください 間違いは発見次第 継続的に直していますが まだ存在する可能性があります 1
カウントデータモデル 労働経済学などで扱うデータには連続変数ではなく離散変数 特に非負の整数の値を取るデータが多い 例えば ある個人による通院の回数などは非負の整数のデータである このようなデータはカウントデータと呼ばれる そのような変数をその性質を考慮しないで たとえば線形回帰分析などを行うと いろいろと不都合なことが起きる 以下ではこのようなデータを扱うのに適したミクロ計量モデルを紹介する 2
ポアソン分布 カウントデータのような非負の整数の値を取る離散型の分布の代表的なものにポアソン分布がある 確率変数 Y がポアソン分布に従うとは Y の確率関数が以下のように与えられる時である Pr(Y = y) = exp( ) y! y, y = 0, 1, 2, ここで λ は正の実数である ( 整数ではないことに注意 ) この分布の期待値は λ で与えられ 分散も λ で与えられる 以下では確率変数 Y がパラメーター λ のポアソン分布に従うことを Y ~ Po(λ) と表すとする 3
ポアソン回帰モデル 前述のポアソン分布においては期待値は λ で一定である より柔軟性のあるモデルにしたい場合 特に説明変数による影響を考えたい場合 この仮定は制約的である 例えば線形回帰分析などは 被説明変数 Y i の期待値を説明変数の線形の関数でモデル化していると考えることができるが 同じような拡張をポアソン分布に対して行ったものがポアソン回帰モデルと呼ばれるモデルである 4
線形回帰分析の時と同様に ポアソン分布の期待値 λ を説明変数の関数として表すことを考える ただし λ は常に正の値をとらないといけないので ポアソン回帰モデルでは それを考慮して λ i = exp(β 0 + β 1 X 1i + + β K X Ki ), のように定式化する このように定式化すると λ i は各説明変数 X ki の値に依存して変わり また X ki の値によらず正の実数をとる λ i の定式化としては 分析の目的に応じて他の定式化を使用することも可能であるが 上記の定期化が非常によく用いられる 5
よって非負の整数の値をとる被説明変数 Y i に対するポアソン回帰モデルは以下で与えられる Pr(Y i = y) = λ i = exp(β 0 + β 1 X 1i + + β K X Ki ) λ i の部分は exp( ) i y! y i, y = 0, 1, 2,. log λ i = β 0 + β 1 X 1i + + β K X Ki と書き直すこともできる 6
最尤法によるポアソン回帰モデルの推定 ポアソン回帰モデルは最尤法で簡単に推定できる β = [β 0, β 1,.., β K ] T とすると 実現値 y 1, y 2,., y n が与えられた下での 対数尤度関数は log L(β) i= 1 λ i = exp(β 0 + β 1 X 1i + + β K X Ki ) によって与えられる n = log Pr( Y = y ) i n n n = + y log y! i i i i i= 1 i= 1 i= 1 i 7
最大化の 1 階の条件は log L( β) n = X + y X = 0, k = 0,..., K ki i i ki k i= 1 i= 1 n である ここで X 0i = 1 とする この連立方程式は β について明示的には解けないので ポアソン回帰モデルの最尤推定値は数値的に求めることになる 8
例題 1 ポアソン回帰モデル Y i ~ Po(λ i ), λ i = exp(β 0 + β 1 X i ) を考える ここで X i は i 番目の観測値に対する説明変数である Y i, i =1,,n は λ i が与えられたもとで互いに独立であるとする 今 n = 3 であり X i と Y i の観測値として X 1 = 1, X 2 = 1, X 3 = 1, y 1 = 5, y 2 = 4, y 3 = 3 が与えられたとする この時 β 0 と β 1 の最尤推定値を求めなさい 9
デュレーションモデル デュレーションモデルとは何らかのイベントが起こる ( または終わる ) までの時間を分析したい場合によく用いられる 例えば 病気が治癒するまでの時間 製品が故障するまでの時間 失業が続く期間 などの分析である デュレーション分析は生存時間分析などとも呼ばれる 10
生存関数 確率変数 T は非負の実数とする ここで T はあるイベントが続いた時間を表すと解釈しよう T の分布関数 F(t) = Pr(T t) はこのイベントが続く時間が t 以下になる確率 言い換えると ( 時点 0 から始まった ) このイベントが時点 t までに終わる確率を表している この時 S(t) = 1 F(t) と定義される関数 S(t) は生存関数と呼ばれる S(t) はこのイベントが続く時間が t より大きくなる確率 つまりこのイベントが時点 t になっても続いている確率 を表している 11
ハザード関数 Δt を任意の小さな値を表すとする このイベントが続く時間が t 以上 t + Δt 以下となる確率は Pr(t < T t + Δt ) = F(t + Δt) F(t) である これは無条件確率であることに注意しよう 今 興味のある確率として このイベントが時点 t まで続いたという条件付きでさらに非常に小さな時間 Δt 続く確率を考えたいとする これは Pr(t < T t + Δt T t ) = と表すことができる F( t + t) F( t) St () 12
ここでこの条件付き確率と Δt の比 Pr( t T t + t T t) F( t + t) F( t) 1 = t t S() t を考え この比が Δt 0 の時にどのように表せるか考える これは条件付き分布関数の点 t における ( 右から極限を取った ) 傾きと解釈できる T の確率密度関数を f(t) とすると この極限は F( t + t) F( t) 1 f ( t) ht ( ) = lim = t 0 t S( t) S( t) と表すことができる これをハザード関数と呼ぶ デュレーション分析ではしばしばこのハザード関数をモデル化する 13
ウェイブル分布 T の分布として ウェイブル分布を考える ウェイブル分布の分布関数は以下で与えられる F(t ; β, α) = 1 exp( (βt) α ) ここで β > 0 である よってその密度関数 生存関数 ハザード関数はそれぞれ f (t ; β, α) = αβ α t α 1 exp( (βt) α ) S(t) = exp( (βt) α ) h(t) = αβ α t α 1 で与えられる 14
ウェイブル分布のハザード関数の性質 t が増加した時 ハザード関数が減少していく場合 負のデュレーション依存があるといい 増加していく場合 正のデュレーション依存があるという またハザード関数が t に依存しないときはデュレーション独立と呼ぶ ウェイブル分布場合は α < 1 負のデュレーション依存 α > 1 正のデュレーション依存 α = 1 デュレーション独立となる 15
対数正規分布 確率変数 Y の分布が対数正規分布であるとは log Y ~ N(μ, σ 2 ) となる分布のことである ウェイブル分布のハザード関数は t についての単調関数であり やや制約的である 実際の分析では ハザード関数が t が大きくなるにつれて 当初は増加するが 次第に減少に転じるような場合がよく観測される そのようなハザード関数はウェイブル分布のハザード関数では表現できないため 他の分布を用いる必要がある そのような分布としてよく用いられるのが対数正規分布である 16
対数正規分布の分布関数 密度関数 生存関数 ハザード関数は 1 F( t) = (log t ) 1 1 f ( t) (log t ) t = 1 1 S( t) = 1 (log t ) = (log t ) 1 ((log t ) / ) ht () = t (log t ) / ( ) で与えられる ここでΦ(.) および (.) はそれぞれ標準正規分布の分布関数と密度関数である 17
デュレーションモデルの最尤推定 T i を i 番目のデュレーションで確率変数とする また t i を i 番目のデュレーションの観測値とする T i は X i という説明変数ベクトルで条件付けした場合 条件付きで独立であるとし その条件付き密度関数は f (t i ; X i, θ) で与えられるとする ここで θ は未知パラメーターベクトルである この時 一般的に尤度関数は n log L( θ) = log f ( t ; X, θ) i= 1 i と書け これを最大化する未知パラメーター θ の値を求めれば最尤推定値が得られる i 18
例えばウェイブル分布において i 番目の分布の β を β = β i = exp(β 0 + β 1 X 1i + + β K X Ki ) とすると 先ほどのウェイブル分布の密度関数より f (t i ; X i, θ) = f(t i ; α, β 0,, β K ) = αβ iα t i α 1 exp( (β i α t iα )) である ここで β iα = exp(αβ 0 + αβ 1 X 1i + + αβ K X Ki ) であるので γ k = αβ k とすれば f (t i ; X i, θ) = f(t i ; α, γ 0,, γ K ) = αβ i t i α 1 exp( (β i t iα )) β i = exp(γ 0 + γ 1 X 1i + + γ K X Ki ) という表現を得ることもできる 19
よってその対数尤度関数は L 0 K = n + i + ti t i i i= 1 i= 1 i= 1 β i = exp(γ 0 + γ 1 X 1i + + γ K X Ki ) log (,,..., ) log log ( 1) log で与えられる n n n 20
例題 2 対数正規分布の分布関数 密度関数 生存関数 ハザード関数が前述のように得られることを確認しなさい 例題 3 デュレーション T i が log T i ~ N(μ, σ 2 ) の対数正規分布に従っている時 T i = t i が観測された時の μ, σ 2 についての対数尤度関数を書きなさい 21