LMEモデルフィッティング R言語編

続き。

LMEモデル推定をR言語を用いて実践してみる。

まずRをインストールし、RGuiを開く。
今回使うパッケージはlme4とlmeTestなので、これらを読み込む。lme4のlmer()関数によってLMEモデル推定を行うが、lme4だけだと固定効果のp値が得られないため、lmeTestを追加で読み込む必要がある。

>Install.packages(“lme4”) #lme4インストール
>library(“lme4”)　#lme4読み込み
>options(repos="http://cran.ism.ac.jp/")
>install.packages("lmerTest")　#lmerTestインストール
>library(“lmerTest”) #lmerTest読み込み

次にデータを読み込んでみる。今回、4人の学生の勉強日数と点数のテストデータをCSVファイルで用意した。

> s_data <- read.csv("score_data.csv")
> s_data
   subject days    Score
1        0   10 97.88387
2        0    8 98.59650
3        0    6 97.65732
4        0    4 96.16080
5        0   10 99.59362
6        0    8 97.77845
7        0   10 98.22498
8        0   10 98.04627
9        0    8 98.25393
10       0   10 97.46261
11       0    6 98.18022
12       0    4 96.99735
13       0   13 97.90047
14       0    6 98.66602
15       0   15 98.52445
16       0   12 97.65170
17       0    8 97.65037
18       0   14 98.28197
19       0   11 97.47814
20       0   10 98.51969
21       0   11 98.76532
22       0    8 97.27131
23       0    6 98.02530
24       0   10 98.52534
25       0   14 99.38336
26       0    4 94.42770
27       0    7 98.01417
28       0   12 98.10540
29       0    4 97.87208
30       1   10 97.82623
31       1   10 96.94224
32       1   14 98.31716
33       1   13 95.45953
34       1   18 96.94513
35       1   11 93.20422
36       1    6 94.51405
37       1    8 91.38044
38       1   16 93.66502
39       1    6 93.81745
40       1   13 97.06315
41       1    8 96.63337
42       1   11 96.66075
43       1    9 96.37080
44       1   13 98.71230
45       1   18 97.57873
46       1    5 94.33944
47       1    7 94.78274
48       1   17 97.75903
49       1   11 96.82348
50       1   17 97.77848
51       1    5 91.73288
52       2   15 90.96240
53       2   12 91.43013
54       2   15 96.27122
55       2   14 93.39732
56       2   18 92.99953
57       2   11 93.39127
58       2    9 91.33970
59       2    5 87.46454
60       2   18 91.80603
61       2    4 90.71793
62       2   10 91.20598
63       2    4 90.24194
64       3    5 90.67530
65       3   17 95.64061
66       3   14 97.77640
67       3   15 96.77373
68       3   14 95.39185
69       3   14 97.48028
70       3    5 91.10883
71       3    4 88.89533
72       3   10 95.65205
73       3    5 92.51790
74       3    6 89.78949
75       3   14 94.84116
76       3   15 96.99177
77       3   13 95.21246

グラフにplotすると

>plot_data <- function(x){
for (i in 1:nrow(x)){
plot(x[i,2],x[i,3],col=x[i,1]+1,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),ann=F)
par(new=T)
}
}
>plot_data(s_data)

各学生の勉強時間（横軸）とテストの点数（縦軸）

さて、ここからが本題であるが、lmer関数でLMEモデル推定を行う。

第一引数にモデル関数、第二引数に解析データの行列を入れる。

 モデル関数の表記は以下のように行う。Yは従属変数ラベル、Xi(i=1,2,…)は説明変数ラベル（切片は省略するが、切片のみを記述する場合は１を入れる）、Ci(i=1,2,…)はクラスタラベルである

Y~X1+X2+X3…+(X1'+X2'+…|C1)+(X1"+X2"+…|C2)+…

例えば、上記のデータで、勉強日数に応じて点数が変化するようなモデルを考える。
勉強日数の係数（傾き）が学生共通で一定（固定効果）で、切片は学生それぞれに対して変わる（変量効果）ような場合、

Score~days+(1|subject)

となる。従って、

lmer(Score~days+(1|subject),data=s_data)

で推定が行える。結果を見るには、

> test_model1 <- lmer(Score~days+(1|subject),data=s_data)
> summary(test_model1)

とすれば見ることができる。
推定結果の見方を説明する。

Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [merModLmerTest]
Formula: Score ~ days + (1 | subject)
   Data: s_data

REML criterion at convergence: 297.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4714 -0.6311  0.1049  0.7079  2.1640 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 8.186    2.861   
 Residual             2.232    1.494   
Number of obs: 77, groups:  subject, 4

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 91.59004    1.51029  3.61000  60.644 1.48e-06 ***
days         0.31828    0.04266 72.11000   7.461 1.53e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
days -0.298

赤文字にした部分が変量効果（Random effects）、固定効果（Fixed effects）についての記述。
変量効果については、切片（Intercept）の変量効果における正規分布の分散(Variance)と標準偏差(Std.Dev)が記してある。
固定効果については、dayの係数と切片について、推定した値（Estimate）とP値（Pr）が記してある。
固定効果についてはfixef()で、それぞれのクラスタにおける変量効果についてはranef()で参照できる。

> fixef(test_model1)
(Intercept)        days 
  91.590043    0.318283 
> ranef(test_model1)
$subject
  (Intercept)
0   3.4623447
1   0.6747299
2  -3.3261432
3  -0.8109314

傾きにも変量効果を加えて再度推定してみる。

> test_model2 <- lmer(Score~days+(days|subject),data=s_data)

plotしてみる

>lmer_plot <- function(x,y){
rmat <- ranef(y)
fmat <- fixef(y)
for (i in 0:max(x[,1])+1){
func_buff <- function(x){(fmat[2]+rmat[[1]][i,2])*x+(fmat[1]+rmat[[1]][i,1])}
plot(func_buff,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),col=i,ann=F)
par(new=T)
}
fix_func <- function(x){fmat[2]*x+fmat[1]}
plot(fix_func,xlim=c(min(x[,2]),max(x[,2])),ylim=c(min(x[,3]),max(x[,3])),lty=2,ann=F)
par(new=T)
}

>lmer_plot(s_data,test_model2)

点線：固定効果モデル
各実線：それぞれの学生の混合効果モデル

また、切片、傾きそれぞれの変量効果の確率分布をplotすると次のようになる。

切片（定数項の係数）の確率分布

傾き（daysの係数）の確率分布