 図1 ランダム切片モデル |
 図2 ランダム切片・傾きモデル |
<< 目次へ
最終更新日 ： 2008.11.19
| Contents
|
経時測定データとは1つの個体から時間の経過とともに複数回にわたって"繰り返し"記録されたデータのことである。ある個体から"繰り返し"データをとったという意味から、経時測定データは繰り返し測定データともいうことがある。
なお、似たような用語で反復測定データというものがあるが、反復測定と繰り返し測定は全く異なるものなので注意が必要である。これは専門家でさえも混同して使っている場合が多いので、参考書を読むときには常に注意して読み勧める必要がある。
反復測定と繰り返し測定を区別するためには、「サンプルサイズはいくつか」ということを自問自答することが重要である。例えば、表1は反復測定のデータ形式を表しており、表2は繰り返し測定のデータ形式を表している。
| A | B | C |
| S1 S2 S3 S4 S5 |
S6 S7 S8 S9 S10 |
S11 S12 S13 S14 S15 |
| A | B | C |
| S1 S2 S3 S4 S5 |
S1 S2 S3 S4 S5 |
S1 S2 S3 S4 S5 |
表中のSは被験者を表しており、Sの添字として番号が振られているが、これは被験者番号である。表1においては、Aという条件の下で5人の異なる被験者が用いられている。同じようにBという条件下でも、Cという条件下でもそれぞれ異なった5人の被験者が割り当てられている。結局、表1におけるサンプルサイズは15である。
これに対して表2はどうだろうか。表2では1人の被験者が3つの異なる条件下において、繰り返してデータが得られるように配置されている。例えば、S1は条件Aだけでなく、条件B、条件Cにおいてもデータが取られている。ここでのサンプルサイズは5である。
このようにその実験に用いられた被験者の人数、すなわちサンプルサイズの大きさを考えれば明らかである。反復測定とは「ある異なる条件下において、条件間で異なる個体から反復してデータを得ること」であり、繰り返し測定とは「ある異なる条件下において、条件間で同じ個体から繰り返してデータを得ること」である。
つまり、極端な話、ある実験において1つの個体から2つ以上のデータがとられた場合、それはもはや独立なデータではないので、擬似反復を疑う必要がある。そして、この擬似反復は容易に見つけることが可能である。これについては、後ほど解析例を示すとともに紹介しよう。
さて、前述した表2のような繰り返し測定データを解析するためには線形混合効果モデルを適用することになる。もし表1のように、独立であるデータである場合は今まで行ってきたような一般線形モデル、ないしは一般化線形モデルで扱ってきたようなモデルを用いればよい。しかしながら、表2のような、いわゆる対応ありのデータを扱う場合は少し工夫が必要となるのである。線形混合効果モデルとは、そういった"少しの工夫"がつけ加えられたモデル解析である。
表1のようなデータセットにおいて、例えば従属変数をY、それを説明する独立変数としてGroupを用いて分散分析モデルとして解析することを考える。このモデル式は次のように表せる。
Y = μ + β*Group + ε ・・・ 式1
μは平均、βはGroupがYを説明する大きさ、εは誤差である。
線形混合効果モデルにおいては、次式のようなモデル式の母数を推定することになる。
Y = μ + β*Group + δ + ε ・・・ 式2
ここに新たなδというパラメータ（母数）が出てきたが、これは個体のランダム効果と呼ばれるものである。簡単に説明してみると、ある個体はA、B、Cと異なる条件下において何らかのデータを測定されるわけだが、個体から得られる測定値には"ランダムなバラつき"をもつと考えられる。例えば、あなたが「長さ10cmの線分をこの紙に10本書いてみてください」といわれて、それを試したときのことを想像してみればよい。おそらく、すべての線分が完全に10cmピッタリに書かれていることはないだろう。9.5cmのものもあれば、11cmくらいのものもあるかもしれないのである。
このように個体から繰り返して測定値を得たときには必ずバラつきが現われる。そのようなことから、これをランダム効果とか変量効果と呼ばれたりする。つまり、線形混合効果モデルとはランダム効果を考慮したモデルであるといえる（より正確にいえば、δはランダム誤差である）。
ここで、とりあえずμ、β、δ、εの4つのパラメータが推定されたとして、上記モデル式に推定値を代入してみよう。
推定されたパラメータがそれぞれμ = 2, β = 3, δ = 0.3, ε= 0.5だとすると、
Y = 2 + 3*Group + 0.3 + 0.5 ・・・ 式3
となる。定数の部分（μ, δ, ε）は1つにまとめられるので、結局、
Y = (2 + 0.3 + 0.5) + 3*Group = 2.8 + 3*Group ・・・ 式4
となる。解析に用いたソフトウェアが出力するのはこの形であろう。すなわち、1つの定数項（切片）と独立変数の係数である。
そして、このように切片にランダム効果を組み込んだモデルのことを特にランダム切片モデルと呼ぶ。さらに、このランダム切片モデルを拡張したモデルとして、係数にもランダム効果を組み込んだモデルとしてランダム切片・傾きモデルと呼ばれるものがある。これは次式によって表される。
Y = μ + (β + λ)*Group + δ + ε ・・・ 式5
先ほど紹介した式2（ランダム切片モデル）に加えて、新たにλが組み込まれているのが分かるだろう。これはGroupの係数にランダム効果を組み込んだものであり、予測する回帰直線（曲線）の傾きがランダムに変化することを意味している。
まとめると、ランダム切片モデルは切片のみがランダムに変化するモデル、ランダム切片・係数モデルは切片と傾きがランダムに変化するモデルであるといえる。少し数式が出てきたものの、要はこの2点だけ理解しておけばよい。より理解しやすいように両モデルを表現している図を載せておこう（図1と図2を参照のこと）。
図1 ランダム切片モデル |
図2 ランダム切片・傾きモデル |
ところで、今まであえて伏せていたのだが、これらのモデルはいずれも「個体のデータごとに異なる回帰直線を当てはめる」というものなのである。とはいっても、線形混合効果モデルではランダム効果を考慮してパラメータを推定しているので、通常の回帰直線の当てはめとは異なる直線の方程式が得られることになる（個々のデータセットについてふつうの回帰分析をしてもダメということ）。
つまり、個体の数が5だとしたら、異なる5つの直線の方程式が得られることになる。表1のデータセットの解析においては、分散分析モデルが用いられたが、これは全てのデータについて（当てはめではないが）1つの直線の方程式が得られた。これに対して表2のデータセットの解析においては、先ほど述べたように5つの直線の方程式が得られるという点でまた大きく異なる。
ここでは実際に仮想データを用いて線形混合効果モデルを扱ってみよう。解析に使用するデータセットとして、文献[4]p2-表1.1とp161-表10.1のデータセットを使わせていただいた。ただし、これから提示するRの出力結果は本文に掲載されているものとは異なることに注意してもらいたい。
余談だが、使用するソフトウェアが異なれば、そのソフトウェアのプログラマーが採用しているアルゴリズムも異なる場合がある。そのため同じ理論に基づいた手法であっても、結果に若干の違いが生じることはしばしば起こりうる。したがって、なぜ異なる出力結果が得られるのかを理解したければ、それぞれのアルゴリズムを理解するしかないといえる。
話を元に戻すが、Rには線形混合効果モデルを扱うための関数が標準（インストール時）では実装されていない。そこでnlmeというパッケージをダウンロードする必要がある。そのnlmeというパッケージの中にlme()という関数が含まれているので、それを使うことにする。似たような関数がいくつかあるので、興味のある人は探してみるとよいだろう。
# nlmeパッケージのダウンロード
install.packages("nlme")
# nlmeの呼び出し
library(nlme)
使用するデータセット： [ p2-表1.1 table11.xls ] [ p161-表10.1 table101.xls ]
# table11.xlsのデータをベクトルで用意する
dat <- c(
117, 119, 122,
134, 123, 122,
133, 126, 117,
133, 128, 120,
126, 117, 115,
117, 114, 108)
# 被験者番号を表すカテゴリカル型変数
subject <- rep(1:6, rep(3, 6))
subject <- as.factor(subject)
# 経過時間を表すカテゴリカル型変数
time <- rep(1:3, 6)
time <- factor(time, levels=c(1, 2, 3), labels=c("before", "1h", "3h"))
# データフレームの作成
mydata <- data.frame(Subject = subject, Time = time, BP = dat)
# まずはdatを行列に変換し、t()を用いて転置行列にする
py <- matrix(dat, ncol=3, byrow=T)
py <- t(py)
# データプロット： matplot()を使用する
matplot(1:3, py, type="l", xlab="経過時間", ylab="血圧", ylim=c(100, 150), xaxt="n")
axis(1, at=1:3, label=c("before", "1h", "3h"))
図3 経過時間と血圧の変化
> model.eins <- lme(BP ~ Time, random=~ 1|Subject, method="ML", data=mydata)
> summary(model.eins)
Linear mixed-effects model fit by maximum likelihood
Data: mydata
AIC BIC logLik
118.2460 122.6978 -54.12299
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 4.519505 3.680529
Fixed effects: BP ~ Time
Value Std.Error DF t-value p-value
(Intercept) 126.66667 2.606616 10 48.59430 0.0000
Time1h -5.50000 2.327771 10 -2.36278 0.0398
Time3h -9.33333 2.327771 10 -4.00956 0.0025
Correlation:
(Intr) Time1h
Time1h -0.447
Time3h -0.447 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.09487990 -0.57828439 -0.03176163 0.53722903 1.79948751
Number of Observations: 18
Number of Groups: 6
ここでlme()の引数について少し解説しておく。モデルの指定（BP ~ Time）はlm()やglm()と全く同じである。異なるのはランダム効果の指定（random=~ 1|Subject）の部分だけである。これはSubjectをランダム効果として組み込んだランダム切片モデルである。つまりこのモデルでは切片のみが異なり、傾きはすべて同じであると仮定されたモデルである。したがって平行ないくつかの直線が引かれるグラフが出来上がる（図4）。
図4 ランダム切片モデルによる予測
図4を描くためのコードは次のようなものになる。若干、U字型に見えるのは横軸（elapsed time）がカテゴリカル型であるためである。横軸が連続型で、もっと細かく値を取るようにすれば完全な直線が引かれることになる（だが今回の横軸の値は3つである）。
# predict()を利用して予測値を行列として格納する
# 描き方は図3の場合と同じである
yosoku <- matrix(predict(model.eins), ncol=3, byrow=T)
matplot(1:3, t(yosoku), type="l", xlab="elapsed time", ylab="blood pressure", ylim=c(100, 150), xaxt="n")
axis(1, at=1:3, label=c("before", "1h", "3h"))
さて、図4の予測値のプロットを見てみると経過時間に伴って血圧が低下していく様子をそれなりに上手く説明しているように見える。しかし、ランダム切片・傾きモデルを当てはめた場合ではどうなるだろうか。試しにランダム切片・傾きモデルとして解析をしてみよう。
> model.zwei <- lme(BP ~ Time, random=~ Time|Subject, method="ML", data=mydata)
> summary(model.zwei)
Linear mixed-effects model fit by maximum likelihood
Data: mydata
AIC BIC logLik
122.9188 131.8226 -51.45942
Random effects:
Formula: ~Time | Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 7.106491 (Intr) Time1h
Time1h 3.436417 -0.837
Time3h 6.280722 -0.774 0.916
Residual 1.747551
Fixed effects: BP ~ Time
Value Std.Error DF t-value p-value
(Intercept) 126.66667 3.272801 10 38.70283 0.0000
Time1h -5.50000 1.892978 10 -2.90547 0.0157
Time3h -9.33333 3.018455 10 -3.09209 0.0114
Correlation:
(Intr) Time1h
Time1h -0.759
Time3h -0.761 0.799
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.98391995 -0.37464131 0.04811629 0.53633945 0.95774013
Number of Observations: 18
Number of Groups: 6
出力結果の詳しい解説は後ほどするとして、とりあえずランダム切片・傾きモデルによる予測値のプロットをしてみよう。そのためのコードは次のように書ける。
# predict()の引数にmodel.zweiとしただけである
yosoku <- matrix(predict(model.zwei), ncol=3, byrow=T)
matplot(1:3, t(yosoku), type="l", xlab="elapsed time", ylab="blood pressure", ylim=c(100, 150), xaxt="n")
axis(1, at=1:3, label=c("before", "1h", "3h"))
図5 ランダム切片・傾きモデルによる予測
図3（原データのプロット）と見比べてみると、どうやら図5（ランダム切片・傾きモデル）の方がより上手く経過時間による血圧の変化を上手く説明しているようである。なぜなら、ある1名の被験者（黒い直線）は時間の経過に伴って血圧が上がっているからである。これを図4（ランダム切片モデル）では上手く説明されていない。つまり、血圧の上昇している被験者のデータを無視していることになるわけだ。
このように予測値のデータプロットを眺めるだけでもどちらが良いモデルであるか見当をつけることはできる。しかし、そうはいってもあくまでも主観（分析者の経験からくる判断力）に大きく頼ることになるので、どちらがより良いモデルであるかを表す指標が欲しいところであろう。
結果からいって、lme()の出力結果の最初の項にモデルの適合度を表す指標が出力されているので、それを参考により良いモデルを選べばよい。あるいは、anova()を利用して両モデルの比較を行うことができる。その場合は尤度比の検定が行われることになる。
> anova(model.eins, model.zwei)
Model df AIC BIC logLik Test L.Ratio p-value
model.eins 1 5 118.2460 122.6978 -54.12299
model.zwei 2 10 122.9188 131.8226 -51.45942 1 vs 2 5.327142 0.3773
AIC（赤池情報基準）、BIC（ベイズ情報基準）の2つはいずれも数値が小さければ小さいほど良いモデルであることを示すものである。だから、これによればmodel.eins（ランダム切片モデル）の方が良いモデルであると判断される。またlogLKとは対数尤度（"たいすうゆうど"）のことであり、これを利用して尤度比検定が行われる。
2つのモデルのlogLKの差を2倍した値（{-51.45942 - (- 54.12299)} * 2 = 5.327142）は、2つのモデルの自由度の差（10 - 5 = 5）のカイ自乗分布に従うことになっている。これをpchisq()で試してみると、
> 1 - pchisq(5.327142, 5) [1] 0.3772727
となり、これはanova()によって返されたp-valueと一致していることが確認できる。ちなみに、私はこの尤度比検定についてはよく知らない。理屈はよく知らないが、p-value < 0.05の場合はModel2の方が良いモデル、p-value > 0.05のときはModel1の方が良いモデルだといえるようである。
[問題] 「図5（ランダム切片・傾きモデル）の方がより上手く経過時間による血圧の変化を上手く説明しているようである」と前述したが、AICの比較ではランダム切片モデルの方がより良いモデルだと判断できる。これについて考察せよ（ヒントは1人の被験者だけ血圧が時間とともに上昇していったという事実にある）。
さて、少し回り道をしたがlme()の出力結果について詳しく解説していこう。
# モデルの適合度
Data: mydata
AIC BIC logLik
118.2460 122.6978 -54.12299
これはすでに説明したとおりである。
# ランダム効果
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 4.519505 3.680529
ランダム効果： これはランダム切片モデルとして、Subjectがランダム効果として切片に影響を与えている。つまり、被験者ごとに異なる切片をもつことになる。
InterceptのStdDev（分散）は被験者間のばらつきを表している。ここで改めてランダム切片モデルを表す式2を見てみよう。Interceptの4.52というのはδの分散のことである。そしてResidualの3.68はεの分散である。
Y = μ + β*Group + δ + ε ・・・ 式2
# 固定効果
Fixed effects: BP ~ Time
Value Std.Error DF t-value p-value
(Intercept) 126.66667 2.606616 10 48.59430 0.0000
Time1h -5.50000 2.327771 10 -2.36278 0.0398
Time3h -9.33333 2.327771 10 -4.00956 0.0025
固定効果： ランダム効果に対してこれは固定効果と呼ばれるものである。固定効果の部分の見方はlm()やglm()で得られる結果と同じである。
Correlation:
(Intr) Time1h
Time1h -0.759
Time3h -0.761 0.799
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.98391995 -0.37464131 0.04811629 0.53633945 0.95774013
Number of Observations: 18
Number of Groups: 6
この部分は特に気にしてみる必要はないだろう。あえて説明するなら、Number of Observations: 18という部分は観測値の数（サンプルサイズではない）、Number of Groups: 6はサンプルサイズである（Groupと書いてあるが）。
以上でmodel.einsについての説明は終わりにするが、つけ加えてmodel.zweiのランダム効果の部分について解説しておこう。
Random effects:
Formula: ~Time | Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 7.106491 (Intr) Time1h
Time1h 3.436417 -0.837
Time3h 6.280722 -0.774 0.916
Residual 1.747551
今度は切片（Subject）のランダム効果に加えて、傾き（Time）にもランダム効果を組み込んでいる。Intercept = 7.12は傾きの分散、Time1hとTime3hの分散はそれぞれ3.44、6.28である。そしてResidual = 1.75は誤差（Residualは残差と訳すが）の分散である。
改めてランダム切片・傾きモデルのモデル式（式5）を見てみよう。式2（ランダム切片モデル）の場合と異なるのは、λの分散があるかどうかという点のみである。λの分散はTime1h = 3.44、Time3h = 6.28である。beforeは分散が0である（つまりベースラインとしている）ので、Time1hとTime3hはbeforeからの偏差であるといえる。
だから例えば、Time1hはbeforeに比べて3.44だけ大きなバラつきをもつということである。同じようにTime3hはbeforeに比べて6.28だけ大きなバラつきをもつということである。
Y = μ + (β + λ)*Group + δ + ε ・・・ 式5
続いてp161-表10.1 table101.xlsのデータセットを解析してみよう。
# とりあえず観測値をdatにベクトルとして格納しておく
dat <- c(
4.0, 4.0, 4.1, 3.6, 3.6, 3.8, 3.1,
4.2, 4.3, 3.7, 3.7, 4.8, 5.0, 5.2,
4.3, 4.2, 4.3, 4.3, 4.5, 5.8, 5.4,
4.2, 4.4, 4.6, 4.9, 5.3, 5.6, 4.9,
4.6, 4.4, 5.3, 5.6, 5.9, 5.9, 5.3,
3.1, 3.6, 4.9, 5.2, 5.3, 4.2, 4.1,
3.7, 3.9, 3.9, 4.8, 5.2, 5.4, 4.2,
4.3, 4.2, 4.4, 5.2, 5.6, 5.4, 4.7,
4.6, 4.6, 4.4, 4.6, 5.4, 5.9, 5.6,
3.4, 3.4, 3.5, 3.1, 3.1, 3.7, 3.3,
3.0, 3.2, 3.0, 3.0, 3.1, 3.2, 3.1,
3.0, 3.1, 3.2, 3.0, 3.3, 3.0, 3.0,
3.1, 3.2, 3.2, 3.2, 3.3, 3.1, 3.1,
3.8, 3.9, 4.0, 2.9, 3.5, 3.5, 3.4,
3.0, 3.6, 3.2, 3.1, 3.0, 3.0, 3.0,
3.3, 3.3, 3.3, 3.4, 3.6, 3.1, 3.1,
4.2, 4.0, 4.2, 4.1, 4.2, 4.0, 4.0,
4.1, 4.2, 4.3, 4.3, 4.2, 4.0, 4.2,
4.5, 4.4, 4.3, 4.5, 5.3, 4.4, 4.4,
3.2, 3.3, 3.8, 3.8, 4.4, 4.2, 3.7,
3.3, 3.4, 3.4, 3.7, 3.7, 3.6, 3.7,
3.1, 3.3, 3.2, 3.1, 3.2, 3.1, 3.1,
3.6, 3.4, 3.5, 4.6, 4.9, 5.2, 4.4,
4.5, 4.5, 5.4, 5.7, 4.9, 4.0, 4.0,
3.7, 4.0, 4.4, 4.2, 4.6, 4.8, 5.4,
3.5, 3.9, 5.8, 5.4, 4.6, 5.3, 5.6,
3.9, 4.0, 4.1, 5.0, 5.4, 4.4, 3.9,
3.1, 3.5, 3.5, 3.2, 3.0, 3.0, 3.2,
3.3, 3.2, 3.6, 3.7, 3.7, 4.2, 4.4,
3.5, 3.9, 4.7, 4.3, 3.9, 3.4, 3.5,
3.4, 3.4, 3.5, 3.3, 3.4, 3.2, 3.4,
3.7, 3.8, 4.2, 4.3, 3.6, 3.8, 3.7,
4.0, 4.6, 4.8, 4.9, 5.4, 5.6, 4.8,
4.2, 3.9, 4.5, 4.7, 3.9, 3.8, 3.7,
4.1, 4.1, 3.7, 4.0, 4.1, 4.6, 4.7,
3.5, 3.6, 3.6, 4.2, 4.8, 4.9, 5.0)
# 群を表すカテゴリカル型
group <- rep(c("T0", "T1", "T2", "T3"), c(63, 70, 56, 63))
group <- factor(group, levels=c("T0", "T1", "T2", "T3"))
# 経過時間を表す連続型
time <- c(1, 3, 5, 7, 9, 11, 13)
time <- rep(time, 36)
# 被験者番号を表すカテゴリカル型
subject <- rep(1:36, rep(7, 36))
subject <- as.factor(subject)
# データフレームの作成
mydata2 <- data.frame(GROUP = group, SUBJECT = subject, TIME = time, KALIUM = dat )
データフレームが作成できたら、まずは原データのプロットを行う。図6は群ごとに原データをプロットしたものである。1つのグラフィックウィンドウに複数の図を書き込むにはpar()関数を用いればよい。引数のmfrow=c(2, 2)というのは、2*2に分割せよという命令である。
# 原データを行列の形に変換 py <- matrix(dat, ncol=7, byrow=T) py <- t(py) # それを転置行列に変換 # グラフィックウィンドウの分割 par(mfrow=c(2, 2)) # 群T0のプロット matplot(c(1,3,5,7,9,11,13), py[,1:9], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T0") matplot(c(1,3,5,7,9,11,13), py[,10:19], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T1") matplot(c(1,3,5,7,9,11,13), py[,20:27], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T2") matplot(c(1,3,5,7,9,11,13), py[,28:36], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T3") # ウィンドウの分割の設定を元に戻す # 元に戻さない限り、ずっと2*2に分割されたままである par(mfrow=c(1, 1))
図6 群ごとに原データをプロット
ここでもランダム切片モデルとランダム切片・傾きモデルの2つを考えてみる。そして、それらをanova()によって比較してみる。
# ランダム切片モデル
> lme.model.eins <- lme(KALIUM ~ TIME + GROUP,
+ random=~ 1|SUBJECT, method="ML", data=mydata2)
# ランダム切片・傾きモデル
> lme.model.zwei <- lme(KALIUM ~ TIME + GROUP,
+ random=~ TIME|SUBJECT, method="ML", data=mydata2)
# 2つのモデルを比較する
> anova(lme.model.eins, lme.model.zwei)
Model df AIC BIC logLik Test L.Ratio p-value
lme.model.eins 1 7 391.1300 415.8360 -188.5650
lme.model.zwei 2 9 363.2297 394.9946 -172.6148 1 vs 2 31.90035 <.0001
# lme.model.zweiの係数表の出力
> summary(lme.model.zwei)
Linear mixed-effects model fit by maximum likelihood
Data: mydata2
AIC BIC logLik
363.2297 394.9946 -172.6148
Random effects:
Formula: ~TIME | SUBJECT
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.38429029 (Intr)
TIME 0.05189712 -0.152
Residual 0.37646814
Fixed effects: KALIUM ~ TIME + GROUP
Value Std.Error DF t-value p-value
(Intercept) 4.177192 0.15175156 215 27.526519 0.0000
TIME 0.039137 0.01059203 215 3.694940 0.0003
GROUPT1 -0.709064 0.20452632 32 -3.466858 0.0015
GROUPT2 -0.393634 0.21629768 32 -1.819873 0.0781
GROUPT3 -0.477967 0.20983957 32 -2.277774 0.0296
Number of Observations: 252
Number of Groups: 36
anova()によって2つのモデルを比較すると、AIC、BICともにlme.model.einsよりlme.model.zweiの方が小さな値になっているので、ランダム切片・傾きモデルの方がより良いモデルであることを示している。
尤度比検定においても、検定等計量（χ2 = 31.90035）についてのp値もp < 0.001となっているので、統計学的にModel2の方が良いモデルであるということがいえる（仮にp > 0.05であるようならばModel1の方が良いモデルであると判断できるようである）。
なおp値を計算したい場合は次のようにして求めることができる。
> 1 - pchisq(31.90035, 9-7) [1] 1.182843e-07
最後にランダム切片・傾きモデルによる予測値のプロットを描いてみよう。
# これでdatの行列と同じ形式で予測値が格納される yosoku <- matrix(predict(lme.model.zwei), ncol=7, byrow=T) t.yosoku <- t(yosoku) # ウィンドウの分割 par(mfrow=c(2, 2)) # T0のプロット matplot(c(1,3,5,7,9,11,13), t.yosoku[,1:9], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T0") matplot(c(1,3,5,7,9,11,13), t.yosoku[,10:19], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T1") matplot(c(1,3,5,7,9,11,13), t.yosoku[,20:27], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T2") matplot(c(1,3,5,7,9,11,13), t.yosoku[,28:36], type="l", xlab="time", ylab="kalium", ylim=c(0, 10), main="T3") # ウィンドウの分割設定を元に戻す par(mfrow=c(1, 1))
図7 ランダム切片・傾きモデルによる予測値のプロット
線形混合効果モデル（関数lmeを使う）と一般線形モデル（関数lmを使う）の違いを述べるとともに、経時測定データ（繰り返し観測データ）を独立なデータとして一般線形モデルで解析すると、どのようなことが生じうるのかを説明する。
改めてp2-表1.1 table11.xlsのデータセットを2つのモデルで解析して比較してみる。2つのモデルとは、ランダム切片モデル（ランダム切片・傾きモデルより良いモデルだと判断されたのであった）と一般線形モデルの分散分析モデルである。
# 分散分析モデルを用いた"誤った"モデル解析
> wrong.model <- lm(BP ~ Time, data=mydata)
> summary(wrong.model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 126.667 2.607 48.594 <2e-16 ***
Time1h -5.500 3.686 -1.492 0.156
Time3h -9.333 3.686 -2.532 0.023 *
# ランダム切片モデルを用いた"正しい"モデル解析
> model.eins <- lme(BP ~ Time, random=~ 1|Subject, method="ML", data=mydata)
> summary(model.eins)
Linear mixed-effects model fit by maximum likelihood
Fixed effects: BP ~ Time
Value Std.Error DF t-value p-value
(Intercept) 126.66667 2.606616 10 48.59430 0.0000
Time1h -5.50000 2.327771 10 -2.36278 0.0398
Time3h -9.33333 2.327771 10 -4.00956 0.0025
まず推定された偏回帰係数を見比べてみると、どちらのモデルでも全く同じであるが、これはたまたま一致したにすぎない。また、各偏回帰係数に対するp値を見比べてみる。そうすると両モデルでp値が異なっているのが確認できるだろう。
例えばTime1hのp値を見てみると、wrong.modelにおいてはp = 0.156と有意ではない。一方でmodel.einsではp = 0.0398と有意である。このように仮に（たまたま）偏回帰係数の値が一致したとしても、それに対するp値が異なることもある。もちろん、正しいのはmodel.einsのほうである。
それではp161-表10.1 table101.xlsのデータセットではどうだろうか。こちらの例でも、両モデルで推定された偏回帰係数は同じであるが、それに対するp値が異なる。例えば、GROUP2についていえば、lme.wrong.modelでは有意であるが、lme.model.zweiでは有意ではない。もちろん、正しいのはlme.model.zweiのほうである。
# 誤ったモデル解析
> lme.wrong.model <- lm(KALIUM ~ TIME + GROUP, data=mydata2)
> summary(lme.wrong.model)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.34826 0.10734 40.508 < 2e-16 ***
TIME 0.03914 0.01010 3.876 0.000136 ***
GROUPT1 -1.06937 0.11135 -9.603 < 2e-16 ***
GROUPT2 -0.48294 0.11776 -4.101 5.59e-05 ***
GROUPT3 -0.68254 0.11425 -5.974 8.02e-09 ***
# ランダム切片・傾きモデル
> lme.model.zwei <- lme(KALIUM ~ TIME + GROUP,
+ random=~ TIME|SUBJECT, method="ML", data=mydata2)
> summary(lme.model.eins)
Fixed effects: KALIUM ~ TIME + GROUP
Value Std.Error DF t-value p-value
(Intercept) 4.348264 0.17148397 215 25.356679 0.0000
TIME 0.039137 0.00697021 215 5.614878 0.0000
GROUPT1 -1.069365 0.22660465 32 -4.719078 0.0000
GROUPT2 -0.482937 0.23964672 32 -2.015202 0.0523
GROUPT3 -0.682540 0.23249147 32 -2.935762 0.0061
このページの最初の方に擬似反復を発見するのは容易であると述べたのを覚えているだろうか。復習しておくと、ある同一の個体から繰り返してデータをとった場合、それはもはや独立なデータではありえない。そして、それをあたかも独立なデータであるとして解析すると誤った結果が得られてしまう（これは前節のlme()とlm()の違いで明らかにした）。繰り返しデータをあたかも独立なデータであるとする（つまり反復測定データとする）ことを擬似反復というのであった。
例えば、先ほどやったlme.wrong.modelの分散分析表を見てみると擬似反復であることが容易に分かるものである。
> summary(aov(lme.wrong.model))
Df Sum Sq Mean Sq F value Pr(>F)
TIME 1 6.176 6.176 15.021 0.0001363 ***
GROUP 3 39.140 13.047 31.733 < 2.2e-16 ***
Residuals 247 101.552 0.411
これによれば全自由度は1 + 3 + 247 = 251となるが、これは誤りである。なぜなら、全自由度はサンプルサイズから1を引いたものになる（n - 1）のであるから、正しくは36 - 1 = 35となるはずだからである（今回の実験で用いられた被験体の数（サンプルサイズ）は36である）。なぜn-1であるかというと、定数項（μ： 母平均）を推定しているからである（1つのパラメータを推定するごとに1つの自由度を失う）。
簡単にいえば、Residuals（誤差、でも正しくは残差と訳される）の自由度がサンプルサイズ（n）よりも大きな値を取るということは在り得ないのである。だから、用いられた被験体の数（サンプルサイズ）が36であるのにもかかわらず、誤差の自由度がそれよりもはるかに大きい247となるのは明らかにおかしいわけである。
したがって、正しい分散分析表はanova()を用いて得る必要がある。summary(aov(lme.model.zwei))とすると、上の分散分析表と全く同じものが返されてしまう。
> anova(lme.model.zwei)
numDF denDF F-value p-value
(Intercept) 1 215 2750.4094 <.0001
TIME 1 215 13.6526 0.0003
GROUP 3 32 4.1210 0.0140
結果としてTIMEもGROUPも有意であるといえるが、それぞれに対するp値の大きさは明らかに異なっている。
擬似反復については文献[1]の第8章、および文献[2]のp15により詳しい解説がなされているので参考にするとよいだろう。
lme()のさらなる使い方を紹介しよう。例えば、切片を0とする（原点を0とする）モデルや2次の項を含めた非線形モデルを解析することも可能である。
例えとして表3のデータを解析してみよう。これはAからFまでの7人の被験者を対象に、60分間走り続けたときの消費カロリー（単位: kcal）の変化量を記録したものである（これは仮想データである）。そしてこのデータをプロットしものが図8である。
| 10 | 20 | 30 | 40 | 50 | 60 | |
| A | 30 | 60 | 120 | 210 | 400 | 550 |
| B | 28 | 54 | 120 | 220 | 372 | 558 |
| C | 45 | 71 | 137 | 243 | 389 | 575 |
| D | 20 | 46 | 112 | 218 | 364 | 550 |
| E | 46 | 72 | 138 | 244 | 390 | 576 |
| F | 34 | 60 | 126 | 232 | 378 | 564 |
図8 時間経過に伴う消費カロリー量の変化
図8を見ると曲線関係があるようであるため、線形モデルを当てはめるのは適切ではないかもしれない。このような場合には非線形モデルを当てはめる必要があるだろう。具体的には経過時間Timeの2次の項をモデルに加えるということである。あるいは原データを対数変換することによって線形モデルを当てはめてもよいかもしれないが、ここでは2次の項を含んだ非線形モデルを当てはめることを考えてみる。
以下に図8のプロットのためのコードと非線形混合効果モデルの解析のためのコードを記す。
# 消費カロリーのデータ
Cal <- c(
30, 60, 120, 210, 400, 550,
28, 54, 120, 220, 372, 558,
45, 71, 137, 243, 389, 575,
20, 46, 112, 218, 364, 550,
46, 72, 138, 244, 390, 576,
34, 60, 126, 232, 378, 564
)
# 被験者を表すカテゴリカル型
Sub <- rep(1:6, rep(6, 6))
Sub <- as.factor(Sub)
# 経過時間を表す連続型
Time <- c(5, 10, 15, 20, 30, 60)
Time <- rep(Time, 6)
# 図8のプロット
y <- matrix(Cal, ncol=6, byrow=T)
matplot(c(5, 10, 15, 20, 30, 60), t(y), type="l", xlab="time", ylab="cal.")
# 非線形ランダム切片モデル
# I(Time * Time)が2次項である
nlme.model.eins <- lme(Cal ~ Time + I(Time * Time), random=~ 1|Sub, method="ML")
# 非線形ランダム切片・傾きモデル
nlme.model.zwei <- lme(Cal ~ Time + I(Time * Time), random=~ Time|Sub, method="ML", data=mydata3)
# 両モデルの比較
> anova(nlme.model.eins, nlme.model.zwei)
Model df AIC BIC logLik Test L.Ratio
nlme.model.eins 1 5 349.9494 357.867 -169.9747
nlme.model.zwei 2 7 353.9494 365.034 -169.9747 1 vs 2 1.325827e-08
p-value
nlme.model.eins
nlme.model.zwei 1
どうやら今回の場合においては、ランダム切片モデルのほうがよいモデルであるといえるようである。そこでnlme.model.einsの係数表を見てみよう。同時に予測値のプロットも見てみる。
> summary(nlme.model.eins)
Linear mixed-effects model fit by maximum likelihood
Random effects:
Formula: ~1 | Sub
(Intercept) Residual
StdDev: 0.001439335 27.18269
Fixed effects: Cal ~ Time + I(Time * Time)
Value Std.Error DF t-value p-value
(Intercept) -89.41254 13.884458 28 -6.439758 0
Time 18.50301 1.153183 28 16.045159 0
I(Time * Time) -0.12608 0.016753 28 -7.526158 0
Correlation:
Number of Observations: 36
Number of Groups: 6
# 予測値のプロット
y.hat <- matrix(predict(nlme.model.eins), ncol=6, byrow=T)
matplot(c(5, 10, 15, 20, 30, 60), t(y.hat), type="l", xlab="time", ylab="cal.")
図9 予測値のプロット
ここでランダム効果の部分を見てみてほしい。
Random effects:
Formula: ~1 | Sub
(Intercept) Residual
StdDev: 0.001439335 27.18269
Sub（被験者）のランダム効果、すなわち切片のバラつきを表す値であるIntercept StdDev = 0.001という点に注目すると、限りなく0に近く、ほとんど被験者間での切片のバラつきがないことを示している。
このモデルはランダム切片モデルなのであるから、最初から傾き（TimeとTime^2に対する係数）は被験者間で同じであることを仮定している。したがって、切片の分散が0に近ければ、全ての曲線はほぼ完全に重なり合う形になり、視覚上は一本の曲線が描かれているようになる。
要するにこの結果から、いずれの被験者も時間の経過に伴う消費カロリーの変化量は同じであるということが分かるのである。もし切片のバラつきが大きければ、各被験者の変化量（つまり曲線の形状）に違いはないが、被験者間での"全体的な差"が認められるわけである。その場合、いくつかの平行な曲線が描かれるようになる。
さらに、もしこれがランダム切片・傾きモデルであるならば、時間経過に伴う消費カロリーの変化量は被験者間で異なるというような結論が得られるわけである。この辺りのことは、数多くのデータを分析して経験を積むことでしか、判断が難しい部分も多いのである。
ところで先ほど行った解析（nlme.model.eins）は不十分なモデルである。なぜなら、理論的にスタート時（0分の時点）では消費カロリー量は0であると考えられるので、切片を0であるとするモデル（原点を通る回帰）を解析する必要があるからである。
# 原点を通る回帰：
# 切片を0と仮定したモデルである
> new.model <- lme(Cal ~ Time + I(Time * Time) + 0, random=~ 1|Sub, method="ML")
> summary(new.model)
Random effects:
Formula: ~1 | Sub
(Intercept) Residual
StdDev: 75.81647 28.45068
Fixed effects: Cal ~ Time + I(Time * Time) + 0
Value Std.Error DF t-value p-value
Time 17.367724 1.103341 29 15.741024 0
I(Time * Time) -0.110993 0.016238 29 -6.835396 0
Number of Observations: 36
Number of Groups: 6
# 図10のプロットのためのコード
y.hat <- matrix(predict(new.model), ncol=6, byrow=T)
matplot(c(5, 10, 15, 20, 30, 60), t(y.hat), type="l", xlab="time", ylab="cal.")
図10 切片を0と仮定したモデルの当てはめ
先ほど説明したように、切片が0であるとしたモデルのほうでは切片の分散（Intercept StdDev = 75.81647）は、ランダム切片モデルよりも大きな値をとっているので、図10を見ても視覚的にいくつかの平行な曲線が描かれているようになっていることが確認できる。
しかしながら、鋭い人は気づいたかもしれないがランダム切片モデルのAIC（349.9494）と比べてみると、切片が0であると仮定されたモデルのAIC（373.8834）のほうが大きくなっている。つまり前者のモデルのほうがよいモデルであることを示している（尤度比検定についてはよく知らない）。
> anova(nlme.model.eins, new.model)
Model df AIC BIC logLik Test L.Ratio
nlme.model.eins 1 5 349.9494 357.8670 -169.9747
new.model 2 4 373.8834 380.2175 -182.9417 1 vs 2 25.93406
p-value
nlme.model.eins
new.model <.0001
このような事態に遭遇すると（いや、実際にはこのようなことは決して少なくない）AICの値を基準にモデルを選択した方がよいのか、それとも理論的には切片を0としたモデルの方が好ましいのか、といったように迷ってしまうかもしれない。だが、これはまさに分析者の判断に委ねられる問題であり、どちらが絶対に正しいということはない。
例えば、似たようなデータセットに次数（今回は2次）を3次、4次と上げていけばデータポイントへの曲線の当てはまり自体は良くなるし、AICの値も小さくなる。実際、次数を限界まで上げていけば、いずれは完全にフィットした理論曲線が求まる。そしてそのモデルについてのAICは0になるだろう。
結局は分析者が論文に自分の考えを書いて、それが査読者や他の読者に認められるかということになってしまうかもしれない。私としては、「AICによれば切片を0としないランダム切片モデルのほうがよいモデルだと判断できるが、その場合だと切片が-89.41となっているのでこれは現実的に考えにくいモデルである。なぜならスタート時において、消費カロリー量がマイナスの値を取ること自体が現実に起こりえないことだからである。したがって、実質科学的な観点から考えても切片を0としたランダム切片モデルのほうを採用すべきであると考えた。」と主張したいところである。
最後に補足としてlme()を使ううえでの典型的なエラーを紹介しておこう。以下は先ほどの切片を0と仮定したランダム切片モデルとは別に、傾きにもランダム効果を組み込んだ、切片を0と仮定したランダム切片・傾きモデルを解析しようとした例である。
# 原点を通る回帰： # ランダム切片・傾きモデル > new.model2 <- lme(Cal ~ Time + I(Time * Time) + 0, random=~ Time|Sub, method="ML") 以下にエラー lme.formula(Cal ~ Time + I(Time * Time) + 0, random = ~Time | : nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (9)
よくワケの分からないエラーメッセージが出てきてしまったが、これはユーザーの関数の使い方が誤っているのではない。簡単にいうとプログラムが「計算できないよ」というものである。パラメータを求めるために「プログラムが解を見つけようと何度も反復して計算する」のであるが、これは反復計算をしたけども、結局、解は見つけられなかったというエラーメッセージである。
"iteration limit reached without convergence"というのは、直訳してみると「反復は限界に達して収束しなかった」というような訳になる。これが、先ほど述べた「プログラムががんばって解を求めようとしたが、結局、解を見つけることができなかった」ということにあたる。これは非線形モデルを扱う際にはけっこう頻繁に遭遇するエラーメッセージであるので覚えておくとよい。
これに対する改善策はないかと言いたくなるかもしれないが、残念ながら「解決策はほぼない」ということを断っておこう。だが、今回の解析例においては原データを対数変換して線形モデルとして扱うことによって解決できるかもしれない（自分で色々とデータを変換してみたり、それに対していくつかの考えられるモデルを当てはめてみるのはよい勉強になるだろう）。
参考文献
[1] Alan Grafen, Rosie Hails著 野間口・野間口訳「一般線形モデルによる生物科学のための現代統計学」共立出版（2007）
[2] Michael J. Crawley著 野間口・菊池「統計学：Rを用いた入門書」共立出版（2008）
[3] B. エヴェリット著 石田ほか訳「RとS-PLUSによる多変量解析」シュプリンガー・ジャパン（2007）
[4] 藤越ほか「経時データ分析」オーム社（2008）
[5] 岡田昌史偏「The R Book -データ解析環境Rの活用事例-」九天社（2004）
[6] Ita Kreft, Jan de Leeuw著 小野寺ほか訳「基礎から学ぶマルチレベルモデル」ナカニシヤ出版（2006）
