2018年12月1日土曜日

日時分の文字値から日付をSAS日付で取り出す

"2018-11-15T10:45"の様な形式の日付の文字変数を,日付部分だけSAS日付に変換する時ってあると思います.
今の例で言うと,"2018-11-15T10:45"文字値を2018-11-15のSAS日付の変数に変えると言うだけです.

以下にプログラムと実行ログを載せているので他に書くこと無いのですが,SAS9.3以下だとおかしくなることがあります.

もともとの文字変数側のlengthが短いと正しくSAS日付に変換できないことがあるようです.
SAS9.4だとこんなことは無く何も気にしなくて大丈夫です.


/*---------- プログラム ----------*/

data HOGE ;
    A = "2018-11-15T10:45" ; output ;
    A = "2018-11--T10:45" ; output ;
 
run ;

data _null_ ;
    set HOGE ;

    if count(A , "-") = 2 then B = input(A,is8601da.) ;
    else putlog A ;
   
    format B yymmdd10. ;
    putlog A B ;

run ;

/*---------- 実行log ----------*/

356
357  data HOGE ;
358      A = "2018-11-15T10:45" ; output ;
359      A = "2018-11--T10:45" ; output ;
360
361  run ;
NOTE: データセットWORK.HOGEは2オブザベーション、1変数です。
NOTE: DATAステートメント処理(合計処理時間):
      処理時間           0.01 秒
      CPU時間            0.00 秒

362
363  data _null_ ;
364      set HOGE ;
365      if count(A , "-") = 2 then B = input(A,is8601da.) ;
366      else putlog A ;
367
368      format B yymmdd10. ;
369      putlog A B ;
370
371  run ;
2018-11-15T10:45 2018-11-15    /*1レコード目をputlog*/
2018-11--T10:45                          /*Aにハイフンが2個以上あるのでelse以下のputlog*/
2018-11--T10:45 .                        /*2レコード目をputlog*/
NOTE: データセットWORK.HOGEから2オブザベーションを読み込みました。
NOTE: DATAステートメント処理(合計処理時間):
      処理時間           0.00 秒
      CPU時間            0.00 秒

2018年11月1日木曜日

MMRMをSASで実行する話_proc mixed


固定効果とランダム効果を持つ混合モデルの推定をするmixedプロシジャの紹介です.
あくまでプロシジャの紹介なので,混合モデルとは,固定効果とは,と言ったところはあまり触れません.
触れるほど知識が私に無いとも言います.

まあ何はともあれデータとプロシジャのプログラムを以下に置きます.
データセットPR2は,PERSONにID,TRTPに群,Yに架空の測定値,AVISITNに測定時点を持っています.
データはSASのonline docから持ってきて改変しています.
SAS/STAT14.1,BASE9.4の環境です.

このデータに対して,測定値Yは群TRTPと時点AVISITNで説明できる,というモデルを想定します.
本来ならmixedプロシジャは欠測があっても実行できますが,今回は欠測なしのデータを用意しています.
群と時点は固定効果としています.このデータの中の水準のみに興味があるからです.
この二つをランダム効果とすると,他にも群や時点があり,
その中からランダムで選ばれたもののみが格納されたデータということになります.
そんなことは無いので固定効果です.

またランダム効果は設定していません.ただしIDであるPERSONはランダム効果とすることは出来ます.
なぜならこのデータは,母集団からランダムで選ばれた人たちから構成されていると考えられるからです.


data PR;
   input PERSON TRTP $ y1 y2 y3 y4;
   y=y1; output;
   y=y2; output;
   y=y3; output;
   y=y4; output;
   drop y1-y4 ;
   datalines;
 1   A   21.0    20.0    21.5    23.0
 2   A   21.0    21.5    24.0    25.5
 3   A   20.5    24.0    24.5    26.0
 4   A   23.5    24.5    25.0    26.5
 5   A   21.5    23.0    22.5    23.5
 6   A   20.0    21.0    21.0    22.5
 7   A   21.5    22.5    23.0    25.0
 8   A   23.0    23.0    23.5    24.0
 9   A   20.0    21.0    22.0    21.5
10   A   16.5    19.0    19.0    19.5
11   A   24.5    25.0    28.0    28.0
12   B   26.0    25.0    29.0    31.0
13   B   21.5    22.5    23.0    26.5
14   B   23.0    22.5    24.0    27.5
15   B   25.5    27.5    26.5    27.0
16   B   20.0    23.5    22.5    26.0
17   B   24.5    25.5    27.0    28.5
18   B   22.0    22.0    24.5    26.5
19   B   24.0    21.5    24.5    25.5
20   B   23.0    20.5    31.0    26.0
21   B   27.5    28.0    31.0    31.5
22   B   23.0    23.0    23.5    25.0
23   B   21.5    23.5    24.0    28.0
24   B   17.0    24.5    26.0    29.5
25   B   22.5    25.5    25.5    26.0
26   B   23.0    24.5    26.0    30.0
27   B   22.0    21.5    23.5    25.0
;
data PR2 ;
    set PR ;
    by PERSON ;
    retain AVISITN 0 ;
    if first.PERSON = 1 then AVISITN = 1 ;
    else AVISITN = AVISITN + 1 ;
