SAS OnDemand for Academics / SAS 9.4 利用手順

非商用目的で 無料で使えるSAS OnDemand for Academicsについて，始め方がわからないという方は，あるいはインストール版のSAS9.4について，インストールの仕方にとまどうって方が多いということで，利用手順をスライドにしてアップしてくれた方がいるので紹介させてくださいませ

https://www.docswell.com/s/NorihiroSuzuki/KNLD4E-SAS_ODA_9.4_installguide#p1

SASユーザー総会 論文集

 1982年～2024年までのSASユーザー総会論文集について，公開場所が分散していたりしていましたので世話人会のほうで一つ所にアーカイブしました

SASユーザー総会論文集アーカイブ

これは高橋 行雄 先生が2014年に当時の世話人会で行った，電子化・公開の取り組みを引き継いだものになります

目録のページには，全演題と，SAS社で公開されている発表スライドへのリンクをつけていってますが，確認しながら手作業でつけていっているので，まだ2000～2016年ぐらいまでです

また目録データはエクセルファイルでダウンロード可能です．
OCR由来なことと，名寄せ処理はおこなっておりませんが，簡単な分析はできると思います

エクセル取り込んでもらって

proc sgplot data=sugi;

where author_cat="筆頭筆者";

 vbar year;

 xaxis values=(1982 to 2024) valueattrs=(size=7);

run;
とかで　演題数の推移など

さて，2025年もユーザー総会について演題募集が始まっておりますのでこちらもぜひ
https://sas-user2025.ywstat.jp/

REGプロシジャでCLASSステートメントがなかったり，交互作用の事前処理が必要な点はHPREGで改善されているよという話

 REGプロシジャにはCLASSステートメントがないため，カテゴリ変数をモデルに含める場合，事前に数値型でダミー変数を用意しないといけなかったり X*Yといった交互作用の記述をmodelステートメント内にかけないことから，例えば

data class1;

 set sashelp.class;

 sexn = choosen(whichc(sex,"男子","女子"),0,1);

 age_height=age*height;;

run;

proc reg data=class1;

   model weight = age height sexn age_height;

run; 

quit; 

などとして

こういう感じに解析したりしてるわけですが

HPREGだと

proc hpreg data=sashelp.class;

   class sex;

   model weight = age height sex age*height;

run; 

quit; 

基本的にはREGプロシジャで，できることを押さえた上で，機能拡張されているような印象ですが，REGにあってHPREGにない機能もあり，似たことのできるGLM,GLMSELECTとの機能比較などは公式リファレンスにあるのでそちらを参照ください．

「PROC HPREG Contrasted with Other SAS Procedures」
https://documentation.sas.com/doc/en/stathpug/15.2/stathpug_hpreg_overview02.htm

テキストをモールス信号に変換する関数

やっぱり モールス信号って業務で良く使いますよね．
プログラマ派遣先にいったら，全ての通信機器を取り上げられ，鉄格子窓のある部屋に閉じ込められて，SASの入ったPCだけ与えられて労働させられたみたいなあるあるの状況の際に

電気の消灯点灯タイミングで外部に助けを求めることができるかどうかが大事です

いちいち調べるの面倒なので，変換関数作ってみました

proc fcmp outlib=work.functions.common;

function text2morse (text $) $ ;

length key1  $1 text  _rvalue rvalue $2000 ;

_text=upcase(compress(text));

 declare dictionary  d1 ;

   d1["A"] = "・－";

   d1["B"] = "－・・・";

   d1["C"] = "－・－・";

   d1["D"] = "－・・";

   d1["E"] = "・";

   d1["F"] = "・・－・";

   d1["G"] = "－－・";

   d1["H"] = "・・・・";

   d1["I"] = "・・";

   d1["J"] = "・－－－";

   d1["K"] = "－・－";

   d1["L"] = "・－・・";

   d1["M"] = "－－";

   d1["N"] = "－・";

   d1["O"] = "－－－";

   d1["P"] = "・－－・";

   d1["Q"] = "－－・－";

   d1["R"] = "・－・";

   d1["S"] = "・・・";

   d1["T"] = "－";

   d1["U"] = "・・－";

   d1["V"] = "・・・－";

   d1["W"] = "・－－";

   d1["X"] = "－・・－";

   d1["Y"] = "－・－－";

   d1["Z"] = "－－・・";

 do i = 1 to length(_text); 

  _key1=char(_text,i);

  _rvalue= d1[_key1];

  rvalue=catx(" ",rvalue,_rvalue);

 end;

  put _text= rvalue=;

  return(rvalue);

 endsub;

