質的データの線型モデル −対数線型モデル・ロジット分析の紹介−  東京大学大学院人文社会系研究科 社会心理学専門分野 江利川 滋 0.線型モデルとは 今までみなさんに説明してきた手法は、一度に多くの変数を扱うため「多変量解 析」と呼ばれています。ごく大雑把にいって多変量解析には“独立変数で従属変数 を説明する”タイプと“似た変数・回答者などをまとめる”タイプの2種類があり ます。後者には、今までに説明した因子分析や主成分分析、対応分析、クラスター 分析などが含まれます。 一方、前者で代表的なのは重回帰分析や判別分析、分散分析などですが、これら の手法が“独立変数で従属変数を説明する”仕方は「線型モデル(liner model) 」 として一般化できます(註1)。線型モデルは、従属変数Yと独立変数X1〜Xpの 間に“Yは適切に重み(パラメータという)θ1〜θpを付けたX1〜Xpの和で表現 される”という(下式のような)関係を仮定するものです。 Y=θ1X1+θ2X2+…θpXp+e ここで、上式のeは、YのうちX1〜Xpでは説明できない部分(誤差)を表しま す。また、X1〜XpやYは実際に測定できた値を用いるので、線型モデルの分析は θ1〜θp(及びe)の値を推測することが主眼となります。 線型モデルとして一括される各手法は、XやYにどんな種類の変数を用いるかに よって区別されます(下表参照)。そして、今回説明するのは、XとYがともに質 的変数の場合の「ロジット分析」です。 独立変数X 従属変数Y 分析手法 -------------------------------------------------- 量的変数 量的変数 → 重回帰分析 量的変数 質的変数 → 判別分析 質的変数 量的変数 → 分散分析 質的変数 質的変数 → ロジット分析 ところが、ロジット分析を説明するには「対数線型モデル(log-liner model)」 を説明しておく必要があります。ロジット分析は対数線型モデルに基づいて行うも のだからです。では、対数線型モデルとは何でしょうか? 註1 このように一般化した線型モデルを、「一般線型モデル(general liner model)」と呼ぶこともあります。また、因子分析や主成分分析、対応分析な ども線型モデルで表現できますが、混乱を避けるため触れていません。 1.対数線型モデルとは 1−1.対数線型モデルの考え方 対数線型モデルは、複数の質的な変数の関連を検討する手法です。“質的変数の 関連の検討”で代表的なのが2元クロス表に対して行うカイ二乗検定ですが、これ は3つ以上の質的変数の関連が検討できません。そこで対数線型モデルの出番です。 質的変数の分析ではクロス表が欠かせません。そこで、まず最初に“p個の選択 肢(またはカテゴリー)を持つ質問A”と“q個の選択肢(またはカテゴリー)を 持つ質問B”の2元クロス表を考えてみます(3元以上の話は後述)。表0-1 はA とBの普通のクロス表で、各セルには該当する回答者の人数(セル度数)が入って います。表の右端と下端には、各行及び各列ごとの合計度数(行度数及び列度数) があり、全体の合計はn..人(総度数)です。また、クロス表を構成するAやBな どの変数を、以下「構成変数」と呼ぶことにします(なお、構成変数という言い方 は私が勝手にそう呼ぶだけですから、他所では使わないように(^^;)。 表0-1 ! b1 b2 … bj … bq ! 計 ! ----+---------------------------+------+ 変数AとBの a1! n11 n12 … n1j … n1q ! n1. ! 実数クロス表 a2! n21 n22 … n2j … n2q ! n2. ! ! ! ! ai! ni1 ni2 … nij … niq ! ni. ! ! ! ! ap! np1 np2 … npj … npq ! np. ! ----+---------------------------+------+ 計 ! n.1 n.2 … n.j … n.q ! n.. ! ----+---------------------------+------+ 表0-1 の各セル度数を総頻度で割ると、表0-2 のようなセル相対度数からなるク ロス表ができます。このとき、0から1の間の値をとるセル相対度数fijは、“あ る回答者を調査したとき、その人のAへの回答がiで、Bへの回答がjである確率 (あるいは、セル(ai,bj) に属するような回答者を得る確率)”と見なすことが できます。 表0-2 ! b1 b2 … bj … bq ! 計 ! ----+---------------------------+------+ 変数AとBの a1! f11 f12 … f1j … f1q ! f1. ! 相対度数 a2! f21 f22 … f2j … f2q ! f2. ! クロス表 ! ! ! ai! fi1 fi2 … fij … fiq ! fi. ! fij=nij/n.. ! ! ! ap! fp1 fp2 … fpj … fpq ! fp. ! ----+---------------------------+------+ 計 ! f.1 f.2 … f.j … f.q ! 1 ! ----+---------------------------+------+ このとき、fijという確率はどのように決まるのかを考えてみましょう。 まず、上の表にはp×q個のセルがあるので、ある回答者は、構成変数のAやB への回答に関係なく、1/pqの確率でセル(ai,bj) に入り得ます。しかし、fij には、“その人のAへの回答がiである確率(=行相対度数fi.)”や“その人の Bへの回答がjである確率(=列相対度数f.j)”も関係すると考えられます。 さらに、“その人のAへの回答がiで、かつBへの回答がjである確率(これを 同時確率pijと呼ぶ)”の影響も考えられます。仮にAへの回答とBへの回答が確 率的に独立な事象ならば、pijはfi.とf.jで説明できる(pij=fi.×f.j)の で考えなくていいのですが、独立でないなら、pijは“fi.やf.jでは説明できな い部分”を表すものとして捉える必要があります。すると、AとBが独立かどうか によって、pijを考える立場と考えない立場の2つがあることになります。 以上から、仮に“AとBは独立でない”という立場に立つと、fijを説明するも のとして1/pq、fi.、f.j、pijの4つが考えられます。これらはみな確率的に 独立ですから、fijは以下のような式で表現されます。 式1 fij=1/pq・fi.・f.j・pij 一方、“AとBは独立である”という立場ではpijは考えませんから、fijを説 明するものは1/pq、fi.、f.jの3つとなり、fijは以下のように表されます。 式2 fij=1/pq・fi.・f.j 式1も式2も右辺の項が積の形なので、ちょっと見づらくなっています。そこで 両辺の対数をとると、式1は式3、式2は式4のようになります。 式3 log(fij)=log(1/pq)+log(fi.)+log(f.j)+log(pij) =u.. + ui. + u.j + uij 式4 log(fij)=log(1/pq)+log(fi.)+log(f.j) =u.. + ui. + u.j なお、u..=log(1/pq) ui.=log(fi.) u.j=log(fj.) uij=log(pij) 式3も式4も明らかに線型モデルで、fijを“構成変数に基づく項の対数”の線 型モデルで説明することから、これを「対数線型モデル」と呼びます(註2)。ま た、対数線型モデルでは一般にuのことを効果項といい、分散分析になぞらえて、 行相対度数や列相対度数に関するuを「主効果」、同時確率に関するuを「交互作 用」と呼びます。 註2 厳密にいうと、一般には、式1や式2(または式3や式4)の左辺はモデ ルの式に基づくfij(の対数)の推定値であり、必ずしもセル相対度数fij (の対数)そのものとは等しくありません。そこで以下では、対数線型モデ ルで推定されるfijの推定値を期待相対度数mijと書くことにします。 式3のように全ての主効果と交互作用の項を取り込んだモデルを「飽和モデル」 といいます。一方、式4のようにいくつかの交互作用の項を(さらには主効果の項 も)削除したモデルを「不飽和モデル」といいます。2つのモデルの違いは変数A と変数Bの関連を独立とみるかどうかの違いですから、どちらのモデルがより元の クロス表に妥当するのかを検討することは、とりもなおさず、AとBが独立かどう かを検討することと等しくなります。 質的変数が3つ以上ある場合も同様で、例えば上記のA、Bに、r個のカテゴリ ーをもつC(Cのカテゴリーを一般にkで表す)を加えると、飽和モデルとして式 5を考えることができます。これは、“A〜Cはどの2つの間にも、また3つ全体 の間にも関連がある”という仮説を表したモデルです。ちなみに、“AとB”や、 “BとC”といった2つの変数の間の交互作用を「一次交互作用」、A〜Cの3つ の変数の間の交互作用を「二次交互作用」といいます。 式5 log(mijk)=log(1/pqr) +log(fi..)+log(f.j.)+log(f..k) +log(pij.)+log(p.jk)+log(pi.k) +log(pijk) =u... ←定数項 + ui.. + u.j. + u..k ←主効果 + uij. + u.jk + ui.k ←一次交互作用 + uijk ←二次交互作用 一方、例えば“A〜Cは、3つ全体の間に関連はなく、さらにAとCの間にも関 連はない(その他の間には関連がある)”という仮説を立てたとすると、これに基 づく不飽和モデルは式6のようになります。 式6 log(mijk)=u... + ui.. + u.j. + u..k + uij. + u.jk そして、式5と式6の妥当性を検討することによって、式5が表現する仮説と式 6が表現する仮説のどちらが妥当であるかを検討します。式6のような、変数の関 連の仕方に関する細かい仮説が検討できる点が、対数線型モデルの利点といえます (註3)。 註3 対数線型モデルでは、効果項の組合せによっていろいろなバリエーション が生まれますが、これをいちいち“二次交互作用とAとCの一次交互作用が ないモデル”などと表記するのは面倒です。そこで、モデルの中で最も高次 の効果項だけを表記することで、特定のモデルを指すことがあります。 例えば、式3〜式6はそれぞれ以下のように表されます。 式3:[AB] / 式4:[A][B] 式5:[ABC] / 式6:[AB][BC] 最も高次の項だけ示してもモデル全体が分かるのは、対数線型モデルが 「階層の原則」に従っているからですが、これについては後述します。 1−2.不飽和モデルの妥当性の検討 ところが、今までの議論をひっくり返すようですが、実は飽和モデルは常に元の クロス表を完全に説明し、mijk=fijkとなることが分かっています。その理由は 飽和モデルには考え得る全ての効果項が含まれているからで、fの値を推測すると いう意味においては、飽和モデルは常に妥当なものなのです。では、なぜ不飽和モ デルなどを考えるのでしょうか? fが完全に予測できるからといって、常に飽和モデルだけ考えればいいわけでも ありません。飽和モデルには、実際は誤差として無視できるような効果項が含まれ ることもありますし、変数が増えるとモデルに含まれる項の数も爆発的に増えるの で、すぐに複雑になってしまいがちです(註4)。 註4 この辺の事情は重回帰分析にも似ています。重回帰分析では、独立変数を 増やせば増やすほどモデルの説明力は上がりますが、逆に無駄な変数も取り 込んでしまい、訳の分からないモデルになりがちです。 従って、飽和モデルと同じくらいの妥当性を持ちつつ(ここが重要)、より項の 数が少ない単純な不飽和モデルを検討する必要があるのです。では、“ある不飽和 モデルが飽和モデルと同じくらいの妥当性を持つかどうか”は、どうやって検討す ればよいのでしょうか? それには尤度比カイ二乗(likelihood ratio chi-square) 検定というものを行い ます。ここで尤度(ゆうど)とは、“対数線型モデルによって各セルの相対度数を 推定したときに、元のクロス表が再現される確率”のようなものだと理解して下さ い。例えば、飽和モデルではmijk=fijkなので、常に元のクロス表と全く同じも のが再現できるため、飽和モデルの尤度Lf は最大です。一方、飽和モデルでは一 般にmijk≠fijkなので、不飽和モデルの尤度Lw はLf より小さくなります。 ここで、Lw/Lfという比(尤度比LR)を考えてみます。飽和モデルが常に妥 当なのは分かっていますが、もしある不飽和モデルも同様に妥当であるなら、Lw はLf と同じくらいの値になり、LRは1に近くなると考えられます。逆に不飽和 モデルが妥当でないならLw は小さくなり、LRは0に近づいていくでしょう。こ のようなLRは以下のような変換によって、Lw垂kfである場合、その分布がカイ 二乗分布に近似されることが分かっています。このとき、この統計量を尤度比カイ 二乗(G^2)と呼びます。 G^2=−2log(LR) G^2はLRが1に近いほど値が小さくなり、LRが0に近いほど値が大きくなり ます。これを利用すると“不飽和モデルが飽和モデルと同じくらいの妥当性を持つ といえるかどうか”を統計的に検討できます。すなわち、G^2が有意に大きくなけ ればLR垂Pとみなして“不飽和モデルは飽和モデルと同じくらい妥当である”と 判断しますが、有意に大きい場合はLR<1なので“不飽和モデルは飽和モデルに 比べて妥当ではない”と判断します。こうした解釈の仕方(有意でない方が良い結 果)は普通の統計的検定の場合と逆なので十分注意して下さい。 繰り返しになりますが、以上のように検定を行なってG^2が有意でなかったとき、 その不飽和モデルは飽和モデルと同じくらいの妥当性を持つと解釈します。そして “同じ妥当性なら項の数が少ない単純なモデルの方がよい”という節約の原理に従 って、不飽和モデルを採用します。不飽和モデルの採用はすなわち、不飽和モデル が表している構成変数間の関連に関する仮説が支持されたことを意味します。 1−3.不飽和モデルの設定と吟味 以上では、線型対数モデルの考え方と、飽和モデル・不飽和モデルの区別、そし て不飽和モデルの妥当性の検討の仕方を説明しました。では、肝心の不飽和モデル はどのように設定すればよいのでしょうか? まず、どの効果項を削ればよいかを考えるわけですが、これを判断するための指 標がSASで出力されるので、それを参考にするのが1つの手です(これは「3. SASでの実行」で説明します)。もちろん、仮説に基づいて不要な効果項を削除 することもあります。 またこれとは別に、効果項の削り方には「階層の原則(hierarchy principle)」 という規則があります。これは“効果項は高次のものから順に削る”というもので、 例えば式5から不飽和モデルをつくる場合、二次交互作用→一次交互作用→主効果 →定数項の順に削らなければならず、二次交互作用を残して主効果を削るといった ことはできません。逆に言うと、高次の交互作用がモデルに含まれるためには、よ り下位の交互作用・主効果が全てモデルに含まれる必要があります。 実際の作業では、高次の効果から順に削り、そのつどG^2が有意でないことを確 認しつつ、残りの効果のうち高次のものからまた削ってみる、ということを繰り返 します。何回も効果項を抜き差しして、ようやく妥当な不飽和モデルが見つかると いうわけです。もちろん、妥当な不飽和モデルが見つからず、飽和モデルを採用せ ざるを得ない場合もあります。 さて、妥当なモデルを採用しました。次に、このモデルの各効果項を詳細に吟味 します。効果項には定数項、主効果、交互作用の3種類があることは前述の通りで すが、以下では順にその解釈の仕方を説明します(説明を簡単にするために、主に 2変数の場合を取り上げます)。 まず定数項は、全てのセルの期待相対度数(の対数)の平均で、全平均ともいい ます。あまり意味がないので、SASの出力では表示されない場合もあります。 1 定数項 :u..=----ΣΣlog(mij) pq i j 次に、一般にiで代表させる変数Aのカテゴリー(i=1,…,p)のうち、あるカテゴ リーs(=as )に関わる主効果us.は、“調査でas に属するような回答者を得る 確率の対数log(msj) を平均したものと全平均u..との差”で表されますが、こ の式を変形すると以下のようになります。 1 as の主効果:us.=----Σlog(msj)−u.. q j 1 1 =----・p・Σlog(msj)−----ΣΣlog(mij) pq j pq i j 1 =----ΣΣlog(msj/mij) pq i j 式の最後に注目するとus.は、“as に属する回答者を得る確率は、Aの他のカ テゴリーに属する回答者を得る確率の何倍か”という見込み(odds)msj/mijの対 数を、表全体について平均したものと考えられます(註5,註6)。 註5 以下では、“as に属する回答者を得る確率は、Aの他のカテゴリーに属 する回答者を得る確率の何倍かという見込み”のことを、単に「as の見込 み」ということにします。 註6 上式からも明らかですが、us を添字ごとに加えた合計は0となります。 この性質はSASで実行した結果の解釈の際、重要になります。 Σui.=0 i us はAのカテゴリー数だけ(すなわちp個)あります。SASではus.の推定 値と標準誤差(推定のブレの幅を示す)が算出されますが、推定値を標準誤差で割 った値は平均0、分散1の標準正規分布に従う統計量になります。これを用いて、 どのカテゴリーについてのus.が有意で、どのカテゴリーについてのus.が有意で ないかが統計的に検討できます。つまり、有意な正(負)のus.から“as に属す る回答者を得る確率ms.は偶然の確率(=1/p) より有意に大きい(小さい)”と 解釈します。 以上の考え方は、一般にjで代表させるBのカテゴリーのうちの1つ、bt に関 する主効果u.tについても全く同様に当てはまります。 1 1 bt の主効果:u.t=----Σlog(mit)−u..=----ΣΣlog(mit/mij) p i pq i j また、 Σu.j=0 j さらに、こうした主効果では説明しきれないlog(mst) の全平均からの偏りが 交互作用ustで、以下のような式になります。 (as,bt) :ust=log(mst)−u..−(us.+u.t) の交互作用 1 =log(mst)−----ΣΣlog(mij) pq i j 1 −----ΣΣ{log(msj/mij)+log(mit/mij)} pq 1 =----ΣΣlog{(mst/mit)/(msj/mij)} pq i j 式の最後に注目するとustは、“B=tのときのas の見込みはB=jのときの as の見込みの何倍か”という見込みの比(以下これを「as のbt での見込比」 と呼ぶ)の対数を、表全体で平均したものと考えられます。ustもp×q個あり、 その推定値を標準誤差で割った値が平均0、分散1の標準正規分布に従うため、値 の大きさを統計的に検定することができます。このとき有意な正(負)のustは、 “as とbt にともに属する回答者を得る確率は、偶然の確率(=1/pq)より大き い(小さい)”と解釈します。 また、交互作用も、添字ごとに加えた合計は0になります。 ΣΣuij=0 i j 以上は構成変数が2つの場合ですが、構成変数が増えるほどより高次の交互作用 が考えられます。例えば、構成変数が3つの場合に存在する二次交互作用を考えて みましょう。これを考えるために、変数AとBにr個のカテゴリーを持つ変数Cを 加えます。また、一般にCのカテゴリーをck で表し、そのうちの1つをcv とし ます。 このとき、主効果us..,u.t.,u,.v や一次交互作用ust.,us.v,u.tv で は説明しきれないlog(mstv)の全平均からの偏りが二次交互作用ustv で、以下 の式のようになります。 (as,bt,cv):ustv=log(mstv)−u... の交互作用 −(us..+u.t.+u..v)−(ust.+us.v+u.tv) 1 (mstv/msjv) (mstk/mijk) =----ΣΣΣlog{ ---------------- / ---------------- } pqr i j k (mitv/mijv) (msjk/mijk) 式の最後に注目するとustv は、“C=vのときのas のbt での見込み比は、 C=kのときのas のbt での見込みの何倍か”という見込比の比の対数を、表全 体で平均したものと考えられます。このustv はp×q×r個あり、推定値を標準 誤差で割った値が平均0、分散1の標準正規分布に従うので、その大きさを統計的 に検定できます。解釈も、有意な正(負)のustv は“as とbt とcv に同時に 属する回答者を得る確率は、偶然の確率(=1/pqr) より大きい(小さい)”と、 以前と同様です。 ここで注意したいのは、主効果が見込み、一次交互作用が見込比、二次交互作用 が見込比の比、…、とある高次の効果は一つ低次の効果の比として表されることで す。この性質は高次交互作用の意味を考える際に役に立ちます。 以上をまとめると、対数線型モデルによる分析の手順は以下のようになります。 1)飽和モデルを設定する。 2)効果項を削っていくつかの不飽和モデルを設定する。 3)飽和モデルとの比較から、採用する不飽和モデルを決定する。 4)不飽和モデルの詳細な吟味を行う。 1−4.対数線型モデルを行う際の注意点 以上の説明から、対数線型モデルはクロス表のセル相対度数の対数を推定する分 析だといえますが、相対度数が0のセルでは対数が取れないため、そのセルの対数 線型モデルは求まりません。そこで、こうした度数0のセルをどう扱うかが、対数 線型モデルの1つの注意点となります。 度数が0になるのは、“たまたまそのセルに属する回答者が得られなかった”場 合と“元から該当者がいないセルである”場合の2つがあり、前者を「標本ゼロ」、 後者を「構造ゼロ」と呼びます。構造ゼロの例としては、(変な例ですが)性別と ガンの種類のクロス表を作ると「男性」×「子宮ガン」のセルは必ず度数0になる といったものが典型的です(いい例があったら教えて下さい)。 SASでは最尤法という手法で対数線型モデルを求めるため、実際は度数0のセ ルでも対数線型モデルが求まりますが、逆に構造ゼロの場合、その旨をプログラム ではっきり指定しないと、全体として不適切な計算をすることになります。従って 標本ゼロと構造ゼロの区別はしっかりして下さい。 また、いくら最尤法を使うからといっても、クロス表の大半が0セルでは問題が あります。そこで、構造ゼロ以外の0セルを極力作らないよう、以下のような努力 が必要になります。 ・回答者を増やす ・むやみに構成変数を増やさない 特に後者については、解釈の容易さも考えると3〜4個の構成変数が限度でしょう。 ただし、このように努力したとしても、構成変数間に関連があるほど分布の偏り が大きくなり、度数の少ないセルが増えてきます。関連を見出そうとして分析する のに、関連があると0セルが増えるというやっかいな問題です。 2.ロジット分析とは 説明が非常に長くなりましたが、簡単にまとめると対数線型モデルとは、クロス 表の構成変数(とその組合せ)による効果項で各セルの相対度数(の対数)を推定 する線型モデルの吟味から、構成変数間の関連を検討するものでした。 普通、重回帰分析や判別分析といった他の多変量解析の線型モデルでは、対数線 型モデルでいう構成変数のうちの1つ(あるいはいくつか)が従属変数で、残りが 独立変数となっています。しかしここで注意したいのは、対数線型モデルではそれ らと違い、セル相対度数(の期待値)が“従属変数”であり、構成変数はいずれも “独立変数“であるという点です。このことからよく、対数線型モデルでは(構成 変数のレベルにおいて)独立変数と従属変数を区別しない、といわれます。 そうは言っても、対数線型モデルである構成変数を従属変数に、残りを独立変数 に見立てることもできないわけではありません。こうした見立てを行うのがロジッ ト分析です(註7)。 註7 ロジット分析で見立てた従属変数のことを、特に「反応変数」と呼んでい ます。これは、ロジット分析の基になる対数線型モデルの従属変数(セル相 対度数の推定値)と混同するのを避けるためです。 例えば、3つの構成変数A(カテゴリー数p個)、B(同q個)、C(同r個) のうち、AとBが独立変数でCが反応変数だとします。またこの3変数のクロス表 では、以下の不飽和モデル[AB][AC][BC]が採用されたとします。 log(mijk)=u...+ui..+u.j.+u..k+uij.+ui.k+u.jk このとき“反応変数Cのカテゴリー間の度数の違いが独立変数A,Bによってど の程度影響されているのか”という問題を、ロジット分析は以下の式で検討します。 P(ck|aibj) logitC(k:v) = log { ----------------- } P(cv|aibj) ここで、P(ck|aibj) とは、“Aのカテゴリーがi(ai) でBのカテゴ リーがj(bj) のときに、Cのカテゴリーがk(ck )である確率”というck についての条件付き確率を表し、“AとBの影響を受けたCのカテゴリーkの相対 度数”であると考えられます。このとき、2つのCのカテゴリー間でこのような条 件付き確率の比(の対数)をとって、カテゴリー間でのAやBの影響の差異を検討 するのです。Cのカテゴリーが2個しかない場合は1回の比較で済みますが、r個 (r≧3)ある場合は(r−1)回の比較が必要となります。 P(ck|aibj) は、対数線型モデルによるセルの期待相対度数とクロス表の 周辺相対度数を用いて、 P(ck|aibj)=mijk/mij. と書けるので、logitC(k:v)は結局、以下のような相対度数の見込みとなります。 mijk/mij. logitC(k:v) = log { -------------- } = log(mijk/mijv) mijv/mij. そこで、上式にA〜Cについての対数線型モデルを代入すると、logitC(k:v)は logitC(k:v) = log(mijk)−log(mijv) = u...+ui..+u.j.+u..k+uij.+ui.k+u.jk −(u...+ui..+u.j.+u..v+uij.+ui.v+u.jv) = [u..k−u..v]+[ui.k−ui.v]+[u.jk−u.jv] = w..k + wi.k + w.jk ただし、 w..k = [u..k−u..v] wi.k = [ui.k−ui.v] w.jk = [u.jk−u.jv] と表され、これをロジットモデルと呼ぶことにします(註8)。このとき、wをロ ジット効果係数とかw係数などといいます。またuについての制約から、 Σwi.k = Σw.jk = 0 i j となります。 註8 この場合、不飽和モデルが採用されているため、上の式のようになります が、仮に飽和モデルが採用された場合、wijk (=uijk−uijv)という項 が追加されます。 w係数は正負の値を取ります。そして、例えばw..k が正(負)の場合、“ck に属する回答者を得る確率はcv に属する回答者を得る確率より大きい(小さい)” と解釈します。また、wi.k の解釈では“Aがiのとき”、w.jk の解釈では“B がjのとき”という条件が付け加わります。 このようにロジットモデルのw係数は、反応変数のカテゴリー間での見込みを説 明変数の水準ごとに分解して示してくれます。しかし、通常の回帰モデルの係数に 比べ、対数線型モデルを基につくるロジットモデルでは、カテゴリーごとに係数の 解釈を行うため複雑になることは避けられません。 また、異なる対数線型モデルから導けるロジットモデルが同じになる可能性もあ ります。今の例の場合、次の不飽和モデル[AC][BC] log(mijk)=u...+ui..+u.j.+u..k+ui.k+u.jk でも全く同じロジットモデルが導けます。逆にいうと、1つのロジットモデルから は基となる対数線型モデルを必ずしも特定できないことを意味します(唯一の例外 は飽和モデルを採用したとき)。これは、ロジットモデルに変換する際に反応変数 の関与しない項は相殺されて残らないためです。 普通、基になる対数線型モデルが違えば効果項の推定値も異なりますから、当然 ロジットモデルのw係数の推定値にもこれが反映されると思われます。従って、同 じロジットモデルによる独立変数の影響力分析でも、w係数がデータの構造(特に 独立変数間の構造)を十分反映しない可能性があります。 これを避けるためには、まず“あてはまり”のよい対数線型モデルを決めて、そ れを基本にロジットモデルを導くという手順が必要です。 3.SASでの実行 お待ちどうさまでした。いよいよ、SASで対数線型モデル分析を行う手順を紹 介します。まず、例題として、以下のデータを分析することにします。これは松田 (1988)で紹介されていたもので、カナダでの犯罪に関する統計資料です(変な例で 済みません)。構成変数は「A:加害者の性別」・「B:加害者の年令」・「C: 殺害事件の分類」で、Cは“加害者が被害者の親族である場合(c1)”、 “強盗 やその他の犯罪を実行中に加えた傷害が死に至った傷害致死(c2)”、 “その他 (c3)” の3カテゴリーです。また以下の表の数字は該当者の人数です。 --------+-----------+----------------------------------------- ! ! C:殺害事件の分類 A:性別 ! B:年令 !----------------------------------------- ! ! 1)親族 2)強盗殺人など 3)その他 --------+-----------+----------------------------------------- 1)男性 ! 1)<20 ! 92 110 144 ! 2)20〜29 ! 175 191 389 ! 3)30〜39 ! 212 85 217 ! 4)≧40 ! 263 34 188 2)女性 ! 1)<20 ! 18 2 8 ! 2)20〜29 ! 68 7 22 ! 3)30〜39 ! 69 4 17 ! 4)≧40 ! 46 3 2 --------+-----------+----------------------------------------- さて、このクロス表に対数線型モデルを当てはめてみます。最初に、飽和モデル [ABC]を検討します。SASで対数線型モデル分析を行うには、CATMODプロシ ジャを用います。プログラムは以下の通りです。 --- SAMPLE PROGRAM IS AFTER THIS LINE ------------------------------------ TITLE 'Canadian Criminal Data from Maxim(1981)'; TITLE2 'Log-liner Model Analysis'; DATA D; INFILE IN; INPUT A B C; PROC CATMOD DATA=D; MODEL A*B*C=_RESPONSE_ / ML NOGLS NOITER NOPROFILE NODESIGN NORESPONSE; LOGLIN A B C A*B A*C B*C A*B*C; RUN; QUIT; -------------------------------------------------------------------------- 最初はタイトルの指定で、なくても構いません。DATAステップでは、A〜Cの3 変数を読み込んでいます。 そしてPROC CATMOD ですが、MODEL ステートメントで、左辺に全部の構成変数を アスタリスク'*' でつないで並べ、右辺には"_RESPONSE_"と書き(これは、飽和モ デルであれ不飽和モデルであれ、共通して必ずこう書きます)、オプションを指定 します。ここで指定しているオプションは、以下のような意味を持っています。 ML 最尤法による効果項の推定値を指定(必須)。 NOGLS 一般化最小二乗法による効果項の推定を行わない(必須)。 NOITER 最尤分析での計算の途中経過を出力させない。 NOPROFILE 構成変数のカテゴリーの組合せとセル度数の情報を出力させない。 NODESIGN 計算に用いる計画行列を出力させない。 NORESPONSE 計算に用いる_RESPONSE_行列を出力させない。 このうち、MLとNOGLS は、対数線型モデル分析を行う際に必ず一緒に指定します。 後は、あってもなくても構いませんが、付けておくと無駄な出力を省けます。この 他にも、以下のようなオプションもあります(まだ他にもあります)。 FREQ 構成変数によるクロス表(セルの値は度数)を出力させる。 PROB 構成変数によるクロス表(セルの値は相対度数)を出力させる。 PRED=FREQ 各セルの期待度数(実際の人数の予測値)を出力させる。 PRED=PROB 各セルの期待相対度数(確率の予測値)を出力させる。 CORRB 効果項の推定値の推定相関行列を出力させる。 COVB 効果項の推定値の推定共分散行列を出力させる。 次に、LOGLINステートメントでモデルを指定します。飽和モデルの場合、主効果 や交互作用を全て書きます。例えば、二次交互作用は“A*B*C” といった書き方に なります。 最後の“QUIT;” はCATMODの会話モードを終了させるものですが、なくても構い ません(詳しくはマニュアルを参照のこと)。 以上を実行すると、以下のような出力が得られます。 -------------------------------------------------------------------------- Canadian Criminal Data from Maxim(1981) Log-liner Model Analysis CATMOD PROCEDURE Response: A*B*C Response Levels (R)= 24 Weight Variable: WT Populations (S)= 1 Data Set: D Total Frequency (N)= 2366 Observations (Obs)= 24 MAXIMUM LIKELIHOOD ANALYSIS OF VARIANCE TABLE Source DF Chi-Square Prob -------------------------------------------------- A 1 448.00 0.0000 B 3 53.94 0.0000 C 2 136.93 0.0000 A*B 3 5.41 0.1438 A*C 2 90.54 0.0000 B*C 6 33.64 0.0000 A*B*C 6 10.74 0.0969 LIKELIHOOD RATIO 0 . . ANALYSIS OF MAXIMUM LIKELIHOOD ESTIMATES Standard Chi- Effect Parameter Estimate Error Square Prob ---------------------------------------------------------------- _RESPONSE_ 1 1.2956 0.0612 448.00 0.0000 〜中略〜 23 0.1031 0.1616 0.41 0.5234 NOTE: _RESPONSE_ = A B C A*B B*C A*C A*B*C -------------------------------------------------------------------------- 最初の飽和モデルの検討でまず注目するのが、「最尤分散分析表」です。まず、 一番下に分析に用いたモデルと飽和モデルとの尤度比カイ二乗値G^2が示されます が、ここでは飽和モデルそのものを分析しているのでG^2は0であり、ピリオドで 示されています。 次に、“Source”のところに各効果項の名前が示され、それぞれに対するカイ二 乗検定の結果が“Prob”に示されています。これが、「1−3.不飽和モデルの設 定と吟味」のところで述べた“効果項削除のための判断指標”で、有意でない効果 項はモデルでの説明力が低いため、削除の候補となります。 この例では“A*B” と“A*B*C” が有意ではありませんが、階層の原則に従って とりあえず二次交互作用だけを削ることにします。ここで注意しておくと、1つの 効果項を抜くだけでモデルが変わり、他の効果項の推定値も変わるので、何を抜く と他がどう変わるかを確認するためにも、高次の効果項から1つづつ削り、そのつ どモデルの妥当性を検討しなければいけません。 MODEL ステートメントとLOGLINステートメントを以下のように変更します。 --- SAMPLE PROGRAM IS AFTER THIS LINE ------------------------------------ PROC CATMOD DATA=D; MODEL A*B*C=_RESPONSE_ / ML NOGLS NOITER NOPROFILE NODESIGN NORESPONSE COVB; LOGLIN A B C A*B A*C B*C; RUN; -------------------------------------------------------------------------- モデルからから“A*B*C” を削り、COVBオプションで効果項の推定共分散行列を 算出させます(推定共分散を出力させる理由は後述します)。再度プログラムを実 行すると、以下のような出力を得ます。 -------------------------------------------------------------------------- Log-liner Model Analysis MAXIMUM LIKELIHOOD ANALYSIS OF VARIANCE TABLE Source DF Chi-Square Prob -------------------------------------------------- A 1 586.95 0.0000 B 3 134.18 0.0000 C 2 175.44 0.0000 A*B 3 20.47 0.0001 A*C 2 140.18 0.0000 B*C 6 174.88 0.0000 LIKELIHOOD RATIO 6 10.40 0.1086 ANALYSIS OF MAXIMUM LIKELIHOOD ESTIMATES Standard Chi- Effect Parameter Estimate Error Square Prob ---------------------------------------------------------------- _RESPONSE_ 1 1.3089 0.0540 586.95 0.0000 2 -0.3773 0.0814 21.49 0.0000 3 0.6119 0.0566 116.80 0.0000 4 0.1963 0.0607 10.44 0.0012 5 0.7756 0.0588 173.91 0.0000 6 -0.9203 0.0930 98.00 0.0000 7 0.0938 0.0813 1.33 0.2486 8 -0.1605 0.0557 8.30 0.0040 9 -0.1476 0.0567 6.79 0.0092 10 -0.6350 0.0573 122.91 0.0000 11 0.4017 0.0900 19.91 0.0000 12 -0.3632 0.0647 31.47 0.0000 13 0.4824 0.0682 50.00 0.0000 14 -0.4166 0.0509 66.86 0.0000 15 0.2999 0.0572 27.47 0.0000 16 0.1382 0.0539 6.58 0.0103 17 -0.0959 0.0687 1.95 0.1627 NOTE: _RESPONSE_ = A B C A*B B*C A*C COVARIANCE MATRIX OF THE MAXIMUM LIKELIHOOD ESTIMATES 1 2 3 4 5 6 -------------------------------------------------------------------------- 1 0.00291874 -.00078459 0.00080129 0.00040080 0.00214837 -.00277864 2 -.00078459 0.00662409 -.00165122 -.00196295 0.00037056 -.00067380 3 0.00080129 -.00165122 0.00320576 -.00024588 0.00063419 -.00060993 4 0.00040080 -.00196295 -.00024588 0.00368957 -.00008031 0.00003979 5 0.00214837 0.00037056 0.00063419 -.00008031 0.00345914 -.00372892 6 -.00277864 -.00067380 -.00060993 0.00003979 -.00372892 0.00864190 7 0.00090223 -.00552948 0.00156065 0.00168742 -.00010966 0.00035643 8 -.00090795 0.00155924 -.00244319 0.00013415 -.00047496 0.00026931 9 -.00038702 0.00169229 0.00013886 -.00253037 0.00012920 -.00017573 10 -.00216223 -.00030237 -.00039531 0.00013858 -.00276962 0.00314174 11 0.00287909 0.00041154 0.00024352 -.00013543 0.00312682 -.00747795 12 0.00021224 -.00114818 0.00004438 0.00030839 0.00044679 0.00006641 13 -.00013627 0.00117493 0.00010378 -.00021026 0.00001092 -.00021664 14 -.00023118 0.00001667 -.00079372 0.00012526 -.00059991 0.00058282 15 0.00015619 0.00012246 0.00099910 -.00011951 0.00057718 -.00117870 16 -.00008056 0.00030589 0.00010680 -.00129852 -.00025882 0.00018794 17 0.00006934 -.00020560 -.00010995 0.00161424 0.00021286 -.00036903 〜以下略〜 -------------------------------------------------------------------------- まずG^2(=10.40) をみると5%水準で有意でないので、ここで分析した不飽 和モデル[AB][AC][BC]は妥当なようです。各効果項のカイ二乗検定の 結果も全部有意なため、これ以上削れる効果項はなさそうです。 そこで、次に各効果項の吟味に移ります。効果項uの推定値は「最尤推定分析」 の表の“Estimate”の欄にあり、対応する標準誤差とカイ二乗値、棄却水準なども 示されています。ところが、この表をみるのは案外面倒なのです。それは、全ての 効果項についての数字が示されているわけではないからです。 上の表では17個の効果項が推定されていますが、どの項目がどの効果項を指して いるかを示すと、以下のようになります(最初の数字は表の“Parameter” に対応 しています)。 1 : u1.. ←Aの主効果 2 : u.1. 3 : u.2. 4 : u.3. ←Bの主効果 5 : u..1 6 : u..2 ←Cの主効果 7 : u11. 8 : u12. 9 : u13. ←A×Bの交互作用 10 : u1.1 11 : u1.2 ←A×Cの交互作用 12 : u.11 13 : u.12 14 : u.21 15 : u.22 16 : u.31 17 : u.32 ←B×Cの交互作用 これをみると、各変数の最後のカテゴリーに対応する主効果や、それが関係する 交互作用が抜け落ちていることが分かります。これは、主効果や交互作用の和が0 になるという性質から、最後のカテゴリーのことはその前までの情報があれは分か ってしまうため、SASではこうした冗長なカテゴリーを除いて計算していること によります。しかし、項の和が0になるために、簡単に残りのカテゴリーについて の推定値は求まります(式だけを示す;数字は後述)。なお、これらの欠けている 数字は手計算で求めるしかありません。 u2..=−u1.. u.4.=−(u.1.+u.2.+u.3.) u..3=−(u..1+u..2) u14.=−(u11.+u12.+u13.) u21.=−u11. u22.=−u12. u23.=−u13. u24.=−u14.=u11.+u12.+u13. u1.3=−(u1.1+u1.2) u2.1=−u1.1 u2.2=−u1.2 u2.3=−u1.3=u1.1+u1.2 u.13=−(u.11+u.12) u.23=−(u.21+u.22) u.33=−(u.31+u.32) u.41=−(u.11+u.21+u.31) u.42=−(u.12+u.22+u.32) u.43=−(u.13+u.23+u.33)=u.11+u.12+u.21+u.22+u.31+u.32 こうして効果項の推定値は求まりますが、さらに、各項の標準誤差も求めなけれ ば効果項の有意性の検定ができません。これも、効果項の和が0になるという性質 に基づいて、他の項から計算します。例えばu..3 の標準誤差は、以下のように求 めます(ここで、COVBオプションで求めた効果項の推定共分散を用いるわけです)。 u..3=−(u..1+u..2) ↓ var(u..3)=var[−(u..1+u..2)] =var(u..1+u..2) =var(u..1)+var(u..2)+2cov(u..1,u..2) ={sd(u..1)}^2+{sd(u)}^2 +2cov(u..1,u..2) =(.0588)^2 + (.0930)^2 + 2*(-.0037) = .0047 ∴ sd(u..3)=普i.0047)=.0686 ただし var(x)はxの分散、sd(x)はxの標準誤差、 cov(x1,x2)はx1とx2の共分散を表す 同様に、u.4. の標準誤差は以下のように求めます。 u.4.=−(u.1.+u.2.+u.3.) ↓ var(u.4.)=var[−(u.1.+u.2.+u.3.)] =var(u.1.+u.2.+u.3.) =var(u.1.+u.2.)+var(u.3.) +2cov[(u.1.+u.2.),u.3.] =var(u.1.)+var(u.2.)+2cov(u.1.,u.2.) +var(u.3.) +2cov(u.1.,u.3.)+2cov(u.2.,u.3.) =var(u.1.)+var(u.2.)+var(u.3.) +2cov(u.1.,u.2.)+2cov(u.1.,u.3.) +2cov(u.2.,u.3.) =(.0314)^2 + (.0566)^2 + (.0607)^2 + 2*(-.0017) + 2*(-.0020) + 2*(-.0002) =.000154 ∴ sd(u.4.)=普i.000154)=.0124 このようにして、SASでは算出されない効果項についても推定値と標準誤差を 手計算で求めて、全ての効果項の値などを整理すると以下のようになります(以下 の項目は順に、効果項の名前/推定値/標準誤差(括弧の中)/標準化推定値(= 推定値/標準誤差);* p<.05)。 Aの主効果              A×Cの交互作用 u1.. 1.309 (0.0540) 24.24 * u1.1 -0.635 (0.0573) -11.08 * u2.. -1.309 (0.0540) -24.24 * u1.2 0.402 (0.0900) 4.47 *                    u1.3 0.233 (0.0663) 3.51 * Bの主効果 u2.1 0.635 (0.0573) 11.08 * u.1. -0.377 (0.0814) -4.63 * u2.2 -0.402 (0.0900) -4.47 * u.2. 0.612 (0.0566) 10.81 * u2.3 -0.233 (0.0663) -3.51 * u.3. 0.196 (0.0607) 3.23 * u.4. -0.431 (0.0124) -34.76 * B×Cの交互作用                     u.11 -0.363 (0.0647) -5.61 * Cの主効果 u.12 0.482 (0.0682) 7.07 * u..1 0.776 (0.0588) 13.20 * u.13 -0.119 (0.0594) -2.00 * u..2 -0.920 (0.0930) -9.89 * u.21 -0.417 (0.0509) -8.19 * u..3 0.145 (0.0686) 2.11 * u.22 0.300 (0.0572) 5.24 * u.23 0.177 (0.0459) 3.86 * A×Bの交互作用           u.31 0.138 (0.0539) 2.56 * u11. 0.094 (0.0813) 1.16 u.32 -0.096 (0.0687) -1.40 u12. -0.161 (0.0557) -2.89 * u.33 -0.042 (0.0532) -0.79 u13. -0.148 (0.0567) -2.61 * u.41 0.642 (0.0609) 10.54 * u14. 0.214 (0.0661) 3.24 * u.42 -0.686 (0.0908) -7.56 * u21. -0.094 (0.0813) -1.16 u.43 0.045 (0.0624) 0.72 u22. 0.161 (0.0557) 2.89 * u23. 0.148 (0.0567) 2.61 * u24. -0.214 (0.0661) -3.24 *  結果の解釈ですが、まず主効果が有意であるというのは、単にカテゴリー間で 人数が有意に異なっていることを示すに過ぎません。例えば、a1.. が有意な正 の値をである(=a2.. が有意な負の値である)ことは、男女半々の状態よりも 有意に男性が多い(女性が少ない)ことを意味しています。他の主効果について も同様で、Bでは20代が多く40歳以上が少ないことが、Cでは親族の殺害が他の よりも多いことが示されています。  次にAとBの交互作用は、20〜30代では女性が多く、40歳以上では男性が多い 傾向を表すと解釈できます。  AとCの交互作用では、親族の殺害は女性の方が多いという、ちょっと興味深 い傾向が読みとれます(この解釈については、以下のロジットモデルで詳述)。  BとCの交互作用では、加害者の年令が高いほど、親族の殺害の割合が増えて いる傾向が分かります(ただし、20代→30代のところで減少傾向がある)。  では、以上の対数線型モデルを基に、AとBを独立変数、Cを反応変数にして ロジットモデルを構成してみましょう。ただし、ロジットモデルのw係数を求め るSASのプロシジャはありません。これも手計算で求めます。  まず親族による殺人は特別な意味を持つと考え、これを基準として相対度数の 見込みをとることにします。すると、2つのロジットモデル logitC(2:1) = log(mij2)−log(mij1) = w..2 + wi.2 + w.j2 logitC(3:1) = log(mij3)−log(mij1) = w..3 + wi.3 + w.j3 を得ます。w係数はu効果項と w..k = u..k−u..1 wi.k = ui.k−ui.1 w.jk = u.jk−u.j1 という関係があるので、計算は容易です。これから、以下のw係数が求まります。 k=2 k=3 ------------------------------------------ w..k w..2=-1.696 w..3=-0.631 wi.k w1.2= 1.037 w1.3= 0.868 w2.2=-1.037 w2.3=-0.868 w.jk w.12= 0.845 w.13= 0.244 w.22= 0.717 w.23= 0.534 w.32=-0.234 w.33=-0.180 w.42=-1.328 w.43=-0.597 ------------------------------------------ c1 (親族の殺害)を基準に、これと他の種類の殺害事件とを比較しているの で、解釈もそれを踏まえて行います。  まずCだけを取り上げて比較した場合(w..k) では、相対的に強盗殺人やそ の他が少ない、すなわち、相対的に親族の殺害が他に比べて多いといえます。  またAとCとの関係(wi.k) では、加害者が男性の場合、親族の殺害よりは 他の殺害事件を起こしている可能性が高く、逆に加害者が女性である場合、親族 の殺害をしている可能性が高いといえます。  しかし、これを単純に“女性の方が親族を殺しやすい”と解釈するのは誤りで す。元のクロス表に戻ると、男女で事件の発生件数自体が大きく異なっており、 男性は各種の殺害事件をたくさん起こしていますが、女性による殺害事件全体が 少ないことが分かります。そして、仮に女性による殺害事件をあったとすると、 それは親族の殺害である可能性が高い、と解釈する方が妥当でしょう。  さらに、BとCの関係(w.jk) では、加害者の年令が上がるに連れ、その他 の殺害事件が相対的に減り、親族の殺害が増える傾向が読みとれます(ただし、 w.13<w.23を例外とする)。  以上のような傾向は、親族の殺害と強盗殺人を比較した場合(c1 vs. c2) に顕著であるといえます。 4.最後に  以上で、対数線型モデルとロジット分析の理論、及びSASでの実行法につい て、最低限のことを説明しました。正直にいって、対数線型モデルは難しいと思 います。理論もそうでしょうし、SASで実行する場合も(非階層的クラスター 分析のように)何回もモデルを変えて繰り返す上に、手計算もかなりしなければ なりません。  しかし、複数の質的データの関連を検討できる、しかもモデルを構築して検討 することができるのは対数線型モデルくらいなもので、これほど便利かつ魅力的 な手法はないと思います。  私の能力では、これ以上、不正確にならない範囲で分かりやすく対数線型モデ ルを説明することはできません。しかし、この説明だけで対数線型モデルの全て を説明したわけではありません。興味のある人は、以下に挙げる参考文献を読ん で勉強してみて下さい。 5.参考文献  対数線型モデルを説明した本はいろいろありますが、独断と偏見により、ラン ク別に(?)参考文献を紹介したいと思います。なおこれらの文献の中には、稲 葉哲郎氏(現・東京大学文学部助手)から示唆されたものが多く含まれています。 稲葉さんに感謝します。 〈初級編〉 ・古谷野亘 1988 数学が苦手な人のための多変量解析ガイド 川島書店.   …いろいろな多変量解析の手法について、“何をする手法なのか”のイメー    ジがわくように、数式を極力使わずに説明している良い入門書です。ただ    し、この本で「対数線型モデル」と呼んでいるのは、私の説明でいうロジ    ットモデルのことなので注意して下さい。 〈中級編〉 ・B.S.エヴェリット(山内光哉監訳 弓野憲一・菱谷晋介訳) 1980 質的データの解析 カイ二乗検定とその展開 新曜社. ・G.J.G.アプトン(池田央・岡太彬訓訳) 1980 調査分類データの解析法   朝倉書店. …いずれも対数線型モデルの基本的な概念を説明しています。内容的には私  の説明と同程度だと思います。 ・A. DeMaris 1992 LOGIT MODELING -Practical Applications-, SAGE.   …これはSAGE社の "Series:Quantitative Applications in the Social Sci- ences"の中の1冊(No.86) で、対数線型モデルやロジット分析、さらに はロジスティック回帰分析までを含めた「ロジットモデル」全般について 非常に分かりやすく説明しています。 ・池田央編 1989 統計ガイドブック 新曜社. …まさに統計のガイドブックで、対数線型モデルに限らず、あらゆる統計的  手法についてコンパクトに説明しています。名前だけ聞いて何をする手法  か分からないとき、ちょっと調べるのに便利です。 〈上級編〉 ・松田紀之 1988 質的情報の多変量解析 朝倉書店.   …私が知っている中で、対数線型モデルの説明で最も詳しい日本語の本だと    思います。今回の私の説明は、この本の抜粋・要約といっても過言ではあ    りません。 ・SASインスティチュートジャパン 1993   SAS/STATソフトウェアユーザーズガイド SAS出版局.   …実は、今までに挙げた本の中で、SASのマニュアルがいちばん読みにく    いです(私にもよく分からないところが多々あります)。 _ 6.覚え書き  このレジュメは本来、東京大学文学部の社会心理学調査実習の資料であり、1994 年から1995年にかけて作成されました。現在は、柴内康文氏(現・東京大学大学院 人文社会系研究科)の好意により、オンライン上で公開されています。  Shigeru Erikawa(c)1996 All Rights Reserved.  このレジュメは、学術目的に限り複製・再配布を許可します。ただし、再配布す る場合は内容を改変しないでください。