run ;

proc mixed data=pr2 method=reml  ;
   class PERSON TRTP AVISITN ;
   model Y = AVISITN TRTP TRTP*AVISITN / s ddfm = kr ;
   repeated AVISITN / type=un subject=PERSON r;
   lsmeans TRTP * AVISITN / pdiff cl ;
   ods output diffs = DIFF ;
   ods output lsmeans = LSMEAN ;
   ods output ConvergenceStatus = CHECK ;
run;


上のプログラムを実行すると何か結果は出てきます.収束もします.
出力結果が多いので紹介は省きます.すみません.

mixedプロシジャのオプションにmethodがありますが,ここで推定方法を指定します.remlは制限付き最尤法を示します.

classステートメントには分類変数を指定します.ここで指定した変数を他のステートメントでも使用できます.

modelステートメントには「説明したい変数=固定効果の変数」の形式で指定します.*で区切ると交互作用を示します.

  オプションのsはsolutionの略で,ddfmには自由度の計算方法を指定します.
  今は「kr」が指定されており,コレは「Kenward-Roger」です.
repeatedステートメントには,反復測定の変数を指定します.今回はAVISITが同じ症例に対しての繰り返しを情報を持っています.
  オプションのtypeには共分散構造を指定します.unは無構造を示します.
  構造の指定を間違えると結果がおかしくなることがあるので,あえて無構造としています.
  ただしこれはモデルが収束しないことがあるので,その際の取り扱いは事前に決めておくのが良いでしょう.
  subjectオプションには被験者を指定します.
lsmeansステートメントで最小二乗法を計算することを指定します.
  オプションのpdiffを指定すると最小二乗平均の差を出力してくれます.
今回はods outputで交互作用の交互作用の最小二乗平均の差(DIFF),交互作用の最小二乗平均(LSMEAN),収束状態(CHECK) を出しています.
ここはお好みで変えてください.

今は反復測定を想定したプログラムの紹介ですが,反復測定をしないときもmixedプロシジャは使えます.
その時はrepeatedステートメントの指定が要らなくなりますね.
また今回はランダム効果を想定してませんが,ランダム効果を考える時は
randomステートメントに対象の変数を指定してください.
おいおい紹介もしたいと思います.

最後におまけでods traceしたoutputの名前を置いておきます.

出力の追加:
-------------
名前 :         ModelInfo
ラベル :       モデルの情報
テンプレート : Stat.Mixed.ModelInfo
パス :         Mixed.ModelInfo
-------------
出力の追加:
-------------
名前 :         ClassLevels
ラベル :       分類変数の水準の情報
テンプレート : Stat.Mixed.ClassLevels
パス :         Mixed.ClassLevels
-------------
出力の追加:
-------------
名前 :         Dimensions
ラベル :       次元の数
テンプレート : Stat.Mixed.Dimensions
パス :         Mixed.Dimensions
-------------
出力の追加:
-------------
名前 :         NObs
ラベル :       オブザベーション数
テンプレート : Stat.Mixed.NObs
パス :         Mixed.NObs
-------------
出力の追加:
-------------
名前 :         IterHistory
ラベル :       反復履歴
テンプレート : Stat.Mixed.IterHistory
パス :         Mixed.IterHistory
-------------
出力の追加:
-------------
名前 :         ConvergenceStatus
ラベル :       収束状態
テンプレート : Stat.Mixed.ConvergenceStatus
パス :         Mixed.ConvergenceStatus
-------------
NOTE: 収束基準は満たされました。
出力の追加:
-------------
名前 :         R
ラベル :       PERSON 1 の推定 R 行列
テンプレート : Stat.Mixed.R
パス :         Mixed.R
-------------
出力の追加:
-------------
名前 :         CovParms
ラベル :       共分散パラメータ推定値
テンプレート : Stat.Mixed.CovParms
パス :         Mixed.CovParms
-------------
出力の追加:
-------------
名前 :         FitStatistics
ラベル :       適合度統計量
テンプレート : Stat.Mixed.FitStatistics
パス :         Mixed.FitStatistics
-------------
出力の追加:
-------------
名前 :         LRT
ラベル :       帰無モデルの尤度比検定
テンプレート : Stat.Mixed.LRT
パス :         Mixed.LRT
-------------
出力の追加:
-------------
名前 :         SolutionF
ラベル :       固定効果の解
テンプレート : Stat.Mixed.SolutionF
パス :         Mixed.SolutionF
-------------
出力の追加:
-------------
名前 :         Tests3
ラベル :       固定効果の Type 3 検定
テンプレート : Stat.Mixed.Tests3
パス :         Mixed.Tests3
-------------
出力の追加:
-------------
名前 :         LSMeans
ラベル :       最小 2 乗平均
テンプレート : Stat.Mixed.LSMeans
パス :         Mixed.LSMeans
-------------

SASで一元配置分散分析

glmプロシジャの紹介です.
ざっくり言うと平均からのズレを見たい時に使うと思います.
実行はSAS9.4,SAS/STAT 14.3で行いました.

