wilcoxon順位和検定の話_npar1wayプロシジャ

対応の無い2群間の比較をする際，wilcoxon順位和検定を実行するときはあると思います．
この検定をSASで実行するには，npar1wayプロシジャのwilcoxonオプションを用います．
この検定は「2群間の比較」であることに注意してください．

univariateプロシジャでも2群の比較をすることはできますが，
univariateプロシジャではwilcoxonの「符号」順位和検定を行うことが出来ます．
この検定は今回のwilcoxon順位和検定と違い，対応のある2群間比較になるので割愛します．

/*---------- testdata ----------*/
data hoge ;
    cat = "A"; var1 = 1 ; output ;
    cat = "A"; var1 = 1 ; output ;
    cat = "A"; var1 = 1 ; output ;
    cat = "A"; var1 = 2 ; output ;
    cat = "A"; var1 = 2 ; output ;
    cat = "A"; var1 = 3 ; output ;
    cat = "A"; var1 = 3 ; output ;
    cat = "A"; var1 = 3 ; output ;
    cat = "A"; var1 = 3 ; output ;
    cat = "A"; var1 = 3 ; output ;
    cat = "A"; var1 = 4 ; output ;
    cat = "A"; var1 = 4 ; output ;
    cat = "A"; var1 = 4 ; output ;
    cat = "A"; var1 = 5 ; output ;
    cat = "A"; var1 = 5 ; output ;

    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 1 ; output ;
    cat = "B"; var1 = 2 ; output ;
    cat = "B"; var1 = 2 ; output ;
    cat = "B"; var1 = 3 ; output ;
    cat = "B"; var1 = 4 ; output ;
    cat = "B"; var1 = 5 ; output ;
    cat = "B"; var1 = 6 ; output ;
    cat = "B"; var1 = 6 ; output ;
run ;

proc npar1way data=HOGE wilcoxon ;
      class CAT ;
      var VAR1 ;
      output out = NPAR ;
run;

正規分布で近似した時の片側P-値は「PR_WIL」
同じ時の両側P-値は「P2_WIL」の変数に格納されています．

結果のデータセットの全体像








おまけですが，npar1wayプロシジャを動かした時の結果の名前の一覧を以下に貼ります．
ods trace on した時に出てくるやつです．
上のプログラムでoutput outで指定すると全ての結果がデータセット化されますが，
結果の名前を指定してods outputすると，データセット化する結果を選ぶことが出来ます．

平均の信頼区間を求める話_proc means

なんかの測定値の平均の95％信頼区間を求めたいときは往々にしてあると思います．

あちこちに載っていますが，そんな時はproc　meansで出せますよ，という話です．
proc univariateでも出せますが，結果をデータセット化するのはproc meansのほうが簡単だと思います．

meansでlclm(信頼区間の下限)，uclm(信頼区間の上限)を指定するだけです．簡単ですね．

/*----- test data -----*/
data aa ;
    a = 1 ; b = 1 ; output ;
    a = 2 ; b = 1 ; output ;
    a = 3 ; b = 1 ; output ;
    a = 4 ; b = 1 ; output ;
    a = 5 ; b = 1 ; output ;
    a = 6 ; b = 1 ; output ;
    a = 9 ; b = 2 ; output ;
    a = 8 ; b = 2 ; output ;
    a = 7 ; b = 2 ; output ;
    a = 5 ; b = 2 ; output ;
    a = 6 ; b = 2 ; output ;
run ;

proc means data = AA noprint nway ;
    var A ;
    class B ;
    output out = HOGE n = lclm = uclm = / autoname ;
run ;

マクロ変数を作成するときに引っかかった話

マクロ変数の値を呼び出すときに、思ってたのと違う値が出てきた時があったのでそれを紹介する話です。

マクロ変数をローカルマクロとして定義した後、同じマクロ変数名でグローバルマクロ変数としてもう一度定義して、そのマクロ変数を呼び出そうとしました。

私のイメージではグローバルマクロ変数として定義したほうの値が取れてくると思ったのですが、先に定義したはずのローカルマクロ変数の値が取れてくる時がありました。
とりあえず以下にプログラムを書いときます。
しかし不思議ですねえ

