説明

隣接カテゴリーロジットモデルにおけるパラメータ推定方法、システム及びプログラム

【課題】 隣接カテゴリーロジットモデルのパラメータ推定において、サンプル数が少ないときでも精度良く安定な推定法、システム及びプログラムを提供する。
【解決手段】 パラメータ推定部23は、入力装置1から与えられたデータに対して罰則付対数尤度関数l(θ)を最大にするパラメータ[外24]を推定し、推定されたパラメータ[外24]を事後確率計算部22へ送る。事後確率計算部は入力装置から与えられたデータに対して、パラメータ推定部から与えられたパラメータ[外24]を用いてyの事後確率を計算し、計算された事後確率を出力装置3へ送る。
[外24]

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、隣接カテゴリーロジットモデルにおけるパラメータ推定方法、システム及びプログラムに関し、特に、パラメータを適切に推定し、サンプル数が少ないときでも安定した推定値を与えることができる、隣接カテゴリーロジットモデルにおけるパラメータの推定方法、システム及びプログラムに関する。
【背景技術】
【0002】
順序カテゴリーデータ(ordered categorical data)は、順序をもったカテゴリーデータであって、医学や生物学あるいは社会学などにおいて非常に良く扱われるデータである。順序カテゴリーデータには、例えば、薬の臨床試験における症状を示す「悪化」、「不変」、「軽度改良」、「中等度改良」及び「著明改良」や、副作用を示す「グレード(Grade)0」、「グレード1」、「グレード2」、「グレード3」、「グレード4」及び「グレード5」などがある。これらは、まとめて多値反応データ(multiple response data)と呼ばれることもある。これに対して、血液型である「A型」、「B型」、「AB型」及び「O型」はデータ間に順序関係がないので、順序カテゴリーデータではない。
【0003】
順序カテゴリーデータに対する解析手法として、複数のカテゴリーを統合したり興味あるカテゴリーだけを抽出したりして、2カテゴリーのデータに作り直し、ロジスティック回帰分析が適用される。例えば、「グレード0」、「グレード1」及び「グレード2」を新しく[カテゴリー0]とし、「グレード3」、「グレード4」及び「グレード5」を新しく[カテゴリー1]とする。しかし、複数のカテゴリーを2カテゴリーに統合するので、本来データが持っている順序の情報が失われてしまう。
【0004】
一方、順序カテゴリーデータを2カテゴリーに統合することなく解析できる統計モデルとして、隣接カテゴリーロジット(adjacent categories logit:ACL)モデルがある(例えば、非特許文献1参照。)。ACLモデルは、順序カテゴリーデータに対する統計モデルとして医学や生物学あるいは社会学などにおいて一般的に広く使用されている。
【0005】
以下、ACLモデルについて説明する。
【0006】
順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を考える。このとき、応答変数Yの値が0となる確率を数10で定義する。
【数10】

【0007】
ここで、θ=(γ,β)、γ=(γ,・・・,γK−1)およびβ=(β,・・・,β)とする。
【0008】
いま、共変量xが与えられた状況下での応答変数Yの条件付確率を数11で定義する。
【数11】

【0009】
このACLモデルのパラメータθの推定法としては、例えば最尤法(以下、Ordinary likelihood method:OL法という)が用いられる。以下、図1及び図2を参照してOL法を用いたACLモデルのパラメータ推定手順、及びデータ解析ついて説明する。
【0010】
図1に従来の隣接カテゴリーロジットモデル解析システムを示す。図示の解析システムは、入力装置1、データ処理装置2、及び出力装置3を備えている。また、データ処理装置2は、通常の最尤法を用いるパラメータ推定部21と事後確率計算部22とを備えている。このデータ解析システムは、図2のフローチャートに示すように動作する。
【0011】
はじめに、入力装置1からデータ(y,x)を入力する(ステップB1)。次に、入力装置1に入力されたデータ(y,x)に対して、データ処理装置2の通常の最尤法を用いるパラメータ推定部21において、数12で示す対数尤度関数l(θ)を最大化するパラメータ[外1]を求める(ステップB2)。
【数12】