quitをつけで実行をしてください.
色々ステートメントがありますが,必須なのはmodelステートメントだけです.
つまり他のステートメントは必要に応じて追加する形です.
またこのページでは以下のステートメントの全てを紹介はしません.すみません.

まず使うデータはSASHELPのFailureです.
これは失敗の数がCOUNTに,失敗した理由がCAUSEに,
失敗したときの方法がPROCESSに格納されています.

今たとえば,プロセスAとBで,どちらの失敗が多いかを調べたいとします.(1要因,2群の比較)
/*1要因,2群の比較*/
proc glm data = SASHELP.FAILURE outstat = OUT plots = none alpha = 0.05;
    class PROCESS ;                                /*群変数を指定*/
    model COUNT = PROCESS / ss3;     /*測定値 = 群 ss3オプションをつけているのでtype3のみ出力*/
    means PROCESS ;                      /*各群の平均値を出力*/
    ods output  FitStatistics = FIT ;   
    ods output  ModelANOVA    = MODEL ; 
    ods output  Means         = MEANS ; 
quit;

--model出力例----
 要因                    自由度          平方和        平均平方       F 値    Pr > F
 Process                      1     44.80000000     44.80000000       1.95    0.1670
--means出力例----
 水準                 ------------Count------------
 Process        N             平均         標準偏差
 Process A     35       4.62857143       6.01510703
 Process B     35       3.02857143       3.12000862

FitStatisticsに適合度統計量,ModelANOVAにanovaのモデル(今回はss3オプションを指定しているのでtype3のみ)ですね.
Meansに各群の値の平均値が出ています.
ods outputしているので,それぞれをデータセットに出力しています.

glmのところのOUTSTATオプションですが,ここに指定した名前のデータセットに結果のMODELを出力します.
ですのでそのデータセットの中身はModelANOVAとほとんど一致します.
違うところはOUTSTATで出力したデータセットには,p-値,F値が丸められていないものが出ます.
ERROR時の平方和も出してくれますが,平均平方は出してくれません.

一方ods outputのModelANOVAではp-値とF値は四捨五入された値が出ます.
こちらはMODELの平均平方もデータセットに出力してくれます.
proc glmのplotsオプションにnoneを指定しているので何もグラフが出ない設定です.
何も指定しないと2群の箱ひげ図が出てきます.
またproc glmのalphaオプションには有意水準を指定します.(0から1の間)
何も指定しなかったら0.05=有意水準95%が指定されます.
ここまでは2群の比較です.(群変数のPROCESSにはAとBの二つが格納されている状態)


3群以上の時は上記に加えて多重比較を行ってくれます.
仮にA,B,C群の3群があった時にAC,AB,BCの比較をします.
3群の比較はデータを変えて,SASHELPのcarsを使います.
このデータセットの内,車種を格納したTYPE,ガソリン1リットルあたりの高速道路走行距離を格納したMPG_Highwayを使います.
/*1要因,3群以上の比較(この例は6群)*/
proc glm data = SASHELP.CARS outstat = OUT plots = none;
    class TYPE ;                     
    model MPG_Highway = TYPE / ss3 ;  
    means TYPE / tukey bon ;           /*オプションにて群間の比較をする際の検定方法を指定(複数可)*/
    ods output  FitStatistics = FIT ;   
    ods output  ModelANOVA    = MODEL ; 
    ods output  CLDiffs       = DIFF ; 
quit;

--means出力例(一部)-----
  有意水準 0.05 で有意に差があることを *** で示しています。

       Type
       比較           平均の差     同時 95% 信頼限界
  Hybrid - Sedan       27.3702      20.4417   34.2987  ***
  Hybrid - Wagon       28.1000      20.8746   35.3254  ***
  Hybrid - Sports      30.5102      23.4133   37.6071  ***
  Hybrid - Truck       35.0000      27.6929   42.3071  ***
  Hybrid - SUV         35.5000      28.4407   42.5593  ***
  Sedan  - Hybrid     -27.3702     -34.2987  -20.4417  ***
  Sedan  - Wagon        0.7298      -1.5701    3.0297

出力結果は基本的には2群のときと同じです.
2群のときの内容に加えて群間の比較を行います.
その内容を,ods outputを使ってDIFFの名前でデータセットにしています.
どの群の比較をしたかをComparison,95%信頼区間をLowerCL/UpperCL,平均の差をDiffrence,
有意であればSignificanceに1を格納されています.有意でなければSignificanceには0が入ります.
output windowでは有意な組み合わせに***が付いてます.
今はmeansステートメントにtukeyオプションを指定していますが,このオプションは複数指定できます.
指定した方法で群間の比較を行うので,複数指定すると複数の方法で比較してくれます.
tukey(tukeyの方法),t(pair t test),bon(bonferroni)などです.

ここまでは1要因の比較でしたが,複数要因の比較もglmプロシジャは行えます.
さっきと同じCARSのデータを用います.群の数が多いのでwhereで少し絞っています.

