VARCHARデータタイプ(可変長文字型)

Viyaとかでサポートしてる新しいCASエンジンだと,いろんなデータ型をサポートしてますが,ご存じ,SAS V9エンジンは,文字型と数値型の2本槍の漢仕様ですね

ところが,そんなV9さんも実はデータステップのセッション中に限りVARCHAR,可変長文字列をサポートすることができます
ステップが終了するとともに,従来の文字型(Char Type)に強制変換されますが,その辺ちょっと紹介

data A;
   length x y z varchar(*);
   x='A'; y='B'; z='あ';
run;

とすると可変長変数が作成できます




ログがなんかいうてますが,気にせず contentsにかけてみます

proc contents data=A;
run;










*を指定すると,文字型に戻った後32767のlengthになるようですね

そういや,ログでmsglevel=iつけろとか言うていたので,つけるとともに
アスタリスクではなくて10としてみると

options msglevel=i;
data B;
   length x y z varchar(10);
   x='A'; y='B'; z='あ';
run;








具体的に変換された変数を教えてくれました.INFOじゃなくてNOTEなんすね


proc contents data=B;
run;











長さは20.



で,だから,SAS9でほんのひととき,VARCHARにできたところで,なにが嬉しいの?
って話ですが

まず一つ目は,SAS忘備録の中の人に教えてもらった話

通常のSASの文字は固定長なので,よくある失敗ですが
||でつなぐと

data C;
   length X Y  $10.;
   x='A'; y='B'; ;
   z=x || y;
run;






余計な空白が入りがちってやつですが

これが可変長であれば

data C;
   length X Y varchar(10);
   x='A'; y='B'; ;
   z=x || y;
run;





このように,実際の長さに応じて調整が入るので,いちいちstripだなんだので空白除去が不要になる



あと,公式でちょっと面白いこと書いてて
https://documentation.sas.com/doc/ja/vdmmlcdc/8.1/nlsref/n01zuyhmhz32ufn1rajvb5xl4462.htm#p1gtkuek0ora5in1k39voups9j0g






へーー,どうやら,通常のSUBSTRやINDEXがVARCHAR型に対しては,KSUBSTR,KINDEXと同じ効用をもつようになるみたい.

早速実験



data A;
length x1 $10.  x2 VARCHAR(10);
x1="あいうえお";
x2="あいうえお";

y1=substr(x1,1,1);
ky1=ksubstr(x1,1,1);
y2=substr(x2,1,1);
ky2=ksubstr(x2,1,1);

run;






おー,ほんまや.通所の文字型に日本語入れてSUBSTRかけると文字化けするけど
VARCHARならKSUBSTRと同じように,文字単位で処理されてる

もしかすると,使い方によっては,活きてくるシチュエーションあるかもしれませんね

MathMLを勉強しているのだけども…

いつも、一応ブログに書くことは、そこそこ、自分で理解してて、そんなに大きな間違いはなかろうってことを書くようにしてるのですが…、今回はちょっと自信なし。むしろ教えてって内容です

MathML(Mathematical Markup Language)に興味をもっていて、EPUB形式の中に統計解析の数式を埋め込むことで、SASと統計解析の電子書籍のテキストが作れないかなぁとか思って、ちみちみ勉強してるのですが、 MathMLの実装がよくわからなくて大苦戦中…
辛いのが、不可解な現象が起きた時に、それがEPUBリーダー起因なのか、私のMathMLの理解不足起因なのか、実装してるSASの問題なのかよくわからないという点

↓Firefoxのepubリーダーアドイン

























↓Google Play books

















そもそも、iPadもKindleも持ってないのに、EPUBとかやってんじゃねえよって話ですが(笑)

ods epub file="XXXX\test.epub" title="MathML Test";

ods escapechar='^';

proc odstext contents="";


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <mn>1</mn>

  <mo>+</mo>

  <mn>2</mn>

  </math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <mn>3</mn>

  <mo>&#xD7;</mo>

  <mn>3</mn>

  <mo>&#xF7;</mo>

  <mn>3</mn>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <mo>|</mo>

  <mrow>

    <mo>-</mo>

    <mi>a</mi>

  </mrow>

  <mo>|</mo>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <mi>y</mi>

  <mo>=</mo>

  <msub>

    <mi>log</mi>

    <mi>a</mi>

  </msub>

  <mi>x</mi>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

 <mrow>

    <mo>(</mo>

    <mfrac>

      <mn>1</mn>

      <mi>x</mi>

    </mfrac>

    <mo>+</mo>

    <mn>1</mn>

    <mo>)</mo>

  </mrow>

  <mrow>

    <mo>(</mo>

    <mi>x</mi>

    <mo>-</mo>

    <mn>1</mn>

    <mo>)</mo>

  </mrow>

</math>}';



