【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

【SAS】hpsvmプロシジャ_サポートベクターマシーン

今回はサポートベクターマシーンを試してみる。


Pythonでの実装コード


Pythonで簡単に実装できる。

#-- 線形SVM
from sklearn.svm import LinearSVC
mod1 = LinearSVC(C=1, random_state=999)
mod1.fit(X_train, y_train)
y_pred1 = mod1.predict(X_test)
 
#-- 非線形SVM
from sklearn.svm import SVC
klist=['poly','rbf','sigmoid']
for k in klist:
    mod1 = SVC(C=1, kernel=k, random_state=999)
    mod1.fit(X_train, y_train)
    y_pred1 = mod1.predict(X_test)


以前、次の記事でも触れている。
【機械学習_Python】サポートベクターマシーン - こちにぃるの日記

SASでの実装コード


ものすごくシンプルな記載。
ここにカーネルやらペナルティやらクロスバリデーションやらを設定できる。

proc hpsvm data = 訓練データ;
    target 目的変数 / order = 順序 level = 変数のタイプ;
    input 特徴量1, ... , 特徴量p;
    code file="hpsvm_fit.sas"; /* 将来データ用に予測アルゴリズム出力 */
run;
 
data _Pred_&cd.;
    set alldata;
    %inc hpsvm_fit; *-- 新しいデータで予測;
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でお試し(モデル作成)


線形SVMの実行。
kernel linear で線形SVM(デフォルトなので書かなくても良い)。
partition で訓練データとバリデーションデータを指定。
select cv=validateset でバリデーションデータを指定。クロスバリデーションで penalty のハイパーパラメータCを決定する。例では0.1~1の間で最適値を探してくれる。

proc hpsvm data=train;
    partition roleVar=selected(train='0' valid='1');
    performance details;
    
    target Survived / order=desc level=nominal;
    input Pclass Sex Age Fare Embarked SibSp Parch;
    
    kernel linear;
    penalty C=0.1 to 1;
    select cv=validateset seed=19;
    
    code file="hpsvm_linear_cv.sas";
run;


線形SVM(10-fold cross validation)の実行。
select cv=random fold=10 で10-fold クロスバリデーション。

proc hpsvm data=train;
    performance details;
    
    target Survived / order=desc level=nominal;
    input Pclass Sex Age Fare Embarked SibSp Parch;
    
    kernel linear;
    penalty C=0.1 to 1;
    select cv=random fold=10 seed=19;
    
    code file="hpsvm_linear_10fcv.sas";
run;


多項式カーネルでの非線形SVM(10-fold cross validation)の実行。
kernel polynom多項式カーネルを指定。
その他、rbf(Gaussカーネル)、sigmoid(シグモイドカーネル)がある。

proc hpsvm data=train;
    performance details;
    
    target Survived / order=desc level=nominal;
    input Pclass Sex Age Fare Embarked SibSp Parch;
    
    kernel polynom / degree=2;
    penalty C=0.1 to 1;
    select cv=random fold=10 seed=19;
    
    code file="hpsvm_polynom_10fcv.sas";
run;


