12月21日，SASの発表会に行ってきた．
年末だけあって何分立て込んでいたがなんとか間に合う時間に退社し，いそいそと会場に向かう．
当日はあいにくの雨，ただでさえ寒いのにさらに冷えるのはほんと勘弁してもらいたい．
途中にある松屋で牛飯をかきこんでから会場に入る．松屋は牛丼ではない，牛飯です．
会場に入って待つことしばし，まあ来るわくるわ参加者が．途中から何人来るかの把握をあきらめていたが，やはり今回は多かった．全部で10人ちょっとか．
とりあえず椅子だけ並べて会長の100万回の管理人さんの挨拶を頂いて，さて皆で発表．
1人5分程度のライトニングセッションなのだがコレが気楽で良い．30分質疑込みとか重すぎてしんどいだけである．
グラフから始まり文字コードなどのバラエティに富んだ発表がたくさんだった．よくもまあそんなネタが出てくるもんである．
終わった後は飲み屋で串を食べ食べ他愛も無い話に花が咲く．私はお酒はダメなのでミックスジュースだった
気づくと終電の時間が近く，食べ過ぎの体に鞭打って駅に走りとても苦しかった．たまたま電車が遅れていたので終電に間に合い無事帰宅と相成った．
毎度場所の選定から飲み会のセットまでやってくださる方には頭が上がらない．
この場を借りて感謝をお伝えします．ありがとうございます．
そろそろ人数が増えてきたので，どこか良い感じの場所があるといいなあ，どこかに親切なS○S社とかいませんかね．
連絡手段もメールはそろそろ限界ですしなあ
2019年1月24日木曜日
2019年1月1日火曜日
文字列の前後のスペースを除いてマクロ変数を作る話
たとえば" test"のような，マクロ変数にしたい文字列の前後に要らないスペースがある時，
スペースを取り除いてマクロ変数にすることがあると思います．
今回の例で言うと"test"と言う値をマクロ変数に格納する場合です．
そんな時はcall symputxが便利ですよ，と言う話です．
call symputxだと前後のスペースを取り除いてマクロ変数にしてくれます．
call symputx("a" , " test")と，call symput("a" , left(trim(" test") ))は同じです．
割と色んな場所で見かける紹介なのですが，あえてcall symputでマクロ変数を作っていてスペースを取り除き忘れたために
意図しない挙動を取ったプログラムを見かけたので私も紹介します．
/*---------- 以下マクロ変数作成時のプログラムの例 ----------*/
data _null_ ;
call symput("a1" , " test") ;
call symput("a2" , left(trim(" test") )) ;
call symputx("a3" , " test") ;
run ;
%put &=a1 ;
%put &=a2 ;
%put &=a3 ;
/*---------- 作成したマクロ変数をlogに表示 ----------*/
39
40 %put &=a1 ;
A1= test
41 %put &=a2 ;
A2=test
42 %put &=a3 ;
A3=test
スペースを取り除いてマクロ変数にすることがあると思います．
今回の例で言うと"test"と言う値をマクロ変数に格納する場合です．
そんな時はcall symputxが便利ですよ，と言う話です．
call symputxだと前後のスペースを取り除いてマクロ変数にしてくれます．
call symputx("a" , " test")と，call symput("a" , left(trim(" test") ))は同じです．
割と色んな場所で見かける紹介なのですが，あえてcall symputでマクロ変数を作っていてスペースを取り除き忘れたために
意図しない挙動を取ったプログラムを見かけたので私も紹介します．
/*---------- 以下マクロ変数作成時のプログラムの例 ----------*/
data _null_ ;
call symput("a1" , " test") ;
call symput("a2" , left(trim(" test") )) ;
call symputx("a3" , " test") ;
run ;
%put &=a1 ;
%put &=a2 ;
%put &=a3 ;
/*---------- 作成したマクロ変数をlogに表示 ----------*/
39
40 %put &=a1 ;
A1= test
41 %put &=a2 ;
A2=test
42 %put &=a3 ;
A3=test
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-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の値に差があると言えそうです．
正直指定できるオプションやらステートメントが多いので，リファレンス見るのが一番ですかね．
ざっくり言うと平均からのズレを見たい時に使うと思います．
実行は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 ;
どんな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だった．結局もったいなくて使えずそのまま放置するんだろうなあ．
しかしよく見ると受付が9:30～10:00，関西から行くにはいささか時間が早いので夜行バスで向かうとする．
結果的にここで夜行バスを選択したのが吉と出た．
迎えた当日，台風「こんにちはー」
SASユーザー総会の時も台風と一緒に東京に行った記憶がよみがえる．
すわ電車運行見送りかと震えるも私は夜行バス．前日に東京に入っていたので大した影響は受けなかった．
建物に入ればこちらのもの．恙無くSDEが始まった．
プロジェクタが何を思ったか赤を青に，青を赤に表示している所為でSASのアイコンが赤い．違和感違和感．
中国の規制当局の話が一番興味深かったが残念英語がイマイチわからない．
うーんこのと聞きつつ，英語の聞き取りは練習したほうがいいなあと実感した．
後はまあ人工な知能をpythonで実装した話の事例である．実に最近のブーム．
RでPPTへの出力，との内容も面白そうではあるがバリデーションの壁は厚く高そう．
うちの大将のポスターセッションも見てきたがやはり大将は話がうまい．
あれくらい話せるようになればなあと思うも，練習あるのみか．
いや人前で話すのは苦手なんだが
終わる頃には風が荒れる空模様だったが雨は無く，後楽園で知り合いと飲んで事故も無くホテルに引っ込む．
どうせ台風で交通機関荒れるからと一泊して本当に良かったと思いつつ12:00頃に就寝．
風が轟くもなんのその，すぐに夢の世界へとさようならし，翌日の昼過ぎに東京に別れを告げた．
自宅の隣家のブロック塀が見事に崩落している様に絶句するも，家は無事で何よりであった．
次は来年のSASユーザー総会まで何も無いかなあ，しかし東京に行くのはお金がかかって仕方が無い．
今回のSASグッズはSASのリングNOTEだった．結局もったいなくて使えずそのまま放置するんだろうなあ．
登録:
投稿 (Atom)