run;

で，使ってみます

options cmplib=work.functions ;

data test;

length x $200.;

x = "SOS";output;

x = "SAS";output;

x = "VVV";output;

run;

data a;

set test;

y = text2morse(x);

run;

proc print data=a noobs;

run;

こっちのけんと さんのおかげで誰でもわかるようになったSOSとかちゃんと変換できてますね
VVVはムスカ大佐が通信してたやつで，本日は晴天なり　みたいな試験信号です

数字や記号は今回無視しています

Excelのガントチャートでイベントの持続期間を可視化

 

👆
こういうのを簡単に作る方法．DDEとかは使わない．SASのRWI(Report Writing Interface)とODS EXCELでやります

テストデータとして，

data wk1;

   SUBJID="AAA";AETERM="事象1"; AESTDY=1;AEENDY=5;output;

   SUBJID="AAA";AETERM="事象2"; AESTDY=3;AEENDY=8;output;

   SUBJID="BBB";AETERM="事象1"; AESTDY=2;AEENDY=30;output;

   SUBJID="BBB";AETERM="事象2"; AESTDY=-2;AEENDY=3;output;

   SUBJID="BBB";AETERM="事象3"; AESTDY=4;AEENDY=.;output;

   SUBJID="CCC";AETERM="事象1"; AESTDY=15;AEENDY=20;output;

run;

要は，開始と終了の相対日付（DY変数）をもっていればなんでもOK
終了日がnullの場合は，継続として最大までバーをひっぱる仕様を想定

以下　コードですが，面白いのがexcelに出力するSAS RWIの部分での指定が，信じられないぐらいシンプルなのと，セル幅が自動でlengthに連動してくれるので細かいメッシュにたいして何の範囲指定や命令を書かなくても全部えー感じにしてくれるところ．
この技術はわりと使える

proc sql noprint;

 select min(min(AESTDY,1)) into:min trimmed from wk1;

 select max(max(AESTDY,AEENDY)) into:max trimmed from wk1;

quit;

%put &=min;

%put &=max;

%let range =%eval( %sysfunc(range(&min,&max))+ (0<&min));

%put &=range;

data wk2;

 set wk1;

 array ar{&range.} $1.;

 if &min < 1 then do;

    st = AESTDY - &min + (AESTDY<1);

    if ^missing(AEENDY) then en = AEENDY - &min + (AEENDY<1);

 end;

 else do; 

    st=AESTDY;

    en=AEENDY;

 end;

 do i = 1 to &range;

    if st <=i <=en then ar{i}="Y";

    else if st=i and missing(en) then ar{i}="Y";  

    else if st<=i and missing(en) then ar{i}="O";  

 end;

run;

%macro hed();

 %do i = &min %to &max;

    if &i ne 0 then ob.format_cell(data: "&i",  style_attr: "background=blue color=white");

 %end;

%mend; 

%macro row();

 %do i = 1 %to %eval(&range);

   if ar&i="Y" then  ob.format_cell(data: ar&i,  style_attr: "background=red color=red");

   else if ar&i="O" then  ob.format_cell(data: ar&i,  style_attr: "background=pink color=pink");

   else ob.format_cell(data: ar&i,  style_attr: "background=white color=white");

 %end;

%mend; 

ods excel options( sheet_name= "GANT" );

data _NULL_;

   set wk2 end=EOF;

   if _N_=1 then do;

       dcl odsout ob();

       ob.table_start();

    *** header ;

       ob.head_start();

       ob.row_start();

            ob.format_cell(data: "ID",  style_attr: "background=darkblue color=white");

            ob.format_cell(data: "Event",  style_attr: "background=darkblue color=white");

            ob.format_cell(data: "STDY",  style_attr: "background=darkblue color=white");

            ob.format_cell(data: "ENDY",  style_attr: "background=darkblue color=white");

            %hed();

       ob.row_end();

       ob.head_end();

   end;

   *** report ;

   ob.row_start();

        ob.format_cell(data: SUBJID);

        ob.format_cell(data: AETERM);

        ob.format_cell(data: AESTDY);

        ob.format_cell(data: AEENDY);

         %row();

   ob.row_end();

   if EOF then ob.table_end();

run;

ods excel close;

