1 共分散分析（線形モデル） - 大学入試センター研究開発部

2010 年 10 月 8 日 (大学入試センター大津起夫 [email protected])
文書の所在 http://www.rd.dnc.ac.jp/~otsu/doc/SyakaitoKenko2010.pdf
平成 22 年度 (2010 年度) 大学院共通講義社会と健康 II (統計解析の基礎)
共分散分析（線形モデル）
1
共分散分析とは，連続変数と離散変数を共に説明変数とし，一つの連続変数を従属変数とする線
形の統計モデルを用いる分析のことである．回帰分析と分散分析 (ANOVA) は，ともに共分散分
析の特殊ケースである．
回帰分析と分散分析はともに，複数の変数（独立変数または説明変数と呼ばれる）によってひと
つの変数（従属変数または被説明変数）の値を予測する手法である．これらの分析の利用目的とし
て，次のようなケースが考えられる．
1. 従属変数の予測を独立変数の値を用いて行う.
2. 特定の独立変数が従属変数にどれだけの影響を与えているかを推定する.
3. 複数の独立変数の従属変数への影響について，それらの交互作用を推定する.
これらのうち，1 については大量のデータがあれば，多くの場合，適切に推定を行うことが可能で
ある（ただし，精度が良いかどうかは個別の状況に依存する）．しかし 2 と 3 とについては，信頼
できる結果を得るためには，研究者が独立変数の値をコントロールすることが重要である．理想的
には無作為割り当てなどを用いた実験 (randomized trial) を計画する必要がある．
疫学や社会科学における多くの調査のように，研究対象への介入を行わず，受動的なデータの観
測を行う場合には，予測は妥当であったとしても，未知の影響要因（交絡変数と呼ばれる）のため
に，計算によって得られる影響の大きさの推定値に偏りが含まれている可能性を否定できない．実
際の分析においては，このような未知の要因が影響を及ぼさないように慎重に分析の枠組みが設定
されるが，大規模な疫学調査においても，これによる結果の曖昧さは残るため，結果の解釈には十
分な注意が必要である．
ここでは，まず回帰分析と分散分析 (ANOVA) について説明し，次にこれらが同一のモデルの
特殊ケースであることを説明する．ついで離散変数と連続変数の両者を説明変数に持つ，共分散分
析 (ANCOVA) について，実例を示して説明する．
1.1
重回帰分析
説明変数が複数ある場合に，回帰分析のモデルを拡張する．y N ×1 は被説明変数であり，X N ×p
は説明変数とする．通常 X の第 1 列は定数 (全て 1) とする．行列 X の各列が１次独立の場合（多
くの場合には成立する）には，p は説明変数の個数になる（定数項をも 1 変数と数える）．
多変量の回帰分析のモデル (重回帰分析) は，行列を用いて記述すると次のようなものである．
y = Xβ + ε
(1)
これは、被説明変数 Y が次の式で表されることを表す．
Y = β1 X1 + · · · + βp Xp + ε
1
(2)
ここで，ε の要素 εi (i = 1, ..., N ) は，平均ゼロ，分散 σ 2 の正規分布に互いに独立に従うものと
する．
教科書や文献により、表記の方法に違いがあるが，ここでは X1 は定数を表す変数であり，常に
１であるとする．
後述の例では Y が新生児の体重であり、X1 , ..., Xp が母親の特性を表す変数（身長、体重、妊娠
期間、喫煙の有無など）である．被験者別に考えると、
yi = β1 xi1 + · · · + βp xip + εi , (i = 1, ..., N )
(3)
と表される．これをまとめて表現したものが 1 である．
残差を最小にする回帰係数 β は，行列（逆行列と転置）を用いて次の式で推定される．ここで
行列の −1 乗の記号は，
「逆行列」を表す．
ˆ = (X T X)−1 X T y
β
(4)
一般的に，ベクトルで表された複数の確率変数の平均について
E(Az) = AE(z), Var(Az) = AVar(z)AT
(5)
ˆ の平均は β であり，分散共分散行列は σ 2 (X T X)−1 となることがわかる．
が成立するので，β
また，1 が本当に成り立っているならば，次の２つが成り立つ．
ˆ の平均は真の値 β に一致する．
1. 最小 2 乗法によって求められた β
ˆ の分散は，観測値の１次式で表され, さらに平均が真の値に一致するもののうちで，最小
2. β
となる．
ˆ で与えられ，また誤差分散の不偏推定量は，
ˆ = Xβ
被説明変数の予測値は，y
ˆ )T (y − y
ˆ )/(N − p) = y − y
ˆ 2 /(N − p)
σ
ˆ 2 = (y − y
(6)
で与えられる．σ
ˆ 2 は，モデルの仮定が正しければ σ 2 × χ2N −p /(N − p) に従う（χ2 (k) は自由度 k
ˆ とσ
のカイ２乗分布を示す）．また，理論的な検討から，β
ˆ 2 とはモデルの仮定の下で統計的に独
立であることが知られている．
X が定数の列を含む場合，次の値を重相関係数の 2 乗 (R2 ) と呼ぶ．R2 は，１変数の場合の r2
と同様に次の式で定義される．
R2 =
ˆ−y
¯
y
¯
y−y
2
2
=1−
ˆ
y−y
¯
y−y
2
2
(7)
ˆ を意味す
R2 は 0 から 1 までの範囲をとる．R2 = 1 は，誤差がゼロであること，すなわち y = y
2
ˆ の分散がゼロ，すなわち yˆi = y¯ (i = 1, ..., N ) を意味する．
る．また，R = 0 は，y
1.2
説明変数の選択
回帰モデルの性質として，多くの説明変数を含むモデルの推定を行なうと，回帰係数の推定精度
が落ちることが知られている．多くの説明変数を用いて R2 を求めると見掛け上大きな値が与えら
ˆ を適用すると予測の性能は概して悪い．もしデー
れることがあるが，別のデータに推定された β
タが十分に大量にあればこのような問題は生じないが，複雑なモデルの推定について注意を払うこ
2
とが不要であるほど大量のデータは普通には利用できない．このような問題が生じるのは，回帰係
数がたまたま目前にあるデータについて日和見的に適合しているためである．これを避けるために
は，モデルの適合度は保ちつつ，極力少数の説明変数によって予測を行なう必要がある．回帰分析
における F 検定は，変数選択のための一つの方法である．
切片 β1 を含む p 個の回帰係数のうち，要素 βq+1 , ..., βp がゼロであるか否かを検討することにし
よう．M0 を βq+1 = · · · βp = 0 と制約を加えたモデル (説明変数の個数が q のモデル) とする．ま
た，M1 を βq+1 , ..., βp が必ずしもゼロでは無いモデルとする．モデル M0 によって推定された y
ˆ 0 とし，M1 による予測値を y
ˆ 1 とする．モデル M0 (より単純なモデルといえる) が正
の予測値を y
しいとの仮定の下で，次が成立する．
ˆ0
y−y
2
ˆ0
y−y
2
∼ σ 2 χ2N −q
(8)
ˆ1
y−y
2
∼
(9)
ˆ1
− y−y
2
∼ σ 2 χ2p−q
σ 2 χ2N −p
(10)
また M0 の仮定の下で，(9) と (10) とは独立である．また，このとき
F =
ˆ0 2 − y − y
ˆ1
y−y
2
ˆ1
y−y
2
×
N −p
p−q
(11)
は，自由度 (p − q, N − p) の F 分布に従う．もし，M0 が正しくなく，M1 が説明のために必要で
あるのならば，(8) の値は M0 が成立する場合よりも大きくなるはずである．そこで，F 値 (11) を
F 分布の値と比較し，もし顕著に大きければ「M0 が成立している」という仮説を棄却することに
する．通常，F 値が分布の上位 5 パーセントに入れば，M0 を棄却することが多い．
また，回帰係数 βj について，次のような検討を加えることが多い．以下で，行列 (X T X)−1 を
C = (cij ) と表記する．上に示したように，回帰係数 βˆj の平均は βj ，分散は σ 2 × cjj である．残
差の分散 σ 2 は未知であり，この推定値として σ
ˆ 2 を用いることにする．モデル M1 の仮定の下で，
βˆj − βj
t = √
σ
ˆ 2 cjj
(12)
は自由度 N − p の t 分布に従う．これを用いて，平均値の信頼区間の場合と同様に，βj の信頼区
間を求めることができる．もし，信頼区間がゼロを含めば，第 j 番目の変数を説明変数に加える必
要は薄いとみなされる．
また次の式で定義される自由度調整済み重相関係数も説明変数の決定に使われることが多い．
√
N −1
∗
(13)
R = 1 − (1 − R2 )
N −p
ここで N は標本数，p はデザイン行列 X のランク (通常は列数，重回帰分析の場合は説明変数の
個数＋ 1)．
この指標は，同一の R については p が大きくなるとともに小さくなる．少ない変数で大きな重
相関係数を持てば，その説明変数の組が望ましいものとする．
より単純で理論的に簡明なモデルの選択方法として，赤池情報量基準 AIC(Akaike’s Information
Criterion) を用いる方法がある．これは AIC とよばれる指標を最小化するモデルを最良モデルと
みなす方法である．誤差分散が未知の回帰モデルの場合，AIC は次の式で与えられる．
ˆ
y−y
N
ˆ
= N log y − y
2
AIC = N log
3
2
+ 2q + 定数
(14)
+ 2q + 定数
(15)
ここで，q はモデルに含まれるパラメータの個数であり，誤差分散未知の回帰モデルの場合には
p + 1 である．
また Schwarz(1978) は，ベイズ推定の立場から次のモデル選択基準 (のちに BIC と呼ばれるよ
うになった) を提案している．
BIC = N log
1.3
ˆ
y−y
N
2
+ q log N + 定数
(16)
回帰係数の解釈
回帰分析は，得られたデータの見掛け上の関係を分析するものであり，計算の手続きにそれ以上
の意味はない．しかし，多くの場合分析によって推定したいのは，変数間の影響関係である．これ
について検討するには，計算手続きだけでなく，データ取得のプロセスや，変数間に想定される影
響構造についての知識が必要になる．
回帰分析の結果（回帰係数）が説明変数が被説明変数に及ぼす影響とは解釈できない場合を考え
てみる．
1. 何らかの観測されない要因があり，これが説明変数と被説明変数の両者に影響を及ぼしてい
る場合．
ここで，自動車の走行速度と停止距離の関係を検討する実験を想定する．実験を低速度から
高速度に順に行ったと仮定する．このとき初め路面が濡れており，時間とともに乾燥してき
たとする．もし，路面の状態が観測されていなければ，速度の停止距離に及ぼす影響は，過
小評価されるだろう．
2. 分析の対象データが，説明変数と被説明変数の両者の影響を受ける変数によって層別または
選別されたものである場合．
大学入学後の成績と学力試験成績および高校の内申書の影響について検討する場合を想定し
よう．合否は学力試験の成績のみによって決まると仮定する．学力試験の成績と高校内申書
の得点の間には，相関があるが，完全に一致するわけではない．ここで，入学後の成績（例
えば教養課程の GPA）と学力試験成績，および内申書成績の相関を検討すると，しばしば
内申書成績との相関が高くなる．これは入学者は学力試験成績が基準以上のものに限定され
ていることの影響があるためである．一方，内申書は直接には合否判定に用いられないため
に，学力試験成績よりは広い範囲のものが含まれる．このため，不合格者までを考慮にいれ
ると（不合格者に大学教育を実際には行えないが），かならずしも高校内申書と学内成績と
の相関が，学力試験と学内成績の相関よりも大きくなるとは限らない（図 1）.
実験においては，条件の割り当ては実験者によって制御され，この割り当ての方法が適切なら他
の要因の影響を受ける可能性は小さい．先の自動車の停止距離の例で考えると，速度は停止距離に
影響を及ぼす要因と関係を持つことがないように，無作為化または事前の計画を用いて設定するこ
とが可能である．特に無作為化を用いる場合には，たとえ未知の要因が存在したとしても，それと
説明変数の関係が断ち切られるので，回帰係数の推定に偏りは生じない．この場合には，回帰係数
を説明変数から被説明変数への影響を表すのものとして解釈することは妥当である．
しかし，社会調査などの実験者による介入を伴わない観察データにおいては，説明変数と被説明
変数の両者に影響を及ぼしている変数は多く考えうる．これらのうちの幾つかは測定されたデータ
の内にあり，説明変数に加えることにより影響を除くことが可能かも知れない．しかしながら，考
4
2
1
0
-2
-1
内申書
-2
-1
0
1
2
学力試験
図 1: 学力試験得点と内申書得点の同時分布図（概念図）
慮すべき変数を全て列挙しているという保証が存在しない．分析者が，主要な変数を全て測定した
と考えたとしても，その他に説明変数と被説明変数の両者に影響を及ぼしている変数を見落してい
る可能性を否定できない．このため，観察データへ回帰分析を適用して得られる結果を影響関係の
推定とみなすことは，変数間の影響構造についての明確な知識が予めなければ困難になる．
1.4
1.4.1
線形モデルとしての分散分析
１元配置の分散分析
初等統計の教科書では，回帰分析と分散分析は別の手法として，取り扱われることが多いが，X
を群への所属を表すダミー変数とすることにより，共通のモデルとして取り扱うことができる．
ここでは，分散分析 (Analysis of Variance, ANOVA) のうち，最も簡単な場合である１元配置
(one-way design) について説明する．回帰分析において説明変数は連続な値をもつもの (身長，
血圧など) であったが，分散分析では離散的な区分が説明変数である．ここで，３種類の教育方法に
よって，生徒の成績に違いがあるか否かを検討することを想定する．３つのクラス j = 1, 2, 3 がそ
れぞれ異なる教育方法の対象となっている．クラス j の生徒の成績を yij , (i = 1, , , , Nj ; j = 1, 2, 3)
とする．また N1 + N2 + N3 = N とする．分散分析のモデルは，各クラス毎にデータが独立に同
一の正規分布に従っていることを仮定する．つまり
yij = µj + εij , (i = 1, ..., Nj ; j = 1, 2, 3)
(17)
であり，εij ∼ N (0, σ 2 ) である．検討すべき課題は µj , (j = 1, 2, 3) が互いに同一か否かである．
5
ここで，次のような行列を X として想定する．










X=









1
..
.
1
1
..
.
1
1
..
.
1
0
..
.
0
1
..
.
1
0
..
.
0
0
..
.
0
0
..
.
0
1
..
.
1




















(18)
最初の列のみが 1 である行は N1 ，1 列目と 2 列目が 1 である行は N2 ，また 1 列目と 3 列目が 1 で
ある行は N3 回繰り返されているものとする．この X を用い，β = (µ1 , µ2 − µ1 , µ3 − µ1 )T と置
くことにより，分散分析のモデル (17) は，
y = Xβ + ε
(19)
と表わされる．ここで，y は y11 , ..., yN1 1 , y12 , ..., yN2 2 , y13 , ..., yN3 3 と yij の値を縦につないだベク
トルである．平均に差があるか否かは，仮説 H0 : β2 = µ2 − µ1 = 0 および β3 = µ3 − µ1 = 0 を検
定すればよい．
上に示した (18) の第 2 列と第 3 列目は，それぞれの行に対応する yij の値が，第 2 群または第 3
群に所属する場合に各々1 となり，それ以外では 0 となる．このように，離散的なカテゴリーへの
所属を 0, 1 で表す変数をダミー変数 (dummy variable) と呼ぶ．
ここで，１元配置の分散分析モデルについて，(18) とは異なる表現が可能であることを示す．(18)
においては，X の第 1 列目は定数項 (すべて 1) を表し，第 2 列目と第 3 列目は各々第 2 グループ
(j = 2) と第 3 グループ (j = 3) への所属を表すダミー変数であった．このようにして構成される
X を X a とし対応する回帰係数ベクトルを β a = (βa0 , βa2 , βa3 )T とする．一般的に，分散分析な
どの離散的な説明変数を含む回帰モデル (一般線形モデル) において，説明変数およびそれらのダ
ミー変数から構成される X のことをデザイン行列 (design matrix) と呼ぶ．ここで示そうとし
ていることは，実質上同等なモデルが複数のデザイン行列によって表すことが可能なことである．
次のような 2 つのデザイン行列を想定しよう．
1. X の第 1 列は第 1 グループへの所属を表すダミー変数とし，第 2 列，第 3 列を X a と同じも
のとし，これを X b と表す．回帰係数ベクトルは β b = (βb1 , βb2 , βb3 )T である．
2. X の第 1 列は定数項，第 2 列は第 1 グループへの所属を表すダミー変数，第 3 列は第 2 グ
ループへの所属を表すダミー変数，さらに第 4 列は第 3 グループへの所属を表すダミー変数
とする．これを X c とする．X c は X b の左側に定数ベクトルを並べたもの，つまり (1|X b )
である．但し，X c に対応する回帰係数のベクトルを β c = (βc0 , βc1 , βc2 , βc3 )T として，次の
制約を満たすものとする．
βc1 + βc2 + βc3 = 0
(20)
ここで，X a によるモデルの表現は，X c をデザイン行列とし，(20) の代わりに βc1 = 0 と制約を
置いたものである．
6
これら 3 つのモデルはいずれも同等なものである．つまり，観測値 y に対して，推定値 yˆa = X a βˆa
が得られるとすると，つねに yˆa = X b βˆb = X c βˆc となる βˆb と βˆc が一通りに決まる．これら 3
∑
つのモデルにおいて，推定される β の値は異なるが，推定値と残差 2 乗和 (yi − yˆi )2 は同一で
ある．ソフトウェアによっては，互いに異なるデザイン行列の表現を利用している場合がある．
1.4.2
帰無仮説の表現
回帰分析においては，
「幾つかの回帰係数がゼロであること」を帰無仮説として取り上げた．分
散分析のように，ダミー変数を説明変数として持つモデルにおいては，より多様な形式の帰無仮説
を取り扱う必要が多い．次の ANOVA モデルを考える．
yij = µj + εij , (i = 1, ..., Nj ; j = 1, 2, 3)
(21)
ここで，通常もっともよく利用される帰無仮説は，µ1 = µ2 = µ3 というものである．デザイン行
列 (13) が設定されているとすると，回帰係数は β1 = µ1 , β2 = µ2 − µ1 , β3 = µ3 − µ1 であり，上の
帰無仮説は H0 : β2 = β3 = 0 として表される．しかし，一般に検定の対象となる帰無仮説は，上
のタイプのものには限らない．例えば，仮説 µ3 − µ2 = 0 は，上のデザイン行列では，β3 − β2 = 0
と表現され，ある i について βi = 0 という形式では表現されない．
上にのべた幾つかの帰無仮説はいずれも，回帰係数の (１つまたは複数の) １次式として表され
る．例えば，β2 = β3 = 0 という仮説は
(
T
Q =
0 1
0 0
0
1
)
(22)
とおくと
QT β = 0
(23)
と表記される．一般的に，QT β = ξ のような回帰係数への線形の制約式で表される帰無仮説を，
線形仮説 (linear hypothesis) と呼ぶ．同じデザイン行列のもとで，仮説 β2 − β3 = 0 は
(
)
QT = 0 1 −1
(24)
とおくと，QT β = 0 として表される．ここで，デザイン行列 X N ×p の各列ベクトルが１次独立で
あれば，QT がどんなものであろうと QT β の値は，データから一意的に推定できる．ところが，X
の列ベクトルが１次従属の場合，つまり rankX < p の場合には，一つのデータ y について QT β
の値を一意的に推定できない場合が生じ得る．これを避けるためには，QT の各行ベクトルが X
の行ベクトルの１次結合として表される必要がある．この条件は，
(
rank
X
QT
)
= rankX
(25)
と表現される．これはまた，KerX ⊆ KerQT とも表される．ここで，KerA とは Ax = 0 となる
x 全体のことをあらわす．
この条件が成り立つことを線形仮説が推定可能 (estimable) であるという．ある線形仮説が推
定可能でないということは，指定されたデザイン行列がその仮説を検定するために適切ではないこ
とを意味する．
7
ここでは，当面 X のランクは p であると仮定する．つまり X の各列は 1 次独立であるものと
する．線形仮説 H0 : QT β = ξ を検定するためには，制約下での推定値と制約無しの下での推定値
ˆ 0 とし，また
を求め，各々について残差の２乗和を求める．帰無仮説の制約下での Y の推定値を y
ˆ 1 とする．Q のランクを k とすると，もし帰無仮説が正しいのなら，
制約無しの推定値を y
ˆ0
y−y
2
ˆ1
− y−y
2
∼ σ 2 χ2 (k)
(26)
である．またこれとは独立な分布として
ˆ1
y−y
2
∼ σ 2 χ2 (N − p)
(27)
が得られる．回帰分析における変数選択の場合と同様に，次の F 値が帰無仮説のもとで自由度
(k, N − p) の F 分布に従うことを用いて検定を行なえる．
F =
1.5
ˆ0 2 − y − y
ˆ1
y−y
ˆ1 2
y−y
2
×
N −p
k
(28)
共分散分析
以上で説明したように，回帰分析と分散分析とは数学上は同一モデルとして表され，回帰係数の
推定や検定もほぼ同一の手順によって行うことができる．これらの手法が異なるのは，モデルの特
徴というよりは，むしろ応用の対象となるデータの性質であることが多い．ここで説明したモデル
は，説明変数が実験者によって設定された値であることを前提としている．分散分析が応用される
場面ではこの前提が当てはまることが多いが，実際に回帰分析の応用，特に社会科学や疫学ではこ
の前提は多くの場合当てはまらない．説明変数が確率的に変動する場合にも，ある種の条件を仮定
すれば，予測値の推定は説明変数が実験者によって設定された場合と同様の手続きで行うことが可
能である．しかし最大の問題は，回帰係数の解釈が不明確になることにある．
上述の線形のモデルにおいて，説明変数に離散値と連続値の両者を含むものは共分散分析 (Anal-
ysis of Covariance, ANCOVA) と呼ばれる．以下に共分散分析の例を示すが，この例におい
てはデータは疫学調査に基づくものであり，研究者によって設定された値ではない．このような場
合，回帰係数が影響の程度を表しているとみなせるためには，説明変数と非説明変数の両者に影響
を及ぼす変数が全てモデルに取り込まれている必要がある．
共分散分析モデルには，次のような特徴がある．
1. 離散値をもつ説明変数 (要因) は X のいくつかの列におけるダミー変数として表現される．
また，これらの交互作用もダミー変数となる．
2. 共分散分析においては，さらに連続値を持つ説明変数と離散値を持つ説明変数の交互作用も
用いられる．これらの交互作用を表す変数は，離散値を表すダミー変数 (0-1 値) に，連続変
数の値をかけたもの，つまりダミー変数が 1 の部分のみ連続変数の値を残し，他の部分はゼ
ロとなる変数によって表される．
ここで X1b と X1c がある要因 (3 要因) を表すダミー変数であるとし，X2 は別の連続値を持
つ説明変数であるとする．これらの各要因に対応する標本を A 群，B 群，C 群と名付ける．
X1b は A 群と B 群における差に対応し，X1c は A 群と C 群との差に対応する．
主効果のみのモデルは，次の式で表される．
E(Y |X) = β0 + β1b X1b + β1c X1c + β2 X2
8
表 1: 出生時体重と母親の喫煙データ (Nolan & Speed, 2000 より)
変数
変数名
内容
出生時体重
bwt
新生児の出生時体重 (オンス = 28.3g)
妊娠期間
gestation
parity
妊娠期間 (日)
母親の受胎時の年齢 (年)
体重
age
height
weight
喫煙
smoke
母親の喫煙 (Smoke)，非喫煙 (Nosmoke)
パリティ
年齢
身長
第 1 子 (First) か否 (Other) か
母親の身長 (インチ)
母親の妊娠前の体重 (ポンド)
ここで，β1b と β1c とは，離散変数であらわされる群間の Y の違いを表し，β2 は各群に共通
な連続変数 X2 と Y との関係を表す．
離散変数と連続変数の交互作用を含むモデルは，
E(Y |X) = β0 + β1b X1b + β1c X1c + β2 X2 + β1b2 X1b X2 + β1c2 X1c X2
となる．ここで，β1b2 が非ゼロであることは，A 群と B 群において X2 が Y に及ぼす効果が
異なることを示す．A 群においては X2 の回帰係数は β2 であり，B 群においては β2 + β1b2
である．β1c2 についても同様のことが当てはまる．
3. 連続変数間の交互作用は，これらを単純に掛算した 2 次式によっては極めて限られた形の効
果しか表現できない．スプライン関数とよばれる区分的多項式 (区分的 3 次式がよく用いら
れる) を用いると，連続変数間の交互作用を柔軟に表現できるが，現在のところ利用可能な
ソフトウェアは限られているが，一般の回帰分析プログラムを工夫して利用することによっ
て分析をおこなうことは可能である．
1.5.1
共分散分析の例
Nolan& Speed (2000) で引用されている新生児の体重と妊娠中の母親の喫煙に関するデータを
分析対象とする．このデータは，1960 年から 1967 年にサンフランシスコにおける健康調査によっ
て得られた (Yerushalmy,1971)．データは 1236 件であるが分析には欠測値のないもの 1174 件を用
いた．1
ここで示されている変数は次の 7 つである．
データの一部を次に示す．
bwt gestation parity age height weight
smoke
1 120
284 First 27
62
100 Nosmoke
2 113
282 First 33
64
135 Nosmoke
3 128
279 First 28
64
115
Smoke
4 123
NA First 36
69
190 Nosmoke
5 108
282 First 23
67
125
Smoke
...
1 データは
http://www.stat.berkeley.edu/users/statlabs/ において公開されている．
9
以下の分析は R-1.5.0 による．分散分析表は説明変数を逐次加えた場合の，残差 2 乗和の減少分
を表している．標本件数のバランスした分散分析とは異なり，残差 2 乗和の減少量は変数の追加順
に依存する．
モデル 1: 5 つの説明変数を用いたモデル．年齢は影響が小さいため説明変数から除いた．
Residuals:
Min
1Q
-57.7164 -10.1500
Median
-0.1594
3Q
9.6885
Max
51.6199
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -80.71321
14.04465 -5.747 1.16e-08 ***
gestation
0.44408
0.02907 15.276 < 2e-16 ***
parityOther -3.28762
1.06281 -3.093 0.00203 **
height
1.15497
0.20473
5.641 2.11e-08 ***
weight
0.04983
0.02503
1.991 0.04672 *
smokeSmoke
-8.39390
0.95117 -8.825 < 2e-16 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.82 on 1168 degrees of freedom
Multiple R-Squared: 0.2579,Adjusted R-squared: 0.2548
F-statistic: 81.2 on 5 and 1168 DF, p-value: < 2.2e-16
AIC は 9823.503
分散分析表
Response: bwt
Df Sum Sq Mean Sq F value
Pr(>F)
gestation
1 65450
65450 261.4292 < 2.2e-16 ***
parity
1
2345
2345
9.3658 0.002261 **
height
1 12554
12554 50.1437 2.462e-12 ***
weight
1
1801
1801
7.1941 0.007418 **
smoke
1 19497
19497 77.8779 < 2.2e-16 ***
Residuals 1168 292412
250
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
10
モデル 2: 年齢以外の 5 つの変数，および妊娠期間と喫煙の交互作用項を用いたモデル．自由度調
整済 R2 は増加，AIC は減少している．いずれもモデルがより適切であることを示唆している．
ただし，実際の残差を検討してみると，非喫煙群で妊娠期間が極端なケース (短い場合と長い場
合) があり，検討が必要と思われる．(図.4)．
Residuals:
Min
1Q
-57.87631 -9.95257
Median
-0.01195
3Q
9.54773
Max
53.49886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)
-59.07260
15.30129 -3.861 0.000119
gestation
0.37020
0.03590 10.313 < 2e-16
parityOther
-3.29621
1.05780 -3.116 0.001877
height
1.14888
0.20377
5.638 2.16e-08
weight
0.04542
0.02494
1.821 0.068855
smokeSmoke
-66.82659
16.83114 -3.970 7.61e-05
gestation:smokeSmoke
0.20970
0.06031
3.477 0.000525
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘
***
***
**
***
.
***
***
’ 1
Residual standard error: 15.75 on 1167 degrees of freedom
Multiple R-Squared: 0.2656,Adjusted R-squared: 0.2618
F-statistic: 70.33 on 6 and 1167 DF, p-value: < 2.2e-16
AIC は 9813.402
分散分析表
Response: bwt
Df Sum Sq Mean Sq F value
Pr(>F)
gestation
1 65450
65450 263.9117 < 2.2e-16 ***
parity
1
2345
2345
9.4547 0.0021550 **
height
1 12554
12554 50.6198 1.952e-12 ***
weight
1
1801
1801
7.2624 0.0071424 **
smoke
1 19497
19497 78.6175 < 2.2e-16 ***
gestation:smoke
1
2999
2999 12.0910 0.0005253 ***
Residuals
1167 289413
248
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
11
120
|
100
|
60
80
babies.nona$bwt
140
160
180
Observed |:Nosmoke, −:Smoke
|
150
200
| | | ||
|
| ||
||||−− |−−
|
|
|
|
||| || −
|−
||| −
− | |−−−
|| |
−| |−|−
|||||||||| |
|
|
| | | | −| | −
−
−
|
|
|
−
−
|
|
|
|
|−
| |−
| | |−−
|−
−−
||−
||−
|−||||| − −−
|||−
||−
−
||−−
|||||−
||−
||||||−
||−
|−
|||||−
|−
|||−
|||| − || |
| || |||| − −
|
−
|
|
|−
|
|
|
|−
||−−
| | || ||−|−|||−
|
||−
−
−
|
|
|
|
−
−
|−−
| −|
−
|−
|−
|
|−−
|
−
−||−
|−
| |||−
|−
|||||−
|−
||||| −
||−
||−
||−
| −
||−
|||−
|||−
||−
|| ||−
|||−
||−
||−
| || |− −| −
|−
||||−
|−
|−
|−||||||−
|
−
−
−−
||−
||||||−
|−
||−
|−
|−
|−−
||||−
−
|||−
|−
| || − || |
−
−
−
−| |||−−−
|
|
|
|
||−
−
−
|
|
|
|
|
|
|
|
|−
|
−
−
−
|−
|−||−
||−
|−
|−
|||−−−
|−
||
|−−
|| −|−−−
|−
| ||−
|−
|−
|||||−
|−
|||−
−
−
−
−−
−
||−−
|−
||−
|−
|||−
|−
||||−
| ||−−
||−−
|−
−−
−
||−
| |−
|−−
||−
|||−−
||−
|−
||−
−−
|−
||−
| |−−−|||−
| | ||
|||−
||−
||−
−
|
|
|
|
|
−
−
|
|
|
−
−
−
−
−
−
−
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
−| −| −−−
−
−−−
− ||−−| |||||
−−
−
−
−
|
|
|
|
−
|
− | |−
−
|−
−
−
|−−
|−
|||−−−
|−
| | −| || −
|−−−
|−−−|−
| −||−
−
−
|
|−
||−−|−||||−−
−
−|
||−
||−−
|−
−− | | −|| −||−−−−−
|||−
|−
−|
||−
|−
−
−
−
−
−
|−
|−
−|−−
| − ||−| −−
||−
−
−
−
−
−
|
−
−
−
||−
−
−||| − −−|−−−
|
−−
−| −−|−
|−|−−− |−|
−−
| | | −−−−−−| |−− −−−
−|− | |−−| −| −
|
|−−−−||
| −− −
− − | || −
−
−
−
−|
| |
−
250
300
|
|
350
babies.nona$gestation
図 2: 妊娠期間と出生時体重 (観測データ)
120
100
80
|
|
| ||||| −
| | || | −−
|
−
|
|
|
|
−
|
|
|
|
| ||| | ||||||||||||−
||| ||||| | − −
||||||
||||−
|||||−
|||||||||||||||||||||||||||−
|−|||−
| −−
|−
||||−
|−
|||||||||||−
−
|||−
|||||||||||||−
|−
|||−
|||−
|||−
−
−
||||||||−
||||−
|
|
−
|
|
|
|
|
|
|
|
|
−
|
|
|
|
|
−
−
|
|
−
−
||−
|−
|−−| −−
||||||−
||−
|−
−
−
−
−
|−
|
−
|−
−
|
−
−
|
−
|||−
|
|
||||||||||||||||||−
−
−
|
−
|
|
|
−
|
−
|
|
|
−
|
−
|
|
−
−
−
||−
−−−−−−− −
−
||−
||−
||−
|||−
−
||−
|−
|−
−
|−
−
|||−
||−
−
||−
−
−
−−
|−
−−
||−
−
−
−
| |||||||| ||||||||||−
|
−
|
−
|
−
|
−
|
−
|
−
−
−
−
|
|
−
−
−
|
−
|
|
|
−
−
−
|
−
−
−
|
−
−
−
−
−
−
−
|−−−−−
−
||−
−−
|−
|−
−−
−
−
|−
|−
||||||||||||||−
|−
−
−−
|−
−
||−
|−
−−
||−
|−
−−
−
|−
−−
||−
−
−
|−
−
−
|−
|−
|−
||−
−
−
|−
|||| |||||−
−
−
|−
||−
−−−−−−
−
−−
−
−
−
−−
−
|−
−−
|−
−
−
−
−−−
−
−−
−−−
||| | ||| ||||−−
−
−
−
−−
−−−−−−−
−
−
−
−
−
| | |−−|−
−
−−
−−
| −−
−−−−−−−−−−
−
|−||| −
−|−
−−| −|−−−−− − −
− −
−
|
|
|
|
|
60
babies.nona$bwt
140
160
180
Estimated (parallel model) |:Nosmoke −:Smoke
150
200
250
300
350
babies.nona$gestation
図 3: 妊娠期間と出生時体重 (モデル 1 の予測値)
12
120
100
|
80
|
|
− |
| −−
− |−|−
||| −
|
−
|
||| |−| −| ||| |
| |||||||||||||||||−
|| −−|
−
| ||||||||−
−
|||||||||−
|| |−
|−
||||−
−
||−
|||||||−
|||−
||||||||||||||||||||||||−
|||||||||−
|−|||||−
|||−
||−
||
|
−
−
|
|
|
|
|
−
|
−
|
|
|
|
|
−
|
|
|
−
−
|−
−
|−−
−
||−
||−
|−
|
−
−
||||||||−
|
|
−
|||−
|
−
|
||−
|
−
−
|
|||||||||||||||||||||||||||−
−
−
−
−
|
|
|
|
|
−
|
|−
|−−|
−
−
|
−
|
|
|
−
|
|
−
|
−
−
−
−
|−
|
−
|
||−
|
||−
−
−
|
|
|
−
|
|
|
||||−
−
|||−
||−
−
−
|
−
| ||||||| ||||||||−
|
−
||−
−−−
−−
|−
||−
−
−
|||−
−−
||−
−
||−
−
|−
|−
−
||−
−
−−
|−− −
−
|−
||−
−
|
|
−
|||||||||||||||||||||||−
|
|
−
−
−
−
−
−
|
−
−
|
−−
−
−
−
−
−
|
|
−
−
|
|
−
−
−
−
−
−
|−
|−
−
|−
|−
|−
−−
−−
−
||−
|−
−
−
| | |||| |||||| ||||−
−
−−
−−−
−
−−−−−
|−
−
−
||−
|−
−−
−−
−
−
−
−−−
| || | | || ||−−|−−−
−−−
−−−
−−
−
−−
−−
−−
−−−−
−
−−
−
−−−−
−
−
| ||| | | −−−−−−−
−
−
−−
|
−−−−−
− −−−−
−−
|−−−−−
−− −− −
−
−
−
−
−
−
−
−− −− −
− −
|
|
60
babies.nona$bwt
140
160
180
Estimated |:Nosmoke −:Smoke
150
200
250
300
350
babies.nona$gestation
図 4: 妊娠期間と出生時体重 (モデル 2 の予測値)
13
２項分布，多項分布とポアソン分布
2
回帰分析や分散分析，共分散分析では従属変数は連続な値をとり，独立変数によって説明される
平均構造と，確率的に変動し正規分布に従う誤差の和であると仮定していた．以下では従属変数が
離散値である場合の分析（多重分割表）について紹介するが，このような場合に用いるモデルで
は，離散値をとる確率分布によって従属変数を表す必要がある．
離散値をとる確率分布として代表的なものには，２項分布，多項分布，ポアソン (Poisson) 分布，
超幾何分布，負の２項分布などがあるが，ここでは最も基本的な２項分布，多項分布とポアソン分
布について説明を加える．また，超幾何分布については，Fisher の正確検定の説明の箇所で説明
する．
2.1
2 項分布 (binomial distribution)
確率変数 X が 2 項分布 Bin(1, p) に従うとは，X = 1 の生じる確率が p であり，X = 0 の生じる
確率が 1 − p であることを意味する．コイン投げで表を 1 とし，裏を 0 と対応させると，Bin(1, p)
は表が p の確率で生じるコインのモデルになる．また，Y が 2 項分布 Bin(n, p) に従うとは，それ
ぞれ独立に Bin(1, p) に従う Xi , (i = 1, ..., n) があり，Y = X1 + · · · + Xn であることを意味する．
このとき Y の確率は
(
Pr(Y = k) =
n
k
)
pk (1 − p)n−k , (k = 0, 1, ..., n)
(29)
となる．繰り返しの数 n が大きくなると, 平均ゼロ，分散１に標準化した Y の分布は次第に標準正
規分布に近づく．確率変数 Y ∼ Bin(n, p) の平均は np，分散は np(1 − p) である．
2.2
多項分布 (multinomial distribution)
2 項分布 Bin(1, p) は X が 1 か 0 のいずれかをとる場合を考えているが，ここで X を長さ m の
ベクトルとし，m 個の要素の内いずれか一つのみが 1 となり，他の要素は 0 であるとする．また，
∑m
それぞれの要素が 1 となる確率を pj とする．ここで j=1 pj = 1 である．これを 1 回の試行とし，
∑n
n 回の独立な試行を行って得られる n 個の m 次元ベクトルを X1 , ..., Xn とし Y = i=1 Xi とす
る．この m 次元ベクトル Y の分布を多項分布とよび Mn(n, (pj )) で表す．
多項分布の確率は
Pr(Y = (y1 , ..., ym )) =
n!
py1 · · · pymm
y1 ! · · · ym ! 1
∑m
(ただし j=1 yj = n) で表される．ベクトル Y の要素 Yj の平均と分散は 2 項分布の場合と同様
にそれぞれ npj および npj (1 − pj ) である．また，Yj と Yk (j = k) の共分散は −npj pk である．
単純な多項分布の例は，サイコロやルーレットである．サイコロの場合には m = 6 であり p1 =
· · · = p6 = 1/6 である（イカサマがなければ）．サイコロを n 回投げて，その結果得られた各々の
目の頻度 (長さ 6 のベクトル）が Y である．
14
2.3
ポアソン分布 (Poisson distribution)
2 項分布 Bin(n, p) に従う確率変数 Y を考える．Y の平均は np であるが，この値を一定に保ち
つつ n を大きくする．np = λ とおくと，p = λ/n である．(29) より，Pr(X = k) は
)n
(
1 − nλ
n × (n − 1) × · · · × (n − k + 1) λk
× k ×(
(30)
)k
k!
n
1− λ
n
であるが，ここで n → +∞ とすると，
Pr(X = k) →
λk
exp(−λ)
k!
(31)
となる．この分布はポアソン (Poisson) 分布と呼ばれるものであり，Po(λ) と表す．
ポアソン分布は，互いに独立な小さい平均を持つ 2 項分布の多数の和を近似していると見なせる
ため，件数 (頻度) の分布のモデルとしてしばしば用いられる．Po(λ) の平均は λ，分散も同じく λ
である．
Poisson lambda=3
0.15
Pr(X=x)
0.05
0.00
0.00
0
2
4
6
8
10
0
X
5
10
15
X
図 5: ２項分布の例
2.4
0.10
0.15
0.10
0.05
Pr(X=x)
0.20
0.20
0.25
Binomial n=10 p=0.6
図 6: ポアソン分布の例
過大分散と過小分散
連続変数の誤差を表す正規分布では，平均が定まっても異なる分散を持つ場合があり，これは分
散パラメータの違いによって表現される．一方，２項分布や多項分布では，平均 (pi ) が定まれば，
必然的に分散も決定される．ポアソン分布についても，類似の性質が成り立つ．
もし，データが同一の母数を持つ２項分布（あるいは多項分布，ポアソン分布）から得られてい
るとの仮定が正しいなら，これは正しく成立するが，実際に得られるデータでは，複数のサンプル
に対応する真の母数が，必ずしも同一ではない場合がありうる．正規分布の場合には，このよう
な母数のバラツキも Y の分散によって吸収することができるので，見かけ上誤差の分散が大きく
なったものとしてパラメータを推定する．しかしながら，被説明変数に２項分布などを用いる場合
には，このような分散の変化を表す自由度がない．これらの現象は過大分散 (over dispersion) と
呼ばれる2 ．離散値の被説明変数を想定する推定手法においては，このようなケースに対応するた
2 場合によっては Y の分散が２項分布から期待されるものよりも小さくなる場合もあり，この場合は過小分散 (under
dispersion) と呼ばれる
15
めの分析手法も開発されている (SAS の GENMOD においては，過大分散パラータの推定を行う
ことができる）.
２×２のクロス表, オッズ比
3
3.1
２×２表の関連性
離散的な値（合格・不合格，性別，現役・既卒など）を持つ変数間の関係を調べる基本的な方法
は，これらの変数間についての分割表を作成することである．ただし，適切な層別を行わないと，
不適切な見かけ上の関係が得られる場合がある．
表 2 は 1973 年におけるカリフォルニア大学バークレイ校の５つの学科での入学志願者数と合
否の集計結果である．このデータは一見すると不可解な情報を与える例としてよく知られている．
Bickel et al. (1975) によれば，1973 年のバークレイ校での大学院の学科またはコースは，全体で
101 あり，志願者の総数は 12,763 名である．Bickel らは，入学許可に男女差別があるのではない
かとの疑問に答える形での分析を行い，単純な分析が誤った結論を導く可能性を指摘している．
ここに示したのは，この入学にかかわるデータのうち，Freedman et al.(1998) に示されている６
つの学科についてのものである3 ．表 2 に示した性別の集計結果を見ると，男子の合格率は 44.5%で
あるのに対し，女子の合格率は 30.4%であり，明らかに女子の合格率が低い．女子が大学院に合格
しづらい傾向があるのだろうか？
以下では，まず２×２の分割表の関係を調べるための統計的な指標とそれらの性質について検討
し，ついで学科別に層別された分割表の検討を行う．
表 2: 1973 年におけるカリフォルニア大学バークレイ校大学院６学科への入学志願と合否の結果
（性別）
実数
率（%）
合格
不合格
合計
合格
不合格
合計
女子
1198
557
1493
1278
2691
1835
44.5
30.4
55.5
69.6
100.0
100.0
合計
1755
2771
4526
38.8
61.2
100.0
男子
Source: Freedman et al. (1998) による．これらは人数の多い６つの学科についてのもの．大学院全体の受
験者数は 12,763 である (参考, Bickel et al., 1975 ).
3.2
２×２分割表の関連性指標
2 つの変数 X1 と X2 とがそれぞれ 2 値 (0 と 1，または Yes と No など) の値をとるものとする．
これらの関係を示すためにいくつかの方法が用いられている．
ふたつの 2 値変数 X1 と X2 のクロス集計されたデータの件数を 2 × 2 の行列 nij , (i = 1, 2; j =
1, 2) で表すものとする．行は X1 に対応し，列が X2 に対応する．また，ni+ によって第 i 行の
和を，n+j によって第 j 列の和を表す．標本の総件数は N とする．また，クロス表の相対頻度を
pij , (i = 1, 2; j = 1, 2) で表し，その周辺相対頻度を pi+ および p+j で表す．
3 このデータはフリーの統計分析ソフトウェアであるＲシステム
16
(Ihaka & Gentleman, 1996) の例題にも含まれている．
より具体的には，p1+ は p11 + p12 であり，p2+ = p21 + p22 である．ここでの例では，それぞれ
男子の志願者とと女子の志願者の相対頻度である．また p+1 = p11 + p21 であり，p+2 = p12 + p22
である．ここではそれぞれ合格者と不合格者の相対頻度となる．
実数
相対頻度（率）
X2 = 1
X2 = 2
行計
合格
不合格
行計
X1 = 1
X1 = 2
n11
n21
n12
n22
n1+
n2+
p11
p21
p12
p22
p1+
p2+
列計
n+1
n+2
N = n++
p+1
p+2
1.0
ここで，行 X1 と列 X2 の間に統計的な関係がないことを，X1 と X2 とが統計的に独立であると
いう．これが成立するのは
p21
p11
=
p12
p22
である場合である．表 2 の例においては，男子と女子で合格率が同じであることを意味する．
２×２表の行と列との関係を示すための指標としては，次のものが用いられるている．
1. オッズ比 (odds ratio)
α=
p11 p22
，
p12 p21
標本から計算されるオッズ比は，α
ˆ = (n11 n22 )/(n12 n21 ) となる．頻度がゼロのセルがあると推
定上の問題が生じる．このような場合には α
ˆ = ((n11 +0.5)(n22 +0.5))/((n12 +0.5)(n21 +0.5))
のように修正した値がしばしば用いられる．
合格率を「賭け率」（オッズ）とみなすと，上の値 α は二つの条件における賭け率の比とな
る．このためこの値 α はオッズ比と呼ばれる．オッズ比は行と列とが独立なら１である．最
小はゼロであり，上限はない．
2. 対数オッズ比．
log α = log
p11 p22
．
p12 p21
オッズ比の対数をとったもの．オッズ比自体よりも扱いやすい性質があるためしばしば利用
される．行と列が独立なら（オッズ比１であり）ゼロとなる．
対数オッズ比の分散は各 nij が大きい場合に次の式で近似される．
[
]
n11 n22
1
1
1
1
Var log
+
+
+
n12 n21
n11
n22
n12
n21
3. 相対リスク (p11 /p21 ). これは疫学での用語である．性別を暴露条件，合格を発病とみなすと
p11 /p21 は，2 つの条件間での発病確率の比になる．もし，発病がまれであるなら，p12 /p22
はほぼ 1 であるので，相対リスクはオッズ比とほとんど同じになる．
4. 4 分点相関係数 (ファイ φ 係数)．2 値データが０または１値をとるものとみなして Spearman
の積率相関係数（一般的な相関係数）を求めると次の式が得られる．
p11 p22 − p12 p21
φ= √
p1+ p2+ p+1 p+2
この値は通常の相関係数と同じ手順で求められるので，便利ではあるが，後述するように周
辺度数の偏りの影響を受けやすい．この値は-1 から 1 の範囲をとる．
17
独立性の検定に用いられる X 2 統計量
X2 =
∑ (nij − N pi+ p+j )2
N pi+ p+j
i,j
は φ2 の N 倍になる．
5. ４分相関係数 (tetrachoric correlation)．
分析対象となる 2 × 2 表の背後に，相関のある 2 変量正規分布を仮定し，2 つの変数がある閾
値によって分割されることにより，4 つに区分されたデータが得られていると解釈する．図 7
は２変量正規分布に基づくデータの例であるが，ここで縦線と横線の位置および分布の相関
を調節することにより，分布から期待される各領域の確率と得られたデータとが適合するよ
うにする．この仮定にもとづいて，データにもっとも良く当てはまる 2 変量正規分布を推定
し，その相関係数を求める．このようにして得られる値を 4 分相関係数と呼ぶ．この値は簡
単な式で求めることができないので，近似式による推定か反復を伴う数値計算が必要になる．
−3
−2
−1
y
0
1
2
3
r=0.3, N=1000
−3
−2
−1
0
1
2
3
x
図 7: ２変量正規分布に基づくデータ例 (母相関は 0.3)
これらの指標にはそれぞれの特徴がある．4 分点相関係数は周辺度数の偏りによって大きく影響
を受けるため，相対頻度が小さいカテゴリーがある場合には適切ではない．
ここで，表 2 について，これらの指標の値を求めてみると次のようになる．
1. オッズ比
1198 × 1278
= 1.84
557 × 1493
値の範囲はゼロから無限大まで．性別と合格率が独立なら１．
18
表 3: 合格者を４分の１にした仮想例
実数
率（%）
合格
不合格
合計
合格
不合格
合計
女子
299
139
1493
1278
1792
1417
16.7
9.8
83.3
90.2
100.0
100.0
合計
438
2771
3209
13.6
86.4
100.0
男子
2. 対数オッズ比
log
1198 × 1278
= 0.610
557 × 1493
求められた値の標準偏差は約 0.064 と推定される. 値の範囲は，負の無限大から正の無限大
までどの値も取りうる．性別と合格率が独立ならゼロ．
3. 相対リスク
1198/2691
0.445
=
= 1.47
557/1835
0.304
4. ４分点相関係数（ファイ係数 φ）
1198 × 1278 − 557 × 1493
φ= √
= 0.1427
(1198 + 557)(1493 + 1278)(1198 + 1493)(557 + 1278)
−１から＋１の範囲の値をとる．性別と合格率が独立ならゼロ．この例では，より正確には
φ = 0.1427318．
ファイ係数から導かれるカイ２乗統計量 X 2 は必ずゼロ以上の値をとる．この場合
X 2 = φ2 × N = 0.14272 × 4526 = 92.21
である．ゼロであるのは，性別と合格率とが独立である場合に限られる．このデータについ
て独立性のカイ 2 乗検定を行うと，上側確率（p 値）は極めてゼロに近い値となり，性別と
合格率が独立ではないことが確認される．
5. ４分相関係数
専用のソフトウェアを用いた反復計算により推定すると，次の値が得られる．
rtetra = 0.230
3.3
関連性指標の特徴
これらの指標にはそれぞれ個性があり，利用にあたって注意を払う必要がある．特に利用される
機会の多い指標は（対数）オッズ比とファイ係数であると思われるが，これらは相反する性格を
持っている．
ここで表 2 のデータに変更を加え合格者数を（ほぼ）４分の１に変更してみる．
表 3 について，指標を求めると表 4 のようになる．
19
表 4: 合格者が４分の１の場合の関連性の指標（仮想例）
1
2
3
4
5
指標
変更後
変更前
オッズ比
1.84
0.610
1.70
1.84
0.610
1.47
対数オッズ比
相対リスク
４分点相関係数（ファイ係数） 0.099
４分相関係数
0.200
0.1427
0.230
このように，オッズ比は，ある行（または列）が一定の比率で減少（増加）しても変化しないが，
他の指標は変化する．特定の列（または行）の比率が小さい場合（例えば合格率が極端に小さい場
合）には，ファイ係数 φ はゼロに近い値になる傾向がある．このため，周辺度数の異なる２×２表
についてファイ係数を求めると，それらの大小が関連性の強さではなく，周辺度数の偏りの程度を
表す場合がある．
表 5 は，2 変量正規分布を各変数の閾値で区分し，(pij , i = 1, 2; j = 1, 2) を求め，4 分点相関係
数 φ とオッズ比 α を計算した結果である．同一の ρ について，φ と α が逆の傾向を示している．φ
は，周辺度数が偏ると小さくなる傾向があるが，逆に α は大きな値を示している．
表 5: 相関係数とオッズ比
母相関
ρ
0.3
0.9
3.4
周辺確率
(p1+ , p+1 )
(0.05, 0.95)
(0.20, 0.80)
(0.20, 0.95)
(0.50, 0.50)
(0.50, 0.80)
(0.50, 0.95)
(0.80, 0.80)
(0.80, 0.95)
(0.95, 0.95)
(0.05, 0.95)
(0.20, 0.80)
(0.20, 0.95)
(0.50, 0.50)
(0.50, 0.80)
(0.50, 0.95)
(0.80, 0.80)
(0.80, 0.95)
(0.95, 0.95)
同時確率
(p11 )
0.0495
0.1809
0.1968
0.2985
0.4336
0.4870
0.6661
0.7704
0.9071
0.0500
0.2000
0.2000
0.4282
0.4981
0.5000
0.7499
0.7988
0.9319
4 分点相関係数
φ
0.043
0.131
0.078
0.194
0.168
0.110
0.163
0.119
0.098
0.053
0.250
0.115
0.713
0.490
0.229
0.687
0.445
0.618
オッズ比
α
5.94
2.78
3.81
2.19
2.38
3.00
2.46
2.99
3.52
極めて大
18398.0
極めて大
35.59
169.97
極めて大
44.85
219.70
90.32
対数オッズ比
log α
1.78
1.02
1.34
0.79
0.87
1.10
0.90
1.08
1.26
−
9.82
−
3.57
5.14
−
3.80
5.39
4.50
仮説からのずれの指標
いま，一般的に 2 × 2 のデータについて仮説が次のように表されるものとする．
H0 : E(nij ) = µij (i, j = 1, 2)
20
(32)
ここで，左辺は該当するセルの頻度の期待値であり，期待頻度 (expected frequency) と呼ばれ
る．また，観測されたデータは nij , i, j = 1, 2 であるとする．
このとき，仮説 H0 からのずれを検討するための指標として，次の 2 つが代表的である．
1. カイ２乗 (χ2 ) 統計量
X2 =
∑ (nij − µij )2
µij
(33)
2. ポアソン分布に基づく対数尤度比 (の 2 倍)
G2 = 2
∑[
(
nij log
nij
µij
)
]
− nij + µij
(34)
いずれの値の分布も仮説が真ならば，各セルのデータの件数が多い場合には I × J 表の場合，自
由度 IJ の χ2 分布で近似される．
また，仮説が頻度ではなく確率 πij を特定する場合には µij = N πij であり，これらの χ2 分布の
∑
∑
自由度は，(IJ − 1) となる．またこのとき (34) は N =
nij = µij が成立し，多項分布に基づ
く対数尤度比 (の 2 倍)
2
G =2
∑
(
nij log
nij
µij
)
=2
∑
(
nij log
nij /N
πij
)
(35)
となる．
この性質を用いて，特定の仮説が成立しているか否かを検定できる．
3.5
独立性の検定
X と Y (行と列) とが独立であることを検定するためには，上に類似の手続きを利用するが，少
し変更を要する．統計的に独立であるということは，各セルの値が特定の固定されたものであるこ
とを必要とするわけではない．課せられている制約はもう少し緩やかなものである．
∑ ∑
∑
∑
ここで標本の周辺度数を i nij = n+j ,
j nij = ni+ ,
i
j nij = n++ とする．また，周
∑ ∑
∑
∑
辺度数の期待値を i µij = µ+j ,
j µij ,
i
j µij = µ++ と表記する．
独立性の仮定は，2 × 2 の分割表においては
µij =
µi+ µ+j
(i, j = 1, 2)
µ++
(36)
が成り立つことに等しい．そこで，検定を行なうには，標本においてこの関係からのずれの大きさ
を測ればよい．
そこで，
µ
ˆij =
ni+ n+j
(i, j = 1, 2)
n++
(37)
とおき，(33),(34) に代入すればよい．ただし，独立性の検定においては，µ
ˆij の値は固定されたも
のではなく，標本の値に応じて変更されるため，χ2 分布の自由度は 1 である．
より一般的に，I × J の 2 重分類表の独立性の検定を行なうには，やはり，(37) によって µ
ˆij を
求める．このとき X 2 ,G2 は漸近的に自由度 (I − 1)(J − 1) の χ2 分布に従う．
X 2 と G2 を比較すると，X 2 の方が χ2 分布への収束が速い．N/IJ < 5 のような条件では，G2
を χ2 分布で近似することは難しい．
21
3.6
Fisher の正確検定 (Fisher’s Exact Test)
標本数 N が小さい場合には，X 2 や G2 の χ2 近似は難しい．この場合には，独立性の仮説のも
とでの正確な分布を用いて検定を行なうことができる．
独立性の仮説が成り立つとき，周辺度数が固定されているとの条件のもとで (つまり n1+ , n+1 の
両者が固定されている)，n11 は超幾何分布 (hyper geometric distribution) に従う．この分布は次
の式で表される．
(
p(n11 ) =
n1+
n11
)(
n2+
n21 = n+1 − n11
(
)
n
n+1
)
(38)
超幾何分布は，次のような確率を表すものと考えればよい．壷のなかに白い玉が n+1 個，赤い
玉が n+2 個含まれている．合わせて n 個である．ここで，壷のなかを見ずに，n1+ 個の玉を取り
出す．取り出したなかに含まれる白い玉の個数が n11 である．
超幾何分布の確率は上の式を用いて計算できるので，この分布と観測された n11 とを比較し，独
立性の仮説を棄却できるかどうか検討すればよい．
ここで，2 × 2 表の第 1 行の比率を p とおくと，n1+ = np となり，独立性の仮定のもとでの n11 の
平均と分散は超幾何分布の性質から次の式で与えられる．これらの式は後述する Cochran-Mantel-
Haenszel 検定で利用する．
n1+ n+1
n
n+1 n+2 n1+ n2+
Var(n11 ) =
(n − 1)n2
E(n11 ) = n+1 p =
3.7
(39)
(40)
シンプソンのパラドックスと層別の分析
ここで表 2 のデータをより詳細に学科別に検討してみる．表 6 は，学科別の内訳を示したもので
ある．各学科別に，オッズ比を求めると６学科のうち，女子の合格率が低いのは，ＣとＥの２つの
みであることがわかる．また，独立性のカイ２乗検定を行うと，顕著な関係がみられるのはＡ学科
のみであり，しかも傾向は女子の合格率の方が高い．これは全体を集計した傾向とは，逆向きのも
のである．
では，なぜこのような相反する傾向が生じたのだろうか．ここで，学科ごとに女子の志願者の比
率と男女合わせての合格率の散布図を図 8 に示す．これを検討すると，女子の志願者の多い学科
は，合格率の低いつまり入学の難しい学科であることがわかる．このため，個別の学科において
は，女子の合格率が低い傾向がなくとも，全体では女子の合格率が男子よりも低くなってしまう．
多重分割表の分析において，全体での傾向と層別での傾向とが相反するこのような逆説的な現象
は，シンプソンのパラドックス (Simpson’s paradox あるいは Yule-Simpson’s paradox) と呼ばれ
ている．
ここでデータの構造を直観的に理解するために，表 6 の内容をグラフに表示したものが，図 9 で
ある．この表示方法は Bertin(1981) によるもので，重み付きマトリックスと呼ばれている．この
場合は，各列の横方向の幅が各学科別・性別の志願者数に比例するようとってある．長方形の高さ
は，合格率と不合格率を表している．横方向の幅がそれぞれのケースの層別の人数に比例している
ので，長方形の面積が該当するセルの人数を表す．また，破線は全体での合格率と不合格率を表し
ており，これを超える部分は表示色を変えて強調してある．カイ 2 乗検定の結果からもわかるよう
22
表 6: 学科別・性別の合否の結果
学科
性別
合否
A
B
C
D
E
F
合計
男子
合格
512
313
353
207
120
205
138
279
53
138
22
351
1198
1493
62
89
63
17
37
202
33
131
28
94
6
24
45
557
19
82
8
68
391
34
244
35
299
24
317
7
1278
30
933
585
918
792
584
714
4526
オッズ比
0.35
0.80
1.13
0.92
1.22
0.83
1.84
φ
X2
-0.14
17.25
-0.02
0.25
0.03
0.75
-0.02
0.30
0.04
1.00
-0.02
0.38
0.61
92.20
(p 値)
0.00
0.61
0.39
0.59
0.32
0.54
0.00
不合格
合格率 (%)
女子
合格
不合格
合格率 (%)
合計
Source: Freedman et al. (1998) p.18, Agresti (2002) p.63). オッズ比が１より大きい（φ がゼロより大き
い）ことは，男子の合格率が女子より大きいことを示す．
に，Ａ学科を除いては男女で大きな合格率の違いはない．また，Ａ学科では女子の方が際立って合
格率の高いことがわかる．
調査データにおける因果推定の難しさ
4
多重分割表の分析を行って因果関係の推論を行う場合には注意すべき点がいくつかある．
4.1
多くの層に対応する推定法
ひとつは，推定上の（数値的な）性質に関わるものである．
ここでは説明しなかったが，多重分割表の構造を分析する手法として最も一般的なものはロジス
ティック回帰および対数線形モデルと呼ばれる手法である．これらのモデルの推定にあたっては，
通常，最尤法という手法を用いている．最尤法はデータの件数が十分あるなら，ある意味において
最善の推定方法であること，つまり最も推定精度の良いことが理論的に保証されている．しかしな
がら，各層に含まれるデータの件数が多くない場合（それぞれ数件などの場合）には，必ずしも適
切な推定値を与えるとは限らない．特に，ケースコントロールスタディとよばれる研究デザイン
（各層には実験群に属する被験者が１名，対照群に属する被験者が１名の計２名しか含まれない）
では，最尤法による推定では著しい偏りが生じる．このようなケースについては，条件付最尤法と
よばれる推定手法を用いるか，またはコクラン・マンテル・ヘンツェル (CMH) 検定などの手法を
用いる必要がある．
23
0.6
0.4
BA
C
E
D
0.2
男女合格率
0.8
1.0
UCBAdmissions
0.0
F
0.0
0.2
0.4
0.6
0.8
1.0
女子志願者率
図 8: 学科別の女子志願者率と合格率（男女込み）
学科別・性別の合格率・不合格率
合格率
図 9: 学科別性別合格率不合格率の重み付きマトリックス表示
24
Ｆ女子
Ｆ男子
Ｅ女子
Ｅ男子
Ｄ女子
Ｄ男子
Ｃ女子
Ｃ男子
Ｂ女子
Ｂ男子
Ａ女子
Ａ男子
不合格率
4.2
Cochran-Mantel-Haenszel 検定
そこで，層別のデータを対象として変数間の独立性の検定を行うための方法 (Cochran-Mantel-
Haenszel 検定，以下では CMH 検定と略す）が開発されている．この方法が検定する帰無仮説は，
「各層においてオッズ比がゼロ」であり，対立仮説は「各層がゼロでないある共通のオッズ比を持
つ」というものである．
CMH 検定の基本的なアイデアは，次のようなものである．まず，各層ごとに n11 の帰無仮説の
もとでの正確な分布（超幾何分布）を考え，その平均と分散を求めると次のようになる．
µ11 = E(n11 ) =
Var(n11 ) =
n1+ n+1
n++
n1+ n2+ n+1 n+2
n2++ (n++ − 1)
(41)
(42)
ついで，各層における帰無仮説のもとでの n11k の期待値 µ11k からのずれの和を求め，これらを
足し合わせる．この値は帰無仮説のもとで近似的に平均ゼロの正規分布に従うとみなせる．また分
散の大きさを求めることもできる．この和の平均と分散は各層のデータが独立に分布することか
ら，各々
E(n11+ ) =
∑
E(n11k ) , Var(n11+ ) =
∑
Var(n11k )
(43)
k
k
となる．さらに，連続修正を加えた次の値が, 帰無仮説（行と列が独立）のもとで近似的に自由度
1 の χ2 分布に従うことを用いて検定を行なう．
CMH =
(|n11+ − E(n11+ )| − 1/2)2
Var(n11+ )
(44)
この近似は，
1. 層の個数が固定されておりそれぞれの層におけるサンプル件数が増加する場合，
2. 個別の層内でのサンプル件数は限られているが，層の数が増加し全体では多くのサンプルが
存在する場合
のいずれにおいても成立する.
CMH 検定は，対立仮説として各層における変数の関係がほぼ同じ傾向を持つものであることを
仮定している．このため，層ごとに関係の方向が異なるデータにはうまく適用できない．
4.3
共通オッズ比の推定
また，層別された 2 × 2 表における共通オッズ比を推定するには，次の値 (Mantel-Haenszel 推
定量) が用いられる．
∑
∑
p11|k p22|k n++k
(n11k n22k /n++k )
= ∑k
θˆMH = ∑k
k (n12k n21k /n++k )
k p12|k p21|k n++k
(45)
ここで，pij|k は k 層全体における ij セルの比率を表す．また，Rk = n11k n22k /n++k ，Sk =
∑
∑
n12k n21k /n++k とすると，θˆMH = R/S = (
Rk )/(
Sk ) となる．
k
k
25
この推定量の対数の分散は次の式で近似される (Agresti, 2002, 6.3.5)．
σ
ˆ 2 (log θˆMH )
=
1 ∑ −1
n++k (n11k + n22k )Rk
2R2
k
1 ∑ −1
+ 2
n++k (n12k + n21k )Sk
2S
k
1 ∑ −1
+
n++k [(n11k + n22k )Sk + (n12k + n21k )Rk ]
2RS
k
4.4
一般化 Cochran-Mantel-Haenszel 検定
Cochran-Mantel-Haenszel 検定を拡張し，I × J の表が K 層あるデータにおける行と列の独立
性の検定を行うことができる. この方法は，一般化 Chchran-Mantel-Haenszel 検定と呼ばれる．２
×２の表において周辺度数が固定された場合の n11 は独立性の仮定 (H0 ) のもとで，超幾何分布に
従う．一方，周辺度数が固定された I × J 表では，(I − 1) × (J − 1) の自由度が存在し，ひとつの
層 (ここでは k で指定する) においては，パラメータ
nk = (n11k , ..., n1,(J−1),k , ..., n(I−1),1,k , ..., n(I−1),(J−1),k )
(46)
によって表現される．
超幾何分布を拡張することにより，行と列が独立である場合の，k 層におけるこれらの値 (46) の
平均 µk と分散共分散行列 V k を数学的に求めることができる．k を省略して表記すると，ベクト
ル要素の平均は µij = ni+ n+j /n++ であり，また共分散は次の式である（ここで δij は i と j が等
しいときは 1 であり，異なるときには 0 となる記号=クロネッカーデルタ）．
Cov(nij , ni j ) =
ni+ (δii n++ − ni + )n+j (δjj n++ − n+j )
n2++ (n++ − 1)
(47)
ここで, ベクトルの和について，次のように表す．
n=
∑
k
nk , µ =
∑
k
µk , V =
∑
Vk
(48)
k
２×２表の場合と類似の形式を持つ統計量を行列とベクトルの積
一般化 CMH = (n − µ)T V −1 (n − µ)
(49)
によって定義する（上の式は有限補正なしの場合）．ベクトル n が近似的に多変量正規分布に従う
ことを用いると, 独立性の仮定のもとで (49) は自由度 (I − 1)(J − 1) の χ2 分布で近似される. こ
れを用いることにより,CMH 統計量が大きければ行と列が独立であるとの仮説を棄却する.
また, 行と列に順序構造がある場合には，行および列のカテゴリーにスコアを与えることにより，
関係性をより敏感に検出するための方法も提案されている (Agresti,2002, 7.5)．
4.5
影響関係の方向
もうひとつの問題は，より概念的な枠組みに関わるものである．上に述べた例では，性別と合格
率の双方に学科が関係しているため，これについて層別することにより，妥当な検討を行うことが
26
できた．では，一般的に，変数 X と変数 Y の関係を調べる場合，それら両者に関係している全て
の変数 Z1 ,...,Zp を列挙し，これらについて同時に層別を行えばよいだろうか．
これについては，影響関係の方向までも含めてより慎重に検討する必要がある．図 10 は，影響
関係についてのいくつかのケースを示したものである．ここで取り上げた例は (i) に示したもので
ある．多重分割表による分析は影響関係の強さを記述することを目的としており，影響関係の向き
についての情報を用いてはいない．しかし，現実の解釈においては，影響関係の向きが重要にな
る．上に示した合格率の例 (i) では，
「学科」について層別することにより，
「性別」(X) の「合否」
(Y ) に及ぼす直接の影響を検討することが可能になる．また (ii) の例のように X と Y とに Z が共
通に影響を与えている場合も，Z について層別を行うことにより，X から Y への直接の影響を推
定することができる．
しかし，常に層別を行うことが妥当であるわけではない．(iii) の例のように影響関係の下流にあ
る変数について層別を行うと，知ろうとしている関係とは異なる見かけ上の関係が X と Y に生じ
るため，分析方法として適切ではない．
ቇ⑼
ᕈ೎
䋿
䋿
วุ
ᕈ೎
(i) 学科による層別（適切）
ᕈ೎
วุ
ቇ⑼
(ii) 性別による層別（適切）
䋿
ቇ⑼
วุ
(iii) 合否による層別（不適切）
図 10: 影響構造と層別の適切さ
4.6
未知の背景要因の影響
また外的な介入の有無は，影響関係の推論の確さと深い関係を持つ. 一般的には，外的介入なし
での推論から決定的な結論を得るのは難しい．以下は，介入を伴わない（大量の）調査データから
27
得られた結論が，他の証拠なしには，明快な解釈が困難であることを示している．
Smith et al. (1992) は，調査データの結果の解釈に，慎重な注意が必要であることを指摘し
た論文である．ここで著者らは，大規模な危険要因についての調査 (MRFIT) を分析することに
よって，まず喫煙と自殺の間に関係が見られることを指摘している．MRFIT(Multiple Risk Factor
Intervention Trial) は米国の 35 歳から 57 歳までの 361,662 人の男性を 12 年間にわたって追跡調
査した研究である．年齢，人種，1 日の喫煙本数，糖尿病治療中であるか，心筋梗塞の有無，社会
経済地位（居住地域から推定）がベースラインデータとしてとられており，これらの効果を考慮し
てもなお喫煙と自殺との間に関係は見られる．また，MRFIT 以外の他の大規模な複数の追跡調査
でも同様の結果が得られている．
しかし，Smith et al.(1992) は，MRFIT について喫煙と他殺による死の間にも同様の関係が見
られることを報告している．喫煙がその薬理作用によって自殺を増加させるという説明は，一見
もっともらしく思えるが，他殺を増加させる原因となるとは考えづらい．喫煙は本当に自殺に影響
を及ぼしているのだろうか？
これについての決定的な答えは現状では分かっていない．しかし，タバコが直接自殺の危険を
高めるとこれらの調査結果から解釈するのは，困難であるように思える．研究者によっては捉えら
れていない何らかの心理的または社会的要因があり，これが喫煙の習慣と自殺傾向の両者に影響を
及ぼしているとみなすことも可能と思われる．
ᐕ㦂
෼౉
䉺䊋䉮
ੱ⒳
䌘
∛᳇
䋿
⥄Ვ
図 11: 喫煙の影響構造
表 7: 喫煙行動と自殺率 (Smith et al.,1992)
人数
自殺件数
自殺件数
1 万人×年
0
228,545
291
1.09
1-19
20-39
40-59
29,333
72,200
27,844
50
166
78
1.47
2.00
2.46
60+
3,740
16
3.78
タバコ/日
28
表 8: リスク要因と自殺の相対リスク (95% 信頼区間) (Smith et al. 1992)
リスク要因
年齢補正済み
完全補正済み相対リスク相対リスク
0
1-19
20-39
1.00
1.36(1.00-1.83)
1.88 (1.55-2.27)
1.00
1.36(1.00-1.84)
1.86(1.54-2.26)
40-59
60+
2.31 (1.79-2.96)
3.44 (2.08-5.69)
2.27(1.77-2.92)
3.33(2.01-5.52)
(78.66,p < 0.0001)
(75.98,p < 0.0001)
1.35(1.13-1.60)
1.33(1.13-1.59)
1.00
1.00
0.85(0.60-1.21)
1.00
0.76(0.53-1.08)
1.00
1.79(1.07-2.99)
1.00
1.73(1.03-2.90)
1.00
1.66(0.99-2.77)
1.00
1.61(0.96-2.70)
1.00
喫煙
(χ2 for trend)
収入
低
他
人種
黒人
他
心筋梗塞
病歴あり
なし
糖尿病
有病
なし
表 9: タバコと他殺による死 (Smith et al. 1992 に基づく)
総数 222, 人種と収入 (居住地域から推定) について補正された
10 万人あたり死亡率の比
リスク要因
他殺による死亡の相対リスク
(一日あたり喫煙本数)
(95% 信頼区間)
0
1.00
1-39
40-
1.71 (1.29-2.28)
2.04 (1.32-3.15)
29
ロジスティック回帰 (Logistic Regression)
5
回帰分析 (分散分析・共分散分析を含む) は，被説明変数が連続な値を持つ場合の方法であるが，
時には 0 − 1 反応または比率を被説明変数とする必要が生じる場合がある．ロジスティック回帰は，
このような分析のための方法である．
線形の回帰分析は，説明変数 x = (x1 , ..., xp ) を固定した場合の Y の平均を次のような式で表す．
E(Y |x) = µx = β0 + β1 x1 + · · · + βp xp
(50)
さらに，Y は µx を平均とする正規分布 N (µx , σ 2 ) に独立に従うと仮定している．
ロジスティック回帰は，上の 2 点を次のように変更したものである．
1) まず，説明変数 x を固定した際，被説明変数 Y は 2 項分布 Bin(mx , πx ) に独立に従うとする．
2) さらに πx = E(Y |x) が，次のような式で表されると仮定する．
logit(πx ) = log
πx
= ηx = β0 + β1 x1 + · · · + βp xp
1 − πx
上の式は，ロジット (logit) 変換と呼ばれる．一方，log
πx =
πx
1−πx
(51)
= ηx の逆関数は，
exp(ηx )
1 + exp(ηx )
(52)
となる．これはロジスティック (logistic) 関数と呼ばれる．ηx の値が −∞ から +∞ まで変化する
とき，πX は 0 から 1 まで変化する．
ˆ を推定すべきかである．線形の回帰分析の場合
次に考えるべきは，標本からどのようにして β
には，最小 2 乗基準により次の式
N
∑
(yi − yˆi )2
(53)
i=1
ˆ を求めた．単純に考えるなら，i 番目の観測値に対応する πi と yi /mi
を最小化することによって，β
を考え，πi と yi /mi の残差 2 乗和を最小化する方法が考えられるが，通常はこのような推定は行
わない．
標準的には最尤法 (maximum likelihood method) と呼ばれる推定方法が用いられる．モデルが
正しくしかもデータが十分にあれば，最尤法は優れた推定方法であることが，理論的に知られて
いる。
ロジスティック回帰の対数尤度（モデルの当てはまりの良さの指標）は，
l(β|y1 , ..., yN ) =
定数 +
N
∑
{yi log πi + (mi − yi ) log (1 − πi )}
(54)
i=1
N
∑
= 定数 +
{yi ηi − mi log (1 + exp(ηi ))}
(55)
i=1
で与えられるので (πi と ηi は β と xi によって値が決まる)，これを最大化する係数 β を数値的に
求め，推定値とする．
5.1
ロジスティック回帰の例
表 10 は，Dobson の教科書の表 8.2 である．このデータは薬品 (二硫化炭素ガス) を与えたあと，
５時間後の昆虫の死亡数を表している．この例では，説明変数は定数項の他には，薬品の投与量
30
(X1 ) のみである．個体数が，2 項分布 Bin(m, π) の試行数 m であり，昆虫の死亡数が被説明変数
Y である．
モデルは次の 2 つの性質を持つ．
1. Yi ∼ Bin(mi , πi ) , (i = 1, ..., N )
2. logit(πi ) = β0 + β1 xi , (i = 1, ..., N )
係数の推定値を表 11 に示す．
表 10: 甲虫の死亡率
投与量
個体数
死亡数
1.6907
59
6
1.7242
1.7552
1.7842
60
62
56
13
18
28
1.8113
1.8369
63
59
52
53
1.8610
1.8839
62
60
61
60
(log10 CS2 mg/l )
表 11: 推定値
最尤推定値
標準誤差
β0
-60.7175
5.1807
β1
34.2703
2.9120
分析すべきデータを yi , mi , xi (i = 1, ..., N ) とする．ここで，yi は i 番目のデータにおける反応
件数 (成功など) であり，mi は i 番目のデータの総数である．また，xi は，説明変数の値を要素と
する p + 1 次元のベクトル (1, xi1 , ..., xip ) とする．
5.2
仮説検定
ある回帰係数がゼロであるとの仮説 H0 : βi1 = · · · = βip−q = 0 を検定するには，尤度比検定を
ˆ とし，制約されたモデルにおける
用いる．つまり，制約されないモデルにおける最尤推定値を β
˜ とすると，帰無仮説 H0 が真なら
最尤推定値を β
ˆ
˜
2l(β|y)
− 2l(β|y)
∼ χ2p−q + O(1/n)
(56)
である．ここで，p は制約されないモデルにおける推定パラメータの数であり，q は制約されたモ
デルにおける推定パラメータの数である．この値を χ2 分布と比較し，通常は上側 5%に入れば帰
無仮説が棄却されるものとする．つまり，パラメータを制限することによる当てはまりの低下程度
が大きければ，これらがゼロではないであろうと結論する．
31
0.6
0.4
0.0
0.2
dob82$Y/dob82$m
0.8
1.0
Logistic Regression
1.70
1.75
1.80
1.85
dob82$x
図 12: ロジスティック回帰分析 (円:観測値, 破線:推定値)
6
項目反応理論 (IRT)
医学研究ではあまり一般的ではないが，心理テストや各種の判定テストなどの分析に有効である
手法に，
「項目反応（応答）理論」(item response theory, IRT) と呼ばれる統計分析法がある．この
手法は，学力試験の品質（難度，学力の識別力）等を一定に保つために心理統計の分野で開発され
てきた．IRT は 1 次元の因子分析 (factor analysis) を，観測データが 2 値である場合に対応させた
ものとみなせる．
IRT は直接には観測されない潜在変数 θ（通常は１次元）の存在を仮定し，観測されている２値変数
Yj , (j = 1, ..., p) （学力試験の場合には正誤反応になる）の期待値（条件付正答率）πj , (j = 1, ..., p)
が，潜在変数のロジスティック関数によって表されると仮定する．
潜在変数 θ は直接は観測されないので，観測されている変数 Yj , (j = 1, ..., p) から間接的に推定
される．通常 θ は，平均０，分散 1 の正規分布に従うことを仮定し，各 Yj の平均（正答率）は，θ
によって定まり，さらに θ の値を固定したもとでは，各 Yj は統計的に独立であることを仮定する．
この最後の性質は「局所独立性」と呼ばれる．
広く用いられているモデル（3 パラメータモデル) では，潜在変数 θ を仮定したもとでの Yj の平
均 πj (θ) を，次のようなロジスティック関数を若干変形した式で表す．
πj (θ)
=
cj +
1 − cj
1 + exp{−aj (θ − bj )}
(57)
多肢選択式の試験問題への解答データでは，cj は選択肢数の逆数（に近い数）を想定している．潜
在変数と項目への応答 Yj との関係を記述する aj , bj , cj , (j = 1, ..., p) は最尤法などを用いて推定さ
れる．ただし，３つのパラメータをすべてデータから推定することは実際にはかなり難しい場合が
多く，cj をゼロと仮定するモデル（2 パラメータモデル) を利用したり，あるいは cj の値をあらか
32
じめ設定しておく場合もある．Rasch モデルは上の式において cj = 0 とし，さらに aj をすべての
j について共通であると仮定したものである．
参考文献
[1] Agresti,A. (2003) カテゴリカルデータ解析入門, サイエンティスト社 (原著 1996, An introduction to
categorical data analysis, Wiley). （原著第２版が 2007 年に出版された. 訳本は初版）
[2] Agresti,A. (2002). Categorical data analysis, Wiley. 前項よりかなり詳しい教科書．710 ページ.
[3] Bertin,J. (1981). Graphics and graphic information-processing, (translated by W.J.Berg and P.Scott)
New York: Walter de Gruyter.
[4] Bickel,P.J., Hammel,E.A., & O’Connell,J.W. (1975). Sex bias in graduate admissions: Data from
Berkeley, Science, 187, 398-404.
[5] Freedman,D., Pisani,R. & Purves,R. (1998). Statistics, 3rd ed., New York: W.W. Norton. (やさし
い初等統計の教科書。専門家の評価が高い。）
[6] Ihaka,R. & Gentleman,R. (1996). R: A language for data analysis and graphics, Journal of Computational and Graphical Statistics, 5, 299–314. ( http://www.r-project.org/) 日本語の解説本が
数種販売されるようになった. SAS のような解説マニュアルが完備しているわけではないが，オンライ
ンマニュアル (英語) が付いている．適切に設定すれば日本語も利用できる．
[7] Pearl,J. (2000). Causality: models, reasoning, and inference. Cambridge University Press. （邦訳
あり）因果関係の確率論的モデルについての独創的な論考．かなり難解．
[8] Smith,G.D., Phillips,A.N. & Neaton,J.D. (1992). Smoking as “independent” risk factor for suicide:
illustration of an artifact from observational epidemiology? Lancet, 340, 709-12.
[9] Dobson,A.J. 田中豊他訳 (1993) 統計モデル入門 – 回帰モデルから一般化線形モデルまでー . 共立出
版. (原著 An Introduction to Generalized Linear Models. Chapman and Hall (1990) )
[10] Draper,D. et al. (1992). Contemporary Statistics,1; Combining Information, American Statistical
Association.
[11] Draper,N.R. & Smith,H. (1998). Applied Regression Analysis, 3rd ed., Wiley. (大部であるが回帰分
析全般についてのわかりやすい解説．森北出版の中村慶一訳は第 1 版．)
[12] Fitzmaurice,G.M., Laird,N.M., & Ware,J.H. (2004) Applied Longitudinal Analysis, Wiley. (継時的
な繰り返し測定や，クラスター化された標本の分析に用いられる手法である混合モデルについての平易
な解説．SAS の proc MIXED や GENMOD の利用法についての説明がある．ハーバード公衆衛生学部
の生物統計の研究者による教科書.)
[13] Gelman,A., & Hill,J. (2006) Data Analysis Using Regression and Multilevel/Hierachical Models,
Cambridge University Press.
Gelman は D.Rubin 門下のベイズ統計の著名な研究者．回帰（ロジステッィク回帰およびポアソン回帰
を含む）と階層モデルについての解説. ソフトウェアとしては BUGS と R の利用法の解説がある. BUGS
は英国のメディカルリサーチカウンシルの生物統計研究グループが開発したベイズ推定のソフトウェア．
[14] Kleinbaum,D.G. & Klein,M. (2002). Logistic regression: A self learning text, 2nd ed. Springer. （わ
かりやすい教科書．ワークブックとして使える．多値反応や相関のある反応の推定方法も扱っている．
）
[15] Nolan,D. & Speed,D. (2000). Stat Labs: Mathematical Statistics Through Applications. Springer.
(統計入門コースのための実データに基づくワークブック)．
[16] 丹後俊郎・高木晴良・山岡和枝 (1996). ロジスティック回帰分析 - SAS を利用した統計解析の実際. 朝
倉書店. （日本語の教科書としては，最も詳しいかも. 分析例はいずれも医学系の題材）.
[17] Yerushalmy,J. (1971). The relationship of parents’ cigarette smoking to outcome of pregnancy –
implications as to the problem of inferring causation from observed associations. American Journal
of Epidemiology, 93, 443–456.
33

Download Report