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;