JP2007513391A - システムの複数の構成要素のうちのサブセットを識別する方法 - Google Patents

Abstract

システムからの少なくとも１つのトレーニングサンプルを使用し、システムから取得されたデータに基づき上記システムの構成要素のうちのサブセットを識別する。本方法は、システムの構成要素とその重み係数との一次結合を取得するステップを含み、重み係数は、既知の特徴を有する少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有する。本方法は、構成要素の一次結合を条件とする、既知の特徴の確率分布モデルを取得するステップと、構成要素の一次結合に係る重み係数の事前分布を取得するステップとを含み、事前分布は、ゼロに近い高い確率密度を有しかつジェフリーズの超事前分布ではない超事前分布を含む。本方法は、事前分布とモデルとを組み合わせて事後分布を生成するステップと、事後分布を最大化する複数の重み係数にてなるセットに基づいて構成要素のうちのサブセットを識別するステップとを含む。

Description

本発明は、システムの複数のサンプルから生成されるデータから、システムの複数の構成要素（又は成分）を識別するための方法及び装置に関し、ここで上記構成要素はシステム内のサンプルの特徴を予測することができるものである。また特に、本発明は、生物学的方法によって生成されるデータから、生物学的システムの複数の構成要素を識別するための方法及び装置に関するが、ただしそれ以外を除外するものではない。ここで、上記構成要素は、生物学的システムに適用されるサンプルに関連づけられる、関心対象の特徴を予測することができるものである。

固有の１つ又は複数の特徴によって分類可能なシステムが、いくつも存在している。本明細書を通じて使用される「システム」という用語は、それからデータ（例えば統計データ）を取得することのできる全てのタイプのシステムを包含するものとされる。このようなシステムの例には、化学システム、財務システム及び地質学的システムがある。システムから取得されるデータを利用してシステムからの複数のサンプルに係る特定の特徴を識別できること、例えば、財務システムの分析を支援して信用のあるグループと信用リスクのあるグループとを識別できることが望ましい。システムから取得されるデータは比較的大量であることが多く、よってそのデータからシステムの複数の構成要素を識別することが望ましい。上記構成要素は、システムからの複数のサンプルに係る特定の特徴を予測するものである。しかしながら、データの量が比較的多い場合には、処理すべき大量のデータが存在することになるので、構成要素を識別することは困難になる可能性がある。その大量のデータのうちのほとんどは、データの取得元である特定のサンプルの特徴を全く示していないか、ほとんど示していない可能性もある。さらに、試験サンプルデータがトレーニングサンプルデータに関して高い可変度を有する場合には、トレーニングサンプルを用いて識別される構成要素は、試験サンプルデータに関する特徴の識別に際して有効でない場合が多い。個々のソースからデータを収集する際の条件を制御することはしばしば困難であるので、例えば、異なる多くのソースからデータが取得されるような状況ではよくこうした事態になる。

これらの問題点が特に顕著であるシステムのタイプの一例は、構成要素が例えば特定の遺伝子又はタンパク質を含む可能性のある生物学的システムである。バイオテクノロジーの最近の進歩は、大規模なシステムスクリーニング及びサンプル分析のための生物学的方法の開発をもたらしている。このような方法には、例えば、ＤＮＡ又はＲＮＡを用いるマイクロアレイ解析と、プロテオミクス（proteomics）解析と、プロテオミクスにおける電気泳動ゲル解析と、高スループットスクリーニング技術とが含まれる。これらのタイプの方法は、検査される各サンプルについて３０，０００個以上もの構成要素を有する可能性のあるデータを生成させる場合が多い。

生物学的システムからのサンプルにおいて、例えば「病気を持っている」及び「病気を持たない」等のグループに分類するというように、関心対象の特徴を識別することができれば、かなり望ましい。これらの生物学的方法の多くは、生物学的システムにおけるサンプルの特徴を予測する診断ツールとして、例えば、組織もしくは体液をスクリーニングして病気を識別する診断ツールとして有用であり、あるいは、例えば医薬化合物の効力を決定するためのツールとしても有用であろう。

今日まで、このようなアプリケーションにおけるバイオテクノロジーアレイ等の生物学的方法の使用は、これらのタイプの方法から生成されるデータが大量であることと、有意義な結果をもたらすデータの効率的なスクリーニング方法がないこととに起因して限定的であった。その結果、既存の方法を用いる生物学的データの解析は時間がかかり、間違った結果になりやすく、また、データから有意義な結果を得ようとすれば大量のコンピュータメモリが必要である。このことは、高速かつ正確なスクリーニングが要求される大規模スクリーニングのシナリオにおいて問題となる。

従って、特に生物学的データを解析する方法、そしてより一般的には、システムからのサンプルに関して関心対象の特徴を予測するためにシステムからのデータを解析する改良された方法を有することが望ましい。

本発明の第１の態様によれば、システムからの少なくとも１つのトレーニングサンプルを使用し、システムから取得されるデータに基づいて上記システムの複数の構成要素のうちのサブセットを識別する方法が提供されていて、上記方法は、
上記システムの複数の構成要素と、上記複数の構成要素の一次結合に係る複数の重み係数との一次結合を取得するステップを含み、上記重み係数は、上記少なくとも１つのトレーニングサンプルを用いて上記システムから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の特徴を有し、
上記既知の特徴の確率分布のモデルを取得するステップを含み、上記モデルは上記複数の構成要素の一次結合を条件とし、
上記複数の構成要素の一次結合に係る重み係数の事前分布を取得するステップを含み、上記事前分布は、ゼロに近い高い確率密度を有する超事前分布（hyperprior）を含み、上記超事前分布はジェフリーズ（Jeffreys）の超事前分布ではないようなものであり、
上記事前分布と上記モデルとを組み合わせて事後分布を生成するステップと、
上記事後分布を最大化する複数の重み係数にてなるセットに基づいて上記複数の構成要素のうちのサブセットを識別するステップとを含む。

本方法は、既知の特徴を有する複数のトレーニングサンプルを利用して、あるトレーニングサンプルの特徴を予測することのできる複数の構成要素のうちのサブセットを識別する。続いて、複数の構成要素のうちのサブセットについての知識は、試験のために、例えば臨床試験のために使用可能であり、それにより、組織サンプルが悪性であるか良性であるか、もしくは腫瘍の重さはどの程度かといった特徴を予測することができ、又は特定の状態を有する患者の推定生存期間を決めることができる。

本明細書を通じて使用される「特徴」という用語は、あるサンプルに関連づけられる任意の応答又は識別可能な特性又は性質を示す。例えば特徴は、特定のサンプルに関するイベントまでの特定の時間である場合もあれば、サンプルのサイズもしくは量であってもよく、又はサンプルを分類するために使用可能なクラスもしくはグループであってもよい。

好適には、上記一次結合を取得するステップは、ベイズの統計的方法を用いて複数の重み係数を推定するステップを含む。

好適には、本方法はさらに、大部分の構成要素は、複数の構成要素のうちのサブセットの一部を形成する構成要素にはなりそうにない、という先験的仮定（アプリオリな仮定）を立てるステップを含む。

先験的仮定は、システムから取得される構成要素が大量に存在する場合において特に適用される。先験的仮定は、本質的には、大部分の重み係数はゼロになるであろう、というものである。先験的仮定に留意して、重み係数が、観測されるデータに与えられる重み係数の事後確率が最大化されるようなものであるように、モデルは構成される。予め決められたしきい値を下回る重み係数を有する構成要素（先験的仮定によるものの大部分はこれになる）は、無視される。このプロセスは、正しい診断構成要素が識別されるまで反復される。従って本方法は、主として、結果的に大部分の構成要素の迅速な除去をもたらす先験的仮定に起因して、高速になる可能性を有する。

好適には、超事前分布は、ゼロに近い事前分布が変更されることを可能にする、１つ又は複数の調節可能なパラメータを含む。

システムの大部分の特徴は典型的には所定の確率分布を提示し、上記特徴の確率分布は、複数のトレーニングサンプルから生成されるデータに基づく複数の統計モデルを用いてモデル化されることが可能である。本発明は、関心対象の特徴又は関心対象の一連の特徴について確率分布をモデル化する統計モデルを使用する。従って、特定の確率分布を有する関心対象の特徴について、その分布をモデル化する適切なモデルが画成される。

好適には、本方法は、少なくとも１つのトレーニングサンプルから取得されるデータに基づいて確率分布を提供する尤度関数の形式の数式を含む。

好適には、尤度関数は、何らかの確率分布を記述するための以前に記述されたモデルに基づく。

好適には、上記モデルを取得するステップは、多項又は２項ロジスティック回帰と、一般化線形モデルと、コックス（Cox）の比例ハザードモデルと、加速度故障モデル（accelerated failure model）と、パラメトリック生存モデルとを含むグループから上記モデルを選択するステップを含む。

第１の実施形態では、尤度関数は多項又は２項ロジスティック回帰に基づく。多項又は２項ロジスティック回帰は、好適には、多項又は２項分布を有する特徴をモデル化する。２項分布は、オン／オフ状態等の２つの可能なクラス又はグループを有する統計分布である。このようなグループの例には、死亡／生存、改良／未改良、抑制（depressed）／非抑制（not depressed）が含まれる。多項分布は２項分布の一般化であり、複数のサンプルの各々について複数のクラス又はグループが可能なものであり、又は言い替えれば、１つのサンプルが複数のクラス又はグループのうちの１つに分類されることが可能なものである。従って、多項又は２項ロジスティック回帰に基づいて尤度関数を定義すれば、あるサンプルを複数の予め定義されたグループ又はクラスのうちの１つに分類することができる複数の構成要素にてなるサブセットを識別することが可能である。これを行うために、複数のトレーニングサンプルは、トレーニングサンプルの予め決められた特徴に基づいて複数のサンプルグループ（又は「クラス」）にグループ化されるが、ここで、各サンプルグループの要素は共通の特徴を有しかつ共通のグループ識別子の割り当てを受ける。尤度関数は、（グループ化された複数のトレーニングサンプルから生成されるデータを組み込んだ）一次結合を条件とする多項又は２項ロジスティック回帰に基づいて定式化される。特徴は、トレーニングサンプルをグループ化すべき際に使用される、所望される任意の分類であることが可能である。例えば組織サンプルを分類するための特徴は、その組織が正常、悪性、良性、白血病細胞、健康細胞であるという場合もあれば、所定の状態を有するかもしくは持たない患者らの血液から複数のトレーニングサンプルが採取されるという場合もあり、又は、正常細胞に比較されるいくつかのタイプの癌のうちの１つの細胞から複数のトレーニングサンプルが採取されるという場合もある。

上記第１の実施形態において、多項又は２項ロジスティック回帰に基づく尤度関数は、次式の形式である。

ここで、ｘ_ｉ ^Ｔβ_ｇは、構成要素の複数の重み係数β_ｇとともにトレーニングサンプルｉの入力データから生成される一次結合であり、ｘ_ｉ ^ＴはＸのｉ番目の行の要素であり、β_ｇはサンプルクラスｇに係る構成要素の複数の重み係数にてなるセットであり、Ｘはｐ個の要素を含むｎ個のトレーニングサンプルからのデータであり、ｅ_ｉｋは本明細書内で後に定義される。

第２の実施形態では、尤度関数は順序分類ロジスティック回帰（ordered categorical logistic regression）に基づく。順序分類ロジスティック回帰は、複数のクラスが特定の順序で存在する（例えば、病気の重さが次第に重くなるか又は次第に軽くなる複数のクラス等の、順序付けられたクラスが存在する）２項又は多項分布をモデル化する。順序分類ロジスティック回帰に基づいて尤度関数を定義すれば、複数の予め定義された順序付けられたクラスのうちの１つであるクラスにサンプルを分類することができる複数の構成要素のうちのサブセットを識別することが可能である。順序付けられたクラスの要素にそれぞれ相当する一連のグループ識別子を定義して、複数のトレーニングサンプルの予め決められた特徴に基づいて当該複数のトレーニングサンプルを順序付けられたクラスのうちの１つにグループ化することにより、尤度関数は、（グループ化された複数のトレーニングサンプルから生成されるデータを組み込んだ）一次結合を条件とする順序分類ロジスティック回帰に基づいて定式化されることが可能である。

上記第２の実施形態では、順序分類ロジスティック回帰に基づく尤度関数は、次式の形式である。

ここで、γ_ｉｋは、ｋ以下の識別子を有するクラスにトレーニングサンプルｉが属する確率である（ここで、順序付けられたクラスの総計はＧである。）。ｒ_ｉは本明細書において後に定義される。

本発明の第３の実施形態では、尤度関数は一般化線形モデルに基づく。一般化線形モデルは、好適には、通常の指数型分布族（regular exponential family of distributions）として分布される特徴をモデル化する。通常の指数型分布族の例には、正規分布、ガウス分布、ポアソン分布、ガンマ分布、及び逆ガウス分布が含まれる。従って、本発明の方法に係るもう１つの実施形態では、予測されるべき特徴をモデル化する一般化線形モデルを定義することによって特に、通常の指数型分布族に属する分布を有するサンプルの予め定義された特性を予測することができる、複数の構成要素のうちのサブセットが識別される。一般化線形モデルを用いて予測される可能性のある特性の例には、例えばサンプルの重さ、サイズ又は他の寸法もしくは量などの指定された分布を呈示する、サンプルに係る任意の量が含まれる。

上記第３の実施形態では、一般化線形モデルは、次式の形式である。

ここで、ｙ＝（ｙ_１，…，ｙ_ｎ）^Ｔであり、ａ_ｉ（φ）＝φ／ｗ_ｉであり、このときｗ_ｉは既知の重み係数にてなる固定されたセットであり、φは単一のスケールパラメータである。この式における他の項は本明細書において後に定義される。

第４の実施形態では、本発明の方法は、ハザードモデルに基づく尤度関数を利用することにより、あるサンプルに係るイベントまでの時間を予測するために使用可能であり、これは好適には、上記イベントがデータ取得時点において発生していないことを条件にイベントまでの時間の確率を推定する。この第４の実施形態では、尤度関数は、コックスの比例ハザードモデル、パラメトリック生存モデル、及び加速度故障回数モデルを含むグループから選択される。コックスの比例ハザードモデルは、イベントまでの時間が、時間に関して制限的仮定を行うことなしに複数の構成要素と構成要素に係る複数の重み係数とのセットに基づいてモデル化されることを可能にする。加速度故障モデルは、複数の生存時間より成るデータのための一般モデルであって、ここでは、構成要素に係る複数の測定値は、時間スケールに対して乗法的に増大するように作用し、よって時間軸に沿って個々に進行する速度（レート）に影響を与えることが仮定される。従って加速度生存モデルは、例えば病気の進行速度に置き換えて解釈されることが可能である。パラメトリック生存モデルは、イベントまでの時間（例えば生存時間）の分布関数が既知の分布によってモデル化されたり、又は特定の（指定された）パラメトリックな定式化を有したりするというものである。一般に使用されている生存分布には、ワイブル（Weibull）分布、指数分布、及び極値分布がある。

上記第４の実施形態では、あるサンプルに係るイベントまでの時間を予測することができる複数の構成要素のうちのサブセットが、コックスの比例標準（proportional standards）モデル、パラメトリック生存モデル、又は加速度生存時間モデルに基づいて尤度を定義することにより識別される。これには、複数のサンプルに関する、サンプルの取得時刻からイベント発生時刻までの経過時間を測定することが含まれる。

上記第４の実施形態では、イベントまでの時間を予測するための尤度関数は、次式の形式である。

ここで、

及び

はモデルパラメータであり、ｙは観測された複数の時刻にてなるベクトルであり、ｃは、ある時間が真の生存時間であるか、それとも打ち切り生存時間（censored survival time）であるかを示す指示子ベクトルである。

上記第４の実施形態では、コックスの比例ハザードモデルに基づく尤度関数は、次式の形式である。

ここで、観測された時刻は昇順で並べられて

で表され、ＺはＸの行の並べ替えであるＮ×ｐ行列を示し、Ｚの行の順序づけは

の順序づけによって導かれた順序づけに対応する。また、

であり、ｚ_ｊはｚのｊ番目の行であり、

は、ｊ番目の順序を有するイベント時刻ｔ_（ｊ）に設定されるリスクである。

尤度関数がパラメトリック生存モデルに基づく第４の実施形態では、尤度関数は、次式の形式である。

ここで、

であり、Λは積分されたパラメトリックハザード関数を示す。

定義される任意のモデルについて、重み係数は、典型的には、ベイズ統計モデルを用いて推定され（コッツ及びジョンソン（Kots and Johnson），１９８３年）、この場合、構成要素に係る複数の重み係数の事後分布であって、尤度関数と事前分布とを組み合わせる事後分布が定式化される。構成要素に係る複数の重み係数は、少なくとも１つのトレーニングサンプルに関して生成されるデータを所与として、複数の重み係数の事後分布を最大化することにより推定される。従って、最大化されるべき目的関数は、先に論じたような特徴に関するモデルに基づく尤度関数と、複数の重み係数の事前分布とから成る。

好適には、事前分布は、次式の形式である。

ここで、ｖは複数のハイパーパラメータにてなるｐ×１ベクトルであり、ｐ（β│ｖ^２）はＮ（０，ｄｉａｇ｛ｖ^２｝）であり、ｐ（ｖ^２）はｖ^２に関する何らかの超事前分布（hyperprior distribution）である。

好適には、超事前分布は、指定された形状及びスケールパラメータを有するガンマ分布を含む。

この超事前分布（これは、好適には本方法の全ての実施形態に関して同じである。）は、異なる表記法を用いて表されることが可能であり、実施形態の詳細な説明（下記参照）では、単に便宜上、特定の実施形態に関して次のような表記法を採用している。