proc glm data = SASHELP.CARS outstat = OUT plots = none;
    where TYPE in("Sedan","SUV","Sports") ;
    class TYPE ORIGIN ;
    model MPG_Highway = TYPE ORIGIN TYPE * ORIGIN  / ss3 ;  
    lsmeans TYPE * ORIGIN / pdiff cl ;
    ods output  FitStatistics = FIT ;   
    ods output  ModelANOVA    = MODEL ; 
    ods output  LSMeanCL      = LSCL ; 
quit;

--model出力例--
   要因            自由度     平方和         平均平方       F 値    Pr > F
   Type            2          3089.294384    1544.647192    96.25   <.0001
   Origin          2          229.974936     114.987468     7.16    0.0009
   Type*Origin     4          34.950460      8.737615       0.54    0.7032

--LSMeanCL出力例--
   Type      Origin            平均           95% 信頼限界
   SUV       Asia         21.680000       20.104362    23.255638
   SUV       Europe       18.700000       16.208698    21.191302
   SUV       USA          20.040000       18.464362    21.615638
   Sedan     Asia         29.968085       29.155512    30.780658
   Sedan     Europe       27.115385       26.223355    28.007414
   Sedan     USA          28.544444       27.714010    29.374878
   Sports    Asia         26.647059       24.736317    28.557800
   Sports    Europe       25.130435       23.487719    26.773151
   Sports    USA          24.222222       21.596159    26.848285

これは出力が多いので絞っています.
MODELにはmodelで指定した組み合わせのモデルが出力されます.ここで有意
LSCLには,組み合わせごとの95%信頼区間が出ます.LSMEANが平均,LowerCLが下限,UpperCLが上限です.
今回はType*Originのp-値が0.05より大きいので有意ではありません.ですので混合効果はありませんね.
Typeの時には有意であるとなっているので,TypeによってMPG_Highwayの値に差があると言えそうです.

正直指定できるオプションやらステートメントが多いので,リファレンス見るのが一番ですかね.

2018年10月1日月曜日

sgplotで使うsymbolの紹介

SGPLOTでグラフを書く時,群毎にグラフにsymbolを出すことがあると思います.
どんなsymbolが使えるかをすぐ忘れてしまうので紹介します.
よく使うのは丸と四角と三角ですかね.下の図に乗せています.

その他のはonline docを参照してください.
https://documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=9.4_3.3&docsetId=grstatproc&docsetTarget=p0i3rles1y5mvsn1hrq3i2271rmi.htm&locale=ja

/*test data*/
data HOGE ;
  do G = 1 to 6 by 1 ;
    do A = 1 to 3 ;
      B = G ;
      output ;
    end ;
  end ;
run ;

proc sgplot data = HOGE noautolegend ;
  styleattrs
    datasymbols=(Circle CircleFilled Square SquareFilled Triangle TriangleFilled)
    datacontrastcolors=(black)
    datalinepatterns=(solid)
  ;
  series X = A Y = B
    / group = G markers
  ;
 
run ;


2018年9月6日木曜日

Phrma SugのSDEに行って来た話

東京でやると言うもんだから早々に行くことを決意,さっさと予約を済ませる.
しかしよく見ると受付が9:30~10:00,関西から行くにはいささか時間が早いので夜行バスで向かうとする.
結果的にここで夜行バスを選択したのが吉と出た.

迎えた当日,台風「こんにちはー」
SASユーザー総会の時も台風と一緒に東京に行った記憶がよみがえる.
すわ電車運行見送りかと震えるも私は夜行バス.前日に東京に入っていたので大した影響は受けなかった.

建物に入ればこちらのもの.恙無くSDEが始まった.
プロジェクタが何を思ったか赤を青に,青を赤に表示している所為でSASのアイコンが赤い.違和感違和感.
中国の規制当局の話が一番興味深かったが残念英語がイマイチわからない.
うーんこのと聞きつつ,英語の聞き取りは練習したほうがいいなあと実感した.
後はまあ人工な知能をpythonで実装した話の事例である.実に最近のブーム.
RでPPTへの出力,との内容も面白そうではあるがバリデーションの壁は厚く高そう.

うちの大将のポスターセッションも見てきたがやはり大将は話がうまい.
あれくらい話せるようになればなあと思うも,練習あるのみか.
いや人前で話すのは苦手なんだが

終わる頃には風が荒れる空模様だったが雨は無く,後楽園で知り合いと飲んで事故も無くホテルに引っ込む.
どうせ台風で交通機関荒れるからと一泊して本当に良かったと思いつつ12:00頃に就寝.
風が轟くもなんのその,すぐに夢の世界へとさようならし,翌日の昼過ぎに東京に別れを告げた.

自宅の隣家のブロック塀が見事に崩落している様に絶句するも,家は無事で何よりであった.
次は来年のSASユーザー総会まで何も無いかなあ,しかし東京に行くのはお金がかかって仕方が無い.
今回のSASグッズはSASのリングNOTEだった.結局もったいなくて使えずそのまま放置するんだろうなあ.

2018年9月1日土曜日

プログラム実行時にログウインドウを常に前面に表示する話

