リファレンスにはないのに、決定木モデリングのプロットを描くdecisiontreeステートメントが追加されている話

以下を実行してみます。
ちなみにGTLのリファレンスには今時点(9.4 fifth ed)で「decisiontree」に関する記載はないです

https://support.sas.com/documentation/cdl/en/grstatug/69747/HTML/default/viewer.htm

proc template;
define statgraph temp;
   begingraph ;
      layout region;
         decisiontree nodeid=NODEID parentid=PARENTID;
      endlayout;
   endgraph;
end;
run;

data test;
NODEID=0;PARENTID=.;output;
NODEID=1;PARENTID=0;output;
NODEID=2;PARENTID=0;output;
NODEID=3;PARENTID=2;output;
NODEID=4;PARENTID=2;output;
NODEID=5;PARENTID=3;output;
NODEID=6;PARENTID=3;output;
run;

proc sgrender data=test template=temp;
run;































うわぁ、なんかでたし。




なんで、こんなことがわかったかというと
9.4で追加された決定木モデリングを行うhpslitプロシジャがods graphics onで作ってくれるプロットがかなり面白い形をしていたので気になりました。
以下のような感じです。


proc hpsplit data=sashelp.iris;
 class Species;
 model Species=SepalLength SepalWidth PetalLength PetalWidth;
 grow entropy;
 prune costcomplexity;
run;



























そこで気になったのでods trace on;にして
2つのプロットが以下のテンプレートで作成されいることを確認し
 HPStat.HPSplit.Graphics.ZoomedClassificationTreePlot
 HPStat.HPSplit.Graphics.WholeClassificationTreePlot

proc template;
source HPStat.HPSplit.Graphics.ZoomedClassificationTreePlot;
source HPStat.HPSplit.Graphics.WholeClassificationTreePlot;
run;

として、(むちゃくちゃパラメータあって複雑なテンプレートなので割愛)

中身をみて初めてdecisiontreeプロットが追加されていることをしったわけです。


ods output WholeTreePlot=out1;
ods output ZoomedTreePlot=out2;

として、流れ込んでいるデータセットの中身と見比べると、どのパラメータが何に効いているかが推測できます。

ただ、プロシジャ内部でdynamicをつかってどのように値を渡しているかは、見えないんですが、どなたか方法知りません?

sgplotで、生存率曲線の中央値地点からX軸に線落とせとか、そういうグラフの特定点から軸に線引くのはdroplineでいけるという話

グラフで、点から両軸に線を落とすのってどうやるんですか?って聞かれました。
その人はreflineで頑張ってましたが、それだと補助線でデータの位置でとめることができないです。
普通にデータまで縦の直線と、横の直線をグラフとして引いてもいいんですが

sgplotにはいつのまにか(多分9.3)からdroplineステートメントという便利な機能があります。

なんか、sgplotにドカンドカン機能追加されてて、嬉しいんですけど、把握しきれないですね。
あと、グラフの呼称って難しいですね。字面だけみても、どういうグラフ描くものなのかわからないし、検索するときも、なんて入れたもんかってね。

たとえば

data q1;
do x=1 to 10;
 y=x**2;
 output;
end;
run;

proc sgplot data=q1;
  series x=x y=y/markers;
  dropline x=x y=y / dropto=both
      lineattrs=(color=blue  pattern=dot);
  yaxis min=0;
run;

とすると、こんな感じ。





















droplineステートメントにx値とy値を指定して、dropto=でどっちの軸に線を落とすかを指定します。
bothは両方で。

あと、まあ、よくあるのが、凝りだしたら終わりのないグラフの定番kaplan-meierで
50%から線落とせってやつ。
面倒なので、ざくっと描くとこんな感じですか?

ods output SurvivalPlot=OUT1 ;
ods output Quartiles=OUT2;
proc lifetest data=SASHELP.BMT ;
 where GROUP="AML-Low Risk";
  time T * STATUS(0);
run;
ods output close;

/*中央値をマクロ変数に*/
proc sql noprint;
select Estimate into:p50
from out2
where Percent=50;
quit;

proc sgplot data=out1 noautolegend;
/*階段 */
step x=time y=Survival /lineattrs=(thickness=3);
/*髭*/
scatter x=time y=Censored/markerattrs=(symbol=plus size=15) ;

dropline x=&p50.  y=0.5  /label="★ここが50%やで" dropto=both lineattrs=(color=blue  pattern=dot thickness=3);

