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;
実行結果

f:id:Prunus1350:20150315204132p:plain

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

f:id:Prunus1350:20150315214440p:plain f:id:Prunus1350:20150315214450p:plain f:id:Prunus1350:20150315214456p:plain

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

f:id:Prunus1350:20150315221910p:plain f:id:Prunus1350:20150315221916p:plain f:id:Prunus1350:20150315221923p:plain

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;

とすることもできます。

まとめ

一般化線形モデルの確率分布とリンク関数の関係はみどりぼんを、ロジスティック回帰の対数尤度関数の書き方はぱじパタを読むことで理解できたので、読んでよかったなと思います。

*1:この例では利点がないですね。

*2:分かり次第記事を更新するつもりです。