確率・統計 (21) ポアソン回帰 (Poisson Regression)

ある事象が発生した回数を記録したデータは非常に多く見かけられ、活用されています。例えば、台風や大雨、地震といった災害の年間発生回数や、店舗にてある商品が一日に売れる個数、ガンなどの病気による年間死亡者数などがそれに該当します。このようなデータは「計数データ (Count Data)」または「度数データ (Frequency Data)」と呼ばれ、季節や月単位といった期間や、商品別・性別・喫煙の有無等の事象によってヒストグラムなどのグラフを使って傾向を観察することはよく行われる推定方法です。一般化線形モデルで計数データを扱う場合は、「ポアソン回帰 (Poisson Regression)」や「対数線形モデル (Log-linear Model)」がよく利用されます。この章では、これらのモデルを「分散分析」と比較しながら紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 分散分析と一般化線形モデル

以前は、一元配置分散分析法を表す線形重回帰式として以下のものを示しました。

y_i,j = α₀ + α₁δ₁ + ... + α_jδ_j + ... + α_p-1δ_p-1 + ε_i,j

ここで、j は事象に対する番号で、その総数は p になります。また、i は各事象が持つ従属変数に対する番号であり、ここではその数を n_j で表します。α_j は回帰係数、δ_j は各事象に対応したダミー変数で、j 番目の事象に属する場合だけ 1 で、その他は 0 になります。よって y_i,J に対応する ( j = J のときの ) 式は

y_i,J = α₀ + α_J + ε_i,J

となります。但し、j = p の場合は該当するダミー変数がないことから、回帰式は

y_i,p = α₀ + ε_i,p

になります。これは、α₀ が基準となって、1 から p - 1 までの事象との差異が α_j で表されると考えると理解しやすいと思います。

最後に、y_i,j が従属変数、ε_i,j が誤差成分で、ε_i,j は正規分布 N( 0, σ² ) に従うと仮定します。この誤差成分は、事象によらずバラツキが一定であると考えるわけです。

この式は、各事象に対して n_j 個存在するので、全事象での式の総数を N とすると、この値は

N = Σ_j{1→p}( n_j )

で求められます。この式を、事象ごとの式がまとめるように順に並べてみると

y_1,1

α₀

α₁

ε_1,1

y_2,1

α₀

α₁

ε_2,1

y_n₁,1

α₀

α₁

ε_n₁,1

y_1,2

α₀

α₂

ε_1,2

y_i,j

α₀

α_j

ε_i,j

y_{n_p-1,p-1}

α₀

α_p-1

ε_{n_p-1,p-1}

y_{n_p,p}

α₀

ε_{n_p,p}

となります。ここで、全成分が 1 からなる n 次元ベクトルを 1_n、全成分が 0 からなる n 次元ベクトルを 0_n で表して

X	=	\|	1_n₁,	1_n₁,	0_n₁,	...	0_n₁	\|
		\|	1_n₂,	0_n₂,	1_n₂,	...	0_n₂	\|
		\|	:	:	:	...	:	\|
		\|	1_{n_p-1},	0_{n_p-1},	0_{n_p-1},	...	1_{n_p-1}	\|
		\|	1_{n_p},	0_{n_p},	0_{n_p},	...	0_{n_p}	\|

α	=	\|	α₀	\|
		\|	α₁	\|
		\|	:	\|
		\|	α_p-1	\|

とし、y を、左辺を並べた N 次元ベクトル、ε を、誤差成分を並べた N 次元ベクトルとすれば、先ほど示した式は

となります。これが一元配置分散分析法における線形重回帰モデル式で、X は「デザイン行列 (Design Matrix)」を意味します。a を α の最尤推定量としたとき、線形重回帰モデル式の正規方程式は

X^TXa = X^Ty

X^TX

1_n₁^T,

1_n₂^T,

...

1_{n_p-1}^T,

1_{n_p}^T

1_n₁,

0_n₁,

...

0_n₁

1_n₁^T,

0_n₂^T,

...

0_{n_p-1}^T,

0_{n_p}^T

1_n₂,

0_n₂,

1_n₂,

...

0_n₂

0_n₁^T,

1_n₂^T,

...

0_{n_p-1}^T,

0_{n_p}^T

...

1_{n_p-1},

0_{n_p-1},

...

1_{n_p-1}

0_n₁^T,

0_n₂^T,

...

1_{n_p-1}^T,

0_{n_p}^T

1_{n_p},

0_{n_p},

...

0_{n_p}

n₁,

n₂,

...

n_p-1

n₁,

...

n₂,

...

n_p-2,

...

n_p-1,

...

n_p-1

X^Ty	=	\|	1_n₁^T,	1_n₂^T,	...	1_{n_p-1}^T,	1_{n_p}^T	\|\|	y₁	\|
		\|	1_n₁^T,	0_n₂^T,	...	0_{n_p-1}^T,	0_{n_p}^T	\|\|	y₂	\|
		\|	0_n₁^T,	1_n₂^T,	...	0_{n_p-1}^T,	0_{n_p}^T	\|\|	:	\|
		\|	:	:	...	:	:	\|\|	y_p-1	\|
		\|	0_n₁^T,	0_n₂^T,	...	1_{n_p-1}^T,	0_{n_p}^T	\|\|	y_p	\|
	=	\|	S_y	\|
		\|	S₁	\|
		\|	S₂	\|
		\|	:	\|
		\|	S_p-1	\|

となります。但し S_j は事象 j に対する従属変数の和 Σ_i{1→n_j}( y_i,j )、S_y は全従属変数の和 Σ_j{1→p}( S_j ) をそれぞれ表します。簡単な考察から、X^TX の逆行列 ( X^TX )^-1 は

( X^TX )^-1	=	\|	1 / n_p,	-1 / n_p,	-1 / n_p,	...	-1 / n_p	\|
		\|	-1 / n_p,	( n₁ + n_p ) / n₁n_p,	1 / n_p,	...	1 / n_p	\|
		\|	:	:	:	...	:	\|
		\|	-1 / n_p,	1 / n_p,	1 / n_p,	...	( n_p-1 + n_p ) / n_p-1n_p	\|

a =

( X^TX )^-1X^Ty

\|	1 / n_p,	-1 / n_p,	-1 / n_p,	...	-1 / n_p	\|\|	S_y	\|
\|	-1 / n_p,	( n₁ + n_p ) / n₁n_p,	1 / n_p,	...	1 / n_p	\|\|	S₁	\|
\|	:	:	:	...	:	\|\|	:	\|
\|	-1 / n_p,	1 / n_p,	1 / n_p,	...	( n_p-1 + n_p ) / n_p-1n_p	\|\|	S_p-1	\|

\|	[ S_y - Σ_j{1→p-1}( S_j ) ] / n_p	\|
\|	[ -S_y + S₁( n₁ + n_p ) / n₁ + Σ_j{2→p-1}( S_j ) ] / n_p	\|
\|	:	\|
\|	[ -S_y + Σ_j{1→p-2}( S_j ) + S_p-1( n_p-1 + n_p ) / n_p-1 ] / n_p	\|

\|	S_p / n_p	\|	≡	\|	m_p	\|
\|	S₁ / n₁ - S_p / n_p	\|		\|	m₁ - m_p	\|
\|	:	\|		\|	:	\|
\|	S_p-1 / n_p-1 - S_p / n_p	\|		\|	m_p-1 - m_p	\|

となります。但し、m_j は事象 j に対する従属変数の平均を表します。

線形重回帰モデルにおける飽和モデルとの対数尤度統計量 D は

D = Σ_i{1→N}( ( y_i - x_i^Ta )² ) / σ²

D	=	{ Σ_j{1→p-1}( Σ_i{1→n_j}( [ y_i,j - ( a₀ + a_j ) ]² ) ) + Σ_i{1→n_p}( ( y_i,p - a₀ )² ) } / σ²
	=	{ Σ_j{1→p-1}( Σ_i{1→n_j}( [ y_i,j - ( m_p + m_j - m_p ) ]² ) ) + Σ_i{1→n_p}( ( y_i,p - m_p )² ) } / σ²
	=	Σ_j{1→p}( Σ_i{1→n_j}( ( y_i,j - m_j )² ) ) / σ²

という結果が得られます。これは、一元配置分散分析において集団内の平方和 S_E を求める式に相当し、自由度 N - p の χ²-分布に「正確に」従います。

回帰係数を一つとし、y_i,j = α₀ + ε_i,j である場合を考えると、線形重回帰モデルの場合における飽和モデルとの対数尤度統計量 D₀ は

D₀ = Σ_i{1→N}( ( y_i - m_y )² ) / σ²

となるのでした (*1-3)。但し、m_y は全従属変数の平均 S_y / N を表します。一元配置分散分析の場合、これが全体の平方和 S_T を意味し、D₀ - D が集団間の平方和 S_C になります。D₀ は自由度 N - 1 の χ²-分布に従うので、D₀ - D は自由度 p - 1 の χ²-分布に従うことになり、

F = [ ( D₀ - D ) / ( p - 1 ) ] / [ D₀ / ( N - p ) ]

は自由度 ( p - 1, N - p ) の F-分布に従うことになります。この値 F を使って検定を行うのが一元配置分散分析でした (*1-4)。

一元配置分散分析における重回帰線形モデル式は、以下のモデル式

y_i,j = α₀ + α₁δ₁ + ... + α_jδ_j + ... + α_p-1δ_p-1 + α_pδ_p + ε_i,j

において α_p = 0 とした場合に相当します。上記の式をそのまま y = Xα + ε の形にすると、

X	=	\|	1_n₁,	1_n₁,	0_n₁,	...	0_n₁,	0_n₁	\|
		\|	1_n₂,	0_n₂,	1_n₂,	...	0_n₂,	0_n₂	\|
		\|	:	:	:	...	:	:	\|
		\|	1_{n_p-1},	0_{n_p-1},	0_{n_p-1},	...	1_{n_p-1},	0_{n_p-1}	\|
		\|	1_{n_p},	0_{n_p},	0_{n_p},	...	0_{n_p}	1_{n_p}	\|

α	=	\|	α₀	\|
		\|	α₁	\|
		\|	:	\|
		\|	α_p-1	\|
		\|	α_p	\|

となって、X の列数と α の次元数は p + 1 となります。X^TX は

X^TX	=	\|	N,	n₁,	n₂,	...	n_p	\|
		\|	n₁,	n₁,	0,	...	0	\|
		\|	n₂,	0,	n₂,	...	0	\|
		\|	:	:	:	...	:	\|
		\|	n_p-1,	0,	0,	...	0	\|
		\|	n_p,	0,	0,	...	n_p	\|

と計算できますが、第一列 ( または第一行 ) は他の列(行)の和に等しいことから線形従属な列(行)ベクトルが存在することになって、X^TX は特異行列であることになります。このままでは正規方程式 X^TXa = X^Ty の解が一意に決まりません。そのため、α_p = 0 として ( α₀ の中に含めることで ) 非特異な行列になるようにしていたわけです。これは「端点制約 (Corner-point Constraint)」と呼ばれます。

X^TXa = X^Tyを連立方程式の形で書き直すと以下のようになります。

Na₀ + n₁a₁ + n₂a₂ + ... + n_pa_p	=	S_y
n₁a₀ + n₁a₁	=	S₁
n₂a₀ + n₂a₂	=	S₂
:		:
n_pa₀ + n_pa_p	=	S_p

このままでは解は一意に決まらないので、n₁a₁ + n₂a₂ + ... + n_pa_p = Σ_j{1→p}( n_ja_j ) = S_y - Nt とすれば、一番目の式から

Na₀ + S_y - Nt = S_y より

a₀ = t

a_j = S_j / n_j - t = m_j - t

となり、a₀ = t = S_y / N = m_y のとき、最初に仮定した式から Σ_j{1→p}( n_ja_j ) = S_y - Nt = 0 で

という結果が得られます。この結果を使って飽和モデルとの対数尤度統計量 D を求めると

D	=	Σ_i{1→N}( ( y_i - x_i^Ta )² ) / σ²
	=	Σ_j{1→p}( Σ_i{1→n_j}( [ y_i,j - ( a₀ + a_j ) ]² ) ) / σ²
	=	Σ_j{1→p}( Σ_i{1→n_j}( ( y_i,j - m_j )² ) ) / σ²

となって、端点制約を用いた結果と等しくなります。これは「零和制約 (Sum-to-zero Constraint)」と呼ばれる手法になります。

次に、以下のような線形重回帰モデル式について検討してみます。

y_i,j,k = μ + α_j + β_k + γ_j,k + ε_i,j,k

μ は従属変数全体の基準値で、α_j はある事象 (これを A とします) ごとの基準との差異を、また β_k は他の事象 (これを B とします) ごとの基準との差異をそれぞれ表し、最後の γ_j,k は事象 A, B の組み合わせによって生じる差異を表します。α_j, β_k は事象 A, B それぞれの「主効果 (Main Effect)」、γ_j,k は事象 A, B の「交互作用効果 (Interaction Effect)」と呼ばれます。添字の i は、事象 ( A, B ) の組み合わせごとの変数に付けられる連番であり、各変数の個数は全て等しく n であるとします。事象 A, B の数をそれぞれ p, q としたとき、連立方程式は

y_1,1,1	=	μ	+	α₁	+	β₁	+	γ_1,1
y_2,1,1	=	μ	+	α₁	+	β₁	+	γ_1,1
:
y_n,1,1	=	μ	+	α₁	+	β₁	+	γ_1,1
y_1,2,1	=	μ	+	α₂	+	β₁	+	γ_2,1
:
y_n,p,1	=	μ	+	α_p	+	β₁	+	γ_p,1
y_1,1,2	=	μ	+	α₁	+	β₂	+	γ_1,2
:
y_i,j,k	=	μ	+	α_j	+	β_k	+	γ_j,k
:
y_n,p,q	=	μ	+	α_p	+	β_q	+	γ_p,q

となって、全部で npq 個の式が得られます。しかし、事象 A, B の組み合わせが等しい従属変数に対する式は等しく、独立した式は高々 pq 個しかありません。それに対して求めるべき未知数は、μ と、α_j が p 個、β_k が q 個、γ_j,k が pq 個で合わせて pq + p + q + 1 = ( p + 1 )( q + 1 ) 個あるので、ここでも零和制約や端点制約を使って一意の解にする必要があります。

<μ>

α_j (p列)

β_k (q列)

γ_j,1 (p列)

γ_j,2 (p列)

...

γ_j,q (p列)

1_n,

0_n,

...

0_n,

1_n,

0_n,

...

0_n,

1_n,

0_n,

...

0_n,

...

0_n,

...

0_n,

...

0_n

1_n,

0_n,

1_n,

...

0_n,

1_n,

0_n,

...

0_n,

1_n,

...

0_n,

...

0_n,

...

0_n,

...

0_n

( np 行 )

...

1_n,

0_n,

...

1_n,

0_n,

...

0_n,

...

1_n,

0_n,

...

0_n,

...

0_n,

...

0_n

1_n,

0_n,

...

0_n,

1_n,

...

0_n,

...

0_n,

1_n,

0_n,

...

0_n,

...

0_n,

...

0_n

...

1_n,

0_n,

...

1_n,

0_n,

...

1_n,

0_n,

...

0_n,

...

0_n,

...

0_n,

...

1_n

α	=	\|	μ	\|
		\|	α₁	\|
		\|	:	\|
		\|	α_p	\|
		\|	β₁	\|
		\|	:	\|
		\|	β_q	\|
		\|	γ_1,1	\|
		\|	:	\|
		\|	γ_p,q	\|

となって、X は npq 行 ( p + 1 )( q + 1 ) 列の行列、α は ( p + 1 )( q + 1 ) 次元のベクトルになります。X の内容がわかりづらいので、次の部分行列(ブロック)

D	=	\|	1_n,	0_n,	...	0_n	\|
		\|	0_n,	1_n,	...	0_n	\|
		\|	:	:	...	:	\|
		\|	0_n,	0_n,	...	1_n	\|

					k 列目
F_k	=	\|	0_n,	...	1_n,	...	0_n	\|
		\|	0_n,	...	1_n,	...	0_n	\|
		\|	:	...	:	...	:	\|
		\|	0_n,	...	1_n,	...	0_n	\|

を使って X を表します。ここで、D は np 行 p 列の行列、F_k は np 行 q 列の行列とします。また、0_np,p を np 行 p 列の零行列とすれば、X は

X	=	\|	1_np,	D,	F₁,	D,	0_np,p,	...	0_np,p	\|
		\|	1_np,	D,	F₂,	0_np,p,	D,	...	0_np,p	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_np,	D,	F_q,	0_np,p,	0_np,p,	...	D	\|

となります。しかし、D を列とする部分行列の各列 ( α_j に対する独立変数成分 ) の和は最左端にある 1_np が並ぶ列に等しくなり、F_k を列とする部分行列の各列 ( β_k に対する独立変数成分 ) の和も同じく 1_np が並ぶ列に等しくなります。そこで、それぞれの最後の列を除き、D を np 行 p - 1 列の、F_k を np 行 q - 1 列の行列とします。このとき、F_q は全ての要素がゼロになります。これは、未知数 α_p, β_q を除外する操作と同じ意味を持ちます。また、この操作によって、D と F_k の全ての列も互いに独立になります。
次に、D が対角上に並んだ後半部の行列 ( γ_j,k に対する独立変数成分 ) の各列を任意の組み合わせで加算することで、その左側に並ぶ 1_np, D, F_k が並ぶ列と等しくすることができるので、D が並ぶ列と独立になるように対角上の最後の D を除外し、F_k が並ぶ列と独立になるように各対角上の D から最後の p 列目を除きます。これらの操作は、γ_j,q ( 1 ≤ j ≤ p ) と γ_p,k ( 1 ≤ k ≤ q - 1 ) を取り除いたことに等しくなります。これでようやく互いの列が独立になります。

除外した列の数は、α_p, β_q の 2 個と γ_j,q の p 個、それに γ_p,k が q - 1 個なので、全部で p + q + 1 個になります。従って未知数の数は pq 個で式の数と等しくなり、解が一意に決まります。具体的に n = 1, p = 3, q = 3 の時の X を例に今までの操作を行うと次のようになります。

X	=	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
		\|	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

2, 3, 4 列目の和が 1 列目になるので 4 列目を除去 ( 1, 2, 3 列目は独立になる )。

→	\|	1,	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

4, 5, 6 列目の和が 1 列目になるので 6 列目を除去 ( 1 から 5 列目までは独立になる )。

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

6, 7, 8 列目の和が 4 列目になるので 8 列目を除去。同様の考え方から 11, 14 列目も除去。

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0	\|

6, 8, 10 列目の和が 2 列目になるので 10 列目を除去。同様の考え方から 11 列目も除去。

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0	\|

X	=	\|	1_np,	D,	F₁,	D,	0_np,p-1,	...	0_np,p-1	\|
		\|	1_np,	D,	F₂,	0_np,p-1,	D,	...	0_np,p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_np,	D,	F_q-1,	0_np,p-1,	0_np,p-1,	...	D	\|
		\|	1_np,	D,	0_np,q-1,	0_np,p-1,	0_np,p-1,	...	0_np,p-1	\|

但し、D は p 列目が除外された np 行 p - 1 列の行列であり、F_k も q 列目が除外された np 行 q - 1 列の行列です。X^TX は

X^TX	=	\|	q1_np^T1_np,	q1_np^TD,	Σ_k{1→q-1}( 1_np^TF_k ),	1_np^TD,	1_np^TD,	...	1_np^TD	\|
		\|	qD^T1_np,	qD^TD,	Σ_k{1→q-1}( D^TF_k ),	D^TD,	D^TD,	...	D^TD	\|
		\|	Σ_k{1→q-1}( F_k^T1_np ),	Σ_k{1→q-1}( F_k^TD ),	Σ_k{1→q-1}( F_k^TF_k ),	F₁^TD,	F₂^TD,	...	0_q-1,p-1	\|
		\|	D^T1_np,	D^TD,	D^TF₁,	D^TD,	0_q-1,p-1,	...	0_q-1,p-1	\|
		\|	D^T1_np,	D^TD,	D^TF₂,	0_q-1,p-1,	D^TD,	...	0_q-1,p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	D^T1_np,	D^TD,	0_q-1,p-1,	0_q-1,p-1,	0_q-1,p-1,	...	D^TD	\|

q1_np^T1_np = npq

1_np^TD = n1_p-1^T

1_np^TF_k = ( 0, ... np, ... 0 ) ( k 番目の要素を np、その他を 0 とする q - 1 次元のベクトル ) より

Σ_k{1→q-1}( 1_np^TF_k ) = np1_q-1^T

D^TD	=	\|	1_n^T1_n,	0,	...	0	\|
		\|	0,	1_n^T1_n,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	1_n^T1_n	\|
	=	\|	n,	0,	...	0	\| = nE_p-1 ( p - 1 の大きさの単位行列 x n )
		\|	0,	n,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	n	\|

					k 列目