Proc ODSTABLEでRTFを作る際に任意の位置で改ページをする方法

Proc OdstableはProc Reportでできる機能をほぼすべて実装可能なうえに，構文がよりシンプルであり,個人的には完全な上位互換だと思っています

ただ，RTFを作成する際の大きな弱点として，Proc ReportのbreakがOdstableにはなく任意の位置での改ページができないことがありました．
2018年にSASユーザー総会で「PROC ODSTABLEを用いた帳票作成」という論文が発表されていますが，そちらでもその点書かれていました．

8年間色々，チャレンジしてきましたが，最近また向き合う機会があり，一応実装できたので共有．

data output;

 set sashelp.class;

 if mod(_N_,3)=0 then break ="Y";

 keep Name Age break;

 run;

↑
break="Y"の直前のレコードまで出力して，Yのレコードから次ページにうつるようにしてみます

ods rtf ;

proc odstable data=output  ;

   column Name Age break ;

   cellstyle break="Y" as data{PRETEXT="(*ESC*)R'\pagebb' " };

   define Name; end;

   define Age;  end;

   define break; print=no; end;

run;

ods rtf close;

なんとたったのこれだけで

みたいに，改ページできていることがわかります

結局種がわかってしまえば簡単で，インラインフォーマット関数に改ページがあることと
print=noでも，RTFに埋め込みさえすれば作動するのでcellstylesのPRETEXTで入れてしまえばいいということでした

開始終了日の情報から持続期間をまとめる話

まあ，医薬なら有害事象の持続とか，薬剤の継続使用を特定のカテゴリでまとめて
切れ間なく続いていた期間をまとめたレコードを合成するってことはあるのですが

これが結構SASだとやりにくい．SQLでもやれるんだけど，レコード一時的にかさむから
環境のパワーないとちょっと辛い時がある

なので結局，ハッシュオブジェクトを使ってしまう

👆適当なサンプルですが，症例ごとにカテゴリで，切れ間なく持続してる期間でまとめて
レコードをおこす

data wk1;

length ID cat term $30.;

format st en yymmdd10.;

ID="001";cat="痛み";term="頭痛";st="01JAN2025"d;en="05JAN2025"d;output;

ID="001";cat="痛み";term="筋肉痛";st="06JAN2025"d;en="10JAN2025"d;output;

ID="001";cat="痛み";term="関節痛";st="12JAN2025"d;en="13JAN2025"d;output;

ID="001";cat="痛み";term="神経痛";st="14JAN2025"d;en="19JAN2025"d;output;

ID="001";cat="胃腸障害";term="胃痛";st="01JAN2025"d;en="08JAN2025"d;output;

ID="001";cat="胃腸障害";term="腸炎";st="06JAN2025"d;en="12JAN2025"d;output;

ID="002";cat="感染症";term="インフルエンザ";st="02JAN2025"d;en="08JAN2025"d;output;

ID="002";cat="感染症";term="感冒";st="01JAN2025"d;en="10JAN2025"d;output;

run;

ここで，setしつつ，setしたもの自体をハッシュオブジェクトにマルチキー許容でいれて
自症例の全検索をかけちゃいます．ハッシュオブジェクトはメモリ内走査なので，特にアルゴリズム工夫なく力技で十分いけます

proc sort data=wk1;

by ID cat st en;

run;

data wk2;

length ID cat $30. cat_st cat_en 8.;

format cat_st cat_en yymmdd10.;

set wk1;

if 0 then set wk1(rename=(st=_st en=_en));

if _N_ = 1 then do;

  declare hash h1(dataset:"wk1(rename=(st=_st en=_en))",multidata:"Y");

  h1.definekey("ID","cat");

  h1.definedata("_st","_en");

  h1.definedone();

end;

cat_st = st;

cat_en = en;

do while(h1.do_over() = 0);

if n(_st,_en) = 2 then do;

  if _st-1<=cat_en<=_en and cat_en <= _en then do;

cat_en = _en;

if _st < cat_st then cat_st = _st; 

  end;

end;

end;

keep ID cat cat_st cat_en;

run;

do_overメソッドはキーがヒットする限り0を返し続けハッシュオブジェクトの中を探しつづけます．なのでdo whileと相性がいいですね．
そのままだと，元レコード数と変わらず集約されてないので，一番包括的なものだけを残します

 

proc sort data=wk2 ;

by ID cat  cat_en cat_st;

run;

 