本明細書で使用しているように、確率分布の尤度関数が多項又は２項ロジスティック回帰に基づく場合、事前分布は次式のように表記される。

ここで、β^Ｔ＝（β_１ ^Ｔ，…，β_Ｇ−１ ^Ｔ）及びτ^Ｔ＝（τ_１ ^Ｔ，…，τ_Ｇ−１ ^Ｔ）であり、ｐ（β_ｇ│τ_ｇ ^２）はＮ（０，ｄｉａｇ｛τ_ｇ ^２｝）であり、Ｐ（τ_ｇ ^２）はτ_ｇ ^２に関する何らかの超事前分布である。

本明細書で使用しているように、確率分布の尤度関数が順序分類ロジスティック回帰に基づく場合、事前分布は次式のように表記される。

ここで、β_１，β_２，…，β_ｎは構成要素の重み係数であり、Ｐ（β_ｉ│ｖ_ｉ）はＮ（０，ｖ_ｉ ^２）であり、Ｐ（ｖ_ｉ）はｖ_ｉに関する何らかの超事前分布である。

本明細書で使用しているように、上記分布の尤度関数が一般化線形モデルに基づく場合、事前分布は次式のように表記される。

ここで、ｖは複数のハイパーパラメータに係るｐ×１ベクトルであり、ｐ（β│ｖ^２）はＮ（０，ｄｉａｇ｛ｖ^２｝）であり、ｐ（ｖ^２）はｖ^２に関する何らかの超事前分布である。

本明細書で使用しているように、上記分布の尤度関数がハザードモデルに基づく場合、事前分布は次式のように表記される。

ここで、ｐ（β^＊│τ）はＮ（０，ｄｉａｇ｛τ^２｝）であり、ｐ（τ）はτに関する何らかの超事前分布である。

事前分布は、可能であればいつでもゼロの重み係数が使用されることを保証する、超事前分布を含む。

ある代替実施形態では、超事前分布は、ｔ_ｉ ^２＝１／ｖ_ｉ ^２が各々独立したガンマ分布を有する逆ガンマ分布である。

あるさらなる代替実施形態では、超事前分布は、（コンテキストに依存して）ｖ_ｉ ^２、τ_ｉ又はτ_ｉ ^２が各々独立したガンマ分布を有するガンマ分布である。

先に論じたように、事前分布と尤度関数とは組み合わされて事後分布を生成する。事後分布は、好適には、
［数１］
ｐ（βφｖ│ｙ）αＬ（ｙ│βφ）ｐ（β│ｖ）ｐ（ｖ）
又は、

の形式である。ここで、

は尤度関数である。

好適には、複数の構成要素のうちのサブセットを識別するステップは、事後分布の確率密度が最大化されるように反復手順を使用するステップを含む。

上記反復手順の実行中において、予め決められたしきい値を下回る値を有する、構成要素に係る複数の重み係数は、好適にはこれらの構成要素の重み係数をゼロに設定することにより除去される。これにより、対応する構成要素は実質上除去されることになる。

好適には、上記反復手順はＥＭアルゴリズムである。

ＥＭアルゴリズムは、事後分布の確率密度を最大化する重み係数を構成要素に与えるように収束する、構成要素の重み係数に係る一連の推定値を生成する。ＥＭアルゴリズムは、Ｅステップ又は期待値計算ステップと、Ｍステップ又は最大化ステップとして知られた２つのステップより成る。Ｅステップでは、観測データを条件とする対数事後関数の期待値が決定される。Ｍステップでは、更新された構成要素の重み係数に係る複数の期待値であってかつ事後分布を増大させる推定値を与えるように、期待される対数事後関数が最大化される。２つのステップは、Ｅステップ及びＭステップの収束が達成されるまで、又は言い替えれば、期待値と、期待される対数事後関数の最大値とが収束するまで交互に行われる。

本発明に係る方法は、測定値の取得先とすることが可能な任意のシステムに適用され得ること、また好適には膨大な量のデータの生成元となるシステムに適用され得ることが想定されている。本発明の方法を適用可能なシステムの例には、生物学的システムと、化学システムと、農業システムと、気象システムと、例えば信用リスク評価システム、保険システム、マーケティングシステム又は企業記録システムを含む財務システムと、電子的システムと、物理的システムと、宇宙物理的システムと、機械的システムとが含まれる。例えば、財務システムでは、サンプルは特定の株式であることが可能であり、構成要素は、企業収益、従業員数、さまざまな都市の降水量、株主数などの、株価に影響を与える可能性のある任意個数のファクタに関して求められる測定値である可能性がある。

本発明の方法は、特に、生物学的システムの分析における使用に適している。本発明の方法は、構成要素の測定可能値を生成する任意の生物学的システムからのサンプルを分類するための複数の構成要素にてなる複数のサブセットであって、構成要素が一意的にラベリングされることが可能なサブセットを識別するために使用可能である。言い替えれば、上記複数の構成要素は、１つの構成要素からのデータを別の構成要素からのデータから区別できるようにラベリングされるか又は編成される。例えば、複数の構成要素は、各構成要素からのデータを空間的位置によって別のものから区別できるように空間的に編成される、例えばアレイ状に編成される場合もあれば、各構成要素は、識別信号又はタグ等の、当該構成要素に関連づけられる何らかの一意的な識別を有する場合もある。例えば、構成要素は、検出可能な識別サインをそれぞれ有する個々のキャリアに拘束される場合もある。識別サインとしては、例えば、量子ドット（例えば、「ローゼンソール，２００１年，ネイチャーバイオテック１９：６２１−６２２（Rosenthal, 2001, Nature Biotech 19: 621-622）」、「ハンほか（２００１年）ネイチャーバイオテクノロジー１９：６３１−６３５（Han et al. (2001) Nature Biotechnology 19: 631-635）」を参照）、蛍光マーカ（例えば、「フーほか（１９９９年）ネイチャーバイオテクノロジー１７：１１０９−１１１１（Fu et al. (1999) Nature Biotechnology 17: 1109-1111）」を参照）、バーコード付きタグ（例えば、「ロックハート及びトラルソン（２００１年）ネイチャーバイオテクノロジー１９：１１２２−１１２３参照（Lockhart and trulson (2001) Nature Biotechnology 19: 1122-1123）」を参照）がある。

ある特に好適な実施形態では、生物学的システムはバイオテクノロジーアレイである。バイオテクノロジーアレイの例には、オリゴヌクレオチドアレイ、ＤＮＡアレイ、ＤＮＡマイクロアレイ、ＲＮＡアレイ、ＲＮＡマイクロアレイ、ＤＮＡマイクロチップ、ＲＮＡマイクロチップ、タンパク質アレイ、タンパク質マイクロチップ、抗体アレイ、化学アレイ、炭水化物アレイ、プロテオミクスアレイ、脂質アレイが含まれる。別の実施形態では、生物学的システムは、例えば、ＤＮＡ又はＲＮＡ電気泳動ゲルと、タンパク又はプロテオミクス電気泳動ゲルと、ビアコア（Biacore）分析等の生体分子間相互作用の分析と、アミノ酸分析と、ＡＤＭＥＴｏｘスクリーニング（例えば、フェレンツ・ダルバシュ及びジェルジ・ドルマーン編，「ハイスループットＡＤＭＥＴｏｘ推定法：インビトロ及びインシリコ的アプローチ」，バイオテクニークスプレス，２００２年（High-throughput ADMETox estimation: In Vitro and In Silico approaches (2002), Ferenc Darvas and Gyorgy Dorman (Eds), Biotechniques Press）を参照）と、タンパク質電気泳動ゲルと、プロテオミクス電気泳動ゲルとを含むグループから選択されることが可能である。

構成要素は、システムに係る測定可能な任意の構成要素であってよい。生物学的システムのケースでは、構成要素は、例えば、遺伝子もしくはその一部、ＤＮＡ配列、ＲＮＡ配列、ペプチド、タンパク質、炭水化物分子、脂質もしくはその混合物、生理学的構成要素、解剖学的構成要素、疫学的構成要素、又は化学的構成要素である可能性がある。

トレーニングサンプルは、サンプルの特徴が既知であるシステムから取得される任意のデータであることが可能である。例えばトレーニングサンプルは、生物学的システムに適用されるサンプルから生成されるデータである可能性がある。例えば、生物学的システムがＤＮＡマイクロアレイであれば、トレーニングサンプルは、そのアレイと、既知の特徴を有する細胞から抽出されたＲＮＡとのハイブリダイゼーション、又は細胞から抽出されたＲＮＡから合成されたｃＤＮＡとのハイブリダイゼーションに続いて当該アレイから取得されるデータである場合もあり、生物学的システムがプロテオミクス電気泳動ゲルであれば、トレーニングサンプルは、システムに適用されるタンパク質又は細胞抽出物から生成される場合もある。

本発明の方法の実施形態は、試験処理剤への反応において雑多な結果を示した検査対象からの試験データを再評価する、又は評価する際に使用され得ることが想定される。こうして、本発明には第２の態様が存在する。

第２の態様は、検査対象を複数の予め定義されたグループのうちの１つに分類することができる、検査対象に係る複数の構成要素のうちのサブセットを識別するための方法を提供し、各グループは試験処理剤への反応によって定義され、上記方法は、
複数の検査対象を試験処理剤にさらし、上記処理剤に対する反応に基づいて上記検査対象を複数の反応グループにグループ化するステップと、
上記検査対象の複数の構成要素を測定するステップと、
統計解析方法を用いて、上記検査対象を反応グループに分類することのできる構成要素のうちのサブセットを識別するステップとを含む。

好適には、上記統計解析方法は本発明の第１の態様に係る方法を含む。

いったん複数の構成要素のうちのサブセットが識別されると、そのサブセットは、試験処理剤に反応する可能性のあるグループとそうでないグループ等の複数のグループに検査対象を分類するために使用可能である。このようにして、本発明の方法は、母集団のうちの一部に有効である可能性のある処理剤の識別を可能にし、また、その試験処理剤に反応する母集団のその一部の識別を可能にする。

本発明の第３の態様によれば、検査対象に係る複数の構成要素のうちのサブセットを識別するための装置が提供されていて、上記サブセットは、上記検査対象を複数の予め定義された反応グループのうちの１つに分類するために使用可能であり、各反応グループは、複数の検査対象を試験処理剤にさらしかつ上記処理剤に対する反応に基づいて上記検査対象を複数の反応グループにグループ化することによって形成され、上記装置は、
上記検査対象に係る測定された複数の構成要素を受信するための入力と、
統計解析方法を用いて、上記検査対象を反応グループに分類するために使用可能な構成要素のうちのサブセットを識別する処理手段とを備える。

好適には、上記統計解析方法は上記第１又は第２の態様に係る方法を含む。

本発明の第４の態様によれば、試験化合物による処理に対して反応するもの又は反応しないものとして検査対象を分類することができる検査対象に係る複数の構成要素のうちのサブセットを識別するための方法が提供されていて、上記方法は、
複数の検査対象を試験化合物にさらし、上記試験化合物に対する各検査対象の反応に基づいて上記検査対象を複数の反応グループにグループ化するステップと、
上記検査対象に係る複数の構成要素を測定するステップと、
統計解析方法を用いて、上記検査対象を反応グループに分類するために使用可能な複数の構成要素のうちのサブセットを識別するステップとを含む。

好適には、上記統計解析方法は上記第１の態様に係る方法を含む。

本発明の第５の態様によれば、検査対象に係る複数の構成要素のうちのサブセットを識別するための装置が提供されていて、上記サブセットは、上記検査対象を複数の予め定義された反応グループのうちの１つに分類するために使用可能であり、各反応グループは、複数の検査対象を化合物にさらしかつ上記化合物に対する反応に基づいて上記検査対象を複数の反応グループにグループ化することによって形成され、上記装置は、
上記検査対象に係る測定された複数の構成要素を受信する入力と、
統計解析方法を用いて、上記検査対象を反応グループに分類することができる複数の構成要素のうちのサブセットを識別する処理手段とを備える。

好適には、上記統計解析方法は本発明の第１又は第２の態様に係る方法を含む。

本発明の上記第２乃至第５の態様において測定される構成要素は、例えば、遺伝子もしくは小塩基多型（small nucleotide polymorphism：ＳＮＰ）、タンパク質、抗体、炭水化物、脂質、又は検査対象の他の任意の測定可能な構成要素であることが可能である。

上記第５の態様の特別な実施形態では、化合物は、医薬化合物であるか、又は医薬化合物と薬剤を受容可能なキャリアとを備えた組成物である。

本発明に係る識別方法は、適切なコンピュータソフトウェア及びハードウェアによって実装されることが可能である。

本発明の第６の態様によれば、システムの複数のサンプルから生成されるデータからシステムに係る複数の構成要素のうちのサブセットを識別するための装置が提供されていて、上記サブセットは試験サンプルの特徴を予測するために使用可能であり、
上記装置は処理手段を備え、上記処理手段は、
上記システムに係る複数の構成要素の一次結合を取得し、上記複数の構成要素の一次結合に係る複数の重み係数を取得するように動作し、上記重み係数の各々は少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の特徴を有し、
第２の特徴の確率分布のモデルを取得するように動作し、上記モデルは上記複数の構成要素の一次結合を条件とし、
上記複数の構成要素の一次結合に係る複数の重み係数に関する事前分布を取得するように動作し、上記事前分布は、ゼロに近い事前確率質量が変更されることを可能にする調節可能な超事前分布を含み、上記超事前分布はジェフリーズの超事前分布ではなく、
上記事前分布と上記モデルとを組み合わせて事後分布を生成するように動作し、
上記事後分布を最大化する構成要素の重み係数を有する複数の構成要素のうちのサブセットを識別するように動作する。

好適には、上記処理手段はソフトウェアを実行するように構成されたコンピュータを備える。

本発明の第７の態様によれば、計算装置によって実行されたときに、本発明の第１の態様に係る方法を上記計算装置に実行させるコンピュータプログラムが提供されている。

上記コンピュータプログラムは、好適なアルゴリズムと、先に論じた本発明の第１又は第２の態様に係る方法ステップとのうちの任意のものを実施することができる。

本発明の第８の態様によれば、本発明の上記第７の態様によるコンピュータプログラムを備えた、コンピュータが読み取り可能な媒体が提供されている。

本発明の第９の態様によれば、システムからのサンプルを検査してサンプルの特徴を識別する方法が提供されていて、
上記方法は、上記特徴の症状を示す複数の構成要素のうちのサブセットについて検査するステップを含み、上記複数の構成要素のうちのサブセットは本発明の第１又は第２の態様に係る方法を用いて決定されている。

好適には、上記システムは生物学的システムである。

本発明の第１０の態様によれば、システムからのサンプルを検査してサンプルの特徴を決定するための装置が提供されていて、上記装置は、本発明の第１又は第２の態様に係る方法に従って識別される構成要素を検査するための手段を備える。

本発明の第１１の態様によれば、計算装置によって実行されたときに、システムからの試験サンプルの特徴を予測するために使用可能なシステムからの構成要素を識別する方法を上記計算装置に実行させるコンピュータプログラムが提供されていて、上記方法においては、
複数のトレーニングサンプルから生成されるデータから複数の構成要素と構成要素に係る複数の重み係数との一次結合が生成され、各トレーニングサンプルは既知の特徴を有し、
ゼロに近い確率質量が変更されることを可能にする調節可能な超事前分布を含む構成要素に係る複数の重み係数の事前分布と、上記一次結合を条件とするモデルとを組み合わせ、上記事後分布を最大化する構成要素に係る複数の重み係数を推定することによって、事後分布が生成される。ここで、上記超事前分布はジェフリーズの超事前分布ではない。

本発明の態様は計算装置によって実装されるが、任意の適切なコンピュータハードウェア、例えばＰＣ又はメインフレームあるいはネットワークで接続されたコンピューティングインフラストラクチャが使用可能であることは認識されるであろう。

本発明の第１２の態様によれば、生物学的システムに係る複数の構成要素のうちのサブセットを識別する方法が提供されていて、上記サブセットは上記生物学的システムからの試験サンプルの特徴を予測することができ、上記方法は、
上記システムに係る複数の構成要素と、上記複数の構成要素の一次結合に係る複数の重み係数との一次結合を取得するステップを含み、上記重み係数の各々は少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の第１の特徴を有し、
第２の特徴の確率分布のモデルを取得するステップを含み、上記モデルは上記複数の構成要素の一次結合を条件とし、
上記複数の構成要素の一次結合に係る複数の重み係数に関する事前分布を取得するステップを含み、上記事前分布は、ゼロに近い確率質量が変更されることを可能にする調節可能な超事前分布を含み、
上記事前分布と上記モデルとを組み合わせて事後分布を生成するステップと、
上記事後分布を最大化する重み係数に基づいて複数の構成要素のうちのサブセットを識別するステップとを含む。

本発明の範囲内にある可能性のある他の任意の実施形態に関わりなく、以下、添付の図面を参照して、本発明の実施形態を単なる例示としてのみ説明する。

本発明の実施形態は、特定のトレーニングサンプルがある特徴を有するか否かを識別するために使用可能な、比較的少数の構成要素を識別する。これらの構成要素はその特徴の「症状」を示し、又はこれらの構成要素は、異なる特徴を有するサンプル間での区別を可能にする。本方法によって選択される構成要素の数は、超事前分布におけるパラメータの選択によって制御されることが可能である。超事前分布は、指定された形状とスケールパラメータとを有するガンマ分布であることが知られている。本質的に本発明の方法は、システムから生成される全てのデータから、特定の特徴の検査に使用可能な比較的少数の構成要素の識別を可能にする。いったん本方法によりこれらの構成要素が識別されると、上記構成要素は、将来に新たなサンプルを評価するために使用可能である。本発明の方法は統計的方法を利用して、特徴を正しく予測するためには不要である構成要素を除去する。

