確率・統計

(19) ロジスティック回帰 (Logistic Regression Model)

前章で、「一般化線形モデル (Generalized Linear Model)」の概要について紹介をしました。「線形重回帰モデル (Linear Multiple Regression Model)」は一般化線形モデルの一つであり、逆に考えれば、線形重回帰モデルを一般化したものが一般化線形モデルであるともいえます。今回は、一般化線形モデルの代表である「ロジスティック回帰 (Logistic Regression)」について紹介したいと思います。

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

1) ロジスティックモデル (Logistic Model)

ある試行が成功・失敗するというような事象を考えます。例えば、サイコロを二つ振ったとき、目の数が一致したら成功とするような事象や、コインを投げて表が出たら成功とする事象などです。この時の確率変数は二つのみであり、例えば

x=1(試行が成功)
=0(試行が失敗)

と表すことができます。このような変数は「二値確率変数 (Binary Random Variable)」と呼ばれます。試行が成功する確率 P(x=1) が π であるとしたとき、当然 0 ≤ π ≤ 1 であり、試行が失敗する確率 P(x=0) は 1 - π になります。この確率はまとめて

πx( 1 - π )1-x

と表すことができます。

このような試行を N 回繰り返したとき、j 回目 ( 1 ≤ j ≤ N ) の試行に対する確率変数 xj は 0, 1 のいずれかの値を取ります。j 回目の試行が成功する ( xj = 1 になる ) 確率を πj としたとき、同時確率は

Πj{1→N}( πjxj( 1 - πj )1-xj )=exp( logΠj{1→N}( πjxj( 1 - πj )1-xj ) )
=exp( Σj{1→N}( logπjxj( 1 - πj )1-xj ) )
=exp( Σj{1→N}( xjlogπj + ( 1 - xj )log( 1 - πj ) ) )
=exp( Σj{1→N}( xj[ logπj - log( 1 - πj ) ] + log( 1 - πj ) ) )
=exp( Σj{1→N}( xjlog( πj / ( 1 - πj ) ) + log( 1 - πj ) ) )

と表されます。これは、h(x) = 1, ηj(π) = log( πj / ( 1 - πj ) ), Tj(x) = xj, A(π) = -Σj{1→N}( log( 1 - πj ) ) としたとき N-母数指数型分布族に属します ( 但し、x, π は N 個の xj, πj からなるベクトルとします )。もし、全ての j に対して πj が等しければ、その値を π としたとき上式は

exp( Σj{1→N}( xjlog( π / ( 1 - π ) ) + log( 1 - π ) ) ) = exp( log( π / ( 1 - π ) )Σj{1→N}( xj ) + N・log( 1 - π ) )

と表すことができます。y = Σj{1→N}( xj ) としたとき、この値は試行が成功した回数を表し、同時確率は

exp( y・log( π / ( 1 - π ) ) + N・log( 1 - π ) )=[ π / ( 1 - π ) ]y( 1 - π )N
=πy( 1 - π )N-y

になります。y = Σj{1→N}( xj ) となる x の値は、N 個の要素から y 個を選ぶときの組み合わせ分だけあるので、y を確率変数とした時の周辺確率は、上記同時確率と組み合わせの数 NCy の積で求められ

NCyπy( 1 - π )N-y

という結果になります。これは、二項分布 BN,π(y) と同じ確率密度関数です。さらに一般的な場合として、K 個のグループがあって、各グループにおいて ni 回の試行を行い、成功した回数が yi 回だとしたとき、その同時確率は

Πi{1→K}( niCyiπiyi( 1 - πi )ni-yi )=exp( logΠi{1→K}( niCyiπiyi( 1 - πi )ni-yi ) )
=exp( Σi{1→K}( logniCyiπiyi( 1 - πi )ni-yi ) )
=exp( Σi{1→K}( logniCyi + yilogπi + ( ni - yi )log( 1 - πi ) ) )
=exp( Σi{1→K}( yilog( πi / ( 1 - πi ) ) + nilog( 1 - πi ) + logniCyi ) )

となり、対数尤度関数 l( π, y ) は

l( π, y ) = Σi{1→K}( yilog( πi / ( 1 - πi ) ) + nilog( 1 - πi ) + logniCyi )

と表されます。


二項分布 BN,π(y) の期待値は E[y] = Nπ であり、N 回の試行で y 回成功した時の割合は p = y / N なので、その期待値は E[p] = E[y] / N = π です。そこで、確率 π を連結関数 g(π) によって

g(π) = xTα

の形のモデル式に当てはめることを検討します。もし、g(π) = π (恒等関数) とした場合、π は x による線形式で表現できることになりますが、x の定義域が有限でなければ π そのものも有限となりません。特に、π は確率を表すので 0 ≤ π ≤ 1 である必要があり、x は有限の値を取らなければなりません。例えば、

π = α0 + α1x

というモデル式に対して、0 ≤ π ≤ 1 を満たすためには 0 ≤ α0 + α1x ≤ 1 でなければならないので、

0 / α1 ≤ x ≤ ( 1 - α0 ) / α1 [ α1 > 0 ]

( 1 - α0 ) / α1 ≤ x ≤ -α0 / α1 [ α1 < 0 ]

という制限が必要になります。そこで、0 ≤ π ≤ 1 を満たすために最適な連結関数として累積分布関数がよく用いられます。すなわち、任意の t に対して f(t) ≥ 0 を満たす関数 (これは確率密度関数が満たすべき条件の一つです) を使い、

P(x) = ∫{-∞→x} f(t) dt

と表します。但し、P(∞) = ∫{-∞→∞} f(t) dt = 1 となる必要があります (確率密度関数はこれも満たす必要があります)。例えば、区間 [ a, b ] 上の一様分布 Pa,b(t) の累積分布関数は

P(x)=0[ x < a ]
=( x - a ) / ( b - a )[ a ≤ x ≤ b ]
=1[ x > b ]

と表されるので、α0 = -a / ( b - a )、α1 = 1 / ( b - a ) とすれば π = α0 + α1x の形になり、連結関数を恒等関数とした時のモデル式と一致します。しかし前述の通り x の値が制限されるため、このモデル式はほとんど使われないようです。代わりに、正規分布 N( μ, σ2 ) を確率密度関数として用いたモデルを考えます。このときの累積分布関数は、

P(x) = ( 1 / ( 2πσ2 )1/2 ) ∫{-∞→x} exp( -( t - μ )2 / 2σ2 ) dt

となり、s = ( t - μ ) / σ とすると ds = dt / σ で、t → -∞ のとき s → -∞、t = x のとき s = ( x - μ ) / σ なので、

P(x)=( 1 / ( 2πσ2 )1/2 ) ∫{-∞→(x-μ)/σ} exp( -s2 / 2 )・σ ds
=( 1 / ( 2π )1/2 ) ∫{-∞→(x-μ)/σ} exp( -s2 / 2 ) ds
=Φ( ( x - μ ) / σ )

と表すことができます。但し、Φ(x) は標準正規分布 N( 0, 1 ) の累積分布関数を表します。P(x) = π とすれば、

Φ-1(π) = -μ / σ + x / σ

なので、α0 = -μ / σ、α1 = 1 / σ とすれば、これは連結関数 g(π) を標準正規分布の逆累積分布関数 Φ-1(π) とした一般化線形モデルとなります。さらに一般化して

Φ-1(π) = xTα

で表される一般化線形モデルを「プロビット・モデル (Probit Model)」といいます。このとき、π = P(x) は

π = P(x) = Φ( xTα ) = ( 1 / ( 2π )1/2 ) ∫{-∞→xTα} exp( -s2 / 2 ) ds

で計算することができます。

図 1-1. プロビット・モデル
正規分布の累積分布関数

プロビット・モデルはモデル式に積分を含むので、いわゆる「初等関数 (Elementary Function)」とは異なる関数です。このモデル式によく似た形状を持った初等関数として、「ロジスティック・モデル (Logistic Model)」または「ロジット・モデル (Logit Model)」があります。ロジスティック・モデルの場合、正規分布の代わりに以下のような関数を使います。

f(t) = et / ( 1 + et )2

s = et とすれば、f(s) = s / ( 1 + s )2、ds / dt = et = s で、t → -∞ のとき s → 0、t = xTα のとき s = exp( xTα ) なので、

P(x) = ∫{-∞→xTα} f(t) dt=∫{0→exp( xTα )} [ s / ( 1 + s )2 ]・[ 1 / s ] ds
=∫{0→exp( xTα )} ( 1 + s )-2 ds
=[ -( 1 + s )-1 ]{0→exp( xTα )}
=1 - 1 / [ 1 + exp( xTα ) ]
=exp( xTα ) / [ 1 + exp( xTα ) ]

と求めることができます。P(x) = 1 - 1 / [ 1 + exp( xTα ) ] より P(x) は単調増加であり、xTα → +∞ のとき P(x) → 1 となることから、f(t) は確率密度関数として成り立っていることになります。P(x) = π とすれば、

[ 1 + exp( xTα ) ]π = exp( xTα ) より

( 1 - π )exp( xTα ) = π
exp( xTα ) = π / ( 1 - π )
xTα = log( π / ( 1 - π ) )

となるので、連結関数を g(π) = log( π / ( 1 - π ) ) とすれば一般化線形モデルとなります。この連結関数は「ロジット関数 (Logit Function)」と一般的に呼ばれています。

図 1-2. ロジスティック・モデル
ロジット関数

その他に、「極値分布 (Extreme Distribution)」と呼ばれる確率密度関数の一つである「標準ガンベル分布 (Standard Gumbel Distribution)」を使った以下のようなモデル式もあります。

f(t) = exp( t - et )

f(t) = et・exp( -et ) より s = et とおくと、ds / dt = s、f(s) = se-s、t → -∞ のとき s → 0、t = xTα のとき s = exp( xTα ) なので、

P(x) = ∫{-∞→xTα} f(t) dt=∫{0→exp(xTα)} [ se-s ]・[ 1 / s ] ds
=∫{0→exp(xTα)} e-s ds
=[ -e-s ]{0→exp(xTα)}
=1 - exp( -exp( xTα ) )

となります。P(x) は単調増加であり、xTα → +∞ のとき P(x) → 1 となることから、極値分布は確率密度関数として成り立ち、P(x) = π とすれば、

exp( -exp( xTα ) ) = 1 - π より

xTα = log( -log( 1 - π ) )

となるので、連結関数を g(π) = log( -log( 1 - π ) ) とすれば一般化線形モデルとなります。この連結関数は「Complementary Log-log 関数」といいます。

図 1-3. 極値分布の累積分布
Complementary Log-log 関数

これまで紹介した累積分布関数は、グラフに表すと全て S 字形の曲線になるため、ギリシャ文字の ς に似た曲線という意味で「シグモイド曲線 (Sigmoid Curve)」と呼ばれます。「シグモイド関数 (Sigmoid Function)」という言葉もありますが、こちらは関数が一意に定義されていて、

ς(x) = eax / ( 1 + eax ) [ = 1 / ( 1 + e-ax ) ]

となります。これは、ロジスティック・モデルにおける累積分布関数の特殊形 (切片がゼロのモデル) と等しくなります。


まずは、今まで紹介した連結関数用のクラスを以下に示します。

/*
  ProbitModelFunc : プロビット・モデル連結関数
*/
class ProbitModelFunc : public LinkFunction_IF
{
  NormalDistribution norm_;

public:

  ProbitModelFunc()
    : norm_( 0, 1 ) {}

  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x < 0.5 )
      return( -binSearch( norm_, 0.5 - x ) );
    else
      return( binSearch( norm_, x - 0.5 ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  { return( exp( pow( (*this)( x ), 2 ) * 0.5 ) * sqrt( 2 * M_PI ) ); }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( norm_.lower_p( y ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Probit Model Function" ); }
};

/*
  LogitFunc : ロジット連結関数
*/
struct LogitFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x == 0 ) return( -INFINITY );
    if ( x == 1 ) return( INFINITY );
    return( log( x / ( 1 - x ) ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    return(
           ( x < 0 || x > 1 ) ? NAN :
           ( ( x == 0 || x == 1 ) ? INFINITY : 1 / ( x * ( 1 - x ) ) )
           );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( 1 - 1 / ( 1 + exp( y ) ) ); }

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

/*
  LoglogFunc : Complementary Log-log 連結関数
*/
struct LoglogFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x == 0 ) return( -INFINITY );
    if ( x == 1 ) return( INFINITY );
    return( log( -log( 1 - x ) ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    return(
           ( x < 0 || x > 1 ) ? NAN :
           ( x == 0 || x == 1 ) ? INFINITY : 1 / ( ( x - 1 ) * log( 1 - x ) )
           );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( 1 - exp( -exp( y ) ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Complementary Log-log Function" ); }
};

プロビット・モデルにおいて、連結関数の値のほかに、導関数と逆関数の値を返す処理が必要になります。プロビット・モデルは

x = Φ(y) = ( 1 / ( 2π )1/2 ) ∫{-∞→y} exp( -t2 / 2 ) dt

に対する逆関数 y = Φ-1(x) が連結関数なので、標準正規分布 N( 0, 1 ) の区間 -∞ ≤ t ≤ y における確率が x になるような y を「二分探索」や「ニュートン-ラフソン法」などを使って求める必要があります。区間 -∞ ≤ t ≤ y における確率は、連続確率密度関数を表すクラス ContDist 用のメンバ関数 lower_p で求めることができます (*1-1)。また、導関数に対しては、

dx / dy = ( 1 / ( 2π )1/2 ) exp( -y2 / 2 ) より

dy / dx = 1 / ( dx / dy ) = ( 2π )1/2exp( y2 / 2 )

となるので、やはり確率が x になるような y を求めてからそれを上式に代入して結果を得ます。最後の逆関数は簡単で、lower_p を使って区間 -∞ ≤ t ≤ y における確率を求めるだけです。

ロジスティック・モデルの場合、y = log( x / ( 1 - x ) ) なので連結関数の値を求めるのは簡単です。導関数は、

z = x / ( 1 - x ) = 1 / ( 1 - x ) - 1

y = logz

としたとき

dz / dx = 1 / ( 1 - x )2

なので、

dy / dx = ( dy / dz )( dz / dx ) = ( 1 / z )[ 1 / ( 1 - x )2 ] = 1 / x( 1 - x )

と求められ、逆関数は P(x) そのものなので

x = ey / ( 1 + ey )

になります。

最後の、Complementary Log-log 関数 y = log( -log( 1 - x ) ) の導関数は、

z = -log( 1 - x )

y = logz

としたとき

dz / dx = 1 / ( 1 - x )

なので、

dy / dx = ( dy / dz )( dz / dx ) = ( 1 / z )[ 1 / ( 1 - x ) ] = 1 / -( 1 - x )log( 1 - x )

と求められ、逆関数は

x = 1 - exp( -ey )

になります。

これらの関数群を、前回作成したスコア法用のサンプル・プログラムに適用すればいいわけですが、一つ問題となるのが π = 0, 1 になる時の処理で、ロジット関数と Complemenmtary Log-log 関数を利用した場合は計算途中で非数 (nan) が発生するため正しい結果が得られません。そのため、0 や 1 に近い値に置き換えて計算する必要があります。プロビット・モデルの場合は計算方式の理由で非数が発生しづらいようです (あくまでサンプル・プログラムの実装での話です)。


二項分布を一般化した場合、対数尤度関数は以下のような式で表されるのでした。

l( π|y ) = Σi{1→K}( yilog( πi / ( 1 - πi ) ) + nilog( 1 - πi ) + logniCyi )

飽和モデルは、パラメータ π を使ったモデル式であり、πi による偏微分が

∂l / ∂πi = yi / πi( 1 - πi ) - ni / ( 1 - πi ) = ( yi - niπi ) / πi( 1 - πi )

と求められることから、∂l / ∂πi = 0 のとき πi = yi / ni なので、対数尤度の最大値を l( π|y )max = l( p|y ) としたとき、

l( p|y ) = Σi{1→K}( yilog( yi / ( ni - yi ) ) + nilog( ( ni - yi ) / ni ) + logniCyi )

です。独立変数 xi = ( xi1, xi2, ... xip )T を使い、g(πi) = xiTα が成り立つような p 個のパラメータ α = ( α1, α2, ... αp )T が存在すれば、πi = g-1(xiTα) より ∂πi / ∂αj = xij / g'(πi) なので、αj による偏微分 uj = ∂l / ∂αj

uj = ∂l / ∂αj=Σi{1→K}( xijyi / g'(πii( 1 - πi ) - xijni / g'(πi)( 1 - πi ) )
=Σi{1→K}( ( yi - niπi )xij / g'(πii( 1 - πi ) )

となります。任意の指数型分布族 exp( η(θ)T(y) - A(θ) + B(y) ) に対し、uj = Σi{1→K}( ( yi - μi )xij / V[yi]g'(μi) ) で表されるのでした。ここでは μi = niπi、V[yi] = niπi( 1 - πi ) であり、g'(μi) = ∂g / ∂μi = ( ∂g / ∂πi )( ∂πi / ∂μi ) = g'(πi) / ni となるので、両式は一致します。yi = nipi と表せば ( つまり、実測値を使って確率を計算した結果を pi とすれば )、上式は

uj = Σi{1→K}( ni( pi - πi )xij / g'(πii( 1 - πi ) )

となるので、もし ni が全て等しければ、その値を N としたとき上式は

uj = NΣi{1→K}( ( pi - πi )xij / g'(πii( 1 - πi ) )

となり、連立方程式 uj = 0 において N は全て除去できるので、ちょうど N = 1 の二項分布 (すなわちベルヌーイ分布) を指数型分布族としたときの方程式と一致します (この場合の最尤推定量は確率そのものになります)。しかし、ni の偏りが大きい場合は、ベルヌーイ分布を使ってスコア法を適用することはできません。そのため、各試行回数に応じて分布が切り替えられるような仕組みが必要になります。


一般化した二項分布に対応したスコア法のサンプル・プログラムを以下に示します。

/*
  ExpFamily_Binomial : 一母数指数型分布族(二項分布)
*/
class ExpFamily_Binomial : public ExpFamily_IF
{
  double n_; // 総試行回数 N

public:

  // コンストラクタ
  ExpFamily_Binomial( unsigned int n )
    : n_( n ) {}

  // A(η)
  virtual double a( double eta ) const
  { return( n_ * log( exp( eta ) + 1 ) ); }

  // 期待値 ( E[y] = A'(η) )
  double average( double eta ) const
  { return( n_ * ( 1.0 - 1.0 / ( exp( eta ) + 1.0 ) ) ); }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  { return( log( mu / ( n_ - mu ) ) ); }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  { return( ( isinf( eta ) ) ? 0 : n_ * exp( eta ) / pow( exp( eta ) + 1, 2 ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Binomial Distribution" ); }
};

/*
  ExpFamily_MultiBinomial : 一母数指数型分布族(複数パラメータの二項分布)

  各 x の総試行回数に応じて式を切り替える
*/
class ExpFamily_MultiBinomial : public ExpFamily_IF
{
  static const double DEFAULT_H_ = 1E-6; // デフォルトの微小値

  vector<double> n_; // 各 x に対応する二項分布の総試行回数

  // 各関数の利用回数(n_.size()の剰余)
  mutable unsigned int aCnt_;      // a()
  mutable unsigned int aveCnt_;    // average()
  mutable unsigned int aveInvCnt_; // aveInv()
  mutable unsigned int varCnt_;    // variance()

  double h_; // 微小値 ( π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値 )

  // π = 0, 1 のときに計算結果がゼロになるのを防ぐためのチェック
  static double checkProb( double p )
  {
    return(
           ( p <= 0 ) ? DEFAULT_H_ :
           ( ( p >= 1 ) ? 1.0 - DEFAULT_H_ : p )
           );
  }

  // η = eta に対する π の値を返す
  static double eta2Pi( double eta )
  { return( checkProb( 1.0 - 1.0 / ( exp( eta ) + 1.0 ) ) ); }

  // π = pi に対する η の値を返す
  static double pi2Eta( double pi )
  { return( log( pi / ( 1.0 - pi ) ) ); }

public:

  // コンストラクタ
  ExpFamily_MultiBinomial( const vector<unsigned int>& n, double h = DEFAULT_H_ )
    : n_( n.size() ), aCnt_( 0 ), aveCnt_( 0 ), aveInvCnt_( 0 ), varCnt_( 0 ), h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
    for ( unsigned int i = 0 ; i < n_.size() ; ++i )
      n_[i] = n[i];
  }

  // A(η)
  virtual double a( double eta ) const
  {
    double ans = n_[aCnt_] * log( exp( eta ) + 1.0 );
    aCnt_ = ( aCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // 期待値 ( E[y] = A'(η) )
  virtual double average( double eta ) const
  {
    double ans = n_[aveCnt_] * eta2Pi( eta );
    aveCnt_ = ( aveCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  {
    double eta = pi2Eta( checkProb( mu / n_[aveInvCnt_] ) );
    aveInvCnt_ = ( aveInvCnt_ + 1 ) % n_.size();
    return( eta );
  }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  {
    double p = eta2Pi( eta );
    double ans = n_[varCnt_] * p * ( 1.0 - p );
    varCnt_ = ( varCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Binomial Distribution (multiple parameters)" ); }
};

/*
  SigmoidFunc : シグモイド型連結関数

  各 x の総試行回数に応じて式を切り替える
*/
class SigmoidFunc : public LinkFunction_IF
{
  static const double DEFAULT_H_ = 1E-6; // デフォルトの微小値

  const LinkFunction_IF& linkFunc_; // 利用する連結関数
  vector<double> n_; // 各 x に対応する二項分布の総試行回数

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

  double h_; // 微小値 ( π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値 )

  // checkProb : π = 0, 1 のときに計算結果がゼロになるのを防ぐためのチェック
  double checkProb( double p ) const
  {
    return(
           ( p <= 0 ) ? h_ :
           ( ( p >= 1 ) ? 1.0 - h_ : p )
           );
  }

public:

  /*
    コンストラクタ

    const LinkFunction_IF& linkFunc : 利用する連結関数
    const vector<unsigned int>& n : 各 x の総試行回数
    double h : 微小値
  */
  SigmoidFunc( const LinkFunction_IF& linkFunc, const vector<unsigned int>& n, double h = DEFAULT_H_ )
    : linkFunc_( linkFunc ), n_( n.size() ), opCnt_( 0 ), dfCnt_( 0 ), invfCnt_( 0 ), h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
    for ( unsigned int i = 0 ; i < n_.size() ; ++i )
      n_[i] = n[i];
  }

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

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    double d = linkFunc_.df( checkProb( x / n_[dfCnt_] ) ) / n_[dfCnt_];
    dfCnt_ = ( dfCnt_ + 1 ) % n_.size();
    return( d );
  }

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

  // 属性を表す文字列
  virtual string ident() const
  { return( linkFunc_.ident() ); }
};

/*
  ScoringMethod_Binomial : 二項分布を確率分布としたスコア法

  連結関数として Probit, Logistic, Complementary Log-log を利用することを想定

  const vector< vector<double> >& x : 独立変数(p個のパラメータのベクトルからなるn個のベクトル)
  const vector<unsigned int>& n : 試行回数
  const vector<unsigned int>& y : 成功した数
  vector<double>& a : 求めた係数
  const LinkFunction_IF& g : 連結関数
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)
  double h : 微小値 (π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値)

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool ScoringMethod_Binomial( const vector< vector<double> >& x, const vector<unsigned int>& n, const vector<unsigned int>& y,
               vector<double>& a, const LinkFunction_IF& g,
               bool verbose, unsigned int maxCount, double threshold, double h )
{
  if ( &n == 0 ) {
    cerr << "n not defined." << endl;
    return( false );
  }
  if ( &y == 0 ) {
    cerr << "y not defined." << endl;
    return( false );
  }
  if ( n.size() != y.size() ) {
    cerr << "n size (" << n.size() << ") and y size (" << y.size() << ") not matched." << endl;
    return( false );
  }
  ExpFamily_MultiBinomial pdf( n ); // 各 x 用の二項分布
  SigmoidFunc sf( g, n );           // 各 x 用の連結関数

  vector<double> yDbl( n.size() ); // y をdouble型に置き換える
  for ( unsigned int i = 0 ; i < n.size() ; ++i ) {
    if ( n[i] == 0 ) {
      cerr << "n[" << i << "] is zero ( n must be greater than zero )." << endl;
      return( false );
    }
    if ( y[i] > n[i] ) {
      cerr << "y[" << i << "] (" << y[i] << ") is more than n[" << i << "] (" << n[i] << ") ( y must be less than n )." << endl;
      return( false );
    }
    yDbl[i] = y[i];
    if ( y[i] == 0 ) yDbl[i] = h * (double)( n[i] );
    if ( y[i] == n[i] ) yDbl[i] -= h * (double)( n[i] );
  }

  return( ScoringMethod( x, yDbl, a, pdf, sf, verbose, maxCount, threshold ) );
}

ExpFamily_Binomial は二項分布を指数型分布族としたクラスで、ni が全て等しければこれをスコア法に用いることができます。ni が各グループに対して異なる場合は、ExpFamily_MultiBinomial クラスを利用します。このクラスの中で、グループごとの試行回数を保持しておき、スコア法の中では各試行回数を順番に使って計算する流れになりますが、スコア法ではメンバ関数の aveInvvariance を一回の計算の中で各グループに対し一回ずつ、合計で独立変数の数 n だけ呼び出しているため、各メンバ関数が呼び出された時にカウンタを 1 ずつ増加させることで切り替えを行っています。スコア法の実装に依存したやり方でありいい方法とは言えませんが、今回は妥協しました。なお、スコア法が反復処理であるため、一回のループでカウンタをリセットする必要があり、n の剰余をカウンタに反映する形で実現しています。

SigmoidFunc は、連結関数オブジェクト ProbitModelFunc, LogitFunc, LoglogFunc を利用して、試行回数を切り替えながらスコア法用のパラメータを求めるためのクラスで。連結関数オブジェクトのメンバ関数も、operator(), df, invf がスコア法の中でそれぞれ一回ずつ呼び出されているので、切り替えの原理は ExpFamily_MultiBinomial クラスと全く同じです。
ProbitModelFunc, LogitFunc, LoglogFunc では、引数や戻り値が πi であることを前提としています。しかし、スコア法の中では μi = niπi が利用されるため、SigmoidFunc の中で補正を行なっています。具体的には、operator()df の引数が μi であるため、ni で除算して連結関数へ渡すようにし、さらに

g'(μi) = g'(πi) / ni

より df では求めた値を ni で除算します。逆関数は

ξi = g(πi) より μi = niπi = nig-1i)

なので、invf で求めた値に ni を掛ければ μi を得ることができます。

ScoringMethod_Binomial 関数は、ExpFamily_MultiBinomialSigmoidFunc を使って、二項分布を指数型分布族としたスコア法を適用するための専用関数で、内部ではスコア法をそのまま利用しています。スコア法とは異なり、各グループの試行回数 n と成功した回数 y をそのまま別々に渡し、n を使って ExpFamily_MultiBinomialSigmoidFunc を初期化した上で、それらを使ってスコア法を行います。


サンプル・プログラムを利用して、文献にあったデータの回帰係数を計算してみたいと思います。下記データは、二硫化炭素 CS2 に 5 時間暴露されたカブトムシの死亡数を、使用したガス濃度ごとに示したものです。

図表 1-4. 二硫化炭素ガス濃度に対するカブトムシの死亡率
iCS2 ガス濃度 xi
(log10mgl-1)
カブトムシの数
ni
死亡数
yi
死亡率
pi=yi/ni
11.69075960.1017
21.724260130.2167
31.755262180.2903
41.784256280.5000
51.811363520.8254
61.836959530.8983
71.861062610.9839
81.883960601.0000
CS2 vs 死亡率グラフ

このデータから、まずは独立変数を CS2 ガス濃度 xi、従属変数を死亡率 pi = yi / ni とし、

g(π) = a0 + a1xi

の連結関数 g(π) にプロビット・モデル、ロジット関数、Complementary Log-log 関数をそれぞれ適用します。また、π が従う指数型分布族を N = 1 のときの二項分布 (つまりベルヌーイ分布) として回帰係数を求めると、結果は以下のようになります (但し、p8 = 1 になるため、1 に近い値として 0.999999 に置き換えて計算します)。

プロビット・モデル : π = -34.80 + 19.65x

ロジスティック・モデル : π = -60.46 + 34.12x

Complementary Log-log 関数 : π = -39.53 + 22.02x

図表 1-5. 実測値と当てはめ値の比較 (ベルヌーイ分布)
ガス濃度カブトムシの数実測値ProbitLogitLoglog
死亡率死亡数死亡率死亡数死亡率死亡数死亡率死亡数
1.6907590.101760.05733.380.05903.480.09475.59
1.7242600.2167130.178910.740.16429.850.187811.27
1.7552620.2903180.378223.450.361422.410.337320.92
1.7842560.5000280.602433.740.603533.800.541330.31
1.8113630.8254520.785949.510.793349.980.757147.70
1.8369590.8983530.902453.240.901953.210.916854.09
1.8610620.9839610.961559.620.954459.170.985461.10
1.8839601.0000600.986859.210.978658.720.999159.95
実測値とモデル式の比較(ベルヌーイ)

グラフを見る限り、Complementary Log-log 関数を使った場合が最も適合しているようです。次に、ScoringMethod_Binomial 関数を各連結関数と組み合わせて処理を行った結果を以下に示します。

プロビット・モデル : π = -34.94 + 19.73x

ロジスティック・モデル : π = -60.72 + 34.27x

Complementary Log-log 関数 : π = -39.57 + 22.04x

図表 1-6. 実測値と当てはめ値の比較 (グループ単位の二項分布)
ガス濃度カブトムシの数実測値ProbitLogitLoglog
死亡率死亡数死亡率死亡数死亡率死亡数死亡率死亡数
1.6907590.101760.05693.360.05863.45750.09475.59
1.7242600.2167130.178710.720.16409.840.188011.28
1.7552620.2903180.378723.480.362122.450.338020.95
1.7842560.5000280.603833.820.605333.900.542330.37
1.8113630.8254520.787549.620.795250.100.758447.78
1.8369590.8983530.903753.320.903253.290.917754.14
1.8610620.9839610.962359.660.955259.220.985761.11
1.8839601.0000600.987159.230.979058.740.999159.95

ベルヌーイ分布を利用した場合と比較すると非常に近い値となっています。しかし、後半部分のデータを ni, yi ともに 1000 倍すると、ScoringMethod_Binomial 関数を利用した場合、結果は以下のように変化します。

プロビット・モデル : π = -40.90 + 23.06x

ロジスティック・モデル : π = -79.74 + 44.83x

Complementary Log-log 関数 : π = -32.59 + 18.26x

図表 1-7. 実測値と当てはめ値の比較 (グループ単位の二項分布)
ガス濃度カブトムシの数実測値ProbitLogitLoglog
死亡率死亡数死亡率死亡数死亡率死亡数死亡率死亡数
1.6907590.101760.02821.670.01911.120.16359.65
1.7242600.2167130.12837.700.08024.810.280416.82
1.7552620.2903180.337420.920.259316.080.439827.27
1.7842560.5000280.598433.510.562331.490.626235.07
1.811363000.825452000.809050970.812451180.80095045
1.836959000.898353000.928554780.931754970.92385451
1.861062000.983961000.978360660.975760490.98166086
1.883960001.000060000.994659680.991259470.99775986
実測値とモデル式の比較(グループ単位の二項分布)

グラフを見るとはっきりと分かるように、前半部分の当てはめは後半部分に比べて非常に悪くなります。uj = Σi{1→K}( ni( pi - πi )xij / g'(πii( 1 - πi ) ) より、ni が小さい成分については影響度が小さくなることが主な理由と考えられます。元々のデータでは ni が互いに非常に近い値だったので、その場合はベルヌーイ分布を利用した近似解でも問題はありませんが、そうでないときは一般化した二項分布で計算する必要があることがこの結果からもわかります。


サンプル・プログラムを利用して uj = 0 ( j = 1, 2, ... p ) を満たす α を求め、y^i = niπ^i = nig-1(xiTα) と計算すれば、π^i は πi の最尤推定量であり、y^i は yi に対する当てはめ値なので、これを対数尤度に代入すれば l( α|y ) の最大値 l( α|y )max = l( a|y ) が得られ、その式は

l( a|y )=Σi{1→K}( yilog( π^i / ( 1 - π^i ) ) + nilog( 1 - π^i ) + logniCyi )
=Σi{1→K}( yilog( y^i / ( ni - y^i ) ) + nilog( ( ni - y^i ) / ni ) + logniCyi )

となります。従って、対数尤度統計量 D は

D = 2[ l( p|y ) - l( a|y ) ]=2[ Σi{1→K}( yilog( yi / ( ni - yi ) ) + nilog( ( ni - yi ) / ni ) + logniCyi )
 - Σi{1→K}( yilog( y^i / ( ni - y^i ) ) + nilog( ( ni - y^i ) / ni ) + logniCyi ) ]
=i{1→K}( yilog( yi( ni - y^i ) / y^i( ni - yi ) ) + nilog( ( ni - yi ) / ( ni - y^i ) ) )
=i{1→K}( yilog( yi / y^i ) - yilog( ( ni - yi ) / ( ni - y^i ) ) + nilog( ( ni - yi ) / ( ni - y^i ) ) )
=i{1→K}( yilog( yi / y^i ) + ( ni - yi )log( ( ni - yi ) / ( ni - y^i ) ) )

となって、yi の実測値と当てはめ値を使えば各モデル式での対数尤度統計量を得ることができます。yi と ni - yi はそれぞれ施行に成功・失敗した回数の実測値であるのに対し、y^i = niπi と ni - y^i はモデル式から得られた確率密度関数による期待値なので、実測値を o、期待値を e で表せば、上式は

2Σ( o・log( o / e ) )

の形で表せることになります。ここでの Σ は、各試行における成功・失敗回数それぞれ計 2K 個の和を取るという意味になります。これを、別の式に置き換えます。

Χ2 = Σ( ( o - e )2 / e )

この式は、「ピアソン・カイ二乗統計量 (Pearson Chi-squared Statistic)」と呼ばれ、「χ2-検定」でも登場しています。o, e にそれぞれの値を代入すると、

Χ2=Σi{1→K}( ( yi - y^i )2 / y^i ) + [ ( ni - yi ) - ( ni - y^i ) ]2 / ( ni - y^i ) )
=Σi{1→K}( [ ( ni - y^i )( yi - y^i )2 + y^i( y^i - yi )2 ] / y^i( ni - y^i ) )
=Σi{1→K}( ni( yi - y^i )2 / niy^i( 1 - πi ) )
=Σi{1→K}( ( yi - niπi )2 / niπi( 1 - πi ) )

となります。f(s) = s・log( s / t ) に対して t のまわりのテーラー展開を使って二次項までの近似式を求めると、

f(t) = t・log( t / t ) = 0

f'(s) = log( s / t ) + s・[ 1 / ( s / t ) ]・( 1 / t ) = log( s / t ) + 1 より f'(t) = 1

f(2)(s) = [ 1 / ( s / t ) ]・( 1 / t ) = 1 / s より f(2)(t) = 1 / t

なので

f(s) ≈ ( s - t ) + ( s - t )2 / 2t

であり、これを D 値に適用すれば

2Σ( o・log( o / e ) )2Σ( ( o - e ) + ( o - e )2 / 2e )
=Σ( [ e + ( o - e ) ]( o - e ) / e )
=Σ( o( o - e ) / e )
=Σi{1→K}( yi( yi - y^i ) / y^i +
 ( ni - yi )[ ( ni - yi ) - ( ni - y^i ) ] / ( ni - y^i ) )
=Σi{1→K}( yi( ni - y^i )( yi - y^i ) / y^i( ni - y^i ) +
 y^i( ni - yi )( y^i - yi ) / y^i( ni - y^i ) )
=Σi{1→K}( [ yi( ni - y^i ) - y^i( ni - yi ) ]( yi - y^i ) / y^i( ni - y^i ) )
=Σi{1→K}( ni( yi - y^i )2 / y^i( ni - y^i ) )
=Σi{1→K}( ( yi - niπi )2 / niπi( 1 - πi ) ) = χ2

となるので、実測値と当てはめ値が非常に近ければ、Χ2 値は、D 値の近似値となります。D が漸近的に自由度 K - p の χ2-分布に従うことから Χ2 もそれに従い、D の代わりに Χ2 値を利用することも可能で、D の方が、値の小さな yi, ni - yi の影響をより強く受けるため Χ2 値の方がよい近似値になるそうです。

D 値、Χ2 値それぞれに対して、「残差 (Residual)」が定義されます。まず、D 値に対する残差は「逸脱度残差 (Deviance Residuals)」と呼ばれ、

di = sign( yi - y^i ){ 2[ yilog( yi / y^i ) + ( ni - yi )log( ( ni - yi ) / ( ni - y^i ) ) ] }1/2

と定義されます。但し、sign(x) は x の正負を表し、実測値と当てはめ値の大小関係によって決まります。di2 は D の和の成分そのものなので、D = Σi{1→K}( di2 ) が成り立ちます。また、Χ2 値に対しては「ピアソン残差 (Pearson Residuals)」

Χi = ( yi - niπi ) / [ niπi( 1 - πi ) ]1/2

で定義され、これも Χ2 = Σi{1→K}( Χi2 ) が成り立ちます。試行回数 N が十分大きければ、二項分布 BN,π(y) に従う確率変数 y に対し、( y - Nπ ) / [ Nπ( 1 - π ) ]1/2 は漸近的に標準正規分布 N( 0, 1 ) に従います(*1-2)。従って、ピアソン残差 Χi は N( 0, 1 ) に漸近的に従い、Χi ≈ di より逸脱度残差 di についてもそれは成り立ちます。よって、それぞれの残差に対して「標準化残差」

rP[i] = Χi / ( 1 - hi )1/2

rD[i] = di / ( 1 - hi )1/2

を定義することができます。但し、hi は「てこ比 (Levarage)」を表します。


対数尤度統計量 D は、飽和モデルとの比較を意味するのに対し、全ての πi が等しいと仮定した場合のモデルである「最小モデル (Minimal Model)」と比較することも行われます。最小モデルは、対数尤度関数において πi = π (一定) とすればよいので、

l( π|y ) = Σi{1→K}( yilog( π / ( 1 - π ) ) + nilog( 1 - π ) + logniCyi )

より ∂l / ∂π を計算すると、

∂l / ∂π=Σi{1→K}( yi / π( 1 - π ) - ni / ( 1 - π ) )
=Σi{1→K}( ( yi - niπ ) / π( 1 - π ) )

なので、∂l / ∂π = 0 のとき π = Σi{1→K}( yi ) / Σi{1→K}( ni ) であり、これを p^ で表せば、l( π|y )max = l( p^|y ) より

C = 2[ l( a|y ) ) - l( p^|y ) ]=2[ Σi{1→K}( yilog( y^i / ( ni - y^i ) ) + nilog( ( ni - y^i ) / ni ) + logniCyi )
 - Σi{1→K}( yilog( p^ / ( 1 - p^ ) ) + nilog( 1 - p^ ) + logniCyi ) ]
=2[ Σi{1→K}( yilog( y^i( 1 - p^ ) / p^( ni - y^i ) ) + nilog( ( ni - y^i ) / ni( 1 - p^ ) ) ) ]
=2[ Σi{1→K}( -yilog( p^( ni - y^i ) / y^i( 1 - p^ ) ) + nilog( ( ni - y^i ) / ni( 1 - p^ ) ) ) ]
=2[ Σi{1→K}( -yilog( ( nip^ / y^i )・[ ( ni - y^i ) / ni( 1 - p^ ) ] ) + nilog( ( ni - y^i ) / ni( 1 - p^ ) ) ) ]
=2[ Σi{1→K}( yilog( y^i / nip^ ) + ( ni - yi )log( ( ni - y^i ) / ni( 1 - p^ ) ) ) ]

が最小モデルとの対数尤度統計量になります。C は「尤度比カイ二乗統計量 (Likelihood Ratio Chi-squared Statistc)」と呼ばれ、

C=2[ l( a|y ) - l( p^|y ) ]
=2{ [ l( a|y ) - l( α|y ) ] - [ l( p^|y ) - l( π|y ) ] - [ l( α|y ) - l( π|y ) ] }

より第一項は自由度 p の、第二項は自由度 1 の χ2-分布に従うので、もし l( α|y ) - l( π|y ) ≈ 0 ならば、C は自由度 p - 1 の χ2-分布に従います。この場合、最小モデルが対象のモデルとよく当てはまることを意味し、さらには切片以外の係数が全てゼロに非常に近いということになります。逆に、C の値がありえないほど大きいということは、係数がゼロではありえない、つまり傾きが意味のあるものであると判断することができます。

線形重回帰モデルにおいて、飽和モデルと最小モデルの間の対数尤度統計量を D0、飽和モデルと関心のあるモデル式の間の対数尤度統計量を D としたとき、

R2 = ( D0 - D ) / D0 = [ Σi{1→N}( ( yi - my )2 ) - Σi{1→N}( ( yi - y^i )2 ) ] / Σi{1→N}( ( yi - my )2 )

が「決定係数」と等しくなることを前章で示しました(*1-3)。これによく似た式として、

R2 = [ l( p^|y ) - l( a|y ) ] / l( p^|y ) = 1 - l( a|y ) / l( p^|y )

を「(マクファデンの) 擬似 R2 値 ( (McFadden's) Pseudo R-Squared)」といい、アメリカの経済学者「ダニエル・マクファデン (Daniel Little McFadden)」によって提唱されています。尤度は 1 以下であることから対数尤度は必ず負値となり、飽和モデルが対数尤度の最大値であったことから考えれば、最小モデルの対数尤度 l( π|y ) は必ず最小となる (つまり、その絶対値は最大となる) ことから、この擬似 R2 値は必ず 1 以下となります。

ここで注意すべき点として、擬似 R2 値の計算において、対数尤度は定数項部分の logniCyi を無視する必要があります。すなわち、この時の尤度は

L( π|y ) = Πi{1→K}( πiyi( 1 - πi )ni-yi )

であり、対数尤度は

l( π|y ) = Σi{1→K}( yilog( πi / ( 1 - πi ) ) + nilog( 1 - πi ) )

で表されます。対数尤度統計量は対数尤度の差で表されるので、定数項は打ち消し合って無視することができるのに対し、擬似 R2 値は定数項の有無によって値が変化します。文献や統計ソフトを確認する限り、擬似 R2 値の計算では定数項は除いてあります。しかし、対数尤度の式には定数項は含まれているので、擬似 R2 値の計算以外では定数項も含めて計算をしています。なお、ni が全て 1 で yi が 0 か 1 のいずれかである場合は、ベルヌーイ分布に従う Σi{1→K}( ni ) 回の独立試行を繰り返すモデルを表し、このときは定数項はゼロになります。


表 1-8. 各統計量の比較
ガス濃度カブトムシの数
ni
実測値 yiProbitLogitLoglog
死亡率死亡数diΧiCidiΧiCidiΧiCi
1.6907590.101761.341.4831.941.281.4132.020.180.1832.83
1.7242600.2167130.750.7718.551.061.1018.270.560.5718.67
1.7552620.290318-1.46-1.4411.49-1.20-1.1811.85-0.80-0.7912.24
1.7842560.500028-1.57-1.590.03-1.59-1.61-0.01-0.63-0.641.06
1.8113630.8254520.750.736.890.610.596.991.291.246.34
1.8369590.898353-0.14-0.1412.80-0.13-0.1312.80-0.52-0.5412.67
1.8610620.9839611.000.8925.961.251.0925.68-0.12-0.1226.46
1.8839601.0000601.250.8829.381.591.1328.880.330.2330.10
総計Σ(ni)p^Σ(yi)DΧ2CDΧ2CDΧ2C
4810.605029110.129.51274.0811.2310.03272.973.453.29280.76

上の表は、各連結関数を使って求めた当てはめ値を元に、逸脱度残差 di、ピアソン残差 Χi の他、尤比度カイ二乗統計量の和の成分 Ci を求めた結果です。表の下に総計として、対数尤度統計量 D、ピアソン・カイ二乗統計量 Χ2、尤比度カイ二乗統計量 C を示しています。D, Χ2Complementary Log-log 関数を使った場合が最も小さくなるため、この関数が最も適合していることが定量的に示されたことになります。これらは自由度 6 のカイ二乗分布に近似的に従うと考えてよいので、上側 5% 点が 12.59 であることから、プロビット・モデルとロジスティック・モデルに対してはあまり当てはまりがよいとは言えません。C は自由度 1 のカイ二乗分布に近似的に従い、求められた値に対する p 値はありえないほど小さいため、死亡率がガス濃度に無関係に全て等しいとする最小モデルは適用できず、傾きは必要であると判断することができます。

表 1-9. 擬似 R2 値の比較
ガス濃度li
ProbitLogitLoglog最小モデル
1.6907-20.30-20.22-19.41-52.24
1.7242-31.64-31.92-31.51-50.19
1.7552-38.42-38.07-37.67-49.91
1.7842-40.05-40.09-39.02-40.08
1.8113-29.46-29.36-30.01-36.35
1.8369-19.41-19.41-19.54-32.21
1.8610-5.62-5.90-5.13-31.58
1.8839-0.78-1.27-0.05-30.15
-185.68-186.24-182.34-322.72
擬似 R20.42460.42290.4350

上表は、各連結関数ごとの擬似 R2 値を計算した結果です。li は、対数尤度を計算するときの和の各成分を表しており、下側に示した合計は対数尤度そのものを意味します。各連結関数に対する対数尤度と最小モデルの対数尤度を使えば擬似 R2 値が得られ、それは表の最も下側に示されています。但し前述の通り、対数尤度の計算において定数項は含んでいないことに注意して下さい。どの連結関数に対しても R2 は 40% を少し超える程度です。

*1-1)確率・統計 (2) 確率空間」の中に「サンプル・プログラム」があり、確率密度関数用クラスが定義されています。lower_p は、連続分布用クラス ContDist の中で定義された純粋仮想関数です。正規分布に対する実装内容は、NormalDistribution クラスのサンプル・プログラムをご覧ください。

*1-2)確率・統計 (5) 正規分布」の「ド・モアブル = ラプラスの定理」を参照。

*1-3)確率・統計 (18) 一般化線形モデル」の「5) 尤度比」を参照。


2) ホズマー・レメショウ検定 (Hosmer-Lemeshow Test)

二項分布を利用した回帰分析は、ある独立変数ベクトル xi に対して二値の観測値が ni 個あって、試行に成功した個数が yi としたときに、pi = yi / ni を元にモデル式を決めるという方式をとっていました。しかし、xi が連続値であれば、全ての試行に対して xi が異なるというのが通常となって、ni はほとんど 1 となり、yi は 0 か 1 のいずれかしかとらなくなります。全ての i に対して ni = 1 の場合、飽和モデルに対する対数尤度 l( p|y ) の和の成分は必ず 0 になるので l( p|y ) = 0 です。l( p|y ) が χ2-分布に漸近的に従うことを利用していることから分かるように、この場合は D や Χ2 値とカイ二乗値との近似性はあまりよくなく、そのまま用いることはできません。そこで、次のような処理を行います。

  1. 一般化線形モデルから yi の当てはめ値 (この場合は確率に対する当てはめ値) π^i を求める。
  2. π^i の大きさに基づいて、当てはめ値を複数のグループに分割する。
  3. 各グループに対する試行回数 ni とその成功回数 yi を求め、その結果を元にピアソン・カイ二乗統計量 Χ2 を計算する。

グループに分割するときの目安としては、通常 10 個程度のグループにすることが多いようです。また、各グループの試行回数 ni はできるだけ等しくなるようにします。このようにして求めたピアソン・カイ二乗統計量 Χ2 を「ホズマー・レメショウ統計量 (Hosmer-Lemeshow Statistic)」といい、Χ2HL で表されます。数値実験の結果から、Χ2HL は自由度を ( グループ数 - 2 ) とするカイ二乗分布に従うことが知られています。

この方法は、二人の生物統計学者「David W. Hosmer」と「Stanley Lemeshow」によって 1980 年に発表されました。


ホズマー・レメショウ検定を行うためのサンプル・プログラムを以下に示します。

/*
  vector用 less 関数オブジェクト
*/
template<class T> struct LessVector
{
  bool operator()( const vector<T>& v1, const vector<T>& v2 )
  { for ( unsigned int i = 0 ; i < v1.size() && i < v2.size() ; ++i )
      if ( v1[i] != v2[i] ) return( v1[i] < v2[i] );
    return( false );
  }
};

/*
  NullCheck : NULLチェック関数

  T& t : 対象の変数へのリファレンス
  string arg : 変数名(出力用)
*/
template<class T> bool NullCheck( T& t, string arg )
{
  if ( &t != 0 ) return( true );

  cerr << arg << " not defined." << endl;
  return( false );
}

/*
  SizeCheck : コンテナクラスのサイズチェック

  T& t : 対象のコンテナクラス
  string arg1 : t の変数名(出力用)
  unsigned int sz : サイズ
  string arg2 : sz の大きさを持つ変数の名称(出力用)
*/
template<class T> bool SizeCheck( T& t, string arg1, unsigned int sz, string arg2 )
{
  if ( t.size() == sz ) return( true );

  cerr << arg1 << " size (" << t.size() << ") and " << arg2 << " size (" << sz << ") not matched." << endl;
  return( false );
}

/*
  CalcPredictiveValue : 予測値の計算

  const vector< vector<double> >& x : 独立変数
  const vector<double>& a : 係数
  vector<double>& y : 計算した予測値
  const LinkFunction_IF& g : 連結関数
*/
void CalcPredictiveValue( const vector< vector<double> >& x, const vector<double>& a, vector<double>& y,
                          const LinkFunction_IF& g )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return;
  if ( ! NullCheck( a, "Coefficient a" ) ) return;
  if ( ! NullCheck( y, "Dependent Variable y" ) ) return;
  if ( ! NullCheck( g, "Link Function g" ) ) return;

  unsigned int n = x.size(); // 独立変数の数
  unsigned int p = a.size(); // パラメータ数

  y.assign( n, 0 );

  for ( unsigned int i = 0 ; i < n ; ++i ) {
    if ( p != x[i].size() ) {
      cerr << "The size of x[" << i << "] is not equal to p ( = " << p << " ). Stop processing." << endl;
      return;
    }
    double d = 0;
    for ( unsigned int j = 0 ; j < p ; ++j )
      d += a[j] * x[i][j];
    if ( &g != 0 )
      y[i] = g.invf( d );
  }
}

/*
  Bin2Cnt : 二値データから計数データに変換する

  vector< vector<double> >& x : 独立変数ベクトル
  const vector<bool>& bin : 二値データ
  vector<unsigned int>& n : 計数データ(試行回数の総数)
  vector<unsigned int>& y : 計数データ(成功回数)
  bool verbose : 冗長モード(ON/OFF)
*/
bool Bin2Cnt( vector< vector<double> >& x, const vector<bool>& bin,
              vector<unsigned int>& n, vector<unsigned int>& y, bool verbose )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( bin, "Binary Data bin" ) ) return( false );
  if ( ! NullCheck( n, "Total Trial Count n" ) ) return( false );
  if ( ! NullCheck( y, "Success Count y" ) ) return( false );

  unsigned int totalCnt = bin.size(); // データ総数
  if ( totalCnt == 0 ) {
    cerr << "bin has no data." << endl;
    return( false );
  }

  if ( ! SizeCheck( x, "x", totalCnt, "bin" ) ) return( false );

  LessVector<double> lessVec; // ベクトル同士の比較関数オブジェクト
  /*
    mapCnt : 計数用Map
     first ... 独立変数ベクトル x
     second ... 総試行回数 n と 成功回数 y のペア
  */
  typedef map< vector<double>, pair<unsigned int, unsigned int>, LessVector<double> > MapCntType;
  MapCntType mapCnt( lessVec );

  // 計数処理
  for ( unsigned int i = 0 ; i < totalCnt ; ++i ) {
    MapCntType::iterator it = mapCnt.find( x[i] );
    if ( it == mapCnt.end() ) {
      mapCnt[x[i]] = pair<unsigned int, unsigned int>( 1, ( bin[i] ) ? 1 : 0 );
    } else {
      ++( ( it->second ).first );
      if ( bin[i] ) ++( ( it->second ).second );
    }
  }

  // x, n, y に計数結果を登録 (各データはクリアされる)
  x.clear();
  n.clear();
  y.clear();
  for ( MapCntType::iterator it = mapCnt.begin() ;
        it != mapCnt.end() ; ++it ) {
    x.push_back( it->first );
    n.push_back( ( it->second ).first );
    y.push_back( ( it->second ).second );
    if ( verbose ) {
      PrintVector( "x : ", x.back() );
      cout << "n = " << n.back() << " ; y = " << y.back() << endl;
    }
  }

  return( true );
}

/*
  HL_ShowPara : ホズマー・レメショウ・テストでのグループ化結果を出力する

  unsigned int grpNo : グループ番号
  const vector< vector<double> >& x : 独立変数ベクトル
  const vector<double>& pi : 確率の当てはめ値
  multimap<double,unsigned int>::const_iterator s, e : 要素番号保持データの範囲(eは最後の要素の次を指す)
  double so, fo : 成功・失敗回数(観測値)
  double se, fe : 成功・失敗回数(期待値)
*/
void HL_ShowPara( unsigned int grpNo,
                  const vector< vector<double> >& x, const vector<double>& pi,
                  multimap<double,unsigned int>::const_iterator s,
                  multimap<double,unsigned int>::const_iterator e,
                  double so, double fo, double se, double fe )
{
  // イテレータの型定義
  typedef multimap<double,unsigned int>::const_iterator MMapCit;

  cout << "*** Group No. = " << grpNo << " ***" << endl;
  MMapCit back; // 末尾のイテレータを保持する
  for ( MMapCit cit = s ; cit != e ; ++cit ) {
    unsigned int i = cit->second; // 要素番号
    std::ostringstream header;
    header << "x[" << i << "] = ";
    PrintVector( header.str(), x[i] );
    back = cit;
  }
  cout << pi[s->second] << " <= pi <= " << pi[back->second] << endl;
  cout << '\t' << "Obs." << '\t' << "Est." << endl;
  cout << "success\t" << so << '\t' << se << endl;
  cout << "failure\t" << fo << '\t' << fe << endl << endl;
}

/*
  HosmerLemeshowTest : ホズマー・レメショウ・テスト

  const vector< vector<double> >& x : 独立変数(表示のみに使用する)
  const vector<unsigned int>& n : 試行回数
  const vector<unsigned int>& y : y の実測値
  const vector<double>& pi : 確率の当てはめ値
  unsigned int grpCnt : グループ数
  bool verbose : 冗長モード(ON/OFF)

  戻り値 : Χ^2HLが得られた ... true ; 引数ミスなど ... false
*/
bool HosmerLemeshowTest( const vector< vector<double> >& x, const vector<unsigned int>& n,
                         const vector<unsigned int>& y, const vector<double>& pi,
                         unsigned int grpCnt, bool verbose )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( n, "Total Trial Count n" ) ) return( false );
  if ( ! NullCheck( y, "Success Count y" ) ) return( false );
  if ( ! NullCheck( pi, "Predictive Probability pi" ) ) return( false );

  // 要素数の取得とチェック
  unsigned int sz = x.size();
  if ( ! SizeCheck( n, "Total Trial Count n", sz, "Independent Variable x" ) ) return( false );
  if ( ! SizeCheck( y, "Success Count y", sz, "Independent Variable x" ) ) return( false );
  if ( ! SizeCheck( pi, "Predictive Probability pi", sz, "Independent Variable x" ) ) return( false );

  if ( grpCnt == 0 ) {
    cerr << "group count must be more than zero." << endl;
    return( false );
  }

  // 計数データの合計 totalCnt を求める
  unsigned int totalCnt = 0;
  for ( unsigned int i = 0 ; i < sz ; ++i ) {
    if ( n[i] < y[i] ) {
      cerr << "y[" << i << "] (" << y[i] << ") is more than n[" << i << "] (" << n[i] << ") ( o must be less than n )." << endl;
      return( false );
    }
    totalCnt += n[i];
  }
  if ( totalCnt == 0 ) {
    cerr << "total count must be more than zero." << endl;
    return( false );
  }
  unsigned int grpSz = totalCnt / grpCnt; // 1 グループあたりの要素数の目安

  // pi と要素番号をペアとする map の作成
  typedef multimap<double,unsigned int> MMap;
  MMap mmapIndex;
  for ( unsigned int i = 0 ; i < sz ; ++i )
    mmapIndex.insert( pair<double,unsigned int>( pi[i], i ) );

  double se = 0, so = 0, fe = 0, fo = 0; // グループごとの成功・失敗回数(実測・期待値)
  unsigned int cnt = 0;   // グループの総数
  unsigned int grpNo = 1; // グループ番号
  double chisqHL = 0;     // 求める X^2HL
  MMap::const_iterator cit = mmapIndex.begin();
  MMap::const_iterator preCit = cit; // グループの範囲の開始番号
  while ( cit != mmapIndex.end() ) {
    unsigned int i = cit->second; // 要素番号を取得
    so += y[i];
    fo += n[i] - y[i];
    se += (double)( n[i] ) * pi[i];
    fe += (double)( n[i] ) * ( 1.0 - pi[i] );
    cnt += n[i];
    ++cit;
    if ( ( ( cnt >= grpSz ) && ( grpNo < grpCnt ) ) || // グループの要素数が目安を超えた(かつ最後のグループではない)
         ( ( cnt > 0 ) && cit == mmapIndex.end() ) ) { // 最後のグループ(全データ集計完了)
      if ( verbose ) HL_ShowPara( grpNo, x, pi, preCit, cit, so, fo, se, fe );
      ++grpNo;
      preCit = cit;
      chisqHL += pow( so - se, 2 ) / se;
      chisqHL += pow( fo - fe, 2 ) / fe;
      so = fo = se = fe = 0;
      cnt = 0;
    }
  }
  cout << "Chi Square HL = " << chisqHL << endl;

  return( true );
}

サンプル・プログラムは大きく三つの関数に分かれていて、得られたデータが二値の場合、同じ独立変数ベクトル x を持つデータごとに総試行回数 n と成功回数 y を計算するための Bin2Cnt、求めた係数 a を使って成功確率の当てはめ値 pi を求めるための CalcPredictiveValue、最後に、ホズマー・レメショウ・テストを行うための HosmerLemeshowTest があります。また、ホズマー・レメショウ・テストを行うためには、独立変数ベクトル x・計数データ ( 総試行回数 n と成功回数 y )・連結関数 g を使って回帰係数 a を求めるための ScoringMethod_Binomial が必要になりますが、これは前節でのサンプル・プログラムで紹介しています。ホズマー・レメショウ・テスト自体は HosmerLemeshowTest 関数だけで実行することができて、その他の関数はサポートに用います。独立変数と二値データがある場合、ロジット関数を連結関数としてホズマー・レメショウ・テストを行うのであれば、以下のような流れで処理することになります。

vector< vector<double> > x; // 独立変数ベクトル
vector<bool> bin;           // 二値データ
vector<unsigned int> n;     // 計数データ(試行回数)
vector<unsigned int> y;     // 計数データ(成功回数)
LogitFunc logit;            // 連結関数(ロジット関数)

vector<double> a;  // 回帰計数
vector<double> pi; // 確率の当てはめ値

const unsigned int grpCnt = 10; // グループの数(ホズマー・レメショウ・テスト用)

  : (データ処理)

Bin2Cnt( x, bin, n, y );                     // 二値データから計数データへ
ScoringMethod_Binomial( x, n, y, a, logit ); // ロジスティック回帰
CalcPredictiveValue( x, a, pi, logit );      // 予測値を計算する
HosmerLemeshowTest( x, n, y, pi, grpCnt );   // ホズマー・レメショウ・テスト

当てはめ値を計算するための関数 CalcPredictiveValue は非常に単純なものですが、様々な場面で利用することができます。実際、前節にてサンプル・データから求められた当てはめ値はこれを利用して計算しています。ここで得られる当てはめ値は確率に対するものなので、n を掛けることで成功確率の当てはめ値が得られます。同様に、

と計算することで成功・失敗回数の実測値と予測値が得られるので、これをいくつかのグループごとに集計すれば、その結果を使って Χ2HL を求めることができます。グループの分け方としては、各グループの試行回数ができるだけ均等になるようにする方がよいので、引数として渡したグループ数 grpCnt で総試行回数 totalCnt を割って 1 グループあたりの試行回数の目安 grpSz を求め、あるグループに対して試行回数の和が grpSz を超えたら次のグループへ切り替えるようにします。


3) 名義ロジスティック回帰 (Nomial Logistic Regression)

前節において、各独立変数 x に対して二項分布 BN,π(y) を指数型分布族とし、連結関数 g(π) = xTα としてロジット関数などの累積分布関数を利用した一般化線形モデルを紹介しました。各従属変数は成功・失敗のいずれかを表す二値であり、同じグループに対する成功回数が二項分布になることを利用してモデル式を定義したわけですが、「成功」と「失敗」という呼び方は便宜的なもので、先の例では「生存率(死亡率)」を意味することになり、その他に性別(男性・女性)などが考えられます。しかし、この分類が三つ以上に分かれる場合を想定すると、二項分布を利用したモデル式は適用できなくなります。そこで、別の確率分布を利用することを検討してみます。

C 個のカテゴリからなる確率変数があり、それぞれのカテゴリの発生する確率を πk ( k = 1, 2, ... C ) とします。但し、Σk{1→C}( πk ) = 1 という制約条件が付加されます。これは、C 個のカテゴリの中で必ずどれか一つのみが発生するということを意味します。このとき、合計 N 回の試行に対して各カテゴリが yk 回発生するときの確率分布は「多項分布 (Multinomial Distribution)」に従い、y = ( y1, y2, ... yC )T, π = ( π1, π2, ... πC )T とすれば

PN,π( y ) = [ N! / Πk{1→C}( yk! ) ]Πk{1→C}( πkyk )

で表されます。

PN,π( y )=[ N! / Πk{1→C}( yk! ) ]exp( log ( Πk{1→C}( πkyk ) ) )
=[ N! / Πk{1→C}( yk! ) ]exp( Σk{1→C}( yklog( πk ) ) )

より、η( π ) = ( log( π1 ), log( π2 ), ... log( πC ) )、h( y ) = N! / Πk{1→C}( yk! ) とすれば上式は

PN,π( y ) = h( y )・exp( yTη( π ) )

となって、母数を複数とする指数型分布族であることを示すことができます。J = 2 のとき、多項分布は二項分布と等しくなり、このときは一母数指数型分布族となることは前に示した通りです。従って、C 個のカテゴリに対しては (C-1)-母数指数型分布族であることが予想できて、実際にそのようになります。これは、yπ の各要素の和に対する制約条件から、母数の数を一つ減らすことができるためです。

N 個の独立変数に対して πi = ( πi1, πi2, ... πiC )T ( i = 1, 2, ... N ) を定義したとき ( 但し Σk{1→C}( πik ) = 1 )、対数尤度 l は

l=Σi{1→N}( log( Pni,πi( yi ) ) )
=Σi{1→N}( yiTη( πi ) + log( h( yi ) ) )
=Σi{1→N}( Σk{1→C}( yiklog( πik ) ) + log ni! - Σk{1→C}( log yik! ) )

となります。但し、yi = ( yi1, yi2, ... yiC ) は i 番目の独立変数に対する各カテゴリの発生回数で、Σk{1→C}( yik ) = ni を満たします。飽和モデルに対しては、∂l / ∂πik = 0 を満たす πik を求めればよいので、

∂l / ∂πik=( ∂ / ∂πiki'{1→N}( log( Pni',πi'( yi' ) ) )
=( ∂ / ∂πikk'{1→C}( yik'log( πik' ) )

がゼロになるような πik を計算すればいいわけですが、πik' 全てが独立した変数ではなく制約条件 Σk'{1→C}( πik' ) = 1 を持つことから、例えば πi1 = 1 - Σk'{2→C}( πik' ) と表され、

∂l / ∂πik=( ∂ / ∂πik )[ yi1log( 1 - Σk'{2→C}( πik' ) ) + Σk'{2→C}( yik'log( πik' ) ) ]
=-yi1 / ( 1 - Σk'{2→C}( πik' ) ) + yik / πik
=yik / πik - yi1 / πi1 = 0

より

πik = yikπi1 / yi1

となります。また、Σk'{2→C}( πik' ) = 1 - πi1 から

1 - πi1 = Σk'{2→C}( yik'πi1 / yi1 ) = ( πi1 / yi1k'{2→C}( yik' )

であり、Σk'{2→C}( yik' ) = ni - yi1 なので

1 - πi1 = πi1( ni - yi1 ) / yi1

で、これを解くと πi1 = yi1 / ni となります。よって、

πik = yikπi1 / yi1 = yik( yi1 / ni ) / yi1 = yik / ni

という結果が得られます。この結果は k = 1 に対しても成り立っているので、k = 1 の場合だけ特別扱いする必要はありません。また、直感的にも理解しやすい結果となっています。

カテゴリが二つの場合、「ロジスティック・モデル」を利用した連結関数は次のような式でした。

log( πi / ( 1 - πi ) ) ≡ log( πi2 / πi1 ) = xiTα ≡ ξi

πi2 ≡ πi はいわば「成功」する確率で、πi1 ≡ 1 - πi が「失敗」する確率を意味します。これを三つ以上のカテゴリの場合に拡張して、ρik ≡ πik / πi1 が、連結関数 log( ρik ) = log( πik / πi1 ) によって p 個の独立変数の成分に対して線形で表されるとします。この時、

log( ρik ) = xiTαk = Σj{1→p}( αkjxij ) ≡ ξik

より

ρik = exp( xiTαk ) = eξik

と表されます。ρik ≡ πik / πi1 の和 Σk'{2→C}( ρik' ) は

Σk'{2→C}( ρik' )=Σk'{2→C}( πik' ) / πi1
=( 1 - πi1 ) / πi1

と計算できることから、

πi1 = 1 / [ 1 + Σk'{2→C}( ρik' ) ]

であり、

πik = ρikπi1 = ρik / [ 1 + Σk'{2→C}( ρik' ) ]

と表すことができます。逆に、πi1 = 1 - Σk'{2→C}( πik' ) より

ρik = πik / [ 1 - Σk'{2→C}( πik' ) ]

なので、πik は ρi2 から ρiC までの C - 1 個の変数を持つ関数で、ρik は πi2 から πiC までの C - 1 個の変数を持つ関数です。対数尤度 l は

l=Πi{1→N}( log( Pni,πi( yi ) ) )
=Σi{1→N}( Σk'{1→C}( yik'log( πik' ) ) + log ni! - Σk'{1→C}( log yik'! ) )
=Σi{1→N}( -yi1log( 1 + Σk'{2→C}( ρik' ) )
 + Σk'{2→C}( yik'log( ρik' / [ 1 + Σl{2→C}( ρil ) ] ) )
 + log ni! - Σk'{1→C}( log yik'! ) )
=Σi{1→N}( -yi1log( 1 + Σk'{2→C}( ρik' ) )
 + Σk'{2→C}( yik'log( ρik' ) - yik'log( 1 + Σl{2→C}( ρil ) ) )
 + log ni! - Σk'{1→C}( log yik'! ) )
=Σi{1→N}( Σk'{2→C}( yik'log( ρik' ) ) - nilog( 1 + Σk'{2→C}( ρik' ) )
 + log ni! - Σk'{1→C}( log yik'! ) )

なので、その αkj による偏微分 ukj = ∂l / ∂αkj

ukj=( ∂ / ∂αkji{1→N}( Σk'{2→C}( yik'log( ρik' ) ) - nilog( 1 + Σk'{2→C}( ρik' ) ) )
=Σi{1→N}( ( yik / ρik )( ∂ρik / ∂αkj )
 - [ ni / ( 1 + Σk'{2→C}( ρik' ) ) ]( ∂ρik / ∂αkj ) )
=Σi{1→N}( ( yik / ρik )( ∂ρik / ∂αkj )
 - { ni / [ 1 + ( 1 - πi1 ) / πi1 ] }( ∂ρik / ∂αkj ) )
=Σi{1→N}( ( yik / ρik - niπi1 )( ∂ρik / ∂αkj ) )
=Σi{1→N}( [ ( yik - niπik ) / ρik ]( ∂ρik / ∂αkj ) )

と表されます。ここで dρik / dξik = eξik = ρik、∂ξik / ∂αkj = xij より

∂ρik / ∂αkj=( dρik / dξik )( ∂ξik / ∂αkj )
=xijρik

なので、ukj

ukj = Σi{1→N}( ( yik - niπik )xij )

と求めることができます。

C = 2 のとき (つまり、二項分布の場合)、uj = ∂l / ∂αj

uj = Σi{1→N}( ( yi - niπi )xij / g'(πii( 1 - πi ) )

で求められるのでした。ロジスティック・モデルでは

g'(πi)=( d / dπi )log( πi / ( 1 - πi ) )
=[ ( 1 - πi ) / πi ]{ [ ( 1 - πi ) + πi ] / ( 1 - πi )2 }
=1 / πi( 1 - πi )

なので、分母は打ち消されて 1 になり、先ほど得られた式に一致します。つまり、二項分布を利用したロジスティック・モデルは、先ほど求めた多項分布によるモデル式に含まれることがわかります。


スコア法の漸化式は

α(m) = α(m-1) - H-1(m-1)u(α(m-1))

で表され、H の jr 行 jc 列目の要素は ∂ujr / ∂αjc = ∂2l / ∂αjr∂αjc となるのでした。ukj は ( C - 1 ) x p 個あるので、αu( α ) は要素数が ( C - 1 ) x p のベクトルで、H は行列数が ( C - 1 ) x p の対称行列です。∂ukj / ∂αk'j' を直接計算してみると、k ≠ k' ならば

∂ukj / ∂αk'j'=( ∂ / ∂αk'j'i{1→N}( ( yik - niπik )xij )
=Σi{1→N}( -nixij( ∂πik / ∂αk'j' ) )
=Σi{1→N}( -nixij・( ∂ / ∂αk'j' ){ ρik / [ 1 + Σl{2→C}( ρil ) ] } )
=Σi{1→N}( -nixijρik{ -1 / [ 1 + Σl{2→C}( ρil ) ]2 }( ∂ρik' / ∂αk'j' ) )
=Σi{1→N}( nixijρik[ 1 + ( 1 - πi1 ) / πi1 ]2xij'ρik' )
=Σi{1→N}( nixijxij'πi12ρikρik' )
=Σi{1→N}( nixijxij'πikπik' )

となりますが、k = k' の場合は

( ∂ / ∂αkj' ){ ρik / [ 1 + Σl{2→C}( ρil ) ] }
={ [ ( 1 + Σl{2→C}( ρil ) ) - ρik ] / [ 1 + Σl{2→C}( ρil ) ]2 }( ∂ρik / ∂αkj' )
=[ ( 1 / πi1 - πik / πi1 ) / ( 1 / πi1 )2 ]xij'ρik
=πi1( 1 - πik )xij'ρik

より

∂ukj / ∂αkj'=Σi{1→N}( -nixij[ πi1( 1 - πik )xij'ρik ] )
=Σi{1→N}( -nixijxij'πi1( 1 - πikik )
=Σi{1→N}( -nixijxij'πik( 1 - πik ) )

となります。多項分布において、niπikπik' は共分散を、また niπik( 1 - πik ) は分散を表します。

α = ( α21, α22, ... α2p, α31, ... , αkj, ... αCp )T

u( α ) = ( u21, u22, ... u2p, u31, ... , ukj, ... uCp )T

とし、H の行と列を、k が等しいものどうしを固める形で構成したとき ( つまり、H の ( k - 2 ) x p + j 行 ( k' - 2 ) x p + j' 列の要素を ∂ukj / ∂αk'j' = ∂2l / ∂ukj∂uk'j' としたとき )、H を ( C - 1 ) x ( C - 1 ) 個の p x p 部分行列からなる分割行列 ( 行列をいくつかの部分行列に区切って表現した行列 ) とみれば、対角成分にあたる部分行列は k = k' となるので後者の成分から構成され、それ以外の部分行列は前者の成分になります。x'j = ( x1j, x2j, ... xNj )T とし、

wikk'=niπikπik'[k ≠ k']
wikk=-niπik( 1 - πik )[k = k']

を対角要素とする対角行列を Wkk' とすれば、∂ukj / ∂αk'j' = x'jTWkk'x'j' で表されます。H を p x p 部分行列 ( ブロック ) に分けたとき、k 行 k' 列めのブロックは XTWkk'X で表され、分割行列は

H =
|XTW22X,XTW23X,...XTW2CX|
|XTW32X,XTW33X,...XTW3CX|
|::...:|
|XTWC2X,XTWC3X,...XTWCCX|
=
|XT,0,...0||W22,W23,...W2C||X,0,...0|
|0,XT,...0||W32,W33,...W3C||0,X,...0|
|::...:||::...:||::...:|
|0,0,...XT||WC2,WC3,...WCC||0,0,...X|

です。スコア法の漸化式は

(m) = (m-1) - u(α(m-1))

と変形できて、H の p x ( k - 2 ) + j 行目の行ベクトル hkjT

hkjT = ( ∂ukj / ∂α21, ∂ukj / ∂α22, ... ∂ukj / ∂α2p, ∂ukj / ∂α31, ... ∂ukj / ∂αk'j', ... ∂ukj / ∂αCp )

なので、

hkjTα=Σk'{2→C}( Σj'{1→p}( ( ∂ukj / ∂αk'j'k'j' ) )
=Σk'{2→C}( Σj'{1→p}( Σi{1→N}( wikk'xijxij'αk'j' ) ) )
=Σk'{2→C}( Σi{1→N}( wikk'xijΣj'{1→p}( xij'αk'j' ) ) )
=Σk'{2→C}( Σi{1→N}( wikk'xijξik' ) )

より、連立方程式の p x ( k - 2 ) + j 行目の式に対する右辺は

(右辺) = Σk'{2→C}( Σi{1→N}( wikk'xijξik' ) ) - ukj

となり、左辺は H を係数行列とする連立方程式なので、これを解くことを αkj が収束するまで繰り返します。

三つ以上のカテゴリに分類できる場合、多項分布とロジスティック・モデルを利用して確率の推定を行うことができることがこれまでの流れでわかりました。この手法は「多項ロジスティック回帰 (Multinomial Logistic Regression)」と呼ばれる中の一種で、分類されたカテゴリが「名義尺度 (Nominal Scale, Categorical Scale)」の場合に利用されることから「名義ロジスティック回帰 (Nominal logistic regression)」ともいいます。


名義ロジスティック回帰を行うためのサンプル・プログラムを以下に示します。

/*
  SizeCheck_Loop : コンテナクラスのサイズチェック(要素もコンテナクラスの場合)

  C& c : コンテナを要素とするコンテナクラス
  string arg1 : t の変数名(出力用)
  unsigned int sz : サイズ
  string arg2 : sz の大きさを持つ変数の名称(出力用)
*/
template<class C> bool SizeCheck_Loop( C& c, string arg1, unsigned int sz, string arg2 )
{
  std::ostringstream oss;
  unsigned int i = 0;
  for ( typename C::const_iterator it = c.begin() ; it != c.end() ; ++it ) {
    oss << arg1 << "[" << i << "]";
    if ( it->size() != sz ) {
      cerr << oss.str() << " size (" << it->size() << ") and " << arg2 << " size (" << sz << ") not matched." << endl;
      return( false );
    }
    ++i;
  }
  return( true );
}

/*
  MultinomialLogit_CalcDiagMatrix : 多項分布モデル用対角行列計算

  const vector< vector<double> >& pi : 確率の当てはめ値
  const vector<double>& ni : 各独立変数ごとの総数
  vector<double>& w : 対角行列 W_kk'
  unsigned int c : カテゴリ数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcDiagMatrix( const vector< vector<double> >& pi, const vector<double>& ni, vector<double>& w,
                                      unsigned int c, unsigned int n )
{
  vector<double>::iterator it = w.begin();
  for ( unsigned int rk = 0 ; rk < c - 1 ; ++rk ) {
    // 対角成分にあたる分割行列を先に処理する
    for ( unsigned int i = 0 ; i < n ; ++i )
      *it++ = -pi[i][rk + 1] * ni[i] * ( 1.0 - pi[i][rk + 1] );
    // その他の分割行列を処理
    for ( unsigned int ck = rk + 1 ; ck < c - 1 ; ++ck )
      for ( unsigned int i = 0 ; i < n ; ++i )
        *it++ = pi[i][rk + 1] * pi[i][ck + 1] * ni[i];
  }
}

/*
  MultinomialLogit_CalcCoef : 多項分布モデル用係数計算

  LinearEquationSystem<double>& s : 係数行列を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
  unsigned int rk, ck : 分割行列の行列番号
  unsigned int rj, cj : 分割行列内の計算対象要素の行列番号
*/
void MultinomialLogit_CalcCoef( LinearEquationSystem<double>& s,
                                const vector<double>& w, const vector< vector<double> >& x,
                                unsigned int c, unsigned int p, unsigned int n,
                                unsigned int rk, unsigned int rj, unsigned int ck, unsigned int cj )
{
  s[rk * p + rj][ck * p + cj] = 0;

  vector<double>::const_iterator cit =
    w.begin() + n * ( ( c * ( c - 1 ) - ( c - rk - 1 ) * ( c - rk ) ) / 2 + ck - rk );

  for ( unsigned int i = 0 ; i < n ; ++i )
    s[rk * p + rj][ck * p + cj] += x[i][rj] * x[i][cj] * *( cit + i );
}

/*
  Multinomial_CalcCoefMatrix : 多項分布モデル用係数行列の計算

  LinearEquationSystem<double>& s : 係数行列を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcCoefMatrix( LinearEquationSystem<double>& s,
                                      const vector<double>& w, const vector< vector<double> >& x,
                                      unsigned int c, unsigned int p, unsigned int n )
{
  for ( unsigned int rk = 0 ; rk < c - 1 ; ++rk ) {
    for ( unsigned int rj = 0 ; rj < p ; ++rj ) {
      // 分割行列もその要素も対角成分
      MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, rk, rj );
      // 分割行列が対角成分
      for ( unsigned int cj = rj + 1 ; cj < p ; ++cj ) {
        MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, rk, cj );
        s[rk * p + cj][rk * p + rj] = s[rk * p + rj][rk * p + cj];
      }
    }
    for ( unsigned int ck = rk + 1 ; ck < c - 1 ; ++ck ) {
      for ( unsigned int rj = 0 ; rj < p ; ++rj ) {
        // 分割行列は非対角で要素は対角成分
        MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, ck, rj );
        s[ck * p + rj][rk * p + rj] = s[rk * p + rj][ck * p + rj];
        // 分割行列もその要素も非対角
        for ( unsigned int cj = rj + 1 ; cj < p ; ++cj ) {
          MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, ck, cj );
          s[ck * p + cj][rk * p + rj] = s[ck * p + rj][rk * p + cj] =
            s[rk * p + cj][ck * p + rj] = s[rk * p + rj][ck * p + cj];
        }
      }
    }
  }
}

/*
  Multinomial_CalcUkj : 多項分布モデル用 u_kj の計算

  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  const vector<double>& ni : 各独立変数の総数 ( n )
  const vector< vector<double> >& pi : 確率の当てはめ値
  unsigned int k : カテゴリ番号
  unsigned int j : 独立変数の番号

  戻り値 : 計算した u_kj の値
*/
double MultinomialLogit_CalcUkj( const vector< vector<double> >& x, const vector< vector<double> >& y,
                                 const vector<double>& ni, const vector< vector<double> >& pi,
                                 unsigned int k, unsigned int j )
{
  unsigned int n = x.size();
  vector<double> d( n );

  for ( unsigned int i = 0 ; i < n ; ++i )
    d[i] = ( y[i][k + 1] - ni[i] * pi[i][k + 1] ) * x[i][j];

  return( sum( d ) );
}

/*
  MultinomialLogit_CalcRSide : 多項分布モデル用連立方程式の右辺の計算

  LinearEquationSystem<double>& s : 右辺を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  const vector<double>& ni : 各独立変数ごとの総数
  const vector< vector<double> >& pi : 確率の当てはめ値
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcRSide( LinearEquationSystem<double>& s, const vector<double>& w,
                                   const vector< vector<double> >& x, const vector< vector<double> >& y,
                                   const vector<double>& ni, const vector< vector<double> >& pi,
                                   unsigned int c, unsigned int p, unsigned int n )
{
  for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      s.ans( k * p + j ) = 0;
      for ( unsigned int k2 = 0 ; k2 < c - 1 ; ++k2 ) {
        unsigned int rk = ( k > k2 ) ? k2 : k; // 分割行列の行番号
        unsigned int ck = ( k > k2 ) ? k : k2; // 分割行列の列番号
        vector<double>::const_iterator cit =
          w.begin() + n * ( ( c * ( c - 1 ) - ( c - rk - 1 ) * ( c - rk ) ) / 2 + ck - rk );
        for ( unsigned int i = 0 ; i < n ; ++i )
          s.ans( k * p + j ) += x[i][j] * log( pi[i][k2 + 1] / pi[i][0] ) * *( cit + i );
      }
      s.ans( k * p + j ) -= MultinomialLogit_CalcUkj( x, y, ni, pi, k, j );
    }
  }
}

/*
  MultinomialLogistic : 名義ロジスティック回帰

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

  n は独立変数ベクトルの数
  p は独立変数ベクトルの要素数
  c はカテゴリの数

  x は p 個のパラメータのベクトルからなる n 個のベクトル
  y は c 個のパラメータのベクトルからなる n 個のベクトル
  a は p 個のパラメータのベクトルからなる c - 1 個のベクトル

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

  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( y, "Occurrence Count y" ) ) return( false );
  if ( ! NullCheck( a, "Coefficient a" ) ) return( false );

  unsigned int n = x.size(); // 独立変数ベクトル x_i の数
  if ( n == 0 ) {
    cerr << "x has no data." << endl;
    return( false );
  }
  if ( ! SizeCheck( y, "Occurrence Count y", n, "Independent Variable x" ) )
    return( false );

  unsigned int p = x[0].size(); // 独立変数ベクトルの要素数
  if ( ! SizeCheck_Loop( x, "Independent Variable x", p, "Independent Variable x[0]" ) )
    return( false );
  if ( p == 0 ) {
    cerr << "The size of x is zero." << endl;
    return( false );
  }
  unsigned int c = y[0].size(); // カテゴリ数
  if ( ! SizeCheck_Loop( y, "Occurrence Count y", c, "Occurrence Count y[0]" ) )
    return( false );
  if ( c == 0 ) {
    cerr << "The size of categories is zero." << endl;
    return( false );
  }

  cout << "N = " << n << " ; p = " << p << " ; c = " << c << endl << endl;

  if ( verbose ) {
    PrintMatrix( "x = ", x );
    cout << endl;
    PrintMatrix( "y = ", y );
    cout << endl;
  }

  vector<double> ni( n ); // 各独立変数の総数 ( n )
  for ( unsigned int i = 0 ; i < n ; ++i )
    ni[i] = sum( y[i] );

  vector< vector<double> > pi( n, vector<double>( c ) ); // 確率の当てはめ値 ( n x c )
  // piは y_ik / ni_i で初期化
  for ( unsigned int i = 0 ; i < n ; ++i )
    for ( unsigned int k = 0 ; k < c ; ++k )
      pi[i][k] = y[i][k] / ni[i];

  LinearEquationSystem<double> s( ( c - 1 ) * p ); // 連立方程式計算用インスタンス
  vector<double> w( n * c * ( c - 1 ) / 2 ); // 対角行列 W_kk'

  // 係数の初期化
  a.resize( c - 1 );
  for ( unsigned int k = 0 ; k < c - 1 ; ++k )
    a[k].assign( p, 0 );

  bool isMatched;   // 収束したか
  unsigned int cnt; // 計算回数
  for ( cnt = 0 ; cnt < maxCount ; ++cnt ) {

    if ( verbose ) {
      cout << "----- cnt = " << cnt << " -----" << endl << endl;
      PrintMatrix( "pi = ", pi );
      cout << endl;
    }

    // 対角行列の計算 W_kk' の計算
    MultinomialLogit_CalcDiagMatrix( pi, ni, w, c, n );
    // 係数行列の計算
    MultinomialLogit_CalcCoefMatrix( s, w, x, c, p, n );
    // 右辺の計算
    MultinomialLogit_CalcRSide( s, w, x, y, ni, pi, c, p, n );

    if ( verbose ) {
      cout << "Equation System :" << endl;
      s.print();
      cout << endl;
    }

    // 連立方程式の計算
    if ( ! GaussianElimination( s ) ) {
      cerr << "Failed to calculate coefficients." << endl;
      return( false );
    }

    // 各係数が収束しているかを確認する
    isMatched = true;
    for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
      for ( unsigned int j = 0 ; j < p ; ++j ) {
        if ( fabs( a[k][j] - s.ans( k * p + j ) ) >= threshold )
          isMatched = false;
        a[k][j] = s.ans( k * p + j );
      }
    }

    if ( verbose ) {
      for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
        std::ostringstream oss;
        oss << "Regression equation : y[" << k << "] = ";
        PrintEquation( oss.str(), a[k] );
      }
      cout << endl;
    }

    if ( isMatched ) break;

    // カテゴリごとの確率計算
    vector<double> rho( c - 1 );
    for ( unsigned int i = 0 ; i < n ; ++i ) {
      for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
        double xi = a[k][0] * x[i][0];
        for ( unsigned int j = 1 ; j < p ; ++j )
          xi += a[k][j] * x[i][j];
        rho[k] = exp( xi );
      }
      pi[i][0] = 1.0 / ( 1.0 + sum( rho ) );
      for ( unsigned int k = 1 ; k < c ; ++k )
        pi[i][k] = rho[k - 1] * pi[i][0];
    }
  }

  if ( cnt < maxCount ) {
    // 係数行列の再計算(=フィッシャー情報行列)
    MultinomialLogit_CalcDiagMatrix( pi, ni, w, c, n );
    MultinomialLogit_CalcCoefMatrix( s, w, x, c, p, n );
    LinearEquationSystem<double> inv( 0 ); // 連立方程式計算用インスタンスの逆行列
    Inverse( s, inv );

    cout << "Estimated regression equation" << endl << endl;
    for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
      std::ostringstream oss;
      oss << " y[" << k << "] = ";
      PrintEquation( oss.str(), a[k] );
      cout << "variance of a[" << k << "] = ( " << -inv[k * p][k * p];
      for ( unsigned int j = 1 ; j < p ; ++j )
        cout << ", " << -inv[k * p + j][k * p + j];
      cout << " )" << endl << endl;
    }

    cout << "Estimated probability" << endl;
    PrintMatrix( "pi = ", pi );
    cout << endl;
  } else {
    cout << "Failed to estimate regression coefficient" << endl << endl;
  }

  return( cnt < maxCount );
}

メインの関数は MultinomialLogistic で、非常に長いプログラムとなっています。最初の半分程度は引数のチェックや必要な変数の初期化などを行っている部分で、その後の for によるループ処理が反復計算を行っている最も重要な個所になります。
計算処理の最初で部分行列 Wkk' の要素を MultinomialLogit_CalcDiagMatrix を使って求めています。Wkk' は N 個の対角成分を持ち、H 内にある部分行列の個数は ( C - 1 )2 個ありますが、明らかに Wkk' = Wk'k であり、Wkk' のうち対角成分とその右上、または左下の部分行列だけに対して計算をすれば充分です。従って、対角成分を保持するのに必要な要素数は N x ( C - 1 ) x [ ( C - 1 ) + 1 ] / 2 = NC( C - 1 ) / 2 個でよく、計算量も半分程度減らすことができます。なお、計算は左上の部分行列から開始して右方向へ処理した後、その下の行の対角成分からまた右方向へ計算するという形で行い、最後に右下の対角成分にあたる部分行列を処理します。
次に、左辺の係数行列にあたる H の要素を MultinomialLogit_CalcCoefMatrix で計算します。XTWkk'X は明らかに対称行列であり、H 内での Wkk' の配置も対称であることから、XTWkk'X の対角成分より右上または左下を計算すれば、XTWkk'X の要素全てだけでなく、H 内で対称の位置にある部分行列 XTWk'kX の全要素も得られたことになります。係数行列の計算では、H の右上側の部分行列に対し、さらに右上の要素のみについて行い、その値を部分行列内の左下だけでなく、H 内での左下の部分行列に対してもコピーしています。これによって、計算量はほぼ 1 / 4 程度になります。
最後に MultinomialLogit_CalcRSide で右辺の計算を行い、連立方程式を解いたら、その解が収束しているかを確認します。なお、連立方程式の解法には「ガウスの消去法 (Gaussian elimination)」を利用しています。


ここでも文献にあったサンプル・データを使って処理を行なってみます。下記データは、車の安全性や装備の嗜好に関するドライバーへの聞き取り調査結果の中で、エアコンとパワーステアリングをどれだけ重視するかを示したものです。

表 3-1. エアコンとパワーステアリングの重要性
性別年齢(歳)反応
C/D B  A 
女性18-232612745
24-409211545
>405144160
男性18-234017865
24-4017151244
>408151841
10594 101300

どれだけ重視するか(反応)は 4 段階評価になっていて、「D : 重要でない」「C : あまり重要でない」「B : 重要」「A : 非常に重要」となっています。但し、表の中では「D : 重要でない」と「C : あまり重要でない」を合算した形にしてあります。このデータに対して、性別と年齢を独立変数ベクトルとし、各反応が独立変数ベクトルそれぞれの中で発生する確率を名義ロジスティック回帰で推定してみたいと思います。ここでは、「D : 重要でない」「C : あまり重要でない」を基準としてこの確率との比率を使います。また、独立変数ベクトルは以下のようにします。

xi1=1(男性)
=0(女性)
xi2=1(年齢 24-40)
=0(それ以外)
xi3=1(年齢 >40)
=0(それ以外)

このとき、正規方程式は

log( πik / πi1 ) = αk0 + αk1xi1 + αk2xi2 + αk3xi3 ( k = 2, 3 )

と表され、πi1 が、基準となる「D : 重要でない」「C : あまり重要でない」の発生確率、πi2, πi3 がそれ以外の「B : 重要」と「A : 非常に重要」に対する発生確率になります。また、18 から 23 歳までの女性は xi1 = xi2 = xi3 = 0 であり、これが他の性別・年齢に対する基準となります。

サンプル・プログラムを使って処理を行った結果を以下に示します (出力を冗長にするため引数 verbose を ON にしています)。

*** Multinomial Logistic Regression ***

N = 6 ; p = 4 ; c = 3

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

y = ( 26, 12, 7 )
    ( 9, 21, 15 )
    ( 5, 14, 41 )
    ( 40, 17, 8 )
    ( 17, 15, 12 )
    ( 8, 15, 18 )

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

pi = ( 0.577778, 0.266667, 0.155556 )
     ( 0.2, 0.466667, 0.333333 )
     ( 0.0833333, 0.233333, 0.683333 )
     ( 0.615385, 0.261538, 0.123077 )
     ( 0.386364, 0.340909, 0.272727 )
     ( 0.195122, 0.365854, 0.439024 )

Equation System :
(-62.6857)x0 + (-31.9524)x1 + (-21.0864)x2 + (-20.2455)x3 + (31.2019)x4 + (12.7686)x5 + (11.0909)x6 + (16.152)x7 = 14.0668
(-31.9524)x0 + (-31.9524)x1 + (-9.88636)x2 + (-9.5122)x3 + (12.7686)x4 + (12.7686)x5 + (4.09091)x6 + (6.58537)x7 = 6.5478
(-21.0864)x0 + (-9.88636)x1 + (-21.0864)x2 + (0)x3 + (11.0909)x4 + (4.09091)x5 + (11.0909)x6 + (0)x7 = -6.10144
(-20.2455)x0 + (-9.5122)x1 + (0)x2 + (-20.2455)x3 + (16.152)x4 + (6.58537)x5 + (0)x6 + (16.152)x7 = 8.43913
(31.2019)x0 + (12.7686)x1 + (11.0909)x2 + (16.152)x3 + (-54.7347)x4 + (-25.8402)x5 + (-18.7273)x6 + (-23.0809)x7 = -2.35318
(12.7686)x0 + (12.7686)x1 + (4.09091)x2 + (6.58537)x3 + (-25.8402)x4 + (-25.8402)x5 + (-8.72727)x6 + (-10.0976)x7 = 7.97945
(11.0909)x0 + (4.09091)x1 + (11.0909)x2 + (0)x3 + (-18.7273)x4 + (-8.72727)x5 + (-18.7273)x6 + (0)x7 = 3.35057
(16.152)x0 + (6.58537)x1 + (0)x2 + (16.152)x3 + (-23.0809)x4 + (-10.0976)x5 + (0)x6 + (-23.0809)x7 = -21.5174

Regression equation : y[0] = -0.594729x0 + -0.381807x1 + 1.13426x2 + 1.58749x3
Regression equation : y[1] = -1.03114x0 + -0.803674x1 + 1.46287x2 + 2.9008x3

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

pi = ( 0.524023, 0.28911, 0.186867 )
     ( 0.23501, 0.40309, 0.3619 )
     ( 0.098186, 0.264972, 0.636842 )
     ( 0.650933, 0.24515, 0.103917 )
     ( 0.349621, 0.409351, 0.241028 )
     ( 0.174038, 0.32061, 0.505352 )

Equation System :
(-63.3591)x0 + (-31.5974)x1 + (-21.4658)x2 + (-20.6163)x3 + (31.7604)x4 + (12.64)x5 + (10.9058)x6 + (16.7676)x7 = 14.4935
(-31.5974)x0 + (-31.5974)x1 + (-10.6384)x2 + (-8.93059)x3 + (12.64)x4 + (12.64)x5 + (4.34126)x6 + (6.64286)x7 = 7.13147
(-21.4658)x0 + (-10.6384)x1 + (-21.4658)x2 + (0)x3 + (10.9058)x4 + (4.34126)x5 + (10.9058)x6 + (0)x7 = -6.14966
(-20.6163)x0 + (-8.93059)x1 + (0)x2 + (-20.6163)x3 + (16.7676)x4 + (6.64286)x5 + (0)x6 + (16.7676)x7 = 8.99699
(31.7604)x0 + (12.64)x1 + (10.9058)x2 + (16.7676)x3 + (-55.4565)x4 + (-24.3506)x5 + (-18.4408)x6 + (-24.1253)x7 = -4.94817
(12.64)x0 + (12.64)x1 + (4.34126)x2 + (6.64286)x3 + (-24.3506)x4 + (-24.3506)x5 + (-8.04907)x6 + (-10.2488)x7 = 6.37968
(10.9058)x0 + (4.34126)x1 + (10.9058)x2 + (0)x3 + (-18.4408)x4 + (-8.04907)x5 + (-18.4408)x6 + (0)x7 = 2.62442
(16.7676)x0 + (6.64286)x1 + (0)x2 + (16.7676)x3 + (-24.1253)x4 + (-10.2488)x5 + (0)x6 + (-24.1253)x7 = -22.8294

Regression equation : y[0] = -0.590797x0 + -0.388125x1 + 1.12827x2 + 1.5877x3
Regression equation : y[1] = -1.03902x0 + -0.813x1 + 1.47805x2 + 2.91668x3

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

pi = ( 0.524195, 0.290344, 0.185461 )
     ( 0.234584, 0.40153, 0.363886 )
     ( 0.0975796, 0.264427, 0.637993 )
     ( 0.652471, 0.245143, 0.102386 )
     ( 0.350992, 0.407527, 0.241481 )
     ( 0.174276, 0.32035, 0.505374 )

Equation System :
(-63.3346)x0 + (-31.5786)x1 + (-21.4374)x2 + (-20.5971)x3 + (31.7196)x4 + (12.5992)x5 + (10.905)x6 + (16.7599)x7 = 14.5863
(-31.5786)x0 + (-31.5786)x1 + (-10.6237)x2 + (-8.92676)x3 + (12.5992)x4 + (12.5992)x5 + (4.33004)x6 + (6.63776)x7 = 7.17966
(-21.4374)x0 + (-10.6237)x1 + (-21.4374)x2 + (0)x3 + (10.905)x4 + (4.33004)x5 + (10.905)x6 + (0)x7 = -6.13127
(-20.5971)x0 + (-8.92676)x1 + (0)x2 + (-20.5971)x3 + (16.7599)x4 + (6.63776)x5 + (0)x6 + (16.7599)x7 = 9.00441
(31.7196)x0 + (12.5992)x1 + (10.905)x2 + (16.7599)x3 + (-55.3536)x4 + (-24.2819)x5 + (-18.4757)x6 + (-24.1063)x7 = -5.07926
(12.5992)x0 + (12.5992)x1 + (4.33004)x2 + (6.63776)x3 + (-24.2819)x4 + (-24.2819)x5 + (-8.05939)x6 + (-10.2488)x7 = 6.25696
(10.905)x0 + (4.33004)x1 + (10.905)x2 + (0)x3 + (-18.4757)x4 + (-8.05939)x5 + (-18.4757)x6 + (0)x7 = 2.62152
(16.7599)x0 + (6.63776)x1 + (0)x2 + (16.7599)x3 + (-24.1063)x4 + (-10.2488)x5 + (0)x6 + (-24.1063)x7 = -22.7995

Regression equation : y[0] = -0.590799x0 + -0.388128x1 + 1.12827x2 + 1.58771x3
Regression equation : y[1] = -1.03908x0 + -0.813018x1 + 1.47811x2 + 2.91675x3

Estimated regression equation

 y[0] = -0.590799x0 + -0.388128x1 + 1.12827x2 + 1.58771x3
 variance of a[0] = ( 0.0806426, 0.0903072, 0.116722, 0.162328 )

 y[1] = -1.03908x0 + -0.813018x1 + 1.47811x2 + 2.91675x3
 variance of a[1] = ( 0.109228, 0.103064, 0.160738, 0.178863 )

Estimated probability
pi = ( 0.524195, 0.290344, 0.185461 )
     ( 0.234584, 0.40153, 0.363886 )
     ( 0.0975796, 0.264427, 0.637993 )
     ( 0.652471, 0.245143, 0.102386 )
     ( 0.350992, 0.407527, 0.241481 )
     ( 0.174276, 0.32035, 0.505374 )

反復処理は三回で収束し、最後に推定結果が出力されます。推定確率は処理時に計算しているので、後で計算する必要はなくそのまま続けて出力されます。その結果は次のようになりました。

表 3-2. 名義ロジスティック回帰結果
性別年齢(歳)反応
123
推定確率当てはめ値実測値ピアソン残差推定確率当てはめ値実測値ピアソン残差推定確率当てはめ値実測値ピアソン残差
女性18-230.524223.59260.49650.290313.0712-0.29480.18558.357-0.465845
24-400.234610.569-0.47900.401518.07210.68960.363916.3715-0.339845
>400.09765.855-0.35330.264415.8714-0.46840.638038.28410.439760
男性18-230.652542.4140-0.37020.245115.93170.26700.10246.6680.521365
24-400.351015.44170.39600.407517.9315-0.69220.241510.63120.421844
>400.17437.1580.31970.320413.13150.51480.505420.7218-0.597641
二乗和0.9971.5971.3333.927

反応の部分は 1 が「D : 重要でない」または「C : あまり重要でない」、2 が「B : 重要」、そして 3 が「A : 非常に重要」になります。各反応に対し、サンプル・プログラムで得られた推定確率とそれによる当てはめ値 ( e )、実測値 ( o )、最後にピアソン残差 ( o - e ) / √e を示してあります。下端の「二乗和」行は、ピアソン残差の二乗和を計算した結果を示しています。右端の列にある「計」は各反応の当てはめ値または実測値の合計ですが、最も右下の数値はピアソン残差の二乗和の合計を表しており、これは「ピアソン・カイ二乗統計量 (Pearson Chi-squared Statistic)」と等しくなります。


当てはめ値から得られる最大対数尤度は

l( π | y )=Σi{1→N}( Σk{1→C}( yiklog( πik ) ) + log ni! - Σk{1→C}( log yik! ) )
=Σi{1→N}( Σk{1→C}( yiklog( πik ) ) ) + K

に推定確率 πik を代入することで求められます。但し、πik に依存しない項は定数項 K で表しています。最大モデルは πik = yik / ni となることは先述した通りです。また最小モデルは、各カテゴリに対する確率が i に依存しない ( すなわち πik = πk ) と考えればよいので、

l( π | y )=Σi{1→N}( Σk{1→C}( yiklog( πk ) ) ) + K
=Σk{1→C}( Σi{1→N}( yik )log( πk ) ) + K
Σk{1→C}( Yklog( πk ) ) + K

より ( 但し、Yk = Σi{1→N}( yik ) としています )、πk による偏微分 ∂l / ∂πk

∂l / ∂πk = Σk'{1→C}( Yk'( ∂ / ∂πk )log( πk' ) )

で求められます。ここでも π1 = 1 - Σk'{2→C}( πk' ) であることから

∂l / ∂πk = Yk / πk - Y1 / π1

となり、これがゼロのとき

πk = Ykπ1 / Y1

となります。Σk{1→C}( πk ) = 1 より

Σk{1→C}( Ykπ1 / Y1 ) = ( π1 / Y1k{1→C}( Yk ) = 1

となるので、

π1 = Y1 / Σk{1→C}( Yk )

と求めることができて、

πk = ( Yk / Y1 )[ Y1 / Σk{1→C}( Yk ) ] = Yk / Σk{1→C}( Yk )

という結果が得られます。つまり、各カテゴリごとの総数を全ての総数で割った値を使えばよいことになります。

名義ロジスティック・モデルを使って求めた推定確率から計算した対数尤度を lα とすると、サンプル・データにおける値は -25.37 になります。また、最大モデルに対する対数尤度 lmax = -23.40、最小モデルの場合 lmin = -64.29 となるので、対数尤度統計量 D と尤度比カイ二乗統計量 C は

D = 2 x ( lmax - lα ) = 3.94

C = 2 x ( lα - lmin ) = 77.84

と求められます。また、擬似 R2 値は

R2 = ( lmin - lα ) / lmin = 1 - lα / lmin

で求められますが、ここでの対数尤度 lmin, lα は定数項を除く必要があるので(*3-1)、その値は lmin = -329.27, lα = -290.35 となって、

R2 = 1 - 290.35 / 329.27 = 0.1182

という結果になります。飽和モデルでのパラメータ数は、性別 X 年齢(三種類) = 6 個の変数それぞれに対して基準カテゴリを除く二種類のカテゴリ(反応)があったので合計 12 個です。名義ロジスティック・モデルでは係数が 8 つあり、最小モデルでは 2 つになります。従って、D は自由度が 12 - 8 = 4、C は自由度が 8 - 2 = 6 の χ2-分布に漸近的に従い、それぞれの p 値は 0.4144, 9.966E-15 となるので、名義ロジスティック・モデルは最小モデルよりもデータに対して有意に適合しており、飽和モデルと比較してもうまくデータを表せていることが示されています。その反面、擬似 R2 値は 0.1182 と低く、このモデルはデータの全変動の 11.82% 程度しか説明できていないことになります。


確率の比 ρik = πik / πi1 は「オッズ (Odds)」と呼ばれます。ギャンブルなどでもよく使われる用語ですが、その場合は勝った時に支払われる量を意味し、日本の競馬などでは払戻金の倍率を表しています。C = 2 の場合、オッズは πi / ( 1 - πi ) となり、その対数はロジット関数そのものになります。
二つのオッズ ρik = πik / πi1 と ρlk = πlk / πl1 の比を「オッズ比(Odds Ratio)」といいます。ρik と ρlk の対数が線形関係

log( ρik ) = xiTαk

log( ρlk ) = xlTαk

にあるとき、オッズ比の対数は

log( ρik / ρlk ) = ( xi - xl )Tαk

となるので、もし αk がゼロベクトルに近い場合、オッズ比は 1 に近い値になるはずです。例えば、単純な例として xi = ( 1, 1 )Txl = ( 1, 0 )T の場合を考えると、

log( ρik ) = αk0 + αk1

log( ρlk ) = αk0

となるので、オッズ比の対数は

log( ρik / ρlk ) = log( ρik ) - log( ρlk ) = αk1

であり、αk1 = 0 ならば ρik / ρlk = 1 です。独立変数の第一成分が定数項、第二成分が要因の有無を表しているとすれば、オッズ比が 1 ならば要因の有無に対して発生確率が影響しないことを表します。また、オッズ比が 1 より大きければ αk1 > 0 であり要因によって確率は上がる傾向に、逆に 1 より小さければ確率は小さくなる傾向にあることも読み取ることができます。

サンプルデータにおける独立変数は次のような構成になっていました。

表 3-3. 独立変数ベクトルの構成
i定数項性別年齢 24-40年齢 >40
11000
21010
31001
41100
51110
61101

この表から、i = 4 の i = 1 に対するオッズ比は性別に対する影響を、また i = 2, 3 の i = 1 に対するオッズ比は年齢に対する影響をそれぞれ表すことがわかります。これらを計算した結果は次のようになります。

表 3-4. オッズ比
反応BA
性別0.6780.444
年齢 24-403.0904.384
年齢 >404.89218.48

性別による影響については、男性の方が女性よりもエアコンやパワーステアリングを重視しない傾向にあることがわかります。また、年齢に関しては、上昇するほど重要視する傾向が強くなることがはっきりと読み取れます。これらの値は定数項以外の係数 αkj を使って exp( αkj ) から求めることもできます。各係数の分散はフィッシャー情報行列から得ることができて、係数が正規分布に漸近的に従うことから 95% 信頼区間は 平均 ± 1.96 x 標準誤差(S.E.) を計算すれば求められ、次のような結果になります。

表 3-5. オッズ比の信頼区間
反応BA
性別[ 0.376, 1.222 ][ 0.236, 0.832 ]
年齢 24-40[ 1.582, 6.037 ][ 1.998, 9.621 ]
年齢 >40[ 2.221, 10.777 ][ 8.067, 42.34 ]

性別による差異については、B に対して区間が 1 を含んでおり、A はそうではないものの、最大値は 1 に近くなっています (実際、信頼度を 99% にすれば 1 を含むようになります)。よって、性別による差異についてははっきりとあると断言することはできません。


*3-1) 多項分布において、N = 1 で y が二値変数 ( 0 または 1 ) からなるベクトルとすれば、P1,π( y ) = Πk( πk ) であり、対数尤度の定数項はゼロになります。二項分布によるモデルにおいて、指数型分布族をベルヌーイ分布と考えた時、定数項が消滅するのと内容は同じです。ロジスティックモデルの場合、指数型分布族は二項分布や多項分布とはせずに、ベルヌーイ分布や、そのカテゴリが三つ以上になった場合の分布と考えるのが通常のようです。


<参考文献>
1. 「一般化線形モデル入門」 Annette J. Dobson 著 (共立出版)
2. Wikipedia

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system