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の値に差があると言えそうです.

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