本発明者らは、複数のトレーニングサンプルから生成されるデータに係る複数の構成要素の一次結合における構成要素に係る複数の重み係数が、トレーニングサンプルの特徴を正しく予測するためには不要である構成要素を除去するような方法で推定され得ることを発見している。その結果、トレーニングサンプルの特徴を正しく予測することのできる複数の構成要素のうちのサブセットが識別される。従って本発明の方法は、大量のデータから、ある特徴を正しく予測することができる比較的少なくかつ制御可能な個数の構成要素を識別することを可能にする。

また本発明の方法は、必要とするコンピュータメモリの使用量が従来技術の方法よりも少ないという優位点を有する。従って、本発明の方法は、例えばラップトップマシン等のコンピュータ上で高速に実行されることが可能である。より少ないメモリの使用により、本発明の方法はまた、本方法が、例えば生物学的データを解析するために複数の構成要素に関する（周辺的（marginal）情報でなく）同時的（joint）情報を用いる他の方法より高速で実行されることを可能にする。

本発明の方法はまた、解析のために複数の構成要素に関する周辺的情報ではなく同時的情報を用いるという優位点をも有する。

次に、マルチクラスのロジスティック回帰モデルに関する第１の実施形態について説明する。

Ａ．マルチクラスロジスティック回帰モデル．
この実施形態に係る方法は、複数のトレーニングサンプルを用いて、上記トレーニングサンプルを複数の予め定義されたグループに分類することのできる複数の構成要素のうちのサブセットを識別する。続いて上記構成要素のうちのサブセットについての知識は、複数のサンプルを疾病クラス等の複数のグループに分類するための試験、例えば臨床試験に使用可能である。例えば、ＤＮＡマイクロアレイの複数の構成要素のうちのサブセットは、複数の臨床サンプルを、例えば健康又は病気等の臨床に関連した複数のクラスにグループ化するために使用可能である。

このようにして、本発明は、特定のトレーニングサンプルがある特定のグループに属するか否かを識別するために使用可能な、好適には少量かつ制御可能な個数の構成要素を識別する。選択される構成要素はそのグループの「症状」を示し、又は、それらの選択される構成要素は複数のグループ間での区別を可能にする。本質的に本発明の方法は、システムから生成される全てのデータから、特定のグループの検査に使用可能な少数の構成要素の識別を可能にする。いったん本方法によりこれらの構成要素が識別されると、上記構成要素は、将来に新たなサンプルをグループに分類する際に使用可能である。本発明の方法は好適には、統計的方法を用いて、上記サンプルが属するグループを正しく識別するためには不要である構成要素を除去する。

複数のサンプルは、予め決められた分類法に基づいて複数のサンプルグループ（又は「クラス」）にグループ化される。この分類法は、トレーニングサンプルがグループ化されるときに使用される、任意の所望の分類法であってよい。例えば分類法は、トレーニングサンプルが白血病細胞からのものかそれとも健康な細胞からのものかという場合もあれば、トレーニングサンプルが、所定の状態を有する患者又は有していない患者の血液から取得されること、あるいは、トレーニングサンプルが、正常な細胞との比較でいくつかのタイプの癌のうちの１つからの細胞によるものであることという場合もある。

ある実施形態では、入力データは、ｎ個のトレーニングサンプル及びｐ個の構成要素が存在する場合のｎ×ｐデータ行列Ｘ＝（ｘ_ｉｊ）に編成される。典型的には、ｐはｎよりずっと大きくなる。

別の実施形態では、データ行列Ｘは、線形予測量ではなく、予測量としてＸの滑らかな関数を得るためにｎ×ｎ核行列Ｋで置換されることが可能である。核行列Ｋの一例は、次式になる。

［数２］
ｋ_ｉｊ＝ｅｘｐ（−０．５＊（ｘ_ｉ−ｘ_ｊ）^ｔ（ｘ_ｉ−ｘ_ｊ）／σ^２）

ここで、ｘの下付き添字は行列Ｘにおける行の番号を示す。理想的には、Ｋの列のうちのサブセットは、これらの滑らかな関数の疎な表現を与えるものが選択される。

各サンプルクラス（グループ）に付随して、トレーニングサンプルがＧ個のサンプルクラスのうちのどれに属するかを示すクラスラベルｙ_ｉが存在する。ここで、ｙ_ｉ＝ｋ，ｋ∈｛１，…，Ｇ｝である。ここでは、要素ｙ_ｉを備えたｎ×１ベクトルをｙ．と記す。ベクトル

を所与とすると、指示子変数を次式のように定義することができる。

（Ａ１）

ある実施形態では、構成要素の重み係数はベイズ統計モデルを用いて推定される（コッツ及びジョンソン，１９８３年を参照）。好適には、上記重み係数は、各トレーニングサンプルから生成されるデータを所与として、上記重み係数の事後分布を最大化することにより推定される。これにより、最大化される目的関数は２つの部分からなる。その第１の部分は尤度関数であり、その第２の部分は複数の重み係数の事前分布であり、これは可能であればいつでもゼロの重み係数が好適であることを保証する。ある好適な実施形態では、尤度関数はマルチクラスのロジスティックモデルから導出される。好適には、尤度関数は次式の確率から計算される。

（Ａ２）
及び、

（Ａ３）

ここで、Ｐ_ｉｇは、入力データＸ_ｉを有するトレーニングサンプルがサンプルクラスｇ内に存在することになる確率であり、ｘ_ｉ ^Ｔβ_ｇは、構成要素の重み係数β_ｇを有するトレーニングサンプルｉからの入力データから生成される一次結合であり、ｘ_ｉ ^Ｔは、Ｘのｉ番目の行の要素であり、β_ｇはサンプルクラスｇの構成要素に係る複数の重み係数のセットである。

典型的には、先に論じたように、構成要素の重み係数は、大部分の構成要素の重み係数がゼロであるという先験的仮定を考慮した方法で推定される。

ある実施形態では、式（Ａ２）における構成要素の重み係数βｇは、大部分の値がゼロであるように推定されるが、それでもなおサンプルは正確に分類されることが可能である

ある実施形態では、構成要素に係る複数の重み係数は、先に言及したベイズモデルにおけるデータを所与として、それらの重み係数の事後分布を最大化することにより推定される。

好適には、構成要素の重み係数は、
（ａ）構成要素の重み係数β_１，…，β_Ｇ−１の階層的事前分布を指定することと、
（ｂ）入力データの尤度関数を指定することと、
（ｃ）上記データが与えられたときの重み係数の事後分布を、（Ａ５）を用いて決定することと、
（ｄ）上記事後分布を最大化する構成要素の重み係数を決定することと
により推定される。

ある実施形態では、パラメータβ_１，…，β_Ｇ−１に関して指定される階層的事前分布は、次式の形式である。

（Ａ４）

ここで、β^Ｔ＝（β_１ ^Ｔ，…，β_Ｇ−１ ^Ｔ）及びτ^Ｔ＝（τ_１ ^Ｔ，…，τ_Ｇ−１ ^Ｔ）であり、ｐ（β_ｇ│τ_ｇ ^２）はＮ（０，ｄｉａｇ｛τ_ｇ ^２｝）であり、ｐ（τ_ｇ ^２）は適切な事前分布である。

ある実施形態では、

である。ここで、ｐ（τ_ｉｇ ^２）は事前分布であり、ｔ_ｉｇ ^２＝１／τ_ｉｇ ^２は独立なガンマ分布を有する。

別の実施形態では、ｐ（τ_ｉｇ ^２）は事前分布であり、τ_ｉｇ ^２が独立なガンマ分布を有する。

ある実施形態では、尤度関数は式（８）における形式の

であり、ｙを所与とするβ及び

の事後分布は、次式になる。

（Ａ５）

ある実施形態では、尤度関数は１階及び２階の導関数を有する。

ある実施形態では、上記１階の導関数は、次のようなアルゴリズムから決定される。

（Ａ６）

ここで、

は、サンプルクラスｇの帰属関係と、クラスｇの確率とをそれぞれ示すベクトルである。

ある実施形態では、上記２階の導関数は、次のようなアルゴリズムから決定される。

（Ａ７）

ここでδ_ｈｇは、ｈがｇに等しければ１であり、そうでなければゼロである。

式Ａ６及び式Ａ７は、次のようにして導出されることが可能である。

（ａ）データの尤度関数は、（Ａ１）、（Ａ２）及び（Ａ３）を用いて、

（Ａ８）
のように書き表すことができる。

（ｂ）式（Ａ６）の対数をとり、また全てのｉについて

であるという事実を用いると、

（Ａ９）
が与えられる。

（ｃ）式（Ａ８）をβｇに関して微分すると、

（Ａ１０）
が与えられる。ここで、

は、サンプルクラスｇの帰属関係と、クラスｇの確率とをそれぞれ示すベクトルである。

（ｄ）式（９）の２階の導関数は、要素、

（Ａ１１）
を有する。ここで、

である。

尤度関数の事後分布を最大化する構成要素の重み係数は、ＥステップとＭステップとを含むＥＭアルゴリズムを用いて特定されることが可能である。

ＥＭアルゴリズムの実行に際しては、Ｅステップは好適には、次式の形式の項を計算するステップを含む。

（Ａ１１ａ）

ここで、

であり、

であれば

である。

好適には、ｐ（β_ｉｇ│τ_ｉｇ ^２）がＮ（０，τ_ｉｇ ^２）でありかつｐ（τ_ｉｇ ^２）が指定された事前分布を有するとき、式（１１ａ）はｔ_ｉｇ ^２＝１／τ_ｉｇ ^２の条件付き期待値を計算することによって演算される。条件付き期待値の明示的な式は、後に提示する。

典型的には、ＥＭアルゴリズムは下記のようなステップを含む。

（ａ）関数

（Ａ１２）
を用いて、構成要素に係る複数の重み係数の事後分布の条件付き期待値を計算することにより、Ｅステップを実行する。ここで、式（８）においてｘ_ｉ ^Ｔβ_ｇ＝ｘ_ｉ ^ＴＰ_ｇγ_ｇであり、

であり、

は

で評価された式（１１ａ）の場合と同様に定義される。この場合のＰ_ｇは、γ_ｇで示されるβ_ｇの非ゼロ要素をＰ_ｇ ^Ｔβ_ｇが選択するように、恒等行列から導出されるゼロ及び１を要素とする行列である。

（ｂ）反復手順を適用してγの関数としてＱを最大化することにより、Ｍステップを実行する。よって、

（Ａ１３）
となる。ここで、α^ｔは０≦α^ｔ≦１であるようなステップ長であり、γ＝（γ_ｇ，ｇ＝１，…，Ｇ−１）である。

式（Ａ１２）は、次のように導出されることが可能である。

観測データｙと、パラメータ推定値のセット

とを所与として、（Ａ５）の条件付き期待値を計算する。

β（及び

）の要素がゼロに設定されるときの、すなわち、ｇ＝１，…，Ｇ−１についてβ_ｇ＝Ｐ_ｇγ_ｇ及び

であるケースについて考察する。

γを含まない項を無視しかつ（Ａ４）、（Ａ５）、（Ａ９）を使用すれば、次式が得られる。

（Ａ１４）

ここで、（Ａ８）において

であり、

であり、

は

で評価された式（Ａ１１ａ）の場合と同様に定義される。

条件付き期待値は、（Ａ４）で与えられる第１の原則から評価され得ることに留意されたい。明示的な式については、後にいくつか述べる。

反復手順は、次のようにして導出されることが可能である。

（１１）で必要とされる導関数を得るためには、まず、

と表記するときに、（Ａ８）、（Ａ９）及び（Ａ１０）から次式が得られることに留意されたい。

（Ａ１５）
及び、

（Ａ１６）

ここで、

であり、かつ、
［数３］
Ｘ_ｇ ^Ｔ＝Ｐ_ｇ ^ＴＸ^Ｔ，ｇ＝１，…，Ｇ−１
（Ａ１７）
である。

ある好適な実施形態では、上記反復手順は、式（Ａ１３）において式（Ａ１６）のブロックの対角要素のみを用いることにより簡略化されることが可能である。するとこれは、ｇ＝１，…，Ｇ−１について、次式を与える。

（Ａ１８）

式（Ａ１８）を変形すると、次式が得られる。

（Ａ１９）

ここで、

である。

Ｙ_ｇの列の数をｐ（ｇ）と書くと、（Ａ１９）は、ｐ（ｇ）×ｐ（ｇ）行列の逆行列の演算を必要とし、これは非常に大規模なものになる可能性がある。これは、ｐ（ｇ）＞ｎのときに、

（Ａ２０）
に注目することにより、ｎ×ｎ行列まで縮小されることが可能である。ここで、Ｚ_ｇ＝Δ_ｇｇ ^２Ｙ_ｇである。好適には、（Ａ１９）はｐ（ｇ）＞ｎのときに使用され、式（Ａ１９）へ（Ａ２０）が代入された形の（Ａ１９）は、ｐ（ｇ）≦ｎのときに使用される。

τ_ｉｇ ^２がジェフリーズ（Jeffreys）の事前分布を有するとき、次式が得られることに留意されたい。

ある実施形態では、ｔ_ｉｇ ^２＝１／τ_ｉｇ ^２は、スケールパラメータがｂ＞０であり形状パラメータがｋ＞０である独立なガンマ分布を有し、よってｔ_ｉｇ ^２の密度は次式になる。

表記を簡単化するために下付き添字を省略し、
［数４］
Ｅ｛ｔ^２｜β｝＝（２ｋ＋１）／（２／ｂ＋β^２）
（Ａ２１）
となることを証明することができる。その手順は、下記の通りである。

を定義する。すると、

になる。

証明．
ｓ＝β^２／２とすると、

になる。ここで、ｕ＝ｔ^２／ｂを代入すると、

が得られる。次に、ｓ’＝ｂｓとして、γ（ｕ，ｌ，ｋ）の式を代入すると、

が得られる。結果は、例えばアブラモビッツ（Abramowitz）及びステガン（Stegun）のラプラス変換表を参照することによって得られる。
条件付き期待値は、
［数５］
Ｅ｛ｔ^２｜β｝＝Ｉ（１，ｂ，ｋ）／Ｉ（０，ｂ，ｋ）
＝（２ｋ＋１）／（２／ｂ＋β^２）
から得られる。

ｋはゼロへ向かい、ｂは無限大へ向かうとき、ジェフリーズの事前分布を用いる場合と同等の結果が得られる。例えば、ｋ＝０．００５及びｂ＝２×１０^５の場合、
［数６］
Ｅ｛ｔ^２｜β｝＝（１．０１）／（１０^−５＋β^２）
となる。

従って、この適正な事前分布により、ジェフリーズの事前分布へ任意に近づくことができる。

このモデルのアルゴリズムは、

を有する。ここで、期待値は上述の方法で計算される。

別の実施形態では、τ_ｉｇ ^２は、スケールパラメータがｂ＞０であり形状パラメータがｋ＞０である独立なガンマ分布を有する。次式が成り立つことを示すことができる。

（Ａ２２）

ここで、γ_ｉｇ＝β_ｉｇ ^２／２ｂであり、Ｋは変形ベッセル関数を表す。式（Ａ２２）において、ｋ＝１のとき、

であり、式（Ａ２２）において、Ｋ＝０．５のとき、

であり、又はこれと等価であるが、

である。

（Ａ．１）の証明．
条件付き期待値の定義から、γ＝β^２／２ｂと書くと、

が得られる。式変形と、簡単化と、ｕ＝τ^２／ｂの代入とにより、（Ａ２２）における第１の式が得られる。
（２２）における積分は、

という結果を用いて評価されることが可能である。ここで、Ｋは変形ベッセル関数を表す。ワトソン（Watson：１９６６年）を参照されたい。
このクラスの要素の例はｋ＝１であり、この場合は、

である。これは、ティブシラニ（Tibshirani）のラッソ技術（Lasso technique：１９９６年）で使用される事前分布に相当する。フィゲイレド（Figueiredo：２００１年）も参照されたい。
ｋ＝０．５の場合は、

になり、又はこれに等価であるが、

になる。ここで、Ｋ_０及びＫ_１は変形ベッセル関数である。アブラモビッツ及びステガン（１９７０年）を参照されたい。これらのベッセル関数を評価するための多項近似式が、アブラモビッツ及びステガン（１９７０年、３７９ページ）に記述されている。上述の各式は、ラッソモデル及びジェフリーズの事前分布モデルとの関連を実証するものである。

当業者には、ｋがゼロに向かいかつｂが無限大に向かうにつれて、事前分布はジェフリーズの特異事前分布（improper prior）に向かうことが認識されるであろう。

ある実施形態では、０＜ｋ≦１かつｂ＞０である事前分布は、ラッソの事前分布とジェフリーズの超事前分布を用いた仕様との間にあるように、ペナルティー的な非ゼロ係数として解釈される場合もある、事前分布のクラスを形成する。

ハイパーパラメータｂ及びｋは、本方法によって選択される構成要素の数を制御するように変更されることが可能である。ｂが固定されたときにｋがゼロに向かうと、選択される構成要素の数はこれに伴って減少可能であり、逆にｋが１へ向かうと、選択される構成要素の数はこれに伴って増加可能である。

ある好適な実施形態では、ＥＭアルゴリズムは下記のように実行される。

１． ｎ＝０，Ｐ_ｇ＝Ｉを設定し、

