カラーリストの出力

SAS忘備録の「カラー名の取得」がとても便利

http://sas-boubi.blogspot.jp/2017/01/blog-post.html

ただ、カラー名だとピンとこないので、ちょっと一手間加えてみる。
以下のコードを実行すれば、その環境で設定されてるSASカラー名一覧が色見本ででるはず。
多分。

filename tempf temp;
proc registry list export= tempf
 usesashelp
 startat='COLORNAMES' ;
 run;

data test;
 infile tempf truncover;
 input text $ 1-100 ;
 color=compress(scan(text,1,"="),'"');
 hex=scan(text,2,"=");
 if _N_>=8;
 dummy="";
 keep color hex dummy;
 run;

proc report data=test nowd;
 column color hex dummy;
 compute dummy;
call define(_col_,'style','style={background='||color||'}');
 endcomp;
run;

以下は、2017/1/22時点のSAS on demandで実行した結果。

color	hex	dummy
AliceBlue	hex: F0,F8,FF
AntiqueWhite	hex: FA,EB,D7
Aqua	hex: 00,FF,FF
Aquamarine	hex: 7F,FD,D4
Azure	hex: F0,FF,FF
Beige	hex: F5,F5,DC
Bisque	hex: FF,E4,C4
Black	hex: 00,00,00
BlanchedAlmond	hex: FF,EB,CD
Blue	hex: 00,00,FF
BlueViolet	hex: 8A,2B,E2
BR	hex: A5,2A,2A
Brown	hex: A5,2A,2A
Burlywood	hex: DE,B8,87
CadetBlue	hex: 5F,9E,A0
Chartreuse	hex: 7F,FF,00
Chocolate	hex: D2,69,1E
Coral	hex: FF,7F,50
CornFlowerBlue	hex: 64,95,ED
Cornsilk	hex: FF,F8,DC
Crimson	hex: DC,14,3C
Cyan	hex: 00,FF,FF
DarkBlue	hex: 00,00,8B
DarkCyan	hex: 00,8B,8B
DarkGoldenrod	hex: B8,86,0B
DarkGray	hex: A9,A9,A9
DarkGreen	hex: 00,64,00
DarkGrey	hex: A9,A9,A9
DarkKhaki	hex: BD,B7,6B
DarkMagenta	hex: 8B,00,8B
DarkOliveGreen	hex: 55,6B,2F
DarkOrange	hex: FF,8C,00
DarkOrchid	hex: 99,32,CC
DarkRed	hex: 8B,00,00
DarkSalmon	hex: E9,96,7A
DarkSeaGreen	hex: 8F,BC,8F
DarkSlateBlue	hex: 48,3D,8B
DarkSlateGray	hex: 2F,4F,4F
DarkSlateGrey	hex: 2F,2F,2F
DarkTurquoise	hex: 00,CE,D1
DarkViolet	hex: 94,00,D3
DeepPink	hex: FF,14,93
DeepSkyBlue	hex: 00,BF,FF
DimGray	hex: 69,69,69
DimGrey	hex: 69,69,69
DodgerBlue	hex: 1E,90,FF
FireBrick	hex: B2,22,22
FloralWhite	hex: FF,FA,F0
ForestGreen	hex: 22,8B,22
Fuchsia	hex: FF,00,FF
G	hex: 00,80,00
Gainsboro	hex: DC,DC,DC
GhostWhite	hex: F8,F8,FF
Gold	hex: FF,D7,00
Goldenrod	hex: DA,A5,20
Gray	hex: 80,80,80
Green	hex: 00,80,00
GreenYellow	hex: AD,FF,2F
Grey	hex: 80,80,80
Honeydew	hex: F0,FF,F0
HotPink	hex: FF,69,B4
IndianRed	hex: CD,5C,5C
Indigo	hex: 4B,00,82
Ivory	hex: FF,FF,F0
Khaki	hex: F0,E6,8C
Lavender	hex: E6,E6,FA
LavenderBlush	hex: FF,F0,F5
LawnGreen	hex: 7C,FC,00
LemonChiffon	hex: FF,FA,CD
LightBlue	hex: AD,D8,E6
LightCoral	hex: F0,80,80
LightCyan	hex: E0,FF,FF
LightGoldenrodYellow	hex: FA,FA,D2
lightGray	hex: D3,D3,D3
LightGreen	hex: 90,EE,90
LightGrey	hex: D3,D3,D3
LightPink	hex: FF,B6,C1
LightSalmon	hex: FF,A0,7A
LightSeaGreen	hex: 20,B2,AA
LightSkyBlue	hex: 87,CE,FA
LightSlateGray	hex: 77,88,99
LightSlateGrey	hex: 77,88,99
LightSteelBlue	hex: B0,C4,DE
LightYellow	hex: FF,FF,E0
Lime	hex: 00,FF,00
LimeGreen	hex: 32,CD,32
Linen	hex: FA,F0,E6
Magenta	hex: FF,00,FF
Maroon	hex: 80,00,00
MediumAquamarine	hex: 66,CD,AA
MediumBlue	hex: 00,00,CD
MediumOrchid	hex: BA,55,D3
MediumPurple	hex: 93,70,DB
MediumSeaGreen	hex: 3C,B3,71
MediumSlateBlue	hex: 7B,68,EE
MediumSpringGreen	hex: 00,FA,9A
MediumTurquoise	hex: 48,D1,CC
MediumVioletRed	hex: C7,15,85
MidnightBlue	hex: 19,19,70
MintCream	hex: F5,FF,FA
MistyRose	hex: FF,E4,E1
Moccasin	hex: FF,E4,B5
NavajoWhite	hex: FF,DE,AD
Navy	hex: 00,00,80
O	hex: FF,A5,00
Oldlace	hex: FD,F5,E6
Olive	hex: 80,80,00
OliveDrab	hex: 6B,8E,23
Orange	hex: FF,A5,00
OrangeRed	hex: FF,45,00
Orchid	hex: DA,70,D6
P	hex: 80,00,80
PaleGoldenrod	hex: EE,E8,AA
PaleGreen	hex: 98,FB,98
PaleTurquoise	hex: AF,EE,EE
PaleVioletRed	hex: DB,70,93
PapayaWhip	hex: FF,EF,D5
Peachpuff	hex: FF,DA,B9
Peru	hex: CD,85,3F
Pink	hex: FF,C0,CB
Plum	hex: DD,A0,DD
PowderBlue	hex: B0,E0,E6
Purple	hex: 80,00,80
Red	hex: FF,00,00
RosyBrown	hex: BC,8F,8F
RoyalBlue	hex: 41,69,E1
SaddleBrown	hex: 8B,45,13
Salmon	hex: FA,80,72
SandyBrown	hex: F4,A4,60
SeaGreen	hex: 2E,8B,57
Seashell	hex: FF,F5,EE
Sienna	hex: A0,52,2D
Silver	hex: C0,C0,C0
SkyBlue	hex: 87,CE,EB
SlateBlue	hex: 6A,5A,CD
SlateGray	hex: 70,80,90
SlateGrey	hex: 70,80,90
Snow	hex: FF,FA,FA
SpringGreen	hex: 00,FF,7F
SteelBlue	hex: 46,82,B4
Tan	hex: D2,B4,8C
Teal	hex: 00,80,80
Thistle	hex: D8,BF,D8
Tomato	hex: FF,63,47
Turquoise	hex: 40,E0,D0
Violet	hex: EE,82,EE
Wheat	hex: F5,DE,B3
White	hex: FF,FF,FF
WhiteSmoke	hex: F5,F5,F5
Yellow	hex: FF,FF,00
YellowGreen	hex: 9A,CD,32

リファレンスにはないのに、決定木モデリングのプロットを描く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年の二人以上の世帯の「第６－１６表　１世帯当たり１か月間の日別支出」を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