作成モデルで予測する。

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


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    /*正診率(Accuracy)*/
                ,calculated TP / (calculated TP + calculated FN) as Recall      /*検出率(Recall) 、類語:感度*/
                ,calculated TN / (calculated FP + calculated TN) as Specificity /*特異度(Specificity)*/
                ,calculated TP / (calculated TP + calculated FP) as Precision   /*適合率(Precision)*/
                ,2 * calculated Recall * calculated Precision / (calculated Recall + calculated Precision) as f1_score /*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_linear_cv);
%_Scores(inds=_Pred_linear_10fcv);
%_Scores(inds=_Pred_polynom_10fcv);


f:id:cochineal19:20210713230951p:plain

今回のいずれのモデルもハイパーパラメータ Cは 0.1 が最適値だったよう。
f1スコアは多項式カーネル0.70 くらい。

参考


SAS Help Center

【SAS】hpforestプロシジャ_ランダムフォレスト

今回はランダムフォレストを試してみる。


Pythonでの実装コード


Pythonでは非常に簡単に実装できる。

from sklearn.ensemble import RandomForestClassifier
mod1 = RandomForestClassifier(n_estimators = 100
                              ,criterion = "gini"
                              ,max_depth = i_integer
                              ,max_features = "auto"
                              ,bootstrap = True
                              ,random_state = 1)


クロスバリデーションも簡単(例は決定木とランダムフォレストをグリッドサーチで実装)。

import pandas as pd
from sklearn import model_selection
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix
 
modfnc1 = [DecisionTreeClassifier(), RandomForestClassifier()]
prm1 = [{"max_depth": stats.randint(1,6)
         ,"criterion":["gini","entropy"]
         ,"random_state":[1]}
        ,{"n_estimators": stats.randint(10,201)
         ,"max_depth": stats.randint(1,6)
         ,"criterion":["gini","entropy"]
         ,"max_features":["auto", "sqrt", "log2"]
         ,"random_state":[1]}]
for modfnc1, prm1 in zip(modfnc1, prm1):
    clf = model_selection.RandomizedSearchCV(modfnc1, prm1)
    clf.fit(X_train, y_train)
    y_pred1 = clf.predict(X_test)
    print(clf.best_estimator_)
    print(confusion_matrix(y_test, y_pred1))
    print(classification_report(y_test, y_pred1))
    print()


SASでの実装コード


ホールドアウトでのモデル作成は次の通り。
任意の訓練データと検証データを用いたい場合、optionに scoreprole = valid を付け、partition statement で指定する。

proc hpforest data = 訓練データ seed = 19 scoreprole = valid
    maxtrees = 100 /* pythonではn_estimators */
    maxdepth = 6   /* pythonではmax_depth */
    ;
 
    partition roleVar=selected(train='0' valid='1');
    performance details;
 
    target 目的変数 / level = 変数のタイプ;
    input 特徴量1, ... , 特徴量p;
 
    ods output FitStatistics      = _Res_FitStatistics;      /* 予測値の出力 */
    ods output VariableImportance = _Res_VariableImportance; /* 変数重要度(特徴量重要度)の出力 */
    save file = "hpforest_fit.bin"; /* 将来データ用に予測アルゴリズム出力 */
run;


SASに任せる場合は trainfraction=訓練データの割合 とすればプロシジャ内でデータ分割してくれる。

proc hpforest data = 訓練データ seed = 19 trainfraction=0.8
    maxtrees = 100 /* pythonではn_estimators */
    maxdepth = 6   /* pythonではmax_depth */
    ;
    ...
run;


なお、hpforest プロシジャでは save file = "出力ファイル名.bin" を付けることでモデルのアルゴリズムを出力できる。
この出力ファイルを hp4score に渡せば、新しいデータで予測ができる。

proc hp4score data = alldata;
    score file="hpforest_fit.bin" out=_Pred_RandomForest;
    id PassengerId selected;
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でお試し(モデル作成)


ランダムフォレストの実行。

proc hpforest data=train seed=19 scoreprole=valid
    maxtrees    = 100 /* pythonではn_estimators */
    maxdepth    = 6   /* pythonではmax_depth */
    vars_to_try = 4 
    ;
 
    partition roleVar=selected(train='0' valid='1');
    performance details;
 
    target Survived / level = binary;
    input Pclass Sex Age Fare Embarked SibSp Parch;
 
    ods output FitStatistics      = _Res_FitStatistics;
    ods output VariableImportance = _Res_VariableImportance;
    save file="hpforest_fit.bin"; 
run;

f:id:cochineal19:20210712224857p:plain:w400

変数重要度を図示する。

title "Feature Importance";
proc sgplot data=_Res_VariableImportance;
    hbar Variable /response=Gini  groupdisplay=cluster categoryorder=respdesc;
run;

f:id:cochineal19:20210712224944p:plain:w400

作成モデルで予測する。

proc hp4score data = alldata;
    score file="hpforest_fit.bin" out=_Pred_RandomForest;
    id PassengerId selected;