の初期値を選ぶ。式（Ａ２２）におけるｂ及びｋの値を選ぶ。例えば、ｂ＝１ｅ７及びｋ＝０は、優れた近似度でジェフリーズの事前分布モデルを与える。これは、ｘ_ｉに関するｌｏｇ（ｐ_ｉｇ／ｐ_ｉＧ）のリッジ（ridge）回帰によって行われる。ここでｐ_ｉｇは、グループｇにおける観測量について１に近い値であるように選択され、そうでなければ、すべての確率の和が１になるという拘束条件のもとで、０より大きい値を有する小さな量であるように選択される。

２． Ｅステップを実行する。すなわち

を評価する。これもまた、ｋ及びｂの値に依存することに留意する。

３． ｔ＝０を設定する。ｇ＝１，…，Ｇ−１について、
ａ） ｐ（ｇ）≧ｎのとき、（Ａ２０）が代入された（Ａ１９）を用いて、δ_ｇ ^ｔ＝γ_ｇ ^ｔ＋１−γ_ｇ ^ｔを計算する。
ｂ） δ^ｔ＝（δ_ｇ ^ｔ，ｇ＝１，…，Ｇ−１）と表すとき、ラインサーチを行って、α^ｔの関数として（１２）を最大化する（又は単に増大させる）

におけるα^ｔの値を見つける。
ｃ）

を設定し、
［数７］
ｔ＝ｔ＋１
を設定する。
ステップ（ａ）及び（ｂ）を収束するまで反復する。
これは、例えばγの関数として流れＱの関数を最大化するγ^＊ｎ＋１を生成する。
ｇ＝１，…，Ｇ−１について、

を決定する。ここで、ε≪１、例えば１０^−５である。Ｐ_ｇを、ｉ∈Ｓ_ｇについてβ_ｉｇ＝０でありかつ

であるように定義する。このステップは、モデルから、小さな係数値を有する変数を除去する。

４． ｎ＝ｎ＋１を設定し、収束するまでは２へ進む。

次に、順序分類ロジスティック回帰に関する第２の実施形態について説明する。

Ｂ．順序付きカテゴリーモデル．
この実施形態に係る方法は、複数のトレーニングサンプルを用いて、試験サンプルがある特定のクラスに属するか否かを決定するために使用可能な複数の構成要素のうちのサブセットを識別する。例えば、マイクロアレイ解析を用いて組織生検サンプルを評価するための遺伝子を識別するためには、昇順又は降順で並んだ病気の重症度のクラスへ予め順序付けられ、例えば正常な組織、良性の組織、局所的な腫瘍及び転移した腫瘍組織として予め順序付けられている組織の一連のサンプルからのマイクロアレイデータが、複数のトレーニングサンプルとして使用され、上記トレーニングサンプルに関連づけられた病気の重症度を示すことができる複数の構成要素のうちのサブセットを識別する。この構成要素のうちのサブセットは、続いて、前もって分類されていない試験サンプルが正常、良性、局所的な腫瘍又は転移した腫瘍として分類されることが可能であるか否かを決定するために使用可能である。従って、複数の構成要素のうちのサブセットは、複数のクラスにてなる順序付けられたセットのうちの特定のクラスに試験サンプルが属するか否かを診断するものとなる。いったん複数の構成要素のうちのサブセットが識別されれば、サンプルがどの順序付けられたクラスに属するかを決定するための将来の診断手順では、上記複数の構成要素のうちのサブセットのみを試験すればよいことは明らかであろう。

本発明の方法は、特に膨大な量のデータの解析に適している。典型的には、試験サンプルから取得される大量のデータ・セットは大幅に変化し、多くの場合、トレーニングサンプルから取得されるものとは著しく異なる。本発明に係る方法は、複数のトレーニングサンプルから生成される膨大な量のデータから複数の構成要素にてなる複数のサブセットを識別することが可能であり、本方法によって識別される複数の構成要素にてなるサブセットは、試験サンプルから生成されるデータが同じクラスに属するトレーニングサンプルから生成されるデータに比べて著しく可変的であっても、次には試験サンプルの分類に使用可能である。従って本発明の方法は、データ品質が悪い場合及び／又は同じ順序のクラスにおけるサンプル間に高いばらつきが存在する場合であってもサンプルを正しく分類する可能性が高い、複数の構成要素のうちのサブセットを識別することができる。

上記構成要素は、その特定の順序のクラスを「予測」する。基本的には、本発明に係る方法は、システムから生成される全てのデータから、トレーニングデータの分類に使用可能な比較的少数の構成要素を識別することを可能にする。いったん本方法によりこれらの構成要素が識別されると、上記構成要素は将来において試験サンプルを分類するために使用可能である。本発明に係る方法は好適には統計的方法を用いて、順序付けされたクラスのうちの要素であるクラスへサンプルを正しく分類するためには必要でない構成要素を除去する。

以下の説明では、Ｎ個のサンプルが存在し、ｙ，ｚ及びμ等のベクトルは、ｉ＝１，…，Ｎについて要素ｙ_ｉ，ｚ_ｉ及びμ_ｉを有する。ベクトルの乗算及び除算は要素に関して定義され、ｄｉａｇ｛・｝は、引数に等しい対角成分を有する対角行列を示す。また、‖・‖はユークリッドノルムを示すために使用される。

好適には、Ｎ個の観測値ｙ_ｉ ^＊が存在する。ここで、ｙ_ｉ ^＊は整数値１，…，Ｇをとる。これらの値は、例えば病気の重症度等の、何らかの方法で順序付けられたクラスを示す。各観測値に付随して、ｎ個の行及びｐ個の列を有する行列Ｘ^＊内に配列される複数の共変量（covariate：変数、例えば遺伝子発現値）のセットが存在する。ここで、ｎはサンプル数であり、ｐは構成要素の個数である。ｘ_ｉ ^＊Ｔという表記は、Ｘ^＊のｉ番目の行を示す。個々の値（サンプル）ｉは、π_ｉｋ＝π_ｋ（ｘ_ｉ ^＊）により与えられるクラスｋに属する確率を有する。

累積確率を定義する。

γ_ｉｋは、単に、ｋ以下のインデックスを有するクラスに観測値ｉが属する確率であることに留意されたい。Ｃを、

で与えられる要素ｃ_ｉｊを備えたＮ×ｐ行列とし、Ｒを、

で与えられる要素γ_ｉｊを備えたｎ×ｐ行列とする。

これらは、行内のＣの列の累積和である。

独立な観測値（サンプル）の場合、データの尤度は、

と書き表すことが可能であり、対数尤度Ｌは、

と書き表すことができる。

これには、ｋ＝２，…，Ｇについて、次のような継続比（continuation ratio：又は逐次ロジット）モデルを採用することができる。

マクロー及びネルダー（McCullagh and Nelder：１９８９年）、マクロー（１９８０年）及びその論考を参照されたい。ここで、次式が成り立つことに注意する。

尤度は、次式の反応ベクトルｙ及び共変量行列Ｘを有するロジスティック回帰の尤度と等価である。

ここで、Ｉ_Ｇ−１は（Ｇ−１）×（Ｇ−１）の恒等行列であり、１_Ｇ−１は１を要素とする（Ｇ−１）×１ベクトルである。ここでｖｅｃ｛ ｝は引数として行列をとり、１行毎にベクトルを形成する。

典型的には、先に論じたように、構成要素に係る複数の重み係数は、大部分の構成要素の重み係数がゼロであるという先験的仮定を考慮した方法で推定される。

フィゲイレド（２００１年）に従って、冗長な変数（共変量）を除去するために、複数のハイパーパラメータにてなるｐ×１ベクトルを導入することによりパラメータβ^＊の事前分布が指定される。

好適には、構成要素の重み係数に関して指定される事前分布は、次式の形式である。

（Ｂ１）

ここで、ｐ（β^＊│ｖ^２）はＮ（０，ｄｉａｇ｛ｖ^２｝）であり、ｐ（ｖ^２）は適切に選ばれた超事前分布である。例えば、

は、適切な形式のジェフリーズの事前分布である。

別の実施形態では、ｐ（ｖ_ｉ ^２）は、ｔ_ｉ ^２＝１／ｖ_ｉ ^２が独立なガンマ分布を有する事前分布である。

別の実施形態では、ｐ（ｖ_ｉ ^２）は、ｖ_ｉ ^２が独立なガンマ分布を有する事前分布である。

シータの要素は、役立つ情報を持たない事前分布を有する。

尤度関数を

と書き表すと、ベイズフレームワークにおいては、ｙを所与とするβ，θ及びｖの事後分布は、次式になる。

［数８］
ｐ（β^＊φｖ│ｙ）αＬ（ｙ│β^＊φ）ｐ（β^＊│ｖ）ｐ（ｖ）
（２）

ｖを，失われたデータのベクトルとして扱うことにより、ＥＭアルゴリズム（デンプスター（Dempster）ほか，１９７７年）等の反復アルゴリズムは、（２）を最大化してβ及びθの最大事後推定値を生成するために使用可能である。上述の事前分布は、最大事後推定値が疎になるようにされる、すなわち多数のパラメータが余分であればβ^＊の多くの要素がゼロになるようにされる。

好適には、下記においてβ^Ｔ＝（θ^Ｔ，β^＊Ｔ）である。

上述の順序付けられたカテゴリーモデルについては、

（１１）

（１２）
であることを証明することができる。ここで、μ_ｉ＝ｅｘｐ（ｘ_ｉ ^Ｔβ）／（１＋ｅｘｐ（ｘ_ｉ ^Ｔβ））及びβ^Ｔ＝（θ_２，…，θ_Ｇ，β^＊Ｔ）である。

複数の構成要素と構成要素に係る複数の重み係数との事後分布を最大化するための反復手順は、例えばデンプスターほか，１９７７年、に記載されているようなＥＭアルゴリズムである。好適には、ＥＭアルゴリズムは下記のように実行される。

１． 超事前分布を選択し、またそのパラメータとして値ｂ及びｋを選択する。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝，φ^（０）及びε＝１０^−５（例えば）を設定する。正則化パラメータκを、１よりずっと大きい値、例えば１００に設定する。これは、以下のＭステップにおいて２階の導関数行列の最初のＧ−１個の対角要素に対して１／κ^２を加算することに相当する。
ｐ≦Ｎであれば、

（Ｂ２）
により初期値β^＊を計算し、ｐ＞Ｎであれば、

（Ｂ３）
により初期値β^＊を計算する。ここで、リッジパラメータλは０＜λ≦１を満足し、ζは小さな値であり、かつζは、リンク関数ｇ（ｚ）＝ｌｏｇ（ｚ／（１−ｚ））がｚ＝ｙ＋ζにおいてうまく定義されているように選ばれる。

２．

を定義し、Ｐ_ｎを、β^（ｎ）の非ゼロ要素γ^（ｎ）が、

を満足するような、ゼロ及び１を要素とする行列であるとする。

であるように
ｗ_β＝（ｗ_βｉ，ｉ＝１，ｐ）を定義し、ｗ_γ＝Ｐ_ｎｗ_βとする。

３．

（１５）
を計算してＥステップを実行する。ここで、Ｌはｙの対数尤度関数であり、

であり、簡単化のために、

であれば

であると定義する。β＝Ｐ_ｎγ及びβ^（ｎ）＝Ｐ_ｎγ^（ｎ）を用いると、（１５）は、

（Ｂ４）
と書くことができる。ここで、β^（ｎ）＝Ｐ_ｎγ^（ｎ）のときにｄ（γ^（ｎ））＝Ｐ_ｎ ^Ｔｄ^（ｎ）と評価される。

４． Ｍステップを実行する。これは、ニュートン＝ラフソン法の反復により次のように実行されることが可能である。γ_０＝γ^（ｎ）を設定し、ｒ＝０，１，２，…についてγ_ｒ＋１＝γ_ｒ＋α_ｒδ_ｒを設定する。ここで、α_ｒは

を保証するようにラインサーチアルゴリズムによって選ばれる。
ｐ≦Ｎの場合、

（Ｂ５）
を用いる。ここで、

及び

である。
ｐ＞Ｎの場合、

（Ｂ６）
を用いる。ここで、Ｖ_ｒ及びｚ_ｒは先に定義した通りである。
γ^＊を、何らかの収束基準、例えば、
［数９］
‖γ_ｒ−γ_ｒ＋１‖＜ε（例えば１０^−５）
が満足されるときのγ_ｒの値であるとする。

５． β^＊＝Ｐ_ｎγ^＊，及び

を定義する。ここで、ε_１は小値の定数、例えば１ｅ−５である。ｎ＝ｎ＋１を設定する。

６． 収束を確認する。ε_２が十分に小さな値であるときに‖γ^＊−γ^（ｎ）‖＜ε_２であれば停止し、そうでなければ上記ステップ２へ進む。

確率の復元．
パラメータβの推定値を取得すると、ｉ＝１，…，Ｎ及びｋ＝２，…，Ｇについて、

を計算する。

好適には、確率を取得するために、帰納法

と、ｉ＝１，…，Ｎについて確率の総和は１になるという事実とを用いる。

ある実施形態では、行ｘ_ｉ ^Ｔを有する共変量行列Ｘは、何らかの核関数κに対してｋ_ｉｊ＝κ（ｘ_ｉ−ｘ_ｊ）のときにｉｊ番目の要素としてｋ_ｉｊを有する行列Ｋで置き換え可能である。この行列は、複数個の１にてなるベクトルによって拡大されることも可能である。以下の表１に核関数のいくつかの例を示す。エブゲニウ（Evgeniou）ほか（１９９９年）を参照されたい。

表１における最後の２つの核関数は、好適には１次元のものである。すなわち、Ｘがただ１つの列を有する場合のものである。これらの核関数の積から多変量バージョンを導出することができる。Ｂ_２ｎ＋１の定義は、デ・ボーア（De Boor：１９７８年）に記載されている。核関数の使用は、共変量Ｘの滑らかな関数（線形変換の場合とは対照的である）である平均値をもたらす。このようなモデルは、データに対して実質的により優れた適合性を与えることができる。

次に、一般化線形モデルに関する第３の実施形態について説明する。

Ｃ．一般化線形モデル．
この実施形態に係る方法は複数のトレーニングサンプルを用いて、サンプルの特性を予測することのできる複数の構成要素のうちのサブセットを識別する。続いて、この構成要素のうちのサブセットについての知識は、関心対象の特性に係る未知の値を予測するための試験、例えば臨床試験に使用可能である。例えば、ＤＮＡマイクロアレイに係る複数の構成要素のうちのサブセットは、例えば血糖レベル、白血球の数、腫瘍の大きさ、腫瘍の成長速度又は生存時間等の、臨床に関連した特性を予測するために使用可能である。

このようにして、本発明は、特定のサンプルの特性を予測するために使用可能な、好適には比較的少数の構成要素を識別する。選択される構成要素は、その特性を「予測」させるものである。超事前分布におけるハイパーパラメータを適正に選ぶことにより、本アルゴリズムは、さまざまなサイズのサブセットを選択できるようになる。本質的に本発明の方法は、システムから生成される全てのデータから、特定の特性を予測するために使用可能な少数の構成要素の識別を可能にする。いったん本方法によりこれらの構成要素が識別されると、上記構成要素は将来において新たなサンプルの特性を予測するために使用可能である。本発明の方法は好適には統計的方法を用いて、上記サンプルの特性を正しく予測するためには不要である構成要素を除去する。

本発明者らは、複数のトレーニングサンプルから生成されるデータに係る複数の構成要素の一次結合に関する構成要素に係る複数の重み係数が、あるトレーニングサンプルの特性を予測するためには不要である構成要素を除去するような方法で推定され得ることを発見している。その結果、トレーニングセットにおける複数の重み係数のサンプルの特性を正しく予測することのできる複数の構成要素のうちのサブセットが識別される。従って本発明の方法は、例えば関心対象の量であるトレーニングサンプルの特性を正しく予測することができる比較的少数の構成要素を、大量のデータから識別することを可能にする。

上記特性は、どのような関心対象の特性であってもよい。ある実施形態では、特性は量又は測定値である。別の実施形態では、これらはあるグループのインデックス番号であってもよく、このとき、複数のサンプルは、予め決められた分類法に基づいて２つのサンプルグループ（又は「クラス」）にグループ化される。この分類法は、複数のトレーニングサンプルがグループ化されるべきときに使用される、所望される任意の分類法であることが可能である。例えば分類法は、トレーニングサンプルが白血病細胞からのものかそれとも健康な細胞からのものかという場合もあれば、トレーニングサンプルが、所定の状態を有する患者又は有していない患者の血液から取得されること、あるいは、トレーニングサンプルが、正常な細胞との比較でいくつかのタイプの癌のうちの１つからの細胞によるものであることという場合もある。別の実施形態では、特性は、特定の患者が少なくとも所定の日数に渡って生存していることを示す打ち切り生存時間である場合もある。別の実施形態では、上記量は、測定可能なサンプルの連続的に可変な任意の特性、例えば血圧であることが可能である。

ある実施形態では、データは、ｉ∈｛１，…，Ｎ｝とするときに、量ｙ_ｉであることが可能である。ここでは、要素ｙ_ｉを備えたｎ×１ベクトルをｙと書き表す。構成要素に係る複数の重み係数（そのうち多数はゼロであると期待される）にてなるｐ×１パラメータベクトルβと、複数のパラメータφ（ゼロであると期待されるわけではない）にてなるｑ×１ベクトルとを定義する。ｑはゼロである可能性もある（すなわち、ゼロであると期待されないパラメータの集合は空である可能性がある）ことに留意されたい。

ある実施形態では、入力データは、ｎ個の試験トレーニングサンプル及びｐ個の構成要素が存在するときのｎ×ｐデータ行列Ｘ＝（ｘ_ｉｊ）に編成される。典型的には、ｐはｎよりずっと大きくなる。

別の実施形態では、データ行列Ｘは、線形予測量ではなく、予測量としてＸの滑らかな関数を得るためにｎ×ｎ核行列Ｋで置換されることが可能である。核行列Ｋの一例は、次式になる。