data wk3;

set wk2;

by ID cat  cat_en cat_st;

if first.cat_en;

run;

どうでしょう？そんな　難しくないと感じるのではないでしょうか

割合の正確な95%信頼区間(Clopper-Pearsonの信頼区間 二項分布に基づく割合の信頼区間)

 多分，SASのFREQに癖があるせいだと思うのですが，割合の95%信頼区間のSASでの処理がどうも苦手という相談を受けましたので

少しまとめました． 

まず，基本的には　1症例1レコードで，対象となるイベントが起きたか起きてないかがY/Nや1/0で変数でもっているケースを考えます

オブザベーション数＝解析対象集団の例数＝分母　の状況です

data wk1;
SUBJID="001";FL1="Y";FL2="N";FL3="Y";;FL4="N";output;
SUBJID="002";FL1="N";FL2="N";FL3="Y";;FL4="N";output;
SUBJID="003";FL1="Y";FL2="Y";FL3="Y";;FL4="N";output;
SUBJID="004";FL1="N";FL2="N";FL3="Y";;FL4="N";output;
SUBJID="005";FL1="Y";FL2="N";FL3="Y";;FL4="N";output;
run;

 

オススメなのは，まずCLASSDATA集計してしまうことです．なぜかというと，YN両方の水準がそろっていないとFREQは計算ができないからです

なので

data clds;
if 0 then set wk1(keep=FL1 FL2 FL3 FL4);
do val ="Y","N";
 FL1=val;
 FL2=val;
 FL3=val;
 FL4=val;
 output;
end;
stop;
drop val;
run;

そしてこれを使って集計します

proc summary data=wk1 classdata=clds;

class FL1 FL2 FL3 FL4;
ways 1;
output out=out1;
run;

 

👆

各変数　N,Yとそのカウント数が_FREQ_に入ります．summaryでways 1に絞ったのは，1水準ごとの結果しか必要ないからです

 そして　FREQプロシジャで

ods output BinomialCLs=out2;

proc freq data=out1;
table FL1 FL2 FL3 FL4/binomial(exact level="Y");
weight _FREQ_/zeros;
run;

とします

ここで　binominalの中のexactはClopperpeasonと書いてもOK．大事なのはlevel=で見たい水準を必ず指定すること．そうじゃないと今回のNのように逆側の割合をとってしまうことがあるので．

そしてzerosオプションをつけないと0の場合の計算ができずにERRORになるので注意．

 ods output BinomialCLs=out2;

により結果は

のようになります．
ods outputをつかわずに

ods output BinomialCLs=out2;
proc freq data=out1;
table FL1 FL2 FL3 FL4/binomial(exact level="Y");
weight _FREQ_/zeros;
output out=a binomial;
run;

 とした場合，

👆

でるっちゃでるんですが，tableステートメントの複数指定には対応していないので，最後のもの(FL4)しかでていません．

なので，対象変数が1つの場合のみですね．そうでなければ　ods outputがいいと思ってます

 あとは適当に成形して

data result1;
length out1 - out3 $200.;
set out2;
out1 = strip(scan(Table,2,":"));
out2 = put( round(proportion * 100,0.1), 8.1 -L);
out3_1 = put( round(LowerCL * 100,0.1),8.1 -L); 
out3_2 = put( round(UpperCL * 100,0.1), 8.1 -L); 
out3 = cats("[",catx(", ",out3_1,out3_2),"]");
keep out1 out2 out3;
run;

 

みたいな．

 

で，次はADAEみたいに，全症例のYNではなくて，発生したイベントのみが積みあがってるデータについて考えてみます

data wk2;
length SUBJID AETERM $20. ;
SUBJID="001";AETERM="頭痛";output;
SUBJID="001";AETERM="腹痛";output;
SUBJID="002";AETERM="倦怠感";output;
SUBJID="002";AETERM="腹痛";output;
SUBJID="002";AETERM="頭痛";output;
SUBJID="003";AETERM="倦怠感";output;
SUBJID="003";AETERM="咳嗽";output;
SUBJID="003";AETERM="頭痛";output;
SUBJID="003";AETERM="腹痛";output;
SUBJID="004";AETERM="頭痛";output;
SUBJID="005";AETERM="頭痛";output;
run;

 

この場合，解析対象集団の全症例にレコードが起きてるわけではないので，別途，解析対象集団の例数(分母になるN)

を求めておき，マクロ変数等に入れておきます 

