JP2005524124A - システムの診断構成要素を識別するための方法および装置 - Google Patents

Abstract

テストサンプルの特徴を予測することのできるシステムの構成要素の部分集合を識別するための方法および装置が記載されている。この方法は、構成要素と構成要素加重との線形結合を発生するステップを含んでおり、この線形結合において、各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が決定される。ある特徴の確率分布に対するモデルが規定され、このモデルは線形結合を条件とし、また、このモデルは２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではない。事前（ｐｒｉｏｒ）分布はゼロに近い高い確率密度を有するハイパープライアを含む線形結合の構成要素加重に対して構成され、また、事前分布とモデルとを組合せることにより、事後分布が生成される。この事後分布を最大にする構成要素加重を有する構成要素の部分集合が識別される。

Description

本発明は、システムのサンプルから発生されたデータから、そのシステム内のサンプルの特徴を予測することのできるシステムの構成要素を識別するための方法および装置に関し、とくに、本発明は、それに限定されるものではないが、生物学的方法によって発生されたデータから、その生物学的システムからのサンプルと関連付けられた重要な特徴を予測することのできる生物学的システムの構成要素を識別するための方法および装置に関する。

種々の重要な特徴に分類されることのできる多くの“システム”が存在する。この“システム”という用語には、本質的に、データが提供されることのできる全てのタイプのシステムが含まれ、それには化学的システム、金融システム（たとえば、個人、グループまたは法人向けクレジットシステム、ローン履歴）、地質学的システムおよびその他多数が含まれる。システムのサンプルから特定の特徴を識別するために［たとえば、金融システムの解析により、金融システム内に存在するグループを識別する（信用が“申し分のない”者および信用貸しは危険である者という非常に簡単な関係で）ことを補助するために］システムから発生されたデータ（たとえば、統計的データ）を使用できることが望ましい。大量の統計的データが存在している場合、そのシステムからのサンプルの特定の特徴を予測するそのデータから構成要素を識別することは困難な作業である。これは一般に、処理すべき大量のデータの大部分が、そのデータが得られた特定のサンプルの重要な特徴を全く、あるいはほとんど示さない可能性があるためである。さらに、テストサンプルデータが訓練サンプルデータに関して高度の可変性を有しているときには、その訓練サンプルデータを使用して識別された構成要素はそのテストサンプルデータに関する特徴の識別には役に立たないことが多い。これは、たとえば、データが多数の異なったソースから得られる状況においてよくあるケースであり、それは、個々のソースのそれぞれからデータが収集される条件を制御することが不可能なことが多いためである。

これらの問題がとくに関係するシステムのタイプの一例は生物学的システムであり、以下の説明はとくに生物学的システムに該当するものである。しかしながら、本発明は生物学的システムによる使用に限定されず、任意のシステムに対して一般的に適用される。

バイオテクノロジーにおける近年の進歩の結果、システムの大規模スクリーニングおよびサンプルの解析に対する生物学的方法が発達してきた。このような方法には、たとえば、ＤＮＡ、ＲＮＡまたは抗体マイクロアレイ解析、プロテオミクス解析、プロテオミクス電気泳動ゲル解析および高度なスループットのスクリーニング技術が含まれる。これらのタイプの方法の結果、テストされたサンプルのそれぞれに対して３０，０００以上の構成要素を有することのできるデータがしばしば発生される。

生物学的システムからのサンプル中の重要な特徴を識別できることが明らかに重要である。たとえば、“病気にかかっている”および“病気にかかっていない”のようなグループを分類すること等が挙げられる。これらの生物学的方法の多くは、生物学的システム中のサンプルの特徴を予測する診断ツール（たとえば、組織または体液をスクリーニングすることにより病気を識別する）として（、あるいは、たとえば、薬剤化合物の効能を決定するツールとして）有用である。

今日、このような適用におけるバイオテクノロジーアレイのような生物学的方法の使用は、これらのタイプの方法から発生された大量のデータと、および意味のある結果を得るためにデータをスクリーニングする効率的な方法の欠如のために制限されている。その結果、従来技術の方法を使用した生物学的データの解析は、データ中の情報を十分に利用することができないか、あるいは時間を要し、誤った陽性および陰性の結果を出し易く、また、データから意味のある結果が得られなければならない場合には大量のコンピュータメモリを必要とする。これは、迅速で正確なスクリーニングが要求される規模スクリーニングシナリオにおいて問題である。

したがって、とくに生物学的データを解析する改善された方法が必要とされており、さらに一般的には、システムからのサンプルに対する重要な特徴を予測するために任意のシステムからのデータを解析する改善された方法が必要とされている。

第１の特徴において、本発明は、テストサンプルの特徴を予測することのできるシステムの構成要素の部分集合を識別する方法を提供し、この方法は、
（ａ）構成要素と構成要素加重との線形結合を発生し、この線形結合において、各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が決定され、
（ｂ）ある特徴の確率分布に対するモデルを規定し、このモデルは線形結合を条件付けし、また、このモデルは２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではなく、
（ｃ）ゼロに近い高い確率密度を有するハイパープライア（ｈｙｐｅｒｐｒｉｏｒ）を含む線形結合の構成要素加重に対して事前（ｐｒｉｏｒ）分布を構成し、
（ｄ）この事前分布とモデルとを組合せて事後（ｐｏｓｔｅｒｉｏｒ）分布を生成し、
（ｅ）この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別するステップを含んでいる。

この方法は、訓練サンプルに対する特徴を予測することのできる構成要素の部分集合を識別するために既知の特徴を有する訓練サンプルを使用する。続いて、たとえば、臨床テストのようなテストに対して構成要素の部分集合情報を使用して、組織サンプルが悪性であるかまたは良性であるか、あるいは腫瘍の重量はどれくらいか等の特徴を予測することができ、あるいは特定の症状を有する患者の推定生存時間を提供することができる。ここで使用されているように、“特徴”という用語は、あるサンプルと関連付けられた任意の応答または識別可能な特性または特質のことである。たとえば、ある特徴は、特定のサンプルに対するある事象までの特定の時間であってもよいし、あるいはあるサンプルの寸法または量、もしくはあるサンプルが分類されることのできるクラスまたはグループであってもよい。

本発明の方法は、ベイズの統計的方法を使用して構成要素の加重を推定する。大量の構成要素がシステムから発生された場合（通常、本発明の方法が有効であるケースである）、本発明は、大部分の構成要素は、ある特徴を予測するための構成要素の部分集合の一部を形成している構成要素ではなさそうであるというアプリオリ仮定を行うことが好ましい。したがって、大部分の構成要素の加重はゼロであるようだという仮定がなされる。この仮定を念頭に置いて加重の事後（ｐｏｓｔｅｒｉｏｒ）確率が最大化されるように構成要素加重を設定したあるモデルが構成される。予め定められたしきい値より下の加重を有する構成要素（アプリオリ仮定によるとそれらの大部分となるであろう）はその適用を免除される。このプロセスは、残った診断構成要素が識別されるまで繰返される。この方法は、主に、大部分の構成要素を迅速に除去するアプリオリ仮定のために迅速である。

あるシステムのほとんどの特徴は典型的に確率分布を示し、ある特徴の確率分布は、訓練サンプルから発生されたデータに基づいた統計モデルを使用してモデル化されることができる。本発明の方法は、１つのまたは一連の重要な特徴に対する確率分布をモデル化する統計モデルを使用する。したがって、特定の確率分布を有する１つの重要な特徴に対して、その分布をモデル化した１つの適切なモデルが規定される。この方法は、線形結合を条件とする任意のモデルを使用することができ、また、訓練サンプルから得られたデータに基づいて確率分布を提供する尤度関数の形態の数学的方程式であることが好ましい。この尤度関数は、ある確率分布を記述するために上述されたモデルに基づいていることが好ましい。１実施形態において、モデルは、多項または二項ロジスティック回帰、一般化された線形モデル、コックスの比例ハザードモデル、加速された故障（ｆａｉｌｕｒｅ）モデル、パラメトリック生存モデル、ｘ²分布モデル、または指数分布モデルからなるグループから選択されたモデルに基づいた尤度関数である。

１実施形態において、尤度関数は多項または二項ロジスティック回帰に基づいている。多項または二項ロジスティック回帰は、多項または二項分布を有する特徴をモデル化することが好ましい。二項分布は、オン／オフ状態のような２つの可能性のあるクラスまたはグループを有する統計的分布である。このようなグループの例には、死亡／生存、改善された／改善されていない、衰弱した／衰弱していない、が含まれる。多項分布は、複数のクラスまたはグループが複数のサンプルのそれぞれに対して可能である、換言すると、あるサンプルが複数のクラスまたはグループの１つに分類されることのできる二項分布を一般化したものである。したがって、多項または二項ロジスティック回帰に基づいて尤度関数を規定することにより、サンプルを複数のグループまたはクラスの１つに分類することのできる構成要素の部分集合を識別することが可能である。これを行うために、訓練サンプルはそれの予め定められた特徴に基づいて複数のサンプルグループ（または“クラス”）にグループ化され、この予め定められた特徴において、各サンプルのメンバーは共通の特徴を有しており、共通のグループ識別子を割当てられる。尤度関数は、線形結合（グループ化された訓練サンプルから発生されたデータを含む）を条件とする多項または二項ロジスティック回帰に基づいて公式化される。特徴は、訓練サンプルがグループ化される任意の所望の分類であってもよい。たとえば、組織サンプルを分類する特徴は、その組織が正常、悪性または良性であることでもよいし、あるいはその細胞が白血病細胞または健康な細胞であることでもよいし、もしくは訓練サンプルがある症状を有する患者または有しない患者のいずれの血液から得られたことでもよく、あるいは訓練サンプルが正常な細胞と比較されるいくつかのタイプの癌の１つからの細胞から得られたことであってもよい。

尤度関数は、

の形態のロジスティック回帰に基づいていることが好ましく、ここで、
ｘ_i ^Tβ_gは、構成要素加重β_gを有する訓練サンプルｉからの入力データから発生された線形結合であり、
ｘ_i ^Tは、Ｘのｉ番目の行に対する構成要素であり、β_gはサンプルクラスｇに対する構成要素加重のセットであり、
訓練サンプルｉがクラスｇのメンバーである場合、ｅ_ig＝１であり、そうでない場合にはｅ_ig＝０であり、
Ｘは、ｐ個の構成要素を含むｎ個の訓練サンプルからのデータである。

別の実施形態において、尤度関数は順序付けられたカテゴリー的ロジスティック回帰に基づいている。この順序付けられたカテゴリー的ロジスティック回帰は、クラスが特定の順序にされている（たとえば、病気の重大度が増加または減少していくクラスのような順序付けられたクラスの）多項分布をモデル化する。順序付けられたカテゴリー的ロジスティック回帰に基づいて尤度関数を規定することにより、あるサンプルをあるクラスに分類することのできる構成要素の部分集合を識別することが可能であり、この場合そのクラスは複数の予め規定された順序付けられたクラスの１つである。各グループ識別子がある順序付けられたクラスのあるメンバーに対応している一連のグループ識別子を規定し、訓練サンプルをその訓練サンプルの予め定められた特徴に基づいてその順序付けられたクラスの１つにグループ化することにより、尤度関数は、線形結合（グループ化された訓練サンプルから発生されたデータを含む）を条件とする順序付けられたカテゴリー的ロジスティック回帰に基づいて公式化されることができる。

尤度関数は、

の形態の順序付けられたカテゴリー的ロジスティック回帰に基づいていることが好ましく、ここで、
γ_ikは、訓練サンプルｉがｋ以下である識別子を有するクラスに属している確率であり（ここで、順序付けられたクラスの合計はＧである）、
ｘ_i ^Tβ^*は、構成要素加重β^* _iを有する訓練サンプルｉからの入力データから発生された線形結合であり、
ｘ_i ^Tは、Ｘのｉ番目の行に対する構成要素であり、
ｒ_ijは、

として定義されており、ここで、
観察ｉがクラスｊ中のものである場合はｃ_ij＝１であり、そうでない場合はｃ_ij＝０である。

本発明の別の実施形態において、尤度関数は一般化された線形モデルに基づいている。一般化された線形モデルは、分布の正規指数族に属する分布を有する特徴をモデル化することが好ましい。正規指数族分布の例には、正規分布、ガウス分布、ポアソン分布、ガンマ分布および逆ガンマ分布が含まれる。したがって、本発明の方法の別の実施形態では、予測されるべき特性をモデル化する一般化された線形モデルを規定することにより、分布の正規指数族の範囲内であるサンプルの予め規定された特性を予測することのできる構成要素の部分集合が識別される。一般化された線形モデルを使用して予測されることのできる特性の例には、たとえば、あるサンプルの重量、寸法、カウント、グループメンバーシップその他のディメンションまたは量あるいは特性のような指定された分布を示すあるサンプルの任意の量が含まれる。

一般化された線形モデルは、

の形態のものであることが好ましく、ここで、
ｙ＝（ｙ₁，・・・，ｙ_n）^Tであり、ｙ_iはｉ番目のサンプルに関して測定された特性であり、
ａ_i（φ）＝φ／ｗ_iであり、ｗ_iは既知の加重の固定されたセットであり、φは単一の
スケールパラメータであり、
関数ｂ（．）およびｃ（．）はNelderおよびWedderburn氏等による文献（1972）において定義されたものであり、
Ｅ｛ｙ_i｝＝ｂ’（θ_i）
Ｖａｒ｛ｙ｝＝ｂ’’（θ_i）ａ_i（φ）＝τ_i ²ａ_i（φ）であることが好ましい。

各観察はコバリエイトｘ_iと線形予測子η_i＝ｘ_i ^Tβとのセットを有していることが好ましい。ｉ番目の観察の平均とそれの線形予測子との間の関係は、リンク関数η_i＝ｇ（μ_i）＝ｇ（ｂ’（θ_i））によって与えられることが好ましい。このリンクの逆数はｈによって表され、それは：
Ｅ｛ｙ_i｝＝ｂ’（θ_i）＝ｈ（η_i）
であることが好ましい。

別の実施形態において、本発明の方法は、ある事象がデータを得た時点では発生していないと仮定してその事象までの時間の確率を推定することが好ましいハザードモデルに基づいた尤度関数を使用することにより、あるサンプルに対して事象までの時間を予測するために使用されてもよい。１実施形態において、尤度関数は、コックスの比例ハザードモデル、パラメトリック生存モデルおよび加速された故障時間モデルからなるグループから選択されたモデルに基づいている。コックスの比例ハザードモデルは、ハザード関数の形式に関して制限的仮定を行わずに、ある事象までの時間が構成要素と構成要素加重のセットに関してモデル化されることを可能にする。加速された故障モデルは、構成要素の測定値が時間スケールに乗法的に作用すると仮定されているので、個体が時間軸に沿って進行する速度に影響を与える生存時間からなるデータに対する一般的なモデルであることが好ましい。したがって、加速された生存モデルは、たとえば、病気の進行速度に関して解釈されることができる。パラメトリック生存モデルは、ある事象までの時間（たとえば、生存時間）に対する分布関数が既知の分布によってモデル化されているか、あるいは特定されたパラメトリック公式を有するものである。一般に使用されている生存分布には、ワイブル指数および極値分布がある。