［数１０］
ｋ_ｉｊ＝ｅｘｐ（−０．５＊（ｘ_ｉ−ｘ_ｊ）^ｔ（ｘ_ｉ−ｘ_ｊ）／σ^２）

ｘの下付き添字は行列Ｘにおける行の番号を示す。理想的には、Ｋの列のうちのサブセットは、これらの滑らかな関数の疎な表現を与えるものが選択される。

典型的には、先に論じたように、構成要素に係る複数の重み係数は、大部分の構成要素の重み係数がゼロであるという先験的仮定を考慮した方法で推定される。

ある実施形態では、構成要素の重み係数に関して指定される事前分布は、次式の形式である。

（Ｃ１）

ここで、ｖは複数のハイパーパラメータにてなるｐ×１ベクトルであり、ｐ（β│ｖ^２）はＮ（０，ｄｉａｇ｛ｖ^２｝）であり、ｐ（ｖ^２）はｖ^２に関する何らかの超事前分布である。

超事前分布の適切な形式は、ジェフリーズの

である。

別の実施形態では、超事前分布ｐ（ｖ^２）は、各ｔ_ｉ ^２＝１／ｖ_ｉ ^２が独立なガンマ分布を有するようなものである。

別の実施形態では、超事前分布ｐ（ｖ^２）は、各ｖ_ｉ ^２が独立なガンマ分布を有するようなものである。

好適には、φに係る情報価値のない事前分布が指定される。

尤度関数は、データ分布のモデルから定義される。好適には、一般に尤度関数は適切な任意の尤度関数である。例えば、尤度関数

は、例えばネルダー及びウェダーバーン（Nelder and Wedderburn：１９７２年）によって記載されているもののような、一般化線形モデル（ＧＬＭ）に適切な形式である可能性があるが、これに制限されるものではない。この場合好適には、尤度関数は、次式の形式である。

（Ｃ２）

ここで、ｙ＝（ｙ_１，…，ｙ_ｎ）^Ｔ及びａ_ｉ（φ）＝φ／ｗ_ｉであり、ｗ_ｉは既知の重み係数にてなる固定されたセットであり、φは単一のスケールパラメータである。

好適には、尤度関数は次のように指定される。次式が所与であるとする。

（Ｃ３）

各観測値は、複数の共変量ｘ_ｉにてなるセットと、線形予測量η_ｉ＝ｘ_ｉ ^Ｔβとを有する。ｉ番目の観測値の平均とその線形予測量との関係は、リンク関数η_ｉ＝ｇ（μ_ｉ）＝ｇ（ｂ’（θ_ｉ））で与えられる。上記リンク関数の逆はｈで表され、すなわち次式になる。

［数１１］
μ_ｉ＝ｂ’（θ_ｉ）＝ｈ（η_ｉ）

スケールパラメータに加えて、一般化線形モデルは、下記の４つの構成要素によって指定されることが可能である。

・尤度関数又は（スケーリングされた）逸脱度関数（deviance function）
・リンク関数
・リンク関数の導関数
・分散関数

一般化線形モデルの共通の例をいくつか次の表に挙げる。

別の実施形態では、リンク関数及び分散関数のみが定義される擬似尤度モデルが指定される。例によっては、このような仕様が上記表内のモデルをもたらす。他の例では、分布は特定されない。

ある実施形態では、ｙを所与とするβ，φ及びｖの事後分布は、次式を用いて推定される。

［数１２］
ｐ（βφｖ│ｙ）αＬ（ｙ│βφ）ｐ（β│ｖ）ｐ（ｖ）
（Ｃ４）

ここで、

は尤度関数である。

ある実施形態では、ｖは失われたデータのベクトルとして扱われることが可能であり、反復手順は、式（Ｃ４）を最大化してβの最大事後推定値を生成するために使用可能である。式（Ｃ１）の事前分布は、最大事後推定値が疎になるようにされる、すなわち多数のパラメータが余分であればβの多くの要素がゼロになるようにされる。

先に述べたように、事後分布を最大化する構成要素の重み係数は、反復手順を用いて決定されることが可能である。好適には、複数の構成要素と構成要素に係る複数の重み係数との事後分布を最大化するための反復手順は、例えばデンプスターほか，１９７７年、に記載されているような、ＥステップとＭステップとを含むＥＭアルゴリズムである。

ＥＭアルゴリズムの実行に際しては、Ｅステップは好適には、次式の形式の項を計算するステップを含む。

（Ｃ４ａ）

ここで、

であり、簡単化のために、

であれば

を定義する。以下、

と記す。同様に、例えばｄ（β^（ｎ））及びｄ（γ^（ｎ））＝Ｐ_ｎ ^Ｔｄ（Ｐ_ｎγ^（ｎ））を定義する。ここで、β^（ｎ）＝Ｐ_ｎγ^（ｎ）であり、Ｐ_ｎはｐ×ｐ恒等行列からβ_ｊ ^（ｎ）＝０である列ｊを削除して得られる。

好適には、Ｐ（β_ｉ│ｖ_ｉ ^２）がＮ（０，ｖ_ｉ ^２）でありかつｐ（ｖ_ｉ ^２）が指定された事前分布を有するとき、式（Ｃ４ａ）はｔ_ｉ ^２＝１／ｖ_ｉ ^２の条件付き期待値を計算することによって演算される。特定の例及び式は、後に提示する。

任意の適切な尤度関数に適する一般的な実施形態では、ＥＭアルゴリズムは下記のようなステップを含む。

（ａ）超事前分布及びそのパラメータの値を選択する。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝を設定してアルゴリズムを初期化し、φ^（０），β^＊を初期化し、εに対して、例えばε＝１０^−５等の値を適用する。

（ｂ）

（Ｃ５）
を定義し、Ｐｎを、β^（ｎ）の非ゼロ要素γ^（ｎ）が、

を満足するような、ゼロ及び１を要素とする行列であるとする。

（ｃ）関数、

（Ｃ６）
を用いて構成要素の重み係数の事後分布の条件付き期待値を計算することにより、推定（Ｅ）ステップを実行する。ここで、Ｌはｙの対数尤度関数である。（Ｃ４ａ）に定義されているようにβ＝Ｐ_ｎγ及びｄ（γ^（ｎ））を用いると、（Ｃ６）は、

（Ｃ７）
と書くことができる。

（ｄ）反復手順を適用してγの関数としてＱを最大化することにより最大化（Ｍ）ステップを実行する。ここで、γ_０＝γ^（ｎ）であり、ｒ＝０，１，２，…に対してγ_ｒ＋１＝γ_ｒ＋α_ｒδｒであり、αｒは、

及び

（Ｃ８）
を保証するようにラインサーチアルゴリズムによって選ばれる。
ここで、（Ｃ４ａ）におけるように、ｄ（γ^（ｎ））＝Ｐ_ｎ ^Ｔｄ（Ｐ_ｎγ^（ｎ））であり、かつβ_ｒ＝Ｐ_ｎγ_ｒに関して、

である。

（ｅ）γ^＊を、何らかの収束基準が満足されるとき、例えば、‖γｒ−γｒ＋１‖＜ε（例えば１０^−５）であるときのγｒの値であるとする。

（ｆ）β^＊＝Ｐ_ｎγ^＊，

を定義する。ここで、ε_１は小値の定数、例えば１ｅ−５である。

（ｇ）ｎ＝ｎ＋１を設定し、φ^{（ｎ＋１）}＝φ^（ｎ）＋κ_ｎ（φ^＊−φ^（ｎ））を選ぶ。ここで、φ^＊は

を満足し、κ_ｎは０＜κ_ｎ≦１であるような減衰係数（damping factor）である。

（ｈ）収束を確認する。ε_２が十分に小さな値であるときに‖γ^＊−γ^（ｎ）‖＜ε_２であれば停止し、そうでなければ上記ステップ（ｂ）へ進む。

別の実施形態では、ｔ_ｉ ^２＝１／ｖ_ｉ ^２は、スケールパラメータがｂ＞０であり形状パラメータがｋ＞０である独立なガンマ分布を有し、そのためｔ_ｉ ^２の密度は次式になる。

次式
［数１３］
Ｅ｛ｔ^２｜β｝＝（２ｋ＋１）／（２／ｂ＋β^２）
が成り立つことは、次のようにして証明することができる。

を定義すると、

になる。

証明．
ｓ＝β^２／２とすると、

になる。ここでｕ＝ｔ^２／ｂを代入すると、

が得られる。次に、ｓ’＝ｂｓとし、γ（ｕ，ｌ，ｋ）の式を代入すると、

になる。結果は、例えばアブラモビッツ及びステガンのラプラス変換表を参照することによって得られる。
条件付き期待値は、
［数１４］
Ｅ｛ｔ^２｜β｝＝Ｉ（１，ｂ，ｋ）／Ｉ（０，ｂ，ｋ）
＝（２ｋ＋１）／（２／ｂ＋β^２）
から得られる。

ｋはゼロへ向かい、ｂは無限大へ向かうとき、ジェフリーズの事前分布を用いる場合と同等の結果が得られる。例えば、ｋ＝０．００５及びｂ＝２×１０^５の場合、
［数１５］
Ｅ｛ｔ^２｜β｝＝（１．０１）／（１０^−５＋β^２）
となる。

従って、この適正な事前分布により、ジェフリーズの超事前分布によるアルゴリズムへ任意に近づくことができる。

別の実施形態では、ｖ_ｉ ^２は、スケールパラメータがｂ＞０であり形状パラメータがｋ＞０である独立なガンマ分布を有する。次式がは証明可能である。

（Ｃ９）

ここでλ_ｉ＝β_ｉ ^２／２ｂであり、Ｋは変形ベッセル関数を示す。これは、次のように証明することができる。

式（Ｃ９）において、ｋ＝１であれば、

である。式（Ｃ９）において、Ｋ＝０．５であれば、

であり、又はこれに等価であるが、

である。

証明．
条件付き期待値の定義から、λ_ｉ＝β_ｉ ^２／２ｂと書くと、

が得られる。式変形と、簡単化と、ｕ＝ｖ_ｉ ^２／ｂの代入とにより、Ａ．１が得られる。
Ａ．１における積分は、

という結果を用いて評価されることが可能である。ここで、Ｋは変形ベッセル関数を表す。ワトソン（１９６６年）を参照されたい。
このクラスの要素の例はｋ＝１であり、この場合は、

である。これは、ラッソ技術、ティブシラニ（１９９６年）で使用される事前分布に相当する。フィゲイレド（２００１年）も参照されたい。
ｋ＝０．５の場合は、

であり、又はこれに等価であるが、

になる。ここで、Ｋ_０及びＫ_１は変形ベッセル関数である。アブラモビッツ及びステガン（１９７０年）を参照されたい。これらのベッセル関数を評価するための多項近似式は、アブラモビッツ及びステガン（１９７０年、３７９ページ）に記載されている。上述の計算の詳細は付録（Appendix）に記されている。

上述の各式は、ラッソモデル及びジェフリーズの事前分布モデルとの関連を実証するものである。

当業者には、ｋがゼロに向かいかつｂが無限大に向かうにつれて、事前分布はジェフリーズの特異事前分布に向かうことが認識されるであろう。

ある実施形態では、０＜ｋ≦１及びｂ＞０である事前分布は、ラッソ事前分布とジェフリーズの超事前分布を用いた元の仕様との間にあるような、ペナルティー的な非ゼロ係数として解釈される場合もある事前分布のクラスを形成する。

別の実施形態では、一般化線形モデルのケースの場合、最大化ステップにおけるステップ（ｄ）は、

をその期待値

で置換することによって推定され得る。これは、データモデルが一般化線形モデルである場合に好適である。

一般化線形モデルでは、期待値

は次のように計算されることが可能である。次式

（Ｃ１０）
から開始する。ここで、Ｘは、ｉ番目の行をｘ_ｉ ^ＴとするＮ×ｐ行列であり、また

（Ｃ１１）
である。このとき、

が得られる。

式（Ｃ１０）及び（Ｃ１１）は、

（Ｃ１２）

（Ｃ１３）
と書くことができる。ここで、

である。

好適には、一般化線形モデルの場合、ＥＭアルゴリズムは下記のステップを含む。

（ａ）超事前分布及びそのパラメータを選ぶ。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝，φ^（０）を設定してアルゴリズムを初期化し、εに対して、例えばε＝１０^−５等の値を適用する。
ｐ≦Ｎであれば、

（Ｃ１４）
によって初期値β^＊を計算し、
ｐ＞Ｎであれば、

（Ｃ１５）
によって初期値β^＊を計算する。ここで、リッジパラメータλは０＜λ≦１を満足し、ζは、小さな値でありかつリンク関数がｙ＋ζでうまく定義されているように選ばれる。

（ｂ）

を定義し、Ｐｎを、β（ｎ）の非ゼロ要素γ（ｎ）が、

を満足するような、ゼロ及び１を要素とする行列であるとする。

（ｃ）関数

（Ｃ１６）
を用いて構成要素の重み係数の事後分布の条件付き期待値を計算することにより、推定（Ｅ）ステップを実行する。ここで、Ｌはｙの対数尤度関数である。β＝Ｐ_ｎγ及びβ^（ｎ）＝Ｐ_ｎγ^（ｎ）を用いると、（Ｃ１６）は、

（Ｃ１７）
と書くことができる。

（ｄ）例えばニュートン＝ラフソン法の反復である反復手順を適用して、γの関数としてＱを最大化することにより最大化（Ｍ）ステップを実行する。ここで、γ_０＝γ^（ｎ）であり、ｒ＝０，１，２，…に対してγ_ｒ＋１＝γ_ｒ＋α_ｒδ_ｒであり、α_ｒは、

であることを保証するようにラインサーチアルゴリズムによって選ばれる。このとき、ｐ≦Ｎについて、

（Ｃ１８）
を使用し、ここで、

であり、下付き添字ｒはこれらの量がμ＝ｈ（ＸＰ_ｎγ_ｒ）で評価されたことを示す。
ｐ＞Ｎに関しては、

（Ｃ１９）
を用いる。ここで、Ｖ_ｒ及びｚ_ｒは先に定義した通りである。

（ｅ）γ^＊を、何らかの収束基準が満足されるとき、例えば‖γ_ｒ−γ_ｒ＋１‖＜ε（例えば１０^−５）であるときのγ_ｒの値であるとする。

（ｆ）β^＊＝Ｐ_ｎγ^＊，

を定義する。ここで、ε_１は小値の定数、例えば１ｅ−５である。ｎ＝ｎ＋１を設定し、φ^ｎ＋１＝φ^ｎ＋κ_ｎ（φ^＊−φ^ｎ）を選ぶ。ここで、φ^＊は

を満足し、κ_ｎは０＜κ_ｎ≦１であるような減衰係数である。ただし、場合によっては、スケールパラメータが既知であることという点に、又はφの更新式を得るためにこの式が明示的に解かれることもある点に注意されたい。

上述の実施形態は、擬似尤度方法を組み込むために拡張されることが可能である（ウェダーバーン（１９７４年）、及びマクロー及びネルダー（１９８３年））。このような実施形態においては、先に詳述したものと同じ反復手順が適切であろうが、Ｌは、先に示したような、また例えばマクロー及びネルダー（１９８３年）の表８．１におけるような擬似尤度によって置換される。ある実施形態では、スケールパラメータφのための変形された更新方法が存在する。これらのモデルを定義するためには、分散関数τ^２、リンク関数ｇ及びリンク関数

の導関数の仕様が必要である。これらが定義されると、上述のアルゴリズムを適用可能である。

擬似尤度モデルの場合の実施形態では、上述のアルゴリズムのステップ５は、

を計算することによってスケールパラメータが更新されるように変形される。ここで、μ及びτはβ^＊＝Ｐ_ｎγ^＊において評価される。好適には、この更新は、モデル内のパラメータ数ｓがＮ未満であれば実行される。Ｎの序数ｓは、ｓがＮよりずっと小さい場合に使用可能である

別の実施形態では、一般化線形モデル及び擬似尤度モデルの双方に関して、行ｘ_ｉ ^Ｔを有する共変量行列Ｘは、何らかの核関数κについてｋ_ｉｊ＝κ（ｘ_ｉ−ｘ_ｊ）とするときに、ｉｊ番目の要素ｋ_ｉｊを有する行列Ｋで置換されることが可能である。またこの行列は、複数の１にてなるベクトルで拡大される場合もある。いくつかの例示的な核関数を以下の表２に示す。エブゲニウほか（１９９９年）を参照されたい。

表２における最後の２つの核関数は、１次元のものである。すなわち、Ｘがただ１つの列を有する場合のものである。これらの核関数の積から多変量バージョンを導出可能である。Ｂ_２ｎ＋１の定義は、デ・ボーア（１９７８年）に記載されている。一般化線形モデル又は擬似尤度モデルの何れにおいても、核関数の使用は共変量Ｘの滑らかな（線形変換の場合とは対照的である）関数である平均値をもたらす。このようなモデルは、データに対して実質的により優れた適合性を与えることができる。

次に、比例ハザードモデルに関する第４の実施形態について説明する。

Ｄ．比例ハザードモデル．
この実施形態に係る方法は、複数のトレーニングサンプルを用いて、定義されたイベント（例えば死亡、回復）が所定の時間期間内に発生する確率に影響を与える可能性のある複数の構成要素のうちのサブセットを識別することができる。トレーニングサンプルはシステムから取得され、トレーニングサンプルの取得からイベント発生までの時間が測定される。イベントまでの時間を複数のトレーニングサンプルから取得されるデータに関連づける統計的方法を用いて、イベントまでの時間の分布を予測することができる複数の構成要素のうちのサブセットが識別され得る。続いてこの構成要素のうちのサブセットについての知識は、例えば臨床試験である試験に使用可能であり、例えば、死亡までの時間又は病気再発までの時間の統計的特徴が予測される。例えば、システムに係る複数の構成要素のうちのサブセットからのデータは、ＤＮＡマイクロアレイから取得される可能性がある。このデータは、例えば、患者の生存時間の期待値又は中央値等の臨床に関連したイベントを予測するために、あるいは所定の症状の発現又は病気の再発を予測するために使用可能である。