何かプログラムを実行するとログが表示されますよね.
プログラムを書く時は拡張エディタを大きくしたいけど実行時にはログ画面を大きくしたい...
そんなとき,あると思います.
パソコンの画面が小さいと,拡張エディタを大きくするとログ画面がエディタの後ろに行ってしまうので,実行時にログ画面をクリックして最前面に持ってくる謎の1手間が必要になります.
画面が十分に大きければ問題ないんですがねえ...
以下のオプションを設定すると,プログラムを実行する時に
自動でログ画面が最前面に出てきてくれる様なります.
実行時に勝手にログ画面が前に出てきてくれるので,ログの確認時に地味に助かっています.
ログ確認後に再度プログラムを修正する時は拡張エディタを選択しないといけませんが.
  1. logウィンドウをクリックして選択する
  2. ツール→オプション→ログ→「表示タブの新しいデータが表示されたときにウインドウを前面にする」
にチェックを入れる
これでOKです.個人的にはなかなか気に入っているオプションです.

2018年8月15日水曜日

プログラムをデバッグする話

データステップにデバッグオプションを付けると,デバッガが起動します.
1stepごとに実行し適宜変数の値を確認したり,指定した箇所まで実行したりすることが出来ます.
excelで言う所のデバッグトレースですかね.

以下のプログラムを実行すると,データステップにdebugオプションが付いているので
データステップデバッガが起動します.

data hoge  ;
    val1 = 1 ; val2 = 2 ; type = "A" ; output ;
    val1 = 1 ; val2 = 2 ; type = "B" ; output ;
    val1 = 1 ; val2 = 2 ; type = "C" ; output ;
run ;

data hoge2 / debug ;
    set hoge ;
    if TYPE = "A" then VAL3 = VAL2 - VAL1 ;
    if TYPE = "B" then VAL3 = VAL2 + VAL1 ;
    if TYPE = "D" then VAL3 = VAL2 * VAL1 ;
run ;





 





  















デバッガが起動すると,debugger log とdebugger sourceの二つのwindowが立ちあがります.

debugger logの下部にある点線の下にコマンドを入力し,デバッグを行います.
debugger sourceの黒線部分が今実行する直前のステートメントです.

ココで出来ることは,以下の感じです.
もちろん他にも機能はありますが,とりあえずこれだけです.
  • 指定した行数分ステートメントを実行(stepコマンド)
  • 指定した変数の値を確認(examineコマンド)
  • ステートメントを指定して,そのステートメントまで実行(breakコマンド)
  • 指定した変数の値が変わるまで実行(watchコマンド)
  • デバッガを終了(quitコマンド)
breakコマンドはexcelで言う所のbreakpointとほぼ同じだと思います.
 デフォルトではエンターキーがstep1に対応しているので,
debugger logでエンターをたたくと1ステートメントだけ実行されます.
以下の画像ですと134行目に黒線があるので,133行目まで実行されている状態です.
この状態でexamineコマンドでVAL3を見ると,値が1になっています.
examineコマンドは黒線がある直前まで実行された段階での変数の値を返してくれます.
ですのでdebugger logにステップ:行132カラム24の時にすべての変数の値を確認すると,VAL3は欠測になります.
(まだ132行目のTYPEが"A"の時にVAL3を演算する処理が実行されていないため.)
 _N_=1のため,今は1レコード目を読みこんでいる状態でデバッグをしていることが分かります.

データステップデバッガも例にもれず,1レコード読みこんで,すべてのステートメントを実行して,
また次のレコードを読みこんで...と進んでいきます.
確認したいレコード番号が大きいとそこまで実行しないといけません.
いちいちエンターキーを叩いて1レコードごとに実行しているとキリが無いので,
watch,breakなどを駆使して確認箇所まで実行するのが良いのではないでしょうか.



















マクロのデバッグにも使える感じですね.

参考資料
<http://www.sas.com/offices/asiapacific/japan/service/help/pdf/lebaseutilref.pdf>
dataステップデバッガの項目

SASユーザー会行ってきました

今年もSASユーザー会行ってきました.
台風がかすっていきましたが,まあ何とか全日程が完了してよかったです.

1日目は台風とともに東京に向かうことに.
前日入りも覚悟していたところでしたが,東にそれたと聞いて当日東京に向かうことに決定.
霧雨の舞う中会場にたどり着き,荷物はびしょびしょ.
雨が風で舞う所為で傘がその機能を果たしてくれず,それはそれはぬれました.
一日目はルームDを中心に発表を聞く.
ODS TABLE,rtf作成時のfont調整等々なかなか便利そうなものばかり.
fontの調整はあの発表をそのままパクってしまおうかな...

セッションが午後のみだったのであっという間に一日目は終了,懇親会へとなだれ込む
雨の中の移動にエネルギーを思ったより使ったのか,ご飯を見るとおなかがみるみる空いてくる.
乾杯の挨拶中に食べたいものにあたりをつけ,乾杯後に早速ご飯をget
周りを見ると先輩が知り合いを連れてきて後輩が群れを成して名刺交換にいそしむ光景がちらほらありましたが,そんなものは知らぬ.おなかが空いている.

台風が明日にでも命中するとの予報だったので,2次会も無く解散しホテルへと引き上げる