あるサンプルに対するある事象までの時間を予測することのできる構成要素の部分集合は、コックスの比例ハザードモデル、パラメトリック生存モデルまたは加速された生存時間モデルに基づいて尤度を規定することにより識別され、この尤度は、複数のサンプルに対してそのサンプルが得られた時間から事象の時間までの経過した時間を測定することを含んでいる。

ある事象までの時間を予測する尤度関数は、以下の形態のものであることが好ましい：

コックスの比例ハザードモデルに基づいた尤度関数は、

の形態のものであることが好ましく、ここで、
ＺはＸの行の再配列であるマトリックスであることが好ましく、ここでＺの行の順序付けは生存時間の順序付けにより誘発された順序付けに対応し、ｄは生存時間を順序付けるために必要とされる同じ置換により検閲削除インデックスを順序付けた結果である。また、Ｚ_jはマトリックスＺのｊ番目の行であり、ｄ_jはｄのｊ番目のエレメントであり、ここで、

であり、Ｒ_j＝｛ｉ：ｉ＝ｊ，ｊ＋１，・・・，Ｎ｝＝ｊ番目の順序付けられた事象時間ｔ(_j)におけるリスクセットである。

パラメトリック生存モデルに基づいた対数尤度関数は、以下の形態のものであることが好ましい：

ｉ番目のサンプルが検閲削除されていない場合はｃ_i＝１、ｉ番目のサンプルが検閲削除されていない場合はｃ_i＝０である。尤度関数のこの形態は、ワイブル指数および極値分布によって共用される。関数λ（．）およびΛ（．）は、AitkinおよびClayton氏等による文献（1980）において定義されたものである。

任意の規定されたモデルに対して、典型的にベイズの統計的モデルを使用して構成要素加重が推定され（KotzおよびJohnson氏等による文献、1983）、その場合、尤度関数と事前分布を組合せる構成要素加重の事後分布が公式化される。構成要素加重は、各訓練サンプルに対して発生されたデータを与えられた加重の事後分布を最大化することによって推定される。したがって、最大化されるべき目標関数は、上述した特徴に対するモデルと加重に対する事前分布とに基づいた尤度関数からなる。

事前分布は、

の形態のものであることが好ましく、
ここで、ｖはハイパーパラメータのｐ×１ベクトルであり、また、ｐ（β|ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、ｐ（ｖ²）はｖ²に対するあるハイパープライア分布である。このハイパープライア分布（この方法の全ての実施形態に対して同じであることが好ましい）は、種々の表記規定を使用して表現されてもよく、好ましい実施形態の詳細な説明（以下、参照）においては、以下の表記規定は単に特定の好ましい実施形態に対する便宜性のために採用されたに過ぎない。

ここで使用されているように、確率分布に対する尤度関数が多項または二項ロジスティック回帰に基づいているとき、事前分布に対する表記は：

また、ｐ（β_g|τ_g ²）はＮ（０，ｄｉａｇ｛τ_g ²｝）であり、ｐ（τ_g ²）はτ_g ²に対するあるハイパープライア分布である。

ここで使用されているように、確率分布に対する尤度関数が順序付けられたカテゴリー的ロジスティック回帰に基づいているとき、事前分布に対する表記は：

ここで、β₁，β₂，・・・，β_nは構成要素加重であり、ｐ（β_i|τ_i）はＮ（０，τ_i ²）であり、ｐ（τ_i）はτ_iに対するあるハイパープライア分布である。

ここで使用されているように、分布に対する尤度関数が一般化された線形モデルに基づいているとき、事前分布に対する表記は：

ここで、ｖはハイパーパラメータのｐ×１ベクトルであり、また、ｐ（β|ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、ｐ（ｖ²）はｖ²に対するある事前分布である。

ここで使用されているように、分布に対する尤度関数がハザードモデルに基づいているとき、事前分布に対する表記は：

ここで、ｐ（β^*|ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、ｐ（ｖ²）はｖ²に対するあるハイパープライア分布である。

事前分布は、可能な限り常に０加重が好ましいことを保証するハイパープライアを含んでいる。
ハイパープライアは、ジェフェリーのハイパープライアである（KotzおよびJohnson氏等の文献、1983）。

上述したように、事後分布を生成するために事前分布と尤度関数が組合せられる。事後分布は以下の形態のものであることが好ましい：

事後分布における構成要素加重は、事後分布の確率密度が最大化されるように反復的な処理手続きで推定されることが好ましい。反復的な処理手続き中、予め定められたしきい値より小さい値を有する構成要素加重は、好ましくはこれらの構成要素をゼロに設定することによって推定される。この結果、対応した構成要素が除去される。

反復的な処理手続きは、ＥＭアルゴリズムであることが好ましい。ＥＭアルゴリズムは、事後分布の確率密度を最大化する所定の構成要素加重に収斂する構成要素加重推定値のシーケンスを処理する。ＥＭアルゴリズムは、Ｅまたは期待ステップおよびＭまたは最大化ステップとして知られている２つのステップから構成されている。Ｅステップにおいて、観察されたデータおよび現在のパラメータ値を条件とする対数−事後関数期待値が決定される。Ｍステップにおいて、尤度を増加させる更新された構成要素加重推定値を求めるために、期待された対数−事後関数が最大化される。２つのステップは、ＥステップおよびＭステップの収斂が達成されるまで、換言すると、期待値および対数−事後関数の最大化された値が収斂するまで交互に実施される。

本発明の方法は、測定値が得られることのできる任意のシステムに対して、好ましくは、非常に大量のデータが発生されるシステムに対して適用可能であると考えられる。本発明の方法が適用可能なシステムの例には、たとえば、クレジットリスク査定システム、保険システム、マーケティングシステムまたは企業記録システム、電子システム、物理的システム、天体物理学的システムおよび機械的システムが含まれる。たとえば、金融システムにおいては、サンプルは特定の株であることができ、構成要素は、企業収益、従業員数、株主の数等の、株価に影響を与える可能性のある任意の数のファクタに対してなされた測定値であることができる。

本発明の方法は、とくに、生物学的システムの解析における使用に適している。本発明の方法は、構成要素に対する測定可能な値を生成し、その構成要素が特有のラベルを付けられる任意の生物学的システムからのサンプルを分類するために、構成要素の部分集合を識別するために使用されることができる。換言すると、構成要素は、１つの構成要素からのデータと別の構成要素からのデータとの区別を可能にするようにラベル付けされるか、あるいは編成される。たとえば、構成要素は、たとえば、各構成要素からのデータと別のものとを空間的な位置で区別することをを可能にするアレイで空間的に編成されてもよいし、あるいは各構成要素は、それと関連付けられた、識別信号またはタグのようなある特有の識別物を有していてもよい。たとえば、構成要素は個々のキャリアに結び付けられてもよく、各キャリアは量子ドット［たとえば、Rosenthal氏による文献（2001，Nature Biotech 19: 621-622）；Han氏他による文献（2001，Nature Biotechnology 19: 631-635）を参照］、蛍光マーカ［たとえば、Fu氏他による文献（1999，Nature Biotech 17: 1109-1111）を参照］、バーコード付きタグ［たとえば、LockhartおよびTrulson氏等による文献（2001，Nature Biotechnology 19: 1122-1123）を参照］のような検出可能な識別シグナチャを有している。

とくに好ましい実施形態において、生物学的システムはバイオテクノロジーアレイである。バイオテクノロジーアレイ［それらの例はSchena氏他による文献（1995，Science 270: 467-470）；Lockhart氏他による文献（1996，Nature Biotechnology 14: 1649）；米国特許第 5,569,5880号明細書に記載されている］の例には、オリゴヌクレオチドアレイ、ＤＮＡアレイ、ＤＮＡマイクロアレイ、ＲＮＡアレイ、ＲＮＡマイクロアレイ、ＤＮＡマイクロチップ、ＲＮＡマイクロチップ、蛋白質アレイ、蛋白質マイクロチップ、抗体アレイ、化学物質アレイ、炭水化物アレイ、プロテオミクスアレイ、脂質アレイが含まれる。別の実施形態において、生物学的システムは、たとえば、ＤＮＡまたはＲＮＡ電気泳動ゲル、蛋白質またはプロテオミクス電気泳動ゲル、Ｂｉａｃｏｒｅ解析、アミノ酸解析、ＡＤＭＥＴｏｘスクリーニング［たとえば、High-throughput ADMETox estimation： In Vitro and In Silico approaches(2002)，Ference Darvas and Gyorgy Dorman(Eds)，Biotechniques Pressを参照］のような生体分子相互作用解析、蛋白質電気泳動ゲルおよびプロテオミクス電気泳動ゲルを含むグループから選択されてもよい。

構成要素は、システムの任意の測定可能な構成要素であることができる。生物学的システムの場合、構成要素は、たとえば、遺伝子またはその一部、ＤＮＡシーケンス、ＲＮＡシーケンス、ペプチド、蛋白質、炭水化物分子、脂質またはその混合物、生理学的構成要素、解剖学的構成要素、疫学的構成要素または化学的構成要素であってもよい。

訓練サンプルは、サンプルの特徴が知られているシステムから得られた任意のデータであることができる。たとえば、訓練サンプルは発生されたデータであってもよい。たとえば、訓練サンプルは、生物学的システムに適用されたサンプルから得られた任意のデータであることができる。たとえば、生物学的システムがＤＮＡマイクロアレイであるとき、訓練サンプルは、既知の特徴を有する細胞から抽出されたＲＮＡとのアレイのハイブリダイゼーションに後続してそのアレイから得られたデータ、または細胞から抽出されたＲＮＡから合成されたｃＤＮＡであることができ、あるいは生物学的システムがプロテオミクス電気泳動ゲルである場合には、訓練サンプルは、そのシステムに適用された蛋白質または細胞抽出から発生されてもよい。

本発明者は、本発明の方法がテスト治療に反応して混合した結果を示した被験者からのテストデータを再推定または推定するときに１実施形態において使用されてもよいと考えている。したがって、第２の特徴において、本発明は、被験者を複数の予め規定されたグループの１つに分類することのできる被験者の構成要素の部分集合を識別する方法を提供し、各グループは、
（ａ）複数の被験者をテスト治療に露出し、その治療に対する反応に基づいて被験者を反応グループにグループ化し、
（ｂ）被験者の構成要素を測定し、
（ｃ）統計的解析方法を使用して被験者を反応グループに分類することのできる構成要素の部分集合を識別するステップを含むテスト治療に対する反応によって規定される。
統計的解析方法は、本発明の第１の特徴による方法であることが好ましい。

構成要素の部分集合が識別されると、その部分集合を使用して、テスト治療に反応する可能性の高いものおよび反応しないもののようなグループに被験者を分類することができる。このようにして、本発明の方法は、人口のほんの一部分に対して有効である可能性がある治療が識別されることを可能にし、また、テスト治療に対して反応するであろう人口のその一部分の識別を可能にする。

第３の特徴において、本発明は、各反応グループが複数の被験者をテスト治療に露出し、被験者をその治療に対する反応に基づいて反応グループにグループ化することにより形成された複数の予め規定された反応グループの１つに被験者を分類することのできる被験者の構成要素の部分集合を識別する装置を提供し、この装置は、
（ａ）被験者の測定された構成要素を受取る手段と、
（ｂ）統計的解析方法を使用して被験者を反応グループに分類することのできる構成要素の部分集合を識別する手段とを備えている。
統計的解析方法は、本発明の第１または第２の特徴による方法であることが好ましい。
第４の特徴において、本発明は、被験者をテスト化合物による治療に反応した者あるいは反応しなかった者として分類することのできる被験者の構成要素の部分集合を識別する方法を提供し、この方法は、
（ａ）複数の被験者を化合物に露出し、その化合物に対する各被験者の反応に基づいて被験者を反応グループにグループ化し、
（ｂ）被験者の構成要素を測定し、
（ｃ）統計的解析方法を使用して被験者を反応グループに分類することのできる構成要素の部分集合を識別するステップを含んでいる。
統計的解析方法は、第１の特徴による方法であることが好ましい。

第５の特徴において、本発明は、各反応グループが複数の被験者を化合物に露出し、被験者をその化合物に対する反応に基づいて反応グループにグループ化することにより形成された複数の予め規定された反応グループの１つに被験者を分類することのできる被験者の構成要素の部分集合を識別する装置を提供し、この装置は、
（ｃ）被験者の測定された構成要素を受取る手段と、
（ｄ）統計的解析方法を使用して被験者を反応グループに分類することのできる構成要素の部分集合を識別する手段とを備えている。
統計的解析方法は、本発明の第１または第２の特徴による方法であることが好ましい。
本発明の第２乃至第５の特徴において測定された構成要素は、たとえば、その被験者の遺伝子または小さいヌクレオチド多型（ＳＮＰｓ）、蛋白質、抗体、炭水化物、脂質、あるいは任意の別の測定可能な構成要素であることができる。

とくに好ましい実施形態において、化合物は、薬剤化合物、または薬剤化合物と薬剤的に許容可能なキャリアとを含む合成物である。

本発明の識別方法は、適切なコンピュータソフトウェアおよびハードウェアによって実施されてもよい。

第６の特徴において、本発明は、テストサンプルの特徴を予測することのできるシステムから発生されたデータからのシステムの構成要素の部分集合を、そのシステムからの複数のサンプルから識別する装置を提供し、この装置は、
（ａ）各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が導かれる、構成要素と構成要素加重との線形結合を発生する手段と、
（ｂ）線形結合を条件とし、また、２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではない、ある特徴の確率分布に対するモデルを規定する手段と、
（ｃ）ゼロに近い高い確率密度を有するハイパープライアを含む線形結合の構成要素加重に対して事前分布を構成する手段と、
（ｄ）この事前分布とモデルとを組合せて事後分布を生成する手段と、
（ｅ）この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別する手段とを具備している。
この装置は、適切にプログラムされた計算装置を備えていてもよい。

第７の特徴によると、本発明は、計算装置にロードされたときにこの計算装置を制御して、本発明の第１の特徴による方法を実施するように構成されているコンピュータプログラムを提供する。

コンピュータプログラムは、上述した本発明の第１または第２の特徴の好ましいアルゴリズムおよび方法ステップの任意のものを実施することができる。

本発明の第８の特徴によると、本発明の第４の特徴によるコンピュータプログラムを提供するコンピュータ読出し可能な媒体が提供される。

本発明の第９の特徴によると、システムからのサンプルをテストしてそのサンプルの特徴を識別する方法が提供され、この方法は、特徴の診断に役立つ構成要素の部分集合をテストするステップを含んでおり、その構成要素の部分集合は本発明の第１または第２の特徴による方法によって決定されている。
システムは生物学的システムであることが好ましい。

本発明の第１０の特徴によると、サンプルの特徴を決定するためにシステムからのサンプルをテストする装置が提供され、この装置は本発明の第１または第２の特徴による方法にしたがって識別された構成要素をテストする手段を備えている。

本発明の第１１の特徴によると、本発明は、計算装置上で実行されたときに、システムからのテストサンプルの特徴を予測することのできるそのシステムからの構成要素を識別する方法でこの計算装置を制御するように構成されており、構成要素と構成要素加重との線形結合は、各訓練サンプルが１つの既知の特徴を有する複数の訓練サンプルから発生されたデータから発生され、ゼロに近い高い確率密度を有するハイパープライアを含む構成要素加重に対する事前分布と、線形結合を条件とするモデルとを組合せることによって事後分布が発生され、それによって事後分布を最大にする構成要素加重を推定し、前記モデルは２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではないコンピュータプログラムを提供する。