このようにして、本発明は、システムのイベントまでの時間の分布を予測するために使用され得る好適には比較的少数の構成要素を識別する。選択される構成要素は、そのイベントまでの時間を「予測」させるものである。本質的に本発明の方法は、システムから生成される全てのデータから、イベントまでの時間を予測するために使用可能な少数の構成要素の識別を可能にする。いったん本方法によりこれらの構成要素が識別されると、上記構成要素は、将来に、新たなサンプルから、システムのイベントまでの時間の統計的特徴を予測するために使用可能である。本発明の方法は好適には、統計的方法を用いて、システムのイベントまでの時間を正しく予測するためには不要である構成要素を除去する。選択されるサブセットのサイズに関する何らかの制御は、モデルにおけるハイパーパラメータを適切に選択することにより達成されることが可能である。

本明細書で使用しているように、「イベントまでの時間」とは、本発明の方法が適用されるサンプルの取得から、イベント発生時刻までの、時間の尺度（単位）を示す。イベントは、観測可能な任意のイベントであることが可能である。システムが生物学的システムである場合、イベントは、例えば、システムに障害が生じるまでの時間、死亡までの時間、特定の１つ又は複数の症状の発現、状態又は病気の発現又は再発、表現型又は遺伝子型の変化、生化学上の変化、有機体又は組織の形態変化、行動の変化、である可能性がある。

サンプルは、先行する１つのイベントまでの複数の時間からの、特定の１つのイベントまでの時間に関連づけられる。イベントまでの時間は、例えば、サンプリングから死亡までの時間が分かっている患者から取得されたデータから決定される時間、言い替えれば「真正の」生存時間であってもよく、また、最後にサンプルが取得された時点では患者は生存していたという情報のみのみを有する患者から取得されたデータから決定される時間、言い替えれば、その特定の患者は少なくとも所定の日数は生存していたことを示す「打ち切りされた」生存時間であってもよい。

ある実施形態では、入力データは、ｎ個の試験トレーニングサンプル及びｐ個の構成要素が存在するときのｎ×ｐデータ行列Ｘ＝（ｘ_ｉｊ）に編成される。典型的には、ｐはｎよりずっと大きくなる。

例として、例えば、Ｎ個の個体（又はサンプル）とその各個体対するｐ個の遺伝子と
が存在するときのマイクロアレイ実験からの、Ｎ×ｐデータ行列Ｘ＝（ｘ_ｉｊ）について考察する。好適には、各個体ｉ（ｉ＝１，２，…，Ｎ）に関連付けられた変数であって、例えば生存時間であるイベントまでの時間を示す変数ｙ_ｉ（ｙ_ｉ≧０）が存在する。また、各個体に関して、その個体の生存時間が真正の生存時間であるかそれとも打ち切りされた生存時間であるかを示す変数が定義されていてもよい。打ち切り指示子をｃ_ｉと表示する。ここで、次式のように定義する。

生存時間ｙ_ｉを備えたＮ×１ベクトルは

と表記され、打ち切り指示子ｃ_ｉを備えたＮ×１ベクトルは

と表記されることが可能である。

典型的には、先に論じたように、構成要素の重み係数は、大部分の構成要素の重み係数がゼロであるという先験的仮定を考慮した方法で推定される。

好適には、構成要素の重み係数に関して指定される事前分布は、次式の形式である。

（Ｄ１）

ここで、β_１，β_２，…，β_ｎは構成要素の重み係数であり、ｐ（β_１│τ_ｊ）はＮ（０，τ_ｉ ^２）であり、ｐ（τ_ｉ）は、ジェフリーズの超事前分布ではない、何らかの超事前分布

である。

ある実施形態では、事前分布はτの逆ガンマ事前分布であり、上記分布におけるｔ_ｉ ^２＝１／τ_ｉ ^２は、スケールパラメータがｂ＞０であり形状パラメータがｋ＞０であるような独立なガンマ分布を有し、そのためｔ_ｉ ^２の密度は次式になる。

次式が成り立つことを証明することができる。

［数１６］
Ｅ｛ｔ^２｜β｝＝（２ｋ＋１）／（２／ｂ＋β^２）
（Ａ）

式Ａは、次のようにして証明することができる。

を定義すると、

になる。

証明．
ｓ＝β^２／２とすると、

になる。ここでｕ＝ｔ^２／ｂを代入すると、

が得られる。次に、ｓ’＝ｂｓとし、γ（ｕ，ｌ，ｋ）の式を代入すると、

になる。結果は、例えばアブラモビッツ及びステガンのラプラス変換表を参照することによって得られる。
条件付き期待値は、
［数１７］
Ｅ｛ｔ^２｜β｝＝Ｉ（１，ｂ，ｋ）／Ｉ（０，ｂ，ｋ）
＝（２ｋ＋１）／（２／ｂ＋β^２）
から得られる。

ｋはゼロへ向かい、ｂは無限大へ向かうとき、ジェフリーズの事前分布を用いる場合と同等の結果が得られる。例えば、ｋ＝０．００５及びｂ＝２×１０^５の場合、
［数１８］
Ｅ｛ｔ^２｜β｝＝（１．０１）／（１０^−５＋β^２）
となる。

従って、この適正な事前分布により、ジェフリーズの超事前分布へと任意に近づくことができる。

このモデルの変更されたアルゴリズムは、
［数１９］
ｂ^（ｎ）＝Ｅ｛ｔ^２｜β^（ｎ）｝^−０．５
を有する。ここで、期待値は上述の方法で計算される。

さらに別の実施形態では、事前分布はτ_ｉｇ ^２のガンマ分布である。好適には、上記ガンマ分布はスケールパラメータｂ＞０及び形状パラメータｋ＞０を有する。

次式が成り立つことを証明することができる。

ここで、γ_ｉ＝β_ｉ ^２／２ｂであり、Ｋは変形ベッセル関数を表す。このクラスのいくつかの特別な要素はｋ＝１であり、この場合、

である。これは、ラッソ技術、ティブシラニ（１９９６年）で使用される事前分布に相当する。フィゲイレド（２００１年）も参照されたい。

ｋ＝０．５の場合は、

であり、又はこれに等価であるが、

になる。ここで、Ｋ_０及びＫ_１は変形ベッセル関数である。アブラモビッツ及びステガン（１９７０年）を参照されたい。これらのベッセル関数を評価するための多項近似式は、アブラモビッツ及びステガン（１９７０年，３７９ページ）に記載されている。

上述の各式は、ラッソモデル及びジェフリーズの事前分布モデルとの関連を実証するものである。

上述の計算の詳細は、次のようなものである。

上述のガンマ事前分布及びγ_ｉ＝β_ｉ ^２／２ｂの場合、

（Ｄ２）
である。ここで、Ｋは変形ベッセル関数を表す。
（Ｄ２）において、ｋ＝２の場合、

である。
（Ｄ２）において、Ｋ＝０．５の場合、

であり、又はこれに等価であるが、

である。

証明．
条件付き期待値の定義から、γ_ｉ＝β_ｉ ^２／２ｂと書くと、

が得られる。式変形と、簡単化と、ｕ＝ｖ_ｉ ^２／ｂの代入とにより、Ａ．１が得られる。
Ａ．１における積分は、

という結果を用いて評価されることが可能である。ここで、Ｋは変形ベッセル関数を表す。ワトソン（１９６６年）を参照されたい。

ある特に好適な実施形態では、ｐ（τ_ｉ）α１／τ_ｉ ^２は、コッツ及びジョンソン（１９８３年）における、ジェフリーズの事前分布である。

尤度関数は、データの分布に基づいてそのデータに適合するモデルを定義する。好適には、尤度関数は、次式の形式である。

ここで、

及び

はモデルパラメータである。上記尤度関数によって定義されるモデルは、システムのイベントまでの時間を予測するための任意のモデルであることが可能である。

ある実施形態では、尤度によって定義されるモデルはコックスの比例ハザードモデルである。コックスの比例ハザードモデルはコックスによって導入され（１９７２年）、好適には生存データの回帰モデルとして使用可能である。コックスの比例ハザードモデルでは、

は、複数の構成要素に関連づけられる（説明的な）複数のパラメータにてなるベクトルである。好適には、本発明の方法は、データＸ，

及び

を所与とするコックスの比例ハザードモデルのパラメータ

からの、節約志向（parsimonious）の選択（及び推定）を提供する。

コックスの比例ハザードモデルの適用は、同じ生存時間に関して、言い替えれば拘束された生存時間（tied survival time）に関してシステムから異なるデータが取得されるような状況においては問題がある可能性がある。よって拘束された生存時間に対しては、一意的な生存時間をもたらす前置処理ステップが実行されてもよい。提案される前置処理は後続のコードを簡単化し、よってコックスの比例ハザードモデルを続いて適用する際における拘束された生存時間に関する懸念を回避する。

生存時間の前置処理は、極めて少量の微小なランダムノイズを付加することによって行われる。好適には、本手順は、拘束された時間にてなる複数のセットを使用し、かつ、拘束された時間にてなるあるセット内の各拘束された時間に対して、ゼロ平均と、ソーティングされた生存時間の間における非ゼロの最小距離に比例した分散とを有する正規分布から引き出されたランダム量を付加する。このような前置処理は、生存時間の過酷な摂動をもたらすことなく、拘束された時間の除去を達成する。

前置処理は、明確に区別される別個の生存時間をもたらす。好適には、これらの時間は、

で表されるように、大きさに関して昇順で順序付けられることが可能である。

Ｚの行の順序づけが、

の順序づけにより導出される順序づけに対応している場合の、Ｘの行の並べ替えであるＮ×ｐ行列をＺで示し、また行列Ｚのｊ番目の行をＺ_ｊで示す。ｄを、

の順序づけに必要とされるものと同じ順列を用いてｃを順序づけした結果であるとする。

拘束された生存時間の前置処理が考慮されかつ生存データ解析に関する標準的文書（例えばコックス及びオークス（Cox and Oakes：１９８４年））が参照された後には、比例ハザードモデルの尤度関数は、好適には、次式で表すことができる。

（Ｄ３）

ここで、

であり、ｚ_ｊはＺのｊ番目の行であり、

は、ｊ番目の順序を有するイベント時刻ｔ_（ｊ）において設定されるリスクである。

尤度の対数（すなわちＬ＝ｌｏｇ（ｌ））は、好適には、次式で表すことができる。

（Ｄ４）

ここで、

である。

モデルはノンパラメトリックであり、ここで、生存分布のパラメトリック形式は指定されず、好適には（リスクセットの決定において）生存時間の順序付けに係る特性のみが使用されるということに留意されたい。これはノンパラメトリックなケースであるので、

は不要である（すなわちｑ＝０）。

本発明の方法の別の実施形態では、尤度関数によって定義されるモデルはパラメトリック生存モデルである。好適には、パラメトリック生存モデルにおいて、

は、複数の構成要素に関連づけられる複数の（説明的）パラメータにてなるベクトルであり、

は、生存密度関数の関数形式に関連づけられる複数のパラメータにてなるベクトルである。

好適には、本発明の方法は、データＸ，

及び

を所与とするときのパラメトリック生存モデルに関する、パラメータ

と、

）の推定とからの、節約志向の選択（及び推定）を提供する。

パラメトリック生存モデルの適用において、生存時間は前置処理を必要とせず、

で示される。パラメトリック生存モデルは、次のように適用される。

生存時間のパラメトリック密度関数を

で表し、その生存関数を

で表す。ここで、

は密度関数のパラメトリック形式に関連するパラメータであり、

，Ｘは先に定義した通りである。ハザード関数は、

と定義される。

好適には、打ち切りされたデータを考慮した対数尤度関数の一般的な定式化は、次式になる。

パラメトリック回帰生存モデルを用いた生存時間データの解析に関する標準的文書を参照すると、使用可能な多数の生存時間分布が存在することがわかる。使用可能な生存分布には、例えばワイブル分布、指数分布又は極値分布が含まれる。

ハザード関数を、

と書くことができれば、

及び

となる。ここで、

は積分されたハザード関数であり、

であり、Ｘ_ｉはＸのｉ番目の行である。

ワイブル分布、指数分布又は極値分布は、直前の段落に提示した形式で書き表すことのできる密度及びハザード関数を有する。

その適用についての詳細は、部分的には、エイトケン及びクレイトン（Aitken and Clayton：１９８０年）のアルゴリズムに依存するが、ユーザは基本となるパラメトリックハザード関数を任意に指定することができる。

エイトケン及びクレイトン（１９８０年）によると、パラメトリック生存モデルをモデリングする好適な尤度関数は、次式になる。

（Ｄ５）

ここで、

である。エイトケン及びクレイトン（１９８０年）は、式（１１）の結果として、ｃ_ｉは平均値μ_ｉを有するポワソン変量として扱われることが可能であり、式（１１）の最後の項は

に依存しない（ただし

に依存する）と述べている。

好適には、

を所与とする

の事後分布は、次式になる。

（Ｄ６）

ここで、

は尤度関数である。

ある実施形態では、

は、失われたデータのベクトルとして扱われることが可能であり、反復手順は、式（Ｄ６）を最大化して

の事後推定値を生成するために使用可能である。式（Ｄ１）の事前分布は、最大事後推定値が疎になるような、すなわち多数のパラメータが余分であれば、

の多くの要素がゼロになるようなものである。

の多くの要素はゼロである、という事前の期待が存在するので、推定は、推定されるβ_ｉの大部分がゼロであり、その他の非ゼロ推定値が生存時間についての適切な説明となるように実行されることが可能である。

マイクロアレイデータのコンテキストにおいては、この実行は、イベント時間に関して適切な説明をもたらす複数の遺伝子にてなる節約志向のセットを識別することへと移行する。

先に述べたように、事後分布を最大化する構成要素の重み係数は、反復手順を用いて決定されることが可能である。好適には、複数の構成要素と構成要素に係る複数の重み係数との事後分布を最大化するための反復手順は、例えばデンプスターほか，１９７７年に記載されているようなＥＭアルゴリズムである。

ベータを含まない項を無視した（Ｄ６）から、ＥＭアルゴリズムのＥステップを調べる場合、次式を計算する必要がある。

（Ｄ７）

ここで、

であり、簡単化のために

であれば

であると定義する。以下、

と記す。同様に、例えばｄ（β^（ｎ））及びｄ（γ^（ｎ））＝Ｐ_ｎ ^Ｔｄ（Ｐ_ｎγ^（ｎ））を定義する。ここで、β^（ｎ）＝Ｐ_ｎγ^（ｎ）であり、Ｐ_ｎは、ｐ×ｐ恒等行列から、β_ｊ ^（ｎ）＝０である列ｊを除去して得られる。

従って、Ｅステップを実行するためには、ｐ（β_ｉ│τ_ｉ ^２）がＮ（０，τ_ｉ ^２）でありかつｐ（τ_ｉ ^２）が先に論じたような指定された事前分布を有するときの、ｔ_ｉ ^２＝１／τ_ｉ ^２の条件付き期待値を計算する必要がある。

ある実施形態では、ＥＭアルゴリズムは下記のようなステップを含む。

１． 超事前分布及びそのパラメータの値、すなわちｂ及びｋを選ぶ。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝を設定してアルゴリズムを初期化し、

を初期化する。

２．

を定義し、Ｐ_ｎを、

の非ゼロ要素

が、

（Ｄ８）
を満足するような、ゼロ及び１を要素とする行列であるとする。

３． 構成要素の重み係数の事後分布の期待値を計算することにより、推定ステップを実行する。これは、関数、

（Ｄ９）
を用いて実行されることが可能である。ここで、Ｌは

の対数尤度関数である。β＝Ｐ_ｎγ及びβ^（ｎ）＝Ｐ_ｎγ^（ｎ）を用いると、

（Ｄ１０）
が得られる。

４． 最大化ステップを実行する。これは、ニュートン＝ラフソン法の反復を用いて下記のように実行されることが可能である。

を設定し、ｒ＝０，１，２，…に関して

とする。ここで、α_ｒは、

を保証するようにラインサーチアルゴリズムによって選ばれ、また、

（Ｄ１１）
である。ここで、

の場合、

である。

を、何らかの収束基準が満足されるとき、例えば、

（例えばε＝１０^−５）のときの

の値であるとする。

５．

を定義する。ここで、ε_１は小値の定数、例えば１０^−５である。ｎ＝ｎ＋１を設定し、

を選ぶ。ここで、

は

を満足し、κ_ｎは０＜κ_ｎ＜１であるような減衰係数である。

６． 収束を確認する。ε_２が十分に小さな値であるときに

であれば停止し、そうでなければ上記ステップ２へ進む。

別の実施形態では、最大化ステップにおけるステップ（Ｄ１１）は、

をその期待値

で置き換えることによって推定されてもよい。

ある実施形態では、ＥＭアルゴリズムは、モデルがコックスの比例ハザードモデルである場合に、事後分布を最大化するために適用される。

モデルがコックスの比例ハザードモデルである場合におけるＥＭアルゴリズムの適用についての説明の助けとするために、「動的重み係数」及びこれらの重み係数に基づく行列を定義することが好適である。上記重み係数は、次式になる。

これらの重み係数に基づく行列は、次式になる。

重み係数行列の観点から、Ｌの１階及び２階の導関数は、次式のように書き表すことができる。