p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <munderover>

    <mi>&#x2211;</mi>

    <mrow>

      <mi>i</mi>

      <mo>=</mo>

      <mn>1</mn>

    </mrow>

    <mi>n</mi>

  </munderover>

  <msub>

    <mi>k</mi>

    <mi>i</mi>

  </msub>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <mo>[</mo>

  <mtable>

    <mtr>

      <mtd>

        <mn>2</mn>

        <mi>x</mi>

        <mo>-</mo>

        <mi>y</mi>

        <mo>=</mo>

        <mn>1</mn>

      </mtd>

    </mtr>

    <mtr>

      <mtd>

        <mi>x</mi>

        <mo>+</mo>

        <mn>2</mn>

        <mi>y</mi>

        <mo>=</mo>

        <mn>8</mn>

      </mtd>

    </mtr>

  </mtable>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <munder>

    <mi>lim</mi>

    <mrow>

      <mi>n</mi>

      <mo>&#x2192;</mo>

      <mi>∞</mi>

    </mrow>

  </munder>

  <mfrac>

    <mn>1</mn>

    <mi>n</mi>

  </mfrac>

  <mo>=</mo>

  <mn>0</mn>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

  <msqrt>

    <mn>2</mn>

  </msqrt>

</math>}';


p '^{mathml <math xmlns="http://www.w3.org/1998/Math/MathML">

<math>

    <mo>(</mo>

    <mtable>

      <mtr>

        <mtd>

          <mn>1</mn>

        </mtd>

        <mtd>

          <mn>2</mn>

        </mtd>

      </mtr>

      <mtr>

        <mtd>

          <mn>-3</mn>

        </mtd>

        <mtd>

          <mn>4</mn>

        </mtd>

      </mtr>

    </mtable>

    <mo>)</mo>

    <mo>(</mo>

    <mtable>

      <mtr>

        <mtd>

          <mn>2</mn>

        </mtd>

        <mtd>

          <mn>2</mn>

        </mtd>

      </mtr>

      <mtr>

        <mtd>

          <mn>0</mn>

        </mtd>

        <mtd>

          <mn>5</mn>

        </mtd>

      </mtr>

    </mtable>

    <mo>)</mo>

  <mo>=</mo>

    <mo>(</mo>

    <mtable>

      <mtr>

        <mtd>

          <mn>2</mn>

        </mtd>

        <mtd>

          <mn>12</mn>

        </mtd>

      </mtr>

      <mtr>

        <mtd>

          <mn>-6</mn>

        </mtd>

        <mtd>

          <mn>14</mn>

        </mtd>

      </mtr>

    </mtable>

    <mo>)</mo>

</math>

</math>}';


p"よくわからないこと";

   list;

      item "<mfenced>が効かない気がする…";

      item "エンティティの参照定義ってどうすんの…?";

      item "連立方程式で{つかいたいのだけど、使うとコメントとみなされてエラーになるのはどういう現象??";

  item "行列の見栄えが悪い…どうしたら…";

    end;


run;




ods epub close;


2変量データに対しての3D Raincloud plot[アニメーションバージョン]

前使ったダミーデータを雑にずらして二峰性にしたのでなんか、ラクダのコブか、戦艦みたいになっちゃいました



3Dグラフの、視点角度をループでずらして何枚も作図しながらanimationでまとめてGIFファイルにしてみました
なんか微妙に縮尺がブレて、ん?みたいになるので、ほんとはなんか調整がいるのかも
底面に配置したのはHPBINプロシジャでの単純なQuantile Binningの10ビンバージョンです
ようするに各変量で、ただ10%ごとに分位点だしてそれでメッシュしただけっす。
2変量に対しての3次元RCPがプロットとして有用かどうかは、ちょっと私自身、まだ実装例がないので、ちょっと保留で。散布図とヒートマップ、BAGPLOTなど2次元のプロットで複合評価するほうが正確に把握できる気もしている… 身も蓋もないけど(笑)

data wk1;

 seed = 15678;

 do i = 1 to 500;

  z1 = rannor(seed);

  z2 = rannor(seed);

  z3 = rannor(seed);

  x = 3*z1 + z2;

  y = 3*z1 + z3;

  output;

 end;

 do i = 1 to 500;

  z1 = rannor(seed)+3;

  z2 = rannor(seed)+3;

  z3 = rannor(seed)+3;

  x = 3*z1 + z2;

  y = 3*z1 + z3;

  output;

 end;

 drop seed;

run;

proc kde data=wk1 ; 

   bivar x y/ out=kde; 

run;

data wk2;

set wk1(in=wk) kde(in=kde);

if wk then do;

    density2=ranuni(777)*0.004;

end;

if kde then do;

    x=value1;

    y=value2;

    if density>10**-5 then density2=density+0.01;

    else delete;

end;

keep x y density2;