data DT_TEST ;
    _AA=1 ; output ;
    _AA=2 ; output ;
    _AA=3 ; output ;
run ;
%macro MCR_TEST ;

    *----------  ローカルマクロを指定;
    %let _aa2=123 ;
    data _NULL ;
        set DT_TEST ;
       
        *----------  グローバルマクロを指定;
        call symputx("_aa2",_N_,"g") ;
    run ;
    *----------  マクロ変数を呼び出し;
    %put &_aa2 ;

%mend ;

*----------  1回目:ローカルマクロとしての値(123)が返る;
%MCR_TEST ;


*----------  2回目以降:グローバルマクロとしての値(3)を返す;
%MCR_TEST ;

whereステートメントでデータセットを抽出するときの話

データセットをwhereステートメントを使って抽出するとき、
コピペミスをして一つのデータステップにwhereステートメントを2回書いてしまいました。
実行したらどうなるのだろうと思って実行すると、

data hoge;
    a = 1; output;
    a = 10; output;
    a = 100; output;
run;

/*これ↓*/

data hogehoge;
    set hoge;
    where a = 1;
    where a = 10;
run;












と、1つ目のwhereが二つ目の式に書き換えられるんですね
結果二つ目のwhereステートメントのみで絞ったのと同じ結果になる、と。

じゃあと思ってsetのオプションにwhere書いたらどうなるのだろうと思い

data hogehoge;
    set hoge(where = (a = 1));
    where a = 10;
run;
とすると、




を実行すると、これはwarningになるのですね

今まで書いたことなかったので知りませんでした。
知ってて当然、なのかもしれませんが。

信頼区間を求めるお話_FREQプロシジャ

変数の値が0,1の二値を取る時に信頼区間を求めるお話です

信頼区間を求めたい変数に、0と1の両方があれば良いのですが
値として0しかない、1しかない時は往々にしてあると思います。
そんな時はweightステートメントにzeroesオプションを付けましょうと言うお話です。

tabelsステートメントの(level='1')は変数の値が1のものをpositiveとして計算してください、の意味です。
これも良く忘れて、二値データのどっちがpositiveだっけ？となる…気がします

/* ---------- ここからプログラム --------- */
*---------- test data/水準0がなく、1しかないデータ;
data hoge;
    a = 1; output;
    a = 1; output;
    a = 1; output;
run;

*---------- 後でweight考慮するために集計;
proc freq data = HOGE  noprint;
    tables a / out = HOGE_FREQ ;
run;

*---------- dummy data;
data DUM;
    a = 0; COUNT = 0; output;
    a = 1; COUNT = 0; output;
run;

*---------- データに無い水準は0件としてデータに持たせる;
data HOGE2;
    merge DUM HOGE_FREQ ;
    by A;
    drop PERCENT;
run;

*---------- 信頼区間/正確な信頼区間はXL_BIN-XU_BIN;
proc freq data = HOGE2  noprint;
    tables a / binomial(level = '1');
    exact binomial;
    weight COUNT / zeroes;
    output out = TEST_OUT;
run;

実行日時を取得する話

月1の更新は維持できると思っていた、そんな時代も今は昔

実行した日時を取得するプログラムを今まで思考停止して使っていましたが
よく見ると複雑な書き方をしていたので書き直してみました。
実行結果は変わらず、対して実行時間も変わらないので完全に自己満足の世界ですね
年月日_日時の形式で実行日時が取得できます。

/*以下が以前私が使っていた日時取得プログラム*/
data _null_;
    length _DAY DAY TIME DATETIME $20;
   
    now = datetime();                                                        
    _DAY = put(datepart(now) , yymmdd10.);                                   
    DAY  = compress(_DAY,"-");                                               
   
    TIME = compress(put( timepart(now) , hhmm5.) , ':');
   
    DATETIME = trim( left(DAY) )||'_'||trim( left(TIME) );
    call symput('today', trim( left(DATETIME) ) );
run;
/*--ここまで以前のもの-----*/