（Ｄ１２）

ここで、Ｋ＝Ｗ−Δ（Ｗ）である。従って、ＥＭアルゴリズムのステップ（２）の一部で説明した変換行列Ｐ_ｎ（式Ｄ８）から、次式が得られることに留意されたい（式Ｄ１１も参照）。

（Ｄ１３）

好適には、上記モデルがコックスの比例ハザードモデルである場合、ＥＭアルゴリズムのＥステップ及びＭステップは下記のようなものになる。

１． 超事前分布及びそのパラメータｂ及びｋを選ぶ。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝を設定する。ｖを、何らかの小さな値ε、例えば．００１に関して、要素

を備えたベクトルであるとする。ｆを、ｌｏｇ（ｖ／ｔ）であると定義する。
ｐ≦Ｎであれば、

により初期値

を計算する。
ｐ＞Ｎであれば、

により初期値

を計算する。ここで、リッジパラメータλは０＜λ≦１を満足する。

２．

を定義する。Ｐ_ｎを、

の非ゼロ要素

が、

を満足するような、ゼロ及び１を要素とする行列であるとする。

３．

を計算してＥステップを実行する。ここで、Ｌは式（８）によって与えられる

の対数尤度関数である。β＝Ｐ_ｎγ及びβ^（ｎ）＝Ｐ_ｎγ^（ｎ）を用いると、

が得られる。

４． Ｍステップを実行する。これは、ニュートン＝ラフソン法の反復を用いて下記のように実行されることが可能である。

を設定し、ｒ＝０，１，２，…に関して

とする。ここで、α_ｒは、

を保証するようにラインサーチアルゴリズムによって選ばれる。
ｐ≦Ｎである場合、

を用いる。ここで、

である。
ｐ＞Ｎである場合、

を用いる。
γ^＊を、何らかの収束基準が満足されるとき、例えば‖γ_ｒ−γ_ｒ＋１‖＜ε（例えば１０^−５）のときにおけるγ_ｒの値であるとする。

５．

を定義する。ここで、ε_１は小値の定数、例えば１０^−５である。このステップは、非常に小さな係数を有する変数を除去する。

６． 収束を確認する。ε_２が十分に小さな値であるとき

であれば停止し、そうでなければｎ＝ｎ＋１を設定して上記ステップ２へ進み、収束が起こるまで手順を反復する。

別の実施形態では、ＥＭアルゴリズムは、上記モデルがパラメトリック生存モデルである場合に、事後分布を最大化すべく適用される。

パラメトリック生存モデルにＥＭアルゴリズムを適用する際には、式（１１）の結果として、ｃ_ｉは、平均値μ_ｉｉを有するポワソン変量として扱われることが可能であり、式（１１）の最後の項はβに依存しない（ただしφに依存する）。

であり、よって問題点をポワソン型の平均値（Poisson-like mean）の対数線形モデルの形で表現することが可能である点に留意されたい。好適には、対数尤度関数の反復的最大化は、

の初期推定値が与えられたときに

の推定値が取得される場合に実行される。次に、

のこれらの推定値を所与として、

の更新された推定値が取得される。本手順は、収束が起きるまで継続される。

先に述べた事後分布の適用に際しては、（固定された

に関して）

（Ｄ１４）
に留意する。

結果的に、式（１１）及び（１２）から、

及び

が得られる。

式（１２）の、パラメトリック生存モデルに関連したバージョンは、次式になる。

（Ｄ１５）

ＥＭアルゴリズムの各Ｍステップの後に

について解くために（下記のステップ５を参照）、好適には、

とする。ここで、０＜κ_ｎ≦１である場合、

は

を満足し、βは以前のＭステップから取得された値に固定される。

パラメータ選択のためのＥＭアルゴリズムを、パラメトリック生存モデル及びマイクロアレイデータのコンテキストにおいて提供することが可能である。好適には、上記ＥＭアルゴリズムは下記の通りである。

１． 超事前分布及びそのパラメータｂ及びｋを選択し、例えばｂ＝１ｅ７及びｋ＝０．５とする。ｎ＝０，Ｓ_０＝｛１，２，…，ｐ｝，

を設定する。ｖを、何らかの小値ε、例えば、．００１に関して、要素

を備えたベクトルであるとする。ｆを、ｌｏｇ（ｖ／Λ（ｙ，φ））であると定義する。
ｐ≦Ｎであれば、

により初期値

を計算する。
ｐ＞Ｎであれば、

により初期値

を計算する。
ここで、リッジパラメータλは０＜λ≦１を満足する。

２．

を定義する。Ｐ_ｎを、

の非ゼロ要素

が、

を満足するような、ゼロ及び１を要素とする行列であるとする。

３．

を計算してＥステップを実行する。ここで、Ｌは

及び

の対数尤度関数である。
β＝Ｐ_ｎγ及びβ^（ｎ）＝Ｐ_ｎγ^（ｎ）を用いると、

が得られる。

４． Ｍステップを実行する。これは、ニュートン＝ラフソン法の反復を用いて下記のように実行されることが可能である。

を設定し、ｒ＝０，１，２，…に関して

とする。ここで、α_ｒは、

を保証するようにラインサーチアルゴリズムによって選ばれる。
ｐ≦Ｎである場合、

を用いる。ここで、

である。
ｐ＞Ｎである場合、

を用いる。
γ^＊を、何らかの収束基準が満足されるとき、例えば‖γ_ｒ−γ_ｒ＋１‖＜ε（例えば１０^−５）のときにおけるγ_ｒの値であるとする。

５．

を定義する。ここで、ε_１は小値の定数、例えば１０^−５である。ｎ＝ｎ＋１を設定し、

を選ぶ。ここで、

は

を満足し、κ_ｎは０＜κ_ｎ＜１であるような減衰係数である。

６． 収束を確認する。ε_２が十分に小さな値であるとき

であれば停止し、そうでなければステップ２へ進む。

別の実施形態では、生存時間はワイブル生存密度関数によって記述される。ワイブルのケースでは、

は好適には１次元であり、かつ、

である。

好適には、αの更新された値を供給するために、各Ｍステップの後に

が解かれる。

コックスの比例ハザードモデルに関して適用されるステップによれば、αを推定することができ、また、生存時間がワイブル分布に従うならば生存時間に関する適切な説明を提供することのできる複数のパラメータにてなる節約志向のサブセットを、

から選択することができる。次に、数値的な例を挙げる。

ここで、以下の限定的ではない例のみを参照して本発明の好適な実施形態について説明する。ただし、以下の例は単に例示的なものであり、いかなる点でも、以上説明した本発明の一般性を限定するものとして解釈されるべきでないことは理解される必要がある。

２０１個のデータポイントと４１個の基礎関数（basis function）に係る全正規回帰（Full normal regression）の例．

ｋ＝０及びｂ＝１ｅ７
正しい４つの基礎関数が以下に識別される。
２ １２ ２４ ３４
推定された分散は０．６７である。

ｋ＝０．２及びｂ＝１ｅ７の場合。
８つの基礎関数が以下に識別される。
２ ８ １２ １６ １９ ２４ ３４
推定された分散は０．６３である。基礎関数の正しいセットはこのセット内に含まれることに留意されたい。

ｋ＝０．２及びｂ＝１ｅ７に対する反復の結果を以下に示す。

［表１］
――――――――――――――――――――――――――――――
EM Iteration: 0 expected post: 2 basis fns 41

sigma squared 0.6004567
EM Iteration: 1 expected post: -63.91024 basis fns 41

sigma squared 0.6037467
EM Iteration: 2 expected post: -52.76575 basis fns 41

sigma squared 0.6081233
EM Iteration: 3 expected post: -53.10084 basis fns 30

sigma squared 0.6118665
EM Iteration: 4 expected post: -53.55141 basis fns 22

sigma squared 0.6143482
EM Iteration: 5 expected post: -53.79887 basis fns 18

sigma squared 0.6155
EM Iteration: 6 expected post: -53.91096 basis fns 18

sigma squared 0.6159484
EM Iteration: 7 expected post: -53.94735 basis fns 16

sigma squared 0.6160951
EM Iteration: 8 expected post: -53.92469 basis fns 14

sigma squared 0.615873
EM Iteration: 9 expected post: -53.83668 basis fns 13

sigma squared 0.6156233
EM Iteration: 10 expected post: -53.71836 basis fns 13

sigma squared 0.6156616
EM Iteration: 11 expected post: -53.61035 basis fns 12

sigma squared 0.6157966
EM Iteration: 12 expected post: -53.52386 basis fns 12

sigma squared 0.6159524
EM Iteration: 13 expected post: -53.47354 basis fns 12

sigma squared 0.6163736
EM Iteration: 14 expected post: -53.47986 basis fns 12

sigma squared 0.6171314
EM Iteration: 15 expected post: -53.53784 basis fns 11

sigma squared 0.6182353
EM Iteration: 16 expected post: -53.63423 basis fns 11

sigma squared 0.6196385
EM Iteration: 17 expected post: -53.75112 basis fns 11

sigma squared 0.621111
EM Iteration: 18 expected post: -53.86309 basis fns 11

sigma squared 0.6224584
EM Iteration: 19 expected post: -53.96314 basis fns 11

sigma squared 0.6236203
EM Iteration: 20 expected post: -54.05662 basis fns 11

sigma squared 0.6245656
EM Iteration: 21 expected post: -54.1382 basis fns 10

sigma squared 0.6254182
EM Iteration: 22 expected post: -54.21169 basis fns 10

sigma squared 0.6259266
EM Iteration: 23 expected post: -54.25395 basis fns 9

sigma squared 0.6259266
EM Iteration: 24 expected post: -54.26136 basis fns 9

sigma squared 0.6260238
EM Iteration: 25 expected post: -54.25962 basis fns 9

sigma squared 0.6260203
EM Iteration: 26 expected post: -54.25875 basis fns 8

sigma squared 0.6260179
EM Iteration: 27 expected post: -54.25836 basis fns 8

sigma squared 0.626017
EM Iteration: 28 expected post: -54.2582 basis fns 8

sigma squared 0.6260166
――――――――――――――――――――――――――――――

２０１個の観測値と１０個の変数とを備えた削減されたデータセットに対して、ｋ＝０及びｂ＝１ｅ７。
正しい基礎関数、すなわち１ ２ ３ ４を与える。ｋ＝０．５及びｂ＝１ｅ７のとき、７個の基礎関数、すなわち１ ２ ３ ４ ６ ８ ９が選択される。反復の記録を以下に示す。このセットもまた正しいセットを含むということに留意されたい。

［表２］
――――――――――――――――――――――――――――――
EM Iteration: 0 expected post: 2 basis fns 10

sigma squared 0.6511711
EM Iteration: 1 expected post: -66.18153 basis fns 10

sigma squared 0.6516289
EM Iteration: 2 expected post: -57.69118 basis fns 10

sigma squared 0.6518373
EM Iteration: 3 expected post: -57.72295 basis fns 9

sigma squared 0.6518373
EM Iteration: 4 expected post: -57.74616 basis fns 8

sigma squared 0.65188
EM Iteration: 5 expected post: -57.75293 basis fns 7

sigma squared 0.6518781
――――――――――――――――――――――――――――――

順序付きカテゴリーの例．
１５個のサンプル及び９６０５個の遺伝子に係るルオ（Luo）の前立腺データ。ｋ＝０及びｂ＝１ｅ７に対して、以下の結果を得る。

［表３］
――――――――――――――――――――――――――――――
misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 2 1 0
3 0 0 4 0
4 0 0 0 4

Identifiers of variables left in ordered categories model
6611
――――――――――――――――――――――――――――――

ｋ＝０．２５及びｂ＝１ｅ７に対して、以下の結果を得る。

［表４］
――――――――――――――――――――――――――――――
misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 3 0 0
3 0 0 4 0
4 0 0 0 4

Identifiers of variables left in ordered categories model
6611 7188
――――――――――――――――――――――――――――――

ここでは、余分なデータの付加により、トレーニングデータは完全に区別されていることに留意されたい。アルゴリズムの反復の記録を以下に示す。

［表５］
――――――――――――――――――――――――――――――
***********************************************
Iteration 1 : 11 cycles, criterion -4.661957

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 9608
――――――――――――――――――――――――――――――

［表６］
――――――――――――――――――――――――――――――
***********************************************
Iteration 2 : 5 cycles, criterion -9.536942

misclassification matrix
fhat
f 1 2
1 23 0
2 1 21
row =true class

Class 1 Number of basis functions in model : 6431
――――――――――――――――――――――――――――――

［表７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 3 : 4 cycles, criterion -9.007843

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 508
――――――――――――――――――――――――――――――

［表８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 4 : 5 cycles, criterion -6.47555

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 62
――――――――――――――――――――――――――――――

［表９］
――――――――――――――――――――――――――――――
***********************************************
Iteration 5 : 6 cycles, criterion -4.126996

misclassification matrix
fhat
f 1 2
1 23 0
2 1 21
row =true class

Class 1 Number of basis functions in model : 20
――――――――――――――――――――――――――――――

［表１０］
――――――――――――――――――――――――――――――
***********************************************
Iteration 6 : 6 cycles, criterion -3.047699

misclassification matrix
fhat
f 1 2
1 23 0
2 1 21
row =true class

Class 1 Number of basis functions in model : 12
――――――――――――――――――――――――――――――

［表１１］
――――――――――――――――――――――――――――――
***********************************************
Iteration 7 : 5 cycles, criterion -2.610974

misclassification matrix
fhat
f 1 2
1 23 0
2 1 21
row =true class

Class 1 : Variables left in model
1 2 3 408 846 6614 7191 8077
regression coefficients
28.81413 14.27784 7.025863 -1.086501e-06 4.553004e-09 -16.25844 0.1412991 -0.04101412

――――――――――――――――――――――――――――――

［表１２］
――――――――――――――――――――――――――――――
***********************************************
Iteration 8 : 5 cycles, criterion -2.307441

misclassification matrix
fhat
f 1 2
1 23 0
2 1 21
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191 8077
regression coefficients
32.66699 15.80614 7.86011 -18.53527 0.1808061 -0.006728619

――――――――――――――――――――――――――――――

［表１３］
――――――――――――――――――――――――――――――
***********************************************
Iteration 9 : 5 cycles, criterion -2.028043

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191 8077
regression coefficients
36.11990 17.21351 8.599812 -20.52450 0.2232955 -0.0001630341

――――――――――――――――――――――――――――――

［表１４］
――――――――――――――――――――――――――――――
***********************************************
Iteration 10 : 6 cycles, criterion -1.808861

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191 8077
regression coefficients
39.29053 18.55341 9.292612 -22.33653 0.260273 -8.696388e-08

――――――――――――――――――――――――――――――

［表１５］
――――――――――――――――――――――――――――――
***********************************************
Iteration 11 : 6 cycles, criterion -1.656129

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
42.01569 19.73626 9.90312 -23.89147 0.2882204

――――――――――――――――――――――――――――――

［表１６］
――――――――――――――――――――――――――――――
***********************************************
Iteration 12 : 6 cycles, criterion -1.554494

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
44.19405 20.69926 10.40117 -25.1328 0.3077712
――――――――――――――――――――――――――――――

［表１７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 13 : 6 cycles, criterion -1.487778

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
45.84032 21.43537 10.78268 -26.07003 0.3209974

――――――――――――――――――――――――――――――

［表１８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 14 : 6 cycles, criterion -1.443949

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
47.03702 21.97416 11.06231 -26.75088 0.3298526

――――――――――――――――――――――――――――――

［表１９］
――――――――――――――――――――――――――――――
***********************************************
Iteration 15 : 6 cycles, criterion -1.415

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
47.88472 22.35743 11.26136 -27.23297 0.3357765

――――――――――――――――――――――――――――――

［表２０］
――――――――――――――――――――――――――――――
***********************************************
Iteration 16 : 6 cycles, criterion -1.395770

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
48.47516 22.62508 11.40040 -27.56866 0.3397475

――――――――――――――――――――――――――――――

［表２１］
――――――――――――――――――――――――――――――
***********************************************
Iteration 17 : 5 cycles, criterion -1.382936

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
48.88196 22.80978 11.49636 -27.79991 0.3424153

――――――――――――――――――――――――――――――

［表２２］
――――――――――――――――――――――――――――――
***********************************************
Iteration 18 : 5 cycles, criterion -1.374340

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.16029 22.93629 11.56209 -27.95811 0.3442109

――――――――――――――――――――――――――――――

［表２３］
――――――――――――――――――――――――――――――
***********************************************
Iteration 19 : 5 cycles, criterion -1.368567

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.34987 23.02251 11.60689 -28.06586 0.3454208

――――――――――――――――――――――――――――――

［表２４］
――――――――――――――――――――――――――――――
***********************************************
Iteration 20 : 5 cycles, criterion -1.364684

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.47861 23.08109 11.63732 -28.13903 0.3462368

――――――――――――――――――――――――――――――

［表２５］
――――――――――――――――――――――――――――――
***********************************************
Iteration 21 : 5 cycles, criterion -1.362068

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.56588 23.12080 11.65796 -28.18862 0.3467873

――――――――――――――――――――――――――――――

［表２６］
――――――――――――――――――――――――――――――
***********************************************
Iteration 22 : 5 cycles, criterion -1.360305

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.62496 23.14769 11.67193 -28.22219 0.3471588

――――――――――――――――――――――――――――――

［表２７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 23 : 4 cycles, criterion -1.359116

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.6649 23.16588 11.68137 -28.2449 0.3474096

――――――――――――――――――――――――――――――

［表２８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 24 : 4 cycles, criterion -1.358314

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.69192 23.17818 11.68776 -28.26025 0.3475789

――――――――――――――――――――――――――――――