/*軸*/
yaxis values=(0 to 1 by 0.2 0.5);
xaxis values=(0 to 3000 by 500 &p50.);
run;


date関数やtime関数等はデータステップとSQLだと考え方が違うから注意という話

明けましておめでとうございます。

新年早々、ちょっとしたクイズです。

まず以下のデータセットがあるとします。

data a;
do i=1 to 10000000;
output;
end;
run;

次に2つのコードを示します。仮に全く同じ時間に実行した場合において
このコードの結果、できる2つのデータセットの内容は一致するでしょうか?しないでしょうか?

/*一つ目*/
data b;
set a;
time=time();
format time time.;
run;

/*2つ目*/
proc sql  ;
create table c as
select i,time() as time  format=time.
from a;
quit;

はい、まあ、問題にするってことは一致しないんですよ。

まず一つ目のデータステップでやった方の最初らへんと最後らへんの中身を見てみましょう。





















若干時間が違いますね。
データステップの場合、関数は、1obsごとに実行されるので、まさに
そのオブザベーションを処理した時間を掴みます。
よくあんのが、凄い時間のかかる処理を深夜にバッチで回したりして、翌朝みたら、
obsの途中から日付変わってるやん!ちくしょーみたいな話ですね。

次はSQLの方ですが、以下の通り




















全部おんなじ時間が入ってます。

なんでかっていうと、SQLプロシジャではdate,time,datetime等の関数は
SQLが実行される前に、日時を取得して、それを一律で関数の結果に返すという仕様を持ってます。
(わざわざ検証してないですが、たぶんtoday,weekday,holidayとか日時依存のやつは全部な気がする)

これをデータステップと同じように、そのobsが抽出された時間にしたければproc sqlにnoconstdatetimeオプションを
つけてやればOK。

proc sql noconstdatetime ;
 create table d as
 select i,time() as time  format=time.
 from a;
quit;


あるいは、実行される全てのSQLをずっとその仕様にしたいというなら
options nosqlconstdatetime;
としてやればOKです。元に戻すには
options sqlconstdatetime;
でOK。

データステップで、全obsに同じ時間を入れたい場合は、1obs目に取得してretainするか、
直前にマクロ変数に入れておいて展開するかになると思うので、こういう処理を書きたいなら
SQLでやった方が楽かもしんないですね。



以上。

さて、このブログも2013年にスタートしたので、今年の秋で丸4年です。
よくもまあ、データステップで何年も書けるなぁって自分で思います。
書きたいことを無責任に書き殴るスタイルだから、続いてるんだと思います。

最初の頃みると、月に37回とか記事あげてるんですね、思い返すと、毎日あんまり寝ずに書き続けてたんで、もう狂ってたとしか思えないですね。

今年も今まで通り、マイペースで書いたり休んだりでいくので宜しくお願いします。

また、よろしければ、是非、SASブログ初めてみませんか?
データステップでも、統計でも、グラフでも、ネタはなんでもいいので、書いてみると案外楽しいですよ。




幾何平均の出し方がちょっとずつ増えてるけど、誰も気づきませんって!という話

9.4からメンテナンスリリースのバージョンによる機能の差がちょっと目立つようになってきた。
プログラムを他の環境に渡して、ゴリゴリ動かすことが想定されているなら
メンテナンスリリースのバージョンまで、移管先とすり合わせておいた方がよくなってきそう。


ちなみに、自分の使っているSASのメンテナンスリリースのバージョンがわからない場合は

proc product_status;
run;

を実行してみるとよいですよ。











9.4_M3と、でてるM3の部分がメンテナンスレベル3を示してます。
だから、BASEのリファレンスをググッて、参照する場合はthird editionとついているやつを見た
方が正確です。
あとSAS/STATのバージョンも重要で14.1とでていますので、BASEにない統計系のプロシジャ
を参照する場合は14.1のものを見てくださいということですね。

9.3→9.4とかの大きなバージョンアップの場合は、多くのユーザーはSASがだしている
what's newを参照してからプログラムしてると思うので(してますよね?)いいんですけど、
メンテナンスリリースは、一応上がるたびに最新に更新してアップされているとはいえ
その辺をくまなくチェックしてるユーザーは少ないかも。

最近ちょっとびっくりしたのが9.4M3で、univariateに幾何平均の指定がこっそりと増えていたこと。
なんでいまさら!しかもメンテナンスリリースレベルのアップでいきなり!
正直に言わせて頂くと…、誰も気づかねぇよ!!