2日目はもはや台風の気配など無く普通に晴れ.
全日程にわたる雨を覚悟し,折り畳みではないかさを持ってきた私は拍子抜けしました.
自分の準備が出来ていなかったので,直前のセッション中に準備をして何とか間に合わせる.
どうもやっつけな準備をしたせいか発表もやっつけに終わってしまった気がするが,終わればいいのだ終われば.

自分の発表を終えるとイマイチ気が乗り切らず,ほどほどにセッションを見てほどほどに切り上げるとする.
東大の中央食堂に乗り込むとそれはそれは綺麗で広大で自分の学生時代との差に悲しみを背負うが,
悲しみも早々に忘れてデザートをパクついて糖分を補給,悠々と会場に戻って最後のセッションを見て今年のユーザー会もお開きに.

今年のノベルティは扇子,拡大鏡,メモ帳の三択でしたが,
扇子以外にもらうものが無かったので実質一択.メモ帳もらっても...ねえ...?
しかしSASユーザー総会,椅子をですね,もう少し増やして欲しいなと思いますね.頼むわ.

会の後は知人との飲み会が新橋であるので急行し,安居酒屋でチャーハン,ジンジャーエール,モツポン酢をいただく.
ジンジャーエールはもう少し辛口のほうが私の好みであったが,所詮安居酒屋.高望みはすまい.
知り合いも変わりなくSASとは全然関係ない話に華を咲かせ,23時頃にホテルに引き上げそのままお布団へ.

9月にPharmaSUGがあるので,次はそのタイミングで東京にいくとしましょうかね.
申請やら報告書がめんどくさいので有給で行きたいが,果たして申請が通るかどうか.
通って欲しいなあ.

2018年6月24日日曜日

SAS飲み会に行ってきたお話

なにやら大層なブログタイトルになってしまいましたが.
単にネタを持ち寄りましょう,という簡単な会に6月23日に行ってきました.

18:30開始,会場が梅田となってはそこそこ早めに会社を出ないと行けません.しかしプログラム書いてると意外と遅くなり,すわ遅刻か.となりながらも退勤.会場へ急ぐ.

道中急いでいたため軽食と飲み物のことをすっかり失念しており,はたと思い至ったときには時すでに遅し,あたりにコンビニはない.これは21:30まで飲まず食わずかと覚悟を決めるが,自販機を発見しお茶を買うことには成功する.軽食は無く空腹と共に人の発表を聞くか...とがっくり肩を落としながらも引き続き会場へと急ぐ.

当日の会場がわかりにくいビルで,ビルにたどり着くまでに同じ道を行ったり来たり.この筋のはずなんだがなぁとウロウロし,なんとかビルを見つけて(なんと2度もそのビルの前を素通りしていた)建物内に入る.
なんと建物の中でも迷子になってしまい,あちこちさ迷い歩く羽目に.壁に並ぶ不審者に注意の張り紙が大変心に響いた.

ほうほうの体で会場へin.全員揃うのをまち,本日のネタお披露目が始まる.

私は今回はsgplotのオプションをいくつか紹介すると言う,目新しさのかけらもない微妙なネタを持っていきました.

私のは手垢のついたネタなのでいいのですが,classdataオプション(proc means)やweightステートメントのzeroオプション(proc freq),果てはフレームワークのようなもの(プログラムの定型部分を自動で行うなど)などと出るわ出るわ盛り沢山ですね.

classdataオプション,存在は知っていましたが使ったことなかったですねー.いっつもダミーデータをmergeしてましたね...一度使ってみましょう.

次は忘年会のシーズンに集まる感じですかね.またその時はどんなネタが出てくるか,今から楽しみです.

最後になりましたが,来てくださった皆様,これなかった皆様,ありがとうございました.

2018年6月15日金曜日

call is8601_convertルーチンの話

is8601_convertコールルーチンは,日付や2つの日付の間の間隔をiso8601フォーマットに変換します.
この日付の間隔のことを,デュレーション値と呼ぶようです.
たとえば始点が4月1日,終点が4月2日のとき,デュレーション値は1日です.
第一引数に変換前の値に何が入るかを指定し,第二引数に変換後の値が何であるかを指定します.
第三引数に変換前の値,第四引数に変換に使う値,第五引数に変換後の値を格納する変数を指定します.

第一引数に'du/dt'と指定すると,第三引数にデュレーション値,第四引数に日付を指定します.
これが逆になるとエラーが返ります.
同様に'dt/du'と指定すると,日付→デュレーション値の順に指定をしないといけません.

とてつもなくわかりにくいので以下にいくつか具体例を示します.

data _null_ ;
  format hoge e8601dn. ;
  /*変換前の値はデュレーション値と日付,変換後に期間の開始日を出力,期間は8週間,期間の終点は2018年1月1日*/
  call is8601_convert('du/dt' , 'start' ,"p8w" ,"2018-01-01" ,  hoge ) ;
  putlog hoge = ;
run ;

hoge=2017-11-06


data _null_ ;
  format hoge e8601dn. ;
  /*変換前の値は日付とデュレーション値,変換後に期間の終了日を出力,期間の開始日は2018年1月30日,期間は8週間*/
  call is8601_convert('dt/du' , 'end' ,"2018-01-30" ,"p8w" , hoge ) ;
