固定効果とランダム効果を持つ混合モデルの推定をする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
-------------
0 件のコメント:
コメントを投稿