9.1でgeomean関数が追加されて、横持の変数の幾何平均がだせるようになり
9.2でttestプロシジャでだせるようになり
9.4M1でsurveymeansが対応して
9.4M3でunivariateが対応して(信頼区間などはでない)…

いや、まどろっこしいから、普通にproc meansに幾何平均がらみの統計キーワード
全部追加してくださいよ


/*9.4M3から univariateでgeomeanが指定可能*/
proc univariate data=sashelp.class;
var weight;
output out=out1 geomean=geomean;
run;


/*9.4M1からsurveymeansでgomean、その他信頼区間も指定可能*/
proc surveymeans data=sashelp.class geomean gmclm gmstderr lgmclm ugmclm;
var weight;
ods output GeometricMeans=out2;
run;


/*9.2からttestでlognormalを指定して幾何平均と信頼区間がだせる*/
proc ttest data=sashelp.class dist=lognormal;
  var weight;
  ods output Statistics=out3;
run;



/*9.1からgeomean関数で幾何平均がだせる*/
proc transpose data=sashelp.class out=temp1;
var weight;
run;
data out4;
set temp1;
geomean=geomean(of COL:);
keep geomean;
run;



/*9.1以前*/
data temp2;
set sashelp.class;
         logval=log(weight);
RUN;
proc means data=temp2 noprint;
  var logval;
  output out=temp3 lclm=lcl mean=mean uclm=ucl;
run;
data out5;
  set temp3;
  lcl_gm = exp(lcl);
  geomean = exp(mean);
  ucl_gm = exp(ucl);
run;

title "univariate";
proc print data=out1 label;
run;









title "surveymeans";
proc print data=out2 label;
run;




title "ttest";
proc print data=out3 label;
run;




title "転置してからgeomean関数";
proc print data=out4 label;
run;










title "logとってからmeansかけてexp";
proc print data=out5 label;
run;



サンタクロースはいつプレゼントを買うのか

せっかくなんでクリスマスネタをやりたいと思いつつも、ツリー書いたり、雪降らせたり
とてもじゃないですが、matsuさんに太刀打ちできそうにないので、別のネタを。

自分の子供はまだ0歳なんで、あんまり関係ないのですが、世の中のお父さんサンタさん達は、いつ頃プレゼントを仕入れに行ってるのかなと疑問に思いました。

どうしましょうね。

日本にはe-statという、政府統計を誰でも自由に利用できる仕組みがあります。
そこで家計調査を選択し、2015年の二人以上の世帯の「第6-16表 1世帯当たり1か月間の日別支出」を2カ月分ほどEXCELで落としてきます。

ちなみに集計表によってはCSVで落とせたり、API叩いて落とせたりするんですが、今回はEXCELを読み込みます。

ダウンロードしたファイルを開いてみると、












わぁ、なんて素敵な神エクセルでしょう。
ユーザーが読み込んで使うかもしれない統計表をこういうレイアウトにするセンスに脱帽です。

データステップで激しく加工して、品目分類を「玩具」に絞ってプロットしてみます
(ちなみに玩具は「テレビゲーム機」「ゲームソフト等」「他の玩具」を含む分類です)

値自体は、子供のいない世帯も含む平均値なので、絶対的にはあまり意味ないです。
また、事前に予約しておいて、支払いは受け取り時だったりするケースや、クレジットの支払いがどうなっているのかなど
細かい話はおいておいて、傾向をざっくり見てみましょう。


libname ex11 "XXXXXXX\raw\a616_11.xls" mixed=yes header=no;
libname ex12 "XXXXXXX\raw\a616_12.xls" mixed=yes header=no;

data wk1;
set ex11.'6-16二人以上の世帯$'n(drop=F43 in=in1)
      ex12.'6-16二人以上の世帯$'n(in=in2)
;
where F9="玩具";
if in1 then month=11;
if in2 then month=12;

run;

proc transpose data=wk1 out=wk2;
  var F11-F42;
  by month;
run;

data wk3;
set wk2;
val=col1;
year=2015;
day=input(substr(_NAME_,2,2),?? best.)-10;
    date=mdy(month,day,year);
wd=choosec(weekday(date),'日','月','火','水','木','金','土',);
label=cats(put(day,z2.),wd);

if date ne .;
format date yymmdds10.;
keep date val label santa;
run;

