マルチレベルデータとは
「通常の」解析*1では、「ある患者のアウトカムは、他の患者のアウトカムから影響を受けない」と考えて解析している。もう少し正確に言えば、アウトカムの変化のうち、説明変数で説明できない誤差の部分に関しては「互いに独立である（mutually independent）」という前提で解析している。
では、「アウトカムの誤差が互いに独立である」とは考えにくい状況とはどういったものだろうか。
1つ目は、アウトカムが類似する傾向がある亜集団（=クラスターcluster、あるいはブロックblock）が存在する場合である。例えば、研究対象集団がいくつかので世帯で構成されているとき、食事や生活習慣などは世帯内で類似する傾向があるだろう。このようなとき、世帯が異なる研究参加者同士は独立と考えてよいが、同一世帯内の参加者を独立と考えるのは無理がある。世帯の他にも、学校や地域などがクラスターになりうる。
2つ目は、同一対象者において反復してアウトカムが測定される場合である。この場合のデータの単位は「人」ではなく「評価時点」になる。同じ人から測定されたアウトカムは類似した傾向を持つと考える方が普通だろう。眼や四肢・指趾といった「臓器単位」あるいは個々の腫瘍といった「病変単位」でアウトカムをみるような研究の場合も、疾患のメカニズムによっては反復測定と同様に考える必要がある。各測定に対して、人がクラスター単位となっていると考えることもできる。
アウトカムの誤差の出どころが、「測定単位（観測単位）」とその上層の「クラスター単位」の複数箇所にわたる、という意味でマルチレベルデータと呼ばれることがある。解析手法の観点からは、以下の用語はいずれも同じものを指していると考えて差し支えない。
- クラスターデータ（clustered data）
- 反復測定データ（repeated-measurement data）= パネルデータ（panel data）
- マルチレベルデータ（multi-level data）
マルチレベルデータの解析
マルチレベルデータを解析する際には、誤差の相関構造や階層性を考慮する必要である。それは以下の理由のため。
データの情報量を過大評価しないため：
アウトカムの誤差に相関があるということは、観測単位間で情報の一部が重複しているということ。これを考慮しないと情報量を水増ししていることになってしまう。
クラスターを考慮した適切なモデルを使うことで検出力アップ：
観測値の誤差の中でクラスター効果として説明できることがあれば、クラスター単位での異質性を考慮することで検出力が上がるかもしれない。
解析する方法としては大きく以下の2種類がある。
- 一般化推定方程式（Generalized Estimating Equation, GEE）
- 混合効果モデル（Mixed-effects model, MEM）
前者は「誤差間の相関構造を定義して解析する方法」で、個々の観測単位（の誤差）間にある相関構造行列を与えることで対応するアプローチ。
後者は「誤差の発生源を階層化する方法」で、上層（例：世帯）が与えられた条件下では下層（例：世帯構成員）の誤差は独立であると考えるアプローチ。
この記事では前者のみ解説する。後者は別の機会に。
一般化推定方程式（GEE）の概要
GEEといってもモデルの基本は一般化線形モデル（Generalized Linear Model, GLM）なので、分布族（family）やリンク関数（link function）についての説明は割愛。
前述のとおり、GEEでは上記に加えて「観測単位の誤差間にどのような相関構造を想定するか」を指定する。代表的な型は以下のとおり。
独立型：independent
観測単位間に無相関（=独立）を仮定するもの。GEEを使おうと思ったときにはこれは使わない。*2
交換可能型：exchangeable
同一クラスター内では同一の相関係数をもった構造（最もよく使われる）。
1次自己相関型：AR(1)
近い時点ほど相関が強くなる構造。同一対象者において3回以上繰り返して測定したときに使われることがある。
非構造化型：unstructured
それぞれのデータの組み合わせごとに個別に相関係数を設定した構造。*3
geepackパッケージで一般化推定方程式（GEE）を使う
サンプルデータの準備
ChickWeightというデータセットを使う。ヒヨコの餌と体重変化の関係性をみたデータらしい（ヒヨコ50匹、測定数578回）。
> head(ChickWeight) Grouped Data: weight ~ Time | Chick weight Time Chick Diet 1 42 0 1 1 2 51 2 1 1 3 59 4 1 1 4 64 6 1 1 5 76 8 1 1 6 93 10 1 1
ヒヨコ毎に体重は類似していると考える方が自然だろう。
geeglm()を使ってモデルを当てはめる
glm()の使い方が分かっていれば、追加で指定する引数は以下の3つくらい。
id：クラスターを識別する変数corstr：相関構造（"independence", "exchangeable", "ar1", "unstructured"）を指定する引数（ユーザーが直接指定したい場合はzcorという引数で指定できるみたい）。std.err：使用する分散推定方法を指定。- "san.se"：sandwich variance estimate（デフォルト）。
- "jack"：approxymate jackknife variance estimate
- "j1s"：1-step jackknife variance estimate
- "fij"：fully iterated jackknife variance estimate
3番目の分散推定方法について少しだけ。 推定された分散は係数βの標準誤差（推定誤差）に使われる。
デフォルトで使われるsandwitch法*4（White-Huber法）。通常は、誤差の分散は観測単位間で等しくないといけない（等分散性の仮定, homoskedasticity）のだが、この方法で推定された残差を使って誤差の分散を補正することで非等分散性（heteroskedasticity）を許容できる。 クラスター数が30以下と少ないときはデフォルトのsandwich法ではバイアスが生じうる。
Jackknife法は再サンプリング法の1つ*5で、元サンプルから反復して計算することで特定の仮定を置くことなく分散（など）を推定する方法。なので計算負荷が高くなりやすい（多分 "j1s"<"jack"<"fij"の順？）。
コードは以下のとおり。
fit <- geeglm(weight ~ Time, data = ChickWeight, id=Chick, family = gaussian(link="identity"), corstr = "exchangeable", std.err = "san.se")
結果を以下のとおり。
> summary(fit) Call: geeglm(formula = weight ~ Time, family = gaussian(link = "identity"), data = ChickWeight, id = Chick, corstr = "exchangeable", std.err = "san.se") Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) 27.845 1.983 197 <2e-16 *** Time 8.726 0.522 280 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation structure = exchangeable Estimated Scale Parameters: Estimate Std.err (Intercept) 1509 259 Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.475 0.0228 Number of clusters: 50 Maximum cluster size: 12
Timeが1増えるとweightは約8.7増えていると解釈できる。
指数変換が必要なとき
リンク関数がlogとかlogitのときは、係数の推定結果を解釈する際に指数変換が必要になる。
方法1：broomパッケージのtidy( )を使う
library(broom) tidy(fit, exponentiate = TRUE, conf.int = TRUE)
> tidy(fit, exponentiate = TRUE, conf.int = TRUE) # A tibble: 2 x 7 term estimate std.error statistic p.value conf.low conf.high <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 1.24e12 1.98 197. 0 25435648463. 6.04e13 2 Time 6.16e 3 0.522 280. 0 2216. 1.71e 4 Warning message: In tidy.geeglm(fit, exponentiate = TRUE, conf.int = TRUE) : Exponentiating coefficients, but model did not use a log or logit link function
ここではlink="identity"だったので指数変換は不要ですよ、とエラーメッセージが表示されている。
方法2：biostat3パッケージのeform( )を使う
library(biostat3) eform(fit, level=0.95, method="Delta")
exp(beta) 2.5 % 97.5 % (Intercept) 1.24e+12 2.54e+10 6.04e+13 Time 6.16e+03 2.22e+03 1.71e+04
おわりに
- geeglmはロバスト分散推定をしたいときにも使える。
参考資料
- 開発者の文献：Halekoh, U.; Højsgaard, S. and Yan, J (2006) The R Package geepack for Generalized Estimating Equations. Journal of Statistical Software, 15, 2, 1-11"
- CoxモデルにGEEアプローチを使いたいときはgeecureパッケージを参照のこと
