910 likes | 1.53k Views
お試しポアソン回帰. 鳥居稔 行動データ科学M 2. データ. A データ 大阪府下の種類別の A の件数 1 日ごとにカウント 月データ 月と地球の間の距離・位相 2001 年度、 365 日分のデータをマージ. 具体的な仮説. 仮説 「満月・新月になるほど(月が地球に近づくほど) A の数が多くなる」. 使用するデータ. 大阪府下の A の数 2001 年度大阪府下で発生した 1 日ごとの件数 月の満月・新月度( fullcres ) 満月・新月のときに値が最も大きくなり、半月のときにもっとも小さくなる連続量. データセット. A のカウント.
E N D
お試しポアソン回帰 鳥居稔 行動データ科学M2
データ • Aデータ • 大阪府下の種類別のAの件数 • 1日ごとにカウント • 月データ • 月と地球の間の距離・位相 • 2001年度、365日分のデータをマージ
具体的な仮説 • 仮説 • 「満月・新月になるほど(月が地球に近づくほど)Aの数が多くなる」
使用するデータ • 大阪府下のAの数 • 2001年度大阪府下で発生した1日ごとの件数 • 月の満月・新月度(fullcres) • 満月・新月のときに値が最も大きくなり、半月のときにもっとも小さくなる連続量
データセット Aのカウント ・・・
A Aのカウント Count of A
Count of A 何か関係がありそうに見えませんか? Count of A
Count of A Count of A
Count of A Count of A
Count of A Count of A
Count of A Count of A
Aの数はポアソン分布 • ポアソン分布 • 二項分布における極限 • n(試行回数)→すごく大(p(発生確率)→すごく少) • ある町の交差点における一日の交通事故数 • 一日におこる間違い電話の数 • 今回のデータ • Aのカウント数→数えれるほど(分母はン百万!) • ポアソン分布ですね!
ポアソン分布 • 確率関数 • パラメータはλで一つ • 確率関数から、λが与えられた下でX=xが起こる確率が理論的に求まる • E(X)=V(X)=λ
ポアソン回帰分析 • 従属変数がポアソン変量 • ポアソン回帰分析を用いるのが適切 • 通常、期待値にlogをとり線形予測するモデル • モデル式log(λ)=Xβ(=β0+β1x1+β2x2+・・・)or λ=exp(Xβ)
*統計手法アラカルト楠元氏発表時の狩野先生のコメントスライドに少し手を*統計手法アラカルト楠元氏発表時の狩野先生のコメントスライドに少し手を 加えたもの ポアソン回帰モデル • ポアソン変量である応答変数を規定する要因を調べるため回帰分析したい • ポアソン分布 • 例 • λに直接回帰すると,λが負の値になることがある Nがわかっているときはその情報も加える
扱う要因 • 従属変数 • 大阪府下のAのカウント数 • 独立変数 • 満月・新月度 • さらに季節も要因として入れておこう • 温かい季節と寒い季節は差がありそう • 4水準の質的変数
モデル式 • log(λ)=β0+βfullcress+ β(季節i)+βi fullcress*季節i • fullcress • 満月・新月度(連続) • 季節 • 4水準の質的変数。3つのダミー変数が作られる。(今回は(SASが自動的に)冬を基準にしている) • fullcress*季節 • 満月新月度と季節の交互作用。季節ごとに満月新月度の影響力が異なるかどうかを検定する。
推定 • 最尤法 • ポアソン確率関数を用いて尤度関数を導く • 尤度を用いてモデルの逸脱度をチェックできる • 基本的にSASのGenmodプロシジャ • 独立変数は連続・質的を問わない
分析 • だいたいの手順 • 飽和モデルからの逸脱度 • サンプルが多いときは微小な差異を検出するのであまり気にすることはない、といわれる • 交互作用・主効果の有意性 • 交互作用が有意でなければ誤差にプール • パラメータ推定値の有意性 • 実際的な興味にもとづくリスク比の算出 • 残差分析 • モデル式とその信頼区間
SASプログラム PROC GENMOD ; class season ; model A = season fullcres / covb dist=poisson link=log type3 obstats; make 'obstats' out=red ; estimate 'summer vs winter' season1 0 1 0 -1 / exp; estimate 'fullcres in Summer VS NON_fullcres in Winter‘ fullcres 7.387 season 1 0 1 0 -1 / exp ; RUN;
結果①飽和モデルからの逸脱 ・2つの指標ともに棄却 ・モデルはあまり良くないなと心に留めつつ、サンプル数が多さが影響していそうなので話を先へ進める
結果②交互作用・主効果 ・非有意なので誤差にプール ・すなわち季節ごとにfullcresの影響は変わらないと考える
結果②交互作用・主効果 ・どちらの要因も有意!有意ですよ奥さん!
結果③推定値の有意性 ・新月・満月度の係数が正で有意 ・春、夏、秋の係数は(冬を基準にしたとき)正で有意
モデル式とその信頼区間 • モデル式 • λ=exp(-0.69+0.1*fullcres+0.5*spring+0.66*summer+0.46*autumn)(spring,summer,autumnは1,0のダミー変数) ^
モデル式とその信頼区間 • モデル式の95%信頼区間を求めたい • exp(Xβ±1.96SE(Xβ)) • Q:1.96というのは正規分布の両側5%点のはず。今はポアソン変量を検討しているのにどうして正規分布が出てくるのか • A:βについての信頼区間を求めたい。βは最尤推定量であり、漸近正規。だから正規分布を仮定する ^ ^ ^ ^
モデル式とその信頼区間 • SE(Xβ)=SE(b0+b1*fullcres+b2*spring+ b3*summer+b4*autumn) =sqrt(var(b1)+var(b2)+var(b3)+var(b4)+var(b5)+2*fullcres*Cov(b0,b1)+2*spring*Cov(b0,b2)+2*summer*Cov(b0,b3)+2*autumn*Cov(b0,b4)+2*fullcres*spring*Cov(b1,b2)+2*fullcres*summer*Cov(b1,b3)+2*fullcres*autumun*Cov(b1,b4))) • Var(X1+X2+・・・+Xn)=ΣVar(Xi)+2ΣCov(Xi,Xj) (i>j)(from 情報処理論課題)を利用 (---最尤推定量の共分散行列(フィッシャー情報行列の逆行列)はcovbオプションで出力)
モデル式とその信頼区間 実はSE(Xβ)はほとんど同じなんです・・ • 春 • (λL,λU)=exp(-0.19+0.1full±(0.02-0.004full)) • 夏 • (λL,λU)=exp(-0.03+0.1full±(0.02-0.004full)) • 秋 • (λL,λU)=exp(-0.23+0.1full±(0.02-0.004full)) • 冬 • (λL,λU)=exp(-0.69+0.1full±(0.02-0.004full)) ^ ^ ^ ^ ^ ^ ^ ^
相対リスク(RelativeRisk)(リスクの比) • 水準間でリスクがどれほど(相対的に)異なるのか? • 例:夏にAがカウントされる確率は冬の何倍か?P夏/P冬=(E(Y夏)/N夏)/ (E(Y冬)/N冬)=λ夏/λ冬(夏の分母と冬の分母を同じと考えます) ^ ^ ^ ^
相対リスク ^ ^ • λ夏/λ冬=exp(b0+b1*fullcres+b3*summer)/exp(b0+b1*fullcres)=exp(b3)同様に、例えば夏と秋の相対リスクはλ夏/λ秋=exp(b0+b1*fullcres+b3)/exp(b0+b1*fullcres+b4) =exp(b3-b4) ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^
SASでestimateステートメント estimate 'summer vs winter' season0 1 0 -1 / exp; ・要因を指定 ・水準を指定し、指定された値を推定せよと命令 --ここでは「夏の推定値 -冬の推定値」を推定せよという命令 ・その推定値にexpをとれという命令 同様に、b3-b4(夏の推定値-秋の推定値)を推定したいときは、 estimate 'summer vs winter' season0 1 –1 0 / exp;
SASでestimateステートメント α=0.05 ・夏と冬の差は有意、夏と秋の差は非有意 ・夏と冬の相対リスクの点推定値は1.93 --夏にAがカウントされる確率は冬の約1.93倍!
夏の満月・新月VS冬の半月 • 最も相対リスクが高いと思われる比較 estimate 'fullcres in Summer VS NON_fullcres in Winter' fullcres 7.387 season1 0 1 0 -1 / exp ; ・冬の半月度最高の日に比べて、夏の満月・新月度最高の日に Aがカウントされる確率は約4倍
結論 • 仮説が実証された! • 満月・新月度が有意 • 月が地球に近くなるほどAがカウントされやすい • 季節の要因は統制 • 冬における観測率が低い • 月はAの観測率に影響する!
え、いいんですか? これで終わっていいんですか?
やることはまだありますよー • 重要な視点が2つ抜けています • 残差分析 • 誤差のふるまいに異常はないですか? • 決定係数 • どれくらいのパーセンテージでモデルがデータを説明できていますか?
残差分析 • 標準化残差 • (実測値-予測値)/標準偏差 • ポアソン回帰の場合 • (y-λi)/sqrt(λi) • ピアソン残差ともいう • x軸に予測値、y軸に標準化残差をとりプロットしてみる ^ ^
残差分析 • 問題 • 一般に、ポアソン回帰において誤差は非正規 • だから標準化残差の扱い方がよくわかんない • 対処法 • anscombe残差を用いる3(y^(2/3)-λi^(2/3))/(2*λi^(1/6)) • 残差が正規分布するような変換がなされている ^ ^
過大評価されている サンプルが多いとい うこと??
決定係数 • Genmodプロシジャにおいては • 予測力を示す指標がない • でも、オーストリア人のMittlboeckという人がポアソン回帰用に決定係数を算出するマクロを作られていました • 平方和に基づく算出法、そしてDevianceに基づく算出法の2通り • Debiance、つまり尤度に基づいたやりかたの方が良いということが違う人の文献で紹介されている • CFIの平均構造版と考えてよいだろう
Count of A Count of A
Count of A Count of A
Count of A Count of A
Count of A Count of A
Count of A & Fullcres Count of A 夏 春 秋 冬