title"品目分類835-837 玩具の支出";
proc sgplot data=wk3 noautolegend;
series x=date y=val;
text x=date y=val text=label;

xaxis values=("01NOV2015"d to "31DEC2015"d by 15 );
run;



































やっぱり12月に入ると、土日を中心に活発に活動してるみたいですね。
年比較すると、昔より購入時期が早くなってるのか?なども見えてくるかもですね。
予約時期はこのデータからだと見えませんが。
様式が面倒すぎてあまりやる気がおきませんが…。

データステップ内でプロシジャを動かして、その結果をキー結合するマクロ

最近思いついた、相当しょうもないアイデア。
書くのは恥だが役に…って類いのコード

例えば、以下のような3つのデータセットがあって

data A;
ID=1;SEX="M";GRP=1;output;
ID=2;SEX="F";GRP=1;output;
ID=3;SEX="M";GRP=3;output;
ID=4;SEX="F";GRP=2;output;
ID=5;SEX="F";GRP=2;output;
run;

data B;
ID=1;SEX="M";WEIGHT=75;output;
ID=2;SEX="F";WEIGHT=55;output;
ID=3;SEX="M";WEIGHT=65;output;
ID=4;SEX="F";WEIGHT=45;output;
ID=5;SEX="F";WEIGHT=85;output;
run;

data C;
GRP=1;HEIGHT=175;output;
GRP=1;HEIGHT=155;output;
GRP=1;HEIGHT=165;output;
GRP=2;HEIGHT=145;output;
GRP=2;HEIGHT=185;output;
run;

BでSEXごとに求めた要約統計量をSEXをキーにしてAに結合、
CでGRPごとに求めた要約統計量をGRPをキーにしてAに結合っていうような
処理を埋め込み系のマクロにして、1ステップ中にいくらでも入れれるようにしたいなと思いました。

ようするにこういうこと
data D;
set A;
%get_summary(dataset=B,key=SEX,var=WEIGHT)
%get_summary(dataset=C,key=GRP,var=HEIGHT)
run;

結果は(画像は途中まで)




で、肝心のマクロの中身は