【0012】
[外1]

【0013】
次に、事後確率計算部22において、推定されたパラメータ[外2]及び入力装置1により入力されたデータ(y,x)を用いて数1及び数2よりyに対する事後確率A(y|x,θ)を計算する(ステップB3)。
【0014】
[外2]

【0015】
全てのデータyの事後確率を計算したかどうかを判定し、全てのyに対する事後確率を計算していれば終了し、計算していなければステップB3に戻る(ステップB4)。
【0016】
出力装置3においては、回帰係数であるβやパラメータγを出力したり、各データ(y,x)に対する事後確率A(y|x,θ)を出力したりする(ステップB5)。回帰係数βより目的変数yに対する説明変数x=(1,x,・・・,x)の寄与の度合いを判断することができる。
【0017】
ここで、OL法を用いたパラメータ推定には、不適切な推定値[外3]を与えたり、サンプル数が少ない場合には収束が不安定となる等の問題点があることが良く知られている。そこで、より安定してパラメータを推定する方法として、条件付最尤法(Conditional maximum likelihood method、以下CML法という)を用いる推定法が提案されている(例えば、非特許文献2及び3参照。)。しかしながら、CML法を用いたパラメータ推定は、ACLモデルの2カテゴリーの場合とみなせるロジスティック回帰モデルにおいてサンプルサイズが大きいときに計算量が膨大となり、実行が困難であることが知られている(非特許文献4及び5参照。)。
【0018】
[外3]

【0019】
【非特許文献1】Agresti, A. (1984) Analysis of Ordinal Categorical Data, New York: John Wiley.
【非特許文献2】Hirji, K.F (1992) Computing exact distribution for polytomous response data, Journal of American Statistical Society, 87, 487-492
【非特許文献3】Agresti, A. (1999) Modeling ordered categorical data: recent advantages and future challenges, Statistics in Medicine, 18, 2191-2207.
【非特許文献4】Bull, S. B., et al, (1997). Jackknife bias reduction for polychotomous logistic regression. Statistics in Medicine, 16, 545-560.
【非特許文献5】Mehta, C. R. et al, (2000) Efficient Monte Carlo methods for conditional logistic regression, Journal of the American Statistical Association, 95, 99-108.
【発明の開示】
【発明が解決しようとする課題】
【0020】
上述したようにOL法を用いたACLモデルのパラメータ推定は、精度が悪く、特にサンプル数が少ないときにパラメータの推定が不安定になるという問題点がある。
【0021】
また、CML法を用いたACLモデルのパラメータ推定は、計算量が非常に膨大になる場合があるという問題点がある。
【0022】
本発明の目的は、ACLモデルのパラメータθの推定において、従来の推定法よりも精度良く、またサンプル数が少ないときにも解の収束が安定するようなパラメータ推定方法、システム及びプログラムを提供することである。
【0023】
また、本発明の他の目的は、CML法を用いたパラメータ推定よりも計算量が少ないパラメータ推定方法、システム及びプログラムを提供することである。
【課題を解決するための手段】
【0024】
本発明は、第1の要旨によれば、入力データを入力するための入力装置と、該入力装置に入力された前記入力データを用いてパラメータを推定するデータ処理装置とを備えた隣接カテゴリーロジットモデル解析システムを用いたパラメータの推定方法において、前記データ処理装置が罰則付最尤法によりパラメータを推定することを特徴とする。
【0025】
具体的には、上記パラメータ推定方法において、順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率が数13で定義され、
【数13】

【0026】
共変量xが与えられた状況下での応答変数Yの条件付確率が数14で定義され、
【数14】

【0027】
前記罰則付最尤法が、数15により表される尤度関数を用いてパラメータθを推定する方法である、
【数15】

