図3は、一実施形態に係る予測モデル構築装置の機能ブロック図である。予測モデル構築装置10は、文書化部1、年代別分類部2及び粒度考慮モデル構築部40を備える。ここで、粒度考慮モデル構築部40はさらにクラスタリング部3、モデル構築部4及び分割幅設定部5を備える。当該各部の処理は以下の通りである。
文書化部1には、予測モデル構築装置10によって予測モデルを構築するための入力としての全データとして、医療データが入力され、当該医療データにおける各対象者の各年代における健康評価内容を文書化して、年代別分類部2に渡す。本発明では、当該入力される全データは、図1の[2]で説明したように現実に入手しうるデータであって、欠損があることが想定されている。
医療データは対象者の健康状態を評価したものであり、例えば健康診断結果のデータ(健診データ)や医師による問診その他診療の結果のデータ(レセプトデータ)等で構成することができ、文書化部1では各対象者の各年代における当該医療データを、当該対象者の当該年代における健康状態を表す特徴ベクトルとして文書化する。具体的には、周知の形式であるバグオブワーズ(Bag of Words)の形式で文書化する。すなわち、健康状態を表す複数m個の単語i1, i2, ,…, imを用意しておき、それらをそれぞれ何回用いるかという単語頻度のm次元ベクトルの形で、文書化する。当該文書化の処理は、後述のクラスタリング部3におけるクラスタリングを可能な形式とするためのものである。
例えば、対象者Xのn歳における健康状態を表す特徴ベクトルをV(X, n)と表記することにすると、図1の[2]に示されている対象者Aの40歳〜43歳の医療データに対しては文書化部1により、V(A, 40), V(A, 41), V(A, 42)及びV(A, 43)が出力されることとなる。従って、文書化部1によりこのようにして、医療データ内に含まれる全ての(X, n)の組について、V(X, n)が得られ、年代別分類部2に渡されることとなる。
なお、本発明においては、年代は上記のように1年ごとに区切られているものとして説明するが、当該区切りの長さは実際に利用する医療データに応じて任意の長さとすることができる。例えば、2年ごととしてもよいし、半年ごととしてもよい。
なお、当該特徴ベクトルV(X, n)を構成する単語頻度における各単語i1, i2, ,…, imは、医療データ内の評価項目(例えば健診データであれば、身長、体重、BMI、血液検査の各項目等)毎に所定の単語を用意しておき、当該評価項目についての対象者Xのn歳時点における評価内容(一般に数値やテキストで与えられている)に応じた頻度を与えることで文書化すればよい。
文書化部1における当該頻度を与えることによる文書化には、評価項目ごとにその評価内容を頻度へと変換する所定の変換式を用意しておき、当該変換式を適用してもよいし、本出願人による特願2013-159323号(数値データ解析装置及びプログラム)、特願2013-163207号(数値データ解析装置及びプログラム)、特願2013-217817号(数値データ文書化装置及びプログラム)などを利用してもよい。
なお、文書化部1による文書化による単語頻度(バグオブワーズ)形式への変換は、前述のようにクラスタリング部3におけるクラスタリングが可能な形式への変換である。後述するように、クラスタリング部3における大きな粒度でのクラスタリングの事前処理として大きな粒度での特徴ベクトルを求める際の処理や、分割幅設定部5の処理等において、単語頻度に変換する前の直接の数値(例えば体重のそのままの値)を利用する場合もある。このような場合の処理のために、文書化部1における単語頻度への変換後も、単語頻度に対して変換前の直接の数値の情報をいわば補助的な情報として紐付けておくものとする。
年代別分類部2は、文書化部1より得た全対象者Xの全n歳における医療データの特徴ベクトルV(X, n)を、年代nごとに分類して、当該年代nごとに分類された医療データD(n)をクラスタリング部3に渡す。
図4は、仮に入力された全医療データが図1の[2]に示すもの(全対象者はA〜F)であったとした場合の、年代別分類部2による出力を模式的に示す図である。図4では、当該年代nに分類された各D(n)を構成している対象者がD(n)の領域内部に表記されている。例えば42歳の全データD(42)を構成しているのは、対象者A, B, Cの当該42歳時点における各データV(A, 42), V(B, 42), V(C, 42)となる。
なお、以上の文書化部1及び年代別分類部2は、図3に示す順とは逆に適用されてもよい。すなわち、予測モデル構築装置10への入力としての文書化される前の医療データを年代別分類部2で年代ごとの各患者の医療データに分類してから、文書化部1で文書化してもよい。いずれの順であっても、次の粒度考慮モデル構築部40におけるクラスタリング部3には同様のデータD(n)が渡される。
粒度考慮モデル構築部40では、年代n別に分類されたデータD(n)を入力として、粒度を考慮した健康の予測モデルを出力する。具体的には、長期的な予測をする際の健康予測モデル(粒度の大きな予測モデル)と、短期的な予測をする際の健康予測モデル(粒度の小さな予測モデル)と、をそれぞれ個別に生成する。
図5は、粒度を考慮した健康の予測モデルの模式的な例を、n=20歳〜28歳の範囲において示す図である。上段側の[1]には粒度の小さな予測モデルの例として、1年単位で健康状態のクラスタを構成するモデルが示されている。なお、当該1年単位の予測モデルは前述の特許文献4の手法により生成可能である。一方、下段側の[2]には粒度の大きな予測モデルの例として、3年単位で健康状態のクラスタを構成するモデルが示されている。
すなわち、n=20歳〜28歳の9年に渡る範囲を、[1]のモデルでは1年毎に分割して時間軸方向に9回分の状態に分けているのに対し、[2]のモデルでは3年毎に分割して時間軸方向に3回分の状態に分けている。当該1年や3年といったような時間軸方向の分割単位が粒度であり、1年と3年とを比較した場合、1年は小さな粒度に、3年は大きな粒度に、相対的な関係において相当することとなる。
図5では3年の粒度のモデルである[2]において例えばn=20〜22歳の範囲でのクラスタがC(20〜22, 1), C(20〜22, 2), C(20〜22, 3)となっており、当該年齢範囲における3つの異なる健康状態にそれぞれ対応している。なお、当該図5の[2]における表記すなわち、n歳〜n+k歳の範囲でのi番目のクラスタをC(n〜n+k, i)と表記することは以降継続して用いることとする。
粒度考慮モデル構築部40では以上、図5にて模式的に説明したように粒度を考慮した健康の予測モデルを生成する。当該生成するため、具体的に各部3,4,5は以下のような処理を行う。
まず、分割幅設定部5は生成する予測モデルの粒度、すなわち、予測モデルの時間単位となる年齢を区切る分割幅を設定する。当該設定は予測モデル構築装置10を利用するユーザ(データ分析担当者等)の設定に従う所定設定としてもよいし、入力される各年代nからなる医療データD(n)全体を解析することで、当該医療データに適した設定としてもよい。
以下、当該設定はユーザ設定に従って小さい粒度については1年単位、大きい粒度についてはN年単位(N≧2)で設定されたものとして説明する。前述の図5の[2]はN=3の例となっている。なお、その他の設定も可能であるが、その他の設定の場合については後述する。
当該設定された分割幅に従ってそれぞれ、クラスタリング部3はクラスタリングを行い、医療データ内の対象者をクラスタに分類する。クラスタリングはまず、小さい粒度である1年単位において実施し、その結果をもとに大きい粒度であるN年単位において実施することができる。具体的にはそれぞれ、以下(1)、(2)の通りである。
(1)1年単位の小さい粒度でのクラスタリングにおいては、n歳の医療データD(n)に含まれる対象者Xがその特徴ベクトルV(X, n)で表現されるものとして、対象者Xをクラスタリングして各クラスタC(n, i)(i=1, 2, …, K(n))を結果として得る。ここでK(n)は当該n歳の医療データD(n)をクラスタリングした際のクラスタ数である。
クラスタリング手法としては、前述のLDA等の潜在トピック分析によるものを用いてもよいし、LDAを拡張した前述のTTMを用いてもよい。なお、TTMを用いる際はLDAを用いる場合とは異なり、クラスタリング可能となる対象が、1年以上の一連の期間において共通してデータが存在する対象者に限定され、これによって時系列影響を考慮して各年齢nにおけるクラスタリングが可能となる。TTMを適用する1年以上の期間は例えば大きい粒度として設定されているN年を用いることができる。
なお、クラスタ数K(n)については、所定数を設定しておいてもよいし、AIC(赤池情報量基準)等で最適な数を決定してもよいし、無限次元LDAや混合ディリクレ過程等を用いることにより、自動で決定するようにしてもよい。
ここで、後述する事項の説明においても必要となるため、LDA等の潜在トピック分析によるクラスタリングについて簡単に説明する。図6は、当該潜在トピック分析によるクラスタリングの結果及びその意味合いを説明するための図である。
図6にて、クラスタリング対象の全データが(n歳のデータの場合であれば)D(n)であり、各対象者X(n歳とする)の特徴ベクトルV(X, n)を行ベクトルとして列挙した行列形式のデータである。潜在トピック分析においては、行列D(n)の各行はバグオブワーズとしての「文書u」であり、その列方向に並ぶ要素は各単語i1, i2, …, imの頻度数である。
本発明においては、各「文書u」が各「対象者X」に相当し、各単語i1, i2, …, imは、前述のように医療データの各評価項目に対応する所定の単語となる。そして、対象者Xに関して各単語i1, i2, …, imの具体的な頻度数を列挙することで特徴ベクトルV(X,n)が与えられる。
LDA等の潜在トピック分析では、文書u(本発明においては各対象者X)をこのようにバグオブワーズ、つまり単語とその出現頻度として取り扱い、文書においてそのトピックを推定する。例えば、「経済」トピックからは、「株価」、「増収増益」・・・といった単語が出現するだろうし、「スポーツ」トピックからは「野球」、「サッカー」といった単語が出現することになる。
これは、観測されたバグオブワーズ表現、つまり単語i(列成分)と文書u(行成分)との関係行列D(n)を、図6に示すように、単語i(列成分)とトピックk(行成分)の関係行列Φと、文書u(行成分)とトピックk(列成分)の関係行列θとの積に分解することを意味している。このトピックkを推定するのがLDA等の潜在トピック分析であり、当該分析の結果として図6に示すような「D(n)=θ×Φ」という行列の積としての分解結果が得られる。
このようにLDAに代表されるトピックモデルでは、各文書が固有のトピック比率を持ち、単語はこのトピック比率に従いトピックを選択したあと、そのトピックに固有の比率で生成されるという仮定をおいている。
なお、各文書u(=対象者X)のクラスタリング結果として解釈する場合には、関係行列θにおける当該文書uの行(行ベクトル)が、当該文書uの各トピックkへの所属確率を与えたものとなっているので、例えば、当該所属確率が最も高いトピックkに対応するクラスタに所属するものとして解釈することができる。
以上、(1)として、クラスタリング部3における小さい粒度である1年単位のクラスタリングを説明した。
(2)クラスタリング部3におけるN年単位の大きい粒度でのクラスタリングは次のようにして実施すればよい。ここで、説明のためにクラスタリング対象のN年の範囲をn歳〜n+N-1歳とする。当該範囲で参照可能なデータはD(n), D(n+1), …, D(n+N-1)である。
ここで、クラスタリング対象となる対象者については、一実施形態では、当該N年に渡る全データD(n), D(n+1), …, D(n+N-1)の全年齢においてデータが存在している対象者のみとしてもよい。別の一実施形態では、クラスタリング対象となる対象者は、当該N年に渡る全データD(n), D(n+1), …, D(n+N-1)のうち少なくとも1つの年齢においてデータが存在しているような対象者としてもよい。
クラスタリングを行うに際して、まず、各対象者Xの当該N年範囲(n歳〜n+N-1歳の範囲)における特徴ベクトルV(X, n〜n+N-1)を定める必要がある。当該N年範囲での特徴ベクトルは、N年範囲内での代表としての要素(代表要素M(X, n〜n+N-1))と、N年範囲内での微視的な変化としての要素(変化要素D(X, n〜n+N-1))と、を含めて構成することができる。すなわち、式で形式的に記すと以下の式(1)の通りである。
V(X, n〜n+N-1)=( M(X, n〜n+N-1), D(X, n〜n+N-1)) …(1)
当該N年範囲内での代表要素M(X, n〜n+N-1)及び変化要素D(X, n〜n+N-1)は、共に1以上の要素からなるベクトル量として構成することができ、それぞれ以下の(2−1)及び(2−2)のようにして定めることができる。
(2−1)代表要素については、例えば以下の式(2)のように、当該対象者XのN年範囲の各年齢n〜n+N-1の特徴ベクトルV(X, n)〜V(X, n+N-1)の平均として算出することができる。
ここで、当該代表の算出においては、以下[1]〜[7]の見出しで示すように種々の実施形態が可能である。当該[1]〜[7]はその一部または全部を任意に組み合わせて適用することができる。
[1] N年範囲のうち一部分のN'年分(N'<N)しかデータが存在しない対象者については、当該データが存在するN'年において平均を算出すればよい。[2] 上記の式(2)では特徴ベクトルすなわち単語頻度の形で平均を取っているので、小数点以下は四捨五入等で整数化すればよい。
[3] 前述のように、式(2)の特徴ベクトルV(X, k)は各単語i1, i2, …, imの頻度という形で与えられており、各単語は元の医療データにおける評価項目に、その頻度は当該評価項目の内容に応じて定まる。評価項目の中には例えば体重のように、数値としてその内容が評価されるものもある。数値として評価される評価項目については、当該元の数値においてN年分の平均を算出してから、当該N年分の平均評価数値を単語頻度の形式に変換してもよい。なお、単語頻度への変換処理は文書化部1と同様の処理により可能である。
[4] 式(2)の特徴ベクトルV(X, k)の要素ごとに、異なる手法で平均を算出してもよい。例えば上記[3]に関して、元の医療データにおいて体重のように数値として評価される要素に関しては元データの数値において平均値を算出してから当該平均値を単語頻度に変換するようにし、元の医療データにおいて数値以外(問診に対する回答としてのテキスト等)で評価される要素に関しては、式(2)のように変換された単語頻度において平均を算出するようにしてもよい。
[5] 式(2)の特徴ベクトルV(X, k)の要素の全部又は任意の一部分のうち、N年全体の平均ではなく直近のN''年(N''<N)の平均(直近データのみを用いるN''=1の場合を含む)として算出されるものがあってもよい。例えば元の医療データが健診データ及びレセプトデータで構成されている場合に、健診データに対応する要素についてはN年全体の平均として算出し、レセプトデータに対応する要素については直近の1年の値を用いるようにしてもよい。
[6] 式(2)では代表M(X, n〜n+N-1)を特徴ベクトルV(X, k)から算出するものとしているが、追加要素として、平均年齢を加えてもよい。すなわち、式で形式的に記すと以下の式(3)のように算出してもよい。
M(X, n〜n+N-1)=(特徴ベクトルV(X, k)から算出される諸要素, 平均年齢) …(3)
当該平均年齢は、当該対象者Xが当該N年範囲において医療データが存在する年齢の平均として算出すればよい。当該平均年齢の利用により、N年範囲の全部に医療データが存在している対象者と、その一部分のみに医療データが存在している対象者とを、クラスタリングの際に区別して扱うことが可能となる。
[7] 上記[6]の追加要素として平均年齢に代えて又は加えて、N年範囲のうち医療データが存在する年齢数(N年すなわちN回のうち何回、医療データがあるかをカウントした数)を用いてもよい。
(2−2)変化要素D(X, n〜n+N-1)については、以下に見出し[11]〜[13]として示すような種々の実施形態によって算出することができる。当該各実施形態は任意に組み合わせて適用可能であり、2つ以上の実施形態を組み合わせる場合は各実施形態で得られる要素を列挙することにより、ベクトル量としての変化要素D(X, n〜n+N-1)を構成することができる。
例えば、実施形態aで変化要素D(X, n〜n+N-1)としてベクトルVaが算出され実施形態bで変化要素D(X, n〜n+N-1)としてベクトルVbが算出される場合、実施形態a,bを組み合わせて適用するのであれば、ベクトル(Va, Vb)によって変化要素D(X, n〜n+N-1)を定めるようにすればよい。
なお、各実施形態による算出においては、変化要素を構成する各要素が実数として算出されることとなるが、文書化部1と同様の処理によって当該実数を単語頻度すなわち整数の形式へと変換したものが最終的に求める変化要素の値となる。当該単語頻度の形式への変換についてはいずれの実施形態においても共通して実施されるので重複する説明は省略し、実数として算出する処理の部分に関して説明する。
[11] 当該対象者XのN年範囲の各年齢n〜n+N-1の特徴ベクトルV(X, n)〜V(X, n+N-1)の各要素につき、その経年変化を回帰直線で表した際の直線の傾きを算出することで、変化要素D(X, n〜n+N-1)の各要素を構成するようにしてよい。
なお、当該特徴ベクトルの各要素は単語頻度という整数の形で与えられているので、前述の[3]と同様に、元の医療データにおいて体重等のように数値で与えられているような要素に関しては、当該元データの数値を対象として回帰直線を求め、その傾きを算出するようにしてもよい。
[12] 当該対象者XのN年範囲の各年齢n〜n+N-1の特徴ベクトルV(X, n)〜V(X, n+N-1)につきそれぞれ、前述の(1)の1年単位における小さい粒度でのクラスタリング結果における関係行列θの値θ(X, n)〜θ(X, n+N-1)を列挙して、その全体を変化要素D(X, n〜n+N-1)としてよい。すなわち、式で形式的に記すと以下の式(4)のようにしてよい。
D(X, n〜n+N-1)=(θ(X, n), θ(X, n+1), …, θ(X, n+N-2), θ(X, n+N-1)) …(4)
ここで、図6で説明したように、例えばθ(X, n)は当該n歳の全データD(n)をLDA等の潜在トピック分析によりクラスタリングした結果「D(n)=θ(n)×Φ(n)(ここではθ及びΦについてn依存を明示してある)」の関係行列θ(n)において、当該対象者X(=文書u)に該当する行ベクトルを抜粋したものであり、当該対象者Xのトピック比率を与えているものである。すなわち、対象者Xに関して入力された特徴ベクトルV(X, n)の、潜在トピック分析によるクラスタリング結果として出力された特徴ベクトルに相当するものがθ(X, n)である。なお、当該変化要素を算出するためのθ(X, n)〜θ(X, n+N-1)は、複数トピックを仮定したLDA等の潜在トピック分析で算出されるものではなく、単一トピックを仮定したk-means等により算出されたものとして求めてもよい。
[13] 当該対象者XのN年範囲の各年齢n〜n+N-1の特徴ベクトルV(X, n)〜V(X, n+N-1)につきそれぞれ、前述の(1)の1年単位における小さい粒度でのクラスタリング結果における関係行列θの値θ(X, n)〜θ(X, n+N-1)及び関係行列Φの値Φ(n)〜Φ(n+N-1)に基づいてトピック値T(X, n)〜T(X, n+N-1)を求め、当該一連のトピック値を変化要素D(X, n〜n+N-1)としてよい。すなわち、式で形式的に記すと以下の式(5)としてよい。
D(X, n〜n+N-1)=(T(X, n), T(X, n+1), …, T(X, n+N-2), T(X, n+N-1)) …(5)
ここで、例えばトピック値T(X, n)はデータD(n)より次の[手順1]〜[手順4]のように求めることができる。その他の年齢n+j (j=1, 2, …, N-1)のトピック値T(X, n+j)もデータD(n+j)より同様に求めることができる。
[手順1] まず、当該n歳の全データD(n)のクラスタリング結果「D(n)=θ(n)×Φ(n)」においては、関係行列Φ(n)をトピックk(図6参照)ごとにその単語頻度を正規化したものとして求めておく。
[手順2] 当該n歳の全データD(n)における平均の関係行列Φmean(n)を求める。ここで、当該平均の関係行列Φmean(n)とは、トピックkの数を1とした際の分解結果「D(n)=θ(n)×Φ(n)」におけるΦ(n)であって行ベクトルの形式を取り、単純に全データD(n)に現れる各単語i1, i2, …, imの頻度を正規化して求めたものに相当する。
[手順3] トピックkに対応する関係行列Φ(n)[k](関係行列Φ(n)の第k行にある行ベクトルに相当;図6参照)と、平均の関係行列Φmean(n)(前述のようにΦ(n)[k]と同サイズの行ベクトルとなる)と、の2乗距離(その他の距離でもよい)を、当該トピックk毎にd(n)[k]として以下の式(6)より求める。
ここで、平均の関係行列Φmean(n)とは[手順2]より当該n歳の全データD(n)における平均の単語比率を表しているので、上記求めたd(n)[k]は、分解結果「D(n)=θ(n)×Φ(n)」のΦ(n)においてトピックk毎に対応する行ベクトルΦ(n)[k]の当該平均単語比率からの乖離を表しており、トピックkの特別性を数値化したものに相当する。
[手順4] トピック値T(X, n)を以下の式(7)によって算出する。なお、θ(X, n)[k]は、[13]にて前述したθ(X, n)のk番目の要素すなわちk番目のトピック比率の値である。当該式(7)の算出はすなわち、対象者Xのクラスタリング結果としての特徴ベクトルθ(X, n)の特別性を、式(6)で求めたトピックk毎の特別性による重みづけ和として数値化することに相当する。当該算出されるトピック値T(X, n)が大きいほど対象者Xのn歳の健康状態が平均から乖離していることを意味する。
以上、変化要素D(X, n〜n+N-1)を算出する各実施形態を[11]〜[13]として説明した。ここで、これらについての補足を述べる。
[11]のみで変化要素を算出する場合は、(1)で説明した1年単位の小さい粒度でのクラスタリング結果は使わないので、当該クラスタリングは省略してもよい。
[12],[13]を適用する場合、N年範囲のうち一部分のN'年分(N'<N)しかデータが存在しない対象者Yについては、データ欠損している年齢n'のデータV(Y, n')を、当該年齢n'のデータV(X, n')が存在するその他の一連の対象者Xのデータより補完したうえで図5のようなクラスタリングを行い、当該補完データを入力としたクラスタリング結果より[12],[13]を適用するようにしてもよい。当該補完処理に関しては、本出願人の特願2014-148739号(データ解析装置及びプログラム;以下、特許文献5とする)に開示の技術を利用することができる。
特許文献5では、例えば対象者Yはn歳のデータV(Y, n)は存在するがn+1歳のデータV(Y, n+1)が存在せず欠損しているというような場合に、n,n+1歳の両データV(X, n),V(X, n+1)が存在するその他の対象者Xのうち当該n歳時点のデータV(X, n)がV(Y, n)に類似している、あるいはクラスタリング結果θ(X, n)がθ(Y, n)に類似しているような一連の対象者Xとその類似度とを求め、欠損データV(Y, n+1)を当該類似度を重みとして一連の類似対象者XのデータV(X, n+1)の線形和等として推定し、補完する。また、欠損が一連の年齢及び対象者に渡っている際は、逐次的に補完することで一連の年齢及び対象者において欠損のない状態を作る。
特許文献5よりも簡略化された手法として、対象者Yが年齢n'においてデータ欠損がある場合は、年齢n'においてデータが存在するその他の対象者Xの平均値によって補完するようにしてもよい。当該補完はV(Y, n')を対象として補完後にクラスタリングを行ってもよいし、補間前のデータでクラスタリングを行ってθ(Y, n')を補完するようにしてもよい。
以上、クラスタリング部3ではN年範囲での大きな粒度での各対象者Xの特徴ベクトルV(X, n〜n+N-1)を定めると、当該特徴ベクトルを用いて各対象者Xをクラスタリングする。当該クラスタリングの手法に関しては、(1)で説明した小さい粒度である1年単位のクラスタリングと同様に、LDA等の潜在トピック分析を利用することができる。
モデル構築部4は、クラスタリング部3で上記(1)のようになされた1年単位の小さい粒度でのクラスタリング結果と、上記(2)のようになされたN年単位の大きい粒度でのクラスタリング結果と、のそれぞれにおいてクラスタ間の遷移確率を求め、クラスタ及び遷移確率を紐付けたものとして各粒度における予測モデルを出力する。
本発明においては特に、N年単位の大きい粒度でのクラスタリング結果による予測モデルにより、前述の特許文献4における1年単位の小さい粒度でのクラスタリング結果よる予測モデルにおける長期予測の発散を防止する、あるいは発散せず所定精度で予測可能な期間を延長することができる。
ここで、当該大きい粒度でのクラスタリングにおいては前述のように、N年範囲での代表要素M(X, n〜n+N-1)のみならず変化要素D(X, n〜n+N-1)も特徴ベクトルとして用いている。従って、大きな粒度での予測モデルを構築しながらも、小さな粒度における微視的変化の情報も考慮することが可能となる。このためN年単位で平均すると同じ健康状態となる対象者同士であっても、N年範囲内において健康状態が改善している、一定している、あるいは悪化しているといった違いを区別してクラスタリングがなされるので、精度の高い予測モデルを得ることができる。
モデル構築部4では具体的にクラスタ間の遷移確率を以下に見出し[21],[22]として示す各実施形態で計算することができる。これら実施形態は、遷移確率を計算しようとしているクラスタ間の類似性を求め、類似性が高いほど遷移確率も高いものとして計算を行う手法として総括することができる。なお、遷移確率を計算するクラスタは1年単位の小さい粒度の場合を例として説明するが、図4で例示したように構造は同じであるので、N年単位の大きい粒度のクラスタ間においても同様に計算することができる。
[21] 一実施形態では、年齢n歳のある健康状態iのクラスタC(n, i)から年齢n+1歳のある健康状態jのクラスタC(n+1,j)への遷移確率P[C(n, i)→C(n+1, j)]を以下の式(8)により求めることができる。なお、式(8)中では積集合の演算子「∩」は、両集合(両データ)に共通して含まれる対象者を抽出する演算子とする。データGの元の数すなわちデータGを構成している対象者の数をnum(G)とする。
すなわち、遷移確率は、式(8)の分子に示されるように、遷移元のクラスタC(n, i)と遷移先のクラスタC(n+1, j)とに共通して含まれる対象者の数に比例するように計算することができる。
なお、式(8)の分母は、当該遷移元のクラスタC(n, i)から遷移可能な、インデクスj=1, 2, …, N(n+1)で指定される全ての遷移先のクラスタC(n+1, j)への遷移確率を足して1となるようにするための規格化項である。
[22] 上記の式(8)のようにクラスタ間で共通の対象者数に基づく手法とは別の手法により、遷移確率を計算する一実施形態として次のようにしてもよい。当該実施形態においては、図6で説明したようなクラスタリングの際に得られる行列分解「θ×Φ」の情報を利用する。
当該実施形態では、各対象者のクラスタへの所属を式(8)の実施形態とは少し異なる解釈で与えたうえで、式(8)と同様にクラスタ間の共通対象者を算出することにより、遷移確率を求める。すなわち、式(8)の実施形態では、クラスタリングにおいて各対象者が、その最大トピック比率を与えるいずれか1つのクラスタに所属しているという解釈のもと、クラスタ間の共通対象者に比例する形で遷移確率を計算した。これに対して、当該実施形態においては、各対象者は、そのトピック比率に比例する寄与率で、それぞれのクラスタに分かれて所属しているという解釈のもとで、式(8)の考え方と同じく、クラスタ間の共通受診者に比例する形で遷移確率を計算する。
図7は、モデル構築部4の当該実施形態における遷移確率の計算を説明するための例を示す図である。図7では、[1]に示すように、クラスタリング部3によるクラスタリング結果の一部分(年代上の一部分)として、2つの年代n,n+1において、C(n, 1), C(n, 2), C(n ,3)及びC(n+1, 1), C(n+1, 2), C(n+1, 3)と分類されている。また、[2-1], [2-2]に示すように、全受診者はA〜Eの5人であり、それぞれ示すようなθ行列の行ベクトルとしてのトピック比率が得られている。年代nについてはトピックの各要素は(k1, k2, k3)であり、年代n+1についてはトピックの各要素は(k'1, k'2, k'3)である。そして、[3-1], [3-2], [3-3]に、以下に説明する遷移確率の計算例が示されている。
なお、当該実施形態においても、図7の[1]に示すように、クラスタリング部3の結果の利用の仕方として、クラスタC(n, 1), C(n, 2), C(n ,3)はそれぞれトピックk1, k2, k3に対応するクラスタであり、C(n+1, 1), C(n+1, 2), C(n+1, 3)はそれぞれトピックk'1, k'2, k'3に対応するクラスタである、とみなす点は、その他の実施形態の場合と共通である。
そして、前述のように当該実施形態においては、例えば図7の[2-1]に関して、トピック比率が(0.8, 0.1, 0.1)である対象者AはクラスタC(n, 1)に「0.8」の割合で所属し、クラスタC(n, 2)に「0.1」の割合で所属し、クラスタC(n ,3)に「0.1」の割合で所属するという解釈を取る。当該解釈のもと、遷移確率の計算は次の[手順11]〜[手順13]のようにして行うことができる。
[手順11] まず、遷移元となる各クラスタに属する対象者の総数は、当該クラスタに対応するトピックのトピック比率の、全ての対象者についての和とする。その計算例は、図7の[3-1]に示す通りであり、例えばトピックk1に対応するクラスタC(n, 1)の所属人数num(C(n, 1))を、示される通り各受診者A〜Eのそれぞれの寄与(当該トピックk1の比率の値として与えられる寄与)の和として、3.7人と計算することができる。
[手順12] 次に、遷移元クラスタから遷移先クラスタへのクラスタ間の遷移人数は、遷移元クラスタに対応するトピックのトピック比率と遷移先クラスタに対応するトピックのトピック比率との積の、全ての受診者についての和とする。その計算例は、図7の[3-2]に示す通りである。例えば、トピックk1のクラスタC(n, 1)からトピックk'1のクラスタC(n+1, 1)への遷移人数num(C(n, 1)→C(n+1, 1))を、[3-2]に示すように1.98人と計算することができる。その他のクラスタ間遷移人数も全く同様に計算することができる。
[手順13] 最後に、式(7)の実施形態と同様の方針、すなわち、クラスタ間の共通対象者数に比例する形で、遷移確率を計算する。ただし、クラスタ所属人数及びクラスタ間遷移人数は、以上述べたトピック比率を用いた手法により計算したものを用いる。その計算例は、図7の[3-3]に示す通りである。当該例では、クラスタC(n, 1)からの全遷移先として3つのクラスタC(n+1, j)(j=1, 2, 3)への遷移確率P[C(n, 1)→C(n+1, j)]がそれぞれ図中に記す通りに計算されている。
なお、確率としての規格化も、特に年代間でデータ欠損(医療データが存在する対象者の過不足)があるような場合は式(7)の実施形態に分母の項の計算と同様にして、ある一つの遷移元クラスタから可能な全ての遷移先クラスタへの遷移確率の和が1となるように、適宜実施すればよい。規格化に関しては、図7の[3-1]〜[3-3]の計算例に示すように、異なる年代間にてデータ欠損がなく、各対象者のトピック比率をその総和が1となるよう予め規格化しておけば、計算される遷移確率も自ずと規格化されることとなる。
なお、図7の例では遷移元の年齢n歳のクラスタ数と遷移先の年齢n+1歳のクラスタ数とが共通の場合を例としたが、クラスタ数が異なっていても全く同様に遷移確率を計算することができる。
以上、粒度考慮モデル構築部40において分割幅設定部5の設定により1年単位の小さい粒度とN年単位の大きい粒度とからなる、2階層での予測モデルを構築する場合を説明した。同様にして、2階層よりも多い任意数の階層の粒度からなる予測モデルを構築することも可能である。
図8は、当該より多い階層の粒度からなる予測モデルの例として、4階層の予測モデルを構築する際の分割幅の例を示す図である。図8の例では、n=20〜60(20歳〜60歳)の範囲の40年間を、[1]では1年単位で40分割し、[2]では5年単位で8分割し、[3]では10年単位で4分割し、[4]では20年単位で2分割することで、4階層の粒度からなる予測モデルが構築される。
このように多階層で予測モデルを構築する場合も、2階層の場合と同様に、粒度の小さい側から順次、クラスタリング部3でクラスタリングを行うと共にモデル構築部4でクラスタ間の遷移確率を計算することで予測モデルを構築していけばよい。
図7の例であれば、最初に、[1]の最小粒度の1年単位の予測モデルと[2]の2番目の粒度の5年単位の予測モデルとは、(N=5とおくことで)以上説明したのと同様にして構築することができる。次に、3番目以降のM番目(M=3,4)の粒度の予測モデルを構築する際は、既に構築されているM-1番目の予測モデルを粒度の小さいモデルとみなし、当該構築しようとするM番目の予測モデルを粒度の大きい予測モデルとみなすことで、以上説明した2段階の予測モデル構築の場合と全く同様にして構築することができる。
すなわち、M番目の予測モデルの粒度(分割幅)をP(M)年と書くことにすると、次のように2段階の予測モデル構築と同様の処理を再帰的に実施することで、当該M番目の予測モデルを構築できる。まず、2段階の粒度の予測モデル構築にて説明したN年範囲の代表要素M(X, n〜n+N-1)に相当するものを当該P(M)年範囲において算出し、変化要素D(X, n〜n+N-1)に相当するものを既にモデル構築されている1段階粒度の低いP(M-1)年単位で変化を算出して求めることにより、クラスタリング部3にて当該P(M)年単位でのクラスタリングを行うことができる。また、モデル構築部4による遷移確率の計算はM段階のいずれの粒度においても同様に行うことができる。こうして、P(M)年単位の予測モデルが構築される。
以上の説明では、分割幅設定部5では各粒度において等しい分割幅を設定する場合を例としたが、ある粒度において異なる分割幅が任意の順番で存在していてもよい。ただし、小さい粒度のM-1番目の分割幅の設定が、大きい粒度のM番目の分割幅の設定に包含されることで、階層性を有するような構造となっていることが好ましい。
図9は、分割幅が各粒度において等しくない場合の例を示す図である。図9の例では、n=20〜50(20歳〜50歳)の範囲の30年間においてM=2段階の粒度からなる予測モデルを[1]の小さい粒度及び[2]の大きい粒度として構築する場合の分割幅が例示されている。
図示するように、大きい粒度の[2]では(a),(b),(c),(d)と4区間に分割され、それぞれ長さが10年、5年、9年、6年である。これに包含される形で、小さい粒度の[1]では(a),(b),(c),(d)がそれぞれ以下のように分割されている。
(a)の10年=1年+2年+2年+3年+2年
(b)の5年=2年+3年
(c)の9年=3年+1年+3年+2年
(d)の6年=3年+3年
なお、最小の粒度における分割幅が1年(医療データの取得単位の例としての1年)よりも長い場合は、当該分割幅におけるクラスタリング部3によるクラスタリング対象としての特徴ベクトルは、2段階の粒度の予測モデル構築にて説明したN年範囲の代表要素M(X, n〜n+N-1)及び変化要素D(X, n〜n+N-1)を求めることで、前述の式(1)の形で定めてもよいし、変化要素D(X, n〜n+N-1)は省略して代表要素M(X, n〜n+N-1)のみの形で定めてもよい。
分割幅設定部5では、図8や図9で説明したような粒度に関しての階層性を有する所定の分割幅の設定を、ユーザ指定に従って実施することができる。ここで、最終的に構築される予測モデルの精度が高くなるような最適な分割幅を自動決定する手法も可能である。以下、当該手法について説明する。当該手法を用いる場合も、最適な分割幅を決定するための候補としての種々の分割幅は、ユーザ指定によって与えておくものとする。
当該自動決定においては、ある分割設定DCが区間1,区間2,…, 区間mであったとした場合に、当該分割設定DCが健康に関する特定の評価指標について適切であるか否かを周知のAIC(赤池情報量基準)を用いた独立性検定によって数値として求め、一連の分割設定DCの中から最小のAIC値(従属性が最も高い)を与えるような分割設定DCを最適なものとして決定する。なお、当該最小のAIC値は、以下詳細を説明するように、「AIC(DM)-AIC(IM)」又は「AIC(DM)」に関する最小値として求めることができる。
当該決定される「最適」な分割設定DCの意味合いは次の通りである。すなわち、当該「最適」な分割設定DCを構成する区間1,区間2,…, 区間mは、当該各区間i(i=1, 2, …, m)内においては健康状態が均質であるものとして決定されるという点に関して、「最適」である。本発明においては特に、当該健康状態が均質であると考えられる一連の年齢を区間i(i=1, 2, …, m)として一括して扱うことで、長期予測の際も精度を保つモデルの構築が可能となる。
AIC値が最小となり最適であるものとして決定された分割設定DCを元に粒度を考慮した予測モデルを構築する際は、以下に見出し[31],[32]として区別して示すような各実施形態が可能である。
[31] 一実施形態では、AIC値が最小の側からi番目となる分割設定をDC[i]とすると、DC[1]を最小粒度として設定した後、2番目の大きさの粒度の分割設定は、2番目以降の分割設定DC[j](j=2, 3, …)のうちDC[1]を包含する関係にあるもの(すなわち、階層性を有するモデルを構築可能なもの)の中から最小のAIC値を有するものを選ぶことで階層モデルを構築することができる。3番目以降の大きさの粒度の分割設定に関しても同様であり、当該選ぶことができなくなった時点で階層モデルの構築を完了すればよい。
例えば、最適なものを決定する候補としての一連の分割設定DCを、図8の[1]の1年区分の区間、[2]の5年区分の区間、[3]の10年区分の区間及び[4]の20年区分の区間、の全4通りとして与えた場合に、算出したAIC値が小さい方から順に[2],[4],[3],[1]であったとする。この場合、構築される階層モデルは、AIC値が最小となる[2]の5年区分の区間を最小粒度とし、2番目の大きさの粒度を[4]の20年区分の区間で与えた2階層のモデルとすることができる。ここで、AIC値が3番目の[3]は10年区分の区間であり、2番目の[4]の20年区分の区間を包含するものではないため、2階層の時点でモデル構築を完了する。
[32] 一実施形態では、2階層で予測モデルを構築し、AIC値が最小となる分割設定DC[1]を粒度の大きい予測モデルとし、粒度の小さいモデルについては1年毎に全て区切る分割設定を用いる、あるいは、2番目以降の分割設定DC[j](j=2, 3, …)のうちAIC値が最小であり且つDC[1]によって包含される分割設定を用いることができる。
以下では一般に区間1,区間2,…, 区間m(実例は例えば図8における[1]〜[4]等となる)として与えられる分割設定DCに関して、その健康に関する特定の評価指標についてのAIC値の算出による独立性検定の詳細を説明する。
具体的には、分割設定DCにおける評価指標の「該当」「未該当」の度数の分布を参照して、赤池情報量基準(AIC)を算出することができる。ここで、独立性検定におけるAIC値の算出のためには、統計分野で周知のクロス集計表を作成する必要があるので、まず当該クロス集計表を作成してから、AIC値を算出する。図10は、当該クロス集計表の作成を説明するための図である。
図10にて(1)は、分割設定DCにおける、各区間1,区間2,…, 区間mの度数を表として示している。ここで、ある区間iが2年以上に渡る場合で、ある対象者Xについて2年分以上(すなわち2回分以上)のn年分のデータが存在する場合、当該対象者Xは1人ではなく、当該データが存在するn人としてカウントするものとする。すなわち、同一対象者Xであっても異なる年齢においてデータが存在すれば、別人として度数をカウントする。当該カウントのもと、区間1,区間2,…, 区間m全体としての対象データのサンプル総数をNとする。従って、区間iの度数をkiとして、「k1+k2+…+km=N」である。
そして、図10の(2)に示すのが、当該対象データの各サンプル(対象者X;以下、統計分野の用語に従い、「サンプル」と称する。前述のように同一対象者Xでも年齢が異なれば別サンプルとなる。)について予め与えられる評価指標が該当するかしないかの区別を、分割設定DCの各区間i(i=1, 2, …, m)に属するサンプルにおいて調べ、統計分野における周知のクロス集計表としたものである。すなわち、(2)は、各区間iにおける評価指標に関して「該当」であるサンプルの度数n1iと「未該当」であるサンプルの度数n2iと、をi=1, 2, …,mに渡って与えることで、2×mのクロス集計表となっている。
なお、(2)の表では、各区間iの度数がkiであり、属するサンプルは評価指標に関して「該当」又は「未該当」のいずれかであるので、周辺度数として「ki=n1i+n2i」の関係がある。また、「h」は対象データ全体のうちの「該当」のサンプル数であり、周辺度数として、「h=n11+n12+…+n1m」の関係がある。さらに、対象データ全体のうちの「未該当」のサンプル数は「N-h」であり、周辺度数として、「N-h=n21+n22+…+n2m」の関係がある。
なお、予め与える評価指標については、健康に関する所定基準とすればよい。例えば、「生活習慣病を発症しているか否か」を評価指標とし、各サンプルすなわち年齢nの対象者Xにつき、元の医療データを参照して当該評価指標に関して「該当」又は「未該当」を判定すればよい。また、前述の決定される最適な分割設定DCを構成する区間i(i=1, 2,…, m)は、当該評価指標に関して「均質」であるものとして決定されることとなる。
図10のようなクロス集計表を参照して、当該分割設定DCのAIC値を、次のいずれかの手法の値として求めることができる。第一手法では、当該クロス集計表に対して従属モデルを適用することにより、以下の[式10]のような従属モデルのAIC値AIC(DM)[ここでDMはDependent Modelの略である]として求める。第二手法では、さらに、当該クロス集計表に対して独立モデルを適用して、以下の[式20]のような独立モデルのAIC値AIC(IM)[ここでIMはIndependent Modelの略である]を求めたうえで、[式30]のように、従属モデルのAIC値から独立モデルのAIC値を引いた差の値として、求める。
なお、[式10]等においてMLL(DM)は、従属モデルにおける最大対数尤度であって、[式10-1]のような値として求めることができる。また、[式20]等において、MLL(IM)は、独立モデルにおける最大対数尤度であって、[式20-2]のような値として求めることができる。なお、上記の各式における文字は、図10のクロス集計表において説明した通りであり、以降説明する各式においても同様である。各モデルにおける最大対数尤度MLL(DM)及びMLL(IM)が上記のように算出されることについては、以下に説明する通りである。
図11は、[式10]及び[式10-2]として示した従属モデルにおける算出を説明するための、図10の(2)のクロス集計表に対応する従属モデルにおける確率の表である。当該表に示されている確率により、以下のように算出がなされる。
まず、従属モデルの確率変数は以下の通りである。
一方、図11に示された2m個の全てが自由に動かせるわけではなく、以下の制約がある。
従って、従属モデルの自由度は2m-1であり、AICの定義(AIC=-2×MLL+2×自由度)より、[式10]の2*(2m-1)の項が得られる。さらに、上記確率変数より対数尤度LLを計算すると、以下のようになる。
上記対数尤度LLを最大にするときの条件は以下である。
上記最大とする条件より、以下が得られる。
上記と同様にして、さらに
等が得られる。そこで、
とすると、
等となるので、それぞれを足すと、
となるから、以下の場合が最尤推定となる。
従って、上記の値をLLに代入することで、その最大値として前述の[式10-2]が得られる。
図12は、[式20]及び[式20-2]として示した独立モデルにおける算出を説明するための、図10の(2)のクロス集計表に対応する従属モデルにおける確率の表である。当該表に示されている確率により、以下のように算出がなされる。
まず、図10の周辺度数kmと、対応する図12の周辺確率qmと、において、以下のような制約がある。
従って、自由に動かせるのはq1〜qm-1とpとであるから、パラメータの自由度は(m-1)+1=mであって、AIC算出の定義より、[式20]の2×mの項が得られる。また、独立モデルの確率変数は以下の通りとなる。
従って、その対数尤度LLは以下の通りとなる。
対数尤度の最大値を与える条件を求めるべく、これをp、q1・・・で偏微分してゼロに等しいとすることにより、以下等の一連の計算ができる。
従って、
となり、また、
とすると、
等となるので、それぞれ足して、
となり、
となるから、最大尤度は
等において得られることとなる。従って、上記の値をLLに代入することで、最大値としての[式20-2]が得られる。
以下、本発明における補足的事項(1)〜(3)を説明する。
(1)粒度考慮モデル構築部40により構築された予測モデルで年齢nの対象者Xに関して実際に将来の健康の予測を実施する際は、当該予測を実施するユーザの設定に応じていずれの粒度で予測を行ってもよいが、長期に渡る精度を保つという観点からは最も粒度の大きな予測モデルを用いることが好ましい。そして、より微視的な健康の変化を確かめたい場合に、適宜粒度の小さい側の予測モデルによる予測結果を参照するという形で、ユーザは解析や検討を行うことが可能である。
ここで、いずれの粒度のモデルで予測する際も、予測結果としてのクラスタ間の遷移(年齢が上がる方向への遷移)は、以下のURL等に開示され周知のビタビアルゴリズムを適用することで、最大確率を与える経路として決定することができる。なお、前述のようにクラスタ間の遷移確率はモデル構築部4において計算されているので、ビタビアルゴリズムの適用が可能である。
(URL)http://www.yobology.info/text/viterbi/viterbi.htm
なお、ビタビアルゴリズムの利用に際しては、状態間の遷移確率から最大確率に対応する経路を計算する他にも、遷移確率を尤度に変換して、最大尤度に対応する経路を計算してもよいし、遷移確率を距離に変換して、最短距離に対応する経路を計算してもよいし、その他の評価指標によって最適経路を計算するようにしてもよい。
また、予測対象となる年齢nの対象者Xについては当該年齢nで属するクラスタ(予測の初期状態に対応するクラスタ)を決定する必要がある。当該年齢nの対象者Xのデータが文書化部1に入力される医療データに含まれていれば、クラスタリング部3の計算により当該属するクラスタは既に求まっているので、再度決定する必要はない。当該年齢nの対象者Xのデータが文書化部1に入力される医療データに含まれていない場合であれば、クラスタリング部3がクラスタリングを行った際と同様の処理により属するクラスタを決定すればよい。
なお、本発明にて構築された予測モデルは、予測モデル構築装置10に入力された当初の欠損を有する医療データ(図1の[2]で説明したような現実に利用可能な医療データ)の観点から見ると、健康状態が似通った対象者同士をチェーン状(鎖状)に連結することにより欠損を補って構築された予測モデルという意味合いを有する。
(2)本発明において、医療データは文書化部1にて文書化するものとして説明した。文書化処理(すなわち、単語頻度形式としてのデータへの変換)を予め行っておくものとすれば、文書化部1は省略されてもよい。
(3)本発明は、コンピュータを予測モデル構築装置10として機能させるプログラムとしても提供可能である。当該コンピュータは、CPU(中央演算装置)、メモリ及び各種I/Fといった周知のハードウェアで構成することができ、当該プログラムを読み込んで実行するCPUが予測モデル構築装置10の各部として機能することとなる。