［表２９］
――――――――――――――――――――――――――――――
***********************************************
Iteration 25 : 4 cycles, criterion -1.357772

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.71017 23.18649 11.69208 -28.27062 0.3476932

――――――――――――――――――――――――――――――

［表３０］
――――――――――――――――――――――――――――――
***********************************************
Iteration 26 : 4 cycles, criterion -1.357407

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.72251 23.19211 11.695 -28.27763 0.3477704

――――――――――――――――――――――――――――――

［表３１］
――――――――――――――――――――――――――――――
***********************************************
Iteration 27 : 4 cycles, criterion -1.35716

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.73084 23.19590 11.69697 -28.28237 0.3478225

――――――――――――――――――――――――――――――

［表３２］
――――――――――――――――――――――――――――――
***********************************************
Iteration 28 : 3 cycles, criterion -1.356993

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.73646 23.19846 11.6983 -28.28556 0.3478577

――――――――――――――――――――――――――――――

［表３３］
――――――――――――――――――――――――――――――
***********************************************
Iteration 29 : 3 cycles, criterion -1.356881

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.74026 23.20019 11.6992 -28.28772 0.3478814

――――――――――――――――――――――――――――――

［表３４］
――――――――――――――――――――――――――――――
***********************************************
Iteration 30 : 3 cycles, criterion -1.356805

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 6614 7191
regression coefficients
49.74283 23.20136 11.69981 -28.28918 0.3478975

1

misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 3 0 0
3 0 0 4 0
4 0 0 0 4

Identifiers of variables left in ordered categories model
6611 7188
――――――――――――――――――――――――――――――――――――

順序付きカテゴリーの例．
１５個のサンプル及び５０個の遺伝子に係るルオの前立腺データ。ｋ＝０及びｂ＝１ｅ７に対して、以下の結果を得る。

［表３５］
――――――――――――――――――――――――――――――
misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 2 1 0
3 0 0 4 0
4 0 0 0 4

Identifiers of variables left in ordered categories model
1
――――――――――――――――――――――――――――――

ｋ＝０．２５及びｂ＝１ｅ７に対して、以下の結果を得る。

［表３６］
――――――――――――――――――――――――――――――
misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 3 0 0
3 0 0 4 0
4 0 0 0 4

Identifiers of variables left in ordered categories model
1 42
――――――――――――――――――――――――――――――

ｋ＝０．２５及びｂ＝１ｅ７のときの反復の記録を以下に示す。

［表３７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 1 : 19 cycles, criterion -0.4708706

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 53
――――――――――――――――――――――――――――――

［表３８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 2 : 7 cycles, criterion -1.536822

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 53
――――――――――――――――――――――――――――――

［表３９］
――――――――――――――――――――――――――――――
***********************************************
Iteration 3 : 5 cycles, criterion -2.032919

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 42
――――――――――――――――――――――――――――――

［表４０］
――――――――――――――――――――――――――――――
***********************************************
Iteration 4 : 5 cycles, criterion -2.132546

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 18
――――――――――――――――――――――――――――――

［表４１］
――――――――――――――――――――――――――――――
***********************************************
Iteration 5 : 5 cycles, criterion -1.978462

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 Number of basis functions in model : 13
――――――――――――――――――――――――――――――

［表４２］
――――――――――――――――――――――――――――――
***********************************************
Iteration 6 : 5 cycles, criterion -1.668212

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 10 41 43 45
regression coefficients
34.13253 22.30781 13.04342 -16.23506 0.003213167 0.006582334 -0.0005943874 -3.557023

――――――――――――――――――――――――――――――

［表４３］
――――――――――――――――――――――――――――――
***********************************************
Iteration 7 : 5 cycles, criterion -1.407871

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 10 41 43 45
regression coefficients
36.90726 24.69518 14.61792 -17.16723 1.112172e-05 5.949931e-06 -3.892181e-08 -4.2906

――――――――――――――――――――――――――――――

［表４４］
――――――――――――――――――――――――――――――
***********************************************
Iteration 8 : 5 cycles, criterion -1.244166

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 10 45
regression coefficients
39.15038 26.51011 15.78594 -17.99800 1.125451e-10 -4.799167

――――――――――――――――――――――――――――――

［表４５］
――――――――――――――――――――――――――――――
***********************************************
Iteration 9 : 5 cycles, criterion -1.147754

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
40.72797 27.73318 16.56101 -18.61816 -5.115492

――――――――――――――――――――――――――――――

［表４６］
――――――――――――――――――――――――――――――
***********************************************
Iteration 10 : 5 cycles, criterion -1.09211

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
41.74539 28.49967 17.04204 -19.03293 -5.302421

――――――――――――――――――――――――――――――

［表４７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 11 : 5 cycles, criterion -1.060238

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
42.36866 28.96076 17.32967 -19.29261 -5.410496

――――――――――――――――――――――――――――――

［表４８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 12 : 5 cycles, criterion -1.042037

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
42.73908 29.23176 17.49811 -19.44894 -5.472426

――――――――――――――――――――――――――――――

［表４９］
――――――――――――――――――――――――――――――
***********************************************
Iteration 13 : 5 cycles, criterion -1.031656

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
42.95536 29.38894 17.59560 -19.54090 -5.507787

――――――――――――――――――――――――――――――

［表５０］
――――――――――――――――――――――――――――――
***********************************************
Iteration 14 : 4 cycles, criterion -1.025738

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.08034 29.47941 17.65163 -19.59428 -5.527948

――――――――――――――――――――――――――――――

［表５１］
――――――――――――――――――――――――――――――
***********************************************
Iteration 15 : 4 cycles, criterion -1.022366

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.15213 29.53125 17.68372 -19.62502 -5.539438

――――――――――――――――――――――――――――――

［表５２］
――――――――――――――――――――――――――――――
***********************************************
Iteration 16 : 4 cycles, criterion -1.020444

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.19322 29.56089 17.70206 -19.64265 -5.545984

――――――――――――――――――――――――――――――

［表５３］
――――――――――――――――――――――――――――――
***********************************************
Iteration 17 : 4 cycles, criterion -1.019349

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.21670 29.57780 17.71252 -19.65272 -5.549713

――――――――――――――――――――――――――――――

［表５４］
――――――――――――――――――――――――――――――
***********************************************
Iteration 18 : 3 cycles, criterion -1.018725

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.23008 29.58745 17.71848 -19.65847 -5.551837

――――――――――――――――――――――――――――――

［表５５］
――――――――――――――――――――――――――――――
***********************************************
Iteration 19 : 3 cycles, criterion -1.01837

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.23772 29.59295 17.72188 -19.66176 -5.553047

――――――――――――――――――――――――――――――

［表５６］
――――――――――――――――――――――――――――――
***********************************************
Iteration 20 : 3 cycles, criterion -1.018167

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.24208 29.59608 17.72382 -19.66363 -5.553737

――――――――――――――――――――――――――――――

［表５７］
――――――――――――――――――――――――――――――
***********************************************
Iteration 21 : 3 cycles, criterion -1.018052

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.24456 29.59787 17.72493 -19.66469 -5.55413

――――――――――――――――――――――――――――――

［表５８］
――――――――――――――――――――――――――――――
***********************************************
Iteration 22 : 3 cycles, criterion -1.017986

misclassification matrix
fhat
f 1 2
1 23 0
2 0 22
row =true class

Class 1 : Variables left in model
1 2 3 4 45
regression coefficients
43.24598 29.59889 17.72556 -19.6653 -5.554354

1

misclassification table
pred
y 1 2 3 4
1 4 0 0 0
2 0 3 0 0
3 0 0 4 0
4 0 0 0 4
Identifiers of variables left in ordered categories model
1 42
――――――――――――――――――――――――――――――

本発明の実施形態に係る方法のフローチャートである。 本発明の実施形態に係る別の方法のフローチャートである。 本発明の実施形態に係る装置のブロック図である。 本発明の実施形態に係るさらなる方法のフローチャートである。 本発明の実施形態に係る追加の方法を示すフローチャートである。 本発明の実施形態に係るさらに別の方法を示すフローチャートである。

Claims (30)

1. システムからの少なくとも１つのトレーニングサンプルを使用し、上記システムから取得されるデータに基づいて上記システムの複数の構成要素のうちのサブセットを識別する方法であって、
   上記システムの複数の構成要素と、上記複数の構成要素の一次結合に係る複数の重み係数との一次結合を取得するステップを含み、上記重み係数は上記少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の特徴を有し、
   上記既知の特徴の確率分布のモデルを取得するステップを含み、上記モデルは上記複数の構成要素の一次結合を条件とし、
   上記複数の構成要素の一次結合に係る重み係数の事前分布を取得するステップを含み、上記事前分布はゼロに近い高い確率密度を有する超事前分布を含み、上記超事前分布はジェフリーズの超事前分布ではないようなものであり、
   上記事前分布と上記モデルとを組み合わせて事後分布を生成するステップと、
   上記事後分布を最大化する複数の重み係数にてなるセットに基づいて上記複数の構成要素のうちのサブセットを識別するステップとを含む方法。
2. 上記一次結合を取得するステップは、ベイズの統計的方法を用いて上記複数の重み係数を推定するステップを含む請求項１記載の方法。
3. 上記複数の構成要素のうちの大部分は上記複数の構成要素のうちのサブセットの一部を形成する構成要素にはなりそうにない、という先験的仮定を立てるステップをさらに含む請求項１又は２記載の方法。
4. 上記超事前分布は、ゼロに近い事前分布が変更されることを可能にする１つ又は複数の調節可能なパラメータを含む、先行する請求項のうちのいずれか１つに記載の方法。
5. 上記モデルは、上記少なくとも１つのトレーニングサンプルから取得されるデータに基づいて確率分布を提供する尤度関数の形式である数式を含む、先行する請求項のうちのいずれか１つに記載の方法。
6. 上記尤度関数は、何らかの確率分布を記述するための先述のモデルに基づく請求項５記載の方法。
7. 上記モデルを取得するステップは、多項又は２項ロジスティック回帰、一般化線形モデル、コックスの比例ハザードモデル、加速度故障モデル、及びパラメトリック生存モデルを含むグループから上記モデルを選択するステップを含む、先行する請求項のうちのいずれか１つに記載の方法。
8. 上記多項又は２項ロジスティック回帰に基づくモデルは、
   

   

   の形式である請求項７記載の方法。
9. 上記一般化線形モデルに基づくモデルは、
   

   

   の形式である請求項７記載の方法。
10. 上記コックスの比例ハザードモデルに基づくモデルは、
    

    

    の形式である請求項７記載の方法。
11. 上記パラメトリック生存モデルに基づくモデルは、
    

    

    の形式である請求項７記載の方法。
12. 上記複数の構成要素のうちのサブセットを識別するステップは、上記事後分布の確率密度が最大化されるように反復手順を使用するステップを含む、先行する請求項のうちのいずれか１つに記載の方法。
13. 上記反復手順はＥＭアルゴリズムである請求項１２記載の方法。
14. 検査対象を複数の予め定義されたグループのうちの１つに分類することができる上記検査対象に係る複数の構成要素のうちのサブセットを識別するための方法であって、各グループは試験処理剤への反応によって定義され、上記方法は、
    複数の検査対象を上記試験処理剤にさらし、上記処理剤に対する反応に基づいて上記複数の検査対象を複数の反応グループにグループ化するステップと、
    上記複数の検査対象に係る複数の構成要素を測定するステップと、
    統計解析方法を用いて、上記複数の検査対象を複数の反応グループに分類することのできる複数の構成要素のうちのサブセットを識別するステップとを含む方法。
15. 上記統計解析方法は請求項１乃至１３のうちのいずれか１つに記載の方法を含む請求項１４記載の方法。
16. 検査対象に係る複数の構成要素のうちのサブセットを識別するための装置であって、上記サブセットは上記検査対象を複数の予め定義された反応グループのうちの１つに分類するために使用可能であり、各反応グループは、複数の検査対象を試験処理剤にさらしかつ上記処理剤に対する反応に基づいて上記複数の検査対象を複数の反応グループにグループ化することによって形成され、上記装置は、
    上記複数の検査対象に係る測定された複数の構成要素を受信するための入力と、
    統計解析方法を用いて、上記複数の検査対象を複数の反応グループに分類するために使用可能な複数の構成要素のうちのサブセットを識別する処理手段とを備える装置。
17. 上記統計解析方法は請求項１乃至１５のうちのいずれか１つに記載の方法を含む請求項１６記載の装置。
18. 試験化合物による処理に対して反応するもの又は反応しないものとして検査対象を分類することができる上記検査対象に係る複数の構成要素のうちのサブセットを識別するための方法であって、上記方法は、
    複数の検査対象を上記試験化合物にさらし、各検査対象の上記試験化合物に対する反応に基づいて上記複数の検査対象を複数の反応グループにグループ化するステップと、
    上記複数の検査対象に係る複数の構成要素を測定するステップと、
    統計解析方法を用いて、上記複数の検査対象を複数の反応グループに分類するために使用可能な複数の構成要素のうちのサブセットを識別するステップとを含む方法。
19. 上記統計解析方法は請求項１乃至１３のうちのいずれか１つに記載の方法を含む請求項１８記載の方法。
20. 検査対象に係る複数の構成要素のうちのサブセットを識別するための装置であって、上記サブセットは上記検査対象を複数の予め定義された反応グループのうちの１つに分類するために使用可能であり、各反応グループは、複数の検査対象を化合物にさらしかつ上記化合物に対する反応に基づいて上記複数の検査対象を複数の反応グループにグループ化することによって形成され、上記装置は、
    上記検査対象に係る測定された複数の構成要素を受信するように動作する入力と、
    統計解析方法を用いて、上記複数の検査対象を複数の反応グループに分類することができる複数の構成要素のうちのサブセットを識別するように動作する処理手段とを備える装置。
21. 上記統計解析方法は請求項１乃至１５のうちのいずれか１つに記載の方法を含む請求項２０記載の装置。
22. システムの複数のサンプルから生成されるデータから上記システムの構成要素のうちのサブセットを識別するための装置であって、上記サブセットは試験サンプルの特徴を予測するために使用可能であり、
    上記装置は処理手段を備え、上記処理手段は、
    上記システムの複数の構成要素の一次結合を取得し、上記複数の構成要素の一次結合に係る複数の重み係数を取得するように動作し、上記重み係数の各々は少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の特徴を有し、
    第２の特徴の確率分布のモデルを取得するように動作し、上記モデルは上記複数の構成要素の一次結合を条件とし、
    上記複数の構成要素の一次結合に係る複数の重み係数に関する事前分布を取得するように動作し、上記事前分布は、ゼロに近い事前確率質量が変更されることを可能にする調節可能な超事前分布を含み、上記超事前分布はジェフリーズの超事前分布ではなく、
    上記事前分布と上記モデルとを組み合わせて事後分布を生成するように動作し、
    上記事後分布を最大化する構成要素の重み係数を有する、複数の構成要素のうちのサブセットを識別するように動作する装置。
23. 上記処理手段は、ソフトウェアを実行するように構成されたコンピュータを備える請求項２２記載の装置。
24. 計算装置によって実行されたとき、上記計算装置に請求項１乃至１３のうちのいずれか１つに記載の方法を実行させるコンピュータプログラム。
25. 請求項２４記載のコンピュータプログラムを備えた、コンピュータが読み取り可能な媒体。
26. システムからのサンプルを検査して上記サンプルの特徴を識別する方法であって、上記特徴の症状を示す複数の構成要素のうちのサブセットについて検査するステップを含み、上記複数の構成要素のうちのサブセットは請求項１乃至１５のうちのいずれか１つに記載の方法を用いて決定されている方法。
27. 上記システムは生物学的システムである請求項２６記載の方法。
28. システムからのサンプルを検査して上記サンプルの特徴を決定するための装置であって、請求項１乃至１５のうちのいずれか１つに記載の方法に従って識別される複数の構成要素を検査するための手段を備える装置。
29. 計算装置によって実行されたときに、システムからの試験サンプルの特徴を予測するために使用可能な上記システムからの複数の構成要素を識別する方法を上記計算装置に実行させるコンピュータプログラムであって、
    複数のトレーニングサンプルから生成されるデータから複数の構成要素と構成要素に係る複数の重み係数との一次結合が生成され、各トレーニングサンプルは既知の特徴を有し、
    事後分布は、ゼロに近い確率質量が変更されることを可能にする調節可能な超事前分布を含む構成要素に係る複数の重み係数の事前分布であって、上記超事前分布はジェフリーズの超事前分布ではない事前分布と、上記一次結合を条件とするモデルとを組み合わせて、上記事後分布を最大化する構成要素に係る複数の重み係数を推定することによって生成されるコンピュータプログラム。
30. 生物学的システムの複数の構成要素のうちのサブセットを識別する方法であって、上記サブセットは上記生物学的システムからの試験サンプルの特徴を予測することができ、上記方法は、
    上記システムの複数の構成要素と、上記複数の構成要素の一次結合に係る複数の重み係数との一次結合を取得するステップを含み、上記重み係数の各々は少なくとも１つのトレーニングサンプルから取得されるデータに基づく値を有し、上記少なくとも１つのトレーニングサンプルは既知の特徴を有し、
    上記既知の特徴の確率分布のモデルを取得するステップを含み、上記モデルは上記複数の構成要素の一次結合を条件とし、
    上記複数の構成要素の一次結合に係る複数の重み係数に関する事前分布を取得するステップを含み、上記事前分布は、ゼロに近い確率質量が変更されることを可能にする調節可能な超事前分布を含み、
    上記事前分布と上記モデルとを組み合わせて事後分布を生成するステップと、
    上記事後分布を最大化する複数の重み係数に基づいて複数の構成要素のうちのサブセットを識別するステップとを含む方法。

Images