D^TF_k	=	\|	0,	...	1_n^T1_n,	...	0	\|
		\|	0,	...	1_n^T1_n,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	1_n^T1_n,	...	0	\|
	=	\|	0,	...	n,	...	0	\|
		\|	0,	...	n,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	n,	...	0	\| より

Σ_k{1→q-1}( D^TF_k )	=	\|	n,	n,	...	n	\|
		\|	n,	n,	...	n	\|
		\|	:	:	...	:	\|
		\|	n,	n,	...	n	\|

					k 列目
F_k^TF_k	=	\|	0,	...	0,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	p1_n^T1_n,	...	0	\| k 行目
		\|	:	...	:	...	:	\|
		\|	0,	...	0,	...	0	\|
	=	\|	0,	...	0,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	np,	...	0	\| k 行目
		\|	:	...	:	...	:	\|
		\|	0,	...	0,	...	0	\| より

Σ_k{1→q-1}( F_k^TF_k )	=	\|	np,	0,	...	0	\| = npE_q-1 ( q - 1 の大きさの単位行列 x np )
		\|	0,	np,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	np	\|

X^TX	= n	\|	pq,	q1_p-1^T,	p1_q-1^T	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|
		\|	q1_p-1,	qE_p-1,	N^T	E_p-1,	E_p-1,	...	E_p-1	\|
		\|	p1_q-1,	N,	pE_q-1	F₁^T,	F₂^T,	...	F_q-1^T	\|
		\|	1_p-1,	E_p-1,	F₁	E_p-1,	0_p-1,	...	0_p-1	\|
		\|	1_p-1,	E_p-1,	F₂	0_p-1,	E_p-1,	...	0_p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_p-1,	E_p-1,	F_q-1	0_p-1,	0_p-1,	...	E_p-1	\|

											<	( p - 1 )( q - 1 ) 列											>
		<1 列>	<	p-1 列		>	<	q-1 列		>	<	p-1 列		>	<	p-1 列		>		<	p-1 列		>
= n	\|	pq,	q,	q,	...	q,	p,	p,	...	p,	1,	1,	...	1,	1,	1,	...	1,	...	1,	1,	...	1	\|	1 行

	\|	q,	q,	0,	...	0,	1,	1,	...	1,	1,	0,	...	0,	1,	0,	...	0,	...	1,	0,	...	0	\|	p-1 行
	\|	q,	0,	q,	...	0,	1,	1,	...	1,	0,	1,	...	0,	0,	1,	...	0,	...	0,	1,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	q,	0,	0,	...	q,	1,	1,	...	1,	0,	0,	...	1,	0,	0,	...	1,	...	0,	0,	...	1	\|

	\|	p,	1,	1,	...	1,	p,	0,	...	0,	1,	1,	...	1,	0,	0,	...	0,	...	0,	0,	...	0	\|	q-1 行
	\|	p,	1,	1,	...	1,	0,	p,	...	0,	0,	0,	...	0,	1,	1,	...	1,	...	0,	0,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	p,	1,	1,	...	1,	0,	0,	...	p,	0,	0,	...	0,	0,	0,	...	0,	...	1,	1,	...	1	\|

	\|	1,	1,	0,	...	0,	1,	0,	...	0,	1,	0,	...	0,	0,	0,	...	0,	...	0,	0,	...	0	\|	p-1 行
	\|	1,	0,	1,	...	0,	1,	0,	...	0,	0,	1,	...	0,	0,	0,	...	0,	...	0,	0,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	0,	0,	...	1,	1,	0,	...	0,	0,	0,	...	1,	0,	0,	...	0,	...	0,	0,	...	0	\|

	\|	1,	1,	0,	...	0,	0,	1,	...	0,	0,	0,	...	0,	1,	0,	...	0,	...	0,	0,	...	0	\|	p-1 行
	\|	1,	0,	1,	...	0,	0,	1,	...	0,	0,	0,	...	0,	0,	1,	...	0,	...	0,	0,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	0,	0,	...	1,	0,	1,	...	0,	0,	0,	...	0,	0,	0,	...	1,	...	0,	0,	...	0	\|

	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|

	\|	1,	1,	0,	...	0,	0,	0,	...	1,	0,	0,	...	0,	0,	0,	...	0,	...	1,	0,	...	0	\|	p-1 行
	\|	1,	0,	1,	...	0,	0,	0,	...	1,	0,	0,	...	0,	0,	0,	...	0,	...	0,	1,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	0,	0,	...	1,	0,	0,	...	1,	0,	0,	...	0,	0,	0,	...	0,	...	0,	0,	...	1	\|

と求められます。但し、N は全要素が 1 の q - 1 行 p - 1 列行列とします。X^TX の逆行列 ( X^TX )^-1 は以下のような結果になります (補足 2)。

( X^TX )^-1	= (1/n)	\|	1,	-1_p-1^T,	-1_q-1^T,	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|
		\|	-1_p-1,	N_p-1 + E_p-1,	N^T,	-( N_p-1 + E_p-1 ),	-( N_p-1 + E_p-1 ),	...	-( N_p-1 + E_p-1 )	\|
		\|	-1_q-1,	N,	N_q-1 + E_q-1,	-( N + F₁^T ),	-( N + F₂^T ),	...	-( N + F_q-1^T )	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₁ )	2( N_p-1 + E_p-1 ),	N_p-1 + E_p-1,	...	N_p-1 + E_p-1	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₂ )	N_p-1 + E_p-1,	2( N_p-1 + E_p-1 ),	...	N_p-1 + E_p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F_q-1 )	N_p-1 + E_p-1,	N_p-1 + E_p-1,	...	2( N_p-1 + E_p-1 )	\|

											<	( p - 1 )( q - 1 ) 列											>
		<1 列>	<	p-1 列		>	<	q-1 列		>	<	p-1 列		>	<	p-1 列		>		<	p-1 列		>
= (1/n)	\|	1,	-1,	-1,	...	-1,	-1,	-1,	...	-1,	1,	1,	...	1,	1,	1,	...	1,	...	1,	1,	...	1	\|	1 行

	\|	-1,	2,	1,	...	1,	1,	1,	...	1,	-2,	-1,	...	-1,	-2,	-1,	...	-1,	...	-2,	-1,	...	-1	\|	p-1 行
	\|	-1,	1,	2,	...	1,	1,	1,	...	1,	-1,	-2,	...	-1,	-1,	-2,	...	-1,	...	-1,	-2,	...	-1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	-1,	1,	1,	...	2,	1,	1,	...	1,	-1,	-1,	...	-2,	-1,	-1,	...	-2,	...	-1,	-1,	...	-2	\|

	\|	-1,	1,	1,	...	1,	2,	1,	...	1,	-2,	-2,	...	-2,	-1,	-1,	...	-1,	...	-1,	-1,	...	-1	\|	q-1 行
	\|	-1,	1,	1,	...	1,	1,	2,	...	1,	-1,	-1,	...	-1,	-2,	-2,	...	-2,	...	-1,	-1,	...	-1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	-1,	1,	1,	...	1,	1,	1,	...	2,	-1,	-1,	...	-1,	-1,	-1,	...	-1,	...	-2,	-2,	...	-2	\|

	\|	1,	-2,	-1,	...	-1,	-2,	-1,	...	-1,	4,	2,	...	2,	2,	1,	...	1,	...	2,	1,	...	1	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-2,	-1,	...	-1,	2,	4,	...	2,	1,	2,	...	1,	...	1,	2,	...	1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-2,	-1,	...	-1,	2,	2,	...	4,	1,	1,	...	2,	...	1,	1,	...	2	\|

	\|	1,	-2,	-1,	...	-1,	-1,	-2,	...	-1,	2,	1,	...	1,	4,	2,	...	2,	...	2,	1,	...	1	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-1,	-2,	...	-1,	1,	2,	...	1,	2,	4,	...	2,	...	1,	2,	...	1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-1,	-2,	...	-1,	1,	1,	...	2,	2,	2,	...	4,	...	1,	1,	...	2	\|

	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|

	\|	1,	-2,	-1,	...	-1,	-1,	-1,	...	-2,	2,	1,	...	1,	2,	1,	...	1,	...	4,	2,	...	2	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-1,	-1,	...	-2,	1,	2,	...	1,	1,	2,	...	1,	...	2,	4,	...	2	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-1,	-1,	...	-2,	1,	1,	...	2,	1,	1,	...	2,	...	2,	2,	...	4	\|

X^Ty	=	\|	1_np^T,	1_np^T,	...	1_np^T,	1_np^T	\|\|	y₁	\|
		\|	D^T,	D^T,	...	D^T	D^T	\|\|	y₂	\|
		\|	F₁^T,	F₂^T,	...	F_q-1^T	0_q-1,np	\|\|	:	\|
		\|	D^T,	0_p-1,np,	...	0_p-1,np	0_p-1,np	\|\|	y_q-1	\|
		\|	0_p-1,np,	D^T,	...	0_p-1,np	0_p-1,np	\|\|	y_q	\|
		\|	:	:	...	:	:	\|
		\|	0_p-1,np,	0_p-1,np,	...	D^T	0_p-1,np	\|

から求めることができます。但し、y_k は事象 B の添字が k である np 個の要素からなるベクトル ( y_1,1,k, y_2,1,k, ... y_n,1,k, y_1,2,k, ... y_n,p,k )^T を表します。X^T の個々の行ブロックに対する計算結果は

Σ_k{1→q}( 1_np^Ty_k ) = Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k ) ) ) ≡ S_y

Σ_k{1→q}( D^Ty_k )	=	( Σ_k{1→q}( Σ_i{1→n}( y_i,1,k ) ), Σ_k{1→q}( Σ_i{1→n}( y_i,2,k ) ), ... Σ_k{1→q}( Σ_i{1→n}( y_i,p-1,k ) ) )^T
	≡	( S_A,1, S_A,2, ... S_A,p-1 )^T ≡ S_A

Σ_k{1→q}( F_k^Ty_k )	=	( Σ_j{1→p}( Σ_i{1→n}( y_i,j,1 ) ), Σ_j{1→p}( Σ_i{1→n}( y_i,j,2 ) ), ... Σ_j{1→p}( Σ_i{1→n}( y_i,j,q-1 ) ) )^T
	≡	( S_B,1, S_B,2, ... S_B,q-1 )^T ≡ S_B

D^Ty_k	=	( Σ_i{1→n}( y_i,1,k ), Σ_i{1→n}( y_i,2,k ), ... Σ_i{1→n}( y_i,p-1,k ) )^T
	≡	( s_1,k, s_2,k, ... s_p-1,k )^T ≡ s_k

a = ( X^TX )^-1X^Ty

= (1/n)

-1_p-1^T,

-1_q-1^T,

1_p-1^T,

...

1_p-1^T

S_y

-1_p-1,

N_p-1 + E_p-1,

N^T,

-( N_p-1 + E_p-1 ),

...

-( N_p-1 + E_p-1 )

S_A

-1_q-1,

N_q-1 + E_q-1,

-( N + F₁^T ),

-( N + F₂^T ),

...

-( N + F_q-1^T )

S_B

1_p-1,

-( N_p-1 + E_p-1 ),

-( N^T + F₁ ),

2( N_p-1 + E_p-1 ),

N_p-1 + E_p-1,

...

N_p-1 + E_p-1

s₁

1_p-1,

-( N_p-1 + E_p-1 ),

-( N^T + F₂ ),

N_p-1 + E_p-1,

2( N_p-1 + E_p-1 ),

...

N_p-1 + E_p-1

s₂

...

1_p-1,

-( N_p-1 + E_p-1 ),

-( N^T + F_q-1 ),

N_p-1 + E_p-1,

...

2( N_p-1 + E_p-1 )

s_q-1

から回帰係数の最尤推定量 a を得ることができます。ここで、

S_y	=	Σ_j{1→p-1}( S_A,j ) + S_A,p
	=	Σ_k{1→q-1}( S_B,k ) + S_B,q
	=	Σ_k{1→q-1}( Σ_j{1→p-1}( s_j,k ) ) + S_A,p + S_B,q - s_p,q

が成り立つことを利用して一行めのブロック単位で計算すると

	( 1, -1_p-1^T, -1_q-1^T, 1_p-1^T, 1_p-1^T, ... 1_p-1^T )( S_y, S_A^T, S_B^T, s₁^T, s₂^T, ... s_q-1^T )^T / n
=	[ S_y - 1_p-1^TS_A - 1_q-1^TS_B + Σ_k{1→q-1}( 1_p-1^Ts_k ) ] / n
=	[ S_y - Σ_j{1→p-1}( S_A,j ) - Σ_k{1→q-1}( S_B,k ) + Σ_k{1→q-1}( Σ_j{1→p-1}( s_j,k ) ) ] / n
=	[ S_y - ( S_y - S_A,p ) - ( S_y - S_B,q ) + ( S_y - S_A,p - S_B,q + s_p,q ) ] / n
=	s_p,q / n ≡ m_p,q

となり、これは、事象 ( A, B ) が ( p, q ) に属する従属変数の平均を意味します。二行目は

	( -1_p-1, N_p-1 + E_p-1, N^T, -( N_p-1 + E_p-1 ), -( N_p-1 + E_p-1 ), ... -( N_p-1 + E_p-1 ) )( S_y, S_A^T, S_B^T, s₁^T, s₂^T, ... s_q-1^T )^T / n
=	( [ -S_y1_p-1 + ( N_p-1 + E_p-1 )S_A + N^TS_B - Σ_k{1→q-1}( ( N_p-1 + E_p-1 )s_k ) ] / n )

で計算することができて、この結果もベクトルになります。その r 行目 ( 1 ≤ r ≤ p - 1 ) は

	[ -S_y + Σ_j{1→p-1}( S_A,j ) + S_A,r + Σ_k{1→q-1}( S_B,k ) - Σ_k{1→q-1}( s_r,k + Σ_j{1→p-1}( s_j,k ) ) ] / n
=	{ -S_y + ( S_y - S_A,p ) + S_A,r + ( S_y - S_B,q ) - [ ( S_A,r - s_r,q ) + ( S_y - S_A,p - S_B,q + s_p,q ) ] } / n
=	( s_r,q - s_p,q ) / n ≡ m_r,q - m_p,q

なので、これは事象 B を q に固定したときの、事象 A に対して r 番目と p 番目の平均差になります。次の三行目は

	( -1_q-1, N, N_q-1 + E_q-1, -( N + F₁^T ), -( N + F₂^T ), ... -( N + F_q-1^T ) )( S_y, S_A^T, S_B^T, s₁^T, s₂^T, ... s_q-1^T )^T / n
=	[ -S_y1_q-1 + NS_A + ( N_q-1 + E_q-1 )S_B - Σ_k{1→q-1}( ( N + F_k^T )s_k ) ] / n

=	[ -S_y + Σ_j{1→p-1}( S_A,j ) + Σ_k{1→q-1}( S_B,k ) + S_B,r - Σ_k{1→q-1}( Σ_j{1→p-1}( s_j,k ) ) - Σ_j{1→p-1}( s_j,r ) ] / n
=	[ -S_y + ( S_y - S_A,p ) + ( S_y - S_B,q ) + S_B,r - ( S_y - S_A,p - S_B,q + s_p,q ) - ( S_B,r - s_p,r ) ] / n )
=	( s_p,r - s_p,q ) / n ≡ m_p,r - m_p,q

となって、事象 A を p に固定したときの、事象 B に対して r 番目と q 番目の平均差になります。四行目以降は q - 1 行のブロックに分かれているので、その R 行目について計算すると、

		1 番目	...	R 番目	...	q-1 番目
	( 1_p-1, -( N_p-1 + E_p-1 ), -( N^T + F_R ),	N_p-1 + E_p-1,	...	2( N_p-1 + E_p-1 ),	...	N_p-1 + E_p-1	)( S_y, S_A^T, S_B^T, s₁^T, ... s_R^T, ... s_q-1^T )^T / n
=	[ S_y1_p-1 - ( N_p-1 + E_p-1 )S_A - ( N^T + F_R )S_B + ( N_p-1 + E_p-1 )s₁ + ... + 2( N_p-1 + E_p-1 )s_R + ... + ( N_p-1 + E_p-1 )s_q-1 ] / n
=	[ S_y1_p-1 - ( N_p-1 + E_p-1 )S_A - ( N^T + F_R )S_B + Σ_k{1→q-1}( ( N_p-1 + E_p-1 )s_k ) + ( N_p-1 + E_p-1 )s_R ] / n

	[ S_y - Σ_j{1→p-1}( S_A,j ) - S_A,r - Σ_k{1→q-1}( S_B,k ) - S_B,R + Σ_k{1→q-1}( Σ_j{1→p-1}( s_j,k ) + s_r,k ) + Σ_j{1→p-1}( s_j,R ) + s_r,R ] / n
=	[ S_y - ( S_y - S_A,p ) - S_A,r - ( S_y - S_B,q ) - S_B,R + ( S_y - S_A,p - S_B,q + s_p,q ) + ( S_A,r - s_r,q ) + ( S_B,R - s_p,R ) + s_r,R ] / n
=	[ ( s_r,R - s_p,R ) - ( s_r,q - s_p,q ) ] / n
=	( m_r,R - m_p,R ) - ( m_r,q - m_p,q )

となります。m_r,R - m_p,R は事象 B を R に固定したときの、m_r,q - m_p,q は事象 B を q に固定したときの、事象 A が r 番目と p 番目の場合の平均差です。事象 A の中のある二つの要素で平均差を求めた時、それが事象 B の選び方でどの程度変化するのかをこの値で見ることができます。

以上の結果から、a = ( m, a₁, ... a_p, b₁, ... b_q, g_1,1, g_2,1, ... g_p,1, g_1,2, ... g_p,q )^T の各係数は以下のように表されることになります。

m = m_p,q

a_j = m_j,q - m_p,q

b_k = m_p,k - m_p,q

g_j,k = ( m_j,k - m_p,k ) - ( m_j,q - m_p,q )

m + a_j + b_k + g_j,k	=	m_p,q + ( m_j,q - m_p,q ) + ( m_p,k - m_p,q ) + [ ( m_j,k - m_p,k ) - ( m_j,q - m_p,q ) ]
	=	m_j,k

D_E	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( ( y_i,j,k - x_{np(k-1)+n(j-1)+i}^Ta )² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( [ y_i,j,k - ( m + a_j + b_k + g_j,k ) ]² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( ( y_i,j,k - m_j,k )² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² - 2y_i,j,km_j,k + m_j,k² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) - 2m_j,kΣ_i{1→n}( y_i,j,k ) + nm_j,k² ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) - nm_j,k² ) ) / σ²

となり、これは「反復のある二元配置分散分析法」での「誤差変動 S_E」に相当します (*1-5)。

次に、「交互作用効果 (Interaction Effect)」γ_j,k を無視した以下のような縮小モデルを考えます。

y_i,j,k = μ + α_j + β_k + ε_i,j,k

このときのデザイン行列 X は、交互作用効果を含んだモデルに対してその係数に対する変数部分を除外したものに等しくなります。すなわち、

X	=	\|	1_np,	D,	F₁	\|
		\|	1_np,	D,	F₂	\|
		\|	:	:	:	\|
		\|	1_np,	D,	F_q-1	\|
		\|	1_np,	D,	0_np,q-1	\|

X^TX	= n	\|	pq,	q1_p-1^T,	p1_q-1^T	\|
		\|	q1_p-1,	qE_p-1,	N^T	\|
		\|	p1_q-1,	N,	pE_q-1	\|

( X^TX )^-1	= (1/n)	\|	( p + q - 1 ) / pq,	-1_p-1^T / q,	-1_q-1^T / p	\|
		\|	-1_p-1 / q,	( 1 / q )( N_p-1 + E_p-1 ),	0_p-1,q-1	\|
		\|	-1_q-1 / p,	0_q-1,p-1,	( 1 / p )( N_q-1 + E_q-1 )	\|

です (補足 3)。α = ( μ, α₁, ... α_p, β₁, ... β_q )^T の最尤推定量 a = ( m, a₁, ... a_p, b₁, ... b_q )^T は次の式を解くことで求められます。

a = ( X^TX )^-1X^Ty	= (1/n)	\|	( p + q - 1 ) / pq,	-1_p-1^T / q,	-1_q-1^T / p	\|\|	S_y	\|
		\|	-1_p-1 / q,	( 1 / q )( N_p-1 + E_p-1 ),	0_p-1,q-1	\|\|	S_A	\|
		\|	-1_q-1 / p,	0_q-1,p-1,	( 1 / p )( N_q-1 + E_q-1 )	\|\|	S_B	\|

	( ( p + q - 1 ) / pq, -1_p-1^T / q, -1_q-1^T / p )( S_y, S_A^T, S_B^T )^T / n