本発明の特徴が計算装置によって実施される場合、たとえば、ＰＣまたは主フレームあるいはネットワーク計算インフラストラクチャ等の、任意の適切なコンピュータハードウェアが使用されてもよいことが認識されるであろう。

第１２の特徴において、本発明は、生物学的システムからのテストサンプルの特徴を予測することのできる生物学的システムの構成要素の部分集合を識別する方法を提供し、この方法は、
（ａ）構成要素と構成要素加重との線形結合を発生し、この線形結合において、各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が決定され、
（ｂ）ある特徴の確率分布に対するモデルを規定し、このモデルは線形結合を条件とし、
（ｃ）ゼロに近い高い確率密度を有するハイパープライアを含む線形結合の構成要素加重に対して事前分布を構成し、
（ｄ）この事前分布とモデルとを組合せて事後分布を生成し、
この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別するステップを含んでいる。

図１は、本発明の１実施形態の予測の成功に関する置換テストの結果を示している。クラスラベルは２００回ランダムに入れ替えられ、各置換に対して解析が繰返された。このヒストグラムは置換の下での予測の成功分布を示している。ｘ軸は正しく分類されたサンプルの数を示し、ｙ軸は頻度を示している。

図２は、本発明の１実施形態の予測の成功に関する置換テストの結果を示している。クラスラベルは２００回ランダムに入れ替えられ、各置換に対して解析が繰返された。このヒストグラムはクラスラベルの置換の下での予測の成功分布を示している。ｘ軸は合計サンプルの割合を示し、ｙ軸（ラムダ）は正しく分類されるケースの割合を示している。

図３は、本発明の方法の１実施形態において使用される一般化された線形モデルに対する曲線のプロットを示している。この方法により選択された５つの構成要素がこのモデル中で使用されたとき、あてはめられた曲線（実線）が生成され、真の曲線（点線）は破線として示されており、その５つの構成要素に基づく２００回の観察（ｘ軸）からのデータ（ｎｆ，ｙ軸）が円として示されている。

図４は、本発明の方法の１実施形態を使用して識別された単一の遺伝子に対するあてはめられた確率のプロットを示している。ｘ軸は遺伝子インデックスを示し、ｙ軸は特定の順序にされたクラスに属するサンプルの確率を示している。ラインは以下のようにクラスを表している：破線＝クラス１、実線＝クラス２、点線＝クラス３、点線と破線を組合せた線＝クラス４。

図５は、本発明によるシステムを構成するために使用されるパーソナルコンピュータの概略図である。

本発明は、特定の訓練サンプルが特定の特徴を有しているか否かを識別するために使用されることのできる、好ましくは最少数の、構成要素を識別する。この最少数の構成要素はその特徴の“診断に役立ち”、あるいは異なった特徴を有するサンプルを区別することを可能にする。本質的に、システムから発生されたデータの全てから、本発明の方法は、特定の特徴をテストするために使用されることのできる最少数の構成要素を識別することを可能にする。この方法によってこれらの構成要素が識別されると、その構成要素は新しいサンプルを評価するために将来使用されることができる。本発明の方法は統計的方法を使用して、特徴を正しく予測するのに必要とされない構成要素を除去する。

本発明者は、訓練サンプルから発生されたデータの構成要素の線形結合の構成要素加重が、訓練サンプルの特徴を正しく予測するのに必要とされない構成要素を消去する方法で推定可能であることを発見した。その結果、訓練サンプルの特徴を正しく予測することのできる構成要素の部分集合が識別される。このようにして、本発明の方法は、ある特徴を正しく予測することのできる比較的少数の構成要素を大量のデータから識別することを可能にする。

また、本発明の方法には、使用する必要のあるコンピュータメモリが、構成要素に関する限界情報ではなく共有（ｊｏｉｎｔ）情報を使用する従来技術の方法よりも少ないという利点がある。したがって、本発明の方法は、たとえば、ラップトップマシンのようなコンピュータ上で迅速に行われることができる。本発明の方法はまた、少ないメモリを使用することにより、たとえば、生物学的データの解析のために構成要素に関する共有情報（限界情報ではなく）を使用する従来技術の方法よりも迅速に行われることを可能にする。

以下、多クラスロジスティック回帰モデルに関する第１の実施形態を説明する。
Ａ．多クラスロジスティック回帰モデル
この実施形態の方法は、訓練サンプルを予め規定されたグループに分類することのできる構成要素の部分集合を識別するために訓練サンプルを使用する。続いて、たとえば、臨床テストのようなテストに対して構成要素の部分集合の情報を使用して、サンプルに病気クラスのようなグループに分類することができる。たとえば、ＤＮＡマイクロアレイの構成要素の部分集合は、臨床サンプルを、たとえば、健康または病気のような臨床的に適切なクラスにグループ化するために使用されることができる。

このようにして、本発明は、ある特定の訓練サンプルがある特定のグループに属しているか否かを識別するために使用されることのできる、好ましくは最少数の構成要素を識別する。最少数の構成要素はそのグループの“診断に役立ち”、あるいはグループの区別を可能にする。本質的に、システムから発生されたデータの全てから、本発明の方法は、特定のグループをテストするために使用されることのできる最少数の構成要素を識別することを可能にする。この方法によってこれらの構成要素が識別されると、その構成要素は新しいサンプルをグループに分類するために将来使用されることができる。本発明の方法は統計的方法を使用して、サンプルが属しているグループを正しく識別するのに必要とされない構成要素を除去することが好ましい。

サンプルは、予め定められた分類に基づいてサンプルグループ（または“クラス”）にグループ化される。この分類は、訓練サンプルがグループ化される任意の所望の分類であってよい。この分類は、たとえば、訓練サンプルが得られたのは白血病細胞または健康な細胞のどちらからかでもよく、あるいは訓練サンプルが得られたのはある症状を有する患者または有しない患者のどちらの血液からかでもよく、もしくは訓練サンプルが得られたのは正常な細胞と比較されるいくつかのタイプの癌の１つからの細胞からであってもよい。

１実施形態において、入力データは、ｎ個の訓練サンプルおよびｐ個の構成要素を有するｎ×ｐデータマトリックスＸ＝（ｘ_ij）に編成される。典型的に、ｐはｎよりはるかに大きくなる。

別の実施形態において、データマトリックスＸは、線形予測子の代りの予測子としてＸの平滑な関数を得るためにｎ×ｎカーネルマトリックスＫによって置換されてもよい。カーネルマトリックスＫの一例はｋ_ij＝ｅｘｐ（−０．５＊（ｘ_i−ｘ_j）^t（ｘ_i−ｘ_j）／σ²）であり、ここでｘの添字はマトリックスＸ中の行の数を示している。理想的には、これらの平滑な関数の希薄な（ｓｐａｒｓｅ）表現を提供するＫの列の部分集合が選択される。カーネルマトリックスＫの別の例は以下の表２に示されている。（表３は少しでも必要か。）
各サンプルクラス（グループ）にはあるクラスラベルｙ_iが関連付けられており、ここでｙ_i＝ｋ，ｋ∈｛１，・・・，Ｇ｝であり、これは、訓練サンプルがＧのサンプルクラスのいずれに属しているかを示す。エレメントｙ_iを有するｎ×ｌベクトルをｙと書く。そのベクトル：

を仮定すると、インジケータ変数が規定されることができる：
ｙ_i＝ｇならば、ｅ_ig＝１
その他場合、ｅ_ig＝０ （１Ａ）

１実施形態において、構成要素加重は、ベイズの統計的モデル（KotzJohnson氏等による1983年の文献を参照されたい）を使用して推定される。この加重は、各訓練サンプルから発生されたデータを与えられた加重の事後分布を最大化することによって推定されることが好ましい。その結果、最大化されるべき目標関数は２つの部分から構成されたものとなる。第１の部分は尤度関数であり、第２の部分は、可能な限り常にゼロ加重が好ましいことを保証する加重に対する事前分布である。好ましい実施形態において、尤度関数は多クラスロジスティックモデルから導出される。尤度関数は、確率から計算されることが好ましい：

ここで、
ｐ_igは、入力データＸ_iを有する訓練サンプルがサンプルクラスｇ中に存在する確率であり、
ｘ_i ^Tβ_gは、構成要素加重β_gを有する訓練サンプルｉからの入力データから発生された線形結合であり、
ｘ_i ^TはＸのｉ番目の行に対する構成要素であり、β_gはサンプルクラスｇに対する構成要素加重のセットである。

上述したように、構成要素加重は典型的に、構成要素加重の大部分がゼロであるアプリオリ仮定を考慮した方式で推定される。

１実施形態において、式（２Ａ）中の構成要素加重β_gは、値の大部分がゼロであるが、そのサンプルが依然として正確に分類されることのできる方式で推定される。

１実施形態において、パラメータβ₁，・・・，β_G-1に対して特定されたプライアは以下の形態のものである：

また、ｐ（β_g|τ_g ²）はＮ（０，ｄｉａｇ｛τ_g ²｝）であり、

はジェフェリーのハイパープライアである（KotzおよびJohnson氏等の文献、1983）。

１実施形態において、尤度関数は式（８Ａ）中の形の

であり、ｙを与えられたβおよび

の事後分布は、

１実施形態において、第１の導関数は以下の式から決定される：

ここで、

は、サンプルクラスｇのメンバーシップおよびクラスｇの確率をそれぞれ示すベクトルである。

１実施形態において、第２の導関数は以下のアルゴリズムから決定される：

式６および式７は、次のように導出されてもよい：
（ａ）式（１Ａ），（２Ａ）および（３Ａ）を使用することにより、データの尤度関数を以下のように書くことができる：

（ｂ）式（８Ａ）の対数をとり、全てのｉに対して

であることを使用することによって、

が得られる。

（ｃ）式（９Ａ）をβ_gに関して微分することにより、次式が得られる：

ここで、

は、サンプルクラスｇのメンバーシップおよびクラスｇの確率をそれぞれ示すベクトルである。

（ｄ）式（９Ａ）の第２の微分はエレメント：

を有しており、ここで、
ｈ＝ｇならば、δ_hg＝１
そうでない場合、δ_hg＝０である。

尤度関数の事後分布を最大化する構成要素加重は、ＥステップおよびＭステップを含むＥＭアルゴリズムを使用して特定されてもよい。

典型的に、ＥＭアルゴリズムは：
（ａ）関数：

を使用して構成要素加重の事後分布の条件付き期待値を計算することによってＥステップを行い、
（ｂ）γの関数としてＱを最大化するために反復的な処理手続きを適用することによりＭステップを行い、このとき、

この場合、α’はステップ長であり、このとき０≦α’≦１である；
β_g＝Ｐ_gγ_g；
ここで、Ｐ_gはゼロと１との、Ｐ^T _gβ_gがβ_gのゼロでないエレメントを選択するようなマトリックスである；
γ＝（γ_g，ｇ＝１，・・・，Ｇ−１）

式（１２Ａ）は以下のように導出されることができる：
パラメータ推定値のセット：

と観察されたデータｙとを与えられた式（５Ａ）の条件付き期待値を計算する。

β（および

）の構成要素がゼロに設定される場合、すなわち、ｇ＝１，・・・，Ｇ−１、β_g＝Ｐ_gγ_gおよび

に対する場合を検討する。ここで、Ｐ_gはゼロと１との、Ｐ^T _gβ_gがβ_gのゼロでないエレメントを選択するようなマトリックスである。以下において、γ＝（γ_g，ｇ＝１，・・・，Ｇ−１）と書く。γ_gは実際にβ_gの構成要素の部分集合であることに注意されたい。それらを使用して、表記を可能な限り単純なものにしておく。

γを含まない項を無視し、式（４Ａ），（５Ａ），（９Ａ）を使用することにより、次式が得られる：

条件付き期待値は、式（４Ａ）を与えられた第１の原理から推定されることができることに注意されたい。

反復的な処理手続きは以下のように得られることができる：
式（１３Ａ）において必要とされる導関数を得るために、最初に、式（８Ａ）、（９Ａ）および（１０Ａ）から以下の式を得ることに注意されたい：

ここで、
ｇ＝ｈの場合、 δ_gh＝１
それ以外の場合、δ_gh＝０
および
Ｘ_g ^T＝Ｐ_g ^TＸ^T ， ｇ＝１，・・・，Ｇ−１ （１７Ａ）

好ましい実施形態において、反復的な処理手続きは、式（１３Ａ）において式（１６Ａ）のブロック対角線上の成分（ｄｉａｇｏｎａｌｓ）だけを使用することによって簡単化されることができる。これによって、ｇ＝１，・・・，Ｇ−１に対して、以下の関係が得られる：

式（１８Ａ）を再構成することにより、次式が得られる：

Ｙ_gの列の数に対してｐ（ｇ）を書くと、式（１９Ａ）は、非常に大きいものとなる可能性のあるｐ（ｇ）×ｐ（ｇ）マトリックスの逆数を必要とする。これは、ｐ（ｇ）＞ｎに対しては、以下のように表すことによりｎ×ｎマトリックスに減少されることができる：

ここで、Ｚ_g＝Δ^1/2 _ggＹ_gである。ｐ（ｇ）＜ｎのときは、式（１９Ａ）が使用され、ｐ（ｇ）≧ｎであるときには、式（１９Ａ）に式（２０Ａ）を代入したものが使用されることが好ましい。

好ましい実施形態において、ＥＭアルゴリズムは以下のように行われる：
１．ｎ＝０、Ｐ_g＝Ｉと設定し、

に対する初期値を選択する。これは、ｘ_iに関するｌｏｇ（ｐ_ig／ｐ_iG）のリッジ（ｒｉｄｇｅ）回帰により行われ、ここでｐ_igは、全ての確率が合計すると１になるという制約を条件として、グループｇにおいて観察されたものに対して１に近くなり、その他の場合は、わずかな量＞０となるように選択される。

２．Ｅステップを行う。すなわち、

を推定する。

３．ｔ＝０を設定する。ｇ＝１，・・・，Ｇ−１に対して、以下を計算する：
（ａ）ｐ（ｇ）≧ｎであるとき、式（１９Ａ）に式（２０Ａ）を代入したものを使用して、δ^t _g＝γ^t+1 _g−γ^t _g
（ｂ）δ^t＝（δ^t _g，ｇ＝１，・・・，Ｇ−１）と書く。式（１２Ａ）をα^tの関数として最大化する（あるいは単に増加させる）

中のα^tの値を見つけるためにラインサーチを行う。
（ｃ）

およびｔ＝ｔ＋１を設定する。
収斂するまでステップ（ａ）乃至（ｃ）を繰り返す。
これによって、たとえば、現在のＱ関数をγの関数として最大化するγ^*n+1が生成される。
ｇ＝１，・・・，Ｇ−１に対して、

を決定する。

ここで、ε＜＜１、たとえば１０^-5である。ｉ∈Ｓ_gおよび

に対してΒ_ig＝０となるようにＰ_gを定める。
このステップにより、小さい係数を有する変数がモデルから除去される。

４．ｎ＝ｎ＋１を設定し、収斂するまで２を行う。

