SASでロジスティック回帰いろいろ
はじめに
SASでロジスティック回帰を行う場合、logisticプロシジャを使うのが一般的ですが、他のプロシジャを使ってもパラメータ推定ができるのでやってみました。
データ作成
* 真のパラメータを設定(これらを推定する) ; %let beta0 = 0.3; %let beta1 = 0.6; %let beta2 = 0.5; data temp; call streaminit(1000); do x1 = -5 to 5 by 0.2; do x2 = -5 to 5 by 0.2; p = 1 / (1 + exp(-(&beta0. + &beta1. * x1 + &beta2. * x2))); if p <= rand("uniform") then t = 1; else t = 0; output; end; end; run; * 作成したデータを描画 ; proc sgplot data = temp; scatter x = x1 y = x2 / group = t markerattrs=(symbol=circlefilled); run;
実行結果
logisticプロシジャ
logisticプロシジャを使った基本的なやり方です。
proc logistic data = temp; model t = x1 x2; run;
実行結果(抜粋)
最尤推定値の分析 | |||||
---|---|---|---|---|---|
パラメータ | 自由度 | 推定値 | 標準誤差 | Wald カイ 2 乗 |
Pr > ChiSq |
Intercept | 1 | 0.2904 | 0.0550 | 27.9265 | <.0001 |
x1 | 1 | 0.6278 | 0.0256 | 600.3999 | <.0001 |
x2 | 1 | 0.5015 | 0.0237 | 447.6218 | <.0001 |
catmodプロシジャ
カテゴリ変数がある場合、logisticプロシジャを使うにはダミー変数化する必要がありますが、catmodプロシジャを使えばカテゴリ変数をそのままmodelステートメントに記述することができます。*1
逆に、連続値の変数はdirectステートメントに記述しておく必要があります。
proc catmod data = temp; direct x1 x2; model t = x1 x2; run; quit;
実行結果(抜粋)
Analysis of Maximum Likelihood Estimates | ||||
---|---|---|---|---|
Parameter | Estimate | Standard Error |
Chi- Square |
Pr > ChiSq |
Intercept | 0.2904 | 0.0550 | 27.93 | <.0001 |
x1 | 0.6278 | 0.0256 | 600.40 | <.0001 |
x2 | 0.5015 | 0.0237 | 447.62 | <.0001 |
genmodプロシジャ
一般化線形モデル用のgenmodプロシジャを使えば、確率分布に二項分布、リンク関数にロジット関数を指定することでロジスティック回帰を行うことができます。
proc genmod data = temp; model t = x1 x2 / dist = binomial link = logit; run;
実行結果(抜粋)
最大尤度パラメータ推定値の分析 | |||||||
---|---|---|---|---|---|---|---|
パラメータ | 自由度 | 推定値 | 標準誤差 | Wald 95% 信頼限界 | Wald カイ 2 乗 |
Pr > ChiSq | |
Intercept | 1 | 0.2904 | 0.0550 | 0.1827 | 0.3981 | 27.93 | <.0001 |
x1 | 1 | 0.6278 | 0.0256 | 0.5776 | 0.6780 | 600.40 | <.0001 |
x2 | 1 | 0.5015 | 0.0237 | 0.4550 | 0.5480 | 447.62 | <.0001 |
Scale | 0 | 1.0000 | 0.0000 | 1.0000 | 1.0000 |
genmodプロシジャ(ベイズ推定)
genmodプロシジャでは、bayesステートメントを追記することによってベイズ推定を行うことができます。
proc genmod data = temp; model t = x1 x2 / dist = binomial link = logit; bayes; run;
この例ではデフォルト設定で使っていますが、事前分布のパラメータなどを指定することもできます。
実行結果(抜粋)
事後要約 | ||||||
---|---|---|---|---|---|---|
パラメータ | N | 平均 | 標準偏差 | パーセント点 | ||
25% | 50% | 75% | ||||
Intercept | 10000 | 0.2904 | 0.0552 | 0.2526 | 0.2905 | 0.3273 |
x1 | 10000 | 0.6292 | 0.0259 | 0.6113 | 0.6286 | 0.6465 |
x2 | 10000 | 0.5026 | 0.0235 | 0.4866 | 0.5024 | 0.5182 |
nlmixedプロシジャ
非線形混合モデル用のnlmixedプロシジャを使えば、二項分布とロジスティック関数を組み合わせることでロジスティック回帰を行うことができます。
proc nlmixed data = temp; parms beta0 = 0 beta1 = 0 beta2 = 0; theta = 1 - logistic(beta0 + beta1 * x1 + beta2 * x2); model t ~ binomial(1, theta); run;
実行結果(抜粋)
Parameter Estimates | |||||||||
---|---|---|---|---|---|---|---|---|---|
Parameter | Estimate | Standard Error | DF | t Value | Pr > |t| | Alpha | Lower | Upper | Gradient |
beta0 | 0.2904 | 0.05497 | 2601 | 5.28 | <.0001 | 0.05 | 0.1826 | 0.3982 | -0.00013 |
beta1 | 0.6278 | 0.02563 | 2601 | 24.50 | <.0001 | 0.05 | 0.5776 | 0.6781 | -0.00137 |
beta2 | 0.5015 | 0.02371 | 2601 | 21.16 | <.0001 | 0.05 | 0.4550 | 0.5480 | -0.0003 |
modelステートメントを
model t ~ binary(theta);
と書いたり、対数尤度関数を記述して
proc nlmixed data = temp; parms beta0 = 0 beta1 = 0 beta2 = 0; theta = 1 - logistic(beta0 + beta1 * x1 + beta2 * x2); llh = t * log(theta) + (1 - t) * log(1 - theta); model t ~ general(llh); run;
とすることもできます。
mcmcプロシジャ
最後にベイズ推定用のmcmcプロシジャを使ってパラメータのベイズ推定を行います。
mcmcプロシジャの書き方はnlmixedプロシジャに似ています。
proc mcmc data = temp nbi = 2000 nmc = 10000; parms (beta0 beta1 beta2) = 0; prior beta0 beta1 beta2 ~ normal(mean=0, var=10000); theta = 1 - logistic(beta0 + beta1 * x1 + beta2 * x2); model t ~ binomial(1, theta); run;
genmodプロシジャのデフォルト設定に合わせてバーンインを2000、サンプリング数を10000に設定しています。
実行結果(抜粋)
Posterior Summaries and Intervals | |||||
---|---|---|---|---|---|
Parameter | N | Mean | Standard Deviation |
95% HPD Interval | |
beta0 | 10000 | 0.2901 | 0.0545 | 0.1820 | 0.3872 |
beta1 | 10000 | 0.6298 | 0.0253 | 0.5807 | 0.6768 |
beta2 | 10000 | 0.5028 | 0.0233 | 0.4571 | 0.5484 |
genmodプロシジャの実行結果と合わせたかったのですが、細かいオプションの使い方がまだ分かっていないせいでサンプリングがちょっと上手く行っていないようです。*2
nlmixedプロシジャでの書き換えと同様に、modelステートメントを
model t ~ binary(theta);
と書いたり、対数尤度関数を記述して
proc mcmc data = temp nbi = 2000 nmc = 10000; parms (beta0 beta1 beta2) = 0; prior beta0 beta1 beta2 ~ normal(mean=0, var=10000); theta = 1 - logistic(beta0 + beta1 * x1 + beta2 * x2); llh = t * log(theta) + (1 - t) * log(1 - theta); model t ~ general(llh); run;
とすることもできます。
まとめ
一般化線形モデルの確率分布とリンク関数の関係はみどりぼんを、ロジスティック回帰の対数尤度関数の書き方はぱじパタを読むことで理解できたので、読んでよかったなと思います。