【0028】
ことを特徴とする。
【0029】
また、本発明は、第2の要旨によれば、入力データを入力するための入力装置と、該入力装置に入力された前記入力データを用いてパラメータを推定するデータ処理装置とを備えた隣接カテゴリーロジットモデル解析システムにおいて、前記データ処理装置が罰則付最尤法により前記パラメータを推定するパラメータ推定部を備えていることを特徴とする。
【0030】
具体的には、上記隣接カテゴリーロジットモデル解析システムにおいて、前記パラメータ推定部内に、順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率を数16で定義し、
【数16】

【0031】
かつ、共変量xが与えられた状況下での応答変数Yの条件付確率を数17で定義し、
【数17】

【0032】
当該パラメータ推定部が、前記罰則付最尤法として、数18により表される尤度関数を用いてパラメータθを推定する、
【数18】

【0033】
ことを特徴とする。
【0034】
さらに、本発明は、第3の要旨によれば、コンピュータに隣接カテゴリーロジットモデルについてパラメータを推定させるプログラムにおいて、罰則付最尤法により前記パラメータを推定するステップをコンピュータに実行させることを特徴とする。
【0035】
具体的には、上記プログラムにおいて、順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率を数19で定義し、
【数19】

【0036】
共変量xが与えられた状況下での応答変数Yの条件付確率を数20で定義し、
【数20】

【0037】
数21により表される罰則付尤度関数を用いてパラメータθを推定する、
【数21】

【0038】
ステップをコンピュータに実行させることを特徴とする。
【発明の効果】
【0039】
本発明の第1の効果は、隣接カテゴリーロジットモデルにおけるパラメータ推定において、対数尤度関数の代わりに罰則付尤度関数を用いてパラメータθの推定を行うようにしたことで、パラメータθを精度良く、またサンプル数が少ないときにも安定して推定することができることである。
【0040】
本発明の第2の効果は、隣接カテゴリーロジットモデルにおけるパラメータ推定において、罰則付最尤法を用いていることにより、通常の最尤法と同じ計算量でパラメータを推定することができること、つまり、条件付最尤法よりも少ない計算量でパラメータを推定することができることである。
【発明を実施するための最良の形態】
【0041】
以下、本発明の実施の形態について図面を参照して詳細に説明する。
【0042】
図3を参照すると、本発明の一実施の形態に係る隣接カテゴリーロジットモデル解析システムは、キーボード等の入力装置1と、プログラム制御により動作するデータ処理装置2と、ディスプレイ装置や印刷装置等の出力装置3とを含む。
【0043】
データ処理装置2は、罰則付最尤法を用いるパラメータ推定部23と、事後確率計算部22とを備えている。データ処理装置2は、例えば、コンピュータプログラムとそれを実行するコンピュータによって実現される。コンピュータプログラムは、フレキシブルディスクやハードディスク等の磁気ディスク、CD−ROMやMO等の光ディスク、あるいは半導体メモリなど、コンピュータ読み取り可能な記録媒体により外部から提供されてもよいし、コンピュータに内蔵する記憶装置に記憶されていてもよい。あるいは、コンピュータにネットワーク接続された外部記憶装置や他のコンピュータ内の記憶装置に記憶されていてもよい。
【0044】
このシステムの解析対象である隣接カテゴリーロジットモデルは、以下のように定義される。
【0045】
順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)が与えられるとき、応答変数Yの値が0となる確率を数22で定義する。
【数22】

【0046】
ここで、θ=(γ,β)、γ=(γ,・・・,γK−1)およびβ=(β,・・・,β)である。
【0047】
いま、共変量xが与えられた状況下での応答変数Yの条件付確率を数23で定義する。
【数23】

【0048】
本実施の形態に係る解析システムは、数24に示す罰則付尤度関数を用いて、ACLモデルに対するパラメータθの推定を行う。なお、このような罰則付尤度関数を用いた最尤法を、以下、PL法(Penalized Likelihood method)という。
【数24】