/*解析対象集団のNが5の場合*/
%let Big_N=5;
 

で，集計したい軸で症例内の重複を除去して

proc sort data=wk2 out=_wk2 nodupkey;
by SUBJID AETERM;
run;

 

それを集計
proc summary data=_wk2 nway;
class AETERM;
output out=out3;
run;

あとは，データステップで計算してあげるだけです．

基本，Clopper-Pearsonの信頼区間の計算式はF分布使うのですが，F分布から変形して求められるベータ分布使うと計算式がグッと短くなるので

そっち使うのが一般的と思います

data result2;
set out3;
if 0 < &Big_N then per = round(divide(_FREQ_ , &Big_N) * 100,0.1);
if per=0 then CI_LOW=0;
if per=100 then CI_HIGH=100;
if per ne 0 then CI_LOW=round((1-betainv(.975,(&Big_N - _FREQ_+1),_FREQ_))*100, 0.1);
if per ne 100 then CI_HIGH=round((1-betainv(.025,(&Big_N - _FREQ_),_FREQ_+1))*100,0.1);
drop _TYPE_;
run; 

F分布で計算すると

data result3;
set out3;
if 0 < &Big_N then per = round(divide(_FREQ_ , &Big_N) * 100,0.1);
if per=0 then CI_LOW=0;
if per=100 then CI_HIGH=100;
if per ne 0 then CI_LOW = round( (1+ (&Big_n - _FREQ_+1) / (_FREQ_ *finv(0.025, 2*_FREQ_, 2*(&Big_n-_FREQ_+1))))**-1*100,0.1);
if per ne 100 then CI_HIGH = round((1 + (&Big_n - _FREQ_) / (( _FREQ_ +1)*finv(0.975, 2*(_FREQ_ + 1), 2*(&Big_n - _FREQ_))))**-1*100,0.1);
drop _TYPE_;
run;

 

同じですね．

FCMP.SOLVEで方程式の実数解をガウス・ニュートン法で求める

 まずSASの公式に乗ってる例から

y = 1 / √x

でy=20となるxの値を求めてみます

グラフにすると

なので，まあ　だいたい値は0.0023～26あたりかなぁって感じですが

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

     y = 1/sqrt(x);

      return(y);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x);

   outargs x  ;

   x=solve("f1",{.}, y , .);

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(x);

 call solver(20,x);

run;

proc print data=wk1;

run;

👆
と，このように，あらかじめ解きたい式をFCMPで指定し，
FCMP.SOLVE関数の第一引数に　解きたい関数名，第三引数に目標値（ここでいうところの左辺，y）を指定します，
第二引数には計算のための5つのパラメータを配列で与えるのですが，{.}とするとデフォルト設定になります
最後の引数は，第一引数にとった関数に追加引数をここで与えることができます

つぎに3次方程式

を解いてみます

グラフにすると

👆どうも実数解は3つありそうです．ぱっとみ，交点は1,2,3らへん

この場合，初期値の設定によって，収束する解が異なるので，コールルーチンをちょっと改造して，パラメータを細かく指定できるようにして実行してみます

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

      return(x**3 -6*x**2 + 11*x -6);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

  /*設定パラメータ */

   array opts[5] initial_value absolute_criterion relative_criterion maximum_iterations solve_status;

   outargs x ,solve_status,y ;

   x=solve("f1"

               ,opts

               , 0

               , .);

   put x= solve_status=;

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(of x);

 absolute_criterion=1.0e-12;

 relative_criterion=1.0e-6;

 maximum_iterations=100;

 solve_status=.;

 y=0;

 initial_value=0.001;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=1.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=2.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

run;

proc print data=wk1;

run;

グラフから，おそらくｘ＝1,2,3あたりに解がありそうなので
初期値(initial_value)を，0.001，1.9, 2.9として回してみると

👆
やはり ｘ=1,2,3が解でした．solve_status=0は収束していることを示します

UpSet plot (アップセット プロット)の実装例

 UpSet-Plotというプロットをよくみかけます．アイテムの組み合わせ集計を，1個での集計，2個のアイテムセットでの集計，3個でのアイテムセットの集計...... n個のアイテムセットの集計　って，バスケット分析のような集計をしたときにそれをグラフ化します．
