土屋 - Researchmap

2011/08/20 第2回心理・医学系研究者のためのデータ
解析環境Rによる統計学の研究会
初心者による初心者の
ための「線形混合モデル」
※著作権に関わりそうな画像類は削除
土屋政雄
本日の内容
• 線形混合モデルって何？
– 名前がいろいろ
• どういう時に使うの？
– 横断と縦断
• どれ位有名なの？
– 推奨してる論文や使われ具合
• 今までの方法はだめなの？
– ANOVAや重回帰との違い
本日の内容
• 線形混合モデルを使うとどんないいこと
があるの？
• 中はどんな風になってるの？
– かんたんにすうしき
• 絵で描いてよ！
– グラフで図示
• Rでやってみてよ！
– データセット変換＆library(HSAUR)
• 論文に，書きたいです
線形混合モデルって何？
• Linear Mixed Model
• ざっくり説明すると，線形モデル
（ANOVAや重回帰）の進化版
• より細やかに
– 式が増えたよ！
• より柔軟に
– 必要な仮定が減ったよ！
本の表紙画像
Multilevel analysis
for applied research
It's just regression!
(イメージ)
Robert Bickel
線形混合モデルって何？
• Mixed models
今までの分析はこれだけ
– 固定効果（fixed effects）とランダム効果*（random
effects）の混ざったモデルのこと。
– 回帰係数（切片や共変量効果）が上位のユニット
でランダムにばらつくものとばらつかないものに
わかれる
– 一般的なマルチレベルモデルに対する特別な
ケースだと考えられるが，同じ意味で使われるこ
とも多い
*「変量効果」とも訳される
Diez Roux, 2002: Glossary for multilevel analysis. JECH;56:588-594
線形混合モデルって何？
ここにランダム効果が
混ざったら，一般化線
形混合モデルになる
一般化線形モデル
（Generalized linear model:GLLM）
一般線形モデル
（General linear model:GLM）
・t検定
・分散分析
・重回帰分析
・共分散分析等
ここに
ランダム効果が
混ぜ込まれた
・ロジスティック
回帰
・ポアソン回帰
・順序回帰等
図．伝統的な統計手法との関係
線形混合モデルって何？
• 名前がいろいろ
– マルチレベル分析（multilevel analysis）
– マルチレベルモデル・多水準モデル
– 階層（線形）モデル（hierarchical [linear] models）
→HLM
– ランダム（効果/係数/回帰）モデル
（random [effects/coefficient/regression] models）
– 共分散要因モデル（covariance components
models）
– 分散要因モデル（variance components models）
Diez Roux, (2002) ；小野寺編訳 2006 『基礎から学ぶマルチレベルモデル』
線形混合モデルって何？
• 名前がいろいろ
– 経験的ベイズモデル（Empirical Bayes models）
– 成長曲線モデル（growth curve models）
– 混合効果モデル（mixed-effect model）
• 間違えやすいもの，似ているもの（土屋の印象）
– 反復測定ANOVAのbetweenとwithinの混合
– 潜在成長曲線モデル→SEM
– 一般化推定方程式（GEE）
Gueorguieva & Krystal (2004)
どういう時に使うの？
• 下位レベルの集団が上位レベルの集団に
入れ子構造になっている時
– 横断でも縦断でも
・・・
学校A
クラスA
クラスB
・・・
・・・
・・・
（例）
・学校
・職場
・医者と患者
クラスC
・・・
どういう時に使うの？
• 「線形混合モデル」と言った時には，特に
縦断デザインで使われることが多いかも
→今日のテーマ
・・・
学校とか病院とか
（例）
・コホート研究
・介入研究
・・・
t1 t2
・・・
tn t1 t2
時点ごとの測定値
・・・
tn
t1 t2
・・・
tn
どれ位有名なの？
• 色んな領域のjournalで解説&紹介論文などが
出ている
– 発達研究・子ども（Boyle & Willms, 2001）
– 看護（Shin, 2009）
– 社会・パーソナリティ心理(Nezlek, 2008; West , 2011)
– 家族・カップル心理（Atkins, 2005）
– 精神生理（Kristjansson et al., 2007）
今までの方法はだめなの？
• 『無作為化試験の分析でANOVAはたやすく誤って
適用される：批判と代替統計アプローチの議論』
Vickers. （2005）Psychosom Med.
• 『ANOVAの地位をゆずれ：反復測定データ分析の
進歩とArchive of General Psychiatryの論文への反映』
Gueorguieva & Krystal (2004). Arch Gen Psychiatry.
• 『消えますよ・・・ほら消えた。：臨床試験
での縦断フォローアップデータの解析に
おける，従来の回帰 vs ランダム効果
回帰モデルの比較』
Nich & Carroll, (1997) J Consult Clin Psychol
ANOVA
is old!
今までの方法はだめなの？
ANOVAの典型的な問題
代替統計アプローチ
F値とp値の報告だけ
群ごとに平均と標準偏差を報告；
平均値の差と95%信頼区間の報告
3つ以上の群の試験の時，「全て
の群は同等であるか？」の仮説
のみへの言及
効果や興味のある群比較の後，係数
または平均値差と95%信頼区間とp値
の報告
ベースラインで取られた指標によベースラインの得点を共変量とした
る，試験への「群」効果，「時間」 ANCOVAを行い，群間の違いのみを
効果，「群×時間の交互作用」の報告する
報告
ベースラインと治療中or後の何回
かで取られた指標による，試験へ
の「群×時間の交互作用」の意
味は不明確
ベースラインの得点を共変量とした
ANCOVAを行い，従属変数の要約
（例：フォローアップ指標の平均）を使
う。または縦断混合効果モデルを使う
Vickers. （2005）Psychosom Med.
今までの方法はだめなの？
エンドポイント
分析
全対象者の完全
データが必要
rANOVA
rMANOVA
Mixed –Effects
Analysis
Y
N*
Y
N
欠損のある対象者
を取り除く影響
サンプルバイ
アス
サンプルバイ
アス
サンプルバイ
アス
非適用
欠損値代入の影響
推定バイアス
推定バイアス
推定バイアス
非適用
異なった時点での
測定
Y
N
N
Y
時間効果の記述
単純
柔軟
柔軟
柔軟
個人の傾向の推定
N
N
N
Y
相関パターンの制
約のある推定
非適用
Y
N
N
時間依存共変量
N
Y
N
Y
実行の簡単さ
とても簡単
簡単
簡単
難しい
計算の複雑さ
低
低
中
高
*欠損値のある者はデータセットから削除される
Gueorguieva & Krystal (2004)
今までの方法はだめなの？
（特にrANOVAとの比較）
• 欠損への対処のため，欠損値のある者を解析
から除外するか，値を代入する必要がある
– 前回の観測値を代入（last observation carried
forward: LOCF）→推定バイアスのリスク
• 全ての対象者が同じ時点でデータを得る必要
がある
Gueorguieva & Krystal (2004)
今までの方法はだめなの？
（特にrANOVAとの比較）
• 球面性の仮定（sphericity or circularity）が必要
– 全ての測定値が同じばらつきで，同一個人内で全
ての2つの連続した測定時点間で同じ相関である
こと
– 測定時点間の時間が長いと相関が弱くなるので，
この仮定が満たされることは少ない
• タイプIエラーの増大→効果の過大評価
– だめなら，Greenhouse-Geisserなどで自由度の修
• 保守的なので効果の検出が難しくなる
Gueorguieva & Krystal (2004)
線形混合モデルを使うとどんな
いいことがあるの？
• 各対象者の利用可能なデータをすべて使える
• ランダムな欠測には影響されない
• 時間の効果を柔軟にモデル化できる
• 現実的で倹約的な分散・共分散パターンが使
える（複合シンメトリー，自己回帰など）
• いくつかの異なったパターンのモデルを比較
することができる
– 適切な分散のパターンを選ぶとより正確な治療効
果やSEの推定につながる
Gueorguieva & Krystal (2004)
線形混合モデルを使うとどんな
いいことがあるの？
• 正確さがupすると必要なサンプルサイズが減
るので，検定力upにつながる
• 時間依存共変量の扱いが簡単
• 査読でつっこまれない！（ Hamer & Simpson, 2009）
– "Mixed models are preferable... when the number
of subjects is sufficiently large and the proportion
of missing data is small enough"
Gueorguieva & Krystal (2004)
中はどんな風になってるの？
基本のすうしき
※本や論文の著者により使われる記号が違うので，混乱しないように！
Yij   0 j  1 j Aij  Rij
Yij: 時点i，個人jのアウトカム指標
Aij：時点i，個人jのレベル1説明変数
β0j：個人jにおける個人特有の切片
β1j：個人jにおける個人特有の傾き
Rij：個人jにおける個人特有の誤差項
Blackwell et al., (2006)
中はどんな風になってるの？
追加
0 j   00  U 0 j
ランダム切片
1 j   10  U1 j
ランダム傾き
γ00: 全ての個人間の平均切片
γ10: 全ての個人間の平均傾き
U0j：レベル2誤差，つまり個人jにおけるγ00から説明でき
ない逸脱
U1j：レベル2誤差，つまり個人jにおけるγ10から説明でき
ない逸脱
Blackwell et al., (2006)
絵で描いてよ！
通常の回帰
ランダム切片のみ
モデル
絵で描いてよ！
ランダム傾きのみ
モデル
ランダム切片＆傾き
モデル
Rでやってみてよ！
• 使用データセット
– パッケージ “HSAUR” のデータセット “BtheB”
– コンピューターによる認知行動療法プログラムの臨
床試験データ（Proudfoot et al., 2003）
– 患者100名,変数8個のデータフレーム
◆Rプログラム
data(“BtheB”, package=“HSAUR”) #データを開く
head(BtheB) #最初の6名分のデータ表示
【R出力】
drug
1 No
2 Yes
3 Yes
4 No
5 Yes
6 Yes
length treatment bdi.pre bdi.2m bdi.4m bdi.6m bdi.8m
>6m
TAU
29
2
2
NA
NA
>6m
BtheB
32
16
24
17
20
<6m
TAU
25
20
NA
NA
NA
>6m
BtheB
21
17
16
10
9
>6m
BtheB
26
23
NA
NA
NA
<6m
BtheB
7
0
0
0
0
【変数】
drug: 抗うつ薬を使用したかどうか（No or Yes）
length: 現在の抑うつエピソードの長さ
（<6m，6カ月未満; >6m，6カ月以上）
treatment：治療群（TAU，通常治療；BtheB，Beat the Blues）
bdi.pre～bdi.8m：（ベック抑うつ尺度の治療前から8カ月フォ
ローアップまでの値）
◆id番号（subject）の付与と解析準備
BtheB$subject <- factor(rownames(BtheB)) # id番号付与
nobs <- nrow(BtheB) #オブザベーション数をnobsに格納
◆データセットを”wide”から”long”へ
BtheB_long <- reshape(BtheB, idvar = "subject",
varying = c("bdi.2m", "bdi.4m", "bdi.6m", "bdi.8m"),
direction = "long")
◆時間変数を作成
BtheB_long$time <- rep(c(2, 4, 6, 8), rep(nobs, 4))
◆idが1~3の者のデータのみ表示
subset(BtheB_long, subject %in% c("1", "2", "3"))
drug
1.2m No
2.2m Yes
3.2m Yes
1.4m No
2.4m Yes
3.4m Yes
1.6m No
2.6m Yes
3.6m Yes
1.8m No
2.8m Yes
3.8m Yes
length treatment bdi.pre subject time bdi
>6m
TAU
29
1
2
2
>6m BtheB
32
2
2
16
<6m
TAU
25
3
2
20
>6m
TAU
29
1
4
2
>6m BtheB
32
2
4
24
<6m
TAU
25
3
4 NA
>6m
TAU
29
1
6 NA
>6m BtheB
32
2
6
17
<6m
TAU
25
3
6 NA
>6m
TAU
29
1
8 NA
>6m BtheB
32
2
8
20
<6m
TAU
25
3
8 NA
同じidの者が4回出現しており，timeがそれぞれ2ヵ
月から8カ月まで変わり，bdiがそれぞれに対応した
値になっている
◆線形混合モデルの適用
library(“lme4”)
#線形混合モデルを行うためのパッケージ呼び出し
#ランダム切片モデル
BtheB_lmer1 <- lmer(bdi ~ bdi.pre + time + treatment +
drug +length + (1 | subject), data = BtheB_long, method
= "ML", na.action = na.omit)
#ランダム傾きモデル
BtheB_lmer2 <- lmer(bdi ~ bdi.pre + time + treatment +
drug + length + (time | subject), data = BtheB_long,
method = "ML", na.action = na.omit)
赤字部分がランダム効果を表す。1は切片，｜の後に
入れ子の上位ユニットを表す変数がくる。ここでは個人
◆モデルの比較
anova(BtheB_lmer1, BtheB_lmer2)
Models:
BtheB_lmer1: bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject)
BtheB_lmer2: bdi ~ bdi.pre + time + treatment + drug + length + (time | subject)
Df
AIC
BIC
logLik
BtheB_lmer1 8 1886.6 1915.7 -935.31
BtheB_lmer2 10 1889.8 1926.2 -934.90
◆線形混合モデルの結果出力
summary(BtheB_lmer1)
Chisq Chi Df
0.8161
2
Pr(>Chisq)
0.665
有意でないのでtime
のランダム傾きを含め
る意味はなさそう
Linear mixed model fit by maximum likelihood
Formula: bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject)
Data: BtheB_long
AIC BIC logLik deviance REMLdev
1887 1916 -935.3
1871 1866
Random effects:
Groups Name
Variance Std.Dev.
subject (Intercept) 48.304 6.9501
Residual
25.128 5.0127
Number of obs: 280, groups: subject, 97
Fixed effects:
Estimate
(Intercept)
5.94371
bdi.pre
0.63819
time
-0.71703
treatmentBtheB -2.37311
drugYes
-2.79786
length>6m
0.25639
Std. Error t value
2.24911 2.643
0.07759 8.225
0.14605 -4.909
1.66365 -1.426
1.71990 -1.627
1.63210 0.157
Computing standard errors:
Mixed-effects ML regression
Group variable: subject
Number of obs
Number of groups
=
=
280
97
Obs per group: min =
avg =
max =
1
2.9
4
参考：Stata
xtmixed bdir bdipre time treatmentr drugr lengthr || subject: , mle
Wald chi2(5)
Prob > chi2
Log likelihood = -935.31206
bdir
Coef.
bdipre
time
treatmentr
drugr
lengthr
_cons
.6381916
-.7170176
2.373078
-2.797836
.2563476
3.738992
Random-effects Parameters
Std. Err.
.0775914
.1460553
1.663748
1.719997
1.63219
4.774566
z
P>|z|
8.23
-4.91
1.43
-1.63
0.16
0.78
0.000
0.000
0.154
0.104
0.875
0.434
=
=
104.17
0.0000
[95% Conf. Interval]
.4861151
-1.003281
-.8878074
-6.168969
-2.942686
-5.618986
.790268
-.4307544
5.633964
.5732964
3.455381
13.09697
Estimate
Std. Err.
[95% Conf. Interval]
sd(_cons)
6.950145
.6088154
5.853704
8.251958
sd(Residual)
5.01274
.2599609
4.528264
5.549048
subject: Identity
LR test vs. linear regression: chibar2(01) =
126.38 Prob >= chibar2 = 0.0000
他にも必要な知識
• 級内相関（intraclass correlation: ICC）
• 推定法
– 最尤法（maximum liklihood： ML)
– 制限付き最尤法(restricted ML: REML)
• 共分散パターン
– unstructured，自己回帰，複合シンメトリーなど
• 中心化
– 例：個々の独立変数－平均値 →平均を0に
論文に，書きたいです
• 記述例：無作為割付比較試験の効果推定
– ANOVAタイプ
• いい例が見つからず・・・
– 回帰タイプ
• 例：Sumathipala et al., 2008 Cognitive-behavioural
therapy v. structured care for medically unexplained
symptoms: randomised controlled trial.
Br J Psychiatry. 193(1):51-9.
• Sumathipala et al., 2008
– Statistical analysis ・・・ analysed the scores from the four fixed
time points (3 months, 6 months, 9 months and 12 months after
randomisation) using a mixed effects model. We included as fixed
effects group allocation, baseline score, time, and the interaction of
group and time. Time was coded as months after the 3-month time
中心化 point, giving values of 0, 3, 6 and 9, so that in the presence of the
interaction the effect of group represents the difference at 3
months. We included the patient as a random effect. We also fitted
models with a random effect of time, with various covariance
patterns, with treating doctor as an effect, and pattern mixture
models to allow for the different drop-out patterns. We report here
the simpler models as they fit as well as any of the more complex
ones. We also examined model residuals. A mixed effects model
was used as this enables effective use of all the information even
from participants who had some missing scores. We used r for the
analysis,35 with the nlme package for fitting the mixed effects
models.36
• Sumathipala et al., 2008
Table 3
Fig 2
文献
• Atkins DC. Using multilevel models to analyze couple and
family treatment data: basic and advanced issues. J Fam
Psychol. 2005 Mar;19(1):98-110.
• Blackwell E, de Leon CF, Miller GE. Applying mixed
regression models to the analysis of repeated-measures data in
psychosomatic medicine. Psychosom Med. 2006;68(6):870-8.
• Boyle MH, Willms JD. Multilevel modelling of hierarchical
data in developmental studies. J Child Psychol Psychiatry.
2001;42(1):141-62.
• Diez Roux AV. A glossary for multilevel analysis. J Epidemiol
Community Health. 2002;56(8):588-94.
• Gueorguieva R, Krystal JH. Move over ANOVA: progress in
analyzing repeated-measures data and its reflection in papers
published in the Archives of General Psychiatry. Arch Gen
Psychiatry. 2004;61(3):310-7
• Hamer RM, Simpson PM. Last observation carried forward
versus mixed models in the analysis of psychiatric clinical trials.
Am J Psychiatry. 2009 Jun;166(6):639-41.
• Kristjansson SD, Kircher JC, Webb AK. Multilevel models for
repeated measures research designs in psychophysiology: an
introduction to growth curve modeling. Psychophysiology.
2007;44(5):728-36.
• 小野寺孝義編訳岩田昇・菱村豊・長谷川孝治・村山航訳
2006 基礎から学ぶマルチレベルモデル入り組んだ文脈から
新たな理論を創出するための統計手法ナカニシヤ出版
（ Kreft IGG, De Leeuw J ： Introducing multilevel
modeling Sage London, 1998）
• Nezlek, J. B. (2008). An Introduction to Multilevel Modeling for
Social and Personality Psychology. Social and Personality
Psychology Compass, 2(2), 842-860.
• Nich C, Carroll K. Now you see it, now you don't: a comparison
of traditional versus random-effects regression models in the
analysis of longitudinal follow-up data from a clinical trial. J
Consult Clin Psychol. 1997;65(2):252-61.
• Shin JH. Application of repeated-measures analysis of variance
and hierarchical linear model in nursing research. Nurs Res.
2009;58(3):211-7.
• Vickers AJ. Analysis of variance is easily misapplied in the
analysis of randomized trials: a critique and discussion of
alternative statistical approaches. Psychosom Med.
2005;67(4):652-5.
• West SG, Ryu E, Kwok OM, Cham H. Multilevel modeling:
current and future applications in personality research. J Pers.
2011;79(1):2-50.

Download Report