以下、順序付けられたカテゴリー的ロジスティック回帰に関する第２の実施形態を説明する。
Ｂ．順序付けられたカテゴリー的モデル
この実施形態の方法は訓練サンプルを使用して、あるテストサンプルが特定のクラスに属しているか否かを決定するために使用されることのできる構成要素の部分集合を識別することができる。たとえば、マイクロアレイ解析を使用して組織生検サンプルを評価するために遺伝子を識別するために、正常な組織、良性の組織、局所化した腫瘍および転移した腫瘍組織のような、前に病気の重大度が増加するまたは減少するクラスに順序付けられた組織からの一連のサンプルからのマイクロアレイデータは、訓練サンプルに関連した病気の重大度を示すことのできる構成要素の部分集合を識別するために訓練サンプルとして使用される。その後、構成要素の部分集合は、以前分類されていなかったテストサンプルが正常、良性、局所化した腫瘍、または転移した腫瘍組織のいずれとして分類されるかを決定するために後続的に使用されることができる。したがって、構成要素の部分集合は、テストサンプルが順序付けられたクラスセット内の特定のクラスに属するかどうかの診断に役立つ。構成要素の部分集合がひとたび識別されると、将来の診断処理手続きにおいて、あるサンプルが順序付けられたクラスのどれに属しているのかを決定するのに、その構成要素の部分集合がテストされるだけでよいことが認識される。

本発明の方法はとくに、非常に大量のデータの解析に適している。典型的に、テストサンプルから得られた大量のデータセットは非常に多種多様であり、また、訓練サンプルから得られたものとは著しく異なっていることが多い。本発明の方法は訓練サンプルから発生された非常に大量のデータからの構成要素の部分集合を識別することが可能であり、また、この方法によって識別された構成要素の部分集合は、その後、テストサンプルから発生されたデータが同じクラスに属している訓練サンプルから発生されたデータと比較して大幅に異なっている場合でも、そのテストサンプルを分類するために使用されることができる。したがって、本発明の方法は、データが品質の低いものであっても、および、または順序付けられた同じクラスのサンプル間に大きいばらつきがある場合であってさえも、サンプルを正しく分類する可能性が高い構成要素の部分集合を識別することができる。

最少数の構成要素は、その特定の順序付けられたクラスに対する“予測的な”ものである。本質的に、システムから発生された全てのデータから、本発明の方法は訓練データを分類するために使用されることのできる最少数の構成要素を識別することができる。この方法によってこれらの構成要素が識別されると、その構成要素は、将来、テストサンプルを分類するために使用されることができる。本発明の方法は統計的方法を使用して、サンプルを順序付けられたクラスのメンバーであるクラスに正しく分類するのに必要でない構成要素を除去することが好ましい。

以下において、Ｎ個のサンプルが存在し、ｙ、ｚおよびμのようなベクトルはｉ＝１，・・・，Ｎに対して構成要素ｙ_i、ｚ_iおよびμ_iを有している。ベクトル乗算および除算が構成要素に関して定義され、ｄｉａｇ｛・｝は、対角線上の成分がアーギュメントに等しい対角行列を表す。ユークリッドのノルムを表すために｜｜・｜｜もまた使用される。

Ｎ個の観察ｙ_iが存在していることが好ましく、ここでｙ_iは整数値１，・・・，Ｇである。これらの値は、たとえば、病気の重大度のようなある方式で順序付けられたクラスを表す。Ｎ行およびｐ列を有するマトリックスＸに配列されたあるセットのコバリエイト（たとえば、遺伝子発現値等の変数）が各観察と関連付けられており、ここでＮはサンプルであり、ｐは構成要素である。表記ｘ_i ^TはＸのｉ番目の行を表す。個体（サンプル）ｉは、π_ik＝π_k（ｘ_i）で与えられるクラスｋに属する確率を有している。

累積確率を規定する：

γ_ikは観察ｉがｋ以下のインデックスを有するクラスに属する確率であることに注意されたい。Ｃを、
観察ｉがクラスｊ中のものである場合、ｃ_ij＝１
そうでない場合、ｃ_ij＝０
により与えられるエレメントｃ_ijを有するｎ×ｐマトリックスとし、Ｒを

により与えられるエレメントｒ_ijを有するｎ×Ｐマトリックスとする。
これらは、行内におけるＣの列の累積的な和である。

独立した観察（サンプル）に関して、データの尤度は以下のように書かれることができる：

また、対数尤度（ｌｏｇ（Ｌ））ｌは以下のように書かれることができる：

ここにおいて、以下のように、ｋ＝２，・・・，Ｇに対して、連続率（ｃｏｎｔｉｎｕａｔｉｏｎ ｒａｔｉｏ）モデルが採用されてもよい：

McCullaghおよびNelder氏等の文献(1989)、McCullagh氏の文献(1980)ならびにこの明細書における説明を参照されたい。

であることに注意されたい。

尤度は、反応ベクトルｙおよびコバリエイトマトリックスＸを有する対数回帰尤度に等しい

ここで、Ｉ_G-1は（Ｇ−１）×（Ｇ−１）識別子マトリックスであり、ｌ_G-1は１の（Ｇ−１）×ｌベクトルである。
ここにおいて、ｖｅｃ｛ ｝はマトリックスをとり、ベクトル行×行を成す。

典型的に、上述したように、構成要素加重は、構成要素加重の大部分がゼロであるというアプリオリ仮定を考慮したやり方で推定される。

Figueiredo氏の文献(2001)にならって、冗長な変数（コバリエイト）を除去するために、ハイパーパラメータ（ｈｙｐｅｒｐａｒａｍｅｔｅｒｓ）のｐ×ｌベクトルを導入することにより、あるプライアがパラメータβ^*に対して特定される。

構成要素加重に対して特定されたプライアは、以下の形のものであることが好ましい：

ここで、KotzおよびJohnson氏等の文献(1983)に示されているように、ｐ（β^*｜ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、

はジェフェリーのプライアである。θ＝（θ₂，・・・θ_G）^Tのエレメントは非有益プライアを有している。

尤度関数を

と書くことにより、ベイズのフレームワークにおいて、ｙを与えられたβ^*、θおよびｖの事後分布は、以下のようになる：

Ｖを失われたデータのベクトルとして処理することにより、式（６Ｂ）を最大化して局所極大であるβ^*およびθの事後推定値を生成するために、ＥＭアルゴリズム［Dempster氏の文献(1977)］のような反復的アルゴリズムが使用可能であることが好ましい。上記のプライアは最大である事後推定値が希薄（ｓｐａｒｓｅ）になる傾向がある、すなわち、多数のパラメータが冗長である場合、β^*の多くの構成要素がゼロになるようなものである。以下においてβ^T＝（θ^T，β^*T）であり、ｄｉａｇ（）は対角行列を表すことが好ましい：
上記の順序付けられたカテゴリーに関して、

であることが証明されることができる。

上述したように、事後分布を最大化する構成要素加重は、反復的な処理手続きを使用して決定されることができる。構成要素および構成要素加重の事後分布を最大化する反復的な処理手続きは、たとえば、Dempster氏の文献(1977)に記載されているようなＥＭアルゴリズムであることが好ましい。ＥＭアルゴリズムは以下のようにして行われることが好ましい：
１．ｎ＝０，Ｓ₀＝｛１，２，・・・，ｐ｝，φ⁽⁰⁾ および（たとえば）ε＝１０^-5を設定する。レギュラリゼーション（regularization）パラメータκを１よりはるかに大きい値、たとえば、１００に設定する。これは、１／κ²を第２の微分マトリックスの第１のＧ−１対角エレメントに以下のＭステップで加算することに対応する。

ｐ≦Ｎである場合、

によって初期値β*を計算し、
ｐ＞Ｎである場合、

によって初期値β*を計算する。ここでリッジパラメータλは、０＜λ≦１を満足させ、ζは小さく、ロジット（ｌｏｇｉｔ）リンク関数ｇがｙ＋ζにおいて明確に定義されるように選択される。

２．

と定義し、Ｐ_nを、β⁽ⁿ⁾のゼロでないエレメントγ⁽ⁿ⁾が
γ⁽ⁿ⁾＝Ｐ_n ^Tβ⁽ⁿ⁾ ， β⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾
γ＝Ｐ_n ^Tβ ， β＝Ｐ_nγ
を満足させるような０と１とのマトリックスとする。

ｉ≧Ｇである場合、ｗ_βi＝１
それ以外の場合、ｗ_βi＝０
であるようにｗ_β＝（ｗ_βi，ｉ＝１，ｐ）と定義し、ｗ_γ＝Ｐ_nｗ_βとする。

３．

を計算することによってＥステップを行う。ここで、ｌはｙの対数尤度関数である。

β＝Ｐ_nγおよびβ⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾を使用することにより、式（１１Ｂ）は以下のように書かれることができる：

４．Ｍステップを行う。これは、以下のようにＮｅｗｔｏｎ Ｒａｐｈｓｏｎ反復により行われることができる。γ₀＝γ⁽ⁿ⁾ を設定し、ｒ＝０，１，２，・・・に対してγ_r+1＝γ_r＋α_rδ_rであり、ここでα_rは、
Ｑ（γ_r+1｜γ⁽ⁿ⁾）＞Ｑ（γ_r｜γ⁽ⁿ⁾）
を保証するようにラインサーチアルゴリズムによって選択される。

ｐ≦Ｎに対しては、以下の式を使用する：

ｐ＞Ｎに対しては、以下の式を使用する：

Ｖ_rおよびｚ_rは前のように定義されている。

γ^*を、ある収斂基準が満足されたとき、たとえば、以下のようなときのγ_rの値とする：
｜｜γ_r−γ_r+1｜｜＜ε（たとえば、１０^-5）
５．

と定義する。ここで、ε₁は、たとえば、ｌｅ−５等の小さい定数である。ｎ＝ｎ＋１を設定する。

６．収斂をチェックする。｜｜γ^*−γ⁽ⁿ⁾｜｜＜ε₂である場合、中止し、ここでε₂は適切に小さく、その他の場合、上記のステップ２に進む。

確率を再生する。

パラメータβの推定値が得られると、ｉ＝１，・・・，Ｎおよびｋ＝２，・・・，Ｇに対して、

を計算する。

確率を得るために、以下の反復が使用されることが好ましい：

また、ｉ＝１，・・・，Ｎに対して、この確率の合計は１になる。

１実施形態において、行ｘ_i ^Tを有するコバリエイトマトリックスＸは、ｉｊ番目のエレメントｋ_ijを有し、あるカーネル関数κに対してｋ_ij＝κ（ｘ_i−ｘ_j）であるマトリックスＫによって置換されることができる。このマトリックスはまた１のベクトルにより増加されることができる。ある例のカーネルが以下の表１に与えられている［Evgeniou氏他による文献(1999)参照］。

カーネル関数 κ（ｘ−ｙ）に対する公式
ガウスのラジアル基礎関数 ｅｘｐ（−｜｜ｘ−ｙ｜｜²／ａ），ａ＞０
逆マルチクアドリック(multiquadric) （｜｜ｘ−ｙ｜｜²＋ｃ²）^-1/2
マルチクアドリック （｜｜ｘ−ｙ｜｜²＋ｃ²）^1/2
薄板スプライン ｜｜ｘ−ｙ｜｜²ⁿ⁺¹
｜｜ｘ−ｙ｜｜²ⁿｌｎ（｜｜ｘ−ｙ｜｜）
マルチレーヤ知覚 適切なθに対して、ｔａｎｈ（ｘ’ｙ−θ）
度数ｄの多項式 （１＋ｘ’ｙ）^d
Ｂスプライン Ｂ_2n+1（ｘ−ｙ）
三角多項式 ｓｉｎ（（ｄ＋１／２）（ｘ−ｙ））/ｓｉｎ（（ｘ−
ｙ）／２）
表１：カーネル関数の例

表１において、最後の２つのカーネルは１次元、すなわち、Ｘが列を１つしか有しない場合に対して１次元であることが好ましい。多変量バージョンはこれらのカーネル関数の積から得られることができる。Ｂ_2n+1の定義はDe Boor氏による文献（1978）に記載されている。カーネル関数を使用した結果、推定された確率はコバリエイトＸの平滑な（線形の変換とは対照的に）関数となる。このようなモデルにより、データに対する実質的に良好な適合が得られる。

以下、一般化された線形モデルに関する第３の実施形態を説明する。
Ｃ．一般化された線形モデル
この実施形態の方法は、あるサンプルの特徴を予測することのできる構成要素の部分集合を識別するために訓練サンプルを使用する。後続的に、構成要素の部分集合の知識は、重要な特徴の未知の値を予測するための、たとえば、臨床テスト等の、テストのために使用されることができる。たとえば、ＤＮＡマイクロアレイの構成要素の部分集合は、たとえば、血液グルコースレベル、白血球細胞数、腫瘍の寸法、腫瘍の成長速度または生存時間等の、臨床的に関連した特徴を予測するために使用されることができる。

このようにして、本発明は、特定のサンプルに対する特徴を予測するために使用されることのできる最少数の構成要素を識別することが好ましい。最少数の構成要素が、その特徴に対する“予測的な”ものである。本質的に、システムから発生された全てのデータから、本発明の方法は、特定の特徴を予測するために使用されることのできる最少数の構成要素を識別することができる。この方法によってこれらの構成要素が識別されると、その構成要素は、将来、新しいサンプルに対する特徴を予測するために使用されることができる。本発明の方法は統計的方法を使用して、サンプルに対する特徴を正しく予測するのに必要でない構成要素を除去することが好ましい。

本発明は、訓練サンプルから発生されたデータの構成要素の線形結合の構成要素加重が、ある訓練サンプルに対する特徴を予測するのに必要とされない構成要素を除去する方式で推定されることができることを見出している。その結果、訓練セット中のサンプルに対する特徴を正しく予測することのできる構成要素の部分集合が識別される。したがって、本発明の方法は、たとえば、重要な量等の、訓練サンプルに対する特徴を正しく予測することのできる比較的少数の構成要素を大量のデータから識別することを可能にする。

特徴は任意の重要な特徴であることができる。１実施形態において、特徴は量または測定値である。別の実施形態において、それらはあるグループのインデックス番号であることができ、ここにおいて、サンプルは予め定められた分類に基づいて２つのサンプルグループ（または“クラス”）にグループ化されている。この分類は、訓練サンプルがグループ化される任意の所望の分類であることができる。たとえば、分類は、訓練サンプルが白血病細胞からのものか、あるいは健康な細胞からのものかであってもよく、もしくは訓練サンプルが得られたのはある症状を有する患者または有しない患者かでもよく、あるいは訓練サンプルが得られたのは正常な細胞と比較されるいくつかのタイプの癌の１つからの細胞からであってもよい。別の実施形態において、特徴は、特定の患者が少なくとも所定の日数生存していることを示す検定された生存時間であってもよい。別の実施形態において、量は、たとえば、血圧等の、測定可能なサンプルの連続的に変化する任意の特徴であってもよい。

１実施形態において、データは量ｙ_iであってもよく、ここでｉ∈｛１，・・・，Ｎ｝。エレメントｙ_iを有するＮ×１ベクトルをｙと書く。構成要素加重のｐ×１パラメータベクトルβ（それらの多くはゼロであると期待される）およびパラメータφのｑ×１ベクトル（ゼロでないと期待される）が規定される。ｑはゼロであることができる（すなわち、ゼロでないと期待されたパラメータのセットは空であることができる）ことに注意されたい。