星とり表の部分がアイテムセットを示して，上のヒストグラムはそのカウント数です．
ここで，たとえば以下のグラフはある仮想の疾患に罹患した際に，患者に発現した症状のプロットです．Aさんが発熱と倦怠感の症状があったい場合，「発熱＆倦怠感」のアイテムセットとして集計されます

👆
ここで注意は
上図のヒストグラム(上側)の一番左は「咳嗽」で，●が一つですが，これは，「咳嗽」だけを症状として，もった人がカウントされます．「咳嗽」と「発熱」がある人は「咳嗽」にはカウントされずに「咳嗽＆発熱」でカウントされます．

では全体として，「咳嗽」を少なくとも1症状として持つ人数は？？となった場合は横に延びてるヒストグラムがそれにあたるという仕組みです

SASコードは以下です（冒頭のグラフ描いたやつです）．一応，データが変わっても，ID , item, itemnameをあわせてもらえれば可変で動くとは思いますが，多少調整いるかもです．

もうちょっと，洗練されてどこかで出そうと思いますが，第一ドラフトということで．
barchartparamで上側のヒストグラム，scatterとhighlowで星取表部分，横のヒストグラムはhighlowで描いてます．

/*テストデータ*/

data test;

call streaminit(123);

do ID = 1 to 100;

  iend=rand("integer",1,3);

  do i = 1 to iend;

    item = rand("integer",1,5);

    itemname = choosec(item,"鼻汁","倦怠感","咽頭痛","咳嗽","発熱");

    output;

  end;

end;

drop i iend;

run;

/*アイテム数*/

proc sql noprint;

 select count(distinct item) into: item_n trimmed from test;

quit;

%put &=item_n;

/*アイテム*/

proc sort data=test(keep=item itemname) out=itemlist nodupkey;

  by  item;

run;

/*アイテムformat*/

data itemformat;

 set itemlist;

 FMTNAME="itemformat";

 START=item;

 LABEL=itemname;

run;

proc format lib=work cntlin=itemformat;

run;

/*個人内での重複除去*/

proc sort data=test out=test1 nodupkey;

  by ID item;

run;

/*アイテム単純集計*/

proc summary data=test1 nway ;

 class item itemname;

 output out=_single_count(drop=_TYPE_ rename=(_FREQ_=single_count) );

run;

data single_count;

 set _single_count;

 single_start=0;

 run;

/*個人内での組み合わせ生成のため転置*/

proc transpose data=test1 out=out1 prefix=item_;

 by ID;

 id item;

 var item;

run;

/*組み合わせ生成*/

data out2;

length comb $200.;

 set out1;

 array ar item_:;

 do over ar;

  if ^missing(ar) then comb=catx("/",comb,ar);

 end;

  count=1;

run;

/*組み合わせ集計*/

proc summary data=out2 nway missing;

 class comb: item_:;

 output out=_comb_count(drop=_TYPE_ rename=(_FREQ_=comb_count) );

run;

 proc sort data=_comb_count;

 by descending comb_count;

run;

/*組み合わせ数と最大カウント*/

proc sql noprint;

 select count(*) into:comb_n from comb_count;

 select max(comb_count) into:comb_max from comb_count;

quit;

/*星取り表のつなぎ線用に最小・最大値*/

data comb_count;

set _comb_count;

 array dummy_{&item_n};

 do i = 1 to &item_n;

    dummy_{i}=i;

 end;

min = min(of item_:);

max = max(of item_:);

x = _N_;

drop i;

run;

/*グラフ用データ完成*/

data graph_data;

 set comb_count single_count;

run;

/*星取表*/

%macro scatterplot();

 %do i = 1 %to &item_n;

    scatterplot x=x y=item_&i /yaxis=y2 markerattrs=(size=12 symbol=circlefilled);

    scatterplot x=x y=dummy_&i /yaxis=y2 datatransparency=0.9 markerattrs=(size=12 symbol=circlefilled);

 %end;

%mend;

