①SASでデータサイエンス100本ノック（構造化データ加工編）

データサイエンティスト協会から2020年6月15日からGitHubで公開してるデータサイエンス100本ノック（構造化データ加工編）ですが

その100本をSASで書いたので，私のブログでレビューしてくださいという，非常に面白い連絡を某Mさんからいただきました

なので，気の向くままに適当にとりあげれればと

005:

【行に対する操作 / 複数条件 / 複数条件に合致する行を抽出する】

レシート明細のテーブル（df_receipt）から売上日（sales_ymd）、顧客ID（customer_id）、商品コード（product_cd）、売上金額（amount）の順に列を指定し、以下の条件を満たすデータを抽出せよ。

顧客ID（customer_id）が"CS018205000001"

売上金額（amount）が1,000以上

/*Mさま*/

data s005;
    format
        sales_ymd
        customer_id
        product_cd
        amount
    ;
    keep
        sales_ymd
        customer_id
        product_cd
        amount
    ;
    set df_receipt;
    where customer_id = 'CS018205000001' & 1000 <= amount;
run;
    
proc print data = s005;
run;

好みの問題で
私は割と
&(and)でwhere抽出条件をつなぐのが好きではなく

where customer_id = 'CS018205000001'; where same 1000 <= amount; としかかかないかなぁ．コメントアウトして条件絞りを変えやすいし

007問題の

• 顧客ID（customer_id）が"CS018205000001"
• 売上金額（amount）が1,000以上2,000以下

where customer_id = 'CS018205000001' & amount between 1000 and 2000;

整数値が前提なら amout in (1000:2000);の方が書きがちですが正しいのはbetweenの気がします

013問題の
顧客テーブル（df_customer）から、ステータスコード（status_cd）の先頭がアルファベットのA～Fで始まるデータを全項目抽出し、10件だけ表示せよ。

where substr(status_cd , 1 , 1) between 'A' and 'F';

↓

where first(status_cd , 1) between 'A' and 'F';

first関数とchark関数が好きで，substrの第三引数が1になるようなケースで私はsubstr使わないかなぁ

末尾ならreverseしてからfirstするかcharで第二引数にlengthいれるかなぁ

とりあえず，ここまで

Nオブザベーション先の値を取得する話①