【0049】
ここで、ξは罰則パラメータ(非負の定数)である。なお、ξ=0のとき、数24は数12に等しい。なお、基準化した共変数を用いる場合、上記罰則項の代わりに数24−1を使用すると、解の偏りは大きくなるが、収束性がさらによくなる。
【数24−1】

以下、図3の解析システムの動作を図4を参照して説明する。
【0050】
まず、入力装置1から順序カテゴリーデータ及び共変量が入力されると(ステップC1)、罰則付最尤法を用いるパラメータ推定部23は、罰則パラメータξを初期化する(ステップC2)。
【0051】
次に、パラメータ推定部23は、入力装置1から与えられたデータに対して、上記数24で示す罰則付対数尤度関数l(θ)を最大にするパラメータ[外4]を推定する(ステップC3)。推定されたパラメータ[外4]は、事後確率計算部22及び出力装置3へ送られる。
【0052】
[外4]

事後確率計算部22は、入力装置1から与えられたデータに対して、罰則付最尤法を用いるパラメータ推定部23から与えられたパラメータ[外5]を用いて、yの事後確率を計算する(ステップC4)。計算された事後確率は出力装置3へ送られる。
【0053】
[外5]

全てのデータyの事後確率を計算したかどうかを判定し、全てのyに対する事後確率を計算していれば終了し、計算していなければステップC4に戻る(ステップC5)。
【0054】
出力装置3においては、回帰係数であるβやパラメータγを出力したり、各データ(y,x)に対する事後確率A(y|x,θ)を出力したりする(ステップC6)。
【0055】
以上のようにして得られた回帰係数βより、目的変数yに対する説明変数x=(1,x,・・・,x)の寄与の度合いを判断することができる。ここで、入力データとして遺伝子データ(説明変数x)と臨床データ(目的変数y)とが与えられたとする。遺伝子データは、例えば、マイクロアレイなどによって得られる遺伝子発現データである。また、臨床データは、例えば、副作用や薬効などの順序カテゴリーデータである。この場合、出力装置3から出力される回帰係数の推定値[外6]より、副作用や薬効などの臨床データに対する各遺伝子の寄与の度合いを判断することができる。つまり、[外6]が正となる遺伝子に対しては、副作用や薬効などの効果を増加させる働きがあることが予想できる。また、[外6]の絶対値の大きさの大小により、各遺伝子の副作用や薬効などに対する寄与の度合いを比較することができる。
【0056】
[外6]

【実施例1】
【0057】
次に、本発明の実施例を、シミュレーションの結果を参照して具体的に説明する。ここでは、OL法とPL法とをそれぞれ用いたパラメータ推定の性能を比較する。なお、本実施例におけるシミュレーションにおいてはx〜N(0,I)とする。ただし、Iは数25で表され、N(0,I)は3次元の標準正規分布とする。
【数25】

【0058】
(シミュレーション1)
シミュレーション1においては、数26で表される確率に従って順序カテゴリーデータy∈{0,1,2}が生成される。
【数26】

【0059】
ここで、γ=1,β=(0,0,1,−3)とする。また、A(y|x,(γ,β))は数27(=数23の特別な場合)で表されるものとする。本モデルはACLモデルに対応している。
【数27】

【0060】
(シミュレーション2)
シミュレーション2においては、数28に従って順序カテゴリーデータy∈{0,1,2}が生成される。
【数28】

【0061】
シミュレーション1及び2のそれぞれにおいて、サンプル数(N)を20,30,50,100と変化させながら1000回ずつ繰り返し解析を行う。ただし、各カテゴリーに属するサンプル数が少なくとも5つあるデータセットのみをシミュレーションに使用する。また、PL法における罰則パラメータξをξ=1.0とする。
【0062】
表1は不適切な推定値が得られたり、解が収束しなかった頻度を示した表である。
【表1】