run;


ods output Mapping=map;

proc hpbin data=wk1  numbin=10 quantile   ;

   input x y;

 run;

data anno2;

set map;

where same ^missing(LB);

if  variable="x" then do;

   function='move'; xsys='1';ysys='2';zsys='2';color="gray";line=2;size=0.1; x=0; y=LB; z=0;output;

   function='draw';   x=100; y=LB; z=0;output;

end;

if  variable="y" then do;

   function='move'; xsys='2';ysys='1';zsys='2';color="gray";line=2;size=0.1; x=LB; y=0; z=0;output;

   function='draw';   x=LB; y=100; z=0;output;

end;

run;


%macro plot;

%do rotate=0  %to 360 %by 10;

goptions reset=all border cback=white htitle=12pt; 

proc g3d data=wk2;

 scatter x*y=density2 /shape="balloon" annotate=anno2  yticknum=2 xticknum=2 size=0.5 rotate=&rotate noneedle

;

run; 

%end;

%mend;

options  ANIMATION=START  ANIMDURATION=0.5   PRINTERPATH=GIF nodate;

ods printer file="xxxxx\test\test.gif" ;

%plot

ods printer close ;

options  ANIMATION=STOP ;











ZIP読み込みの際のZIP内ファイル名

論文読んで実行していただいた方からのご指摘で,ダイレクトにzipを読み込む際に
各レコードがzip内からのどのレコードから来たかの識別で,filenameオプションで取得できるとか書いてしまったのですが,それはzip名を取得する方法でした…
Direct Zip Readingのオプションの説明が正しくないです
それもよく使うのでごっちゃになってしまってました… 恥ずかしい

うぅ… 本当すみません.せめてスライドだけでも差し替えができないかは問い合わせてみますしかし,こういったご指摘は大変ありがたいです

とりあえず正しい方法だけ 

a.csvとb.csvがあって





それをzというzipファイルに格納した場合

FILENAME indata ZIP "xxxx\z.zip" member="*";

data sample;

length fname zipname memname csvname $200.;

memname = ' ';

INFILE indata dlm="," filename=fname memvar=memname end=done ;

do while(^done);

 input a b;

 zipname=fname;

 csvname=memname;

 output;

end;

run;



ちょっとmemvarオプションがわかりにくいのですが,ブランクにするとzip内の全部の
ファイルが読み込み対象となり,読み込み開始時点でファイルの名前が入ります
zip内の全ファイルの走査が終わるまでinputを繰り返すという作業のためendとの併用で
ループになります


あぁ…,後輩には散々,見返せとか検証しろとか言っておいて,自分の論文でこれとは酷い

2変量データに対しての3D Raincloud plot(未完成版)

絶対SASでやらない方がいい気がする
3Dグラフの機能については,もうサポートする気がないって宣言してるレベルで無理.
たぶん,ヒートマップとかもそうですけど,2次元に上手く落とし込めるでしょ?みたいな主張なのかなぁ
SG系での3次元プロットの機能が旧GPLOT系より悪くなってるんだからどうしようもない

ただ,無理ゲーには挑みたくなるのが,魔界塔士の時代からの人のサガなので…

SAS社のサイトに載ってた適当なテストデータから

data wk1;

 seed = 15678;

 do i = 1 to 1000;

  z1 = rannor(seed);

  z2 = rannor(seed);

  z3 = rannor(seed);

  x = 3*z1 + z2;

  y = 3*z1 + z3;

  output;

 end;

 drop seed;

run;

2変量に対しての確率密度もKDEで問題なく計算できます
あとの流れは通常のRain cloudと一緒なんですけど
proc G3Dが一種類のプロットしか描けずに,組み合わせができないので 雲の部分をSURFACE(面プロット)で
描くのを諦めて,雨と同じ散布図の集合で表現するしかない… そうすると空を雨が覆ってしまうので
truncateの有無にかかわらず,確率密度関数が極小のところは消してしまうとスッキリ(いいのか?しらん)

proc kde data=wk1 ; 

   bivar x y/ out=kde; 

run;

data wk2;

set wk1(in=wk) kde(in=kde);

if wk then do;

density2=0.1+ranuni(777)*0.01;

end;

if kde then do;
          x=value1;

         y=value2;

     if density>10**-5 then density2=density+0.13;

else delete;

end;

keep x y density2;

run;

proc g3d data=wk2;

 scatter x*y=density2 /shape="balloon" size=0.5 noneedle;

run;



















Z軸に意味ないから軸けしたいけどGPLOTのAXISステートメントがZ軸に効かんという謎制約….NOAXISは3軸とも消すか消さないかを選択できるという,野性的すぎる大胆な仕様…
あとは,これだと,1変量の時の箱髭にあたる,位置の指標がないので
2変量用箱ひげ図的な用途で箱と同じくTukeyさまが編み出したBAGPLOTとかを配置してみたいけど,マイナー過ぎてSASに実装されてない(JMPにはあるのに…)
手で実装するにはしんどいし,実装できたとして,このG3Dに組み込む術は多分ANNOTATEしかない…