Nオブザベーション前の値を保持して取得する場合，retain使ったり，lag関数使えばよいのですが(lag関数の注意点参照 https://sas-boubi.blogspot.com/2021/05/lag.html)

data A;

set sashelp.class;

pre1_age=lag(age);

pre2_age=lag2(age);

keep name age pre:;

run;

これを逆にしようと思うと，1obs読んで1obs出力するSASの仕組み的にはつらいものがあります．
もう一回sort順逆にして前のobsをとれば一応は先をとったのと同じことにはなりますが
データが巨大だったりすると，みだりにソートやステップを重ねることが難しかったりします

data B ;

set sashelp.class(drop=height weight);

_DSID = open("sashelp.class(keep=age rename=(age=_age))") ;

call missing(_age);

call set(_DSID) ;

rc = fetchobs(_DSID, _N_+1) ;

if rc ne 0 then call missing(_age);

run;

open関数でsetと同じデータセットを指定して，IDづけして
setしてる親元の_N_に+1した数で，任意のobsをブッコ抜くfetchobsでもってくる．
欲しい変数のみ残していて，リネームしているので結果的にリネーム結合と同義になる
ハッシュオブジェクトと同じくPDV外からきた変数は，PDVでの値のクリーンアップがかからないので，retainされてしまう．そこでリターンコードrcを処理にからめる
rcはfetchobsが成功，つまり1obs先がある場合のみ0なので，0以外の場合は初期化
SCL関数をデータステップ内で使った場合，終了時に自動でcloseがかかるのでclose関数は省略

これをマクロで汎用化するなら

%macro next_look(master=,var=,postn= );

if 0 then set &master(keep=&var rename=(&var=next&postn._&var));

length &master.&var.&postn. 8.;

retain &master.&var.&postn.;

if missing(&master.&var.&postn.) then do;

&master.&var.&postn. = open("&master(keep=&var rename=(&var=next&postn._&var))") ; 

call set(&master.&var.&postn.) ;

end;

rc = fetchobs(&master.&var.&postn., _N_+&postn) ;

if rc ne 0 then call missing(next&postn._&var);

drop &master.&var.&postn. rc ;

%mend;

こんなの

data C ;

set sashelp.class;

%next_look(master=sashelp.class,var=name,postn=1);

%next_look(master=sashelp.class,var=age,postn= 1);

%next_look(master=sashelp.class,var=age,postn= 2);

keep name age next:;

run;

でも，優秀な後輩からの指摘で，それなら単純にpoint setの方が速度でるんじゃねえっすか？
ってことで

data z;

set sashelp.class;

__n = _n_ + 1;

if 1 <= __n <= nof then set sashelp.class(keep = age rename = (age = next_age)) point = __n nobs = nof;

else call missing(next_age);

keep name age next_age;

run;

確かにねー．

で，巨大なデータにしてみたり，色んなシチュエーションで実験してみたけど
fetchobsもパフォーマンス悪くはないけど，常にpoint setの方がややいい感じ．

簡素性や，SCL関数とか下手に使ってない観点からしてもset pointの方がいいですね

いやー，できる後輩がいると辛いですね(笑)

次回は，いつになるかわからないですが，条件付きで先を見る話です
（○○時間経過以内の特定の値をとるとか，そのobsの先に特定条件を満たすものがあった場合のみ取得する　みたいなものには単純な割り当て法だとつらい）

通常のMergeとDS2のMergeは別物という話

当初DS2プロシジャがリリースされた時，なんと文法からMergeステートメントが廃されてました．データの結合は，内部埋め込みできるFedSQLやHash Packageでも使いなさいという思い切った仕様でした

個人的には，左側外部結合の多い医薬分野の自分のプログラムの範囲では，SASにMergeステートメントは必要ない，基本ハッシュオブジェクトがあればいい，SQLの方が性にあるというならjoinする人がいても良いという過激思想なので(人に押し付けるつもりはないですが…笑)

当時は思い切ったことをしたものだと感心しました．
まあ，DS2の場合，基本的に内部処理のベースがSQLで行われているようで，byの挙動などもSQLのorder by,group byに準じていたりするので，SAS固有の特殊な挙動をするMergeステートメントを親和させれたなったんだろうなぁと推察してました

ただ，やはりちょっと思い切り過ぎたのでしょうかね，メンテナンスリリース3あたりで，突如DS2にもMergeが搭載されました．旧来のユーザーからの要望？

ただ，やっぱり，本来仕組みの違うものに，組み込んだので，ちょっと注意が必要です．
まあ，基本的なケースでは同じ結果が返るので，使わない方がとまでは言わないのですが，知った上で気を付けて使ってくださいって話です

data T1;

length K 8. A B $20.;

K=101;A="Cat";B="Apple";output;

K=101;A="Pig";B="Banana";output;

K=102;A="Duck";B="Fig";output;

run;

data T2;

length K 8. A C $20.;

K=101;A="Cow";C="Red";output;

K=102;A="";C="Yellow";output;

K=103;A="";C="Blue";output;

run;

旧来のデータステップのMergeであれば

data A1;

merge T1 T2;

by K;

run;

後続のデータが，先行の共通で変数バッティングするところを上書いていく，いかにもSASのマージっぽい結果ですね

一方，DS2のmergeはというと

proc DS2;

data A2 (overwrite=yes);

   method run();

      merge T1 T2;

      by K;

   end;

enddata;

run;

quit;

（ライブラリやSODの出力ビューが自動更新されないので，製品SASの場合はライブラリを切り替えて戻ってきて見るか，SODならライブラリ開いて確認できます）

ほぼ同じと見せかけて，Cの2行目が欠損なことに注目です．
T2に101のデータが1つしかないからこうなってます

SASの場合はkey一致分，全てに値がつくのですが，見方を変えると同key 前obsから値がretainされているとも解釈できます

そして，DS2のM5でMergeにretainオプションが追加され，これを使うと

proc DS2;

data A3 (overwrite=yes);

   method run();

      merge T1 T2 / retain;

      by K;

   end;

enddata;

run;

quit;

おっ，今度はCの2行目に値がはいった！
しかし実はよく見ると，これも旧来のデータステップの結果と違うのです
Aの2行目，1行目の値がそのままretainされています

というわけで，今までのmergeと同じ出力結果をだすのは中々に難しいのですね
挑戦してみると面白いかもですが，SQLで今回のデータステップのmergeと同じ結果だすの結構難しいはず

ただ，まあそれぞれの結果のどれが正しいという話ではないのですね．旧来の結果に見慣れてるので，DS2のマージはおかしいって言いたくなるのですが，それもあくまで仕様どおりにでてるだけということですね

RMSTのグラフを描いてみた話

とても反抗的で手の付けられない後輩が，SASユーザー総会でRMSTの発表をしてました．
RMSTオプションは良いなぁと思うのですが，ods gprahicsででてくるKM plotは特にいつも通りのものしかでません

Pythonのlifelinesライブラリとかだと割と簡単で

from matplotlib import pyplot as plt

from lifelines.utils import restricted_mean_survival_time

from lifelines.plotting import rmst_plot

bmt=impsas(folder,"bmt")

ix=bmt['Group']=='ALL'

T, E = bmt['T'], bmt['Status']

time_limit = 1500

kmf_exp = KaplanMeierFitter().fit(T[ix], E[ix], label='exp')

rmst_plot(kmf_exp, t=time_limit)

それっぽいものがでて，視覚的にもわかりやすいのですが，これをSASでやるのは面倒です

すでにSASoneさんの方で作図例があるのですが
http://sasonediver.blog.fc2.com/blog-entry-253.html 

annotation使わずに，描けないかなぁと思いました

Number at Risk表つくってる場合，通常，境界時間τは，そんな中途半端な時点にならずにat Riskだしてるポイントと被ることが多いので
それなら，信頼区間バンドつけるときの下限を0にしちゃえばいいんじゃねってことで

%let tau=1500;

ods output Survivalplot=SurvivalPlotData;

ods output RMST=RMST(rename=(Stratum=StratumNum));

proc lifetest data=sashelp.bmt  rmst(tau=&tau) plots=survival(atrisk=0 to 2500 by 500 );

    time T*STATUS(0);

    strata group;

run; 

data SurvivalPlotData_1;

  set SurvivalPlotData;

  if 0 then set RMST;

  if _N_=1 then do;

    declare hash h1(dataset:"RMST");

    h1.definekey("StratumNum");

    h1.definedata("Estimate");

    h1.definedone();

  end;

  retain _Survival;

  if ^missing(Survival) then _Survival=Survival;

  markerchar=" ";

  if ^missing(Censored) then hige=Censored+0.03;

  if time<=&tau then limit_time=time;

  h1.find();

  if tAtRisk eq &tau then do;

    tau=&tau;

    rmst=cats("RMST=",round(Estimate,0.1));

    ry=_Survival -0.05;

  end;

run;

proc sgplot data=SurvivalPlotData_1 noautolegend;

styleattrs datacontrastcolors=(black black black) datacolors =(blue yellow red);

  step x=time y=survival /CURVELABEL group=stratum lineattrs=(thickness=3);

  scatter x=time y=censored /NOERRORCAPS yerrorupper=hige errorbarattrs=(pattern=1 thickness=3)  markerchar=markerchar GROUP=stratum;

  band  x=limit_time upper=_survival lower=0/ group=stratum type=step transparency=0.7;

  text x=tau y=ry text=rmst/textattrs=(size=15);

  xaxistable atrisk / x=tatrisk class=stratum location=outside colorgroup=stratum valueattrs=(size=10);

  refline tau/axis=x;

  yaxis label="Survival Probability" min=0 offsetmax=0.1;

  xaxis label="Day" min=0 offsetmin=0.01 ;

run;

おー，なんかそれっぽい