【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;
作成モデルで予測する。
%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);
今回の 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);
今回のいずれのモデルもハイパーパラメータ Cは 0.1 が最適値だったよう。
f1スコアは多項式カーネルで 0.70 くらい。
参考
【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;
変数重要度を図示する。
title "Feature Importance";
proc sgplot data=_Res_VariableImportance;
hbar Variable /response=Gini groupdisplay=cluster categoryorder=respdesc;
run;
作成モデルで予測する。
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);
今回の 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回帰もやってみる。
【目次】
- 使用データ
- データ分割
- LASSO回帰 Hold-out法
- LASSO回帰 Cross-validation法
- ElasticNet回帰 Hold-out法
- ElasticNet回帰 Cross-Validation法
- モデル評価
- 参考
使用データ
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;
LASSO回帰 Cross-validation法
クロスバリデーションでは model statement で selection=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;
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;
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;
今回は 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);
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の土台があり、パラメタ調整で Lasso ～ ElasticNet ～ Ridge を使い分ける。
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（両親・子供の数）
欠損値処理は 前の記事 と同じ（前回から SibSp と Parch 追加）。
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（乗船料金） が残ったみたい。
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（年齢） が残った。
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での構築モデルが良好だった。
なお、正則化に関わる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_estimater（ num_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を組み合わせる（その１） | ポテサラ
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