LMEモデルフィッティング MATLAB編

線形混合効果（LME）モデル推定をMATLABで行う。

2013b以上じゃないとLMEモデルが扱えないので注意。（非線形混合効果（NLME）モデルなら2012からある）

LMEモデルについては過去記事参照のこと

統計初心者でも2分で雰囲気理解できるLME（線形混合）モデル

http://foslave.blogspot.jp/2014/12/2lme.html

LMEモデル推定 R言語編

http://foslave.blogspot.jp/2014/12/lme.html

---

以下、やり方

①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