１実施形態において、入力データはＮ個のテスト訓練サンプルとｐ個の構成要素とを有するＮ×ｐデータマトリックスＸ＝（ｘ_ij）に編成される。典型的に、ｐはＮよりはるかに大きい。

別の実施形態において、データマトリックスＸは、線形予測子の代りの予測子としてＸの平滑な関数を得るためにＮ×ＮカーネルマトリックスＫによって置換されることができる。カーネルマトリックスＫの一例はｋ_ij＝ｅｘｐ（−０．５＊（ｘ_i−ｘ_j）^t（ｘ_i−ｘ_j）／σ²）であり、ここでｘの添字はマトリックスＸ中の行番号を示す。理想的には、これらの平滑な関数のスパース表現を提供するＫの構成要素の部分集合が選択される。

典型的に、上述したように、構成要素加重は、構成要素加重の大部分がゼロであるというアプリオリ仮定を考慮した方式で推定される。

１実施形態において、構成要素加重に対して特定されたプライアは以下の形態のものである：

また、ｐ（β|ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、

はジェフェリーのプライアである（KotzおよびJohnson氏等の文献、1983）。φに対する非有益プライアが特定されることが好ましい。

尤度関数は、データの分布に基づいてそのデータに適合するモデルを規定する。尤度関数は、一般化された線形モデルから得られることが好ましい。たとえば、尤度関数：

は、たとえば、ＮｅｌｄｅｒおよびＷｅｄｄｅｒｂｕｒｎ氏等の文献(1972)に記載されているもののような、一般化された線形モデル（ＧＬＭ）に対して適した形態であることができる。尤度関数は、以下の形態のものであることが好ましい：

ここで、ｙ＝（ｙ₁，・・・，ｙ_n）^Tおよびａ_i（φ）＝φ／ｗ_iであり、ｗ_iは既知の加重の固定されたセットであり、φは単一のスケールパラメータである。

尤度関数は、以下のように特定されることが好ましい：

各観察はコバリエイトのセットｘ_iおよび線形予測子η_i＝ｘ_i ^Tβを有している。ｉ番目の観察の平均とそれの線形予測子との間の関係は、リンク関数η_i＝ｇ（μ_i）＝ｇ（ｂ’（θ_i））によって与えられる。このリンクの逆数はｈによって表される。すなわち、
μ_i＝＝ｂ’（θ_i）＝ｈ（η_i）
である。

スケールパラメータに加えて、一般化された線形モデルは４つの構成要素によって特定されることができる：
・尤度または（スケールされた）逸脱関数
・リンク関数
・リンク関数の微分
・分散関数。

一般化された線形モデルのいくつかの共通した例は、以下の表２に与えられている。

別の実施形態において、尤度関数は多クラスロジスティックモデルから得られる。

別の実施形態において、リンク関数および分散関数だけが規定される擬似尤度モデルが特定される。いくつかの例において、このような仕様は結果的に上記の表中のモデルを生じさせる。別の例では、分散は特定されない。

１実施形態において、ｙを与えられたβφおよびｖの事後分布は、以下の式を使用して推定される：

１実施形態において、ｖは失われたデータのベクトルとして処理されてもよく、反復的な処理手続きは局所極大であるβの事後推定値を生成するように式（２Ｃ）を最大化するために使用されてもよい。式（５Ｃ）のプライアは、最大である事後推定値が希薄になるようなもの、すなわち、非常に多くのパラメータが冗長である場合にβの多くの構成要素がゼロになるようなものである。

上述したように、事後分布を最大化する構成要素加重は、反復的処理手続きを使用して決定されることができる。構成要素および構成要素加重の事後分布を最大化する反復的処理手続きは、たとえば、Dempster氏他の文献(1977)に記載されているようなＥＭアルゴリズムである。

１実施形態において、ＥＭアルゴリズムは以下のステップを含んでいる：
（ｃ）ｎ＝０，Ｓ0＝｛１，２，・・・，ｐ｝を設定することによりアルゴリズムを初期化し、φ⁽⁰⁾ ，β*を初期化し、たとえばε＝１０^-5のようなεに対する値を適用する；
（ｄ）ｉεＳ_nのときは、β_i ⁽ⁿ⁾ ＝β_i ^*，
それ以外のときは、β_i ⁽ⁿ⁾ ＝０ （５Ｃ）
と定義し、ｐnを、β(n)のゼロでないエレメントγ(n)が
γ⁽ⁿ⁾＝Ｐ_n ^Tβ⁽ⁿ⁾ ， β⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾
γ＝Ｐ_n ^Tβ ， β＝Ｐ_nγ
を満足させるような０と１とのマトリックスとする。
（ｅ）以下の関数を使用して構成要素加重の事後分布の条件付き期待値を計算することにより推定（Ｅ）ステップを行う：

ここで、ｌはｙの対数尤度関数である。

β＝Ｐ_nγおよびβ⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾を使用することにより、この式（６Ｃ）は以下のように書かれることができる：

（ｆ）反復的な処理手続きを適用してＱをγの関数として最大化することにより、最大化（Ｍ）ステップを行う。このとき、γ₀＝γ⁽ⁿ⁾ であり、ｒ＝０，１，２，・・・に対してγ_r+1＝γ_r＋α_rδrであり、ここでα_rは、
Ｑ（γ_r+1｜γ⁽ⁿ⁾，φ⁽ⁿ⁾）＞Ｑ（γ_r｜γ⁽ⁿ⁾，φ⁽ⁿ⁾）
を保証するようにラインサーチアルゴリズムによって選択され、また、

（ｇ）γ*を、ある収斂基準が満足されたとき、たとえば、以下のようなときのγrの値とする：
｜｜γr−γr+1｜｜＜ε（たとえば、１０^-5）；
（ｈ）

と定義する。ここで、ε₁は、たとえば、ｌｅ−５等の小さい定数である。
（ｉ）ｎ＝ｎ＋１を設定し、φ⁽ⁿ⁺¹⁾＝φ⁽ⁿ⁾＋κ_n（φ^*−φ⁽ⁿ⁾）を選択し、ここでφ^*は、

を満足させ、κ_nは０＜κ_n≦１となるような減衰定数である；
（ｊ）収斂をチェックする。｜｜γ^*−γ⁽ⁿ⁾｜｜＜ε₂である場合、ε₂は適切に小さく、停止し、その他の場合、上記のステップ（ｂ）に進む。

別の実施形態において、最大化ステップにおけるステップ（ｄ）は、

をその期待値：

で置換することにより推定されることができる。これは、データのモデルが一般化された線形モデルであるときに好ましい。

一般化された線形モデルについて、期待値：

は以下のように計算されることができる：

ここで、ＸはＮ×ｐマトリックスであり、ｉ番目の行がｘ_i ^Tであり、

これは、以下のように書かれることができる：

ＥＭアルゴリズムは、以下のステップを含んでいることが好ましい：
（ａ）ｎ＝０，Ｓ0＝｛１，２，・・・，ｐ｝，φ(0) を設定することによりアルゴリズムを初期化し、たとえばε＝１０^-5のようなεに対する値を適用し、
ｐ≦Ｎである場合、

によって初期値β*を計算し、
ｐ＞Ｎである場合、

によって初期値β*を計算する。ここでリッジパラメータλは、０＜λ≦１を満足させ、ζは小さく、リンク関数ｇがｙ＋ζにおいて明確に定義されるように選択される。

（ｂ）ｉεＳ_nのときは、β_i ⁽ⁿ⁾ ＝β_i ^*，
それ以外のときは、β_i ⁽ⁿ⁾ ＝０
と定義し、ｐnを、β(n)のゼロでないエレメントγ(n)が
γ⁽ⁿ⁾＝Ｐ_n ^Tβ⁽ⁿ⁾ ， β⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾
γ＝Ｐ_n ^Tβ ， β＝Ｐ_nγ
を満足させるような０と１とのマトリックスとする。
（ｃ）以下の関数を使用して構成要素加重の事後分布の条件付き期待値を計算することにより推定（Ｅ）ステップを行う：

ここで、ｌはｙの対数尤度関数である。

β＝Ｐ_nγおよびβ⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾を使用することにより、この式（１６Ｃ）は以下のように書かれることができる：

（ｄ）たとえば、Ｎｅｗｔｏｎ Ｒａｐｈｓｏｎ反復等の、反復的な処理手続きを適用してＱをγの関数として最大化することにより、最大化（Ｍ）ステップを行う。このとき、γ₀＝γ⁽ⁿ⁾ であり、ｒ＝０，１，２，・・・に対してγ_r+1＝γ_r＋α_rδ_rであり、ここでα_rは、
Ｑ（γ_r+1｜γ⁽ⁿ⁾，φ⁽ⁿ⁾）＞Ｑ（γ_r｜γ⁽ⁿ⁾，φ⁽ⁿ⁾）
を確実にするようにラインサーチアルゴリズムによって選択され、また、
ｐ≦Ｎである場合、

添字ｒはこれらの量がμ＝ｈ（ＸＰ_nγ_r）で推定されることを示す。

ｐ＞Ｎである場合、

Ｖ_rおよびｚ_rは前のように定義されている。

γ*を、ある収斂基準が満足されたとき、たとえば、以下のようなときのγ_rの値とする：
｜｜γ_r−γ_r+1｜｜＜ε（たとえば、１０^-5）
１）

と定義する。ここで、ε₁は、たとえば、ｌｅ−５等の小さい定数である。

ｎ＝ｎ＋１を設定し、φⁿ⁺¹＝φⁿ＋κ_n（φ^*−φⁿ）を選択し、ここでφ^*は、

を満足させ、κ_nは０＜κ_n≦１となるような制動定数である。いくつかの場合において、スケールパラメータは知られており、あるいはこの式は明示的に解かれることができるため、φに対する更新方程式が得られることに注意されたい。

上記の実施形態は、Wedderburn氏による文献(1974)ならびにMcCullaghおよびNelder氏等による文献(1983)に記載されている擬似尤度方法を含むように拡張されることができる。このようなの実施形態において、詳細に上述された同じ反復的な処理手続きは適切であるが、しかし尤度が上記および、たとえば、McCullaghおよびNelder氏等による文献(1983)中の表８．１に示されている擬似尤度によって置換される。１実施形態において、スケールパラメータφに対する修正された更新方法が存在する。これらのモデルを規定するには、分散関数τ²、リンク関数ｇ、およびリンク関数の導関数：

の仕様が必要である。

これらがひとたび規定されると、上記のアルゴリズムは適用されることができる。１実施形態において、擬似尤度モデルに対して、上記のアルゴリズムのステップ５は、スケールパラメータが以下を計算することにより更新されるように修正される：

ここで、μおよびτはβ^*＝Ｐ_nγ^*で推定される。この更新は、モデル中のパラメータの数ｓがＮより小さいときに行われることが好ましい。Ｎ−ｓの除数は、ｓがＮよりはるかに小さいときに使用されることができる。

別の実施形態において、一般化された線形モデルおよび擬似尤度モデルの両者に対して、行ｘ_i ^Tを有するコバリエイトマトリックスＸは、ｉｊ番目のエレメントｋ_ijを有し、あるカーネル関数κに対してｋ_ij＝κ（ｘ_i−ｘ_j）であるマトリックスＫによって置換されることができる。このマトリックスはまた１のベクトルにより増加されることができる。ある例のカーネルが以下の表３に与えられている［Evgeniou氏他による文献(1999)参照］。

カーネル関数 κ（ｘ−ｙ）に対する公式
ガウスのラジアル基礎関数 ｅｘｐ（−｜｜ｘ−ｙ｜｜²／ａ），ａ＞０
逆マルチクアドリック(multiquadric) （｜｜ｘ−ｙ｜｜²＋ｃ²）^-1/2
マルチクアドリック （｜｜ｘ−ｙ｜｜²＋ｃ²）^1/2
薄板スプライン ｜｜ｘ−ｙ｜｜²ⁿ⁺¹
｜｜ｘ−ｙ｜｜²ⁿｌｎ（｜｜ｘ−ｙ｜｜）
マルチレーヤ知覚 適切なθに対して、ｔａｎｈ（ｘ’ｙ−θ）
度数ｄの多項式 （１＋ｘ’ｙ）^d
Ｂスプライン Ｂ_2n+1（ｘ−ｙ）
三角多項式 ｓｉｎ（（ｄ＋１／２）（ｘ−ｙ））/ｓｉｎ（（ｘ−
ｙ）／２）
表３：カーネル関数の例

表３において、最後の２つのカーネルは１次元である、すなわち、Ｘが列を１つしか有しない場合に対して１次元である。多変量バージョンはこれらのカーネル関数の積から得られることができる。Ｂ_2n+1の定義はDe Boor氏による文献（1978）に記載されている。一般化された線形モデルまたは擬似尤度モデルのいずれかにおいてカーネル関数を使用した結果、推定された確率はコバリエイトＸの平滑な（線形の変換とは対照的に）関数である平均値が得られる。このようなモデルにより、データに対する実質的に良好な適合が得られる。

以下、比例ハザードモデルに関する第４の実施形態を説明する。
Ｄ．比例ハザードモデル
この実施形態の方法は、規定された事象（たとえば、死、回復等）がある時間期間内に発生する確率に影響を与えることのできる構成要素の部分集合を識別するために訓練サンプルを使用することができる。訓練サンプルはシステムと、その訓練サンプルが得られたときからその事象が発生したときまで測定された時間とから得られる。統計的方法を使用して、複数の訓練サンプルから得られたデータにより時間を事象に関連付けることにより、事象までの時間の分布を予想することのできる構成要素の部分集合が識別されることができる。後続的に、その構成要素の部分集合の情報は、たとえば、死までの時間または病気の再発までの時間の統計的特徴を予測するために、たとえば、臨床テスト等のテストに対して使用されることができる。たとえば、システムの構成要素の部分集合からのデータは、ＤＮＡマイクロアレイから得られることができる。このデータは、たとえば、期待されるまたは中間的な患者の生存時間のような臨床的に関連した事象を予測するために、あるいはある症状の始まりまたは病気の再発を予測するために使用されることができる。

このようにして、本発明は、あるシステムのある事象までの時間の分布を予測するために使用されることのできる最少数の構成要素を識別することが好ましい。この最少数の構成要素が、そのある事象に対する時間の“予測的な”ものである。本質的に、本発明の方法は、システムから発生された全てのデータから、ある事象までの時間を予測するために使用されることのできる最少数の構成要素を識別することができる。この方法によってこれらの構成要素が識別されると、その構成要素は、将来、あるシステムのある事象までの時間の統計的特徴を新しいサンプルから予測するために使用されることができる。本発明の方法は統計的方法を使用して、あるシステムのある事象までの時間を正しく予測するのに必要でない構成要素を除去することが好ましい。

この明細書において使用されているように、“ある事象までの時間”とは、本発明の方法が適用されるサンプルを獲得したときからある事象の時間までの時間の尺度のことである。ある事象はどのような観察可能な事象であってもよい。システムが生物学的システムである場合、事象は、たとえば、あるシステムの故障までの時間、死までの時間、特定の症状の始まり、病気の始まりまたは再発、表現型または遺伝子型の変化、生化学的特徴の変化、組織体または組織の形態学的特徴の変化、行動の変化であることができる。

