線形混合効果(LME)モデル推定をMATLABで行う。
2013b以上じゃないとLMEモデルが扱えないので注意。(非線形混合効果(NLME)モデルなら2012からある)
LMEモデルについては過去記事参照のこと
以下、やり方
①CSVファイルをテーブルでインポートする。
②Statistics Toolboxのfitlme関数でLMEモデル近似を行う。
>> lme1 = fitlme(allchangeblog,'change_score_ratio~days+(days|subject)')
fitlmeの第一引数は先ほどインポートしたデータテーブル名、第二引数がモデル関数である。
モデル関数の書き方はRのlmerと全く同じで、
従属変数 ~ 説明変数 + (説明変数|グループ変数) + …
である。
切片は省略されるので、もし切片のみを変量効果に加えるなら、1を入力する。
従属変数 ~ 説明変数 + (1|グループ変数) + …
出力は以下のようになる。AIC(赤塚情報量基準)や対数尤度等も出してくれる。
選択した線形混合効果モデルの近似方法: ML
モデル情報:
観測の数 309
固定効果の係数 2
変量効果の係数 20
共分散パラメーター 4
公式:
Linear Mixed Formula with 2 predictors.
モデルの適合性に関する統計量:
AIC BIC LogLikelihood Deviance
-1384.5 -1362.1 698.27 -1396.5
固定効果の係数 (95% 信頼区間):
Name Estimate SE
'(Intercept)' 0.10148 0.01004
'days' -0.003319 0.00039318
tStat DF pValue Lower
10.108 307 6.3802e-21 0.081726
-8.4414 307 1.2587e-15 -0.0040927
Upper
0.12124
-0.0025453
変量効果の共分散パラメーター (95% 信頼区間):
グループ: subject (10 水準)
Name1 Name2 Type
'(Intercept)' '(Intercept)' 'std'
'days' '(Intercept)' 'corr'
'days' 'radius' 'std'
Estimate Lower Upper
0.029357 0.01751 0.04922
-0.85901 -0.99775 0.67225
0.00069051 0.00016489 0.0028916
グループ: 誤差
Name Estimate Lower Upper
'Res Std' 0.023834 0.021973 0.025851
③固定効果と変量効果の値を取得する。
変量効果はグループごとに順番に出力され、一列のベクトルで返されるので注意。
>> F = fixedEffects(lme1) %固定効果の取得
F =
0.1015
-0.0033
>> [R,Rnames] = randomEffects(lme2) %変量効果の取得
R =
-0.0599
0.0012
-0.0271
0.0005
0.0296
-0.0004
-0.0003
-0.0003
0.0224
-0.0004
0.0140
-0.0003
-0.0163
0.0002
-0.0002
0.0003
0.0441
-0.0009
-0.0064
0.0001
Rnames =
Group Level Name
_________ _____ _____________
'subject' '1' '(Intercept)'
'subject' '1' 'radius'
'subject' '2' '(Intercept)'
'subject' '2' 'radius'
'subject' '3' '(Intercept)'
'subject' '3' 'radius'
'subject' '4' '(Intercept)'
'subject' '4' 'radius'
'subject' '5' '(Intercept)'
'subject' '5' 'radius'
'subject' '6' '(Intercept)'
'subject' '6' 'radius'
'subject' '7' '(Intercept)'
'subject' '7' 'radius'
'subject' '8' '(Intercept)'
'subject' '8' 'radius'
'subject' '9' '(Intercept)'
'subject' '9' 'radius'
'subject' '10' '(Intercept)'
'subject' '10' 'radius'
④推定した関数をプロットする。
>> h = gscatter(days,change_score_ratio,subject) %データをグループごとにプロット
>> hold on %グラフの重ね合わせ設定
>> colors = get(h,'Color'); %プロットした色情報を取得
>> for I = 1:10 %グループごとにループを回す
>> fitted_model = @(t)((F(2)+R(I*2))*t+F(1)+R(I*2-1)) %関数を作成
>> plot((1:20),fitted_model(1:20),'Color',colors{I}) %プロット
>> end