/*--ここから新しく作ったもの--------*/
data _null_;
    length DAY TIME $10;
    DAY  = put(today() , yymmddn8.);
    TIME = compress(put(time() ,tod5.) , ,"dk");
    call symputx("_today" , catx("_" , DAY , TIME) );
run;

複数のグラフを1枚にまとめる話_今回は4枚のグラフを一つに

何枚かのグラフを一枚の紙に出力したいときのお話です。
グラフをまとめるのはgreplayプロシジャを使います。
GTLさんには座っててもらいましょう。

以下のプログラムは
・1ページに4枚のグラフを並べて表示
・IDごとに改ページ
・一つのIDの中でもグラフ4枚ごとに改ページ
したRTFファイルを出力するものです。

一番下のods rtfで指定した場所にhoge.rtfが出てきます。

以下の3枚の画像が出力イメージです。






--ここからプログラム-----

/* ---------- testdata --------- */
data test;
    id = 1; cd = 1 ; time = 1; value = 1; output;
    id = 1; cd = 1 ; time = 2; value = 1; output;
    id = 1; cd = 1 ; time = 3; value = 1; output;
    id = 1; cd = 1 ; time = 4; value = 1; output;
    id = 1; cd = 1 ; time = 5; value = 1; output;
    id = 1; cd = 2 ; time = 4; value = 2; output;
    id = 1; cd = 2 ; time = 5; value = 2; output;
    id = 1; cd = 2 ; time = 6; value = 2; output;
    id = 1; cd = 2 ; time = 7; value = 2; output;
    id = 1; cd = 2 ; time = 8; value = 2; output;
    id = 1; cd = 3 ; time = 1; value = 3; output;
    id = 1; cd = 3 ; time = 2; value = 3; output;
    id = 1; cd = 3 ; time = 3; value = 3; output;
    id = 1; cd = 3 ; time = 4; value = 3; output;
    id = 1; cd = 4 ; time = 4; value = 4; output;
    id = 1; cd = 4 ; time = 5; value = 4; output;
    id = 1; cd = 4 ; time = 6; value = 4; output;
    id = 1; cd = 5 ; time = 2; value = 5; output;
    id = 1; cd = 5 ; time = 6; value = 5; output;
    id = 1; cd = 5 ; time = 7; value = 5; output;
    id = 1; cd = 5 ; time = 9; value = 5; output;
    id = 2; cd = 1 ; time = 1; value = 1; output;
    id = 2; cd = 1 ; time = 2; value = 1; output;
    id = 2; cd = 1 ; time = 3; value = 1; output;
    id = 2; cd = 1 ; time = 4; value = 1; output;
    id = 2; cd = 2 ; time = 3; value = 2; output;
    id = 2; cd = 2 ; time = 4; value = 2; output;
    id = 2; cd = 2 ; time = 5; value = 2; output;
    id = 2; cd = 3 ; time = 1; value = 3; output;
    id = 2; cd = 3 ; time = 2; value = 3; output;
    id = 2; cd = 3 ; time = 3; value = 3; output;
    id = 2; cd = 3 ; time = 4; value = 3; output;
run;

/* ---------- 出力 --------- */
  goptions gunit=pct;a
  goptions
      reset   = all
      gsfmode = replace
      xmax    = 6 in
      ymax    = 9 in
      vsize   = 5 in
      hsize   = 9 in
  ;
*-========================================;
*-グラフ出力位置指定(4枚)*;
  proc greplay gout=work.gseg tc=work.tmplt nofs;
    tdef PK1
   
    1 / llx=0 lly=50                                                            /*left upper panel*/
        ulx=0 uly=100
        lrx=50 lry=50
        urx=50 ury=100
   
    2 / llx=50 lly=50                                                           /*right upper panel*/
        ulx=50 uly=100
        lrx=100 lry=50
        urx=100 ury=100
   
    3 / llx=0 lly=0                                                             /*left lower panel*/
        ulx=0 uly=50
        lrx=50 lry=0
        urx=50 ury=50
   
    4 / llx=50 lly=0                                                            /*right lower panel*/
        ulx=50 uly=50
        lrx=100 lry=0
        urx=100 ury=50
   ;
  run;
  quit;