【0063】
ここで不適切な推定値とは、例えば数29で示されるように、推定されたパラメータの下限が−10を下回り、又は上限が10を超えることをいう。
【数29】

【0064】
表1から、サンプル数及びシミュレーションタイプに関係なくPL法がOL法よりも安定に適切な推定値を与えることが分かる。特にシミュレーション2においては、サンプル数が少ない(n=20,30)場合においては、OL法では不適切な推定値を与える頻度が30%を超えるにもかかわらず、PL法では1%程度の頻度である。
【0065】
図5乃至図8にシミュレーション1(N=50)において推定された回帰係数[外7]及び[外8]の分布を示す。なお、β=0,β=0,β=1,β=−3,γ=1の直線はそれぞれのパラメータの真の値を示している。図5乃至図8を見て分かるように、罰則の項を尤度関数に導入することによって、PL法によって推定されたパラメータ[外7]と[外8]の相関がOL法によって推定されたパラメータ[外7]と[外8]の相関よりも小さくなっていることがわかる。即ち、PL法によって推定されたパラメータの方が、OL法によって推定されたパラメータよりも安定している。
【0066】
[外7]

【0067】
[外8]

【0068】
図9乃至図13にシミュレーション1(N=20,30,50,100)において推定された回帰係数[外9]及び[外10]の分布を箱型図で示す。図9乃至図13から分かるように、サンプル数(N)が大きくなるにつれてOL法、PL法ともに推定値の分布が安定することがわかる。また、サンプル数が少ない場合、OL法よりもPL法の方が推定値の分布が安定していることが分かる。特に図14のN=20,30の場合を比較すると、サンプル数が少ない場合においてもPL法の方がより高精度にβを推定できることが分かる。
【0069】
[外9]

【0070】
[外10]

【図面の簡単な説明】
【0071】
【図1】従来の隣接カテゴリーロジットモデル解析システムの構成を示すブロック図である。
【図2】図1の解析システムの動作を説明するための流れ図である。
【図3】本発明の一実施の形態に係る隣接カテゴリーロジットモデル解析システムの構成を示すブロック図である。
【図4】図3の解析システムの動作を説明するための流れ図である。
【図5】本発明の実施例(シミュレーション1:N=50)により求められた[外11]と[外12]の相関を表す相関図である。
【0072】
[外11]

[外12]

【図6】本発明の実施例(シミュレーション1:N=50)により求められた[外13]と[外14]の相関を表す相関図である。
【0073】
[外13]

[外14]

【図7】本発明の実施例(シミュレーション1:N=50)により求められた[外15]と[外16]の相関を表す相関図である。
【0074】
[外15]

[外16]

【図8】本発明の実施例(シミュレーション1:N=50)により求められた[外17]と[外18]の相関を表す相関図である。
【0075】
[外17]

[外18]

【図9】本発明の実施例(シミュレーション1:N=20,30,50,100)により求められた[外19]の分布を示す図である。
【0076】
[外19]

【図10】本発明の実施例(シミュレーション1:N=20,30,50,100)により求められた[外20]の分布を示す図である。
【0077】
[外20]

【図11】本発明の実施例(シミュレーション1:N=20,30,50,100)により求められた[外21]の分布を示す図である。
【0078】
[外21]

【図12】本発明の実施例(シミュレーション1:N=20,30,50,100)により求められた[外22]の分布を示す図である。
【0079】
[外22]

【図13】本発明の実施例(シミュレーション1:N=20,30,50,100)により求められた[外23]の分布を示す図である。
【0080】
[外23]

【符号の説明】
【0081】
1 入力装置
2 データ処理装置
3 出力装置
21 通常の最尤法を用いるパラメータ推定部
22 事後確率計算部
23 罰則付最尤法を用いるパラメータ推定部

