前回に引き続き、バイオリンプロットの作図方法を解説します。
前回はカーネル密度推定を実行し、作図データセットを作成しました。
今回はGTLで実際にバイオリンプロットを作成します。
公式ブログで紹介された方法をためす
公式ブログでは、クラス変数ごとに作図を実行するsgpanelプロシジャを使用し、
カーネル密度関数のグラフをを中心線に対して左右対称に2つ作図することでバイオリンプロットを作図する方法を紹介しています。
このため密度の符号を反転させた新規変数を作成する必要があります。
密度関数の作図は信頼区間を作図などに使われるbandplotステートメントを使用します。
今回はsgplotの代わりにGTLを使って作図します。
kdeで推定した密度をそのままプロットすると観察データの範囲外にも密度が表示されます。普観察データの範囲内で密度を作図するので、
作図前にデータを抽出する必要があります。
今回は観察されたコレステロール濃度の最大値と最小値に範囲内でカーネル密度を作図します。やり方はいろいろありますが、sqlで最小値と最大値を算出し、この値を基に
作図データを抽出します。
バイオリンプロットの塗りと枠線の設定は、bandplotステートメントのfillattrsオプションとoutlineオプションで指定できます。
data chol;
set sashelp.heart;
where deathcause ne "";
keep deathcause cholesterol;
proc sort; by deathcause;
run;
*KDE;
proc kde data=chol;
univar cholesterol / out=density;
by deathcause;
run;
*最大値と最小値を求める;
proc means data=chol nway;
var cholesterol;
class deathcause;
output out=stat min=min max=max;
run;
*作図用変数作成;
data density;
merge density stat;
by deathcause;
mirror=-density;
if min<=value<=max;
run;
proc template;
define statgraph violine;
begingraph;
layout datalattice columnvar=deathcause / headerlabeldisplay=value headerborder=false border=false columnheaders=bottom
rowaxisopts=(linearopts=(tickvaluesequence=(start=0 end=600 increment=100)))
columnaxisopts=(display=(line) offsetmin=0.1 offsetmax=0.1);
layout prototype / walldisplay=none;
bandplot y=value limitlower=mirror limitupper=density /
display=(fill outline)
fillattrs=(color=cx5975a4)
outlineattrs=(color=cx555555 pattern=solid);
endlayout;
endlayout;
endgraph;
end;
run;
ods graphics / width=30cm height=15cm;
proc sgrender data=density template=violine;
label value="Cholesterol";
run;
うーん。この方法ですとprototypeレイアウトで作図する必要があるんですよね・・・prototypeレイアウトはdatalatticeレイアウトに設置できる専用のレイアウトなのですが、
使用できる作図ステートメントが限られる上に内部で制御構文や数値演算ができないデメリットがあるので、このような複雑なグラフではあまり使いたくありません。
特に箱ひげ図を作図するboxplotステートメントが使えないのは致命的ですね。
またこの方法だと各バイオリンプロットは独立したプロットエリアに作図しているため、平均値を線で結ぶことはできません。
別の方法で作図方法を検討することにします。
改良版
prototypeレイアウトの使用はできれば避けたいので、overlayレイアウトでの作図を試みます。
ここで問題となるのはX軸の変数です。離散軸の中に密度をプロットするための線形軸を入れ子にして配置することはできないため、
ダミーのX軸変数を用意し、カテゴリを等間隔に並ぶようにバイオリンプロットを平行移動させます。
どれくらい移動させるかですが、定数にすると密度分布に応じて定数を調整する手間が発生するので、最大密度の定数倍だけバイオリンプロットを右へ平行移動させることにしました。
今回は密度の最大値の２.5倍だけバイオリンプロットを平行移動させます。
X軸はこのままだとダミーの数値が入った目盛ができてしまうため、tickvalueformatオプションとtickvaluelistオプションを併用して、各カテゴリ名を目盛として表示させます。
tickvalueformatで指定するフォーマットは作図用データセットのグループ名とダミーX軸変数の値から自動生成するようにしています。
マクロ化すれば任意のデータのバイオリンプロットを作図できるようになるはずです。
*バイオリンプロット同士の間隔;
%let cat_param=2.5;
*作図データセット;
data chol;
set sashelp.heart;
where deathcause ne "";
x=deathcause;
keep x cholesterol;
proc sort; by x;
run;
*最大値と最小値を求める;
proc means data=chol nway;
var cholesterol;
class x;
output out=stat min=min max=max;
run;
*KDE;
proc kde data=chol;
univar cholesterol / out=density;
by x;
run;
*KDE結果から最大密度を取得し、プロット間隔を算出;
proc sql noprint;
select max(density)*&cat_param. into : cat_width
from density;
quit;
proc sort data=chol(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;
*テンプレート;
proc template;
define statgraph violine;
begingraph;
layout overlay /
xaxisopts=(label="Deathcause" linearopts=(tickvalueformat=cat. tickvaluelist=(&valuelist.))
offsetmin=0.05 offsetmax=0.05);
bandplot y=value limitupper=dens limitlower=dens2 /
group=dum_x
display=(fill outline)
fillattrs=(color=cx3375A3)
outlineattrs=(color=black pattern=solid);
endlayout;
endgraph;
end;
run;
ods graphics / width=30cm height=15cm;
ods listing gpath="画像出力パス";
proc sgrender data=density2 template=violine;
run;うまく行きました。バイオリンプロットの配置を手動で調整しているので操作が煩雑なのですが、併用できるステートメントが増えるので汎用性が向上します。
バイオリンプロットは他の言語では簡単に作図できるのでSASも早く機能追加してほしいものです。
