前回に引き続き、バイオリンプロットを作図方法を検証します。
バイオリンプロットは箱ひげ図とともに作図されることがあります。
RのggplotやpythonのSeabornではviolinplotメゾッドのオプションで追加表示させることができますが、GTLの場合は素直にboxplotステートメントを
追加することで同様な作図が実施可能です。
ただしSAS公式ブログに取り上げられているdatalatticeレイアウトを使う場合は、boxplotステートメントは併用不可のため、手動で統計量を算出して、箱ひげ図の各パーツを
polygonplotステートメントとvectorplotステートメントで描画する必要があります。非常にめんどくさいのお勧めできません。
前回ご紹介した改良版であればboxplotステートメントが使えるので、前回のコードを微修正するだけで作図可能になります。
datalatticeレイアウトは使用不可
まず最初に思いつくのは、箱ひげ図のステートメントを追加することです。
SAS公式ブログのコードににboxplotステートメント（sgplotの場合のvboxに相当)を追加すればOKかと思いきや…これでは動きません。
実はdatalatticeレイアウト内では内部で統計量を算出するステートメントは使用できません。
boxplotは内部で箱ひげ図に必要な統計量を自動で計算するため、今回のテンプレートでは使用できずエラーとなります。
ではあらかじめ統計量を算出してboxplotparmステートメントを利用して作図するのはどうでしょう？
残念ですがこれも作図不可です。
boxplotparmステートメント自体は利用可能ですが、データセットは１つしか読み込めないため、統計量のデータセットと推定した密度のデータセットをマージする必要があります。
この時両者のレコード数は異なるため、統計量と統計量名を格納した変数は必ず欠損となるレコードが発生します。
boxplotparmステートメントのstat引数で指定する変数は、欠損値を許容しません。そのため上記の方法ではとエラーとなります。
なのでSAS公式ブログで公開されているdatalatticeレイアウトを使用した方法ではboxplotステートメントを使用することは不可能だと思われます。
そのため考えられる方法はpolygonplotで箱ひげ図の各パーツを描画する方法ですがとても煩雑です。作図手順は以下のようになります。
- KDEを実行する。
- 箱ひげ図に必要な統計量を算出する。
- 箱ひげ図の図形を作図するのに必要な座標データを作成する。
- 2と3で作成したデータをマージし、作図データセットを作成する。
- bandplotで密度を作図し、vectorplot、polygonplot、scatterplotで箱ひげ図を作図する。
要は箱ひげ図の各要素の頂点座標を求めて各種ステートメントで１つずつ作図するわけですね。
出来なくはないですが、かなり複雑になるためお勧めできません。
作図コード
前回紹介した改良版のほうがずっと簡単です。作図用データセットに箱ひげ図作図用変数をマージし、boxplotステートメントを追加するだけです。今回はマクロ化して使いまわせるようにしました。
%macro violin(ds= , x=, y=, xlabel=, ylabel=, cat_param=2.5);
*ds: プロットに使用するデータセット
x:カテゴリ変数(フォーマットを当てる)
y:応答変数（KDEで密度を推定したい変数)
xlabel:X軸ラベル
ylabel:y軸ラベル
cat_param: 各カテゴリの間隔を定義する。数値を大きくするとプロットの間隔が広がる;
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 violin;とりあえず箱ひげ図の体裁はデフォルトのものを使用します。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の設定を変えたい場合もあるかもしれないのでマクロの引数を増やして設定を変更できるようにするのもありかと思います。
箱ひげ図は複雑な作図なのでステートメントが使えたほうがはるかに楽ちんです。オプションで体裁を制御できますし汎用性は最初に紹介したコードよりも上だと思います。