【特許請求の範囲】
【請求項1】
入力データを入力するための入力装置と、該入力装置に入力された前記入力データを用いてパラメータを推定するデータ処理装置とを備えた隣接カテゴリーロジットモデル解析システムを用いたパラメータの推定方法において、
前記データ処理装置が罰則付最尤法によりパラメータを推定することを特徴とする隣接カテゴリーロジットモデルのパラメータ推定方法。
【請求項2】
請求項1に記載のパラメータ推定方法において、
順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率が数1で定義され、
【数1】

共変量xが与えられた状況下での応答変数Yの条件付確率が数2で定義され、
【数2】

前記罰則付最尤法が、数3により表される尤度関数を用いてパラメータθを推定する方法である、
【数3】

ことを特徴とするパラメータ推定方法。
【請求項3】
入力データを入力するための入力装置と、該入力装置に入力された前記入力データを用いてパラメータを推定するデータ処理装置とを備えた隣接カテゴリーロジットモデル解析システムにおいて、
前記データ処理装置が罰則付最尤法により前記パラメータを推定するパラメータ推定部を備えていることを特徴とする隣接カテゴリーロジットモデル解析システム。
【請求項4】
請求項3に記載の隣接カテゴリーロジットモデル解析システムにおいて、
前記パラメータ推定部内に、
順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率を数4で定義し、
【数4】

かつ、共変量xが与えられた状況下での応答変数Yの条件付確率を数5で定義し、
【数5】

当該パラメータ推定部が、前記罰則付最尤法として、数6により表される尤度関数を用いてパラメータθを推定する、
【数6】

ことを特徴とする隣接カテゴリーロジットモデル解析システム。
【請求項5】
請求項3又は4に記載された隣接カテゴリーロジットモデル解析システムにおいて、
前記パラメータ推定部が推定したパラメータを利用して、前記入力データについて、事後確率計算を行う事後確率計算部をさらに備えていることを特徴とする隣接カテゴリーロジットモデル解析システム。
【請求項6】
請求項3又は4に記載された隣接カテゴリーロジットモデル解析システムにおいて、
前記パラメータ推定部が推定したパラメータを外部出力する出力装置をさらに備えていることを特徴とする隣接カテゴリーロジットモデル解析システム。
【請求項7】
請求項5に記載された隣接カテゴリーロジットモデル解析システムにおいて、
前記パラメータ推定部が推定したパラメータ及び前記事後確率計算部が計算した事後確率を外部出力する出力装置をさらに備えていることを特徴とする隣接カテゴリーロジットモデル解析システム。
【請求項8】
コンピュータに隣接カテゴリーロジットモデルについてパラメータを推定させるプログラムにおいて、
罰則付最尤法により前記パラメータを推定するステップをコンピュータに実行させるプログラム。
【請求項9】
請求項8に記載のプログラムにおいて、
順序カテゴリーデータとしてY∈{0,・・・,K}、共変量としてx=(1,x,・・・,x)を与えたときの応答変数Yの値が0となる確率を数7で定義し、
【数7】

共変量xが与えられた状況下での応答変数Yの条件付確率を数8で定義し、
【数8】

数9により表される罰則付尤度関数を用いてパラメータθを推定する、
【数9】

ステップをコンピュータに実行させることを特徴とするプログラム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11】
image rotate

【図12】
image rotate

【図13】
image rotate


【公開番号】特開2006−24040(P2006−24040A)
【公開日】平成18年1月26日(2006.1.26)
【国際特許分類】
【出願番号】特願2004−202377(P2004−202377)
【出願日】平成16年7月8日(2004.7.8)
【新規性喪失の例外の表示】特許法第30条第1項適用申請有り 平成16年5月29日 九州大学主催の「九州大学21世紀COEプログラム柳川堯先生退官記念シンポジウム バイオ統計学最近の展開」において文書をもって発表
【出願人】(000004237)日本電気株式会社 (19,353)
【出願人】(504136568)国立大学法人広島大学 (924)
【出願人】(504264160)
【Fターム(参考)】