%macro get_summary(dataset=,key=,var=);
%let qkey  = %sysfunc( tranwrd( %str("&key") , %str( ) , %str(",") ) );
length &var &var._SUM &var._MEAN &var._StdDev &var._MIN &var._MAX 8.;
call missing(of &var &var._SUM &var._MEAN &var._StdDev &var._MIN &var._MAX );
if _N_=1 then do;
rc=dosubl("proc summary data=&dataset nway;class &key;var &var.;
                         output out=temp sum= mean= std= min= max=/autoname");
declare hash h&sysindex.(dataset:"temp");
h&sysindex..definekey(&qkey);
h&sysindex..definedata("&var._SUM","&var._MEAN","&var._StdDev","&var._MIN"
                                         ,"&var._MAX  ");
h&sysindex..definedone();
rc=dosubl("proc delete data=temp;run;");
end;
if h&sysindex..find() ne 0 then call missing(of &var._SUM &var._MEAN &var._StdDev &var._MIN &var._MAX);
drop rc &var. ;
%mend get_summary;

知ってる人からしたら、なんだdosubl関数使うのか、そりゃできるでしょうよって馬鹿にされそう。

けど、まあ、当然に感じる人が大半かもしれないけど、ハッシュに格納する直前で、このタイミングでデータセット作って間に合うっていう感覚は個人的には斬新。

dosubl関数は、call executeと違って、データステップの途中で即時的に処理が開始できる特性があります。
なのでdosubl関数でデータセット作って、そのままハッシュに入れて、すぐ消してるんですね。小賢しい…
だから見かけ上は1ステップでも、速くはないです(mergeよりかは速いけど)。
あらかじめプロシージャ回してデータセット作ってからハッシュに格納するのと同じです。

本来これはDS2でやるか、ハッシュ内でループで計算してから繋ぐ書き方でやる処理なんですね…。

ちなみに、思いつきでパパッと書いたマクロなんで、使うならちゃんと検証してくださいね。
多分、穴はあると思います。まず、再帰ができない、ようするにset対象自身を指定できないとかね。
あと複数変数指定にも対応してないですね

あ、今書いてて思ったんですけど、この記事、SAS忘備録の
「データステップ内でプロシジャを実行する。」の焼き増しですね。
http://sas-boubi.blogspot.jp/2016/05/blog-post_17.html


DATEKEYSプロシジャを使って日付にイベントを関連付ける(なんか凄い応用効く予感)

SASにはholidayname関数、holidaycount関数、その他holidayから始まる
たくさんの関数があります。それらは日付や期間を与えることで祝日の名前や祝日の数を
返したりする機能を持っているのですが…

試しに11月11日を与えて、そこに定義されている祝日の名前を書きだしてみます

data _null_;
dt='11NOV2016'd;
do count= 1 to holidaycount(dt);
name=holidayname(dt,count);
put dt yymmdds10. +2  count= name= ;
end;
run;






と、なんか3つでてきましたね
ちなみに
VETERANS : 退役軍人の日 11月11日
VETERANSUSG : 退役軍人の日(米国政府の祝日)
VETERANSUSPS : 退役軍人の日(米国郵政公社の祝日)
らしいです。

要するに、対応しているのが米国の祝日なんですね。
日本じゃ使い道ね~、工数計算とかスケジュールをSASで作るときに使えたら楽なのに~ってずっと思ってました。
ところが9.4で試用版として導入されたproc datekeysを使うことで話は変わってきます。

proc datekeysは結構深いプロシジャのようなのですが、今回は一番単純に使ってみます

2016年の祝日を直打ちで設定(データセットからも流し込める)して
それを使って5月の祝日判定をしてみましょう

まずdatekeysプロシジャを使って、イベント日の定義情報をデータセット化します

proc datekeys;
   /*元日*/ datekeydef New_Years_Day= '1JAN2016'd / pulse=day;
   /*成人の日*/ datekeydef Coming_of_Age_Day= '11JAN2016'd / pulse=day;
   /*建国記念の日*/ datekeydef National_Foundation_Day= '11FEB2016'd / pulse=day;
   /*春分の日*/ datekeydef Vernal_Equinox_Day= '20MAR2016'd / pulse=day;
   /*振替休日*/ datekeydef Substitute_Holiday= '21MAR2016'd / pulse=day;
   /*昭和の日*/ datekeydef Showa_Day= '29APR2016'd / pulse=day;
   /*憲法記念日*/ datekeydef Constitution_Memorial_Day= '3MAY2016'd / pulse=day;
   /*みどりの日*/ datekeydef Greenery_Day= '4MAY2016'd / pulse=day;
   /*子供の日*/ datekeydef Childrens_Day= '5MAY2016'd / pulse=day;
   /*海の日*/ datekeydef Marine_Day= '18JUL2016'd / pulse=day;
   /*山の日*/ datekeydef Mountain_Day= '11AUG2016'd / pulse=day;
   /*敬老の日*/ datekeydef Respect_for_the_Aged_Day= '19SEP2016'd / pulse=day;
   /*秋分の日*/ datekeydef Autumnal_Equinox_Day= '22SEP2016'd / pulse=day;
   /*体育の日*/ datekeydef Health_Sports_Day= '10OCT2016'd / pulse=day;
   /*文化の日*/ datekeydef Culture_Day= '3NOV2016'd / pulse=day;
   /*勤労感謝の日*/ datekeydef Labour_Thanksgiving_Day= '23NOV2016'd / pulse=day;
   /*天皇誕生日*/datekeydef Emperor_Akihitos_Birthday= '23DEC2016'd / pulse=day;

   datekeydata out=MyHolidays ;
run;

次にeventds=オプションで定義データセットを指定します。ここでnodefaultsも追加しておくと
デフォルトの定義(=アメリカの祝日)を外すことができます。日米両方の祝日をだしたければnodefaultsはいりません。

options eventds=(nodefaults MyHolidays);
data MAYHOLI;
length dt 8. name $50.;
do dt='1MAY2016'd  to '31MAY2016'd ; 
call missing(name);
if holidaycount(dt) = 0 then output;
else do count= 1 to holidaycount(dt);
name=holidayname(dt,count);
output;
end;
end;
format dt yymmdds10.;
keep dt name;
run;






















というわけです(画像は途中まで)。
ちなみにイベントデータは、週単位や月単位、あるいは幅を持って設定もできますし、
リファレンス見てるともっと色々できそうです。

これ、祝日とか判定する以外にもなんか使えそうじゃないですか?

特に、開始と終了日を与えてその間の祝日回数を数えるHOLIDAYCK関数とか応用できそうじゃないですか?

頭のいい人は何か使い道を編み出して発表してください