proc template ;

  define statgraph upset;

      begingraph ;

       /*2*2レイアウト*/

          layout lattice/ rowdatarange = union

                             columnweights = (0.25 0.75)

                             rowweights = (0.6 0.4)

                             columns = 2 rows = 2 

                             pad=(top=0 bottom=0 left=0 right=0) 

                             ; 

              /*余白 */

              layout overlay / pad=(top=0 bottom=0 left=0 right=0)

                                     outerpad=(top=0 bottom=0 left=0 right=0) 

                                     xaxisopts  = (display =none);

              endlayout;

             

              /*ヒストグラム（barchartparm）*/

              layout overlay /pad=(top=0 bottom=0 left=0 right=0) 

                                   outerpad=(top=0 bottom=0 left=0 right=0)

                                    walldisplay=none

                                    xaxisopts  = (display =none linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)))

                                    yaxisopts  = (display= (line tickvalues) linearopts = (tickvaluesequence = (start=0 end=&comb_max. increment=5))) ;

                                   barchartparm category=x response=comb_count/datalabel=comb_count;

              endlayout;

              /*横ヒストグラム(highlowplot)*/

              layout overlay /pad=(top=0 bottom=0 left=0 right=0)

                                    outerpad=(top=0 bottom=0 left=0 right=0)

                                    walldisplay=none

                                    xaxisopts  = (display= none reverse=true)

                                    yaxisopts = (display =none linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)))

                                    y2axisopts = (display =(tickvalues) tickvalueattrs=(size=9) linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)));

                                    highlowplot y=item high=single_count low=single_start 

                                    /yaxis=y2 group=item type=bar barwidth=0.5 fillattrs=(color=VLIGB) outlineattrs=(color=VLIGB) 

                                     highlabel=single_count labelattrs=(color=black)

                                    ;

              endlayout;

              /*星取表(scatterplot とhighlowplot)*/

              layout overlay  /pad=(top=0 bottom=0 left=0 right=0)

                                     outerpad=(top=0 bottom=0 left=0 right=0)

                                     xaxisopts  = (display =none linearopts = (tickvaluesequence = (start=1 end=&comb_n. increment=1)))

                                     yaxisopts  = (display= (line ) linearopts = (viewmin=1 viewmax=&item_n  ))

                                     y2axisopts  = (display=none linearopts = (viewmin=1 viewmax=&item_n  ))

                                     walldisplay=none;

                                     %scatterplot();

                                     highlowplot x=x low=min high=max/yaxis=y2;

              endlayout;

           endlayout;

      endgraph;

  end;

run;

* グラフ作成実行 ;

proc sgrender data=graph_data template=upset;

format item itemformat.;

run;

IMLを使わずFCMPで行列計算⑦ コレスキー分解

proc fcmp;

   array mat1[3,3] 2 2 3 2 4 2 3 2 6;

   array mat2[3,3] / nosym;

   call chol(mat1, mat2, 0);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("MAT2", mat2 ) ;

run;

proc print data=MAT1;

run;

proc print data=MAT2;

run;

mat2の転置行列を作って，mat2とmat2tの行列積をつくればmat1に戻ることが確認できる
つまり，コレスキー分解は　下三角行列とその転置の積に分解する操作

 

IMLを使わずFCMPで行列計算⑥ 行列式の算出

 proc fcmp;

   array mat1[2,2] (2 1 3 4);

   call det(mat1, result);

   put result=;

quit;

result = 5

2次正方行列なので

２ * 4  - 1 * 3 = 5

IMLを使わずFCMPで行列計算⑤ 単位行列の生成

proc fcmp;

  array mat1[3,3] /nosym;

   call identity(mat1);

   rc = write_array("RESULT", mat1 ) ;

quit;

 

IMLを使わずFCMPで行列計算④ 転置行列

 proc fcmp;

   array mat1[3,2] (3,5,0,9,1,2);

   array trans_mat1[2,3] /nosym;

  /*call transposeでmat1の転置行列をtrans_mat1に格納 */

   call transpose(mat1,trans_mat1);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("TRANS_MAT1", trans_mat1 ) ;

quit;

proc print data=MAT1;

run;

proc print data=TRANS_MAT1;

run;

IMLを使わずFCMPで行列計算③ 逆行列

 proc fcmp;

   array mat1[3,3] (3,5,0,9,1,2,8,7,9);

   array inv_mat1[3,3] /nosym;

   array result[3,3] /nosym;

  /*call invでmat1の逆行列をinv_mat1に格納 */

   call inv(mat1,inv_mat1);

  /*mat1と逆行列の行列積をresultに格納*/

   call mult(mat1, inv_mat1, result);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("INV_MAT1", inv_mat1 ) ;

   rc = write_array("RESULT", result ) ;

quit;

proc print data=MAT1;

run;

proc print data=INV_MAT1;

run;

proc print data=RESULT;

run;

浮動小数点誤差で-0とかってなってますが，逆行列かけたら　単位行列になるのが確認できる