ある事象までの以前の時間からのサンプルがある事象までの特定の時間と関連付けられる。ある事象までの時間は、たとえば、サンプリングしたときから死までの時間が知られている患者から得られたデータ、換言すると“本当の”生存時間と、およびサンプルが最後に得られたときに生きていたことが唯一の情報である患者から得られたデータ、換言すると、特定の患者が少なくとも所定の日数生存していることを示す“検閲削除された”生存時間とから決定された時間であってもよい。

１実施形態において、入力データはＮ×ｐマトリックスＸ＝（ｘ_ij）に編成され、ここで、Ｎはテスト訓練サンプルであり、ｐは構成要素である。典型的に、ｐはＮよりはるかに大きい。

たとえば、マイクロアレイ実験からのＮ×ｐデータマトリックスＸ＝（ｘ_ij）を考慮する。ここで、Ｎは個人（またはサンプル）であり、各個人に対して同じｐの遺伝子が存在する。たとえば、生存時間等のある事象までの時間を示す変数ｙ_i（ｙ_i≧０）が各個人ｉ（ｉ＝１，２，・・・，Ｎ）と関連付けられることが好ましい。各個人に対して、その個人の生存時間が本当の生存時間であるか、あるいは検閲削除された生存時間であるかを示す変数が規定されることができることも好ましい。ｃ_iのような検閲削除インジケータは、
ｙ_iが検閲削除されていない場合、ｃ_i＝１、
ｙ_iが検閲削除された場合、ｃ_i＝０
を表す。

生存時間ｙ_iを有するＮ×１ベクトルは、

と書くことができ、検閲削除インジケータｃ_iを有するＮ×１ベクトルは、

と書くことができる。

典型的に、上述したように、構成要素加重は、その構成要素加重の大部分がゼロであるというアプリオリ仮定が何を考慮したかで推定される。

構成要素加重に対して特定されたプライアは以下の形態のものである：

ここで、β₁，β₂，・・・，β_nは構成要素加重であり、Ｐ（β_i｜τ_i）はＮ（０，τ_i ²）であり、Ｐ（τ_i）α１／τ_i ²はジェフェリーのプライアである（KotzおよびJohnson氏等の文献、1983）。

尤度関数は、データの分布に基づいてデータに適合したモデルを規定する。尤度関数は以下の形態のものであることが好ましい：

尤度関数によって規定されたモデルは、あるシステムのある事象までの時間を予測する任意のモデルであってもよい。
１実施形態において、尤度関数によって規定されたモデルは、コックスの比例ハザードモデルである。コックスの比例ハザードモデルはコックス氏によって導入され(1972年)、好ましいことに、生存データに対する回帰モデルとして使用されることができる。コックスの比例ハザードモデル：

は、構成要素と関連付けられた（説明的な）パラメータのベクトルである。本発明の方法は、データ：

を与えられたコックスの比例ハザードモデルに対して、パラメータ：

からの非常に少数のものの選択（および除去）を行うことが好ましい。

コックスの比例ハザードモデルの適用は、同じ生存時間に対して種々のデータがシステムから得られる状況において、換言すると、同点の生存時間が生じた場合に問題となる可能性が高い。同点の生存時間は予備処理ステップを実施され、それによって特有の生存時間にされることができる。提案された予備処理により、コックスの比例ハザードモデルが後続的に適用されたときに同点の生存時間に関する問題が回避されるため、必然的に続いて発生するアルゴリズムは簡単なものになる。

生存時間の予備処理は、非常に小量のランダム雑音を追加することにより適用される。この処理手続きは、同点の時間のセットを選び、ゼロ平均と、分類された生存時間の間の最小の非ゼロ距離に比例した分散とを有する正規分布から引き出されたランダムな量を同点の時間のセット内の各同点時間に付加することであることが好ましい。このような予備処理により、生存時間の厳しい摂動を生ぜずに同点時間の除去が行われる。

予備処理により、異なった生存時間が発生する。好ましいことに、これらの時間は、以下に表されるように大きさが増加していく順に順序付けられることができる：

Ｘの行の再配列であるＮ×ｐマトリックスをＺで表し、ここで、Ｚの行の順序付けは、

の順序付けにより誘発された順序付けに対応し、また、マトリックスＺのｊ番目の行はＺ_jで表される。ｄを、

を順序付けるために必要とされる同じ置換によりｃを順序付けた結果とする。

同点の生存時間に対する予備処理を考慮し、生存データ解析に関する標準テキスト（たとえば、CoxおよびOakes氏等の文献，1984）への参照を行った後、確率ハザードモデルに対する尤度関数を以下のように書くことができることが好ましい：

尤度の対数（すなわち、ｌ＝ｌｏｇ（Ｌ））は以下のように書くことができることが好ましい：

ここで、
ｊ＜ｉである場合、ζ_i,j＝０
ｊ≧ｉである場合、ζ_i,j＝１
このモデルは、生存分布のパラメトリック形態が特定されておらず、生存時間の通常の特性だけが使用されることが好ましい（リスクセットの決定において）ため、非パラメトリックであることに注意されたい。これは非パラメトリックなケースであるため、

は必要とされない（すなわち、ｑ＝０）。

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

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

は生存密度関数の関数形態に関連付けられたパラメータのベクトルであることが好ましい。好ましくは、本発明の方法は、データ：

を与えられたパラメトリック生存モデルに対して、パラメータ：

および

の推定値からの非常に少数のものの選択（および除去）を行うことが好ましい。

パラメトリック生存モデルを適用するとき、生存時間は予備処理を必要とせず、

として表される。パラメトリック生存モデルは次のように適用される：
生存時間のパラメトリック密度関数を、

で表し、その生存関数を、

で表し、ここで、

は密度関数のパラメトリック形態に関連したパラメータであり、

は上記に規定されたものである。ハザード関数は以下のように規定される：

対数尤度関数の一般的な公式は、検閲削除されたデータを考慮して、以下のように表されることが好ましい：

パラメトリック回帰生存モデルを介した生存時間データの解析に関する標準的なテキストへの参照により、使用されることのできる生存時間分布の集りが示される。使用されることのできる生存分布には、たとえば、ワイブル、指数または極値分布が含まれる。

ハザード関数が

として書かれることができる場合、

ワイブル、指数または極値分布は、このようにして表された形態で書かれることのできる密度およびハザード関数を有している。

詳細な適用は部分的に、AitkenおよびClayton氏等の文献(1980)に記載されているアルゴリズムに依存するが、しかしながら、それによってユーザは任意のパラメトリック基礎ハザード関数を指定することが可能になる。

AitkinおよびClayton氏等の文献(1980)によると、パラメトリック生存モデルをモデル化する好ましい尤度関数は以下のように表される：

AitkinおよびClayton氏等の文献(1980)には、式（５Ｄ）の結果として、ｃ_iは平均μ_iを有するポアソン変量として処理されてもよく、また、式（１１Ｄ）中の最後の項は、

に依存しない（もっとも、それは

に依存する）ことが示されている。

の事後分布は以下のように表されることが好ましい：

１実施形態において、

は失われたデータのベクトルとして処理され、式（６Ｄ）を最大化して、

の事後推定値を生成するために反復的な処理手続きが使用されることができる。式（１Ｄ）のプライアは、極大である事後推定値が希薄になる、すなわち、多数のパラメータが冗長である場合に、

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

の多くの構成要素がゼロであるプライア期待値が存在するため、推定は、推定されたβ_iの大部分がゼロであり、残りのゼロでない推定値が生存時間の適切な説明を提供するように行われる。

マイクロアレイデータの状況においては、この動作の結果、その事象時間に対する適切な説明を提供する非常に少数の遺伝子のセットが識別されることとなる。

上述したように、事後分布を最大化する構成要素加重は、反復的な処理手続きを使用して決定されることができる。構成要素および構成要素加重の事後分布を最大化するための反復的な処理手続きは、たとえば、Dempster氏他による文献(1977)に記載されているようなＥＭアルゴリズムであることが好ましい。

１実施形態において、ＥＭアルゴリズムは以下のステップを含んでいる：
１．ｎ＝０，Ｓ₀＝｛１，２，・・・，ｐ｝を設定することによりアルゴリズムを初期化し、

を初期化する
２． ｉεＳ_nのときは、β_i ⁽ⁿ⁾ ＝β_i ^*，
それ以外のときは、β_i ⁽ⁿ⁾ ＝０
と定義し、Ｐ_nを、

ような０と１とのマトリックスとする。
３．構成要素加重の事後分布の期待値を計算することにより推定ステップを行う。これは以下の関数を使用して行われる：

４．最大化ステップを行う。これは、以下のようにＮｅｗｔｏｎ Ｒａｐｈｓｏｎ反復を使用して行われることができる：

を設定し、ｒ＝０，１，２，・・・に対して

であり、ここでα_rは、以下を保証するようにラインサーチアルゴリズムによって選択される：

５．

と定義する。ここでε₁は、たとえば、１０^-5等の小さい定数である。ｎ＝ｎ＋１を設定し、

を選択する。ここで、

また、κ_nは、０＜κ_n＜１であるような減衰定数である。

６．収斂をチェックする。

である場合、中止し、ここでε₂は適切に小さく、その他の場合、上記のステップ２に進む。

別の実施形態において、最大化ステップの中のステップ（４）は、

をその期待値：

で置換することにより推定されることができる。

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

モデルがコックスの比例ハザードモデルであるときのＥＭアルゴリズムの適用の説明を助けるために、“動的加重”およびこれらの加重に基づいたマトリックスを規定することが好ましい。加重は以下のように表される：

これらの加重に基づいたマトリックスは以下のように表される：

加重のマトリックスに関して、ｌの第１および第２の微分を以下のように書いてもよい：

ここで、Ｋ＝Ｗ^**−Δ（Ｗ^*）である。したがって、ＥＭアルゴリズムのステップ（２）の一部（式７Ｄ）（式１０Ｄも参照）として示された変換マトリックスＰ_nから、

となることに注意しなければならない。

モデルがコックスの比例ハザードモデルであるとき、ＥＭアルゴリズムのＥステップおよびＭステップは、以下のように行われることが好ましい：
１．ｎ＝０，Ｓ₀＝｛１，２，・・・，ｐ｝を設定する。ｖを、構成要素Ｖ_iを有するベクトルとする。ここで、たとえば０．００１のような、ある小さいεに対して、
ｃ_i＝１である場合、Ｖ_i＝１−ε，
ｃ_i＝０である場合、Ｖ_i＝ε
である。ｆはｌｏｇ（ｖ／ｔ）であると定義する。

ｐ≦Ｎである場合、

によって初期値：

を計算し、
ｐ＞Ｎである場合、

によって初期値：

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

２． ｉ∈Ｓ_nの場合、β_i ⁽ⁿ⁾ ＝β_i ^*，
それ以外の場合、β_i ⁽ⁿ⁾ ＝０
と定義する。
Ｐ_nを、

ような０と１のマトリックスとする。

３．以下を計算することによりＥステップを行う：

β＝Ｐ_nγおよびβ⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾を使用することにより、以下の式が得られる：

４．Ｍステップを行う。これは、以下のようにＮｅｗｔｏｎ Ｒａｐｈｓｏｎ反復により行われることができる。

を設定し、ｒ＝０，１，２，・・・に対して、

ここで、α_rは、

を保証するようにラインサーチアルゴリズムによって選択される。
ｐ≦Ｎに対しては、以下の式を使用する：

ｐ＞Ｎに対しては、以下の式を使用する：

γ*を、ある収斂基準が満足されたとき、たとえば、
｜｜γ_r−γ_r+1｜｜＜ε（たとえば、１０^-5）
であるときのγ_rの値とする。

５．

と定義する。ここで、ε₁は、たとえば、１０^-5等の小さい定数である。このステップにより、非常に小さい係数を有する変数が除去される。

６．収斂をチェックする。

である場合、中止し、ここでε₂は適切に小さく、その他の場合には、ｎ＝ｎ＋１を設定し、上記のステップ２に進んで、収斂が発生するまで処理手続きを繰返す。

別の実施形態において、ＥＭアルゴリズムは、モデルがパラメトリック生存モデルであるときの事後分布を最大化するために適用される。

パラメトリック生存モデルにＥＭアルゴリズムを適用したとき、式（５Ｄ）の結果として、ｃ_iは平均μ_iを有するポアソン変量として処理されてもよく、また、式（５Ｄ）中の最後の項はβに依存しない（もっとも、それはφに依存する）。

であり、したがって、ポアソンのような平均に対する対数・線形モデルで問題を表すことができることに注意されたい。対数尤度関数の反復的な最大化が行われ、

の初期推定値が与えられた場合に

の推定値が得られることが好ましい。それ故、

のこれらの推定値が与えられた場合には、

の更新された推定値が得られる。この処理手続きは、収斂が発生するまで続けられる。

上述した事後分を適用することにより、

であることが認識される。
結果的に、式（１１Ｄ）および（１２Ｄ）から、

パラメトリック生存モデルに関連した式（１２Ｄ）の別の形は以下のとおりである：

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

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

とし、ここで、０＜κ_n≦１に対して

βは前のＭステップから得られた値に固定されていることが好ましい。

パラメトリック生存モデルおよびマイクロアレイデータの状況においてパラメータ選択のためのＥＭアルゴリズムを提供することができる。このＥＭアルゴリズムは以下のようなものであることが好ましい：
１．

を設定する。ｖを、構成要素Ｖ_iを有するベクトルとする。ここで、たとえば０．００１のような、ある小さいεに対して、
ｃ_i＝１である場合、Ｖ_i＝１−ε，
ｃ_i＝０である場合、Ｖ_i＝ε
である。ｆはｌｏｇ（ｖ／Λ（ｙ，φ））であると定義する。
ｐ≦Ｎである場合、

によって初期値：

を計算し、
ｐ＞Ｎである場合、

によって初期値：

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

２． ｉ∈Ｓ_nの場合、β_i ⁽ⁿ⁾ ＝β_i ^*，
それ以外の場合、β_i ⁽ⁿ⁾ ＝０
と定義する。
Ｐ_nを、

ような０と１のマトリックスとする。

３．以下を計算することによりＥステップを行う：

β＝Ｐ_nγおよびβ⁽ⁿ⁾＝Ｐ_nγ⁽ⁿ⁾を使用することにより、以下の式が得られる：

４．Ｍステップを行う。これは、以下のようにＮｅｗｔｏｎ Ｒａｐｈｓｏｎ反復により行われることができる。

を設定し、ｒ＝０，１，２，・・・に対して、

ここで、α_rは、

を保証するようにラインサーチアルゴリズムによって選択される。

ｐ≦Ｎに対しては、以下の式を使用する：

ｐ＞Ｎに対しては、以下の式を使用する：

γ*を、ある収斂基準が満足されたとき、たとえば、
｜｜γ_r−γ_r+1｜｜＜ε（たとえば、１０^-5）
であるときのγ_rの値とする。

５．

と定義する。ここで、ε₁は、たとえば、１０^-5等の小さい定数である。ｎ＝ｎ＋１を設定し、

を選択し、ここで、

κ_nは０＜κ_n＜１となるような減衰定数である。

６．収斂をチェックする。

