前回に引き続き、バイオリンプロットを作図方法を検証します。
バイオリンプロットは箱ひげ図とともに作図されることがあります。
RのggplotやのSeabornではviolinplotメゾッドのオプションで追加表示させることができますが、GTLの場合は素直にboxplotステートメントを
追加することで同様な作図が実施可能です。
ただしSAS公式ブログの要はdatalatticeレイアウトを使う場合は、boxplotステートメントは併用不可のため、
手動で統計量を算出して、箱ひげ図の各パーツをpolygonplotステートメントとvectorplotステートメントで描画する必要があります。非常にめんどくさいのお勧めできません。
前回ご紹介した改良版であればboxplotステートメントが使えるので、前回のコードを微修正するだけで作図可能になります。
datalatticeレイアウトを使用する場合（非推奨）
まず最初に思いつくのは、箱ひげ図のステートメントを追加することです。
前回のテンプレートにboxplotステートメントを追加すればOKかと思いきや…これでは動きません。
実はdatalatticeレイアウト内では内部で統計量を算出するステートメントは使用できません。
boxplotは内部で箱ひげ図に必要な統計量を自動で計算するため、今回のテンプレートでは使用できずエラーとなります。
ではあらかじめ統計量を算出してboxplotparmステートメントを利用して作図するのはどうでしょう？
残念ですがこれも作図不可です。
boxplotparmステートメント自体は利用可能ですが、データセットは１つしか読み込めないため、統計量のデータセットと推定した密度のデータセットをマージする必要があります。この時両者のレコード数は異なるため、統計量と統計量名を格納した変数は必ず欠損となるレコードが発生します。
boxplotparmステートメントのstat引数で指定する変数は、欠損値を許容しません。そのため上記の方法ではとエラーとなります。
なのでこのテンプレートでboxplotステートメントを使用することは不可能だと思われます。
そのため考えられる方法ははdatalatticeレイアウト上で、polygonplotで箱ひげ図の各パーツを描画する方法ですがとても煩雑です。
一応作図方法を掲載します。
作図法
作図手順は以下の通りです。
- KDEを実行する。
- 箱ひげ図に必要な統計量を算出する。
- 箱ひげ図の図形を作図するのに必要な座標データを作成する。
- 2と3で作成したデータをマージし、作図データセットを作成する。
- bandplotで密度を作図し、vectorplot、polygonplot、scatterplotで箱ひげ図を作図する。
要は箱ひげ図の各要素の頂点座標を求めて各種ステートメントで１つずつ作図するわけですね。
非常にめんどくさいですが、各パーツの体裁を微調整できるし、マクロ化することも可能です。
ここでは前回使用したデータセットを用いて箱ひげ図とバイオリンプロットを作図します。
KDEを実行する
前回と同じ方法でKDEを実行し、密度を推定します。
箱ひげ図に必要な統計量を算出する
箱ひげ図に必要な統計量は、4分位数であればmeansプロシジャなどで算出可能なのですが、外れ値の抽出はできません
外れ値を検出するプロシジャはないらしいので本来はマニュアルで算出する必要がありますが、公式ヘルプページに箱ひげ図の作図データセットを作成できるマクロが公開されているのでこれを使用します。ていうかuniverateかmeansプロシジャに実装しなさいよ・・・
マクロの説明通りに実行するとデータセットboxdataに算出結果が出力されます。また要約統計量が格納されたデータセットsummaryも出力されます。
boxdataはそのままboxplotparmステートメントでの作図に使用できます。カテゴリ変数deathcauseは自動的に変数xに変換されます。
data chol;
set sashelp.heart;
where deathcause ne "";
keep deathcause cholesterol;
proc sort; by deathcause;
run;
*統計量の算出;
%boxcompute(indsn=chol,x=deathcause,y=cholesterol);作図に各四分位数の値が必要になりますので、boxdataとsummaryをマージしておきます。
data boxdata2;
merge boxdata summary;
by x;
run;
箱ひげ図の図形を作図するのに必要な座標データを作成する
統計量から箱ひげ図の図形描画に必要な座標を計算します。座標の詳細は以下の通りです。dは箱の幅を決定する変数です。体裁に応じて調節しましょう。箱の頂点の座標を出力します。idは図形識別変数で、polygonplotステートメントには必須です。noは密度データセットとマージするのに必要なキー変数です。
*箱の幅;
%let boxwidth = 0.002;
*箱の座標;
data box;
length id $20;
set boxdata2;
if stat="Q1" then do;
id="box";
no=1;
box_x=&boxwidth.;
box_y=value;output;
id="box";no=2;
box_x=-&boxwidth.;
box_y=value;output;
end;
if stat="Q3" then do;
id="box";
no=3;
box_x=-&boxwidth.;
box_y=value;output;
id="box";no=4;
box_x=&boxwidth.;
box_y=value;output;
end;
keep x no id box_x box_y;
run;平均値は散布図マーカーで描画します。このときダミーのx変数dum_xを用意します。
*平均値の座標;
data meanData;
set boxdata2;
where stat ="MEAN";
mean=value;
dum_x=0;
no=1;
keep x no mean dum_x;
run;ひげと中央値は線分として表現します。線分の始点と終点を求め、キー変数noを生成します。
*中央値とひげの線分;
data medianData;
set boxdata2;
if stat ="MEDIAN" then do;
line_x1=-&boxwidth.;
line_y1=value;
line_x2=&boxwidth.;
line_y2=value;
no=1; output;
end;
if stat="MAX" then do;
line_x1=0;
line_y1=Q3;
line_x2=0;
line_y2=value;
no=2; output;
end;
if stat="MIN" then do;
line_x1=0;
line_y1=Q1;
line_x2=0;
line_y2=value;
no=3; output;
end;
keep no line_x1 line_y1 line_x2 line_y2 x;
proc sort; by x no;
run;外れ値も平均値のデータセットと同様です。キー変数noを作成し連番を振っておきます。
OUTLIERとFAROUTLIERでマーカーを変更できるようにするため、外れ値の区分を格納した変数も用意します。
*外れ値;
data outliner;
set boxdata2;
where stat in ("FAROUTLIER","OUTLIER") and value^=.;
retain no;
by x;
if first.x then no=1;
else no+1;
outliner=value;
*作図用ダミー変数;
dum_x=0;
outlinertype=stat;
keep x no stat dum_x outlinertype outliner;
run;作図データセットを作成する
作成した箱ひげ図用データセットと実行結果データセットをマージします。
マージする前に密度データが観察データの範囲内になるように抽出した後、キー変数noを作成し、x毎に連番を振ります。
proc sort data=density;
by deathcause value;run;
*観察データの最大値最小値;
proc means data=chol nway;
var cholesterol;
class deathcause;
output out=st(rename=(deathcause=x))
min=min max=max;
run;
data density2;
merge density(rename=(deathcause=x)) st;
by x;
density=round(density,1E-8);
mirror=-density;
if min<=value<=max; run;
*連番付与;
data density3;
set density2;
retain no;
by x;
if first.x then no=1;
else no+1;
run;
*箱ひげ図用データとマージ;
data plot;
merge density3 box meandata mediandata outliner;
by x no;
run;作図
proc template;
define statgraph violine2;
begingraph;
discreteattrmap name="out";
value "FAROUTLIER" / markerattrs=(symbol=diamondfilled size=10 color=cxE1812C);
value "OUTLIER" / markerattrs=(symbol=circlefilled size=10 color=cxE1812C);
enddiscreteattrmap;
discreteattrvar attrmap="out" var=outlinertype attrvar=_out;
layout datalattice columnvar=x / headerlabeldisplay=value headerborder=false border=false columnheaders=bottom
rowaxisopts=(linearopts=(tickvaluesequence=(start=0 end=600 increment=100)) label="Cholesterol")
columnaxisopts=(display=(line) label="Deathcause");
layout prototype / walldisplay=none;
*バイオリンプロット;
bandplot y=value limitupper=density limitlower=mirror/
display=(fill outline)
fillattrs=(color=cx3375A3)
outlineattrs=(color=black);
*箱;
polygonplot id=id x=box_x y=box_y / display=(fill outline)
fillattrs=(color=white) outlineattrs=(thickness=1 color=black);
*中央値、ひげ;
vectorplot x=line_x2 y=line_y2 xorigin=line_x1 yorigin=line_y1 /arrowheads=false lineattrs=(thickness=1 color=black) ;
*外れ値;
scatterplot x=dum_x y=outliner / group=_out
;
*平均値;
scatterplot x=dum_x y=mean /markerattrs=(size=15 symbol=circlefilled color=red);
endlayout;
endlayout;
endgraph;
end;
run;
ods graphics / width=30cm height=15cm noborder;
proc sgrender data=plot template=violine2;
run;一応できましたが、冒頭でも書いた通りマジでお勧めしません。計算がめんどくさすぎます。
改良版
前回紹介した改良版のほうがずっと簡単です。作図用データセットに箱ひげ図作図用変数をマージし、boxplotステートメントを追加するだけです。今回はマクロ化して使いまわせるようにしました。
%macro violin(ds= , x=, y=, xlabel=, ylabel=, cat_param=2.5);
*ds: プロットに使用するデータセット
x:カテゴリ変数(文字列化フォーマットを当てる)
y:応答変数（KDEで密度を推定したい変数)
xlabel:X軸ラベル
ylabel:y軸ラベル;
data dat;
set &ds.;
where &x. ne "";
x=&x.;
y=&y.;
keep x y;
proc sort; by x;
run;
*最大値と最小値を求める;
proc means data=dat nway;
var y;
class x;
output out=stat min=min max=max;
run;
*KDE;
proc kde data=dat;
univar y / out=density;
by x;
run;
*KDE結果から最大密度を取得し、プロット間隔を算出;
proc sql noprint;
select max(density)*&cat_param. into : cat_width
from density;
quit;
proc sort data=dat(keep=x) nodupkey out=cat; by x; run;
*ダミー軸変数、軸用フォーマットの作成;
data dummy;
length tickvalue $200;
set cat end=eof;
retain tickvalue ;
dum_x=_N_*round(&cat_width.,0.001);
if _N_=1 then tickvalue =strip(put(dum_x,best.));
*表示する目盛りを指定するため、リストを生成;
else tickvalue = catx(" ",tickvalue, strip(put(dum_x,best.)));
start=dum_x;
end=dum_x;
if eof then call symputx("valuelist", tickvalue);
fmtname="cat";
label=x;
run;
*目盛り用フォーマット作成;
proc format cntlin=dummy;
run;
*作図用変数の作成,ダミー軸変数の値だけプロットを右へ平行移動させる;
data density2;
merge density dummy(keep=x dum_x) stat;
by x;
label value="Cholesterol";
dens=density+dum_x;
dens2=-density+dum_x;
format dum_x cat.;
if min<=value<=max;
run;
*箱ひげ図用データセットを作成;
data dat2;
merge dat dummy(keep=x dum_x);
by x;
rename dum_x=dum_x2;
run;
data plot ;
merge density2 dat2;
run;
proc template;
define statgraph violine;
begingraph;
layout overlay /
xaxisopts=(label="&xlabel." linearopts=(tickvalueformat=cat. tickvaluelist=(&valuelist.))
offsetmin=0.05 offsetmax=0.05)
yaxisopts=(label="&ylabel.");
bandplot y=value limitupper=dens limitlower=dens2 /
group=dum_x
display=(fill outline)
fillattrs=(color=cx3375A3)
outlineattrs=(color=black pattern=solid);
boxplot x=dum_x2 y=y / capshape=none boxwidth=0.1 spread=true;
endlayout;
endgraph;
end;
run;
ods graphics / width=30cm height=15cm;
ods listing gpath="グラフ出力パス";
proc sgrender data=plot template=violine;
run;
%mend;とりあえず箱ひげ図の体裁はデフォルトのものを使用します。boxplotステートメントのオプションで
体裁を変更できます。
試しにsashelpのテストデータで作図してみましょう。
%violin(ds=sashelp.heart, x=deathcause, y=cholesterol , xlabel=deathcause, ylabel=cholesterol);
%violin(ds=sashelp.cars, x=type, y=mpg_city , xlabel=type, ylabel=mpg);
カテゴリの並び順がアルファベット順になってしまうので、任意の順番にしたい場合はソート変数を指定するコードの追加が必要ですね。またKDEの設定を変えたい場合もあるかもしれないのでマクロの引数を増やして設定を変更できるようにするのもありかと思います。
箱ひげ図は複雑な作図なのでステートメントが使えたほうがはるかに楽ちんです。オプションで体裁を制御できますし汎用性は最初に紹介したコードよりも上だと思います。