ビニング①

デジタルバイオマーカーの話するときに,入れそびれたネタを少し.
ビニング(binning)という処理があります
例えば,1~100までのデータがある時,1-10,11-20,21-30...のようにカテゴライズして数えあげてヒストグラム書いたりしますね
これは連続量を離散化してるわけですが,ヒストグラムは元のデータの特徴をつかむ時に使われるように,データの持つ特徴を損なわないまま階層化して,見通しをよくしているともいえます.1-10のような区間・階層を「ビン」といいます
また1000個のデータを10個のビンにはめると,データ量は一気に圧縮されます
また,良くも悪くも,外れ値や区間内のバラつきをビンは飲み込んでくれます

こういった性質から機械学習の前処理などでも用いられる手法です

さてSASの話

data test;

   length id 8;

   call streaminit("1234");

   do id=1 to 500;

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

      output;

   end;

   do id=501 to 1000;

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

      output;

   end;

run;



こういったデータあったとします.
分布としては以下の感じ















SASにおけるビニングは porc HPBINで実施します

ods output mapping=bin;
proc hpbin data=test numbin=10  bucket   ;
   input val;
 run;

numbinでビンの数を指定します. 

バケットビニング(デフォルト)

デフォルトなのでbucketはつけなくても同じ


















<0.4096..のビンには11データ
0.4096~4.517…のビンには56データというように ビニングがおきます
ビニングの結果はデータセット化もできます

バケットビニングとは
最大値ー最小値をビン数でわって出した値を,最小値から足して区切っていく方法です
ビン化イメージとしては














こんな感じ,ビン幅は等間隔ですね

次はQUANTILEビニング,ビン数の分だけ分位点つくるタイプです

proc hpbin data=test numbin=10 quantile   ;
   input val;
 run;



















分位点なので,境界値ぴったりとかはおいておいて,基本度数と比率はビン間で等しくなりますが,代わりにビン幅が不均等になり
データが疎のところは広く,密のところはせまくなります














あとはウィンザライズドビニング

proc hpbin data=test output=out numbin=10  winsor winsorrate=0.05   ;
   input val;
 run;

これはウィンザライズド平均を使う方法で,ウィンザライズド平均というのはトリム平均っぽいやつで,トリムが両端をカットするのに対して,ウィンザライズド平均は,端に当たらない値で,極致をLOCF上書きしちゃうイメージ
0 ,51, 52, 53 ,100というデータなら 51,51,52,53,53にして平均とる感じ
winsoerateで両端をどこまでとるか(なん%分までか)を定義

































あとはPSEUDO-QUANTILEってのも利用できますが,QUANTILEの結果に似ます
今回は割愛.

さて,次は,ビン化の性能評価とかを考えていきましょう だいぶ先になるかもだけど

HPDMDBでデータを要約データベースを得る

もともとSAS  enterprise miner にproc DMDB というものがあって(今は廃盤?)
それがHP(High Performance)として復活して,通常のSAS STATに13.2ぐらいから導入されたって流れですかね

Data Mining DataBase (DMDB)を作るプロシジャということなんですが,まあ,とてもシンプルなプロシジャで

 proc hpdmdb data=sashelp.cars 

     classout=cout  varout=vout1 ;

     var _numeric_  ;

     class _character_;

 run;

こうすると
CLASSOUTで出力したデータセット中は,指定した文字変数を全部FREQしたような
要約データセットができてます(形として立て積みの使いやすい形)
















でVAROUOTで出したデータセットの中は















こんな感じ.

まあ,通常のプロシジャのCLASSステートメントとちょっと違うのと
byも効かないので,層別集計とかには使いにくい

どちらかというと,初見のデータとか,大きなデータとかをとりあえずつっこんで
全体を把握するための要約データセットを作るみたいな(それをData Mining Databaseと一般的にいうのかどうか知らないですが笑
).

あと,臨床系ならSDTMとかADaM作った後に,とりまFREQみたいなことしてざっとQCしたりすることもあるかもしれませんが,カテゴリ値集約リストつくるなら,FREQやSUMMARY使うよりこっちの方がいいかもですよ.maxlevelオプションもあるようだし

で,いちおう,この子はHPプロシジャなのでPERFORMANCEステートメントとかあって,そっちのチューニングができるようになってるので,環境によってはかなりの性能を引き出せる余地があるかもです.ただ,普通のSingleマシーンモードで大きめのデータに使っても,そこそこ早い印象です

デバイスデータとかの1次クリーニングとかにもそこそこ便利ですよ


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;