である場合、中止し、ここでε₂は適切に小さく、その他の場合には、ステップ２に進む。

別の実施形態において、生存時間はワイブル生存密度関数によって表される。ワイブルの場合、

は１次元であり、

であることが好ましい。

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

が解かれることが好ましい。コックスの比例ハザードモデルに対して適用されたステップに後続して、αが推定されることができ、また、生存時間がワイブル分布に従ったものである場合にその生存時間を適切に説明することのできる

からのパラメータの非常に少数の部分集合が選択されることができる。

本発明の特徴および利点は、以下の例の説明から明らかになるであろう。

例
例１：対数回帰モデルを使用する前立腺癌に対する２グループ分類
組織を前立腺および非前立腺グループに分類することのできる遺伝子の部分集合を識別するために、Luo氏他の文献(2001)で報告され、解析されたマイクロアレイデータセットは、二項ロジスティック回帰がモデルとして使用される本発明の方法を使用して解析された。このデータセットは６５００の人の遺伝子に関するマイクロアレイを含んでいる。研究対象には、前立腺癌を有することが知られている１６人の被験者と、良性の前立腺細胞過形成である９人の被験者が含まれている。しかしながら、表現を単に簡潔にするために、５０の遺伝子が解析のために選択された。この５０の全ての遺伝子の遺伝子発現率（行）および２５人の患者（列）は表４に示されている。

以下、この方法を適用した結果を示す。モデルはＧ＝２クラスを有しており、５０の遺伝子の全てをそのモデル中の前立腺変数（構成要素または基礎関数）として開始された。２１回の反復後（以下を参照）、アルゴリズムは、完全な分類を与える２つの遺伝子（表５の番号３６および４７）を発見した。この結果がデータセット中で利用可能な多数の遺伝子（変数）による人為結果だったのか否かを決定するために、クラスラベルがランダムに置換され、後続的にアルゴリズムが適用される置換テストがランされた。これは２００回繰返された。図１は、正しく分類された場合の数のヒストグラムを示している。実際のデータセットに対する１００％の正確さはこの置換分布の末尾にあり、０．０１５のｐ値を有する。これは、この結果が偶然によるものではないことを示唆している。

以下に交換されていないデータに対する反復の詳細を示す：

例２：大きいデータセットおよび二項ロジスティック回帰モデルを使用する２グループ分類
組織を異なった臨床タイプのリンパ腫に分類することのできる遺伝子の部分集合を識別するために、Alizadeh,A.A.氏他の文献［Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling,Nature 403:503-511(2000)］で報告され、解析されたデータセットは、二項ロジスティック回帰がモデルとして使用される本発明の方法を使用して解析された。

データセットには、ｎ＝４０２６の遺伝子およびｎ＝４２のサンプルが存在している。以下において、ＤＬＢＣＬとは“びまん性大Ｂ細胞リンパ腫”のことである。サンプルは、ＧＣ Ｂ様ＤＬＢＣＬ（２１個のサンプル）および活性化されたＢ様ＤＬＢＣＬ（２１個のサンプル）の２つの疾病タイプに分類されている。このセットを使用して、異なった疾病タイプの診断に役立つ遺伝子を迅速に発見するための上記の方法の使用法を説明する。

以下において、この方法を適用した結果が与えられている。モデルはＧ＝２クラスを有しており、全ての遺伝子をそのモデル中の潜在的変数（基礎関数）として開始された。２０回の反復後、アルゴリズムは、以下の誤分類（表５）を与える番号１２８１および１３１２（GENE3332XおよびGENE3258X）の２つの遺伝子を発見し、全体的な分類成功率は９８％であった。この例は、ラップトップマシンで約２０秒で稼動した。

表 ５

予測されたクラス１ 予測されたクラス２
真のクラス１ ２０ １
真のクラス２ ０ ２１

結果がデータセット中で利用可能な多数の遺伝子（変数）による人為結果だったのか否かを決定するために、クラスラベルがランダムに置換され、後続的にアルゴリズムが適用される置換テストがランされた。これは１０００回繰返された。図２は、正しく分類された場合の割合（ラムダ）のヒストグラムを示している。実際のデータセットに対する９７．６％の正確さはこの置換分布の末尾にあり、０．０１３のｐ値を有する。これらの観察は、この結果が偶然によるものではないことを示唆している。

例３：多グループ分類
サンプルを多数のクラスの１つに分類することのできる遺伝子を識別するために、Yeoh氏他の文献［Cancer Cell v1 : 133-143(2002)］で報告され、解析されたデータセットは、尤度が多項ロジスティック回帰に基づいて使用される本発明の方法を使用して解析された。Yeoh氏他の文献に記載されているものと同じ予備処理が適用されている。これは、以下から構成されている：
・８つのアレイ：ＢＣＲ．ＡＢＬ．Ｒ４，ＭＬＬ．Ｒ５，Ｎｏｒｍａｌ．Ｒ４，Ｔ．ＡＬＬ．Ｒ７，Ｔ．ＡＬＬ．Ｒ８，Ｈｙｐｅｒｄｉｐ．５０．２Ｍ．３，Ｈｙｐｏｄｉｐ．２Ｍ．３，およびＨｙｐｏｄｉｐ．２Ｍ．２をドロップする
・各アレイの平均反応値を２５００に設定する
・しきい値化する；１００より小さい４５０００個の値に設定された４５０００個を越える値が１に設定される
・存在が０．０１未満の遺伝子が除去される；これは合計１６０７個の遺伝子である
・最大値と最小値との間の差が１００未満であった遺伝子（１６０４個）が除去される

予備処理後、ｎ＝１１００５個の遺伝子およびｎ＝２４８個のサンプルが存在している。サンプルは６つの疾病タイプに分類されている：
１．ＢＣＲ−ＡＢＬ；
２．Ｅ２Ａ−ＰＢＸ１；
３．Ｈｙｐｅｒｄｉｐ＞５０；
４．ＭＬＬ；
５．Ｔ−ＡＬＬおよび
６．ＴＥＬ−ＡＭＬ１

このセットは、異なった疾病タイプの診断に役立つ遺伝子を迅速に発見するための方法の使用法を示すために使用された。この方法を適用した結果を以下に示す。モデルはＧ＝６のクラスを有しており、全ての遺伝子をそのモデル中の潜在的変数（基礎関数）として開始された。２０回の反復の後、アルゴリズムは、以下の１０個の遺伝子がクラスを分けることを発見した：

１５重の相互確認（ｃｒｏｓｓ ｖａｌｉｄａｔｉｏｎ）から、分類成功率が９４％である以下の誤分類表（表６）が得られた：

交換テスト（クラスラベルを交換する）は、０．９４％の相互確認されたエラーレートが非常に重要であること（ｐ＝０．００）を示した。

例４：一般化された線形モデルを使用する標準回帰
この例は、この方法が一般化された線形モデルフレームワークでどのようにして実施されることができるかを示す。この例は、２００の観察および４１の変数（基礎関数）に関する標準回帰問題である。真の曲線はエラー（または雑音）を伴なって観察され、変数の一部だけに依存することが知られている。反応は連続し、正規分布される。これらのデータは、一般化された線形モデル変数選択のために本発明のアルゴリズムを使用して解析される。

これは、以下を有する一般化された線形モデルである：
リンク関数：ｇ（μ）＝μ
リンク関数の偏微分：

分散関数：τ²＝１
スケールパラメータφ＝σ²
逸脱（尤度関数）：

σ²に対する更新公式は

であり、ここでμ_i ^*はアルゴリズムのステップ５においてβ^*で推定された平均である。

以下にアルゴリズムの出力を示す。
ＥＭ反復：１ 期待事後分布：-55.45434 基礎関数 ４１ シグマの２乗 0.5607509
ＥＭ反復：２ 期待事後分布：-43.96193 基礎関数 ４１ シグマの２乗 0.5773566
ＥＭ反復：３ 期待事後分布：-48.87198 基礎関数 ３９ シグマの２乗 0.5943395
ＥＭ反復：４ 期待事後分布：-52.79632 基礎関数 ３１ シグマの２乗 0.6072137
ＥＭ反復：５ 期待事後分布：-55.18578 基礎関数 ２８ シグマの２乗 0.6161707
ＥＭ反復：６ 期待事後分布：-56.5303 基礎関数 ２３ シグマの２乗 0.6224545
ＥＭ反復：７ 期待事後分布：-57.47589 基礎関数 １７ シグマの２乗 0.626674
ＥＭ反復：８ 期待事後分布：-58.0566 基礎関数 １５ シグマの２乗 0.6293923
ＥＭ反復：９ 期待事後分布：-58.41912 基礎関数 １３ シグマの２乗 0.6315789
ＥＭ反復：１０ 期待事後分布：-58.6923 基礎関数 １１ シグマの２乗 0.633089
ＥＭ反復：１１ 期待事後分布：-58.88766 基礎関数 １０ シグマの２乗 0.6343793
ＥＭ反復：１２ 期待事後分布：-59.05261 基礎関数 １０ シグマの２乗 0.635997
ＥＭ反復：１３ 期待事後分布：-59.24126 基礎関数 ９ シグマの２乗 0.6381456
ＥＭ反復：１４ 期待事後分布：-59.47668 基礎関数 ９ シグマの２乗 0.640962
ＥＭ反復：１５ 期待事後分布：-59.7677 基礎関数 ９ シグマの２乗 0.6443392
ＥＭ反復：１６ 期待事後分布：-60.10277 基礎関数 ９ シグマの２乗 0.6477088
ＥＭ反復：１７ 期待事後分布：-60.44193 基礎関数 ９ シグマの２乗 0.6508144
ＥＭ反復：１８ 期待事後分布：-60.7684 基礎関数 ９ シグマの２乗 0.6539145
ＥＭ反復：１９ 期待事後分布：-61.09251 基礎関数 ９ シグマの２乗 0.6565873
ＥＭ反復：２０ 期待事後分布：-61.38427 基礎関数 ８ シグマの２乗 0.6589498
ＥＭ反復：２１ 期待事後分布：-61.65061 基礎関数 ８ シグマの２乗 0.6615976
ＥＭ反復：２２ 期待事後分布：-61.92217 基礎関数 ８ シグマの２乗 0.664281
ＥＭ反復：２３ 期待事後分布：-62.17683 基礎関数 ７ シグマの２乗 0.6663748
ＥＭ反復：２４ 期待事後分布：-62.37402 基礎関数 ７ シグマの２乗 0.6679655
ＥＭ反復：２５ 期待事後分布：-62.51645 基礎関数 ７ シグマの２乗 0.6689011
ＥＭ反復：２６ 期待事後分布：-62.59567 基礎関数 ６ シグマの２乗 0.6689011
ＥＭ反復：２７ 期待事後分布：-62.6151 基礎関数 ６ シグマの２乗 0.6690962
ＥＭ反復：２８ 期待事後分布：-62.61717 基礎関数 ６ シグマの２乗 0.6691031
ＥＭ反復：２９ 期待事後分布：-62.61739 基礎関数 ５ シグマの２乗 0.6691035

アルゴリズムは、４１の基礎ベクトル（変数）の中の５つを含むモデルにより収斂する。図３には、アルゴリズムにより選択された５つの変数（基礎関数）を有するモデルに対するあてはめられた曲線（実線）、真の曲線（破線）、および雑音の多いデータのプロットが示されており、ここでｙ変数はｎｆで示されている。

例５：一般化された線形モデルを使用する小さい線形回帰の例
この例は例４に類似しているが、しかし簡潔にするために、少数の変数（１０）が使用される。これによって、全データセットが一覧表にされることが可能になる（表７参照）。従属変数は最初の４つの変数だけの関数であり、残りの変数は雑音である。

データは、アイデンティティリンク、一定の分散および正常な反応を有する一般化された線形モデルとして解析された。１２回の反復後、アルゴリズムは、予測情報を有することが知られている４つの変数を含み、６つの雑音変数を全て廃棄した解に収斂した。

例６：リンパ種生存解析
この例は、http://llmpp.nih.gov/lymphoma/data.shtml からの実際の生存データを使用する。

参考文献は、Alizadeh AA氏他の文献(2000)［Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature 403(6769):503-11］である。

データは４０２６個の遺伝子および４０個のサンプル（個人）に対するデータからなるマイクロアレイデータであり、生存時間および検閲削除インジケータが各サンプルに対して利用可能である。結果は、アルゴリズムを使用し、コックスの比例ハザードモデルを実施して解析された。

アルゴリズムは３つの遺伝子（遺伝子：３７９７Ｘ，３３０２Ｘ，３５６Ｘ）を生存時間と関連付けられたものとして選択したことに注意されたい。

例７：減少されたリンパ腫生存解析
文書による証明を完全なものにするために、ランダムに選択された４７個の遺伝子と全データセットの解析で意味のあるものとして識別された３個の遺伝子とを含む５０個の遺伝子が、Alizadeh氏他の文献に記載されている遺伝子の部分集合に基づいて選択された例もまた示されている。データは以下の表９に示されている。

データは、コックスの比例ハザード生存モデルを含むアルゴリズムの別形式を使用して解析された。２２回の反復後、全セットに対する解からの２つの遺伝子を含む５つの遺伝子が選択された。全ての結果（反復履歴を含む）を以下に示す：

例８：パラメトリックハザードによる生存解析
データは、

からの1694w.datである。これは、黒色腫の生存に関するデータである。ｎ＝２５５の個人が存在し、その中の１００人は生存時間を検閲削除されている。各個人は４つのコバリエイト、すなわち、治療、シックネス（ｔｈｉｃｋｎｅｓｓ）、年齢および性別を有している。方法を示すために、４０００個のダミーの遺伝子がこのデータセットに追加され、その結果４００４列と２５５行を有するデータマトリックスが得られた。故意に、４０００個の“遺伝子”は生存時間と関連付けられていない。そのほとんどが情報を持たない４００４個の可能性のある予測子から重要な変数を識別することは、アルゴリズム的に難題である。ハザード関数に対してパラメトリックワイブルモデルを使用して、データが解析された。

アルゴリズムは変数：年齢に関してのみ選択を行った。擬似遺伝子変数は全て迅速に廃棄された。ワイブル形状パラメータは０．６８と推定された。

例９：前立腺癌に対する順序付けられたカテゴリー的解析
この例は、Dhanasekaran氏他による文献(2001)からのものである。

および

における補足ファイルもまた参照されたい。

９６０５個の遺伝子と共に１５個のサンプル（個人）が存在する。失われた値は、行平均＋列平均マイナス総（ｇｒａｎｄ）平均によって置換された。４つの順序付けられたカテゴリー（Ｇ＝４）、すなわち、
１．ＮＡＰ正常な
２．ＢＰＨ良性の
３．ＰＣＡ局所化した
４．ＭＥＴ転移した
カテゴリーが存在する。

アルゴリズムは、１つの誤分類を除いて全ての個人を正しく分類することのできる１つの遺伝子（遺伝子番号６６１１、それらの取得ＩＤ Ｒ３１６７９）を発見した。

ＥＭアルゴリズムからの反復は以下のとおりである：

あてはめられた確率のプロットは以下の図６に示されている。ラインは以下のようにクラスを表している：破線＝クラス１、実線＝クラス２、点線＝クラス３、点線と破線を組合せた線＝クラス４。観察（インデックス）１乃至３はクラス２中にあり、観察４乃至７はクラス１中にあり、観察８乃至１１はクラス３中にあり、観察１２乃至１５はクラス４中にある。

