平均の信頼区間を求める話_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;

annotateマクロの%systemの話

gplotなりsgplotなりグラフを作る時に、annotateマクロを使用することはあると思います。
数あるannotateマクロの中に、%system(hoge,hoge,hoge);があります。今回はそのお話
前にあげた<http://dengonmemo.blogspot.jp/2016/09/sas_29.html>と関連しています。
私も正しく理解しているかは怪しいところなので、間違っているかもしれませんです。

このマクロはannotateで文字やら線やらをグラフ内に書くときに、書く場所の設定をするマクロです。
そもそもsasで出したグラフには、data Area、Graphics Area、output Areaの3つがあります。
その三つを図で書くとだいたい下記のイメージです。(それぞれの斜線部分がそのAreaに対応)

data area(実際にグラフの値があるエリア)
|-----------------------|
|   |-----------------|   |
|   |                       |   |
|   |  ↑ /////////////    |   |
|   | ｜/////////////    |   |
|   | ｜/////////////    |   |
|   | ｜－－－→  |   |
|   |----------------  |   |
|                               |
|---------------------- |

Graphics Area(画面全体)
|---------------------|
|/////////////////////////|
|/|-----------------|///|
|/|/////////////////////|///|
|/|/ ↑/////////////////|///|
|/|/｜////////////////|///|
|/|/｜////////////////|///|
|/|/｜－－－→//|///|
|/|/////////////////////|///|
|/|------------------|///|
|//////////////////////////|
|----------------------|

Output Area(全体からheaderとfooterを除いたところ)
|-----------------------|
|                               |
| |-------------------|   |
| |///////////////////////|   |
| |/↑////////////////////|   |
| |/｜//////////////////|   |
| |/｜//////////////////|   |
| |/｜－－－→////|   |
| |///////////////////////|   |
| |-------------------|   |
|                               |
|----------------------- |

%systemの引数に1,2を指定するとData Area、3,4でGraphics Area、5,6でoutput Areaに
以降で使用するannotateマクロが出力されるようになります。
また1,3,5を指定すると百分率、2,4,6を指定すると座標で、以降のannotateマクロでの座標が指定できるようになります。

例えば%system(2,2,3)で、グラフのdata Areaの座標でannotateマクロを指定できます。
この下に%line( 1 , 30 , 1 , 28.5   , black, 1  , 0.15   );とすると
[1,30]から[1，28.5]に線を引きます。
%system(1,2,3)と記載した際、第一引数が1なので出力が百分率扱いなので、
この下のannotateマクロの指定で100より大きい値は指定できません。(100%が最大値)

また以下の様に一つのデータステップで複数回%systemは指定できます。
指定のたびに設定が切り替わります。
data anno;
    %system(2,2,3);
    *-----lineのX,Yは座標で指定できる;
    %line(~~);

　 *-----lineのXは百分率(%)、Yは座標で指定できる;
    %system(1,2,2);

    %line(~~);
run;

%systemの第三引数は私も良くわかっていないので、今は扱っていません。
どなたか教えていただければ…