=	( [ ( p + q - 1 ) / pq ]S_y - 1_p-1^TS_A / q - 1_q-1^TS_B / p ) / n
=	[ ( 1 / q + 1 / p - 1 / pq )S_y - Σ_j{1→p-1}( S_A,j ) / q - Σ_k{1→q-1}( S_B,k ) / p ] / n
=	[ ( 1 / q + 1 / p - 1 / pq )S_y - ( S_y - S_A,p ) / q - ( S_y - S_B,q ) / p ] / n
=	S_A,p / nq + S_B,q / np - S_y / npq
≡	m_A,p + m_B,q - m_y

となります。ここで m_A,p, m_B,q, m_y はそれぞれ、p 番目の事象 A、q 番目の事象 B、全従属変数の平均を意味します。二行目のブロックは

	( -1_p-1 / q, ( N_p-1 + E_p-1 ) / q, 0_p-1,q-1 )( S_y, S_A^T, S_B^T )^T / n
=	[ -S_y1_p-1 + ( N_p-1 + E_p-1 )S_A ] / nq

	[ -S_y + Σ_j{1→p-1}( S_A,j ) + S_A,r ] / nq
=	[ -S_y + ( S_y - S_A,p ) + S_A,r ] / nq
=	m_A,r - m_A,p

	( -1_q-1 / p, 0_q-1,p-1, ( N_q-1 + E_q-1 ) / p )( S_y, S_A^T, S_B^T )^T / n
=	[ -S_y1_q-1 + ( N_q-1 + E_q-1 )S_B ] / np

=	[ -S_y + Σ_k{1→q-1}( S_B,k ) + S_B,r ] / np
=	[ -S_y + ( S_y - S_B,q ) + S_B,r ] / np
=	m_B,r - m_B,q

m = m_A,p + m_B,q - m_y

a_j = m_A,j - m_A,p

b_k = m_B,k - m_B,q

m + a_j + b_k = m_A,j + m_B,k - m_y

D_AB	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( [ y_i,j,k - ( m + a_j + b_k ) ]² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( [ y_i,j,k - ( m_A,j + m_B,k - m_y ) ]² ) ) ) / σ²

となります。交互作用効果を含んだモデルとの対数尤度比 D_AB - D_E は

D_AB - D_E	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( [ y_i,j,k - ( m_A,j + m_B,k - m_y ) ]² ) - [ Σ_i{1→n}( y_i,j,k² ) - nm_j,k² ] ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( -2( m_A,j + m_B,k - m_y )Σ_i{1→n}( y_i,j,k ) + n( m_A,j + m_B,k - m_y )² + nm_j,k² ) ) / σ²
	=	nΣ_k{1→q}( Σ_j{1→p}( m_j,k² - 2( m_A,j + m_B,k - m_y )m_j,k + ( m_A,j + m_B,k - m_y )² ) ) / σ²
	=	nΣ_k{1→q}( Σ_j{1→p}( ( m_j,k - m_A,j - m_B,k + m_y )² ) ) / σ²

となって、これは二元配置分散分析における「交互作用 S_RC」に該当します (*1-5)。

さらに事象 B に関する変動分 β_k まで除外したモデル式とデザイン行列は

y_i,j,k = μ + α_j + ε_i,j,k

X	=	\|	1_np,	D	\|
		\|	1_np,	D	\|
		\|	:	:	\|
		\|	1_np,	D	\|

X^TX	= n	\|	pq,	q1_p-1^T	\|
		\|	q1_p-1,	qE_p-1	\|

( X^TX )^-1	= (1/nq)	\|	1,	-1_p-1^T	\|
		\|	-1_p-1,	N_p-1 + E_p-1	\|

です。α = ( μ, α₁, ... α_p )^T の最尤推定量 a = ( m, a₁, ... a_p )^T は

a = ( X^TX )^-1X^Ty	= (1/nq)	\|	1,	-1_p-1^T	\|\|	S_y	\|
		\|	-1_p-1,	N_p-1 + E_p-1	\|\|	S_A	\|

m = [ S_y - Σ_j{1→p-1}( S_A,j ) ] / nq

a_r = [ -S_y + Σ_j{1→p-1}( S_A,j ) + S_A,r ] / nq ( 2 ≤ r ≤ p )

m + a_j = S_A,j / nq = m_A,j ( 2 ≤ j ≤ p )

D_A	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( [ y_i,j,k - ( m + a_j ) ]² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( ( y_i,j,k - m_A,j )² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² - 2y_i,j,km_A,j + m_A,j² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) - 2m_A,jΣ_i{1→n}( y_i,j,k ) + nm_A,j² ) ) / σ²
	=	[ Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) ) ) - 2Σ_j{1→p}( m_A,jΣ_k{1→q}( s_j,k ) ) + nqΣ_j{1→p}( m_A,j² ) ] / σ²
	=	[ Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) ) ) - 2Σ_j{1→p}( m_A,jS_A,j ) ) + nqΣ_j{1→p}( m_A,j² ) ] / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) - nm_A,j² ) ) / σ²

となります。最小モデルに対する対数尤度統計量を D_T とすると、

D_T	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k - m_y )² ) ) ) / σ²
	=	Σ_k{1→q}( Σ_j{1→p}( Σ_i{1→n}( y_i,j,k² ) - nm_y² ) ) / σ²

D_T - D_A	=	Σ_k{1→q}( Σ_j{1→p}( [ Σ_i{1→n}( y_i,j,k² ) - nm_y² ] - [ Σ_i{1→n}( y_i,j,k² ) - nm_A,j² ] ) ) / σ²
	=	nΣ_k{1→q}( Σ_j{1→p}( m_A,j² - m_y² ) ) / σ²
	=	nqΣ_j{1→p}( m_A,j² - m_y² ) ) / σ²

となり、これは事象 A での ( 事象 B の差異を考慮しない ) 変動を表します。同様に、事象 B での ( 事象 A の差異を考慮しない ) 変動は

npΣ_k{1→q}( m_B,k² - m_y² ) ) / σ²

であり、これらは二元配置分散分析における「列間・行間変動 S_R, S_C」を意味します。

分散分析は、一般化線形モデルにおける尤度比を使った検定と同等であることが、これらの結果からわかります。その内容についても比較的理解がしやすく、それが広く利用されている理由の一つなのではないかと思います。但し、連結関数は恒等関数であり、利用する確率分布も正規分布であることから、それ以外のモデルを考慮しなければならない場合は分散分析をそのまま利用することはできません。

2) ポアソン回帰 (Poisson Regression)

分散分析の場合、独立変数が事象をダミー変数で表したものであるのに対し従属変数は連続値になるのが通常で、モデル式との誤差が平均ゼロの正規分布に従うと仮定した上で検定を行います。また、正規分布の分散は、どの従属変数に対しても一定であるとみなします。しかし、従属変数が「計数データ (Count Data)」または「度数データ (Frequency Data)」と呼ばれる、ある事象の発生回数であった場合、正規分布よりも「ポアソン分布 (Poisson Distribution)」がよく用いられます (*2-1)。ポアソン分布は、

f(μ|y) = μ^ye^-μ / y!

で表される離散確率密度関数で、一母数指数型分布族なので一般化線形モデルを適用することができます。確率変数の y は事象の発生した回数を表し、その平均発生回数が μ になります。すなわち、E[y] = μ です。また、ポアソン分布は分散が平均と等しいという特徴があります。よって、V[y] = μ ということになります。μ は割合や比率で表されます。例えば、一定時間内に発生する事象の平均回数や、一定の試行回数あたりで発生する事象の平均回数などで表現されます。特に後者の場合、試行回数を n で表せば μ = nθ とすることができて、さらに

μ = ne^xα より

log μ = log n + xα

となって、連結関数は対数関数になります。なお、θ を指数関数にすることで、任意の x に対して μ は正の数であることが保証されます。

ある独立変数 x_i に対してその従属変数 y_i の期待値が E[y_i] = μ_i = n_iθ_i であるとします。これらを上式にあてはめると

g( μ_i ) = log( μ_i / n_i ) = x_iα

g'( μ_i ) = [ 1 / ( μ_i / n_i ) ]・( 1 / n_i ) = 1 / μ_i

であり、また V[y_i] = μ_i なので、1 / V[y_i]g'(μ_i)² = μ_i となって、これがフィッシャー情報行列 Φ = X^TWX における対角行列 W の対角成分となります。すでにポアソン分布用のクラス ExpFamily_Poisson と対数関数を連結関数とした LogFunc クラスは「確率・統計 (18) 一般化線形モデル」において用意されているので、これらを使ってスコア法 ScoringMethod 関数を呼び出せば処理をすることができます。

l( μ \| y )	=	log( Π_i{1→N}( μ_i^y_ie^-μ_i / y_i! ) )
	=	Σ_i{1→N}( y_ilog μ_i - μ_i - log y_i! )

となります。ここで、N は独立変数・従属変数の数を表します。飽和モデルの場合は μ を使ったモデル式を考えればよいので

∂l / ∂μ_i = y_i / μ_i - 1 = 0

より μ_i = y_i、すなわち推定量は実測値と等しくなります。あるモデル式から得られた当てはめ値を m_i とすれば、

l( m | y ) = Σ_i{1→N}( y_ilog m_i - m_i - log y_i! )

なので、飽和モデルの対数尤度関数の最大値を l( y | y ) で表した時、飽和モデルとの対数尤度統計量 D は

D = 2[ l( y \| y ) - l( m \| y ) ]	=	2[ Σ_i{1→N}( y_ilog y_i - y_i - log y_i! ) - Σ_i{1→N}( y_ilog m_i - m_i - log y_i! ) ]
	=	2Σ_i{1→N}( y_ilog( y_i / m_i ) - ( y_i - m_i ) )