例１０：前立腺癌−選択された遺伝子に対する順序付けられたカテゴリー的解析
この例は、データセットが５０個の選択された遺伝子に減少されていることを除いて、例９のものと同一である。これらの遺伝子の１つは例９において重要であると認められた遺伝子であり、その他のものはランダムに選択された。この例の目的は、完全に表にされたデータセット（表１０）に基づいて例証を行うことである。

失われた値は、行平均＋列平均マイナス総（ｇｒａｎｄ）平均によって置換された。４つの順序付けられたカテゴリー（Ｇ＝４）、すなわち、
１．ＮＡＰ正常な
２．ＢＰＨ良性の
３．ＰＣＡ局所化した
４．ＭＥＴ転移した
カテゴリーが存在する。

アルゴリズムは、例９の遺伝子６６１１（取得 Ｒ３１６７９）に等価である１つの予測遺伝子（表１０の遺伝子１）を発見した。予測の成功は、当然ながら、例９のものと同じであった（それは同じ単一の遺伝子に基づいていたので）。

例１１：方法を使用するための装置
図５を参照とすると、本発明の実施形態による方法を実施するのに適したパーソナルコンピュータ20が示されている。コンピュータ20は、ハードディスクデータ記憶装置21上に記憶されたソフトウェアタプログラムの命令の下で動作する。コンピュータ20はさらに、プロセッサ22、メモリ23、表示スクリーン24、プリンタ25、ならびにマウス26およびキーボード27のような入力装置を備えている。コンピュータはインターネット28へのネットワーク接続27のような通信手段、またはデータのダウンロードまたは収集および共用を容易にするデータ収集手段28を有していてもよい。

データ収集手段はデータをシステムから収集し、あるいはダウンロードする。コンピュータはソフトウェアに埋込まれた操作手段を備えており、この操作手段はマウス26およびキーボード27と通信して、ユーザが本発明の実施形態による方法をデータに関して実施することを可能にする。システムは、本発明の実施形態による方法を実施するソフトウェアに埋込まれた手段と、およびグラフィックを生成する手段とを備えている。この方法が実施された後、出力は表示スクリーン24上にグラフィック的に示され、および、またはプリンタ25でプリントされることができる。

上記の例では、生物学的システムに関して本発明の構成を説明してきた。前に述べたように、本発明は、サンプルの特徴が予測されることを必要とする任意の“システム”に適用されることができる。システムの例には、化学的システム、農業システム、気象システム、たとえば、クレジットリスク査定システムを含む金融システム、保険システム、マーケティングシステム、または企業記録システム、電子システム、物理的システム、天文物理学的システムおよび機械的システムが含まれる。

当業者に明らかであるように、修正および変形は本発明の技術的範囲内であると考えられる。

参照文献

本発明の１実施形態の予測の成功に関する置換テストの結果を示すヒストグラム。 本発明の１実施形態の予測の成功に関する置換テストの結果を示すヒストグラム。 本発明の方法の１実施形態において使用される一般化された線形モデルに対する曲線のプロット。 本発明の方法の１実施形態を使用して識別された単一の遺伝子に対するあてはめられた確率のプロット。 本発明によるシステムを実施するために使用されるパーソナルコンピュータの概略図。

Claims (35)

1. テストサンプルの特徴を予測することのできるシステムの構成要素の部分集合を識別する方法において、
   （ａ）構成要素と構成要素加重との線形結合を発生し、この線形結合において、各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が導かれ、
   （ｂ）ある特徴の確率分布に対するモデルを規定し、このモデルは線形結合を条件付けし、また、このモデルは２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではなく、
   （ｃ）ゼロに近い高い確率密度を有するハイパープライア（ｈｙｐｅｒｐｒｉｏｒ）を含む線形結合の構成要素加重に対して事前（ｐｒｉｏｒ）分布を構成し、
   （ｄ）この事前分布とモデルとを組合せて事後（ｐｏｓｔｅｒｉｏｒ）分布を生成し、
   （ｅ）この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別するステップを含んでいる方法。
2. モデルは、多項または二項ロジスティック回帰、一般化された線形モデル、コックスの比例ハザードモデル、およびパラメトリック生存モデルを含むグループから選択されたモデルに基づいた尤度関数である請求項１記載の方法。
3. モデルは、多項または二項ロジスティック回帰に基づいた尤度関数である項１または２記載の方法。
4. ロジスティック回帰は、多項または二項分布を有する特徴をモデル化する請求項２または３記載の方法。
5. 構成要素の部分集合はロジスティック回帰を規定することによってサンプルを複数の予め規定されたグループの１つに分類することが可能であり、このロジスティック回帰はサンプルを複数のサンプルグループにグループ化することを含んでおり、各サンプルグループは共通のグループ識別子を有している請求項１乃至４のいずれか１項記載の方法。
6. ロジスティック回帰は、
   

   

   の形態のものであり、ここで、
   ｘ_i ^Tβ_gは、構成要素加重β_gを有する訓練サンプルｉからの入力データから発生された線形結合であり、
   ｘ_i ^Tは、Ｘのｉ番目の行に対する構成要素であり、β_gはサンプルクラスｇに対する構成要素加重のセットであり、
   訓練サンプルｉがクラスｇのメンバーである場合、ｅ_ig＝１であり、そうでない場合にはｅ_ig＝０であり、
   Ｘは、ｐ個の構成要素を含むｎ個の訓練サンプルからのデータである請求項１乃至５のいずれか１項記載の方法。
7. 構成要素の部分集合は、ロジスティック回帰を規定することによってサンプルをクラスに分類することができ、このクラスは複数の予め規定された順序付けられたクラスの１つであり、前記ロジスティック回帰は、各グループ識別子が順序付けられたクラスのメンバーに対応している一連のグループ識別子を規定し、サンプルを順序付けられたクラスの１つにグループ化することを含んでいる請求項１または２記載の方法。
8. ロジスティック回帰は、
   

   

   の形態のものであり、ここで、
   

   

   は、訓練サンプルｉがｋ以下である識別子を有するクラスに属している確率であり（ここで、順序付けられたクラスの合計はＧである）、
   

   

   は、構成要素加重：
   

   

   を有する訓練サンプルｉからの入力データから発生された線形結合であり、
   Ｘは、ｐ個の構成要素を含むｎ個の訓練サンプルからのデータであり、
   ｘ_i ^Tは、Ｘのｉ番目の行に対する構成要素であり、
   ｒ_ijは、
   

   

   として定義されており、ここで、
   観察ｉがクラスｊ中のものである場合はｃ_ij＝１であり、そうでない場合はｃ_ij＝０である請求項７記載の方法。
9. モデルは、一般化された線形モデルに基づいた尤度関数である請求項１または２記載の方法。
10. 一般化された線形モデルは、分布の正規指数族として分布された特徴をモデル化する請求項９記載の方法。
11. 分布の正規指数族は、正規分布、ガウス分布、ポアソン分布、指数分布、ガンマ分布、ｘ²分布、および逆ガンマ分布からなるグループから選択される請求項１０記載の方法。
12. 構成要素の部分集合は、一般化された線形モデルを規定することによってサンプルの予め規定された特性を予測することができ、この一般化された線形モデルは予測される特性をモデル化することを含んでいる請求項１または２記載の方法。
13. 一般化された線形モデルは、
    

    

    の形態のものであり、ここで、
    ｙ＝（ｙ₁，・・・，ｙ_n）^Tであり、ｙ_iはｉ番目のサンプルに関して測定された特性であり、
    ａ_i（φ）＝φ／ｗ_iであり、ｗ_iは既知の加重の固定されたセットであり、φは単一の
    スケールパラメータであり、
    関数ｂ（．）およびｃ（．）はNelderおよびWedderburn氏等による文献（1972）において定義されたものであり、
    Ｅ｛ｙ_i｝＝ｂ’（θ_i）
    Ｖａｒ｛ｙ｝＝ｂ’’（θ_i）ａ_i（φ）＝τ_i ²ａ_i（φ）であり、
    ここで、各観察はコバリエイト（ｃｏｖａｒｉａｔｅ）ｘ_iと線形予測子η_i＝ｘ_i ^Tβとのセットを有している請求項９または１０記載の方法。
14. モデルは、コックスの比例ハザードモデル、パラメトリック生存モデルおよび加速された生存時間モデルからなるグループから選択されたモデルに基づいた尤度関数である請求項１または２記載の方法。
15. 構成要素の部分集合は、コックスの比例ハザードモデル、パラメトリック生存モデルまたは加速された生存時間モデルに基づいて尤度を規定することにより、あるサンプルに対するある事象までの時間を予測することができ、この尤度は、複数のサンプルに対して、そのサンプルが得られた時間からその事象の時間までに経過した時間を測定することを含んでいる請求項１記載の方法。
16. コックスの比例ハザードモデルは、
    

    

    の形態のものであり、ここで、
    Ｘは、ｐ個の構成要素を含むｎ個の訓練サンプルからのデータであり、
    Ｚは、Ｘの行の再配列であるマトリックスであり、ここで、Ｚの行の順序付けは生存時間の順序付けにより誘発された順序付けに対応し、
    ｄは、生存時間を順序付けるために必要とされる同じ置換により検閲削除インデックスを順序付けた結果であり、
    Ｚ_jはマトリックスＺのｊ番目の行であり、ｄ_jはｄのｊ番目のエレメントであり、
    

    

    Ｒ_j＝｛ｉ：ｉ＝ｊ，ｊ＋１，・・・，Ｎ｝＝ｊ番目の順序付けられた事象時間ｔ(_j)におけるリスクセットである請求項１４記載の方法。
17. パラメトリックハザードモデルは、
    

    

    の形態のものであり、ここで、
    

    

    ｉ番目のサンプルが検閲削除されていない場合はｃ_i＝１、ｉ番目のサンプルが検閲削除されていない場合はｃ_i＝０であり、
    関数λ（．）およびΛ（．）は、AitkinおよびClayton氏等による文献（1980）において定義されたものであり、
    ｘ_jはＸのｉ番目の行であり、Ｘはｐ個の構成要素を含むｎ個の訓練サンプルからのデータである請求項１４記載の方法。
18. 事前分布は、
    

    

    の形態のものであり、
    ここで、ｐ（β|ｖ²）はＮ（０，ｄｉａｇ｛ｖ²｝）であり、
    ｖはハイパーパラメータであり、
    ｐ（ｖ²）はハイパープライア分布である請求項１乃至１７のいずれか１項記載の方法。
19. ハイパープライアは、
    

    

    の形態のジェフェリーのプライア（ｐｒｉｏｒ）である請求項１乃至１８のいずれか１項記載の方法。
20. 事後分布は、
    

    

    の形態のものであり、
    ここで、
    

    

    は尤度関数である請求項１乃至１９のいずれか１項記載の方法。
21. 事後分布は、反復的な処理手続きを使用して最大化される請求項１乃至２０のいずれか１項記載の方法。
22. 反復的な処理手続きは、ＥＭアルゴリズムである請求項２１記載の方法。
23. システムは、生物学的システムである請求項１乃至２２のいずれか１項記載の方法。
24. 生物学的システムは、バイオテクノロジーアレイである請求項２３記載の方法。
25. バイオテクノロジーアレイは、ＤＮＡアレイ、蛋白質アレイ、抗体アレイ、ＲＮＡアレイ、炭水化物アレイ、化学物質アレイ、脂質アレイからなるグループから選択される請求項２４記載の方法。
26. 各グループがテスト治療に対する反応によって規定されている複数の予め規定されたグループの１つに部分集合を分類することのできる被験者の構成要素の部分集合を識別する方法において、
    （ｄ）複数の被験者をテスト治療に露出し、そのテスト治療に対する反応に基づいて被験者を反応グループにグループ化し、
    （ｅ）被験者の構成要素を測定し、
    （ｆ）請求項１乃至２５のいずれか１項記載の方法を使用して被験者を反応グループに分類することのできる構成要素の部分集合を識別するステップを含んでいる方法。
27. 構成要素は、遺伝子、小さいヌクレオチド多型（ＳＮＰｓ）、蛋白質、抗体、炭水化物、脂質からなるグループから選択される請求項２６記載の方法。
28. システムから発生されたデータからのテストサンプルの特徴を、そのシステムからの複数のサンプルから予測することのできるシステムの構成要素の部分集合を識別する装置において、
    （ａ）各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が導かれる、構成要素と構成要素加重との線形結合を発生する手段と、
    （ｂ）線形結合を条件付けし、また、２クラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではない、ある特徴の確率分布に対するモデルを規定する手段と、
    （ｃ）ゼロに近い高い確率密度を有するハイパープライアを含む線形結合の構成要素加重に対して事前分布を構成する手段と、
    （ｄ）この事前分布とモデルとを組合せて事後分布を生成する手段と、
    （ｅ）この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別する手段とを具備している装置。
29. 計算装置にロードされたときにこの計算装置を制御して、請求項１乃至２７のいずれか１項記載の方法を実行するように構成されているコンピュータプログラム。
30. 請求項１乃至２７のいずれか１項記載の方法により実施される請求項２９記載のコンピュータプログラム。
31. 請求項２９または３０記載のコンピュータプログラムを提供するコンピュータ読出し可能な媒体。
32. システムからのサンプルをテストしてそのサンプルの特徴を識別する方法において、
    特徴の診断に役立つ構成要素の部分集合をテストするステップを含んでおり、
    構成要素の部分集合は請求項１乃至２７のいずれか１項記載の方法によって決定されている方法。
33. サンプルの特徴を決定するためにシステムからのサンプルをテストする装置において、
    請求項１乃至２７のいずれか１項記載の方法にしたがって識別された構成要素をテストする手段を備えている装置。
34. 計算装置上で実行されたときに、システムからのテストサンプルの特徴を予測することのできるそのシステムからの構成要素を識別する方法でこの計算装置を制御するように構成されており、構成要素と構成要素加重との線形結合は、各訓練サンプルが１つの既知の特徴を有する複数の訓練サンプルから発生されたデータから発生され、ゼロに近い高い確率密度を有するハイパープライアを含む構成要素加重に対する事前分布と、線形結合に条件付けするモデルとを組合せることによって事後分布が発生され、それによって事後分布を最大にする構成要素加重を推定し、前記モデルは２つのクラス反応に対する二項分布を、線形結合とこの反応の期待値とをリンクするプロビット関数と組合せたものではないコンピュータプログラム。
35. 生物学的システムからのテストサンプルの特徴を予測することのできる生物学的システムの構成要素の部分集合を識別する方法において、
    （ａ）構成要素と構成要素加重との線形結合を発生し、この線形結合において、各訓練サンプルが１つの既知の特徴を有している複数の訓練サンプルから発生されたデータから各構成要素に対する値が決定され、
    （ｂ）ある特徴の確率分布に対するモデルを規定し、このモデルは線形結合に条件付けし、
    （ｃ）ゼロに近い高い確率密度を有するハイパープライアを含む線形結合の構成要素加重に対して事前分布を構成し、
    （ｄ）この事前分布とモデルとを組合せて事後分布を生成し、
    （ｅ）この事後分布を最大にする構成要素加重を有する構成要素の部分集合を識別するステップを含んでいる方法。

Images