【SAS】treeboostプロシジャ_Gradient Boosting Tree(勾配ブースティング木)

今回は Gradient Boosting Tree(勾配ブースティング木)を試してみる。

【目次】


決定木を学習器としたアンサンブル学習であり、ブースティングと言われる手法。
※XGBoost ではないのであしからず。

Pythonでの実装コード


Pythonでは scikit learn で実装できる。
sklearn.ensemble.GradientBoostingClassifier — scikit-learn 0.24.2 documentation

from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import classification_report, confusion_matrix
 
mod1 = GradientBoostingClassifier(n_estimators=100, learning_rate=1.0, max_depth=6, random_state=19)
mod1.fit(train_X, train_y)
pred_y = mod1.predict(valid_X)
 
print(confusion_matrix(valid_y, pred_y))
print(classification_report(valid_y, pred_y))



SASでの実装コード


SASでは treeboost プロシジャで実装できる。
任意の検証データを用いたい場合、assess validata = バリデーションデータ で指定する。

特徴量は input 特徴量 / level=尺度 で名義尺度(nominal)、順序尺度(ordinal)、間隔尺度(interval)を指定する。 何も指定しない場合、文字型は名義尺度、数値型は間隔尺度と見なされる。
本プロシジャでは Huber型 M-推定 を行うようで極端な外れ値の影響を受けにくい。

オプションはたくさんあるので、詳細は プロシジャガイド を参照。

proc treeboost data=train(where=(selected=0))
                iterations = 1000 /* pythonではn_estimators */
                maxdepth = 6      /* pythonではmax_depth */
                seed = 19
                ; 
 
    target 目的変数 / level=binary;
    input 名義尺度の特徴量 / level=nominal;
    input 順序尺度の特徴量 / level=ordinal order=ascending;
    input 間隔尺度の特徴量 / level=interval;
    
    assess validata = valid_dataset;
    subseries best;
    
    code file="gradboost.sas";
    save fit=_gb_fit importance=_gb_importance;
quit;
run;


SASでお試し(対象データ、データ加工)


今回もKaggleチュートリアルタイタニックデータを用いる。

データ加工は以前と同じなのでコードをたたんでいる。
【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰(LASSO) - こちにぃるの日記

★コードを見る場合ここをクリック★

/*---------------------------------------------------
    欠損値処理、変数作成
---------------------------------------------------*/
*-- 数値型:欠損値を中央値に置き換え;
proc stdize data=train out=train method=median reponly;
    var age;
run;

*-- カテゴリ型:欠損値を最頻値に置き換え;
data train;
    set train;
    Embarked = coalescec(Embarked,"S");
run;

/*---------------------------------------------------
    データ分割
---------------------------------------------------*/
proc surveyselect data=train rate=0.2 out=train outall method=srs seed=19;
run;
data train(drop=selected) test(drop=selected);
    set train;
    if selected=1 then output test; else output train;
run;

proc surveyselect data=train rate=0.2 out=train outall method=srs seed=19;
run;

proc format;
    value rollf 0="train-data" 1="validation-data" 2="test-data";
run;
data alldata;
   set train test(in=in1);
   format selected rollf.; *-- フォーマット適用;
   if in1 then selected=2; *-- テストデータに番号振る;
run;


SASでお試し(モデル作成)


Gradient Boosting Tree(勾配ブースティング木)の実行。

proc treeboost data=train(where=(selected=0))
                iterations = 1000 /* pythonではn_estimators */
                maxdepth = 6      /* pythonではmax_depth */
                seed = 19
                ; 
 
    target Survived / level=binary;
    input Sex Embarked / level=nominal;
    input Pclass SibSp Parch / level=ordinal order=ascending;
    input Age Fare / level=interval;
    
    assess validata = train(where=(selected=1));
    subseries best;
    
    code file="gradboost.sas";
    save fit=_gb_fit importance=_gb_importance;
quit;
run;


変数重要度を図示する。

title "Feature Importance of validation dataset";
proc sgplot data=_gb_importance;
    hbar name /response=vimportance  groupdisplay=cluster categoryorder=respdesc;
run;

f:id:cochineal19:20210725010400p:plain:w400

作成モデルで予測する。

%macro _Pred(cd=);
    data _Pred_&cd.;
        set alldata;
        %inc &cd.;
        if P_Survived1 > 0.5 then U_Survived = 1; else U_Survived = 0;
    run;
%mend _Pred;
%_Pred(cd=gradboost);


SASでお試し(モデル評価)


評価指標を自作し、スコアを算出する。

%macro _Scores(inds=, tlabel=Survived, pred=U_Survived, group=selected);
    proc sql;
        create table &inds._Scores as
            select
                 %if &group.^="" %then %do; &group. , %end;
                 count(*) as N
                ,sum(case when &tlabel.=1 and &pred.=1 then 1 else 0 end) as TP
                ,sum(case when &tlabel.=1 and &pred.=0 then 1 else 0 end) as FN
                ,sum(case when &tlabel.=0 and &pred.=1 then 1 else 0 end) as FP
                ,sum(case when &tlabel.=0 and &pred.=0 then 1 else 0 end) as TN
                ,(calculated TP + calculated TN) / count(*)      as Accuracy    label="正診率(Accuracy)"
                ,calculated TP / (calculated TP + calculated FN) as Recall      label="検出率(Recall) 、類語:感度"
                ,calculated TN / (calculated FP + calculated TN) as Specificity label="特異度(Specificity)"
                ,calculated TP / (calculated TP + calculated FP) as Precision   label="適合率(Precision)"
                ,2 * calculated Recall * calculated Precision / (calculated Recall + calculated Precision) as f1_score label="F1スコア"
            from &inds.
             %if &group.^="" %then %do; group by &group. %end;
        ;
    quit;
    
    title "&inds.";
    proc print data=&inds._Scores; run;
    title "";
%mend _Scores;
%_Scores(inds=_Pred_gradboost);

f:id:cochineal19:20210725010328p:plain

今回の f1スコアは 0.754 くらい。
この前、ランダムフォレストが 0.709 くらいだったのでだいぶ良い予測性能。
【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰(LASSO) - こちにぃるの日記

※毎度のことハイパーパラメータチューニングはしていないのであしからず。

参考


https://support.sas.com/documentation/solutions/emtmsas/93/emprcref.pdf

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。