ページ

2016年3月19日土曜日

【R】 対数線形モデルを使った分割表分析(前編)

分割表分析については、カイ自乗検定AICなど基礎的なことを解説してきたが、本コラムではいよいよ対数線形モデルを取り扱う。今までと何が違うか。2元分割表より高次の分割表(高次の交互作用)を取り扱える、ということだ(これらは互いに異質な分析方法というのではなく、自然な拡張になっている、と理解してほしい)。その過程で、条件付き独立やグラフィカル・モデリングとの関係も併せて紹介していこう。
(何故か、wikipedia日本版には対数線形モデルの項目が無かった。英語版でもlog linear modelではなくlog linear analysisとなっているので注意されたし。)

対数線形モデルを3元分割表に基づいて解説していこう。質的尺度を持つ変数A,B,Cがあったとして、これらの3元分割表を構成する。一つ一つのセルの確率を、

(1)
\[
\log p_{ijk} = u + u_{i}^{a} + u_{j}^{b} + u_{k}^{c} + u_{ij}^{ab} + u_{jk}^{bc} + u_{ki}^{ca} + u_{ijk}^{abc}
\]
と8種のパラメータで表現する。これを対数線形モデルにおける「飽和モデル」(saturated model)あるいは「Fullモデル」と呼ぶ。(記号の上付き文字は累乗ではなく、記号節約のための単なる規約であると見てもらいたい。)

一般に対数線形モデルというと、8種のパラメータのうちいくつかを0と置いたもの(制約をかけたもの)ということになるが、このコラムでは制約のかけ方に一定のルールを課す。階層型モデル族、というものだ。



階層型モデル族

階層型モデル族というのは、「低次のパラメータを0と置いた」場合に、「対応する高次のパラメータも同時に0と置かなければならない」というルールに従うことだ。例えば、飽和モデルで、

\[
u_{ij}^{ab} \equiv 0
\]

と2次パラメータに制約をかけたものとしよう。すると、同時にijを含む3次のパラメータも以下のように制約をかけなければならない。

\[
u_{ijk}^{abc} \equiv 0
\]

出来上がるモデルは次のような6種のパラメータのモデルになる。

(2)
\[
\log p_{ijk} = u + u_{i}^{a} + u_{j}^{b} + u_{k}^{c} + u_{jk}^{bc} + u_{ki}^{ca}
\]

次にkiを含む項も消してみる。この場合は対応する高次項は既に消去されているので、

(3)
\[
\log p_{ijk} = u + u_{i}^{a} + u_{j}^{b} + u_{k}^{c} + u_{jk}^{bc}
\]

さらに、jkを含む項を消去する。やはり高次項はない。

(4)
\[
\log p_{ijk} = u + u_{i}^{a} + u_{j}^{b} + u_{k}^{c}
\]

この最後に得られたモデルは「主効果モデル(main effect model)」であるとか、独立モデルであるとか、Nullモデルであるとか、呼ばれる。これは変数A、B、Cが完全に独立であることを意味しており、これまで解説してきた2元分割表で独立性の判定の有無に使ってきたものである。

なお、番外編であるが、

\[
\log p_{ijk} = u
\]

というものも考えられる。これは、主効果さえ存在せず、どのセルも定数である(一様分布)という特殊ケースである。

条件付き独立と階層型モデル族

何故このような階層型分布族を考えるのか、また、何故セル確率を対数変換しているのか、その動機は独立性、条件付き独立性の概念に関係がある。

飽和モデルの全体を指数変換してみる。

\[
p_{ijk} = \tilde{u} \cdot \tilde{u}_{i}^{a} \cdot \tilde{u}_{j}^{b} \cdot \tilde{u}_{k}^{c} \cdot \tilde{u}_{ij}^{ab} \cdot \tilde{u}_{jk}^{bc} \cdot \tilde{u}_{ki}^{ca} \cdot \tilde{u}_{ijk}^{abc}
\]

ここで、
\[
\tilde{u} = \exp(u)
\]

である。ここで、iとjという二つに着目して、全体をiを含む項とjを含む項を完全に分離できるかを考えてみる。すぐに、

\[
\tilde{u}_{ij}^{ab}, \; \tilde{u}_{ijk}^{abc}
\]
の二つがあるので、分離不可能、と分かるだろう。

では2番目のモデルについて同様に考えてみよう。

\[
\begin{eqnarray} p_{ijk} & = & \tilde{u} \cdot \tilde{u}_{i}^{a} \cdot \tilde{u}_{j}^{b} \cdot \tilde{u}_{k}^{c} \cdot \cdot \tilde{u}_{jk}^{bc} \cdot \tilde{u}_{ki}^{ca} \\ & = & \left( \tilde{u} \cdot \tilde{u}_{i}^{a} \cdot \tilde{u}_{k}^{c} \cdot \tilde{u}_{ki}^{ca} \right) \left( \tilde{u}_{j}^{b} \cdot \tilde{u}_{jk}^{bc} \right) \end{eqnarray}
\]

こうして、全体の確率はiを含む項とjを含む項に分離できることが分かる(kは両方にあっても良い)。これは変数Cを固定すると、変数Aと変数Bが独立となることを示しており、「条件付き独立」と表現される。記号的に表現すると、
\[
A \perp \! \! \! \perp B | C
\]

補足:「カテゴリカルデータのモデル分析」の中で坂元慶行先生が「対数線形モデルは無駄が多い」なる批評をしているのは、たぶん階層型モデルに限定しない対数線形モデル全体に向けられた言葉であり、全否定する意図はなかったものと理解している。

グラフと階層型モデル族

対数線形モデルの階層型モデル族の良いところは、モデルが変数間の関係を示す無向グラフと一対一の関係になることだろう。グラフというのは頂点と辺からなる図形のことで、頂点が変数、辺(の欠如)が条件付き独立関係を示す。

Rでグラフを扱う準備として、管理者権限で以下のようにパッケージをいくつか導入しておく。

install.packages("devtools")
require(devtools)
source("http://bioconductor.org/biocLite.R")
biocLite("Rgraphviz")
biocLite("RBGL")
biocLite("gRim")
biocLite("gRain")

試しに飽和モデルを表現してみる。結果は冒頭に示した図のようになる。

library(gRbase)
library(Rgraphviz)
ug0 <- ug(~a:b:c)
plot(ug0)
2番目のモデルを表示するには、
plot(removeEdge("a","b",ug0))
結果は右図のようになる。

まとめ

今回は対数線形モデルがどのようなものであるのかを説明したが、実際に分析するところまで進まなかった。次回、具体的な例題を用いてR上で対数線形モデルによる分析例を示す。

なお、参考文献として二つ挙げておきたい。どちらも分かりやすく、お勧めである。
  1. 日本品質管理学会テクノメトリックス研究会 (編集)、グラフィカルモデリングの実際
  2. Søren Højsgaard, David Edwards, Steffen Lauritzen、Graphical Models with R (Use R!)