run;
 
data _Pred_RandomForest;
    set _Pred_RandomForest;
    if P_Survived1 > 0.5 then U_Survived = 1; else U_Survived = 0;
run;


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    /*正診率(Accuracy)*/
                ,calculated TP / (calculated TP + calculated FN) as Recall      /*検出率(Recall) 、類語:感度*/
                ,calculated TN / (calculated FP + calculated TN) as Specificity /*特異度(Specificity)*/
                ,calculated TP / (calculated TP + calculated FP) as Precision   /*適合率(Precision)*/
                ,2 * calculated Recall * calculated Precision / (calculated Recall + calculated Precision) as f1_score /*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_RandomForest);

f:id:cochineal19:20210712225319p:plain

今回の f1スコアは 0.709 くらい。
この前 hpgenselect プロシジャで作成したモデル(Lasso ロジスティック回帰)よりは良い予測性能。
【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰(LASSO) - こちにぃるの日記

ハイパーパラメータはあまりいじってないので、ここをいじれば性能は良くなるかもしれない。
なお、forest プロシジャではハイパーパラメータを自動調整してくれる autoturn statement があるようだが、hpforest プロシジャに同様の機能があるかは調べられていない。
SAS Help Center

参考


SAS Help Center
SAS Help Center
https://www.lexjansen.com/mwsug/2016/AA/MWSUG-2016-AA20.pdf
https://support.sas.com/content/dam/SAS/support/en/books/free-books/machine-learning-with-sas-special-collection.pdf

【SAS】glmselectプロシジャ_Lasso回帰、ElasticNet回帰

前回 hpgenselect プロシジャで罰則付きロジスティック回帰を実行してみた。
【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰(LASSO) - こちにぃるの日記

目的変数が連続値の時のLASSO回帰をやってなかったのでまとめてみる。
ついでにElasticNet回帰もやってみる。

【目次】


使用データ


SASHELP にあるベースボールデータを使用する。選手のヒット数、エラー数などの成績から年俸(対数)を予測してみる。

*-- データ;
data ads;
    set sashelp.baseball;
    if cmiss(of n:, of Cr:, logSalary, Division, League) = 0;
    if Division="West"   then Division_West  =1; else Division_West  =0; *One-hot encoding;
    if League="American" then League_American=1; else League_American=0; *One-hot encoding;
run;


データ分割


ベースボールデータの2割をテストデータとする。
ホールドアウト法用に残り8割を 7 : 3 = 訓練データ : バリデーションデータ に細分する。

*-- データ分割;
proc surveyselect data=ads rate=0.2 out=ads outall method=srs seed=19;
run;
data train(drop=selected) test(drop=selected);
    set ads;
    if selected=1 then output test; else output train;
run;
 
proc surveyselect data=train rate=0.3 out=train outall method=srs seed=19;
run;
 
proc format;
    value rollf 0="train-data" 1="validation-data" 2="test-data";
run;
data alldata;
    format selected selected2 rollf.; *-- フォーマット適用;
    set train test(in=in1);
    if in1 then do;
         selected =2; *-- テストデータに番号振る;
         selected2=2;
    end;
    else selected2=0;
run;


LASSO回帰 Hold-out法


ホールドアウト法においては、partition statement で訓練データと検証データを指定できる。
なお、model statement の option に stb を入れると偏回帰係数を出してくれる。

*-- ラッソ回帰;
ods graphics on;
proc glmselect data=train plots=all seed=19;
    partition roleVar=selected(train='0' validate='1');
    model logSalary = nAtBat nHits nHome nRuns nRBI nBB nOuts nAssts nError
                      CrAtBat CrHits CrHome CrRuns CrRbi CrBB
                      Division_West League_American
                      / details=all stats=all
                        selection=lasso(choose=validate stop=none) stb;
    code file="glmselect_lasso.sas";
run;
ods graphics off;
 
data pred_lasso;
    format selected selected2 logSalary P_logSalary;
    set alldata;
    %include glmselect_lasso;
run;

f:id:cochineal19:20210711162118p:plain:w300

LASSO回帰 Cross-validation法


クロスバリデーションでは model statementselection=lasso(choose=cv) cvmethod=random(10) とする。
random(k) で k-fold cross-validation になる。なお、seedも指定可能。

*-- ラッソ回帰 CV;
ods graphics on;
proc glmselect data=train plots=all seed=19;
    model logSalary = nAtBat nHits nHome nRuns nRBI nBB nOuts nAssts nError
                      CrAtBat CrHits CrHome CrRuns CrRbi CrBB
                      Division_West League_American
                      / details=all stats=all
                        selection=lasso(choose=cv stop=none) cvmethod=random(10) stb;
    score out=score_lasso_cv predicted residual;
    code file="glmselect_lasso_cv.sas";
run;
ods graphics off;
 
data pred_lasso_cv;
    format selected2 logSalary P_logSalary;
    set alldata;
    %include glmselect_lasso_cv;
run;

f:id:cochineal19:20210711162331p:plain:w300

ElasticNet回帰 Hold-out法


ElasticNetもLassoと基本は同じ。
L2正則化項の係数を指定でき、グリッドサーチも可能。selection=elasticnet(choose=validate l2search=grid showl2search) とする。
showl2search を入れると最終決定値を表示できる。

もちろん、係数を直接指定することもできる。例:selection=elasticnet(choose=validate l2=0.5)

*-- ElasticNet回帰 ;
ods graphics on;
proc glmselect data=train plots=all seed=19;
    partition roleVar=selected(train='0' validate='1');
    model logSalary = nAtBat nHits nHome nRuns nRBI nBB nOuts nAssts nError
                      CrAtBat CrHits CrHome CrRuns CrRbi CrBB
                      Division_West League_American
                      / details=all stats=all
                        selection=elasticnet(choose=validate l2search=grid showl2search stop=none) stb ;
    score out=score_elasticnet predicted residual;
    code  file="glmselect_elasticnet_l2grid.sas";
run;
ods graphics off;
 
data pred_elasticnet_l2grid;
    format selected selected2 logSalary P_logSalary;
    set alldata;
    %include glmselect_elasticnet_l2grid;
run;

f:id:cochineal19:20210711163023p:plain:w250 f:id:cochineal19:20210711163109p:plain:w300

ElasticNet回帰 Cross-Validation法


クロスバリデーションについても Lassoと同じ。

*-- ElasticNet回帰;
ods graphics on;
proc glmselect data=train plots=all seed=19;
    model logSalary = nAtBat nHits nHome nRuns nRBI nBB nOuts nAssts nError
                      CrAtBat CrHits CrHome CrRuns CrRbi CrBB
                      Division_West League_American
                      / details=all stats=all
                        selection=elasticnet(choose=cv l2search=grid showl2search stop=none) cvmethod=random(10) stb;
    score out=score_elasticnet predicted residual;
    code  file="glmselect_elasticnet_cv_l2grid.sas";
run;
ods graphics off;
 
data pred_elasticnet_cv_l2grid;
    format selected selected2 logSalary P_logSalary;
    set alldata;
    %include glmselect_elasticnet_cv_l2grid;
run;

f:id:cochineal19:20210711163316p:plain:w250 f:id:cochineal19:20210711163352p:plain:w300

今回は L2=0 でLasso回帰と等価になった。

モデル評価


モデル評価指標を複数並べたいので自作する。

*-- 評価指標;
%macro _Scores(inds=, act=logSalary, pred=P_logSalary, group=);
    proc sql;
        create table &inds._Scores as
            select
                 %if &group.^="" %then %do; &group. , %end;
                 count(*) as N
                ,sum(abs(&pred. - &act.)) / count(*) as MAE label="平均絶対誤差(MAE:Mean Absolute Error)"
                ,sum((&pred. - &act.)**2) / count(*) as MSE label="平均二乗誤差(MSE:Mean Squared Error)"
                ,sqrt(calculated MSE)                as RMSE label="MSEの平方根(RMSE:Root Mean Squared Error)"
                ,sum(abs((&pred. - &act.) / &act.)) / count(*) as MAPE label="平均絶対誤差(MAE:Mean Absolute Error)"
                ,sum(((&pred. - &act.) / &act.)**2) / count(*) as MSPE label="平均二乗パーセント誤差(MSPE:Mean Squared Percentage Error)"
                ,sqrt(calculated MSPE) as RMSE as RMSPE label="平均二乗パーセント誤差の平方根(RMSPE:Root Mean Squared Percentage Error)"
            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_lasso,                group=selected);
%_Scores(inds=pred_lasso_cv,             group=selected2);
%_Scores(inds=pred_elasticnet_l2grid,    group=selected);
%_Scores(inds=pred_elasticnet_cv_l2grid, group=selected2);

f:id:cochineal19:20210711163511p:plain:w500

RMSEで見ると今回はクロスバリデーションで作ったモデルが良好な予測性能。

参考


https://support.sas.com/documentation/onlinedoc/stat/131/glmselect.pdf

【SAS】hpgenselect プロシジャ_罰則付きロジスティック回帰(LASSO)

特徴量選択でSTEPWISE法などのほか、LASSO(L1正則化)を視野に入れる人もいる(だろう)。
LASSOについては以前 Python の記事で触れた。
cochineal19.hatenablog.com
今回は「SASで使ったみたい」&「ロジスティック回帰で使ってみたい」という視点でまとめる。(需要は低いかもしれないが、)SAS機械学習に踏み込める!という点は面白い。

【目次】


Pythonでの実装コード


Pythonでは sklearn でめちゃくちゃ簡単に実装できる。
使用するプログラム言語に制約がないならPythonが一番使いやすい。

from sklearn.linear_model import LogisticRegression
model1 = LogisticRegression(C=lambda, penalty="l1", solver='liblinear')
model1.fit(train_X, train_y)
model1.predict(valid_X)


Rでの実装コード


Rでも glmnet で簡単に実行できる。
RではElasticNetの土台があり、パラメタ調整で LassoElasticNetRidge を使い分ける。

library(glmnet)

vals <- c("特徴量1","特徴量2",...,"特徴量p")
train_X <- data.matrix(train[vals])
train_y <- train$目的変数
valid_X <- data.matrix(valid[vals])
valid_y <- valid$目的変数
 
model1 <- glmnet(train_X, train_y, family=binomial(link=logit), alpha=1, lambda=lamda)
coef(model1)
pred1 <- predict.glmnet(model1, valid_X, type="response")
 
    #-- alpha:  1=Lasso, 0=Ridge、0~1の間=ElasticNet
    #-- lambda: PythonのCと同じ。
    #--    lamda=cv.lasso$lambda.min でMSE最小のLassoモデル
    #--    lamda=cv.lasso$lambda.1se でSE最小のLassoモデル


SASでの実装コード


ここから本題のSASでの実装。
難しいのかなと思ったが、hpgenselect プロシジャで比較的簡単に実装できる。

*-- モデル作成;
proc hpgenselect data=train lassorho=0.8 lassosteps=20;
    partition roleVar=selected(train='0' validate='1');
    model 目的変数(event="1")=特徴量1, ..., 特徴量p / distribution=binary link=logit;
    selection method=lasso(choose=validatestop=none) details=all;
    performance details;
    code file="MyCode.sas" residual; *-- 予測用のPGM生成;
    ods output ParameterEstimates=_Coef_lasso; *-- 回帰係数;
    output data=_Pred_lasso role=ROLE pred=PRED resid=RESID;
run;
 
*-- 構築モデルで予測(hpgenselectプロシジャのcodeステートメントで出力したPGMを実行);
data test_pred;
   set test;
   %include MyCode;
run;


hpgenselect プロシジャは、バリデーションデータを用いた評価に対応でき、partition statement で指定する。

また、新しいデータで予測したい場合、code statement で作成モデルのPGMコード(SASファイル)を出力し、それをdataステップで実行( %include )すれば良い。
この構成は Python, R と異なり面白い。

なお、hpgenselect のガイド を見ると Model statement の option に NOCENTER があり、Requests that continuous main effects not be centered and scaled とあるので、おそらく内部処理で標準化してくれている。

もし標準化できていないのなら stdize プロシジャで標準化することもできる。

proc stdize data=train out=train2;
    var 特徴量;
run;


SASでお試し(前処理)


お試しに、Kaggleチュートリアルタイタニックデータで分析してみる。

  • 目的変数: Survived(死亡有無)

  • 特徴量: Pclass(乗客の階級)Sex(性別)Age(年齢)Fare(乗船料金)Embarked(乗船港)SibSp(兄弟・配偶者数)Parch(両親・子供の数)


欠損値処理は 前の記事 と同じ(前回から SibSpParch 追加)。
stdize プロシジャでゼロ、中央値、平均値など色々埋めれるようで便利。

*-- 数値型:欠損値を中央値に置き換え;
proc stdize data=train out=train method=median reponly;
    var age;
run;
 
*-- カテゴリ型:欠損値を最頻値に置き換え;
data train;
    set train;
    Embarked=coalescec(Embarked, "S");
run;


データ分割は surveyselect プロシジャで実装できる。シード指定も簡単。 今回は訓練データ、バリデーションデータ、テストデータに分けてみる。

*-- 訓練データとテストデータに 8:2 で分割;
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;
 
*-- 訓練データを訓練:バリデーション=8:2 に再分割;
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でお試し(LASSO)


お待ちかねのLASSO。

*-- ラッソを使ったロジスティック回帰;
proc hpgenselect data=train;
    partition roleVar=selected(train='0' validate='1');
    class Pclass(ref="1") Sex(ref="male") Embarked(ref="S") / param=ref;
    model Survived(event="1")=Pclass Sex Age Fare Embarked SibSp Parch / distribution=binary link=logit;
    selection method=lasso(choose=validate stop=none) details=all;
    performance details;
    code file="hpgenselect_lasso.sas" residual; *-- 予測用のPGM生成;
    *ods output ParameterEstimates=_Coef_lasso; *-- 回帰係数;
    *output data=_Pred_lasso role=ROLE pred=PRED resid=RESID;
run;
data Prediction_lasso;
   set alldata;
   %include hpgenselect_lasso; *-- 予測用PGM実行(構築モデルで予測);
run;


特徴量は Sex(性別)Age(年齢)Fare(乗船料金) が残ったみたい。

f:id:cochineal19:20210710160738p:plain:w250

SASでお試し(STEPWISE)


hpgenselect はSTEPWISE法にも対応しているので、比較がてら実行してみる。

proc hpgenselect data=train;
    partition roleVar=selected(train='0' validate='1');
    class Pclass(ref="1") Sex(ref="male") Embarked(ref="S") / param=ref; 
    model Survived(event="1") = Pclass Sex Age Fare Embarked SibSp Parch SibSp*Parch / distribution=binary link=logit;
    selection method=stepwise(slentry=0.05 slstay=0.05);
    performance details;
    code file="hpgenselect_stepwise.sas" residual; *-- 予測用のPGM生成;
    *ods output ParameterEstimates=_Coef_stepwise; *-- 回帰係数;
    *output data=_Pred_stepwise role=ROLE pred=PRED resid=RESID;
run;
data Prediction_stepwise;
   set alldata;
   %include hpgenselect_stepwise; *-- 予測用PGM実行(構築モデルで予測);
run;


STEPWISEでは、Pclass(乗客の階級)Sex(性別)Age(年齢) が残った。

f:id:cochineal19:20210710160928p:plain:w400

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


hpgenselect は評価指標を算出してくれないので自作してみる。
なお、二値分類問題の評価指標は、以前の記事で解説している。
【統計】ロジスティック回帰分析(ホスマー・レメショウ検定、ROC、AUC等) - こちにぃるの日記

%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    /*正診率(Accuracy)*/
                ,calculated TP / (calculated TP + calculated FN) as Recall      /*検出率(Recall) 、類語:感度*/
                ,calculated TN / (calculated FP + calculated TN) as Specificity /*特異度(Specificity)*/
                ,calculated TP / (calculated TP + calculated FP) as Precision   /*適合率(Precision)*/
                ,2 * calculated Recall * calculated Precision / (calculated Recall + calculated Precision) as f1_score /*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_lasso);
%_Scores(inds=_Pred_stepwise);


今回、F1スコアで評価すると、わずかだが LASSOでの構築モデルが良好だった。
f:id:cochineal19:20210710161015p:plain

なお、正則化に関わるLASSOのパラメタ(lassorho lassosteps)はデフォルトのままでいじってない。ここをいじるとさらに良くなるかもしれない。
グリッドサーチやランダムサーチが簡単にできるとありがたいが、どうなのだろうか。

参考


https://support.sas.com/documentation/onlinedoc/stat/142/hpgenselect.pdf
https://www.mwsug.org/proceedings/2017/AA/MWSUG-2017-AA02.pdf
https://www.lexjansen.com/sesug/2016/SD-253_Final_PDF.pdf
Solved: HPGENSELECT - LASSO- LOGISTIC - SAS Support Communities
sklearn.linear_model.LogisticRegression — scikit-learn 0.24.2 documentation
https://cran.r-project.org/web/packages/glmnet/glmnet.pdf

【SASオンデマンド】CSVファイル等を自分のPCからアップロードする

地味に困ったのでメモ(きっと他にも困る人がいるはず)。

自分のPCからアップロード


SASオンデマンド上で

  • [サーバーファイルとフォルダ] > [odaws##-XXXX##] > [ファイル(ホーム)] を選択。

  • ツールバーの [アップロード] コマンド(上矢印のアイコン)が有効になるのでクリック。
    (もしくは [ファイル(ホーム)] を右クリックして、ポップアップから [ファイルのアップロード] を選択)

  • [ファイルのアップロード] 画面が表示されるのでアップロード。


Work領域にインポート


アップロード後、次のコードでWork領域へ。
エンコーディングはデータに合わせる。例示はCSVファイル)

filename file1 '/home/ユーザー名/ファイル名.csv' encoding='utf-8';
proc import datafile=file1
            out=work.ads
            dbms=csv;
            getnames=yes;
run;


参考


Importing Data into SAS OnDemand for Academics - SAS Tutorials - LibGuides at Kent State University

【Kaggle】Titanicデータを XGBoost, LightGBM, CatBoost で分析してみる

GBDT(Gradient Boosting Decision Tree)を使ってみたかったのでメモ。
KaggleのTitanicデータを予測してみる。

scikit-learnに準拠した model.fit(データ) で記載。
あまりパラメタいじらず&Cross Validateせずでも Catboost は Score:0.78 が出た。
n_estimaternum_boost_round )はでっかい値にして、early_stopping_rounds で切った。

前処理


前回記事と同じ。Embarkedはダミー変数に変更。

import numpy as np 
import pandas as pd
from scipy import stats
from sklearn import model_selection
 
#-- import
train = pd.read_csv("../input/titanic/train.csv")
test = pd.read_csv("../input/titanic/test.csv")
###--- 欠損処理
train["Age"] = train["Age"].fillna(train["Age"].median()) #--中央値で補完
train["Embarked"] = train["Embarked"].fillna(train["Embarked"].mode().iloc(0)) #-- 最頻値で補完 
test["Age"] = test["Age"].fillna(test["Age"].median()) #--中央値で補完
test["Fare"] = test["Fare"].fillna(test["Fare"].median()) #--中央値で補完
 
###--- カテゴリの数値化
train["Sex"] = train["Sex"].map({'male':0, 'female':1})
train = pd.get_dummies(train, columns=['Embarked'])
test["Sex"] = test["Sex"].map({'male':0, 'female':1})
test = pd.get_dummies(test, columns=['Embarked'])
###--- データ分割
Xval = ["Pclass", "Sex", "Age", "Fare", "Embarked_C", "Embarked_Q", "Embarked_S"]
df_y = train["Survived"].values
df_X = train[Xval].values
train_X, valid_X, train_y, valid_y = model_selection.train_test_split(df_X, df_y, random_state=19)
test_X = test[Xval].values


XGBoost (eXtreme Gradient Boosting)


XGBoost Python Package — xgboost 1.5.0-dev documentation

import xgboost as xgb
from sklearn.metrics import accuracy_score
 
#-- モデル
xgb = xgb.XGBClassifier(max_depth = 10
                        ,learning_rate = 0.01
                        ,n_estimators = 100000
                        ,objective = "binary:logistic"
                        ,eval_metric = "logloss"
                        ,random_state = 19
                       )
 
#-- 学習
xgb.fit(train_X, train_y
        ,eval_set = [(valid_X, valid_y)]
        ,early_stopping_rounds = 100, verbose = 10
       )
 
#-- バリデーションデータで予測&正解率算出
pred = xgb.predict(valid_X)
print('score', accuracy_score(valid_y, pred))
 
#-- テストデータで予測
pred_y = xgb.predict(test_X)


Score
- Valid: 0.8565022421524664
- Test: 0.75837

LightGBM


lightgbm.LGBMClassifier — LightGBM 3.2.1.99 documentation

import lightgbm as lgb
from sklearn.metrics import accuracy_score
 
#-- モデル
lgb = lgb.LGBMClassifier(n_estimators = 100000, objective = 'binary')
 
#-- 学習
gbm.fit(train_X, train_y
        ,eval_set = [(valid_X, valid_y)]
        ,early_stopping_rounds = 100, verbose = 10
       )
 
#-- バリデーションデータで予測&正解率算出
pred = gbm.predict(valid_X, num_iteration = gbm.best_iteration_)
print('score', accuracy_score(valid_y, pred))
 
#-- テストデータで予測
pred_y = gbm.predict(test_X, num_iteration = gbm.best_iteration_)


Score
- Valid: 0.8565022421524664
- Test: 0.76315


CatBoost


Usage examples - CatBoost. Documentation

from catboost import CatBoostClassifier
from sklearn.metrics import accuracy_score
 
#-- モデル
cbc = CatBoostClassifier(num_boost_round = 10000
                         ,depth = 10
                         ,learning_rate = 0.01 
                         ,eval_metric = 'Accuracy'
                         ,loss_function = 'Logloss'
                        )
 
#-- 学習
cbc.fit(train_X, train_y
        ,eval_set = [(valid_X, valid_y)]
        ,early_stopping_rounds = 100, verbose = 10
       )
 
#-- バリデーションデータで予測&正解率算出
pred = cbc.predict(valid_X)
print('score', accuracy_score(valid_y, pred))
 
#-- テストデータで予測
pred_y = cbc.predict(test_X)


Score
- Valid: 0.8565022421524664
- Test: 0.78229


参考


【Python覚書】XGBoostで多値分類問題を解いてみる | ポテサラ
【Python覚書】アンサンブル学習:XGBoost、LightGBM、CatBoostを組み合わせる(その1) | ポテサラ
XGBoost論文を丁寧に解説する(1) - Qiita
XGBoost論文を丁寧に解説する(2): ShrinkageとSubsampling - Qiita
GBDT系の機械学習モデルであるXGBoost, LightGBM, CatBoostを動かしてみる。 - Qiita
GBDT系の機械学習モデルのパラメータチューニング奮闘記 ~ CatBoost vs LightGBM vs XGBoost vs Random Forests ~ その1 - Qiita
XGBoostのScikit-Learn APIでearly stoppingを利用する - おおかみ山
XGBoostパラメータのまとめとランダムサーチ実装 - Qiita
なぜn_estimatorsやepochsをパラメータサーチしてはいけないのか - 天色グラフィティ
Python: CatBoost を使ってみる - CUBE SUGAR CONTAINER
PythonでCatBoostの解説 – S-Analysis
Pythonでデータ分析:XGboost - データサイエンティスト(仮)
Pythonでデータ分析:Catboost - データサイエンティスト(仮)
Pythonでcatboostを使ってみる - Qiita

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