Raincloud plot「クロスオーバー試験」用バージョン

 あんまり医薬に特化しすぎた話は好まないのですが，少しだけ．

臨床試験というものには色んなデザインがあるのですが，クロスオーバー試験っていうのがあって，1人の参加者がAの薬を服薬前に測定，服薬後に測定．その薬の効果が抜けるであろう十分な時間がたった後（ウォッシュアウト期間とかいう），Bの薬を服薬前に測定，服薬後に測定　みたいなイメージで

1人の参加協力者から，複数薬剤に対しての投与前後データが取得できるみたいな感じです
服薬順はA→B，B→Aにランダムで決めたりします．なにが嬉しいのって言うと，薬の効果を個人ごとに比較できるので，効果への個人差をキャンセルしうるってノリです

まあ，それは本題ではないのでおいておいて

1人の人からDrug-AAのper-postの値，Placeboのper-postの値をとりだして
薬剤ごとの推移を↑こんな感じでプロットしたりします

例えば，血中の何かを下げる薬だったとして，Drug-AAの方が下がり傾向が強いようになんとなく見えます（実際そういう風にデータ作ってます）


これをRaincloud　plotの応用で描いてみるとどうでしょうって話です

個別の推移が見れると共に，その背後の集団の分布そのものが投与前後でなんか動いてるのが視覚的にわかりやすい気がします

まぁ，ただ，実際のデータではここまでクッキリでないことも多いので，プロットを重ねてしまうことで，かえって見にくいとかもあるかもしれません

まあその辺は，試行錯誤ということで．

Raincloud　plotは，カーネル密度推定で得た確率密度関数の推定を雲にして（バイオリンプロット），それと四分位点と平均を表現した箱髭図を傘にして，データそのもののプロットを雨にする．データ間に対応関係がある場合は線でつなぐこともある（線は雷）ってコンセプトがおおまかに決まってるだけで，アレンジというか，状況やデータに応じて，カスタマイズするところが面白いなぁって感じてます

%let outpath=XXXXXX;

/*Test dataset*/

data test;

call streaminit(1080);

    do id=1 to 25;

     do Timepoint=1 to 2;

       do  TREAT="Drug-AA","Placebo";

          if TREAT="Drug-AA" then do;

               select(timepoint);

                    when(1) val=rand("normal",30,5);

                    when(2) val=rand("normal",25,4);

                end;

            end;

           if TREAT="Placebo" then do;

             select(timepoint);

                when(1) val=rand("normal",30,5);

                when(2) val=rand("normal",31,3);

             end; 

            end;

            output;

        end;

      end;

    end;

run;

/* calculation----kernel density estimation*/

proc sort data = test ;

 by TimePoint TREAT;

run;

proc kde data = test ;

 univar val /out=kde ;

 by TimePoint TREAT;

run;

/*dataset-fix*/

data wk2;

set kde(in=ina)

     test(in=inb)

;

by Timepoint ;

retain timecount 0;

conv=choosen(timepoint,-1,1);

if ina then do;

 if timepoint=1 then density1= conv*density + conv*0.2;

 if timepoint=2 then density2= conv*density + conv*0.2;

end;

if inb then do;

  box_x= conv*0.17;

  series_x=conv*0.1;

end;

TREAT_ID=cats(TREAT,ID);

run;

/*Plot-define*/

ods graphics on /   width=700 height=700;

ods html path="&outpath" file="test.html";

proc template ;

  define statgraph RCP ;

     begingraph;

              layout overlay

                / yaxisopts  = (display=none offsetmin=0.03 offsetmax=0.03 )

                   xaxisopts  = ( display=(label tickvalues ) label ='TimePoint' linearopts=(tickvaluelist = ( -0.2 0.2 ) tickdisplaylist=('Pre' 'Post') ))

;

                   bandplot y=value  limitupper=-0.2  limitlower=density1

                        / group=TREAT  fillattrs=(transparency=0.4) tip=none  ;

                   bandplot y=value  limitupper=density2  limitlower=0.2

                        / group=TREAT  fillattrs=(transparency=0.4) tip=none;

                   boxplot y = val x =box_x 

                        /boxwidth = 1 group=TREAT groupdisplay=cluster clusterwidth=0.1 name="box";

                   seriesplot x=series_x y=val  

                        /display=all group=TREAT_ID

                         markerattrs=(symbol=circlefilled size=8 transparency=0.4) 

                         lineattrs=(thickness=1 pattern=solid ) 

                        linecolorgroup=TREAT markercolorgroup=TREAT;

                    discretelegend "box"/location=inside valign=top halign=right;

              endlayout;   

     endgraph;

  end;

run;

/*Plot-submit*/

proc sgrender data=wk2 template=RCP ;

run;

ods html close;

Raincloud plot「経時的繰り返し測定データ」用バージョン

多群で経時反復測定のデータをどのように視覚的に把握するかという点においても
Raincloud plotはよいと思われます．梅雨時なので季節にもあってます