goptions gunit=pct ;
%macro grepout(_indt = , _wh = );
*---------- annotate macroを有効にする;
%annomac;

*---------- 各グラフの出力を一旦抑制/まとめたグラフ以外は出力不要;
goptions nodisplay;

*---------- 種類をマクロ変数に格納;
data FDAT;
    set &_indt.(where = (&_wh.)) end = _EOF;
    if _eof = 1 then do;
        CD4 = ceil(CD/4);
        call symputx("_n" , CD);
        call symputx("_n4" , CD4);
    end;
    drop CD4;
run;

*---------- CDの種類だけGPLOTを実行;
%do i = 1 %to &_n ;

    *---------- nameを取得/グラフにnameを表示する;
    data _NULL_;
        set TEST(where = (CD = &i. and &_wh.));
        call symputx("_name" , CD);
        call symputx("_id" , ID);
    run;

    ods escapechar = '^' ;

    data anno1;
        set TEST(where = (CD = &i.));
        length TEXT $80;
        %dclanno;
        %system(3,3,3);
 
        *---------- 軸ラベル;
        %label(20, 95  , "CDが&_name.のグラフ" , black, 0, 0, 4 , 'Arial'  , 5);
        %label(80, 95  , "IDが&_id.のグラフ" , black, 0, 0, 4   , 'Arial'  , 5);
    run;

    proc gplot data = FDAT(where = (CD = &i.)) anno=anno1;
 
        plot  VALUE  * TIME   / nolegend vaxis = axis1 haxis = axis2 skipmiss ;      
        symbol1 mode = include v = triangle h = 2 c = black  i = join l = 21 w = 1.5;

        axis1                                                                   /*y軸目盛*/
          label=none
          offset=(2,2)  minor=none major=(w=1.5 h=0.8)
          length=60 width=1 order=(0 to 6 by 1)
          value=(font='Times New Roman' h=4 )
          origin=(10.3, 30);

        axis2                                                                   /*x軸目盛*/
          label=none
          offset=(2, 2)  minor= none  major=(w=1.5 h=0.8)
          length=80 width=1 order=( 1 to 9 by 1 )
          value=(font='Times New Roman' h=4 )
          origin=(10.3, 30);
    run;
    quit;
%end;

*---------- グラフの出力を再開_4グラフ1ページのものだけを出力;

goptions hsize = 9 in display;

*-出力;
    %macro greplay/*(path=, out=)*/;
      %put Number of Parameter is &_n;
      %put Number of Page is &_n4.;
      %do i = 1 %to &_n4. %by 1;
       
        %if %eval(&i) = 1 %then %let j = 1;
                          %else %let j = %eval(&j + 4);
        %put i is &i.;
        %put j is &j.;
        proc greplay igout=work.gseg gout=work.gseg tc=work.tmplt nofs;
          template PK1;
          %if &_n4. = 1  %then %do;
            treplay 1:gplot
                    2:gplot1
                    3:gplot2
                    4:gplot3
            ;
          %end;
          %if &_n4. >= 2 and &i = 1 %then %do;
            treplay 1:gplot
                    2:gplot1
                    3:gplot2
                    4:gplot3
            ;
          %end;

          %if &_n4 >= 2 and &i >= 2  %then %do;
            treplay 1:gplot%eval(&j - 1)
                    2:gplot%eval(&j)
                    3:gplot%eval(&j+1)
                    4:gplot%eval(&j+2)
            ;
          %end;
        quit;
      %end;
    %mend greplay;
    %greplay

    *---------- 出力したグラフをリセット;
    proc catalog  catalog = work.Gseg kill  ;
        run;
   
    quit ;
%mend grepout;

ods listing close;
ods rtf file = "hogehoge\hoge.rtf";

%grepout(_indt = TEST , _wh = %nrstr(ID = 1) );
%grepout(_indt = TEST , _wh = %nrstr(ID = 2) );

ods rtf close;
ods listing;