putlog hoge = ;
run ;

hoge=2018-03-27


data _null_ ;
  length hoge $20 ;
  format hoge $n861b. ;
  /*変換前の値に日付と日付,変換後はデュレーション値,期間の始点は1月30日,終点は3月14日*/
  call is8601_convert('dt/dt' , 'du' ,"2018-01-30T12:30" ,"2018-03-14T14:45" , hoge ) ;
  putlog hoge = ;
run ;

hoge=P1M15DT2H15M

2018年5月23日水曜日

ods rtfでのファイル名の上限の話

何やら身辺が立て込んで仕方がないです.

ods rtfではパス名とファイル名の合計文字数は180文字が上限のようですね.
windows上の制限は256文字なので,上限いっぱい使ってもまだ余裕があるみたいですが.

ただパス名込の文字数なので,76文字以上のフォルダに移動するとwindowsでの制限に引っかかりますね

180文字を超えた名前を指定すると,
無効なファイル名です,致命的なODSエラーが発生しました.のerrorが返ります.


2018年4月9日月曜日

putlogでlogに任意の値を出力する話

putlogをご存知でしょうか.logに任意の文字列を出力するステートメントです.

これを使ってデータの中のイレギュラーデータを(ある程度)見つけようと言う試みです.
ためしに下のデータセットがあるとして,変数aには欠測は本来入らないと仮定してください.
その変数aに欠測が入ってるかどうかと,欠測の時のほかの変数の値をlogに出力するのが以下のプログラムです.
私は最近知っていたく感動したのですが,普通にご存知の方が多いのでしょうかね...?

putlogで指定した変数の値が出力できるので,それを利用した形になります.
自分で指定した条件に該当するレコードがなかった場合,logには何も出てきません.

"E" "RROR"とputlogに指定しているのは,ログには"ERROR"と出したいが,
ログに出す条件のレコードが何もないときはERRORとログに残したくないためです.
涙ぐましいですね.この小細工が.

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

data hoge ;
 a = 1 ; b = 1 ; c = 6 ; output ;
 a = 1 ; b = 2 ; c = 5 ; output ;
 a = . ; b = 3 ; c = 4 ; output ;
 a = 1 ; b = 4 ; c = 3 ; output ;
 a = . ; b = 5 ; c = 2 ; output ;
 a = 1 ; b = 6 ; c = 1 ; output ;
run ;

data _null_ ;
  set hoge ;

 /* -- 想定の値以外(今回は欠測)があるときにlogに出力 -- */
  if a = . then do ;
    putlog "E" "RROR_A" ;
    putlog b = c = ;
  end ;

run ;


2018年4月4日水曜日

ルール間違ってた

「pandemic」というボードゲームにリベンジしてきました.
致命的なルールの勘違いをしていたことに気づいてしまいました.
4戦全敗した後でですが...
規則の確認を雰囲気でやってはいけない好例です.

このゲーム,ワクチンを作って病原菌を根絶(盤上から赤青黒黄の四色の駒を完全に取り除くこと)したら勝ちだと思っていたのですが,ワクチンを4種作った段階で勝ちなのですね.そら勝てないね.

ワクチンは一枚目の画像の版下部にある,なんか四角い四色の駒ですね.
一枚目の画像は赤のワクチン駒を獲得する前に時間切れで負けです.















2枚目は一応勝ちの状態です.盤下部に赤青黄黒の四角いワクチン駒を獲得できています.
しかしこれは参考程度です...手放しで喜べない勝ちなのです...
本来ランダムに決まるプレイヤの役職を,事前に選んでプレイしましたからね...
普通に勝利は次回に持ち越しですね...


2018年3月3日土曜日

グラフのラベルに改行を入れる話

SASでグラフを出す時に,軸のラベルや軸目盛の値の改行をするにはannotateをするしかない,と思っていました.
ところがどっこいaxisステートメントの指定の仕方で改行できるようです.
軸ラベルを改行させてほしいとリクエストが来るたびに,相手のことを呪っていたわけですが
これからはそんな呪いからも解放されますね.

sgplotで軸目盛のラベルの改行は出来たのですが,軸そのものへのラベルの改行のやり方が分からないのが次なる課題ですね...

/*---------- testデータ ----------*/

data hoge;
    call streaminit(11);
    do A=1 to 5;
        Y = int(rand('uniform') *10 );
        output;
    end;
run;

*---------- gplotでの場合 ;
proc gplot data = hoge;
    plot Y * A / vaxis = axis1 haxis = axis2 ;
    symbol1 c = black V = dot I = join;

    *---------- y軸目盛;
    axis1
      order = (0 to 5 by 1)
      origin=(10, 30)
      ;

    *---------- x軸目盛;
    axis2
      order=( 1 to 5 by 1 )
      origin=(10, 30)
      label=("x軸" j = c "ラベル")
      value=(tick = 1 "1" j = c "時間" j=c "経過"
             tick = 2 "2" j = c "時間" j=c "経過"
             tick = 3 "3" j = c "時間" j=c "経過"
             tick = 4 "4" j = c "時間" j=c "経過"
             tick = 5 "5" j = c "時間" j=c "経過"
            )
    ;
 run;