前回のコードからの改良として
パネル分割をやめて，1枚にして，それにともなって
BOXPARAMをBOXにして，事前計算を省略．さらにグラフ間の位置取りの調整も
基本的にtemplate内のevalにつっこんで，データステップでは，各時点間のスペースの具合を調整する用のところのみ残した





 %let outpath=XXXX;

/*Test dataset*/

data test;

call streaminit(1080);

do Timepoint=1 to 5;

    GROUP="A";

    do id=1 to 400;

            select(timepoint);

                when(1) val=rand("normal",30,5);

                when(2) val=rand("normal",33,4);

                when(3) val=rand("normal",35,5);

                when(4) val=rand("normal",34,6);

                when(5) val=rand("normal",34,4);

            end;

            output;

    end;

    GROUP="B";

    do id=1 to 400;

            select(timepoint);

                when(1) val=rand("normal",18,5);

                when(2) val=rand("normal",20,4);

                when(3) val=rand("normal",22,5);

                when(4) val=rand("normal",24,5);

                when(5) val=rand("normal",21,3.5);

            end; 

        output;

    end;

    GROUP="C";

    do id=1 to 400;

            select(timepoint);

                when(1) val=rand("normal",40,3);

                when(2) val=rand("normal",39,4);

                when(3) val=rand("normal",42,5);

                when(4) val=rand("normal",41,5);

                when(5) val=rand("normal",41,3.5);

            end; 

        output;

    end;

end;

run;

 

/* calculation----kernel density estimation*/

proc kde data = test ;

 univar val /out=kde ;

 by TimePoint GROUP;

run;

/*dataset-fix*/

data wk2;

set kde(in=ina)

     test(in=inb)

;

by Timepoint ;

retain timecount 0;

if first.Timepoint then timecount+0.3;

run;

/*Plot-define*/

ods graphics on /imagemap=on tipmax=1000000  ANTIALIASMAX=1000000;

ods html path="&outpath" file="test.html";

proc template ;

  define statgraph RCP ;

     begingraph;

              layout overlay

                / yaxisopts  = (display=none offsetmin=0.03 offsetmax=0.03 ) ;

                   bandplot x=value limitupper=eval(density-timecount)  limitlower=eval(0-timecount) 

                        / group=GROUP  fillattrs=(transparency=0.4) tip=none;

                   boxplot y = val x = eval(-0.03-timecount)

                        / orient   = horizontal group=GROUP  groupdisplay=cluster clusterwidth=0.2 

                           fillattrs=(transparency=0.4) connect=mean display=(caps fill mean median outliers notches  connect)

                            boxwidth = 1;

                   scatterplot x=val y=eval(-0.12 + 0.05*cdf('NORMAL', rannor(1234)) + coalesce(0, val)-timecount)

                       / group=group jitter=auto jitteropts=(axis=Y width=1) markerattrs=(symbol=circle size=5 transparency=0.4)

                    rolename=(tip1=ID tip2=VAL tip3=GROUP) 

                    tip=(tip1 tip2 tip3)

                    tiplabel=(tip1="ID" tip2="Value" tip3="GROUP")

                    ; 

              endlayout;   

     endgraph;

  end;

run;

/*Plot-submit*/

proc sgrender data=wk2 template=RCP ;

run;

ods html close;

SASで錯視，SAS視シリーズ①

data wk1;

GRP=1;x=2;y=2;output;

GRP=1;x=5;y=2;output;

GRP=2;x=2;y=1;output;

GRP=2;x=5;y=1;output;

run;

proc sgplot data=wk1 noautolegend;

series x=x y=y/group=GRP lineattrs=(color=black thickness=3);

xaxis values=(0 to 7) display=none;

yaxis values=(0 to 3) display=none;

run; 

同じ長さにみえますが

矢羽根をたせば

data wk1;

GRP=1;x=2;y=2;output;

GRP=1;x=5;y=2;output;

GRP=1.1;x=1.5;y=2.3;output;

GRP=1.1;x=2;y=2;output;

GRP=1.2;x=1.5;y=1.7;output;

GRP=1.2;x=2;y=2;output;

GRP=1.3;x=5;y=2;output;

GRP=1.3;x=5.5;y=2.3;output;

GRP=1.4;x=5;y=2;output;

GRP=1.4;x=5.5;y=1.7;output;

GRP=2;x=2;y=1;output;

GRP=2;x=5;y=1;output;

GRP=2.1;x=2;y=1;output;

GRP=2.1;x=2.5;y=1.3;output;

GRP=2.2;x=2;y=1;output;

GRP=2.2;x=2.5;y=0.7;output;

GRP=2.3;x=4.5;y=1.3;output;

GRP=2.3;x=5;y=1;output;

GRP=2.4;x=4.5;y=0.7;output;

GRP=2.4;x=5;y=1;output;

run;

proc sgplot data=wk1 noautolegend;

series x=x y=y/group=GRP lineattrs=(color=black thickness=3);

xaxis values=(0 to 7) display=none;

yaxis values=(0 to 3) display=none;

run;

