特徴量選択で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（乗船料金） が残ったみたい。

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）はデフォルトのままでいじってない。ここをいじるとさらに良くなるかもしれない。
グリッドサーチやランダムサーチが簡単にできるとありがたいが、どうなのだろうか。