gplotでの出力イメージ












*---------- sgplotでの場合 ;
proc sgplot data = hoge  ;
    series x = A y = Y  ;

    xaxis
      type  = discrete
      fitpolicy=splitalways
      splitchar="*"
      label="x軸ラベル"
      values=(1 2 3 4 5)
      valuesdisplay=("1*時間*経過" "2*時間*経過" "3*時間*経過" "4*時間*経過" "5*時間*経過")
    ;   
run ;

sgplotでの出力イメージ

2018年2月19日月曜日

勝利を知りたい

この記事はsas関係ないです.ボードゲームです.

私趣味でボードゲームをたまにやるのですが,こないだぼこぼこに負けました.
「パンデミック」というプレイヤで協力して目的を達成するゲームがあるのですが,
そこでまあ見事に2戦2敗キメましたのでその供養です.

4人のプレイヤで協力して世界に蔓延る 4種類の病原菌を根絶することを目指すゲームです.詳しくは適当に調べてください.
私はそこで「ほかのプレイヤを動かせる」という固有の能力を持っていたのですが,あまり活用できませんでしたね…指定した離れた土地にプレイヤを一人ワープさせる力をもっと使いこなせていれば敗北とは違った結果が得られたかもしれません.

次やったときは確実に病原菌を根絶してくれよう覚悟してろ


敗北の盤面がこちら

2018年2月14日水曜日

割り算に欠測とか0が混ざってる時の話

割り算する時に,分子が0の時や分母が0とか欠測の時はどうなりますか,と聞かれました.

割った結果は何となくわかるのですが,logになんと出力されるのかが分からなかったので
実際に実行してlogを見たと言うだけの話です.

data DATA ;
    a = 1 ; b = 2 ; output ;
    a = 0 ; b = 2 ; output ;
    a = 1 ; b = 0 ; output ;
    a = . ; b = 2 ; output ;
    a = 1 ; b = . ; output ;
run ;

data DIVISION ;
    set DATA ;
    c = a / b ;
run ;


実行log





結果















分母が0の時は,0による割り算がありましたってnoteが出ていますね.
欠測がどっちかに入っている時は,欠測値を含んだ計算により...と言うよく見るやつが出ています.
その下の以下の箇所で演算式を計算できなかった...と言うのは0による割り算に対して出ています.
分子が0で分母に数が入っている場合は,結果が0として出ます.









2018年1月27日土曜日

SASオフ会に行ってきました

露骨にブログ記事を稼いでいくスタイルを見せつけています.
特に内容はないです.

先日(1月26日)にデータステップ100万回の管理人さん含めて6人でオフ会しました
せっかくの飲み会…ということで一人一つsasネタを持ち寄って軽く紹介しあいました
ユーザー会より気楽にネタ作ってポンと発表できるのはいい感じでしたね…

私が持って行ったネタの資料はswayに転がしておきました
こちらから閲覧できます.
pptの資料をswayに変換しただけなのでイマイチ見にくいかもしれません

次の開催は4月か5月頃ですかねえ.年度末は忙しいですし.

2018年1月25日木曜日

Wilcoxonの符号付き順位和検定_proc univariate

SASでWilcoxonの符号付き順位和検定は,proc univariateを実行することで出来ます.
こちらはWilcoxonの順位和検定と異なり,対応のあるデータに対して行います.

univariateプロシジャに,対応前後の差の値(変化量の様なもの)を格納した変数を指定します.
以下の仮のデータセットで言う所の「CHG」に当たります.(after-beforeで導出)

以下のプログラムではTEST_UNIのデータセットにunivariateプロシジャの結果を出力しています.
Wilcoxonの順位和検定のp-値は「_PROBS_」に格納されています.

今回はnoprintを指定することでoutput画面への出力を抑制していますが,
outputでは「位置の検定」の項目に結果が出されています.


/*---------- プログラム ----------*/

data TEST ;
 CAT = 1 ;  BEFORE = 1 ; AFTER = 3 ; CHG = 2  ; output ;
 CAT = 1 ;  BEFORE = 2 ; AFTER = 6 ; CHG = 4  ; output ;
 CAT = 1 ;  BEFORE = 3 ; AFTER = 1 ; CHG = -2 ; output ;
 CAT = 1 ;  BEFORE = 4 ; AFTER = 3 ; CHG = 1  ; output ;
 CAT = 2 ;  BEFORE = 5 ; AFTER = 2 ; CHG = -3 ; output ;
 CAT = 2 ;  BEFORE = 6 ; AFTER = 1 ; CHG = -5 ; output ;
 CAT = 2 ;  BEFORE = 7 ; AFTER = 6 ; CHG = -1 ; output ;
 CAT = 2 ;  BEFORE = 8 ; AFTER = 5 ; CHG = 3  ; output ;
run ;

proc univariate data = TEST outtable = TEST_UNI noprint ;
 class CAT ;
    var CHG  ;
run ;