同じ長さと認識しにくくなります
これがミュラー・リヤー錯視です
https://en.wikipedia.org/wiki/M%C3%BCller-Lyer_illusion


もう一個

data wk2;

do x=1.2 to 10.2 by 2;

do y=1.5 to 10.5  by 2;

output;

end;

end;

do x=0.8 to 8.8 by 2;

do y=0.5 to 8.5 by 2;

output;

end;

end;

run;

proc sgplot data=wk2;

scatter x=x y=y /markerattrs=(symbol="squarefilled" size=43 color=black);

refline 0 1 2 3 4 5 6 7 8 9 10/lineattrs=(thickness=4 color=gray);

xaxis values=(0 to 10) offsetmin=0 offsetmax=0 display=none;

yaxis values=(0 to 10) offsetmin=0.009 offsetmax=0.009 display=none;

run;

■を除いた線は，まっすぐな平行線だと思いますか？
scatterのblackをwhiteにでもしてもらえばすぐに確認できますが

案の定，ただの平行線です

まっすぐな線が，まっすぐに見えなくなる．これ，結構錯視の中でも私が好きなやつで

カフェウォール錯視っていいます
https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%95%E3%82%A7%E3%82%A6%E3%82%A9%E3%83%BC%E3%83%AB%E9%8C%AF%E8%A6%96

Raincloud plot「序章」

OpenResearchの，こちらの論文を読んで凄く感銘をうけました

Raincloud plots: a multi-platform tool for robust data visualization 

https://wellcomeopenresearch.org/articles/4-63

Raincloud plotへの並々ならぬ愛もなんですが，データビジュアライゼーションへの真摯な知見もですし，なによりWebコミュニティの集合知を最大限に活用しており，非常にいいなと．

だけど，まあ，しょうがないんですけど，SASで描いたバージョンはないのです(笑)

ということで，今年の某総会の発表ネタの1本はRainCloud PlotにSASバージョンの彩りを加えるにしようかなとか思ってます

個別の値にマウスオンしたら，どのサンプルで，どんな値かもバッチリ表示

まだまだここから磨いていかねばと思ってますがとりあえず．
datalatticeとprototypeでカテゴリ分割してますが，次は時系列で，代表値推移を線で結んで雷にみたてたバージョンもやりたいので，それだと今のやり方だとまずいですからね
また考えようかなと思ってます

%let outpath=XXXXX;

/*Test dataset*/

data test;

call streaminit(1080);

GROUP="A";

do id=1 to 200;

        val=rand("normal",10,5);

        output;

end;

do id=201 to 400;

        val=rand("normal",30,5);

        output;

end;

GROUP="B";

do id=1 to 200;

        val=rand("normal",18,5);

        output;

end;

do id=201 to 400;

        val=rand("normal",22,5);

        output;

end;

run;

 

/* calculation----kernel density estimation*/

proc kde data = test ;

 univar val /out=kde ;

 by GROUP;

run;

/* calculation----for box-and-whisker*/

ods _all_ close;

ods output SGPlot=box;

proc sgplot data=test;

 vbox val/group=group;

run;

ods html;

/*dataset-fix*/

data wk2;

set kde(in=ina)

     box(in=inb where=(^missing(BOX_VAL_GROUP_GROUP____Y)))

     test(in=inc)

;

call streaminit(0615);

if ina then low=0;

if missing(BOX_VAL_GROUP_GROUP___ST) then BOX_VAL_GROUP_GROUP___ST="DUMMY";

if ^missing(BOX_VAL_GROUP_GROUP___GP) then group=BOX_VAL_GROUP_GROUP___GP;

if inb then dummy_x=-0.01;

if inc then do;

    dummy_y=-0.05;

    random=rand("uniform")*0.01;

    if ranuni(777)<0.5 then dummy_y=dummy_y - random;

    else dummy_y=dummy_y + random;

end;

run;

/*Plot-define*/

ods graphics on /imagemap=on tipmax=5000;

ods html path="&outpath" file="test.html";

proc template ;

  define statgraph RCP ;

     begingraph;

          entrytitle "Rain Cloud Plot";

           layout datalattice rowvar=group/  columnaxisopts=(label="Value") rowaxisopts=(display=none);

              layout prototype;

                   bandplot x=value limitupper=density  limitlower=low / display=all ;

                   boxplotparm y=BOX_VAL_GROUP_GROUP____Y x=dummy_x   stat=BOX_VAL_GROUP_GROUP___ST /boxwidth=0.3 orient = horizontal  ;

                   scatterplot x=val y=dummy_y/ jitter=auto jitteropts=(axis=Y width=1) markerattrs=(symbol=circle size=8 transparency=0.4)

                    rolename=(tip1=ID tip2=VAL) 

                    tip=(tip1 tip2)

                    tiplabel=(tip1="ID" tip2="Value")

                    ; 

              endlayout;   

           endlayout;

     endgraph;

  end;

run;

/*Plot-submit*/

proc sgrender data=wk2 template=RCP ;

run;

ods html close;