u_j = Σ_i{1→N}( ( y_i - μ_i )x_ij / V[y_i]g'(μ_i) )

で表され、指数型分布属をポアソン分布、連結関数を対数関数とした場合、

V[y_i] = μ_i

g'(μ_i) = 1 / μ_i

u_j = Σ_i{1→N}( ( y_i - μ_i )x_ij )

となるので、特に定数項に対するスコア統計量は Σ_i{1→N}( y_i - μ_i ) です。スコア統計量がゼロになるときの各係数が最尤推定量になることから、特に定数項に対しては

Σ_i{1→N}( y_i - m_i ) = 0

が成り立ちます。従って、(モデル式が定数項を含むなら) 対数尤度統計量 D は

D = 2Σ_i{1→N}( y_ilog( y_i / m_i ) )

と単純化することができます。なお、ロジスティック回帰の時と同様に「逸脱度残差」を定義すると、

d_i = SIGN * { 2 * [ y_ilog( y_i / m_i ) - ( y_i - m_i ) ] }^1/2

となります。但し、SIGN は y_i - m_i の符号を意味します ( 補足 5 )。

最小モデルに対する対数尤度関数は、θ_i が全て共通の値 θ であると考えて

l( θ | y ) = Σ_i{1→N}( y_ilog n_iθ - n_iθ - log y_i! )

Σ_i{1→N}( y_i / θ - n_i ) = 0 より

θ = Σ_i{1→N}( y_i ) / Σ_i{1→N}( n_i )

となるので、θ の最尤推定量は y_i の和と n_i の和の比率となります。これを t として、最小モデルの対数尤度関数の最大値を l( t | y ) で表した時、最小モデルとの対数尤度統計量 ( = 「尤度比カイ二乗統計量 C」 ) は、

C = 2[ l( m \| y ) - l( t \| y ) ]	=	2Σ_i{1→N}( y_ilog m_i - m_i - log y_i! ) - 2Σ_i{1→N}( y_ilog n_it - n_it - log y_i! )
	=	2Σ_i{1→N}( y_ilog( m_i / n_it ) - ( m_i - n_it ) )

Σ_i{1→N}( m_i - n_it )	=	Σ_i{1→N}( m_i ) - Σ_i{1→N}( n_i )Σ_i{1→N}( y_i ) / Σ_i{1→N}( n_i )
	=	Σ_i{1→N}( m_i - y_i )

より、モデル式に定数項が含まれることを前提として、C は

C = 2Σ_i{1→N}( y_ilog( m_i / n_it ) )

と単純化することができます。また、「(マクファデンの) 擬似 R² 値」は [ l( t | y ) - l( m | y ) ] / l( t | y ) で求められるので、

R² = -Σ_i{1→N}( y_ilog( m_i / n_it ) ) / Σ_i{1→N}( y_ilog n_it - n_it - log y_i! )

線形重回帰分析において、観測した従属変数 y_i の予測値を y^_i とし、観測値 y_i と予測値 y^_i の平均をそれぞれ m_y, m^_y としたとき、

予測値の平方和 S^_y = Σ( ( y^_i - m^_y )² )
残差の平方和 S_ε = Σ( ( y_i - y^_i )² )
観測値の平方和 S_y = Σ( ( y_i - m_y )² )

とすれば、S^_y が自由度 p - 1 の、S_ε が自由度 N - p の、S_y が自由度 N - 1 のカイ二乗分布にそれぞれ従います。また、S^_y = S_y - S_ε が成り立つことから、S^_y は自由度 p - 1 のカイ二乗分布に従うことになります。しかし、正規分布には「局外パラメータ」として標準偏差 σ があり、そのままでは検定をすることができないので、v^_y ≡ S^_y / ( p - 1 ) と v_T = S_T / ( N - p ) の比率 F₀ が自由度 ( p - 1, N - p ) の F-分布に従うことを利用して検定を行います。重回帰分析が、誤差項が従う確率分布として利用する正規分布の分散を全て等しいと仮定しているのはそのためです。また、( S_y - S_ε ) / S_y = S^_y / S_y は「決定係数」または「寄与率」を表し、これは予測値がどの程度モデル式を反映しているかを見るための指標となるのでした。これらの結果は通常、分散分析表で表されます (*2-2)。
一般化線形モデルでは、S_ε が飽和モデルと線形モデルの間の対数尤度統計量 D と等しく、S_y が飽和モデルと最小モデルの間の対数尤度統計量 D₀ と等しいのでした。また、寄与率は ( D₀ - D ) / D₀ となり、これは線形モデルが飽和モデルと最小モデルの対数尤度差を指標としてどの程度の割合でモデル式を表せているのかを表していると考えることができます (*2-3)。

ポアソン回帰を行うためのサンプル・プログラムを示します。

/*
  LogFunc_Poisson : 対数連結関数(ポアソン回帰用)

  定数項として試行回数を加える。試行回数は各試行に応じて切り替える
*/
class LogFunc_Poisson : public LinkFunction_IF
{
  const vector<double>& ni_; // 試行回数

  // 各関数の利用回数(n_.size()の剰余)
  mutable unsigned int opCnt_;   // operator()()
  mutable unsigned int invfCnt_; // invf()

public:

  // コンストラクタ
  LogFunc_Poisson( const vector<double>& ni )
    : ni_( ni ), opCnt_( 0 ), invfCnt_( 0 ) {}

  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    double d = log( x ) - log( ni_[opCnt_] );
    opCnt_ = ( opCnt_ + 1 ) % ni_.size();
    return( d );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  { return( 1 / x ); }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  {
    double d = ni_[invfCnt_] * exp( y );
    invfCnt_ = ( invfCnt_ + 1 ) % ni_.size();
    return( d );
  }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Logarithm Function ( for Poisson Regression )" ); }
};

/*
  PoissonRegression : ポアソン回帰

  const vector<double>& y : 発生回数(従属変数)
  const vector<double>& ni : 試行回数
  const vector< vector<double> >& x : 独立変数
  vector<double>& a : 求めた係数
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool PoissonRegression( const vector<double>& y, const vector<double>& ni,
                        const vector< vector<double> >& x, vector<double>& a,
                        bool verbose, unsigned int maxCount, double threshold )
{
  cout << "*** Poisson Regression ***" << endl << endl;

  // NULL のチェック ( x, y, a は ScoringMethod で行われるので省略
  if ( ! NullCheck( ni, "Total Count ni" ) ) return( false );

  // x と ni のサイズ比較
  if ( ! SizeCheck( ni, "Total Count ni", x.size(), "Independent Variable x" ) )
    return( false );

  if ( verbose ) {
    PrintVector( "ni = ", ni );
    cout << endl;
  }

  // 連結関数の作成
  LogFunc_Poisson g( ni );

  // 指数型分布族(ポアソン分布)
  ExpFamily_Poisson pdf;

  return( ScoringMethod( x, y, a, pdf, g, verbose, maxCount, threshold ) );
}

総試行回数は各試行によって変わるので、連結関数に含まれる定数項 log n_i も切り替える必要があります。そこで、以前「(19) ロジスティック回帰 (Logistic Regression Model)」の章で紹介した、二項分布を利用したスコア法のサンプル・プログラムと同様の手法を使い、各メンバ関数が呼び出された時にカウンタを 1 ずつ増加させることで定数項が切り替わるようにします。スコア法の関数 ScoringMethod は、計算に必要な LinkFunction_IF のメンバ関数 operator(), df, invf をそれぞれ 1 ループに一度しか呼び出さないのでこのような手法が可能になりますが、前にも書いた通り実装に依存した方法なのであまりいいやり方とは言えません。内容は対数関数用の連結関数とほとんど変わりませんが、定数項を使うため以下のような式になります。

ξ_i = g( μ_i ) = log μ_i - log n_i

dξ_i / dμ_i = g'( μ_i ) = 1 / μ_i

μ_i = g^-1( ξ_i ) = n_ie^ξ_i

定数項を計算するための試行回数 n_i はインスタンス化をするときに配列のまま渡すようにしています。カウンタは opCnt_ と invfCnt_ で、それぞれ operator(), invf 用に利用します。df の計算には n_i は含まれないので専用のカウンタは不要です。これらのメンバ関数が呼び出されたら、計算後にカウンタを 1 増加させます。カウンタが n_i の個数以上にならないように、個数との剰余をカウンタに代入するようにしています。これによって、カウンタが循環するような形になります。

サンプル・プログラムを使って、参考文献にあるデータをポアソン回帰で解析してみたいと思います。下表は、イギリスの男性医師の冠動脈心疾患による死亡者数を 1951 年から 10 年間調査し、喫煙の有無と年齢ごとに分類した結果です。総数はその事象の調査対象者の人数を示しています。

グラフの死亡者数は、10 万人を総数とした時の比率で表されています。グラフから、死亡者数は年齢に対して二次のオーダーで増加しているように見えます。また、喫煙者と非喫煙者の間でわずかな差異が見られます。喫煙・非喫煙の有無を x_i1 = 1 (喫煙) または 0 (非喫煙)、年齢を x_i2 として 35～44 から 75～84 まで順に 1 から 5 まで連番で表します。また、死亡者数が年齢に対して二次のオーダーで増加していることを表すため、x_i3 = x_i2² として二次項を追加します。最後に、年齢と喫煙・非喫煙の相関を表す項として x_i4、定数項を x_i0 とします。x_i = ( x_i0, x_i1, x_i2, x_i3, x_i4 )^T は 5 次元のベクトルで、X = ( x₁, x₂, ... x₁₀ ) はデザイン行列です。X の具体的な値は

X	=	\|	1,	1,	1,	1,	1	\|
		\|	1,	1,	2,	4,	2	\|
		\|	1,	1,	3,	9,	3	\|
		\|	1,	1,	4,	16,	4	\|
		\|	1,	1,	5,	25,	5	\|
		\|	1,	0,	1,	1,	0	\|
		\|	1,	0,	2,	4,	0	\|
		\|	1,	0,	3,	9,	0	\|
		\|	1,	0,	4,	16,	0	\|
		\|	1,	0,	5,	25,	0	\|

であり、W の対角行列は μ_i なので、これらからフィッシャー情報行列 Φ を得ることができます。

これらを使ってスコア法により係数の最尤推定量を計算すると以下の結果が得られます ( verbose = true にして冗長モードで出力しています )。

*** Poisson Regression ***

ni = ( 52407, 43248, 28612, 12663, 5317, 18790, 10673, 5710, 2585, 1462 )

*** Scoring Method ***

Exponential Family of Distribution : Poisson Distribution
Link Function : Logarithm Function ( for Poisson Regression )

N = 10; p = 5

x = ( 1, 1, 1, 1, 1 )
    ( 1, 1, 2, 4, 2 )
    ( 1, 1, 3, 9, 3 )
    ( 1, 1, 4, 16, 4 )
    ( 1, 1, 5, 25, 5 )
    ( 1, 0, 1, 1, 0 )
    ( 1, 0, 2, 4, 0 )
    ( 1, 0, 3, 9, 0 )
    ( 1, 0, 4, 16, 0 )
    ( 1, 0, 5, 25, 0 )

y = ( 32, 104, 206, 186, 102, 2, 12, 28, 28, 31 )

----- cnt = 0 -----

mu = ( 32, 104, 206, 186, 102, 2, 12, 28, 28, 31 )
g'(mu) = ( 0.03125, 0.00961538, 0.00485437, 0.00537634, 0.00980392, 0.5, 0.0833333, 0.0357143, 0.0357143, 0.0322581 )
w = ( 0.03125, 0.00961538, 0.00485437, 0.00537634, 0.00980392, 0.5, 0.0833333, 0.0357143, 0.0357143, 0.0322581 )

Equation System :
(731)x0 + (630)x1 + (2489)x2 + (9353)x3 + (2112)x4 = -3563.5
(630)x0 + (630)x1 + (2112)x2 + (7828)x3 + (2112)x4 = -3068.66
(2489)x0 + (2112)x1 + (9353)x2 + (37601)x3 + (7828)x4 = -11428.8
(9353)x0 + (7828)x1 + (37601)x2 + (158753)x3 + (31080)x4 = -41233.4
(2112)x0 + (2112)x1 + (7828)x2 + (31080)x3 + (7828)x4 = -9696.75

Regression equation : y = -10.7355x0 + 1.39899x1 + 2.35751x2 + -0.19629x3 + -0.297892x4

----- cnt = 1 -----

mu = ( 29.7732, 106.94, 207.951, 182.681, 102.819, 3.54945, 11.8203, 25.0376, 30.3064, 30.9487 )
g'(mu) = ( 0.0335872, 0.00935107, 0.00480882, 0.00547401, 0.00972583, 0.281734, 0.0846, 0.03994, 0.0329964, 0.0323116 )
w = ( 0.0335872, 0.00935107, 0.00480882, 0.00547401, 0.00972583, 0.281734, 0.0846, 0.03994, 0.0329964, 0.0323116 )

Equation System :
(731.827)x0 + (630.164)x1 + (2490.6)x2 + (9357.26)x3 + (2112.33)x4 = -3570.14
(630.164)x0 + (630.164)x1 + (2112.33)x2 + (7822.47)x3 + (2112.33)x4 = -3068.6
(2490.6)x0 + (2112.33)x1 + (9357.26)x2 + (37626.3)x3 + (7822.47)x4 = -11440.4
(9357.26)x0 + (7822.47)x1 + (37626.3)x2 + (158935)x3 + (31044)x4 = -41256.5
(2112.33)x0 + (2112.33)x1 + (7822.47)x2 + (31044)x3 + (7822.47)x4 = -9704.44

Regression equation : y = -10.7912x0 + 1.44048x1 + 2.37629x2 + -0.197665x3 + -0.307437x4

----- cnt = 2 -----

mu = ( 29.5862, 106.813, 208.197, 182.827, 102.579, 3.41625, 11.5448, 24.7469, 30.2303, 31.0696 )
g'(mu) = ( 0.0337996, 0.00936217, 0.00480315, 0.00546966, 0.0097486, 0.292718, 0.0866192, 0.040409, 0.0330794, 0.0321858 )
w = ( 0.0337996, 0.00936217, 0.00480315, 0.00546966, 0.0097486, 0.292718, 0.0866192, 0.040409, 0.0330794, 0.0321858 )

Equation System :
(731.009)x0 + (630.001)x1 + (2489.02)x2 + (9353.05)x3 + (2112)x4 = -3564.38
(630.001)x0 + (630.001)x1 + (2112)x2 + (7820.31)x3 + (2112)x4 = -3067.31
(2489.02)x0 + (2112)x1 + (9353.05)x2 + (37611)x3 + (7820.31)x4 = -11430.2
(9353.05)x0 + (7820.31)x1 + (37611)x2 + (158868)x3 + (31028.7)x4 = -41232.5
(2112)x0 + (2112)x1 + (7820.31)x2 + (31028.7)x3 + (7820.31)x4 = -9702.87

Regression equation : y = -10.7918x0 + 1.44097x1 + 2.37648x2 + -0.197677x3 + -0.307548x4

Estimated regression equation : y = -10.7918x0 + 1.44097x1 + 2.37648x2 + -0.197677x3 + -0.307548x4
variance of a = ( 0.202529, 0.138504, 0.043238, 0.000748941, 0.0094156 )

独立変数の一つ x_j に対して、j 以外の事象が全て同一なときの、x_j = 1 での期待値と x_j = 0 での期待値の比率は「率比 (Rate Ratio)」と呼ばれ、それを RR で表すと

RR = E[y|x_j = 1] / E[y|x_j = 0] = e^α_j

となります。e^α_j は θ に対して乗法的に作用するので、α_j ではなく率比で表した方が意味を解釈しやすくなります。率比を計算した結果は以下のようになります。

年齢による差異を調整したときの喫煙の有無に対する死亡率は、率比から 4.22 倍と推定することができます。但し、x_i4 に対する率比が 1 より小さく、係数も負数であることから、その差は年齢とともに小さくなっていくことになります。

求めた係数から得られる期待値と残差を求めた結果は以下のようになります。

表の中で、Χ_i はピアソン残差、d_i は逸脱度残差、C_i は尤比度カイ二乗統計量の和の成分をそれぞれ表し、Χ² = Σ_i( Χ_i² ) がピアソン・カイ二乗統計量、D = Σ_i( d_i² ) が飽和モデルとの対数尤度統計量、C = Σ_i( C_i ) が尤比度カイ二乗統計量をそれぞれ表します。

Χ² も D も小さな値を示していることから、このモデル式への当てはめは非常によいことになります。飽和モデルのパラメータ数は 10 であり、仮定したモデル式のパラメータ数は 5 なので、Χ² と D は自由度 5 のカイ二乗分布に近似的に従うことになり、その上側 5% 点は 11.07 であることからはるかに小さな値を示しています。C はかなり大きな値であり、自由度 5 - 1 = 4 のカイ二乗分布に従うことから p 値を計算すると極端に小さな値となるので、定数項以外の係数がゼロであるという仮定は成り立ちません。擬似 R² 値は 0.943 であり、このことからもモデル式の当てはめが非常によいと判断できます。

3) 対数線形モデル (Log-linear Model)

ポアソン回帰では、各独立変数に対応する従属変数が従う確率分布のパラメータ μ_i が互いに制約を受けないことを前提としていました。しかし、全体の母数や各条件の母数にあらかじめ制約があるなどの理由でポアソン回帰をそのまま利用することができない場合があります。例えば、総試行回数が N 回で、その中である条件 A_i ( 1 ≤ i ≤ m ) に当てはまる回数が y_i 回であったと仮定します。これは、N 回の試行の中で A_i が y_i 回発生した時の確率密度を意味するので、A_i の発生確率を θ_i としたとき、多項分布

P_N,θ( y ) = N!Π_i{1→m}( θ_i^y_i / y_i! )

で表すことができます。ここで、μ_i = Nθ_i として各試行がポアソン分布

P_{μ_i}( y_i ) = μ_i^y_ie^-μ_i / y_i!

Π_i{1→m}( P_{μ_i}( y_i ) ) = exp( -Σ_i{1→m}( μ_i ) )Π_i{1→m}( μ_i^y_i / y_i! )

となりますが、総試行回数 Σ_i{1→m}( y_i ) が N となる確率分布はやはりポアソン分布に従い (補足 4)、

P_μ( N ) = μ^Ne^-μ / N!

と表すことができます。但し、μ = Σ_i{1→m}( μ_i ) です。従って、総試行回数が N という条件下での条件付き確率 P( y | N ) は

P( y \| N )	=	Π_i{1→m}( P_{μ_i}( y_i ) ) / P_μ( N )
	=	[ e^-μΠ_i{1→m}( μ_i^y_i / y_i! ) ] / ( μ^Ne^-μ / N! )
	=	N!Π_i{1→m}( μ_i^y_i / μ^y_iy_i! )
	=	N!Π_i{1→m}( ( μ_i / μ )^y_i / y_i! )

となって、μ_i / μ = θ_i とすれば多項分布と一致します。このとき、μ = Σ_i{1→m}( μ_i ) より Σ_i{1→m}( θ_i ) = 1 が成り立ちます。y_i の期待値 E[y_i] は

E[y_i] = Nθ_i = μ_i

g(μ_i) = log( μ_i / N ) = log( θ_i ) = x_iα

が成り立ち、これはポアソン回帰での連結関数において n_i = N とした場合に相当します。

ポアソン回帰では n_i も θ_i も未知数であり、これらを決定するのは連結関数です。ポアソン回帰の場合、連結関数は

g(μ_i) = log( μ_i / n_i ) = log μ_i - log n_i

g(μ_i) = log( μ_i / N ) = log μ_i - log N

となります。両者の差異は定数項のみであり、前者は各事象に対して定数項を決めた上で方程式を解くのに対し、後者は定数項を全て等しく、つまり連結関数は共通であるとした上で方程式を解いています。

デザイン行列 X を次のような形で表します。各行は一つの事象を表すので行数は m になります。また、端点制約により最後の m 番目の係数は定数項の係数の中に含めてしまい、列数は定数項を含め m 個となります。

X	=	\|	1,	1,	0,	...	0	\|
		\|	1,	0,	1,	...	0	\|
		\|	:	:	:	...	:	\|
		\|	1,	0,	0,	...	1	\|
		\|	1,	0,	0,	...	0	\|

これは、一元配置分散分析において y_i = 1 にした場合に相当します。回帰係数の最尤推定量を a = ( a₀, a₁, ... a_m-1 ) としたとき、i 行目の式は

a₀ + a_i ( i < m )
a₀ ( i = m )

a₀ = log ( μ_m / N ) = log θ_m
a_i = ( log μ_i - log N ) - ( log μ_m - log N ) = log ( μ_i / μ_m )

が成り立ち、率比 RR = e^a_i がそのまま μ_i の基準カテゴリ ( この場合 i = m ) との比率を表し、Ne^a_ie^a₀ が μ_i の推定量となります。

W は対角成分を μ_i とした対角行列で表すことができるので、スコア法における漸化式の左辺の係数行列 X^TWX は、

X^TWX

=	\|	1,	1,	...	1,	1	\|\|	μ₁,	0,	...	0,	0	\|\|	1,	1,	0,	...	0	\|
	\|	1,	0,	...	0,	0	\|\|	0,	μ₂,	...	0,	0	\|\|	1,	0,	1,	...	0	\|
	\|	0,	1,	...	0,	0	\|\|	:	:	...	:	:	\|\|	:	:	:	...	:	\|
	\|	:	:	...	:	:	\|\|	0,	0,	...	μ_m-1,	0	\|\|	1,	0,	0,	...	1,	\|
	\|	0,	0,	...	1,	0	\|\|	0,	0,	...	0,	μ_m	\|\|	1	0,	0,	...	0	\|

=	\|	μ₁,	μ₂,	...	μ_m-1,	μ_m	\|\|	1,	1,	0,	...	0	\|
	\|	μ₁,	0,	...	0,	0	\|\|	1,	0,	1,	...	0	\|
	\|	0,	μ₂,	...	0,	0	\|\|	:	:	:	...	:	\|
	\|	:	:	...	:	:	\|\|	1,	0,	0,	...	1,	\|
	\|	0,	0,	...	μ_m-1,	0	\|\|	1,	0,	0,	...	0	\|

=	\|	Σ_i{1→m}( μ_i ),	μ₁,	μ₂,	...	μ_m-1	\|
	\|	μ₁,	μ₁,	0,	...	0	\|
	\|	μ₂,	0,	μ₂,	...	0	\|
	\|	:	:	:	...	:	\|
	\|	μ_m-1,	0,	0,	...	μ_m-1	\|

z_i = g(μ_i) + ( y_i - μ_i )g'(μ_i)

z = ( z₁, z₂, ... z_m )^T

z_i = log( μ_i / N ) + ( y_i - μ_i ) / μ_i

X^TWz

=	\|	μ₁,	μ₂,	...	μ_m-1,	μ_m	\|\|	z₁	\|
	\|	μ₁,	0,	...	0,	0	\|\|	z₂	\|
	\|	0,	μ₂,	...	0,	0	\|\|	:	\|
	\|	:	:	...	:	:	\|\|	z_m-1	\|
	\|	0,	0,	...	μ_m-1,	0	\|\|	z_m	\|

=	\|	Σ_i{1→m}( μ_iz_i )	\|
	\|	μ₁z₁	\|
	\|	μ₂z₂	\|
	\|	:	\|
	\|	μ_m-1z_m-1	\|

なので、回帰係数の最尤推定量を a = ( a₀, a₁, ... a_m-1 )^T とすれば

	a₀Σ_i{1→m}( μ_i ) + Σ_i{1→m-1}( a_iμ_i )
=	a₀μ_m + Σ_i{1→m-1}( ( a₀ + a_i )μ_i )
=	Σ_i{1→m}( μ_iz_i )

a₀μ_i + a_iμ_i = ( a₀ + a_i )μ_i = μ_iz_i より a₀ + a_i = z_i

a₀μ_m + Σ_i{1→m-1}( μ_iz_i ) = Σ_i{1→m}( μ_iz_i ) より

a₀ = z_m

a₀ + a_i = z_i より

a_i = z_i - z_m

という結果が得られますが、一元配置分散分析での回帰係数が、n_i = 1 ならば

a₀ = y_m

a_i = y_i - y_m

多項分布を仮定した方式では、事象を複数としても成り立ちます。例えば、二つの事象 A, B の組み合わせに対して発生回数 y_j,k の同時確率が多項分布に従うと仮定することができます。さらに、それぞれの事象の発生回数に制約を設ける考え方もあります。事象 A がカテゴリ A_j ( 1 ≤ j ≤ p ) に分類できたとして、その発生回数の総計がそれぞれ n_j に固定されている場合、カテゴリ A_j の同時確率が多項分布

P_{n_j,θ}( y ) = n_j!Π_k{1→q}( θ_j,k^y_j,k / y_j,k! )

に従うと仮定します。q は事象 B の分類数を表し、Σ_k{1→q}( θ_j,k ) = 1、Σ_k{1→q}( y_j,k ) = n_j です。このとき、全事象に対する同時確率は

Π_j{1→p}( P_{n_j,θ}( y ) ) = Π_j{1→p}( n_j!Π_k{1→q}( θ_j,k^y_j,k / y_j,k! ) )

となります。これは「積多項分布 (Product Multinomial Distribution)」と呼ばれます。この場合の連結関数は、μ_j,k = n_jθ_j,k より

g(μ_j,k) = log( μ_j,k / n_j ) = log( θ_j,k ) = x_j,kα

二つの事象 A, B の組み合わせに対して多項分布を仮定したとき、事象 A, B の中であるカテゴリが発生する確率が互いに独立ならば、事象 A のカテゴリ A_j が発生する確率を θ_A,j、事象 B のカテゴリ B_k が発生する確率を θ_B,k としたとき、

θ_j,k = θ_A,jθ_B,k

μ_j,k = Nθ_A,j・Nθ_B,k / N ≡ μ_A,jμ_B,k / N

g(μ_j,k) = log( μ_A,jμ_B,k / N ) = log μ_A,j + log μ_B,k - log N

g(μ_j,k) = log( μ_j,k / N ) = log μ_j,k - log N

なので、この両者を比較することで独立性の推定・検定を行うことができます。これはよく利用される手法です。

X	=	\|	1_p,	D,	F₁,	D,	0_p,p-1,	...	0_p,p-1	\|
		\|	1_p,	D,	F₂,	0_p,p-1,	D,	...	0_p,p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_p,	D,	F_q-1,	0_p,p-1,	0_p,p-1,	...	D	\|
		\|	1_p,	D,	0_p,q-1,	0_p,p-1,	0_p,p-1,	...	0_p,p-1	\|

となります。ここで、D は行列数が p の単位行列 E_p から第 p 列を除外した p 行 p - 1 列の行列、F_k は k 列目の全要素が全て 1 で、残りは 0 の p 行 q - 1 列行列になります。回帰係数の最尤推定量を a = ( a₀, a₁, ... a_p-1, b₁, ... a_q-1, g_1,1, ... g_p-1,1, g_1,2, ... g_p-1,q-1 )^T とした時、j 行目のブロックの k 行目の要素にあたる式は

a₀ + a_j + b_k + g_j,k ( j < p, k < q )
a₀ + a_j ( j < p, k = q )
a₀ + b_k ( j = p, k < q )
a₀ ( j = p, k = q )

log( μ_p,q / N ) = a₀ より μ_p,q = Ne^a₀
log( μ_j,q / N ) = a₀ + a_j より μ_j,q = Ne^a₀e^a_j
log( μ_p,k / N ) = a₀ + b_k より μ_p,k = Ne^a₀e^b_k
log( μ_j,k / N ) = a₀ + a_j + b_k + g_j,k より μ_p,k = Ne^a₀e^a_je^b_ke^g_j,k

W は、μ_j,k を対角成分 ( 但し、μ_1,1, μ_2,1, ... μ_p,1, μ_1,2 ... の順に並べます ) とした行列数 pq の対角行列になるので、μ_1,k から μ_p,k までを対角成分とする行列数 p の対角行列を Μ_k とすれば、

X^TWX

1_p^T,

...

1_p^T,

1_p^T

Μ₁,

0_p,

...

0_p,

0_p

1_p,

F₁,

0_p,p-1,

...

0_p,p-1

D^T,

...

D^T,

D^T

0_p,

Μ₂,

...

0_p,

0_p

1_p,

F₂,

0_p,p-1,

...

0_p,p-1

F₁^T,

F₂^T,

...

F_q-1^T,

0_q-1,p

...

D^T,

0_p-1,p,

...

0_p-1,p,

0_p-1,p

0_p,

...

Μ_q-1,

0_p

1_p,

F_q-1,

0_p,p-1,

...

0_p-1,p,

D^T,

...

0_p-1,p,

0_p-1,p

0_p,

...

0_p,

Μ_q

1_p,

0_p,q-1,

0_p,p-1,

...

0_p,p-1

...

0_p-1,p,

...

D^T,

0_p-1,p

1_p^TΜ₁,

1_p^TΜ₂,

...

1_p^TΜ_q-1,

1_p^TΜ_q

1_p,

F₁,

0_p,p-1,

...

0_p,p-1

D^TΜ₁,

D^TΜ₂,

...

D^TΜ_q-1,

D^TΜ_q

1_p,

F₂,

0_p,p-1,

...

0_p,p-1

F₁^TΜ₁,

F₂^TΜ₂,

...

F_q-1^TΜ_q-1,

0_q-1,p

...

D^TΜ₁,

0_p-1,p,

...

0_p-1,p,

0_p-1,p

1_p,

F_q-1,

0_p,p-1,

...

0_p-1,p,

D^TΜ₂,

...

0_p-1,p,

0_p-1,p

1_p,

0_p,q-1,

0_p,p-1,

...

0_p,p-1

...

0_p-1,p,

...

D^TΜ_q-1,

0_p-1,p

=	\|	Σ_k{1→q}( 1_p^TΜ_k1_p ),	Σ_k{1→q}( 1_p^TΜ_kD ),	Σ_k{1→q-1}( 1_p^TΜ_kF_k ),	1_p^TΜ₁D,	1_p^TΜ₂D,	...	1_p^TΜ_q-1D	\|
	\|	Σ_k{1→q}( D^TΜ_k1_p ),	Σ_k{1→q}( D^TΜ_kD ),	Σ_k{1→q-1}( D^TΜ_kF_k ),	D^TΜ₁D,	D^TΜ₂D,	...	D^TΜ_q-1D	\|
	\|	Σ_k{1→q-1}( F_k^TΜ_k1_p ),	Σ_k{1→q-1}( F_k^TΜ_kD ),	Σ_k{1→q-1}( F_k^TΜ_kF_k ),	F₁^TΜ₁D,	F₂^TΜ₂D,	...	F_q-1^TΜ_q-1D	\|
	\|	D^TΜ₁1_p,	D^TΜ₁D,	D^TΜ₁F₁	D^TΜ₁D,	0_p-1,	...	0_p-1	\|
	\|	D^TΜ₂1_p,	D^TΜ₂D,	D^TΜ₂F₂	0_p-1,	D^TΜ₂D,	...	0_p-1	\|
	\|	:	:	:	:	:	...	:	\|
	\|	D^TΜ_q-11_p,	D^TΜ_q-1D,	D^TΜ_q-1F_q-1	0_p-1,	0_p-1,	...	D^TΜ_q-1D	\|

1_p^TΜ_k = ( μ_1,k, μ_2,k, ... μ_p,k ) より

1_p^TΜ_k1_p = Σ_j{1→p}( μ_j,k ) ≡ μ_B,k
Σ_k{1→q}( 1_p^TΜ_k1_p ) = Σ_k{1→q}( μ_B,k ) ≡ μ ( = N )

1_p^TΜ_kD = ( μ_1,k, μ_2,k, ... μ_p-1,k )

Σ_k{1→q}( 1_p^TΜ_kD )	=	( Σ_k{1→q}( μ_1,k ), Σ_k{1→q}( μ_2,k ), ... Σ_k{1→q}( μ_p-1,k ) )
	≡	( μ_A,1, μ_A,2, ... μ_A,p-1 )

1_p^TΜ_kF_k = ( 0, 0, ... μ_B,k, ... 0 ) ( k 列のみ要素を持つ行ベクトル )
Σ_k{1→q-1}( 1_p^TΜ_kF_k ) = ( μ_B,1, μ_B,2, ... μ_B,q-1 )

D^TΜ_k	=	\|	μ_1,k,	0,	...	0,	0	\|
		\|	0,	μ_2,k,	...	0,	0	\|
		\|	:	:	...	:	:	\|
		\|	0,	0,	...	μ_p-1,k,	0	\| より

D^TΜ_k1_p = ( μ_1,k, μ_2,k, ... μ_p-1,k )^T
Σ_k{1→q}( D^TΜ_k1_p ) = ( μ_A,1, μ_A,2., ... μ_A,p-1 )^T

D^TΜ_kD	=	\|	μ_1,k,	0,	...	0	\|
		\|	0,	μ_2,k,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	μ_p-1,k	\|

Σ_k{1→q}( D^TΜ_kD )	=	\|	μ_A,1,	0,	...	0	\|
		\|	0,	μ_A,2,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	μ_A,p-1	\|

D^TΜ_kF_k	=	\|	0,	0,	...	μ_1,k,	...	0	\|
		\|	0,	0,	...	μ_2,k,	...	0	\|
		\|	:	:	...	:	...	:	\|
		\|	0,	0,	...	μ_p-1,k,	...	0	\| ( k 列のみ要素を持つ行列 )

Σ_k{1→q-1}( D^TΜ_kF_k )	=	\|	μ_1,1,	μ_1,2,	...	μ_1,q-1	\|
		\|	μ_2,1,	μ_2,2,	...	μ_2,q-1	\|
		\|	:	:	...	:	\|
		\|	μ_p-1,1,	μ_p-1,2,	...	μ_p-1,q-1	\|

F_k^TΜ_k	=	\|	0,	0,	...	0	\|
		\|	0,	0,	...	0	\|
		\|	:	:	...	:	\|
		\|	μ_1,k,	μ_2,k,	...	μ_p,k	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	0	\| より

F_k^TΜ_k1_p = ( 0, 0, ... μ_B,k, ... 0 )^T ( k 行のみ要素を持つ列ベクトル )
Σ_k{1→q-1}( F_k^TΜ_k1_p ) = ( μ_B,1, μ_B,2, ... μ_B,q-1 )^T

F_k^TΜ_kD	=	\|	0,	0,	...	0	\|
		\|	0,	0,	...	0	\|
		\|	:	:	...	:	\|
		\|	μ_1,k,	μ_2,k,	...	μ_p-1,k	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	0	\|

Σ_k{1→q-1}( F_k^TΜ_kD )	=	\|	μ_1,1,	μ_2,1,	...	μ_p-1,1	\|
		\|	μ_1,2,	μ_2,2,	...	μ_p-1,2	\|
		\|	:	:	...	:	\|
		\|	μ_1,q-1,	μ_2,q-1,	...	μ_p-1,q-1	\|

F_k^TΜ_kF_k	=	\|	0,	...	0,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	μ_B,k,	...	0	\|
		\|	:	...	:	...	:	\|
		\|	0,	...	0,	...	0	\| ( k 行 k 列のみ要素を持つ行列 )

Σ_k{1→q-1}( F_k^TΜ_kF_k )	=	\|	μ_B,1,	0,	...	0	\|
		\|	0,	μ_B,2,	...	0	\|
		\|	:	:	...	:	\|
		\|	0,	0,	...	μ_B,q-1	\|

X^TWX =	\|	μ,	μ_A,1,	μ_A,2,	...	μ_A,p-1,	μ_B,1,	μ_B,2,	...	μ_B,q-1,	μ_1,1,	μ_2,1,	...	μ_p-1,1,	μ_1,2,	...	μ_p-1,q-1	\|

	\|	μ_A,1,	μ_A,1,	0,	...	0,	μ_1,1,	μ_1,2,	...	μ_1,q-1,	μ_1,1,	0,	...	0,	μ_1,2,	...	0	\|
	\|	μ_A,2,	0,	μ_A,2,	...	0,	μ_2,1,	μ_2,2,	...	μ_2,q-1,	0,	μ_2,1,	...	0,	0,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	...	:	\|
	\|	μ_A,p-1,	0,	0,	...	μ_A,p-1,	μ_p-1,1,	μ_p-1,2,	...	μ_p-1,q-1,	0,	0,	...	μ_p-1,1,	0,	...	μ_p-1,q-1	\|

	\|	μ_B,1,	μ_1,1,	μ_2,1,	...	μ_p-1,1,	μ_B,1,	0,	...	0,	μ_1,1,	μ_2,1,	...	μ_p-1,1,	0,	...	0	\|
	\|	μ_B,2,	μ_1,2,	μ_2,2,	...	μ_p-1,2,	0,	μ_B,2,	...	0,	0,	0,	...	0,	μ_1,2,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	...	:	\|
	\|	μ_B,q-1,	μ_1,q-1,	μ_2,q-1,	...	μ_p-1,q-1,	0,	0,	...	μ_B,q-1,	0,	0,	...	0,	0,	...	μ_p-1,q-1	\|

	\|	μ_1,1,	μ_1,1,	0,	...	0,	μ_1,1,	0,	...	0,	μ_1,1,	0,	...	0,	0,	...	0	\|
	\|	μ_2,1,	0,	μ_2,1,	...	0,	μ_2,1,	0,	...	0,	0,	μ_2,1,	...	0,	0,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	...	:	\|
	\|	μ_p-1,1,	0,	0,	...	μ_p-1,1,	μ_p-1,1,	0,	...	0,	0,	0,	...	μ_p-1,1,	0,	...	0	\|
	\|	μ_1,2,	μ_1,2,	0,	...	0,	0,	μ_1,2,	...	0,	0,	0,	...	0,	μ_1,2,	...	0	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	...	:	\|
	\|	μ_p-1,q-1,	0,	0,	...	μ_p-1,q-1,	0,	0,	...	μ_p-1,q-1,	0,	0,	...	0,	0,	...	μ_p-1,q-1	\|

という結果が得られます。μ_j,k = 1 とすれば、これは ( n = 1 のときの ) 二元配置分散分析での左辺の係数行列 X^TX と等しくなります。各行の式は、1 行目、2 行目からの p - 1 行分、その次の q - 1 行分、残りの行と分けた時、

μa₀ + Σ_j{1→p-1}( μ_A,ja_j ) + Σ_k{1→q-1}( μ_B,kb_k ) + Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kg_j,k ) )
μ_A,ja₀ + μ_A,ja_j + Σ_k{1→q-1}( μ_j,kb_k ) + Σ_k{1→q-1}( μ_j,kg_j,k ) [ 1 ≤ j ≤ p - 1 ]
μ_B,ka₀ + Σ_j{1→p-1}( μ_j,ka_j ) + μ_B,kb_k + Σ_j{1→p-1}( μ_j,kg_j,k ) [ 1 ≤ k ≤ q - 1 ]
μ_j,ka₀ + μ_j,ka_j + μ_j,kb_k + μ_j,kg_j,k [ 1 ≤ j ≤ p - 1, 1 ≤ k ≤ q - 1 ]

z_j,k	=	g(μ_j,k) + ( y_j,k - μ_j,k )g'(μ_j,k)
	=	log( μ_j,k / N ) + ( y_j,k - μ_j,k ) / μ_j,k

z_k = ( z_1,k, z_2,k, ... z_p,k )^T

z = ( z₁, z₂, ... z_q )^T

X^TWz	=	\|	1_p^TΜ₁,	1_p^TΜ₂,	...	1_p^TΜ_q-1,	1_p^TΜ_q	\|\|	z₁	\|
		\|	D^TΜ₁,	D^TΜ₂,	...	D^TΜ_q-1,	D^TΜ_q	\|\|	z₂	\|
		\|	F₁^TΜ₁,	F₂^TΜ₂,	...	F_q-1^TΜ_q-1,	0_q-1,p	\|\|	:	\|
		\|	D^TΜ₁,	0_p-1,p,	...	0_p-1,p,	0_p-1,p	\|\|	z_q-1	\|
		\|	0_p-1,p,	D^TΜ₂,	...	0_p-1,p,	0_p-1,p	\|\|	z_q	\|
		\|	:	:	...	:	:	\|
		\|	0_p-1,p,	0_p-1,p,	...	D^TΜ_q-1,	0_p-1,p	\|

Σ_k{1→q}( 1_p^TΜ_kz_k )	=	Σ_k{1→q}( Σ_j{1→p}( μ_j,kz_j,k ) )
	≡	S_z
Σ_k{1→q}( D^TΜ_kz_k )	=	( Σ_k{1→q}( μ_1,kz_1,k ), Σ_k{1→q}( μ_2,kz_2,k ), ... Σ_k{1→q}( μ_p-1,kz_p-1,k ) )^T
	≡	( S_A,1, S_A,2, ... S_A,p-1 ) ≡ S_A
Σ_k{1→q-1}( F_k^TΜ_kz_k )	=	( Σ_j{1→p}( μ_j,1z_j,1 ), Σ_j{1→p}( μ_j,2z_j,2 ), ... Σ_j{1→p}( μ_j,q-1z_j,q-1 ) )^T
	≡	( S_B,1, S_B,2, ... S_B,q-1 ) ≡ S_B
D^TΜ_kz_k	=	( μ_1,kz_1,k, μ_2,kz_2,k, ... μ_p,kz_p-1,k )^T [ 1 ≤ k ≤ q - 1 ]

という結果が得られます。μ_j,k = 1 かつ z_j,k = y_j,k ならば、これも二元配置分散分析と全く同じ式になります。それぞれの正規方程式は

μa₀ + Σ_j{1→p-1}( μ_A,ja_j ) + Σ_k{1→q-1}( μ_B,kb_k ) + Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kg_j,k ) ) = S_z --- (1)
μ_A,ja₀ + μ_A,ja_j + Σ_k{1→q-1}( μ_j,kb_k ) + Σ_k{1→q-1}( μ_j,kg_j,k ) = S_A,j [ 1 ≤ j ≤ p - 1 ] --- (2)
μ_B,ka₀ + Σ_j{1→p-1}( μ_j,ka_j ) + μ_B,kb_k + Σ_j{1→p-1}( μ_j,kg_j,k ) = S_B,k [ 1 ≤ k ≤ q - 1 ] --- (3)
μ_j,ka₀ + μ_j,ka_j + μ_j,kb_k + μ_j,kg_j,k = μ_j,kz_j,k [ 1 ≤ j ≤ p - 1, 1 ≤ k ≤ q - 1 ] --- (4)

	μ_A,ja₀ + μ_A,ja_j + Σ_k{1→q-1}( μ_j,kb_k ) + Σ_k{1→q-1}( μ_j,kg_j,k )
=	Σ_k{1→q-1}( μ_j,k( a₀ + a_j ) ) + μ_j,q( a₀ + a_j ) + Σ_k{1→q-1}( μ_j,kb_k ) + Σ_k{1→q-1}( μ_j,kg_j,k )
=	Σ_k{1→q-1}( μ_j,k( a₀ + a_j + b_k + g_j,k ) ) + μ_j,q( a₀ + a_j )
=	Σ_k{1→q-1}( μ_j,kz_j,k ) + μ_j,q( a₀ + a_j )
	μ_B,ka₀ + Σ_j{1→p-1}( μ_j,ka_j ) + μ_B,kb_k + Σ_j{1→p-1}( μ_j,kg_j,k )
=	Σ_j{1→p-1}( μ_j,kz_j,k ) + μ_p,k( a₀ + b_k )

	Σ_j{1→p-1}( Σ_k{1→q-1}( μ_j,kz_j,k ) + μ_j,q( a₀ + a_j ) )
=	Σ_j{1→p-1}( Σ_k{1→q-1}( μ_j,kz_j,k ) ) + ( μ_B,q - μ_p,q )a₀ + Σ_j{1→p-1}( μ_j,qa_j )

	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) + μ_p,k( a₀ + b_k ) )
=	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + ( μ_A,p - μ_p,q )a₀ + Σ_k{1→q-1}( μ_p,kb_k )

2Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + ( μ_A,p + μ_B,q - 2μ_p,q )a₀

+ Σ_j{1→p-1}( μ_j,qa_j ) + Σ_k{1→q-1}( μ_p,kb_k )

μa₀	=	Σ_j{1→p}( Σ_k{1→q}( μ_j,ka₀ ) )
	=	Σ_j{1→p}( Σ_k{1→q-1}( μ_j,ka₀ ) ) + μ_B,qa₀
	=	Σ_j{1→p-1}( Σ_k{1→q-1}( μ_j,ka₀ ) ) + ( μ_A,p + μ_B,q - μ_p,q )a₀

Σ_j{1→p-1}( μ_A,ja_j )	=	Σ_j{1→p-1}( Σ_k{1→q}( μ_j,ka_j ) )
	=	Σ_j{1→p-1}( Σ_k{1→q-1}( μ_j,ka_j ) + μ_j,qa_j )
	=	Σ_j{1→p-1}( Σ_k{1→q-1}( μ_j,ka_j ) ) + Σ_j{1→p-1}( μ_j,qa_j )
Σ_k{1→q-1}( μ_B,kb_k )	=	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kb_k ) ) + Σ_k{1→q-1}( μ_p,kb_k )

	μa₀ + Σ_j{1→p-1}( μ_A,ja_j ) + Σ_k{1→q-1}( μ_B,kb_k ) + Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kg_j,k ) )
=	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,ka₀ + μ_j,ka_j + μ_j,kb_k + μ_j,kg_j,k ) )
	+ ( μ_A,p + μ_B,q - μ_p,q )a₀ + Σ_j{1→p-1}( μ_j,qa_j ) + Σ_k{1→q-1}( μ_p,kb_k )
=	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + ( μ_A,p + μ_B,q - μ_p,q )a₀
	+ Σ_j{1→p-1}( μ_j,qa_j ) + Σ_k{1→q-1}( μ_p,kb_k )

	Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + ( μ_A,p + μ_B,q - μ_p,q )a₀
	+ Σ_j{1→p-1}( μ_j,qa_j ) + Σ_k{1→q-1}( μ_p,kb_k )
	- 2Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) - ( μ_A,p + μ_B,q - 2μ_p,q )a₀
	- Σ_j{1→p-1}( μ_j,qa_j ) - Σ_k{1→q-1}( μ_p,kb_k )
=	-Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + μ_p,qa₀

	S_z - Σ_j{1→p-1}( S_A,j ) - Σ_k{1→q-1}( S_B,k )
=	Σ_k{1→q}( Σ_j{1→p}( μ_j,kz_j,k ) ) - Σ_j{1→p-1}( Σ_k{1→q}( μ_j,kz_j,k ) ) - Σ_k{1→q-1}( Σ_j{1→p}( μ_j,kz_j,k ) )
=	Σ_k{1→q}( μ_p,kz_p,k ) - Σ_j{1→p}( Σ_k{1→q-1}( μ_j,kz_j,k ) )
=	-Σ_k{1→q-1}( Σ_j{1→p}( μ_j,kz_j,k ) - μ_p,kz_p,k ) + μ_p,qz_p,q
=	-Σ_k{1→q-1}( Σ_j{1→p-1}( μ_j,kz_j,k ) ) + μ_p,qz_p,q

となり、これは二元配置分散分析において a₀ = m_p,q となることに対応します。これを (2), (3) に代入すれば

Σ_k{1→q-1}( μ_j,kz_j,k ) + μ_j,q( z_p,q + a_j ) = Σ_k{1→q}( μ_j,kz_j,k )
Σ_j{1→p-1}( μ_j,kz_j,k ) + μ_p,k( z_p,q + b_k ) = Σ_j{1→p}( μ_j,kz_j,k )

a_j = z_j,q - z_p,q
b_k = z_p,k - z_p,q

で、二元配置分散分析において a_j = m_j,q - m_p,q、b_k = m_p,k - m_p,q となることに対応し、(4) より

g_j,k	=	z_j,k - a₀ - a_j - b_k
	=	z_j,k - z_p,q - ( z_j,q - z_p,q ) - ( z_p,k - z_p,q )
	=	( z_j,k - z_p,k ) - ( z_j,q - z_p,q )

となるので、二元配置分散分析において g_j,k = ( m_j,k - m_p,k ) - ( m_j,q - m_p,q ) となることに対応します。これらの最尤推定量を使って求められる対数尤度統計量は、二元配置分散分析における誤差変動 D_E に対応します。このように、計数データに多項分布や積多項分布と対数関数を組み合わせたモデルは分散分析と非常によく似た性質を持っていますが、右辺には μ_j,k が存在してこの中にも係数 α が (しかも線形式ではない状態で) 含まれているので、分散分析のときのように連立方程式を解くことはできません。従って、スコア法を利用して係数を求める必要があります。これは、交互作用効果 g_j,k を含む独立変数を除外した縮小モデル (加法モデル) でも同様です。

log θ_j,k = α₀ ( 定数 )

となるので θ_j,k が全て等しいことになります。しかし、θ_j,k の総和は 1 なので、θ_j,k = 1 / pq ( α₀ = -log pq ) であり、従って μ_j,k = N / pq になります。

ポアソン回帰を含め、今まで紹介したモデルはすべて連結関数に対数関数を用いてきました。計数データの場合、線形式の値に関係なく μ は正の数となるので、対数関数を用いることは非常に都合がいいことになります。特に、指数型分布族にポアソン分布を仮定すれば回帰方程式も比較的シンプルな形にすることができます。このような形式の一般化線形モデルは「対数線形モデル (Log-linear Model)」と呼ばれます。

前節で示したサンプル・プログラム PoissonRegression を使って、参考文献にあったデータを処理してみたいと思います。このデータは、悪性黒色腫と呼ばれる皮膚がんの患者数を、三つの部位 ( 頭頸部・体幹部・四肢 ) と四つの病型 ( ハッチンソン黒色斑・表在拡大型黒色腫・結節型黒色腫・不明 ) の組み合わせごとに分類した結果です。患者数は 400 名と固定されているので、各数は N = 400 とした多項分布に従うと解釈することができます。

表 3-1. 悪性黒色腫の病型・部位別計数
病型	部位
病型	頭頸部	体幹部	四肢	計
ハッチンソン黒色斑	22	2	10	34
表在拡大型黒色腫	16	54	115	185
結節型黒色腫	19	33	73	125
不明	11	17	28	56
計	68	106	226	400

まずは、このデータに飽和モデルを適用して係数の最尤推定量を計算してみます ( verbose = true にして冗長モードで出力しています )。

*** Poisson Regression ***

ni = ( 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400 )

*** Scoring Method ***

Exponential Family of Distribution : Poisson Distribution
Link Function : Logarithm Function ( for Poisson Regression )

N = 12; p = 12

x = ( 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0 )
    ( 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0 )
    ( 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 )
    ( 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 )
    ( 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0 )
    ( 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 )
    ( 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0 )
    ( 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1 )
    ( 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 )
    ( 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 )
    ( 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 )
    ( 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 )

y = ( 22, 2, 10, 16, 54, 115, 19, 33, 73, 11, 17, 28 )

----- cnt = 0 -----

mu = ( 22, 2, 10, 16, 54, 115, 19, 33, 73, 11, 17, 28 )
g'(mu) = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )
w = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )

Equation System :
(400)x0 + (68)x1 + (106)x2 + (34)x3 + (185)x4 + (125)x5 + (22)x6 + (2)x7 + (16)x8 + (54)x9 + (19)x10 + (33)x11 = -846.361
(68)x0 + (68)x1 + (0)x2 + (22)x3 + (16)x4 + (19)x5 + (22)x6 + (0)x7 + (16)x8 + (0)x9 + (19)x10 + (0)x11 = -212.734
(106)x0 + (0)x1 + (106)x2 + (2)x3 + (54)x4 + (33)x5 + (0)x6 + (2)x7 + (0)x8 + (54)x9 + (0)x10 + (33)x11 = -254.754
(34)x0 + (22)x1 + (2)x2 + (34)x3 + (0)x4 + (0)x5 + (22)x6 + (2)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -111.295
(185)x0 + (16)x1 + (54)x2 + (0)x3 + (185)x4 + (0)x5 + (0)x6 + (0)x7 + (16)x8 + (54)x9 + (0)x10 + (0)x11 = -302.987
(125)x0 + (19)x1 + (33)x2 + (0)x3 + (0)x4 + (125)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (19)x10 + (33)x11 = -264.4
(22)x0 + (22)x1 + (0)x2 + (22)x3 + (0)x4 + (0)x5 + (22)x6 + (0)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -63.8093
(2)x0 + (0)x1 + (2)x2 + (2)x3 + (0)x4 + (0)x5 + (0)x6 + (2)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -10.5966
(16)x0 + (16)x1 + (0)x2 + (0)x3 + (16)x4 + (0)x5 + (0)x6 + (0)x7 + (16)x8 + (0)x9 + (0)x10 + (0)x11 = -51.502
(54)x0 + (0)x1 + (54)x2 + (0)x3 + (54)x4 + (0)x5 + (0)x6 + (0)x7 + (0)x8 + (54)x9 + (0)x10 + (0)x11 = -108.134
(19)x0 + (19)x1 + (0)x2 + (0)x3 + (0)x4 + (19)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (19)x10 + (0)x11 = -57.8935
(33)x0 + (0)x1 + (33)x2 + (0)x3 + (0)x4 + (33)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (0)x10 + (33)x11 = -82.3336

Regression equation : y = -2.65926x0 + -0.934309x1 + -0.498991x2 + -1.02962x3 + 1.41273x4 + 0.958255x5 + 1.72277x6 + -1.11045x7 + -1.03803x8 + -0.256957x9 + -0.411711x10 + -0.294961x11

----- cnt = 1 -----

mu = ( 22, 2, 10, 16, 54, 115, 19, 33, 73, 11, 17, 28 )
g'(mu) = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )
w = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )

Equation System :
(400)x0 + (68)x1 + (106)x2 + (34)x3 + (185)x4 + (125)x5 + (22)x6 + (2)x7 + (16)x8 + (54)x9 + (19)x10 + (33)x11 = -846.361
(68)x0 + (68)x1 + (0)x2 + (22)x3 + (16)x4 + (19)x5 + (22)x6 + (0)x7 + (16)x8 + (0)x9 + (19)x10 + (0)x11 = -212.734
(106)x0 + (0)x1 + (106)x2 + (2)x3 + (54)x4 + (33)x5 + (0)x6 + (2)x7 + (0)x8 + (54)x9 + (0)x10 + (33)x11 = -254.754
(34)x0 + (22)x1 + (2)x2 + (34)x3 + (0)x4 + (0)x5 + (22)x6 + (2)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -111.295
(185)x0 + (16)x1 + (54)x2 + (0)x3 + (185)x4 + (0)x5 + (0)x6 + (0)x7 + (16)x8 + (54)x9 + (0)x10 + (0)x11 = -302.987
(125)x0 + (19)x1 + (33)x2 + (0)x3 + (0)x4 + (125)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (19)x10 + (33)x11 = -264.4
(22)x0 + (22)x1 + (0)x2 + (22)x3 + (0)x4 + (0)x5 + (22)x6 + (0)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -63.8093
(2)x0 + (0)x1 + (2)x2 + (2)x3 + (0)x4 + (0)x5 + (0)x6 + (2)x7 + (0)x8 + (0)x9 + (0)x10 + (0)x11 = -10.5966
(16)x0 + (16)x1 + (0)x2 + (0)x3 + (16)x4 + (0)x5 + (0)x6 + (0)x7 + (16)x8 + (0)x9 + (0)x10 + (0)x11 = -51.502
(54)x0 + (0)x1 + (54)x2 + (0)x3 + (54)x4 + (0)x5 + (0)x6 + (0)x7 + (0)x8 + (54)x9 + (0)x10 + (0)x11 = -108.134
(19)x0 + (19)x1 + (0)x2 + (0)x3 + (0)x4 + (19)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (19)x10 + (0)x11 = -57.8935
(33)x0 + (0)x1 + (33)x2 + (0)x3 + (0)x4 + (33)x5 + (0)x6 + (0)x7 + (0)x8 + (0)x9 + (0)x10 + (33)x11 = -82.3336

Regression equation : y = -2.65926x0 + -0.934309x1 + -0.498991x2 + -1.02962x3 + 1.41273x4 + 0.958255x5 + 1.72277x6 + -1.11045x7 + -1.03803x8 + -0.256957x9 + -0.411711x10 + -0.294961x11

Estimated regression equation : y = -2.65926x0 + -0.934309x1 + -0.498991x2 + -1.02962x3 + 1.41273x4 + 0.958255x5 + 1.72277x6 + -1.11045x7 + -1.03803x8 + -0.256957x9 + -0.411711x10 + -0.294961x11
variance of a = ( 0.0357143, 0.126623, 0.0945378, 0.135714, 0.0444099, 0.0494129, 0.272078, 0.694538, 0.197819, 0.121752, 0.192954, 0.138539 )

このデータは列 j が 3 個、行 k が 4 個で、基準となるカテゴリは ( j, k ) = ( 3, 4 ) としています。デザイン行列は 3 x 4 = 12 行、1 + ( 3 - 1 ) + ( 4 - 1 ) + ( 3 - 1 ) x ( 4 - 1 ) = 12 列で、実際の要素はサンプル・プログラムの出力結果で確認することができます。指数型分布族として使用するポアソン分布のパラメータ ni は全て 400 として計算しています。各独立変数 x0 から x11 はそれぞれ

x0 ... 定数項
x1, x2 ... 部位に対する差異を表すダミー変数 ( j = 1, 2 )
x3, x4, x5 ... 病型に対する差異を表すダミー変数 ( k = 1, 2, 3 )
x6 ～ x11 ... 病型に対する差異を表すダミー変数 ( ( j, k ) = ( 1, 1 ), ( 2, 1 ), ( 1, 2 ), ( 2, 2 ), ( 1, 3 ), ( 2, 3 ) )

を表します。回帰式を使って ( j, k ) = ( 3, 4 ) すなわち基準カテゴリ ( 部位 = 四肢 ; 病型 = 不明 ) の期待値を求めると

log θ_3,4 = -2.65926 より θ_3,4 = e^-2.65926 = 0.0700

μ_3,4 = 400 x 0.0700 = 28

となって観測値と一致します。( j, k ) = ( 1, 1 ) ( 部位 = 頭頸部 ; 病型 = ハッチンソン黒色斑 ) では

log θ_3,4 = -2.65926 - 0.934309 - 1.02962 + 1.72277 = -2.90 より θ_1,1 = e^-2.90 = 0.0550

μ_1,1 = 400 x 0.0550 = 22

であり、やはり観測値と一致します。このように、飽和モデルでは期待値が観測値と完全に一致します。飽和モデルの回帰係数は

a₀ = z_p,q
a_j = z_j,q - z_p,q
b_k = z_p,k - z_p,q
g_j,k = ( z_j,k - z_p,k ) - ( z_j,q - z_p,q )

log( μ_j,k / N ) = a₀ + a_j + b_k + g_j,k = z_j,k [ 1 ≤ j ≤ p - 1 ; 1 ≤ k ≤ q - 1 ]
log( μ_p,k / N ) = a₀ + b_k = z_p,k [ 1 ≤ k ≤ q - 1 ]
log( μ_j,q / N ) = a₀ + a_j = z_j,q [ 1 ≤ j ≤ p - 1 ]
log( μ_p,q / N ) = a₀ = z_p,q

z_j,k = log( μ_j,k / N ) + ( y_j,k - μ_j,k ) / μ_j,k

log( μ_j,k / N ) = log( μ_j,k / N ) + ( y_j,k - μ_j,k ) / μ_j,k

なので、μ_j,k = y_j,k が成り立ち、スコア法の結果と一致します。

次に、加法モデルを適用して係数の最尤推定量を計算すると以下の結果が得られます ( verbose = true にして冗長モードで出力しています )。

*** Poisson Regression ***

ni = ( 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400 )

*** Scoring Method ***

Exponential Family of Distribution : Poisson Distribution
Link Function : Logarithm Function ( for Poisson Regression )

N = 12; p = 6

x = ( 1, 1, 0, 1, 0, 0 )
    ( 1, 0, 1, 1, 0, 0 )
    ( 1, 0, 0, 1, 0, 0 )
    ( 1, 1, 0, 0, 1, 0 )
    ( 1, 0, 1, 0, 1, 0 )
    ( 1, 0, 0, 0, 1, 0 )
    ( 1, 1, 0, 0, 0, 1 )
    ( 1, 0, 1, 0, 0, 1 )
    ( 1, 0, 0, 0, 0, 1 )
    ( 1, 1, 0, 0, 0, 0 )
    ( 1, 0, 1, 0, 0, 0 )
    ( 1, 0, 0, 0, 0, 0 )

y = ( 22, 2, 10, 16, 54, 115, 19, 33, 73, 11, 17, 28 )

----- cnt = 0 -----

mu = ( 22, 2, 10, 16, 54, 115, 19, 33, 73, 11, 17, 28 )
g'(mu) = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )
w = ( 0.0454545, 0.5, 0.1, 0.0625, 0.0185185, 0.00869565, 0.0526316, 0.030303, 0.0136986, 0.0909091, 0.0588235, 0.0357143 )

Equation System :
(400)x0 + (68)x1 + (106)x2 + (34)x3 + (185)x4 + (125)x5 = -846.361
(68)x0 + (68)x1 + (0)x2 + (22)x3 + (16)x4 + (19)x5 = -212.734
(106)x0 + (0)x1 + (106)x2 + (2)x3 + (54)x4 + (33)x5 = -254.754
(34)x0 + (22)x1 + (2)x2 + (34)x3 + (0)x4 + (0)x5 = -111.295
(185)x0 + (16)x1 + (54)x2 + (0)x3 + (185)x4 + (0)x5 = -302.987
(125)x0 + (19)x1 + (33)x2 + (0)x3 + (0)x4 + (125)x5 = -264.4

Regression equation : y = -2.56211x0 + -1.09147x1 + -0.717314x2 + 0.0371766x3 + 1.22812x4 + 0.802184x5

----- cnt = 1 -----

mu = ( 10.7517, 15.6303, 32.0253, 35.375, 51.4267, 105.369, 23.1055, 33.5897, 68.8228, 10.3593, 15.0599, 30.8566 )
g'(mu) = ( 0.0930087, 0.0639782, 0.0312253, 0.0282685, 0.0194452, 0.00949043, 0.0432798, 0.029771, 0.0145301, 0.0965315, 0.0664014, 0.032408 )
w = ( 0.0930087, 0.0639782, 0.0312253, 0.0282685, 0.0194452, 0.00949043, 0.0432798, 0.029771, 0.0145301, 0.0965315, 0.0664014, 0.032408 )

Equation System :
(432.372)x0 + (79.5915)x1 + (115.707)x2 + (58.4074)x3 + (192.171)x4 + (125.518)x5 = -971.159
(79.5915)x0 + (79.5915)x1 + (0)x2 + (10.7517)x3 + (35.375)x4 + (23.1055)x5 = -240.006
(115.707)x0 + (0)x1 + (115.707)x2 + (15.6303)x3 + (51.4267)x4 + (33.5897)x5 = -298.474
(58.4074)x0 + (10.7517)x1 + (15.6303)x2 + (58.4074)x3 + (0)x4 + (0)x5 = -194.829
(192.171)x0 + (35.375)x1 + (51.4267)x2 + (0)x3 + (192.171)x4 + (0)x5 = -339.025
(125.518)x0 + (23.1055)x1 + (33.5897)x2 + (0)x3 + (0)x4 + (125.518)x5 = -270.734

Regression equation : y = -2.53886x0 + -1.19039x1 + -0.754493x2 + -0.375803x3 + 1.19571x4 + 0.802959x5

----- cnt = 2 -----

mu = ( 6.59568, 10.1992, 21.689, 31.751, 49.0982, 104.409, 21.4382, 33.1511, 70.4968, 9.60437, 14.8517, 31.5827 )
g'(mu) = ( 0.151614, 0.0980466, 0.0461063, 0.0314951, 0.0203673, 0.00957772, 0.0466456, 0.030165, 0.014185, 0.104119, 0.0673322, 0.0316629 )
w = ( 0.151614, 0.0980466, 0.0461063, 0.0314951, 0.0203673, 0.00957772, 0.0466456, 0.030165, 0.014185, 0.104119, 0.0673322, 0.0316629 )

Equation System :
(404.867)x0 + (69.3893)x1 + (107.3)x2 + (38.4839)x3 + (185.258)x4 + (125.086)x5 = -888.832
(69.3893)x0 + (69.3893)x1 + (0)x2 + (6.59568)x3 + (31.751)x4 + (21.4382)x5 = -207.459
(107.3)x0 + (0)x1 + (107.3)x2 + (10.1992)x3 + (49.0982)x4 + (33.1511)x5 = -273.184
(38.4839)x0 + (6.59568)x1 + (10.1992)x2 + (38.4839)x3 + (0)x4 + (0)x5 = -132.198
(185.258)x0 + (31.751)x1 + (49.0982)x2 + (0)x3 + (185.258)x4 + (0)x5 = -323.928
(125.086)x0 + (21.4382)x1 + (33.1511)x2 + (0)x3 + (0)x4 + (125.086)x5 = -267.755

Regression equation : y = -2.53707x0 + -1.20087x1 + -0.757068x2 + -0.491626x3 + 1.195x4 + 0.802962x5

----- cnt = 3 -----

mu = ( 5.82349, 9.07662, 19.3515, 31.4541, 49.0251, 104.522, 21.2528, 33.125, 70.6231, 9.52124, 14.84, 31.6392 )
g'(mu) = ( 0.171718, 0.110173, 0.0516756, 0.0317923, 0.0203977, 0.00956734, 0.0470527, 0.0301886, 0.0141597, 0.105028, 0.0673853, 0.0316064 )
w = ( 0.171718, 0.110173, 0.0516756, 0.0317923, 0.0203977, 0.00956734, 0.0470527, 0.0301886, 0.0141597, 0.105028, 0.0673853, 0.0316064 )

Equation System :
(400.254)x0 + (68.0516)x1 + (106.067)x2 + (34.2516)x3 + (185.002)x4 + (125.001)x5 = -873.139
(68.0516)x0 + (68.0516)x1 + (0)x2 + (5.82349)x3 + (31.4541)x4 + (21.2528)x5 = -202.634
(106.067)x0 + (0)x1 + (106.067)x2 + (9.07662)x3 + (49.0251)x4 + (33.125)x5 = -268.744
(34.2516)x0 + (5.82349)x1 + (9.07662)x2 + (34.2516)x3 + (0)x4 + (0)x5 = -117.854
(185.002)x0 + (31.4541)x1 + (49.0251)x2 + (0)x3 + (185.002)x4 + (0)x5 = -323.173
(125.001)x0 + (21.2528)x1 + (33.125)x2 + (0)x3 + (0)x4 + (125.001)x5 = -267.366

Regression equation : y = -2.53704x0 + -1.20103x1 + -0.757096x2 + -0.498964x3 + 1.195x4 + 0.802962x5

----- cnt = 4 -----

mu = ( 5.78016, 9.01024, 19.2105, 31.45, 49.025, 104.525, 21.25, 33.125, 70.625, 9.52, 14.84, 31.64 )
g'(mu) = ( 0.173006, 0.110985, 0.0520548, 0.0317965, 0.0203978, 0.00956709, 0.0470588, 0.0301887, 0.0141593, 0.105042, 0.0673854, 0.0316056 )
w = ( 0.173006, 0.110985, 0.0520548, 0.0317965, 0.0203978, 0.00956709, 0.0470588, 0.0301887, 0.0141593, 0.105042, 0.0673854, 0.0316056 )

Equation System :
(400.001)x0 + (68.0002)x1 + (106)x2 + (34.0009)x3 + (185)x4 + (125)x5 = -872.262
(68.0002)x0 + (68.0002)x1 + (0)x2 + (5.78016)x3 + (31.45)x4 + (21.25)x5 = -202.428
(106)x0 + (0)x1 + (106)x2 + (9.01024)x3 + (49.025)x4 + (33.125)x5 = -268.492
(34.0009)x0 + (5.78016)x1 + (9.01024)x2 + (34.0009)x3 + (0)x4 + (0)x5 = -116.992
(185)x0 + (31.45)x1 + (49.025)x2 + (0)x3 + (185)x4 + (0)x5 = -323.166
(125)x0 + (21.25)x1 + (33.125)x2 + (0)x3 + (0)x4 + (125)x5 = -267.361

Regression equation : y = -2.53704x0 + -1.20103x1 + -0.757096x2 + -0.498991x3 + 1.195x4 + 0.802962x5

Estimated regression equation : y = -2.53704x0 + -1.20103x1 + -0.757096x2 + -0.498991x3 + 1.195x4 + 0.802962x5
variance of a = ( 0.0197819, 0.0191306, 0.0138587, 0.0472681, 0.0232625, 0.0258571 )

加法モデルでは交互作用効果が除外されるので、デザイン行列の列数は 1 + ( 3 - 1 ) + ( 4 - 1 ) = 6 列になります。回帰式を使って ( j, k ) = ( 3, 4 ) すなわち基準カテゴリ ( 部位 = 四肢 ; 病型 = 不明 ) の期待値を求めると

log θ_3,4 = -2.53704 より θ_3,4 = e^-2.53704 = 0.0791

μ_3,4 = 400 x 0.0791 = 31.6

という結果が得られます。他の組み合わせについても期待値を求めると以下のようになります。

表 3-2. 悪性黒色腫の病型・部位別計数 (期待値)
病型	部位
病型	頭頸部	体幹部	四肢	計
ハッチンソン黒色斑	5.78	9.01	19.21	34
表在拡大型黒色腫	31.45	49.02	104.53	185
結節型黒色腫	21.25	33.12	70.62	125
不明	9.52	14.84	31.64	56
計	68	106	226	400

部位 = 頭頸部 ; 病型 = ハッチンソン黒色斑を例にとると、回帰式は

y_1,1 = -2.53704 - 1.20103 - 0.498991 = -4.237061 より
μ_1,1 = 400 x e^-4.237061 = 5.78

この期待値は、各行・列の和の積を総和で割った結果と一致します。例えば ( j, k ) = ( 1, 1 ) ( 部位 = 頭頸部 ; 病型 = ハッチンソン黒色斑 ) の場合

34 x 68 / 400 = 5.78

μa₀ + Σ_j{1→p-1}( μ_A,ja_j ) + Σ_k{1→q-1}( μ_B,kb_k ) = S_z --- (4)
μ_A,ja₀ + μ_A,ja_j + Σ_k{1→q-1}( μ_j,kb_k ) = S_A,j [ 1 ≤ j ≤ p - 1 ] --- (5)
μ_B,ka₀ + Σ_j{1→p-1}( μ_j,ka_j ) + μ_B,kb_k = S_B,k [ 1 ≤ k ≤ q - 1 ] --- (6)

となりますが、右辺の S_A,j, S_B,k, S_z は ( a_p = b_q = 0 とみなして )

S_A,j	=	Σ_k{1→q}( μ_j,kz_j,k )
	=	Σ_k{1→q}( μ_j,klog( μ_j,k / N ) + ( y_j,k - μ_j,k ) )
	=	Σ_k{1→q}( μ_j,k( a₀ + a_j + b_k ) + ( y_j,k - μ_j,k ) )
	=	μ_A,ja₀ + μ_A,ja_j + Σ_k{1→q-1}( μ_j,kb_k ) + Σ_k{1→q}( y_j,k - μ_j,k )
S_B,k	=	μ_B,ka₀ + Σ_j{1→p-1}( μ_j,ka_j ) + μ_B,kb_k + Σ_j{1→p}( y_j,k - μ_j,k )
S_z	=	Σ_k{1→q}( Σ_j{1→p}( μ_j,kz_j,k ) )
	=	Σ_k{1→q}( Σ_j{1→p}( μ_j,klog( μ_j,k / N ) + ( y_j,k - μ_j,k ) ) )
	=	Σ_k{1→q}( Σ_j{1→p}( μ_j,k( a₀ + a_j + b_k ) + ( y_j,k - μ_j,k ) ) )
	=	μa₀ + Σ_j{1→p-1}( μ_A,ja_j ) + Σ_k{1→q-1}( μ_B,kb_k ) + Σ_k{1→q}( Σ_j{1→p}( y_j,k - μ_j,k ) )

Σ_k{1→q}( y_j,k - μ_j,k ) = 0
Σ_j{1→p}( y_j,k - μ_j,k ) = 0
Σ_k{1→q}( Σ_j{1→p}( y_j,k - μ_j,k ) ) = 0

Σ_k{1→q}( μ_j,k ) = Σ_k{1→q}( y_j,k )
Σ_j{1→p}( μ_j,k ) = Σ_j{1→p}( y_j,k )
Σ_k{1→q}( Σ_j{1→p}( μ_j,k ) ) = Σ_k{1→q}( Σ_j{1→p}( y_j,k ) )

Σ_k{1→q}( y_j,k ) = Y_A,j
Σ_j{1→p}( y_j,k ) = Y_B,k

Σ_j{1→p}( Y_A,j ) = N
Σ_k{1→q}( Y_B,k ) = N
Σ_k{1→q}( Σ_j{1→p}( y_j,k ) ) = N

Σ_k{1→q}( μ_j,k ) = Σ_k{1→q}( Y_A,jY_B,k / N ) = Y_A,j
Σ_j{1→p}( μ_j,k ) = Σ_j{1→p}( Y_A,jY_B,k / N ) = Y_B,k

という結果が得られ、先ほど得られた関係式を満たします。

飽和モデルの最大対数尤度を l(μ) とすると、それは各組み合わせに対するポアソン分布において μ_j,k = y_j,k として、その積の対数を求めればよいので

l(μ)	=	log( Π_k{1→q}( Π_j{1→p}( y_j,k^y_j,ke^-y_j,k / y_j,k! ) ) )
	=	Σ_k{1→q}( Σ_j{1→p}( y_j,klog y_j,k - y_j,k - log y_j,k! ) )

から求めることができます。実際に計算を行うと、l(μ) = -29.56 という結果が得られます。また加法モデルでは μ_j,k = Y_A,jY_B,k / N なので、最大対数尤度を l(a,b) で表すと

l(a,b) = Σ_k{1→q}( Σ_j{1→p}( y_j,klog( Y_A,jY_B,k / N ) - Y_A,jY_B,k / N - log y_j,k! ) )

から得られ、l(a,b) = -55.45 となります。従って、両者の差 ΔD は

ΔD = 2[ -29.56 - ( -55.45 ) ] = 51.78

です。飽和モデルの係数は 12 個あるので、l(μ) は自由度 12 の χ²分布に従い、一方、加法モデルの係数は 6 個なので自由度 6 の χ²分布に従います。従って、ΔD は自由度 12 - 6 = 6 のχ²分布に従い、ΔD = 51.78 のときの上側確率は 2.06E-9 と非常に小さく、加法モデルによる当てはめはよくない、すなわち各組み合わせが互いに独立であるという仮定が棄却されることになります。

加法モデルによる期待値と実測値を並べた表を以下に示します。() 内は実測値に対する差異を表しています。

表 3-3. 悪性黒色腫の病型・部位別計数の実測値と加法モデルによる期待値
病型	部位
	頭頸部		体幹部		四肢		計
	実測値	期待値	実測値	期待値	実測値	期待値	計
ハッチンソン黒色斑	22	5.78 (-16.22)	2	9.01 (+7.01)	10	19.21 (+9.21)	34
表在拡大型黒色腫	16	31.45 (+15.45)	54	49.02 (-4.98)	115	104.53 (-10.47)	185
結節型黒色腫	19	21.25 (+2.25)	33	33.12 (+0.12)	73	70.62 (-2.38)	125
不明	11	9.52 (-1.48)	17	14.84 (-2.16)	28	31.64 (+3.64)	56
計	68		106		226		400

実測値との差異が大きいのは、病型が「ハッチンソン黒色斑」と「表在拡大型黒色腫」の個所で、「ハッチンソン黒色斑」では頭頸部の部分で期待値に比べ頻発しており、「表在拡大型黒色腫」は相対的に頭頸部の頻度が小さくなっています。従って少なくともこの二つの病型については、部位との相関関係がある可能性があります。

事象ごとの発生回数が多項分布に従うとして検定を行う手法として「(ピアソンの)χ²-検定」を以前に紹介しました。下表のような、二種類の事象 A, B ごとの発生回数に対し、ピアソン・カイ二乗統計量 χ² = Σ_k{1→q}( Σ_j{1-p}( ( y_j,k - Y_A,jY_A,k / N )² / ( Y_A,jY_B,k / N ) ) ) が自由度 ( p - 1 )( q - 1 ) の χ²-分布に従うことを利用して、二つの事象 A, B の独立性を検定する手法です。

表 3-4. 二次元分割表
	A₁	A₂	...	A_j	...	A_p	計
B₁	y₁₁	y₂₁	...	y_j1	...	y_p1	Y_B,1
B₂	y₁₂	y₂₂	...	y_j2	...	y_p1	Y_B,2
:	:	:	...	:	...	:	:
B_k	y_1k	y_2k	...	y_jk	...	y_pk	Y_B,k
:	:	:	...	:	...	:	:
B_q	y_1q	y_2q	...	y_jq	...	y_pq	Y_B,q
計	Y_A,1	Y_A,2	...	Y_A,j	...	Y_A,p	N

加法モデルにおいて、実測値 y_j,k に対する期待値は Y_A,jY_B,k / N となるのでした。ピアソン・カイ二乗統計量は、実測値 o_i と期待値 e_i について

χ² = Σ_i( ( o_i - e_i )² / e_i )

で求められるので、加法モデルでのピアソン・カイ二乗統計量は χ²-検定でのそれと一致します。つまり、上記の例において独立性の検定を行った場合、χ²-検定とほぼ同じ結果が得られることになります。但し、対数線形モデルではピアソン・カイ二乗統計量だけではなく対数尤度統計量も検定に利用することができるところが異なります。ちなみに、ピアソン・カイ二乗統計量は表 3-3 の ()内にある実測値と期待値の差と、期待値そのものを使えば計算できて、その値は χ² = 65.81 となります。

対数線形モデルの場合、更に複雑な検定を行うこともできます。文献にある別のサンプル・データを以下に示します。

このデータは、胃潰瘍・十二指腸潰瘍の発症とアスピリン使用有無の関係を調査した時のデータです。ここでは、胃潰瘍・十二指腸潰瘍の発症有無ごとに 4 つのグループがあり、それぞれのサンプル数が固定となります。従って、発症の有無を表す添字 j、胃潰瘍・十二指腸潰瘍のいずれかを表す添字を k、アスピリンの使用有無を表す添字を l とし、それぞれの発生回数を y_j,k,l としたとき、4 つのグループのサンプル数 n_j,k は

n_j,k = Σ_l{1→2}( y_j,k,l )

Π_k{1→2}( Π_j{1→2}( n_j,k!Π_l{1→2}( θ_j,k,ly_j,k,l / y_j,k,l! ) ) )

但し、Σ_l{1→2}( θ_j,k,l ) = 1

log( E[y_j,k,l] ) = κ + αx_α + βx_β + γx_γ + λx_λ + η₁x_η1 + η₂x_η2

とし、κ は定数項、x_α は発症の有無を表す二値変数、x_β は胃潰瘍 ( = 0 ) と十二指腸潰瘍 ( = 1 ) を表すカテゴリ変数、x_γ は発症の有無と胃潰瘍・十二指腸潰瘍の相互作用 (十二指腸潰瘍の発症ありのみ 1 )、x_λ はアスピリン使用の有無、x_η1 と x_η2 はそれぞれ発症の有無及び胃潰瘍・十二指腸潰瘍とアスピリン使用有無の交互作用を表します。デザイン行列は以下のようになります。

X	=	\|	1,	0,	0,	0,	0,	0,	0	\|
		\|	1,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	0	\|
		\|	1,	1,	1,	1,	0,	0,	0	\|
		\|	1,	0,	0,	0,	1,	0,	0	\|
		\|	1,	1,	0,	0,	1,	1,	0	\|
		\|	1,	0,	1,	0,	1,	0,	1	\|
		\|	1,	1,	1,	1,	1,	1,	1	\|

また、従属変数ベクトル y = ( 62, 39, 53, 49, 6, 25, 8, 8 ) です。これらのデータを使ってスコア法によって処理すると以下の回帰式と回帰係数の分散値が得られます。

Estimated regression equation : y = -0.149896x0 + -0.260289x1 + 0.0726564x2 + 0.112942x3 + -1.82193x4 + 1.14288x5 + -0.700046x6
variance of a = ( 0.0165437, 0.0375758, 0.0325526, 0.0683647, 0.0948408, 0.123955, 0.119735 )

回帰式から得られる予測値と実測値の差は以下のようになります。

表 3-6. モデル 1 の計数予測値 *()内は実測値との差
部位		アスピリン		総計
部位		使用なし	使用あり	総計
胃潰瘍	なし	58.53(-3.47)	9.47(3.47)	68
胃潰瘍	あり	42.47(3.47)	21.53(-3.47)	64
十二指腸潰瘍	なし	56.47(3.47)	4.53(-3.47)	61
十二指腸潰瘍	あり	45.53(-3.47)	11.47(3.47)	57

このモデル式から x_η2 ( 胃潰瘍・十二指腸潰瘍とアスピリン使用有無の交互作用 )、x_η1 ( 発症の有無とアスピリン使用有無の交互作用 )、x_λ ( アスピリン使用の有無 ) の順で一つずつ変数を除去してスコア法にて回帰係数を求めると、以下のような結果になります。

■ 胃潰瘍・十二指腸潰瘍とアスピリン使用有無の交互作用を除去 (モデル 2)

Estimated regression equation : y = -0.11488x0 + -0.203573x1 + 9.31899e-16x2 + 3.55945e-11x3 + -2.10587x4 + 1.12505x5
variance of a = ( 0.0156496, 0.0343737, 0.0310993, 0.0642679, 0.0801242, 0.12179 )

■ 発症の有無とアスピリン使用有無の交互作用を除去 (モデル 3)

Estimated regression equation : y = -0.208255x0 + 1.4785e-10x1 + -5.77581e-11x2 + -1.50665e-10x3 + -1.46306x4
variance of a = ( 0.0156318, 0.0303302, 0.0310988, 0.064267, 0.0262006 )

■ アスピリン使用の有無を除去 (モデル 4)

Estimated regression equation : y = -0.693147x0 + -3.48683e-12x1 + -3.41539e-12x2 + 3.44528e-12x3
variance of a = ( 0.0147058, 0.0303308, 0.0310993, 0.0642681 )

回帰式から得られる予測値と実測値の差は以下の通りです。

表 3-7. 各モデル式での計数予測値 *()内は実測値との差
部位		アスピリン		総計
部位		使用なし	使用あり	総計
モデル 2
胃潰瘍	なし	60.62(-1.38)	7.38(1.38)	68
胃潰瘍	あり	46.55(7.55)	17.45(-7.55)	64
十二指腸潰瘍	なし	54.38(1.38)	6.62(-1.38)	61
十二指腸潰瘍	あり	41.45(-7.55)	15.55(7.55)	57
モデル 3
胃潰瘍	なし	55.22(-6.78)	12.78(6.78)	68
胃潰瘍	あり	51.97(12.97)	12.03(-12.97)	64
十二指腸潰瘍	なし	49.53(-3.47)	11.47(3.47)	61
十二指腸潰瘍	あり	46.28(-2.72)	10.72(2.72)	57
モデル 4
胃潰瘍	なし	34.00(-28.00)	34.00(28.00)	68
胃潰瘍	あり	32.00(-7.00)	32.00(7.00)	64
十二指腸潰瘍	なし	30.50(-22.50)	30.50(22.50)	61
十二指腸潰瘍	あり	28.50(-20.50)	28.50(20.50)	57

l	=	Σ_l{1→2}( Σ_k{1→2}( Σ_j{1→2}( log( μ_j,k,l^y_j,k,l e^-μ_j,k,l / y_j,k,l! ) ) ) )
	=	Σ_l{1→2}( Σ_k{1→2}( Σ_j{1→2}( y_j,k,llog μ_j,k,l - μ_j,k,l - log y_j,k,l! ) ) )

に対して求めた予測値 μ_j,k,l を代入して各モデルごとの差を計算すれば対数尤度統計量を得ることができます。各モデル式での回帰係数の数と対数尤度は以下の通りです。なお、飽和モデルでは μ_j,k,l = y_j,k,l なので、実測値をそのまま当てはめれば求めることができます。

表 3-8. 各モデル式での対数尤度
モデル式	回帰係数の数	対数尤度
飽和モデル	8	-19.81
モデル 1	7	-22.95
モデル 2	6	-25.08
モデル 3	5	-30.70
モデル 4	4	-83.16

発症の有無とアスピリン使用有無の交互作用を調べるにはモデル 2 とモデル 3 の対数尤度統計量を計算すればよいので、その値を D_η1 としたとき

D_η1 = 2( -25.08 + 30.70 ) = 11.25

となります。回帰係数の数の差は 1 なので D_η1 は自由度 1 のカイ二乗分布に従うことになり、11.25 に対する上側確率 (p値) は 0.001 より小さな値なので、発症がアスピリンの使用に影響を受けていると判断することができます。同様に、胃潰瘍・十二指腸潰瘍とアスピリン使用有無の交互作用を調べるにはモデル 2 とモデル 3 の対数尤度統計量を計算すればよく、その値を D_η2 とすれば

D_η2 = 2( -22.95 + 25.08 ) = 4.26

です。この値も自由度 1 のカイ二乗分布に従い、p 値が 0.039 となることからアスピリンの使用に対して胃潰瘍と十二指腸潰瘍の間でわずかながら差異があると判断できます。しかし、モデル 1 と飽和モデルの間の統計尤度統計量を D₀ としたとき、

D₀ = 2( -19.81 + 22.95 ) = 6.28

であり、飽和モデルに対してパラメータがわずか一つしか少なくないにもかかわらず比較的大きな値となっています。これは、このモデル式の当てはめがあまりよくないことを意味しています。

今回は、計数データに対して有効な一般化線形モデルとして「ポアソン回帰」と「対数線形モデル」と紹介し、特に分散分析と比較したときの相違点について細かく見てきました。こうして調べてみると、重回帰分析からの一連の推定・検定が一般化線形モデルという枠組みの中で一つにまとまることがわかってきて、なかなか興味深いです。

補足 1) 一元配置分散分析のデザイン行列の二乗 X^TX の逆行列

以下の行列 V に対する逆行列を M として、M を求める方法を検討してみます。

V = X^TX	=	\|	N,	n₁,	n₂,	...	n_p-1	\|
		\|	n₁,	n₁,	0,	...	0	\|
		\|	n₂,	0,	n₂,	...	0	\|
		\|	:	:	:	...	:	\|
		\|	n_p-2,	0,	0,	...	0	\|
		\|	n_p-1,	0,	0,	...	n_p-1	\|

対角成分が 1 でそれ以外は 0 になるような行列になる必要があることから、二行目以降の行ベクトル v_r^T ( 2 ≤ r ≤ p ) と M の列ベクトル m_c ( 1 ≤ c ≤ p ) の積は

v_r^Tm_c = δ_rc

である必要があります。但し δ_rc は「クロネッカーのデルタ (Kronecker Delta)」で、r = c のときのみ 1 でそれ以外は 0 になります。M の r 行 c 列の成分を m_rc として積を計算すると、

v_r^Tm_c = n_r-1m_1c + n_r-1m_rc = n_r-1( m_1c + m_rc )

なので、r ≠ c ならば m_rc = -m_1c、r = c ならば m_rc = ( 1 - n_r-1m_1c ) / n_r-1 を満たす必要があります。

v₁^Tm₁ = Nm₁₁ + Σ_r{2→p}( n_r-1m_r1 )

となりますが、N = Σ_r{1→p}( n_r ) であり、r ≠ 1 ならば m_r1 = -m₁₁ なので右辺は n_pm₁₁ であり、これが 1 に等しいことから

m₁₁ = 1 / n_p, m_r1 = -1 / n_p ( r ≥ 2 )

となることがわかります。ところで V は対称行列なので、その逆行列 M も対称行列です (*n1-1)。従って、

m_1r = m_r1 = -1 / n_p ( r ≥ 2 )

m_rc = -m_1c = 1 / n_p ( r ≠ c ; r ≥ 2 )

m_rr = ( 1 - n_r-1m_1r ) / n_r-1 = ( 1 + n_r-1 / n_p ) / n_r-1 = ( n_r-1 + n_p ) / n_r-1n_p

M	=	\|	1 / n_p,	-1 / n_p,	-1 / n_p,	...	-1 / n_p	\|
		\|	-1 / n_p,	( n₁ + n_p ) / n₁n_p,	1 / n_p,	...	1 / n_p	\|
		\|	-1 / n_p,	1 / n_p,	( n₂ + n_p ) / n₂n_p,	...	1 / n_p	\|
		\|	:	:	:	...	:	\|
		\|	-1 / n_p,	1 / n_p,	1 / n_p,	...	1 / n_p	\|
		\|	-1 / n_p,	1 / n_p,	1 / n_p,	...	( n_p-1 + n_p ) / n_p-1n_p	\|

補足 2) 二元配置分散分析のデザイン行列の二乗 X^TX の逆行列

二元配置分散分析の場合、X^TX は以下のようになるのでした。

X^TX	= n	\|	pq,	q1_p-1^T,	p1_q-1^T	\|	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|
		\|	q1_p-1,	qE_p-1,	N^T	\|	E_p-1,	E_p-1,	...	E_p-1	\|
		\|	p1_q-1,	N,	pE_q-1	\|	F₁^T,	F₂^T,	...	F_q-1^T	\|

		\|	1_p-1,	E_p-1,	F₁	\|	E_p-1,	0_p-1,	...	0_p-1	\|
		\|	1_p-1,	E_p-1,	F₂	\|	0_p-1,	E_p-1,	...	0_p-1	\|
		\|	:	:	:	\|	:	:	...	:	\|
		\|	1_p-1,	E_p-1,	F_q-1	\|	0_p-1,	0_p-1,	...	E_p-1	\|

但し、1_n は全要素が 1 の n 次元ベクトル、E_n は大きさ n の単位行列、0_n は大きさ n の正方零行列、F_k は k 列目のみが 1 で残りはゼロの p - 1 行 q - 1 列行列、N は全要素が 1 の q - 1 行 p - 1 列行列をそれぞれ表します。行列方向に q + 2 個ずつのブロックを持ったこの行列を 3 対 q - 1 個ずつのブロックに分け (上記の仕切り線を参照)、

X^TX	= n	\|	P,	Q^T	\|
		\|	Q,	E_(p-1)(q-1)	\|

とします。ここで、P は p + q - 1 の大きさの正方行列、Q は ( p - 1 )( q - 1 ) 行 p + q - 1 列の行列です。k ≠ 0 である任意の実数 k に対して、M が逆行列を持つならば明らかに ( kM )^-1 = (1/k)M^-1 なので、各ブロックの行列数が等しいように、X^TX の逆行列を以下のように定義します。

( X^TX )^-1	= (1/n)	\|	A,	B^T	\|
		\|	B,	C	\|

PA + Q^TB = E_p+q-1	--- (1)
QB^T + C = E_(p-1)(q-1)	--- (2)
QA + B = 0_{(p-1)(q-1),p+q-1}	--- (3)
PB^T + Q^TC = 0_{p+q-1,(p-1)(q-1)}	--- (4)

が成り立ちます。但し、0_r,c は r 行 c 列の零行列です。(3) より B = -QA を (1) に代入すると

PA - Q^TQA = ( P - Q^TQ )A = E_p+q-1

Q^TQ =

\|	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|\|	1_p-1,	E_p-1,	F₁	\|
\|	E_p-1,	E_p-1,	...	E_p-1	\|\|	1_p-1,	E_p-1,	F₂	\|
\|	F₁^T,	F₂^T,	...	F_q-1^T	\|\|	:	:	:	\|
					\|	1_p-1,	E_p-1,	F_q-1	\|

\|	( q - 1 )1_p-1^T1_p-1,	( q - 1 )1_p-1^T,	Σ_k{1→q-1}( 1_p-1^TF_k )	\|
\|	( q - 1 )1_p-1,	( q - 1 )E_p-1,	Σ_k{1→q-1}( F_k )	\|
\|	Σ_k{1→q-1}( F_k^T1_p-1 ),	Σ_k{1→q-1}( F_k^T ),	Σ_k{1→q-1}( F_k^TF_k )	\|

\|	( p - 1 )( q - 1 ),	( q - 1 )1_p-1^T,	( p - 1 )1_q-1^T	\|
\|	( q - 1 )1_p-1,	( q - 1 )E_p-1,	N^T	\|
\|	( p - 1 )1_q-1,	N,	( p - 1 )E_q-1	\|

P - Q^TQ	=	\|	p + q - 1,	1_p-1^T,	1_q-1^T	\|
		\|	1_p-1,	E_p-1,	0_p-1,q-1	\|
		\|	1_q-1,	0_q-1,p-1,	E_q-1	\|

となります。これはちょうど、補足 1 の行列にて n_i = 1 ( 1 ≤ i ≤ p + q - 1 ) とした場合に相当し、その逆行列は

( P - Q^TQ )^-1

=	\|	1,	-1,	-1,	...	-1	\|
	\|	-1,	2,	1,	...	1	\|
	\|	-1,	1,	2,	...	1	\|
	\|	:	:	:	...	:	\|
	\|	-1,	1,	1,	...	2	\|

=	\|	1,	-1_p-1^T,	-1_q-1^T	\|
	\|	-1_p-1,	N_p-1 + E_p-1,	N^T	\|
	\|	-1_q-1,	N,	N_q-1 + E_q-1	\|

となって ( N_n は大きさが n で全要素が 1 の正方行列とします )、これが A と等しくなります。これを (3) に代入すれば

B = -QA =

-	\|	1_p-1,	E_p-1,	F₁	\|\|	1,	-1_p-1^T,	-1_q-1^T	\|
	\|	1_p-1,	E_p-1,	F₂	\|\|	-1_p-1,	N_p-1 + E_p-1,	N^T	\|
	\|	:	:	:	\|\|	-1_q-1,	N,	N_q-1 + E_q-1	\|
	\|	1_p-1,	E_p-1,	F_q-1	\|

\|	F₁1_q-1,	-E_p-1 - F₁N,	-F₁( N_q-1 + E_q-1 )	\|
\|	F₂1_q-1,	-E_p-1 - F₂N,	-F₂( N_q-1 + E_q-1 )	\|
\|	:	:	:	\|
\|	F_q-11_q-1,	-E_p-1 - F_q-1N,	-F_q-1( N_q-1 + E_q-1 )	\|

\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₁ )	\|
\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₂ )	\|
\|	:	:	:	\|
\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F_q-1 )	\|

C =

E_(p-1)(q-1) - QB^T

E_(p-1)(q-1) -	\|	1_p-1,	E_p-1,	F₁	\|\|	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|
	\|	1_p-1,	E_p-1,	F₂	\|\|	-( N_p-1 + E_p-1 ),	-( N_p-1 + E_p-1 ),	...	-( N_p-1 + E_p-1 )	\|
	\|	:	:	:	\|\|	-( N + F₁^T ),	-( N + F₂^T ),	...	-( N + F_q-1^T )	\|
	\|	1_p-1,	E_p-1,	F_q-1	\|

E_(p-1)(q-1) -	\|	-( E_p-1 + 2N_p-1 ),	-( N_p-1 + E_p-1 ),	...	-( N_p-1 + E_p-1 )	\|
	\|	-( N_p-1 + E_p-1 ),	-( E_p-1 + 2N_p-1 ),	...	-( N_p-1 + E_p-1 )	\|
	\|	:	:	...	:	\|
	\|	-( N_p-1 + E_p-1 ),	-( N_p-1 + E_p-1 ),	...	-( E_p-1 + 2N_p-1 )	\|

\|	2( N_p-1 + E_p-1 ),	N_p-1 + E_p-1,	...	N_p-1 + E_p-1	\|
\|	N_p-1 + E_p-1,	2( N_p-1 + E_p-1 ),	...	N_p-1 + E_p-1	\|
\|	:	:	...	:	\|
\|	N_p-1 + E_p-1,	N_p-1 + E_p-1,	...	2( N_p-1 + E_p-1 )	\|

( X^TX )^-1	= (1/n)	\|	1,	-1_p-1^T,	-1_q-1^T,	1_p-1^T,	1_p-1^T,	...	1_p-1^T	\|
		\|	-1_p-1,	N_p-1 + E_p-1,	N^T,	-( N_p-1 + E_p-1 ),	-( N_p-1 + E_p-1 ),	...	-( N_p-1 + E_p-1 )	\|
		\|	-1_q-1,	N,	N_q-1 + E_q-1,	-( N + F₁^T ),	-( N + F₂^T ),	...	-( N + F_q-1^T )	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₁ )	2( N_p-1 + E_p-1 ),	N_p-1 + E_p-1,	...	N_p-1 + E_p-1	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F₂ )	N_p-1 + E_p-1,	2( N_p-1 + E_p-1 ),	...	N_p-1 + E_p-1	\|
		\|	:	:	:	:	:	...	:	\|
		\|	1_p-1,	-( N_p-1 + E_p-1 ),	-( N^T + F_q-1 )	N_p-1 + E_p-1,	N_p-1 + E_p-1,	...	2( N_p-1 + E_p-1 )	\|

											<	( p - 1 )( q - 1 ) 列											>
		<1 列>	<	p-1 列		>	<	q-1 列		>	<	p-1 列		>	<	p-1 列		>		<	p-1 列		>
= (1/n)	\|	1,	-1,	-1,	...	-1,	-1,	-1,	...	-1,	1,	1,	...	1,	1,	1,	...	1,	...	1,	1,	...	1	\|	1 行

	\|	-1,	2,	1,	...	1,	1,	1,	...	1,	-2,	-1,	...	-1,	-2,	-1,	...	-1,	...	-2,	-1,	...	-1	\|	p-1 行
	\|	-1,	1,	2,	...	1,	1,	1,	...	1,	-1,	-2,	...	-1,	-1,	-2,	...	-1,	...	-1,	-2,	...	-1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	-1,	1,	1,	...	2,	1,	1,	...	1,	-1,	-1,	...	-2,	-1,	-1,	...	-2,	...	-1,	-1,	...	-2	\|

	\|	-1,	1,	1,	...	1,	2,	1,	...	1,	-2,	-2,	...	-2,	-1,	-1,	...	-1,	...	-1,	-1,	...	-1	\|	q-1 行
	\|	-1,	1,	1,	...	1,	1,	2,	...	1,	-1,	-1,	...	-1,	-2,	-2,	...	-2,	...	-1,	-1,	...	-1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	-1,	1,	1,	...	1,	1,	1,	...	2,	-1,	-1,	...	-1,	-1,	-1,	...	-1,	...	-2,	-2,	...	-2	\|

	\|	1,	-2,	-1,	...	-1,	-2,	-1,	...	-1,	4,	2,	...	2,	2,	1,	...	1,	...	2,	1,	...	1	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-2,	-1,	...	-1,	2,	4,	...	2,	1,	2,	...	1,	...	1,	2,	...	1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-2,	-1,	...	-1,	2,	2,	...	4,	1,	1,	...	2,	...	1,	1,	...	2	\|

	\|	1,	-2,	-1,	...	-1,	-1,	-2,	...	-1,	2,	1,	...	1,	4,	2,	...	2,	...	2,	1,	...	1	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-1,	-2,	...	-1,	1,	2,	...	1,	2,	4,	...	2,	...	1,	2,	...	1	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-1,	-2,	...	-1,	1,	1,	...	2,	2,	2,	...	4,	...	1,	1,	...	2	\|

	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|

	\|	1,	-2,	-1,	...	-1,	-1,	-1,	...	-2,	2,	1,	...	1,	2,	1,	...	1,	...	4,	2,	...	2	\|	p-1 行
	\|	1,	-1,	-2,	...	-1,	-1,	-1,	...	-2,	1,	2,	...	1,	1,	2,	...	1,	...	2,	4,	...	2	\|
	\|	:	:	:	...	:	:	:	...	:	:	:	...	:	:	:	...	:	...	:	:	...	:	\|
	\|	1,	-1,	-1,	...	-2,	-1,	-1,	...	-2,	1,	1,	...	2,	1,	1,	...	2,	...	2,	2,	...	4	\|

補足 3) 二元配置分散分析のデザイン行列の二乗 X^TX の逆行列 (加法モデル)

二元配置分散分析のデザイン行列から交互作用効果を除いた場合、以下のような行列になるのでした。

X^TX	= n	\|	pq,	q1_p-1^T,	p1_q-1^T	\|
		\|	q1_p-1,	qE_p-1,	N^T	\|
		\|	p1_q-1,	N,	pE_q-1	\|

この逆行列を以下のような成分で表します。各ブロックの行列数は、X^TX における同じ行・列のものと等しいとします。

( X^TX )^-1	= (1/n)	\|	a,	b_p-1^T,	c_q-1^T	\|
		\|	b_p-1,	S_p-1,	T^T	\|
		\|	c_q-1,	T,	U_q-1	\|

a はスカラーであり、b_p-1 は次元数 p - 1、c_q-1 は次元数 q - 1 のベクトルです。S_p-1 は大きさが p - 1、U_q-1 は大きさが q - 1 の正方行列であり、T は q - 1 行 p - 1 列の行列です。この二つの行列の積が単位行列になることから、各ブロックの積和を求める式が 9 つ得られ、以下のような連立方程式が得られます。

pqa + q1_p-1^Tb_p-1 + p1_q-1^Tc_q-1 = 1	--- (1)
pqb_p-1^T + q1_p-1^TS_p-1 + p1_q-1^TT = 0_p-1^T	--- (2)
pqc_q-1^T + q1_p-1^TT^T + p1_q-1^TU_q-1 = 0_q-1^T	--- (3)

qa1_p-1 + qb_p-1 + N^Tc_q-1 = 0_p-1	--- (4)
q1_p-1b_p-1^T + qS_p-1 + N^TT = E_p-1	--- (5)
q1_p-1c_q-1^T + qT^T + N^TU_q-1 = 0_p-1,q-1	--- (6)

pa1_q-1 + Nb_p-1 + pc_q-1 = 0_q-1	--- (7)
p1_q-1b_p-1^T + NS_p-1 + pT = 0_q-1,p-1	--- (8)
p1_q-1c_q-1^T + NT^T + pU_q-1 = E_q-1	--- (9)

b_p-1, c_q-1 の j 番目の成分をそれぞれ b_j, c_j とすると、(1) より

pqa + qΣ_j{1→p-1}( b_j ) + pΣ_j{1→q-1}( c_j ) = 1 --- (1')

qa + qb_r + Σ_j{1→q-1}( c_j ) = 0 --- (4')

qΣ_j{1→p-1}( b_j ) - pqb_r = 1

が全ての r について成り立ちます。特に r ≠ s に対して上式が成り立つなら

qΣ_j{1→p-1}( b_j ) - pqb_r = 1
qΣ_j{1→p-1}( b_j ) - pqb_s = 1

を辺々引けば容易に b_r = b_s が成り立つことがわかります。従って、

qΣ_j{1→p-1}( b_j ) - pqb_r = q( p - 1 )b_r - pqb_r = -qb_r = 1

より b_p-1 = -1_p-1 / q となります。(1), (7) に対して同様の操作を行えば c_q-1 = -1_q-1 / p という結果が得られ、a = ( p + q - 1 ) / p であることもわかります。この結果を (2) に代入すれば、

-p1_p-1^T + q1_p-1^TS_p-1 + p1_q-1^TT = 0_p-1^T

-p + qΣ_i{1→p-1}( s_i,c ) + pΣ_i{1→q-1}( t_i,c ) = 0 --- (2')

-N_p-1 + qS_p-1 + N^TT = E_p-1

-1 + qs_r,c + Σ_i{1→q-1}( t_i,c ) = δ_rc --- (5')

となります。但し、s_r,c は S の、t_r,c は T の r 行 c 列の成分、δ_rc は「クロネッカーのデルタ (Kronecker Delta)」で、r = c のときのみ 1 でそれ以外は 0 になります。(2') - p x (5') より

qΣ_i{1→p-1}( s_i,c ) - pqs_r,c = -pδ_rc --- (A)

qΣ_i{1→p-1}( s_i,c ) - pqs_s,c = -pδ_sc
qΣ_i{1→p-1}( s_i,c ) - pqs_r,c = -pδ_rc

pq( s_r,c - s_s,c ) = p( δ_rc - δ_sc )

となるので、r ≠ c かつ r ≠ c ならば s_r,c = s_s,c であり、r = c ならば

pq( s_r,c - s_s,c ) = p より s_r,c = s_s,c + 1 / q

となって、対角成分は非対角成分より 1 / q だけ大きいことを意味します。(A) に対し、r ≠ c のとき、前半の和の中には一つだけ対角成分があるので

q[ ( p - 1 )s_r,c + 1 / q ] - pqs_r,c = 0 より

s_r,c = 1 / q

であり、対角成分は 2 / q です。よって、S_p-1 = ( 1 / q )( N_p-1 + E_p-1 ) という結果が得られ、(3) と (6) に同様の操作を行うことで U_q-1 = ( 1 / p )( N_q-1 + E_q-1 ) となります。b_p-1 と S_p-1 を (8) に代入すると

	-( p / q )1_q-11_p-1^T + ( 1 / q )N( N_p-1 + E_p-1 ) + pT
=	-( p / q )N + [ ( p - 1 ) / q ]N + ( 1 / q )N + pT
=	pT = 0_q-1,p-1

( X^TX )^-1	= (1/n)	\|	( p + q - 1 ) / pq,	-1_p-1^T / q,	-1_q-1^T / p	\|
		\|	-1_p-1 / q,	( 1 / q )( N_p-1 + E_p-1 ),	0_p-1,q-1	\|
		\|	-1_q-1 / p,	0_q-1,p-1,	( 1 / p )( N_q-1 + E_q-1 )	\|

補足 4) ポアソン分布の和の分布

二つの確率変数 x, y がそれぞれ互いに独立に確率分布 P_x(x), P_y(y) に従うとき、確率変数 x, y の和 u = x + y は畳み込み積分

P_x(u) * P_y(u)	=	∫{-∞→∞} p_x( u - v )p_y(v) dv	(連続分布)
	=	Σ_v{-∞→∞}( p_x( u - v )p_y(v) )	(離散分布)

となるのでした。よって、二つの確率変数 x, y がそれぞれ互いに独立にポアソン分布

P_λ₁(x) = e^-λ₁λ₁^x / x!

P_λ₂(y) = e^-λ₂λ₂^y / y!

P_λ₁(u) * P_λ₂(u)	=	Σ_v{0→u} p_λ₁( u - v )p_λ₂(v)
	=	Σ_v{0→u}( [ e^-λ₁λ₁^u-v / ( u - v )! ]・[ e^-λ₂λ₂^v / v! ] )
	=	e^{-(λ₁+λ₂)}Σ_v{0→u}( λ₁^u-vλ₂^v / ( u - v )!v! )
	=	e^{-(λ₁+λ₂)}Σ_v{0→u}( [ u! / ( u - v )!v! ]λ₁^u-vλ₂^v / u! )
	=	[ e^{-(λ₁+λ₂)} / u! ]Σ_v{0→u} _uC_vλ₁^u-vλ₂^v
	=	e^{-(λ₁+λ₂)}( λ₁ + λ₂ )^u / u!

となって、これはポアソン分布 P_λ₁+λ₂(u) そのものです。この結果から、ポアソン分布は「再生性(Reproductive Property)」を持つことがわかり、上記操作を繰り返せば明らかなように、複数の確率変数 x_i がポアソン分布 P_{λ_i}(x_i) に従うとき、その和 x = Σ_i( x_i ) の確率密度関数は

P_λ(x) = e^-λλ^x / x!

但し、λ = Σ_i( λ_i )

補足 5) 逸脱度残差が実数になることの証明

ポアソン分布を利用した場合「逸脱度残差」の平方根の中身は次のようになります。

y_ilog( y_i / m_i ) - ( y_i - m_i )

ここで y_i は実測値、m_i はモデル式から得られた当てはめ値をそれぞれ表します。この値が正値でなければ逸脱度残差は実数になりませんが、これが成り立つことは以下のように証明できます。

上式を m_i の関数 f( m_i ) とし ( y_i は定数とみなします )、導関数 f'( m_i ) を計算すると、

f'( m_i ) = -y_i / m_i + 1

よって、m_i = y_i のとき f( m_i ) は極値をとります。二階導関数 f''( m_i ) は

f''( m_i ) = y_i / m_i²

なので y_i > 0 ならば常に正値となり、f( m_i ) は極小値です。f( m_i ) = 0 より y_i > 0 のとき f( m_i ) ≥ 0 が成り立ちます。ポアソン分布では y_i > 0 なので、逸脱度残差は必ず実数となることになります。

X	=	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
		\|	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0	\|

X	=	\|	1,	1,	1,	1,	1	\|
		\|	1,	1,	2,	4,	2	\|
		\|	1,	1,	3,	9,	3	\|
		\|	1,	1,	4,	16,	4	\|
		\|	1,	1,	5,	25,	5	\|
		\|	1,	0,	1,	1,	0	\|
		\|	1,	0,	2,	4,	0	\|
		\|	1,	0,	3,	9,	0	\|
		\|	1,	0,	4,	16,	0	\|
		\|	1,	0,	5,	25,	0	\|

	喫煙の有無 (x_i1)	年齢 (x_i2)	年齢² (x_i3)	年齢 x 喫煙の有無 (x_i4)
係数の最尤推定量 a	1.441	2.376	-0.198	-0.308
a の分散	0.139	0.0432	7.49E-4	9.42E-3
a の標準偏差	0.372	0.208	0.0274	0.0970

年齢	喫煙	実測値	期待値	Χ_i	d_i	C_i
35～44	有	32	29.58	0.444	0.438	237.28
45～54		104	106.81	-0.272	-0.273	33.05
55～64		206	208.19	-0.152	-0.152	57.74
65～74		186	182.82	0.235	0.234	211.23
75～84		102	102.57	-0.0565	-0.0566	157.22
35～44	無	2	3.41	-0.766	-0.830	132.16
45～54		12	11.54	0.135	0.134	31.34
55～64		28	24.74	0.655	0.641	0.6039
65～74		28	30.23	-0.405	-0.410	20.05
75～84		31	31.07	-0.0125	-0.0125	52.75
				Χ²	D	C
				1.550	1.635	933.43

X	=	\|	1,	0,	0,	0,	0,	0,	0	\|
		\|	1,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	0	\|
		\|	1,	1,	1,	1,	0,	0,	0	\|
		\|	1,	0,	0,	0,	1,	0,	0	\|
		\|	1,	1,	0,	0,	1,	1,	0	\|
		\|	1,	0,	1,	0,	1,	0,	1	\|
		\|	1,	1,	1,	1,	1,	1,	1	\|

=	\|	1,	-1,	-1,	...	-1	\|
	\|	-1,	2,	1,	...	1	\|
	\|	-1,	1,	2,	...	1	\|
	\|	:	:	:	...	:	\|
	\|	-1,	1,	1,	...	2	\|

X	=	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
		\|	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0	\|

X	=	\|	1,	1,	1,	1,	1	\|
		\|	1,	1,	2,	4,	2	\|
		\|	1,	1,	3,	9,	3	\|
		\|	1,	1,	4,	16,	4	\|
		\|	1,	1,	5,	25,	5	\|
		\|	1,	0,	1,	1,	0	\|
		\|	1,	0,	2,	4,	0	\|
		\|	1,	0,	3,	9,	0	\|
		\|	1,	0,	4,	16,	0	\|
		\|	1,	0,	5,	25,	0	\|

X	=	\|	1,	0,	0,	0,	0,	0,	0	\|
		\|	1,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	0	\|
		\|	1,	1,	1,	1,	0,	0,	0	\|
		\|	1,	0,	0,	0,	1,	0,	0	\|
		\|	1,	1,	0,	0,	1,	1,	0	\|
		\|	1,	0,	1,	0,	1,	0,	1	\|
		\|	1,	1,	1,	1,	1,	1,	1	\|

=	\|	1,	-1,	-1,	...	-1	\|
	\|	-1,	2,	1,	...	1	\|
	\|	-1,	1,	2,	...	1	\|
	\|	:	:	:	...	:	\|
	\|	-1,	1,	1,	...	2	\|

X	=	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
		\|	1,	0,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

		\|	1,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
		\|	1,	0,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	1,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0,	1	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0,	0,	0	\|

→	\|	1,	1,	0,	1,	0,	1,	0,	0,	0	\|
	\|	1,	0,	1,	1,	0,	0,	1,	0,	0	\|
	\|	1,	0,	0,	1,	0,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	1,	0,	0,	1,	0	\|
	\|	1,	0,	1,	0,	1,	0,	0,	0,	1	\|
	\|	1,	0,	0,	0,	1,	0,	0,	0,	0	\|

	\|	1,	1,	0,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	1,	0,	0,	0,	0,	0,	0	\|
	\|	1,	0,	0,	0,	0,	0,	0,	0,	0	\|

X	=	\|	1,	1,	1,	1,	1	\|
		\|	1,	1,	2,	4,	2	\|
		\|	1,	1,	3,	9,	3	\|
		\|	1,	1,	4,	16,	4	\|
		\|	1,	1,	5,	25,	5	\|
		\|	1,	0,	1,	1,	0	\|
		\|	1,	0,	2,	4,	0	\|
		\|	1,	0,	3,	9,	0	\|
		\|	1,	0,	4,	16,	0	\|
		\|	1,	0,	5,	25,	0	\|

X	=	\|	1,	0,	0,	0,	0,	0,	0	\|
		\|	1,	1,	0,	0,	0,	0,	0	\|
		\|	1,	0,	1,	0,	0,	0,	0	\|
		\|	1,	1,	1,	1,	0,	0,	0	\|
		\|	1,	0,	0,	0,	1,	0,	0	\|
		\|	1,	1,	0,	0,	1,	1,	0	\|
		\|	1,	0,	1,	0,	1,	0,	1	\|
		\|	1,	1,	1,	1,	1,	1,	1	\|

=	\|	1,	-1,	-1,	...	-1	\|
	\|	-1,	2,	1,	...	1	\|
	\|	-1,	1,	2,	...	1	\|
	\|	:	:	:	...	:	\|
	\|	-1,	1,	1,	...	2	\|