JP7365496B2

JP7365496B2 - 脳機能結合相関値のクラスタリング装置、脳機能結合相関値のクラスタリングシステム、脳機能結合相関値のクラスタリング方法、脳機能結合相関値の分類器プログラム、脳活動マーカー分類システムおよび脳機能結合相関値のクラスタリング分類器モデル

Info

Publication number: JP7365496B2
Application number: JP2022514450A
Authority: JP
Inventors: 雄人柏木; 智磯徳田; 雄史高原; 光男川人; 歩山下; 宙人山下; 雄希酒井; 潤一郎吉本
Original assignee: ATR Advanced Telecommunications Research Institute International
Current assignee: ATR Advanced Telecommunications Research Institute International
Priority date: 2020-04-06
Filing date: 2021-04-02
Publication date: 2023-10-19
Anticipated expiration: 2041-04-02
Also published as: CN115484864A; WO2021205996A1; US20230107263A1; EP4134884A4; JPWO2021205996A1; EP4134884A1

Description

この発明は、複数の装置において脳機能画像法により計測される脳機能結合相関値のパターンをクラスタリングする技術、より特定的には、脳機能結合相関値のクラスタリング装置、脳機能結合相関値のクラスタリングシステム、脳機能結合相関値のクラスタリング方法、脳機能結合相関値の分類器プログラム、脳活動マーカー分類システムおよび脳機能結合相関値のクラスタリング分類器モデルに関する。この出願は、２０２０年４月６日出願の日本出願特願２０２０－０６８６６９号に基づく優先権を主張し、前記日本出願に記載された全ての記載内容をここに参照により援用する。

（データ駆動型のクラスタリング方法）
近年の人工知能技術、特に、データ駆動型人工知能技術の進展により、音声認識、翻訳、画像認識などの分野で、一部、人間の能力に匹敵し、あるいは、一部の領域では、人間の能力を凌駕するような応用も実現されている（たとえば、特許文献１）。

医療技術の分野においても、画像診断などで、深層学習などの機械学習が利用されることが増えてきた。深層学習は、多層のニューラルネットワークを用いた機械学習であり、画像認識の分野では、深層学習の１つである畳込みニューラルネットワーク（Convolutional Neural Network、以下、「ＣＮＮ」という）を用いた学習方法が従来の方法に比べて非常に高い性能を示すことが知られている（たとえば、特許文献２）。

たとえば、大腸がんの内視鏡による画像診断などでは、診断の正確度が人間を上回るような診断機器が実用化されるに至っている（非特許文献１）。

ただし、これらの人工知能技術は、機械学習の分類上は、正解データと入力データ（たとえば、画像データ）との組を大量に準備して、これを入力として、人工知能に学習処理をさせるという、いわゆる「教師あり学習」の範疇に入っているものがほとんどである。

一方で、データ駆動型人工知能の応用用途としては、与えられたデータを、その特徴量に基づいて、いくつかのクラスタに分類するというタスクの実行もある。この場合は、正解データの存在しない、いわゆる「教師なし学習」や、少量の「正解ラベル付きの学習データ」による学習と大量の「正解ラベルなし学習データ」による学習とを組み合わせた「半教師あり学習」などが知られている（たとえば、特許文献３）。

たとえば、特許文献３では、「半教師あり学習は、比較的少ないラベル付きデータとラベル無しデータとに基づいて学習を行う学習方法であり、例えばラベル付きデータ（状態データＳと判定データＬを含む教師データＴ）を用いて分類を行う学習モデルを生成し、該学習モデルとラベル無しデータ（状態データＳ）とを用いて、該学習モデルに対して追加で学習を行うことで学習の精度を向上させるブートストラップ法や、ラベル付きデータ及びラベル無しデータのデータ分布に基づいてグループ分けすることで分類器としての学習モデルを生成するグラフベースアルゴリズム等を含む。」とされている。ただし、この例にもあるとおり、「半教師あり学習」においては、教師データが存在するのは、少数の学習データであり、これで、まずは分類器を生成した後、大量の「正解ラベル無し学習データ」用いることで、この分類器自身に対して再学習することなどが前提となっていることになる。

（バイオマーカー）
以下では、人工知能技術による判別やクラスタリングを適用する分野として、医学分野を例にとる。

生体内の生物学的変化を定量的に把握するため、生体情報を数値化・定量化した指標のことを「バイオマーカー」と呼ぶ。

ＦＤＡ(米国食品医薬品局)はバイオマーカーの位置づけを、「正常なプロセスや病的プロセス、あるいは治療に対する薬理学的な反応の指標として客観的に測定・評価される項目」としている。また疾患の状態や変化、治癒の程度を特徴づけるバイオマーカーは、新薬の臨床試験での有効性を確認するためのサロゲートマーカー(代用マーカー)として使われる。血糖値やコレステロール値などは、生活習慣病の指標として代表的なバイオマーカーである。尿や血液中に含まれる生体由来の物質だけでなく、心電図、血圧、ＰＥＴ（positron emission tomography）画像、骨密度、肺機能なども含まれる。またゲノム解析やプロテオーム解析が進んできたことによって、ＤＮＡやＲＮＡ、生体蛋白等に関連したさまざまなバイオマーカーが見出されている。

バイオマーカーは、疾患にかかった後の治療効果の測定だけでなく、疾患を未然に防ぐための日常的な指標として疾患の予防に、さらに副作用を回避した有効な治療法を選択する個別化医療への応用が期待されている。

たとえば、肺疾患に対して、遺伝子情報を用いて、疾患への罹患の可能性を判断するためのバイオマーカーについて開示がある（特許文献４）。特許文献４では、「バイオマーカー」または「マーカー」とは、「その生体システムの生理学的な状態の特徴を示すものとして、客観的に測定することが可能な生物学的な分子である」とされている。そして、この特許文献４では、「通常、バイオマーカー測定値は典型的にはタンパク質またはポリペプチドである発現産物の定量的な測定に関する情報である。本発明は、バイオマーカー測定値をＲＮＡ（翻訳前）レベルでまたはタンパク質レベル（翻訳後修飾も含まれ得る）で決定することを想定している」とされている。そして、特許文献４では、このようなバイオマーカー測定値に対する「分類システム」として使用する分類器としては、決定木、ベイジアン分類器、ベイジアン信念ネットワーク、ｋ-最近傍法、事例ベース推論、およびサポートベクトルマシンなどが例示されている。

一方で、神経・精神疾患の場合、現状の診断は、ＤＳＭ－５（精神障害の診断・統計マニュアル（Diagnostic and Statistical Manual of Mental Disorders）第５版）などに基づいて、いわゆる症候に基づく診断となっていることもあり、生化学的もしくは分子遺伝学的観点から客観的な指標として利用可能な分子マーカーなども研究されているものの、検討段階というべき状況である。

もっとも、ＮＩＲＳ（Near-infraRed Spectroscopy）技術を用いて、生体光計測により計測されたヘモグロビン信号から特徴量に応じて、統合失調症、うつ病などの精神疾患について分類を行う疾患判定システムなどについては報告がある（特許文献５）。

（脳活動に基づくバイオマーカー）
一方で、いわゆる画像診断の分野においては、上述したような「生物学的な分子」というバイオマーカーの概念とは異なり、「画像バイオマーカー」と呼ばれるものも存在する。たとえば、脳神経領域での分子イメージングに、ＰＥＴを用いて、神経伝達機能や受容体機能解析を行おうとする試みもある。

さらに、核磁気共鳴映像法（ＭＲＩ:Magnetic Resonance Imaging）では、血流量の変化に応じて、検出される信号に変化が現れることを用いて、外部刺激等に対する脳の活動部位を視覚化することも可能である。このような核磁気共鳴映像法を、特に、ｆＭＲＩ（functional ＭＲＩ)と呼ぶ。

ｆＭＲＩでは、装置としては通常のＭＲＩ装置に、さらに、ｆＭＲＩ計測に必要なハードおよびソフトを装備したものが使用される。

ここで、血流量の変化がＮＭＲ信号強度に変化をもたらすのは、血液中の酸素化および脱酸素化ヘモグロビンは磁気的な性質が異なることを利用している。酸素化ヘモグロビンは反磁性体の性質があり、周りに存在する水の水素原子の緩和時間に影響を与えないのに対し、脱酸素化ヘモグロビンは常磁性体であり、周囲の磁場を変化させる。したがって、脳が刺激を受け、局部血流が増大し、脱酸素化ヘモグロビンが変化すると、その変化分をＭＲＩ信号として検出する事ができる。このような被験者への刺激は、たとえば、視覚による刺激や聴覚による刺激、あるいは所定の課題（タスク）の実行等が用いられることが一般的である。

そして、脳機能研究においては、微小静脈や毛細血管における赤血球中の脱酸素化ヘモグロビンの濃度が減少する現象（ＢＯＬＤ効果）に対応した水素原子の核磁気共鳴信号（ＭＲＩ信号）の上昇を測定することによって脳の活動の測定が行われている。

このように、ｆＭＲＩ装置により計測される脳活動を反映した血中酸素濃度依存信号をＢＯＬＤ信号（Blood Oxygen Level Dependent Signal）と呼ぶ。

特に、人の運動機能に関する研究では、被験者に何らかの運動を行わせつつ、上記ｆＭＲＩ計測によって脳の活動を測定することが行われている。

ところで、ヒトの場合、非侵襲的な脳活動の計測が必要であり、この場合、ｆＭＲＩデータから、より詳細な情報を抽出できるデコーディング技術が発達してきている。特に、ｆＭＲＩが脳におけるボクセル単位（volumetric pixel : voxel）で脳活動を解析することで、脳活動の空間的パターンから、刺激入力や認識状態を推定することが可能となっている。

さらに、このようなデコーディング技術を発展させた技術として、特許文献６には、神経・精神疾患に対して、脳機能画像法による「診断バイオマーカー」を実現するための脳活動解析方法が開示されている。この方法では、健常群、患者群において測定された安静時機能結合的ＭＲＩのデータから、それぞれの被験者について、所定の脳領域間の活動度の相関行列（脳機能結合パラメータ）を導出する。被験者の疾患／健常ラベルを含む被験者の属性と相関行列とについて正則化正準相関解析により特徴抽出が行われる。正則化正準相関解析の結果に基づいて、スパースロジスティック回帰(Sparse Logistic Regression: SLR)による判別分析によりバイオマーカーとして機能する判別器が生成される。このような機械学習の技術により、安静時のｆＭＲＩデータから導き出される脳領野間の結合にもとづいて神経疾患の診断結果を予測可能であることが示された。しかも、その予測性能の検証は１つの施設において計測された脳活動のみではなく、他の施設で計測された脳活動に対しても、ある程度の汎化が可能であることが示された。

さらに、このような「診断バイオマーカー」について、汎化性能を一層向上させるための技術改良もなされている（特許文献７）。

また、最近では、米国のヒューマンコネクトームプロジェクトのように、大規模な脳画像データを得て共有することは、基礎的な神経科学研究と、精神疾患の診断および治療のような臨床応用の間のギャップを埋めることにとって重要な意味を持つ、と認識されている（非特許文献２）。

２０１３年には、日本の国立研究開発法人日本医療研究開発機構が、８つの研究所が２，２３９のサンプルおよび５つの疾患を含む複数サイトの安静時の機能的磁気共鳴(安静時の機能ＭＲＩ)データを収集し、公にＳＲＰＢＳ（Strategic Research Program for Brain Sciences，https://www.amed.go.jp/program/list/01/04/001_nopro.html）の複数サイトの複数疾患のデータベース(https://bicr-resource.atr.jp/decnefpro/)を通じてそれらを共有するという、デコーディッドニューロフィードバック(DecNef)プロジェクトを組織した。このプロジェクトは、完全に独立したコホートに汎化することができる、いくつかの精神疾患の、安静時の機能的なコネクティビティ(安静時の機能結合ＭＲＩ)に基づくバイオマーカーを同定している。

このように健常群と疾患群の診断については一定の成果が得られつつある。ところで、疾患群の中において、たとえば、一般に「うつ病」と診断されている患者群も、実は、複数のサブタイプに分かれていることが知られている。たとえば、通常の「抗うつ薬」の投与で寛解に向かう患者群がある一方で、寛解しにくい「治療抵抗性」の患者群があることなどが知られている。

このような「うつ病」の患者に対して、上述したような「脳機能結合パラメータ」に対して、データ駆動型人工知能によるクラスタリングを適用して、分類しようとする試みも存在し、一定の傾向が存在することを示した文献も存在する（非特許文献３、４）。

ただし、このような疾患群のサブタイプを分類する手法を実用化するためには、当該疾患群について、大規模データが必要となる。ところが、大規模に脳画像データを収集することは、健常者についても、そして、特に、患者に対しては容易ではない。

そのため、大規模なデータ収集をするために、複数のサイトでの計測を実施すると、各計測サイトでの計測データのサイト間差が問題となる。上述した非特許文献４でも、多施設の大量の計測データに対するクラスタリングの「汎化（generalization）」は、将来の課題であると言及されている。

たとえば、上述した非特許文献３では、うつ患者が、４つのサブタイプに層別化され、ＴＭＳ(transcranial magnetic stimulation)に対する治療反応性の違いがあることなどを指摘された。ところが、脳機能結合指標を発見する過程で、うつ症状データを２度使用しており過学習のために、うつ症状との関連性に統計学的有意性が確認できず層別化の安定性も悪いことが、別の文献で指摘されている（非特許文献５）。

したがって、たとえば、うつ病に関しては、独立した検証データにおける層別化の精度確認は未実施というのが現状である。

一方で、たとえば、複数の計測サイトで、ＭＲＩ計測をした場合の計測データのサイト間差を評価するために、複数のサイトへ多数の参加者が移動しつつ測定を受けるという、いわゆる「トラベリングサブジェクト（旅行被験者）」を採用することにより、安静時の機能的なコネクティビティに対する測定バイアスの効果を調査するという試みも実施されている（非特許文献６、非特許文献７）。

いずれにしても、ｆＭＲＩデータから被験者の属性の分類を行う場合には、機械学習では過学習の問題を避けるために、一人の被験者を除いて検証用として用いる交差検証法：leave-one-subject-out cross validationやデータを１０分割し、１０分の９で学習し、残りの１０分の１で検証を行う10-fold cross validationを用いて分類器の評価を行うことが多い。しかし、単一施設から得られた少数のサンプルに対して機械学習を適用すると、予測のインフレーションを起こす危険が、精神医学分野でも近年認識されるようになってきている。

少数のデータに対する機械学習では、学習用データにおける特定の施設のｆＭＲＩ装置や測定方法、実験者、参加者群などに存在する特定の傾向、あるいはノイズに対して過学習してしまう可能性が高い。

例えば、脳の解剖画像から自閉スペクトラム症を判別する分類器は、開発に使われた英国の学習用データには感度も特異度も９割以上の高性能を示すが、日本人のデータでは５割になってしまうことが報告されている例もある。このことから、学習用データとは全く異なる施設と被験者群からなる独立検証コホートで検証していない分類器は、科学的にも実用的にも意義が少ないといえる。

上述したような計測サイト間のサイト間差を補償するための「ハーモナイゼーション法」については、本件の出願人による報告もある（非特許文献８）。

特許文献１から特許文献７、及び非特許文献１から非特許文献８の各々の記載の全体をここに参照により援用する。

再表２０１８／１４７１９３号公報（国際公開WO2018/147193）特開２０１９－１９８３７６号公報特開２０２０－０２４１３９号公報特表２０１９－５１６９５０号公報（国際公開WO2017/162773）再表２００５／０２５４２１号公報（国際公開WO2005/025421）特開２０１５－６２８１７号公報特開２０１７－１９６５２３号公報

国立研究開発法人日本医療研究開発機構平成30年12月10日プレスリリース「AIを搭載した内視鏡診断支援プログラムが承認―医師の診断補助に活用へ―」 https://www.amed.go.jp/news/release_20181210.html Glasser MF, et al. The Human Connectome Project's neuroimaging approach. Nat Neurosci 19, 1175-1187 (2016). Andrew T Drysdale, Logan Grosenick, Jonathan Downar, Katharine Dunlop, Farrokh Mansouri, Yue Meng1, Robert N Fetcho, Benjamin Zebley, Desmond J Oathes, Amit Etkin, Alan F Schatzberg, Keith Sudheimer, Jennifer Keller, Helen S Mayberg, Faith M Gunning, George S Alexopoulos, Michael D Fox, Alvaro Pascual-Leone, Henning U Voss, BJ Casey, Marc J Dubin & Conor Liston, "Resting-state connectivity biomarkers define neurophysiological subtypes of depression", nature medicine, VOLUME 23, NUMBER 1, JANUARY 2017 Tomoki Tokuda, Junichiro Yoshimoto,, Yu Shimizu, Go Okada, Masahiro Takamura, Yasumasa Okamoto, Shigeto Yamawaki, Kenji Doya, "Identification of depression subtypes and relevant brain regions using a data-driven approach", SCIENTIFIC REPORTS | (2018) 8:14082 | DOI:10.1038/s41598-018-32521-z Richard Dinga, Lianne Schmaal, Brenda W.J.H.Penninx, Marie Josevan Tol, Dick J.Veltman, Laura van Velzen, Maarten Mennes, Nic J.A.van der Wee, Andre F.Marquand, "Evaluating the evidence for biotypes of depression: Methodological replication and extension of Drysdale et al. (2017)", NeuroImage: Clinical 22 (2019) 101796 Noble S, et al. Multisite reliability of MR-based functional connectivity. Neuroimage 146, 959-970 (2017). Pearlson G. Multisite collaborations and large databases in psychiatric neuroimaging advantages, problems, and challenges. Schizophr Bull 35, 1-2 (2009). Ayumu Yamashita, Noriaki Yahata, Takashi Itahashi, Giuseppe Lisi, Takashi Yamada, Naho Ichikawa, Masahiro Takamura, Yujiro Yoshihara, Akira Kunimatsu, Naohiro Okada, Hirotaka Yamagata, Koji Matsuo, Ryuichiro Hashimoto, Go Okada, Yuki Sakai, Jun Morimoto, Jin Narumoto, Yasuhiro Shimada, Kiyoto Kasai, Nobumasa Kato, Hidehiko Takahashi, Yasumasa Okamoto, Saori C Tanaka, Mitsuo Kawato, Okito Yamashita, and Hiroshi Imamizu，"Harmonization of resting-state functional MRI data across multiple imaging sites via the separation of site differences into sampling bias and measurement bias."，PLOS Biology. DOI: 10.1371/journal.pbio.3000042, http://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3000042

以上のように、機能的磁気共鳴画像法などの脳機能画像法による脳活動の解析の神経・精神疾患の治療への応用を考えた場合、たとえば、上述したようなバイオマーカーとして、脳機能画像法による脳活動の解析は、非侵襲的な機能マーカーとして、診断法の開発、根本治療を実現するための創薬に向けた標的分子の探索・同定などへの応用も期待される。

たとえば、これまで、精神疾患に対しては、遺伝子を用いた実用的なバイオマーカーの完成には至っておらず、それゆえに、薬物の効果判定などが困難であるため、治療薬の開発も困難であった。

脳活動の計測データに基づいて、機械学習により診断マーカーとしての判別器（識別器）や層別化マーカとしての分類器を生成し、これをバイオマーカーとして実用化させるためには、１つの施設において計測された脳活動に対して機械学習により生成されたバイオマーカーの予測精度の向上が必要である。また、このようにして生成されたバイオマーカーが、他の施設において計測された脳活動に対しても、汎化可能であることが必要である。

すなわち、脳活動の計測データに基づく、機械学習により疾患を識別するための判別器や疾患のサブタイプへの分類器を構築する際には、主な２つの課題が存在する。

第１の課題は、小標本サイズの問題である。

被験者の人数であるデータ量Ｎが、計測された脳活動計測データの次元Ｍよりはるかに小さいので、判別器のパラメーターは、容易に訓練データに過剰適合することになる。

この過剰適合のために、構築された判別器は、新しくサンプリングされたテスト・データについては、非常に貧弱なパフォーマンスを示すことになる。これは、これらのテスト・データは判別器を訓練するのに使用されていないものだからである。

したがって、判別器の望ましい汎化に対して本質的な特徴量だけを識別し利用するために特徴量選択や次元削減を適切に導入する必要がある。

第２の課題は、構築された判別器が、訓練データが集められたサイトとは異なる撮像サイトで走査されたＭＲＩデータに対しても、よいパフォーマンスを維持する場合のみ、判別器が臨床的に有用で、科学的に信頼できるということである。

これは撮像サイトにわたる、いわゆる汎化能力である。

精神疾患と関連する大規模な脳画像データを集める際には、脳画像は、１つのサイトで取得できるデータの量には制限があるため、複数のサイトから取得する必要がある。
しかしながら、複数サイトでの脳画像データを取得する場合、サイト間差は最も大きな障壁となる。

すなわち、臨床応用では、異なるサイト撮像されたデータに、特定のサイトで取得したデータを使用して訓練された判別器を汎化することができないことはしばしば観測されている。

このため、上述したヒューマンコネクトームプロジェクトでは、これまでのところ、単一サイトで、単一のスキャナを用いた計測を行うこととされている。

この発明は、上記のような問題点を解決するためになされたものであって、その目的は、複数施設から取得された脳計測データに基づいて、所定の属性を有する被験者のクラスタリングを実行する脳機能結合相関値のクラスタリング装置、脳機能結合相関値のクラスタリングシステム、脳機能結合相関値のクラスタリング方法及び脳機能結合相関値のクラスタリング分類器モデルを提供することである。

この発明の他の目的は、脳活動計測による分類器マーカーを実現するための脳機能結合相関値の分類器プログラムおよび脳活動マーカー分類システムを提供することである。

この発明の１つの局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリングを実行するための脳機能結合相関値のクラスタリング装置であって、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者について、脳活動の計測値に基づいて、クラスタリングの処理を実行するための計算処理システムを備え、計算処理システムは、記憶装置と演算装置とを含み、演算装置は、ｉ）各複数の被験者について、所定の複数の脳領野ペア間の脳活動の時間相関をそれぞれ表す複数の脳機能結合相関値に基づく特徴量を記憶装置に格納し、ｉｉ）記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するように構成され、演算装置は、識別器モデルを生成する機械学習において、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成し、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択し、演算装置は、さらに、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングしてクラスタの分類器を生成する。

好ましくは、脳機能結合相関値のクラスタリング装置は、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受け、計算処理システムは、各複数の被験者について複数の脳機能結合相関値について、計測サイトの測定バイアスを除去するように補正することで、補正された調整値を特徴量として記憶装置に格納するハーモナイゼーション算出手段を含む。

好ましくは、機械学習による識別器の生成の処理は、複数の学習用サブサンプルに対してそれぞれ複数の識別器サブモデルを生成し、複数の識別器サブモデルを統合して識別器モデルを生成するアンサンブル学習である。

好ましくは、属性は、所定の精神疾患であるとの診断結果のラベルで表現され、クラスタリングは、データ駆動による機械学習により、第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理である。

好ましくは、演算装置は、機械学習による識別器の生成において、ｉ）調整値を、機械学習用のトレーニングデータセットと検証用のテストデータセットとに分割し、ｉｉ）トレーニングデータセットに対して、アンダーサンプリングおよびサブサンプリングを所定数だけ実行して、所定数の学習用サブサンプルを生成し、ｉｉｉ）学習用サブサンプルごとに識別器サブモデルを生成し、ｉｖ）識別器サブモデルの出力を統合して、属性の有無に対する識別器モデルを生成する。

好ましくは、機械学習による識別器の生成の処理は、外側の交差検証と内側の交差検証を有する入れ子構造の交差検証であり、演算装置は、入れ子構造の交差検証の処理において、ｉ）外側の交差検証をＫ分割交差検証として、調整値を、機械学習用のトレーニングデータセットと検証用のテストデータセットとに分割し、ｉｉ）トレーニングデータセットに対して、アンダーサンプリングおよびサブサンプリングを所定数だけ実行して、所定数の学習用サブサンプルを生成し、ｉｉｉ）Ｋ分割交差検証の各ループにおいて、内側の交差検証で、ハイパーパラメータを調整して、学習用サブサンプルごとに識別器サブモデルを生成し、ｉｖ）識別器サブモデルに基づいて、属性の有無に対する識別器モデルを生成する。

好ましくは、機械学習による識別器の生成の処理は、特徴量選択を伴う機械学習法であり、クラスタリングのための特徴量の選択においては、和集合に属する特徴量の重要度を、識別器サブモデルを生成する際に特徴量が選択される頻度のランキングにより決定する。

好ましくは、機械学習による識別器の生成の処理は、ランダムフォレスト法であり、クラスタリングのための特徴量の選択においては、和集合に属する特徴量の重要度は、各特徴量についてランダムフォレスト法においてジニ不純度に基づいて算出される重要度である。

好ましくは、機械学習による識別器の生成の処理は、Ｌ２正則化による機械学習法であり、クラスタリングのための特徴量の選択においては、和集合に属する特徴量の重要度を、Ｌ２正則化により算出される識別器サブモデルにおける特徴量の重みに基づくランキングにより決定する。

好ましくは、記憶装置は、複数の計測サイトにおいて共通に計測対象となる複数の旅行被験者に対して、各旅行被験者の予め定められた複数の脳領域について脳活動を計測した結果を予め格納しており、演算装置は、複数の脳領域のペアについての脳活動の時間相関を表す脳機能結合行列の所定の要素を各旅行被験者について算出し、一般化線形混合モデル法を用いることで、機能結合行列の所定の要素ごとに、複数の計測サイトおよび複数の旅行被験者にわたる当該要素の平均に対する各計測サイトにおける固定効果として、測定バイアスを算出する。

好ましくは、演算装置は、複数の計測サイト以外の計測サイトで計測された対象者についての計測データに基づいて、サブタイプへの分類処理を実行する。

この発明の他の局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリングを実行するための脳機能結合相関値のクラスタリングシステムであって、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者についての脳活動を時系列で計測するために、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置と、複数の被験者について、脳活動の計測値に基づいて、クラスタリングの処理を実行するための計算処理システムとを備え、計算処理システムは、記憶装置と演算装置とを含み、演算装置は、ｉ）各複数の被験者について、所定の複数の脳領野ペア間の脳活動の時間相関をそれぞれ表す複数の脳機能結合相関値に基づく特徴量を記憶装置に格納し、ｉｉ）記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するように構成され、演算装置は、識別器モデルを生成する機械学習において、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成し、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択し、演算装置は、さらに、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングしてクラスタの分類器を生成する。

好ましくは、計算処理システムは、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受け、各複数の被験者について複数の脳機能結合相関値について、計測サイトの測定バイアスを除去するように補正することで、補正された調整値を特徴量として記憶装置に格納するハーモナイゼーション算出手段を含む。

この発明のさらに他の局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリング処理を計算処理システムが実行するための脳機能結合相関値のクラスタリング方法であって、計算処理システムは、記憶装置と演算装置とを含み、演算装置が、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を記憶装置に格納するステップと、演算装置が、記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを備え、識別器モデルを生成する機械学習を実行するステップは、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、さらに、演算装置が、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを備える。

この発明のさらに他の局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳機能結合相関値の分類器プログラムであって、分類器プログラムは、各クラスタの確率分布のモデルに基づき、コンピュータが、入力データが最大事後確率を有するクラスタへ分類する分類機能を有し、計算処理システムは、記憶装置と演算装置とを含み、計算処理システムは、クラスタリング処理に基づく分類器プログラムの生成処理において、演算装置が、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を記憶装置に格納するステップと、演算装置が、記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、識別器モデルを生成する機械学習を実行するステップは、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、さらに、演算装置が、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを実行する。

好ましくは、計算処理システムは、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受信するステップと、各複数の被験者について複数の脳機能結合相関値について、計測サイトの測定バイアスを除去するように補正することで、補正された調整値を特徴量として記憶装置に格納するハーモナイゼーションを実行するステップとを実行する。

この発明のさらに他の局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳活動マーカー分類システムであって、脳活動マーカー分類システムは、各クラスタの確率分布のモデルに基づき、コンピュータが、入力データが最大事後確率を有するクラスタへ分類する分類機能を有し、計算処理システムは、記憶装置と演算装置とを含み、計算処理システムは、クラスタリング処理に基づく脳活動マーカー分類システムの生成処理において、演算装置が、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を記憶装置に格納するステップと、演算装置が、記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、識別器モデルを生成する機械学習を実行するステップは、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、さらに、演算装置が、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを実行する。

好ましくは、計算処理システムは、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受信するステップと、各複数の被験者について脳活動の時間相関を表す複数の脳機能結合相関値について、計測サイトの測定バイアスを除去するように補正することで、補正された調整値を特徴量として記憶装置に格納するハーモナイゼーションを実行するステップとを実行する。

この発明のさらに他の局面に従うと、対象者の脳活動の計測結果に基づいて、対象者において、少なくとも１つの所定の属性を有する対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳機能結合相関値のクラスタリング分類器モデルであって、クラスタリング分類器モデルは、学習データに含まれる対象者を特徴づける特徴量群を分割したビューごとに、各ビューに含まれる特徴量の情報と、各ビューの対象者のクラスタごとの確率密度関数を特定する情報とに基づいて、入力データに対して算出される確率密度関数の値に応じて、入力データが最大事後確率を有するクラスタへ分類する機能を有し、計算処理システムは、記憶装置と演算装置とを含み、計算処理システムは、クラスタリング処理に基づくクラスタリング分類器モデルの生成処理において、演算装置が、所定の属性を有する第１の被験者群と所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を記憶装置に格納するステップと、演算装置が、記憶装置に格納された特徴量に基づいて、属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、識別器モデルを生成する機械学習を実行するステップは、第１の被験者群と第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、さらに、演算装置が、選択されたクラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、第１の被験者群をクラスタリングして、特徴量をビューへ分割し、各ビューにおける対象者のクラスタごとの確率密度関数を生成するステップを実行する。

この発明によれば、複数の施設で計測される脳の活動の計測データに対して、各施設における測定バイアスを調整および補正することが可能である。これにより、複数の施設での計測データに基づく、脳機能結合相関値を調整してクラスタリングすることが可能となる。

また、この発明によれば、複数の施設で計測される脳の活動の計測データをハーモナイズして、疾患群のサブタイプを客観的に判断することが可能な脳機能結合相関値のクラスタリング装置、脳機能結合相関値のクラスタリングシステム、脳機能結合相関値のクラスタリング方法、脳機能結合相関値の分類器プログラムおよび脳活動マーカー分類システムを実現することができる。

複数の計測サイトに設置されたＭＲＩ計測システムにより計測されたデータについて、ハーモナイゼーション処理を説明するための概念図である。被験者の脳の関心領域（ＲＯＩ:Region of Interest）について、安静時の機能結合の相関を示す相関行列を抽出する手続きを示す概念図であり、ｒｓｆＭＲＩの時系列データを示す図である。被験者の脳の関心領域について、安静時の機能結合の相関を示す相関行列を抽出する手続きを示す概念図であり、相関行列を示す図である。「測定パラメータ」の内容の例を示す概念図である。「被験者属性データ」の内容の例を示す概念図である。各計測サイトに設置されるＭＲＩ装置１００．ｉ（１≦ｉ≦Ｎｓ）の全体構成を示す模式図である。データ処理部３２のハードウェアブロック図である。相関行列から、診断マーカーとなる判別器を生成する過程およびクラスタリング処理について説明する概念図である。計算処理システム３００の構成を説明するための機能ブロック図である。計算処理システム３００の構成を説明するための機能ブロック図である。アンサンブル学習による疾患識別器を生成するための機械学習の手続きを説明するためのフローチャートである。学習用のデータセット（データセット１）のデモグラフィック特性を示す図である。独立検証のデータセット（データセット２）のデモグラフィック特性を示す図である。全撮像サイトについて、学習用データセットに対するＭＤＤの予測性能（出力の確率分布）を示す図である。各撮像サイトについて、学習用データセットに対するＭＤＤの予測性能（識別器の出力の確率分布）を示す図である。独立した検証データセットにおけるＭＤＤの識別器の出力の確率分布を示す図である。各撮像サイトについて、独立検証データセットに対するＭＤＤの識別器の出力の確率分布を示す図である。特徴量を選択して、教師なし学習によりクラスタリングを行う処理を説明するためのフローチャートである。「特徴量選択を伴う学習処理」により、複数個（たとえば、Ｎch個）の特徴量が存在する場合に、特徴量の選択が実施される概念を示す図である。特徴量選択を伴う学習処理により１つの識別器を生成する際に最終的に選択される特徴量を示す概念図である。アンダーサンプリングおよびサブサンプリングの処理を複数回実施して、識別器を生成される際に、特徴量が選択される様子を示す概念図である。特徴量に依存して複数のクラスタの分け方が存在する場合を説明するための概念図である。複数の特徴量により、複数の対象が特徴づけられる場合のクラスタリングの概念を説明するための概念図である。複数の特徴量により、複数の対象が特徴づけられる場合のクラスタリングの概念を説明するための概念図である。多重クラスタリングを説明するための概念図である。多重共クラスタリングを説明するための概念図である。「多重共クラスタリング」において、１つのビューの中に、異なる種類の確率分布の確率モデルが想定されている場合を示す概念図である。多重共クラスタリングの学習方法の概要を説明するためのフローチャートである。多重共クラスタリングの学習方法におけるベイズ推定のグラフ表現を示す図である。２つに分割されたデータセットのうちデータセット１を示す図である。２つに分割されたデータセットのうちデータセット２を示す図である。各データセットにおいてクラスタリングを実施する概念を説明する概念図である。被験者データに対する多重共クラスタリングの例を示す概念図である。データセット１およびデータセット２に対して、実際に多重共クラスタリング処理を実施した結果を示す図である。データセット１とデータセット２において、それぞれのビューに割り当てられた脳機能結合（ＦＣ）の個数を示す表である。クラスタリングの類似度（層別化の汎化性能）の評価手法を説明するための概念図である。ＡＲＩについて説明するための概念図である。ＡＲＩについて説明するための概念図である。データセット１とデータセット２について、それぞれのビューに対するＡＲＩを算出した結果を表形式で示す図である。クラスタリング１とクラスタリング１´の間、および、クラスタリング２とクラスタリング２´の間の類似度の評価結果を示す図であり、図３３Ａに対応するパーミュテーションテスト（Permutation Test）の結果を示す図である。クラスタリング１とクラスタリング１´において、それぞれビュー１の各クラスタに割り当てられる被験者数の分布を示す表である。サイト間を移動しながら計測を受ける移動被験者によるサイト間差の評価手法を説明するための概念図である。被験者ａのｂ番目の機能的なコネクティビティの表現を説明するための概念図である。ハーモナイゼーションのために、測定バイアスを算出する処理を説明するためのフローチャートである。新たな計測サイトが追加された場合のハーモナイゼーション処理のための測定バイアスの算出処理を説明するための概念図である。ハーモナイゼーション処理に使用する、多疾患データベースのデータセットを示す概念図である。ハーモナイゼーション処理に使用する、トラベリングサブジェクトのデータセットを示す概念図である。ＳＲＰＢＳの多疾患データセットの内容を示す図である。ＳＲＰＢＳの多疾患データセットの内容を示す図である。ＳＲＰＢＳの多疾患データセットの内容を示す図である。ＳＲＰＢＳの多疾患データセットの内容を示す図である。ＳＲＰＢＳの多疾患データセットの内容を示す図である。各計測サイトにおける撮像プロトコルを示す図である。各計測サイトにおける撮像プロトコルを示す図である。各計測サイトにおける撮像プロトコルを示す図である。各計測サイトにおける撮像プロトコルを示す図である。各計測サイトにおける撮像プロトコルを示す図である。主成分分析によるサイト間差および疾患効果の視覚化を示すグラフである。階層的クラスタリング分析による樹状図である。各要因についての寄与サイズを示す図である。ハーモナイゼーションプロセスの影響を視覚化した図であり、図４２と対比される図である。データ収集、推定処理および対象者の脳活動計測を、分散して処理する場合の一例を示す機能ブロック図である。トラベリングサブジェクトの巡回方式を示す図である。クラスタリング分類器の構成を示す概念図である。

以下、本発明の「脳機能結合相関値のクラスタリング装置」、「脳機能結合相関値のクラスタリング方法」等を説明するために、複数の脳活動計測装置からなる計測システムで計測した被験者（精神疾患の患者を含む）の脳機能結合画像データに対する人工知能技術による「クラスタリング」を例として説明することとする。

そこで、以下、本発明の実施の形態の計測システム、より具体的にはＭＲＩ計測システムの構成について、図に従って説明する。なお、以下の実施の形態において、同じ符号を付した構成要素および処理工程は、同一または相当するものであり、必要でない場合は、その説明は繰り返さない。

また、本実施の形態では、複数の施設に設置された「脳活動計測装置」、より具体的には、「ＭＲＩ装置」によって、脳の複数の領野間の脳活動を時系列で計測し、これらの領野間の時間相関（「脳機能結合」と呼ぶ）のパターンに基づいて、複数施設に汎化可能に、特定の疾患の被験者を、さらに複数のグループ（サブグループ）に分類するものとして本発明を説明する。

なお、これも特に限定されないが、「特定の疾患」については、「大うつ病」を例として説明することとする。ただし、以下の説明にも示されるように、本発明は、被験者の「脳機能結合相関値」をデータ駆動による分類する技術に関するものであって、被験者の疾患は、「大うつ病」に限定されず、他の疾患であってもよい。また、被験者の「脳機能結合相関値」のパターンにより分類される被験者の属性であれば、必ずしも、疾患である必要もなく、他の属性であってもよい。

そして、このような「ＭＲＩ計測システム」では、複数の「ＭＲＩ装置」が、複数の異なる施設に設置されており、後述するように、これらの計測施設（計測サイト）間での計測に対するサイト間差を、測定機器を原因とする測定バイアスと、計測サイトにおける被験者のポピュレーションによる差（標本バイアス）とを独立に評価する。その上で、各計測サイトにおける計測値について、測定バイアスの効果を除いてサイト間差を補正する処理を実施することで、計測サイト間での測定結果に対する調和処理（ハーモナイゼーション）を実現する。そして、ハーモナイゼーション後の脳機能結合値について、特定疾患の診断ラベルを教師データとするアンサンブル学習を用いて、「特徴量選択」を行った後に、教師なし学習によるクラスタリングを実行することで、被験者属性（たとえば、精神疾患のサブタイプ）の分類を実行するものとして、説明する。

［実施の形態１］
図１は、複数の計測サイトに設置されたＭＲＩ計測システムにより計測されたデータについて、クラスタリング（層別化）処理を説明するための概念図である。

図１を参照して、計測サイトＭＳ．１～ＭＳ．Ｎｓ（Ｎｓ：サイト数）には、それぞれ、ＭＲＩ装置１００．１～１００．Ｎｓが設置されているものとする。

また、計測サイトＭＳ．１～ＭＳ．Ｎｓにおいては、それぞれ、被験者群ＰＡ．１～ＰＡ．Ｎｓの計測が行われる。被験者群ＰＡ．１～ＰＡ．Ｎｓの各々には、少なくとも２つ以上に分類される群、たとえば、患者群と健常者群とが含まれるものとする。また、患者群としては、特に限定されないが、たとえば、精神疾患の患者、より特定的には、「大うつ病の患者」の群が相当するものとする。

そして、各計測サイトにおいては、原則として、ＭＲＩ装置の仕様上、可能な範囲で、統一した計測プロトコルで被験者の計測が実施されるものとする。

ここで、特に限定されないが、計測プロトコルとしては、たとえば、以下のような内容が規定されているものとする。

１）頭部のスキャンを実行する方向
たとえば、頭部の後側（posterior：以下、“Ｐ”と略記する）から前側（anterior：以下、“Ａ”と略記する）に向かう方向（以下、「Ｐ→Ａ方向」と呼ぶ）と、逆の方向、すなわち、前側から後側に向かう方向（以下、「Ａ→Ｐ方向」と呼ぶ）にスキャンを行うのかを規定する必要がある。状況によっては、両方の場合のスキャンを行うことを規定することもあり得る。

ＭＲＩ装置によって、デフォルトの方向が違っていたり、あるいは、両方を任意に設定できない場合もあり得る。

スキャンの方向は、たとえば、画像としての「ひずみ方」を規定する可能性があり、プロトコルとして、条件を設定する。

２）脳構造画像の撮像条件
いわゆるスピンエコー法により、「Ｔ１強調画像」や「Ｔ２強調画像」のいずれか、または両方を撮像する条件を設定する。

３）脳機能画像の撮像条件
ｆＭＲＩ（functional Magnetic Resonance Imaging）法により、「安静時」の被験者の脳機能画像を撮像する条件を設定する。

４）拡散強調画像の撮像条件
拡散強調画像（ＤＷＩ：diffusion (weighted) image）を撮像するか否か、そしてその条件について設定する。

拡散強調画像とは、ＭＲＩ撮像のシーケンスの一種で、水分子の拡散運動を画像化したものである。通常用いられるスピンエコー法のパルス系列では拡散による信号の減衰は無視できるが、大きな傾斜磁場が長時間にわたって印加されると、その間の各磁化ベクトルの移動によって生じる位相のずれが無視できなくなり、拡散が活発な領域ほど低信号として現れる、ことを利用するものである。

５）ＥＰＩ歪みを画像処理で補正するための撮像
たとえば、ＥＰＩ歪みを画像処理で補正するための１つの方法として、「フィールドマップ法」が知られており、空間歪みの補正に対する撮像の条件を設定する。

フィールドマップ法は、多重エコー時間によってＥＰＩ画像を収集し、これらのＥＰＩ画像に基づいてＥＰＩ歪みの量を算出する。フィールドマップ法を適用して、新規の画像に含まれるＥＰＩ歪みを補正することができる。異なるエコー時間による同一の解剖学的構造の一組の画像を前提として、ＥＰＩ歪みを計算して、画像の歪みを補正することが可能である。

たとえば、以下の公知文献１に、「フィールドマップ法」が開示されている。公知文献１の記載の全体をここに参照により援用する。

公知文献１：特開2015-112474号明細書
なお、計測プロトコルには、上記の条件の中から適宜、必要となるシーケンス部分を抜き出してもよいし、必要に応じて、他のシーケンスやその条件を追加してもよい。

再び、図１を参照して、各計測サイトＭＳ．１～ＭＳ．Ｎｓにおいて、計測を行う被験者として採用することを、「被験者をサンプリングする」と呼び、各計測サイトでのサンプリングの偏りにより生じる、計測値のサイト間差の原因を、「標本バイアス」と呼ぶ。

たとえば、上述した例では、従来の診断基準で、「大うつ病」と診断される患者には、実は、いくつかのサブタイプが含まれることが知られている。

典型的には、「メランコリー型」、「非定型」、「季節型」、「産後」うつ病などがある。また、「２種類以上の作用機序の異なる抗うつ薬による十分な治療により十分に改善せず、中等症以上の症状が続く場合」は、「治療抵抗性うつ病」と呼ばれ、うつ病の１０～２０％と推定されるとの報告もある。すなわち、一般に、「大うつ病」と診断される患者群は、決して、均質なものではないことが知られている。ただし、このようなサブタイプを、客観的な計測データに基づいて、分類する手法は、これまでは、必ずしも実用化に至っていない。

各計測サイトでは、当該計測サイトの病院に来院する患者の地域性に由来する性向の偏りや、当該病院での診断の傾向等々のさまざまの要因により、一口に、「大うつ病」の患者といったとしても、その患者群に含まれるサブタイプの分布が均一であるとは、必ずしも言えない。この結果、各計測サイトの患者群には、そのサブタイプの分布に、偏りが生じることが、むしろ一般的であり、その結果、上述した「標本バイアス」が生じると考えられる。

また、「健常群」と呼ばれる被験者の群であっても、その中には、複数のサブタイプが存在することが一般であり、その点では、「健常群」においても、「標本バイアス」が存在することになる。

また、ＭＲＩ装置１００．１～１００．Ｎｓも、各計測サイトで決して同一計測特性を有するＭＲＩ装置が使用されているとは言えない。

たとえば、ＭＲＩ装置の製造メーカ、ＭＲＩ装置の型番、ＭＲＩ装置の静磁場強度、ＭＲＩ装置における（送）受信コイルのコイル数（チャンネル数）、等々、ＭＲＩ装置の条件、ＭＲＩ装置の計測条件などの測定条件に応じて、計測データのサイト間差異が生じ得る。このような測定条件により生じるサイト間差異のことを「測定バイアス」と呼ぶ。

仮に、同一のＭＲＩ装置の製造メーカの同一の型番のＭＲＩ装置であっても、装置の固有な個性により、完全に同一の計測特性を実現するとも、必ずしも言えない。

ここで、（送）受信コイルは、計測される信号のＳＮ比を向上させることを目的として、「マルチアレイコイル」が採用されることが一般的である。「受信コイルのコイル数」とは、マルチアレイコイルを構成する「エレメントコイル」の個数を指す。個々のエレメントコイルの感度を高め，この出力を束ねることで受信感度の向上を図るものである。

そして、特に限定されないが、後に説明するようなハーモナイゼーション手法により、本実施の形態では、「標本バイアス」と「測定バイアス」とを独立に評価することを可能とする。

再び、図１を参照して、各計測サイトＭＳ．１～ＭＳ．Ｎｓのそれぞれからの計測関連データＤＡ１００．１～ＤＡ１００．Ｎｓは、データセンター２００内の記憶装置２１０に集積され、格納される。

ここで、「計測関連データ」には、各計測サイトでの「測定パラメータ」と、各計測サイトにおいて計測された「患者群データ」および「健常群データ」とが含まれる。

さらに、「患者群データ」および「健常群データ」には、それぞれ、各々の被験者に対応して「患者のＭＲＩ計測データ」および「健常者のＭＲＩ計測データ」とが含まれる。

以下では、このような「計測関連データ」について説明する。

図２Ａおよび図２Ｂは、被験者の脳の関心領域について、安静時の機能結合の相関を示す相関行列を抽出する手続きを示す概念図である。図２ＡはｒｓｆＭＲＩの時系列データを示す図であり、図２Ｂはその相関行列を示す図である。

ここで、図１において、「患者群データ」および「健常群データ」における、「患者のＭＲＩ計測データ」および「健常者のＭＲＩ計測データ」には、少なくとも、以下のようなデータが含まれる。
ｉ）相関行列のデータを算出するための時系列の「脳機能画像データ」、および／または、相関行列のデータ自身
すなわち、図１における計算処理システム３００が、記憶装置２１０内に格納されたデータに基づいて、後述するような脳活動バイオマーカを算出する際の基礎となるデータとして、これらのデータが使用される。

ここで、相関行列のデータについては、時系列の「脳機能画像データ」に基づいて、各計測サイトで計算された後に、記憶装置２１０内に格納され、計算処理システム３００が、この記憶装置２１０内の相関行列のデータに基づいて、脳活動バイオマーカを算出する、という構成とすることができる。

あるいは、時系列の「脳機能画像データ」が記憶装置２１０内に格納され、計算処理システム３００が、この記憶装置２１０内の「脳機能画像データ」に基づいて、相関行例のデータを算出し、さらに、脳活動バイオマーカを算出する、という構成とすることも可能である。

したがって、「患者のＭＲＩ計測データ」および「健常者のＭＲＩ計測データ」の各々については、少なくとも、相関行列のデータを算出するための時系列の「脳機能画像データ」、および、相関行列のデータ自身のいずれか一方が含まれることになる。

ｉｉ）被験者の構造画像データ、および、拡散強調画像データ
なお、特に限定されないが、ＥＰＩ歪みを画像処理で補正する処理については、各計測サイトで演算処理がされた後に、記憶装置２１０内にデータが格納される構成とすることができる。

また、特に限定されないが、個人情報保護の観点から、記憶装置２１０内へデータが格納される前に、各計測サイトにおいて、匿名加工処理を実行する構成とすることが可能である。ただし、匿名加工処理については、計算処理システム３００を運営する主体が、法規制上、個人情報を扱うことが認められている場合などは、計算処理システム３００において匿名加工処理を実行する構成としてもよい。

図２Ａおよび図２Ｂに戻って、図２Ａに示すように、リアルタイムで測定した安静時のｆＭＲＩのｎ個（ｎ：自然数）の時刻分のｆＭＲＩデータから、各関心領域の平均的な「活動度」を算出し、図２Ｂに示すように、脳領域間（関心領域間）の機能的なコネクティビティ（「活動度の相関値」）の相関行列を算出する。

（脳領域の区画化（パーセレーション：Parcellation））
機能的なコネクティビティは、各参加者に対して２つの脳領域間の安静時の機能ＭＲＩの血中酸素濃度に依存する（ＢＯＬＤ）信号の時間的な相関性として計算される。

ここで、関心領域としては、上述したようにＮr領域を考えているので、相関行列における独立な非対角成分は、対称性を考慮すると、
Ｎr×（Ｎr－１）／２（個）
ということになる。

関心領域の設定の方法としては、以下のような方法が想定される。

方法１）“解剖学的な脳領野に基づいて、関心領域を定義する。”
ここでは、脳活動バイオマーカーのために、たとえば、関心領域として１４０の領域を採用する。

つまり、この方法では、ＲＯＩについては、Brain Sulci Atlas （BAL）に含まれる１３７のＲＯＩに加え、 Automated Anatomical Labeling Atlas の小脳（左右）および虫部のＲＯＩを用いる。これら合計１４０のＲＯＩ間の機能的結合ＦＣを特徴量として用いる。

ここで、Brain Sulci Atlas （BAL）およびAutomated Anatomical Labeling Atlasについては、以下に開示がある。以下の公知文献２及び公知文献３の各々の記載の全体をここに参照により援用する。

公知文献２：Perrot et al., Med Image Anal, 15(4), 2011
公知文献３：Tzourio-Mazoyer et al., Neuroimage, 15(1), 2002
このような関心領域としては、たとえば、以下のような領域である。

背内側前頭前皮質（DMPFC）
前頭前野腹内側部（VMPFC）
前帯状皮質（ACC）
小脳虫部、
左視床、
右下頭頂葉、
右尾状核、
右中後頭葉、
右中帯状皮質
ただし、採用する脳の領野は、このような領域に限定されるものではない。

たとえば、対象とする神経・精神疾患に応じて、選択する領域を変更してもよい。

方法２） “機能的なコネクティビティを、全体の脳をカバーする機能的な脳地図の脳領域に基づいて、定義する。”

ここで、このような機能的な脳地図の脳領域については、以下のような文献にも開示がされており、特に限定されないが、たとえば、２６８のノード（脳領域）から成る、という構成とすることができる。以下の公知文献４から公知文献７の各々の記載の全体をここに参照により援用する。

公知文献４：Noble S, et al. Multisite reliability of MR-based functional connectivity. Neuroimage 146,959-970 (2017).
公知文献５：Finn ES, et al. Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nat Neurosci 18, 1664-1671 (2015).
公知文献６：Rosenberg MD, et al. A neuromarker of sustained attention from whole-brain functional connectivity. Nat Neurosci 19, 165-171 (2016).
公知文献７：Shen X, Tokoglu F, Papademetris X, Constable RT. Groupwise whole-brain parcellation from resting-state fMRI data for network node identification. Neuroimage 82, 403-415 (2013).

方法３）表面ベースの方法
脳領域のパーセレーションについては、ヒューマンコネクトームプロジェクト(ＨＣＰ)スタイルのマルチモダリティなイメージング（ミエリン・タスク・ファンクショナル）を用いることにより、脳を脳溝に沿ってシート状に変換して作成された脳地図に基づく「表面ベースの方法」で、データを分析することも可能である。

このようなパーセレーション法のためには、以下のサイトに開示されるようなツールボックスを使用することができる（ciftifyツールボックスバージョン2.0.2）。https://edickie.github.io/ciftify/#/
このツールボックスは、使用するデータを分析することを、（たとえば、ＨＣＰパイプラインに必要なＴ２強調画像を欠いている場合でも）ＨＣＰ類似の表面ベースのパイプラインの中で可能にするものである。

そして、方法３の解析では、関心領域（ＲＯＩ）として、以下の公知文献８に開示される３７９の表面ベースの区画（皮質の３６０の区画＋皮質下の１９の区画）を使用する。以下の公知文献８の記載の全体をここに参照により援用する。
公知文献８：Glasser, M.F., Coalson, T.S., Robinson, E.C., Hacker, C.D., Harwell, J., Yacoub, E., et al. (2016). A multi-modal parcellation of human cerebral cortex. Nature 536(7615), 171-178. doi: 10.1038/nature18933.
したがって、ＢＯＬＤ信号の時間的変化は、これらの３７９個の関心領域（ＲＯＩ）から抽出される。

さらに、以下の公知文献９に開示されるような解剖学的な自動標識付（ＡＡＬ）およびNeurosynth（http://neurosynth.org/locations/）の使用により、重要なＲＯＩの解剖学的名称、およびＲＯＩを含む内在性の脳ネットワークの名称を特定することができる。以下の公知文献９の記載の全体をここに参照により援用する。

公知文献９：Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., et al. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage 15(1), 273-289. doi: 10.1006/nimg.2001.0978.

方法４）データドリブンに脳の領域を定めた方法
下記の公知文献１０にも開示があるように、事前情報（脳地図）なしに位相のそろったボクセルからネットワークを新規に同定する手法であり、「Canonical ICA」や「辞書学習」などと呼ばれる手法である。以下の公知文献１０の記載の全体をここに参照により援用する。

公知文献１０：Kamalaker Dadi，Mehdi Rahim，Alexandre Abraham，Darya Chyzhyk，Michael Milham，Bertrand Thirion，Gael Varoquaux，“Benchmarking functional connectome-based predictive models for resting-state fMRI”，Preprint submitted to NeuroImage，October 31，2018.

以下では、基本的に「方法３」の表面ベースの脳地図の脳領域に基づいて機能的なコネクティビティを定義する方法を用いるものとして説明する。

また、相関値の計算としては、たとえば、タンジェント法や部分相関法のような機能的なコネクティビティの計測に対するいくつかの候補がある。

しかしながら、以下では、特に限定されないが、ピアソンの相関係数を使用するものとして説明する。

前処理されたＢＯＬＤ信号の可能なノード組の各々の時間経過の間についてのフィッシャーのｚ変換後のピアソンの相関係数が計算され、要素がそれぞれ２つのノード間のコネクティビティの強度を表わす、３７９×３７９の対照的な脳機能結合マトリックスを構築するために使用される。

図３Ａおよび図３Ｂは、「測定パラメータ」および「被験者属性データ」の内容の例をそれぞれ示す概念図である。

「被験者属性データ」は、図１における「患者群データ」または「健常群データ」において、「患者のＭＲＩ計測データ」および「健常者のＭＲＩ計測データ」に、それぞれ関連づけて格納されているものとする。

図３Ａに示すように、計測サイトを識別するためのサイトＩＤと、サイト名と、測定パラメータを識別するための条件ＩＤと、計測装置に関する情報と、計測条件に関する情報とを含む。
「測定パラメータ」は、「計測装置に関する情報」と「計測条件に関する情報」とを含む。
「計測装置に関する情報」には、各計測サイトにおいて被験者の脳活動を計測するためのＭＲＩ装置のメーカ名と、型番と、（送）受信コイル数とを含む。
なお、「計測装置に関する情報」は、これらのものに限られず、たとえば、静磁場強度、シム調整後の磁場の均一性、等々の他の計測装置の性能を表す指標などを含んでいてもよい。

「計測条件に関する情報」には、画像再構成時の位相エンコードの方向（Ｐ→ＡまたはＡ→Ｐ）、画像種類（Ｔ１強調、Ｔ２強調、拡散強調など）、撮像シーケンス（スピンエコーなど）、撮像中の被験者の開眼／閉眼、等々の情報が含まれる。
「計測条件に関する情報」も、これらのものに限られるものではない。

図３Ｂに示されるように、「被験者属性データ」は、被験者の識別ができないように仮名制御した被験者仮ＩＤ，その被験者を計測した際の測定条件を表す条件ＩＤ、被験者の属性情報とを含む。

そして、「被験者の属性情報」としては、被験者の性別、年齢、健常または疾患のいずれかを示すラベル、医師による被験者の診断病名、被験者への投薬履歴、診断履歴などの情報が含まれる。
なお、「被験者の属性情報」は、必要に応じて、たとえば、計測サイトにおいて、匿名加工処理が行われているものとする。

たとえば、年齢や性別などについては、準識別子(同一属性)のデータがｋ件以上になるようにデータを変換するなどの処理を行うことで、個人が特定される確率をｋ分の１以下に低減し特定を困難にする「ｋ-匿名性」を維持するように加工することができる。ここで、「準識別子」とは、「年齢」「性別」「居住地」など単体では特定できないが、組み合わせることで個人の特定が可能になる属性をいう。
また、投薬履歴や診断履歴は、日付のランダム化やシフト化（相対化）など、必要に応じて、匿名化のための処理を行う。

そして、以下では、「患者のＭＲＩ計測データ」および「健常者のＭＲＩ計測データ」で、上述した方法で、各被験者の脳領野間の各々について経時的な活動の相関として算出される機能的結合のことを、各領野間について総称して「機能的なコネクティビティ（Connectivity）」（省略する際には、“ＦＣ”と表記）と呼ぶことにする。脳領野ごとに、機能的なコネクティビティを区別する必要がある場合は、後述するように、下添え字を付加して区別するものとする。

（ＭＲＩ装置の構成）
図４は、各計測サイトに設置されるＭＲＩ装置１００．ｉ（１≦ｉ≦Ｎｓ）の全体構成を示す模式図である。
図４においては、１番目の計測サイトのＭＲＩ装置１００．１を、例示的に、詳しく説明している。他のＭＲＩ装置１００．２～１００．Ｎｓについても、基本的な構成は同様である。

図４に示すように、ＭＲＩ装置１００．１は、被検者２の関心領域に制御された磁場を付与してＲＦ波を照射する磁場印加機構１１と、この被検者２からの応答波（ＮＭＲ信号）を受信してアナログ信号を出力する受信コイル２０と、この被検者２に付与される磁場を制御するとともにＲＦ波の送受信を制御する駆動部２１と、この駆動部２１の制御シーケンスを設定するとともに各種データ信号を処理して画像を生成するデータ処理部３２とを備える。
なお、ここで、被検者２が載置される円筒形状のボアの中心軸をＺ軸にとりＺ軸と直交する水平方向にＸ軸および鉛直方向にＹ軸を定義する。

ＭＲＩ装置１００．１は、このような構成であるので、磁場印加機構１１により印加される静磁場により、被検者２を構成する原子核の核スピンは、磁場方向（Ｚ軸）に配向するとともに、この原子核に固有のラーモア周波数でこの磁場方向を軸とする歳差運動を行う。

そして、このラーモア周波数と同じＲＦパルスを照射すると、原子は共鳴しエネルギーを吸収して励起され、核磁気共鳴現象（ＮＭＲ現象；Nuclear Magnetic Resonance）が生じる。この共鳴の後に、ＲＦパルス照射を停止すると、原子はエネルギーを放出して元の定常状態に戻る緩和過程で、ラーモア周波数と同じ周波数の電磁波（ＮＭＲ信号）を出力する。

この出力されたＮＭＲ信号を被検者２からの応答波として受信コイル２０で受信し、データ処理部３２において、被検者２の関心領域が画像化される。

磁場印加機構１１は、静磁場発生コイル１２と、傾斜磁場発生コイル１４と、ＲＦ照射部１６と、被検者２をボア中に載置する寝台１８とを備える。

被験者２は、寝台１８に、たとえば、仰臥する。被験者２は、特に限定されないが、たとえば、プリズムメガネ４により、Ｚ軸に対して垂直に設置されたディスプレイ６に表示される画面を見ることができる。このディスプレイ６の画像により、必要に応じて、被験者２に視覚刺激を与えることもできる。なお、被験者２への視覚刺激は、被験者２の目前にプロジェクタにより画像が投影される構成であってもよい。
このような視覚刺激は、被験者にニューロフィードバックを行う場合は、フィードバック情報の提示に相当する。

駆動部２１は、静磁場電源２２と、傾斜磁場電源２４と、信号送信部２６と、信号受信部２８と、寝台１８をＺ軸方向の任意位置に移動させる寝台駆動部３０とを備える。

データ処理部３２は、操作者（図示略）から各種操作や情報入力を受け付ける入力部４０と、被検者２の関心領域に関する各種画像および各種情報を画面表示する表示部３８と、表示部３８による表示を制御する表示制御部３４と、各種処理を実行させるプログラム・制御パラメータ・画像データ（構造画像等）およびその他の電子データを記憶する記憶部３６と、駆動部２１を駆動させる制御シーケンスを発生させるなどの各機能部の動作を制御する制御部４２と、駆動部２１との間で各種信号の送受信を実行するインタフェース部４４と、関心領域に由来する一群のＮＭＲ信号からなるデータを収集するデータ収集部４６と、このＮＭＲ信号のデータに基づいて画像を形成する画像処理部４８と、ネットワークとの間で通信を実行するためのネットワークインタフェース５０を備える。

また、データ処理部３２は、専用コンピュータである場合の他、各機能部を動作させる機能を実行する汎用コンピュータであって、記憶部３６にインストールされたプログラムに基づいて、指定された演算やデータ処理や制御シーケンスの発生をさせるものである場合も含まれる。以下では、データ処理部３２は、汎用コンピュータであるものとして説明する。

静磁場発生コイル１２は、Ｚ軸周りに巻回される螺旋コイルに静磁場電源２２から供給される電流を流して誘導磁場を発生させ、ボアにＺ軸方向の静磁場を発生させるものである。このボアに形成される静磁場の均一性の高い領域に被検者２の関心領域を設定することになる。ここで、静磁場発生コイル１２は、より詳しくは、たとえば、４個の空芯コイルから構成され、その組み合わせで内部に均一な磁界を作り、被験者２の体内の所定の原子核、より特定的には水素原子核のスピンに配向性を与える。

傾斜磁場発生コイル１４は、Ｘコイル、ＹコイルおよびＺコイル（図示省略）から構成され、円筒形状を示す静磁場発生コイル１２の内周面に設けられる。
なお、傾斜磁場の均一性を向上させるためにシムコイル（図示せず）が設けられ、「シム調整」が実施される。
これらＸコイル、ＹコイルおよびＺコイルは、それぞれＸ軸方向、Ｙ軸方向およびＺ軸方向を順番に切り替えながら、ボア内の均一磁場に対し傾斜磁場を重畳させ、静磁場に強度勾配を付与する。Ｚコイルは励起時に、磁界強度をＺ方向に傾斜させて共鳴面を限定し、Ｙコイルは、Ｚ方向の磁界印加の直後に短時間の傾斜を加えて検出信号にＹ座標に比例した位相変調を加え（位相エンコーディング）、Ｘコイルは、続いてデータ採取時に傾斜を加えて、検出信号にＸ座標に比例した周波数変調を与える（周波数エンコーディング）。

この重畳される傾斜磁場の切り替えは、制御シーケンスに従って、Ｘコイル、ＹコイルおよびＺコイルにそれぞれ異なるパルス信号が傾斜磁場電源から出力されることにより実現される。これにより、ＮＭＲ現象が発現する被検者２の位置を特定することができ、被検者２の画像を形成するのに必要な三次元座標上の位置情報が与えられる。

ここで、上述のように、３組の直交する傾斜磁場を用いて、それぞれにスライス方向、位相エンコード方向、および周波数エンコード方向を割り当ててその組み合わせにより様々な角度から撮影を行える。たとえば、一般にＸ線ＣＴ装置で撮像されるものと同じ方向のトランスバーススライスに加えて、それと直交するサジタルスライスやコロナルスライス、更には面と垂直な方向が３組の直交する傾斜磁場の軸と平行でないオブリークスライス等について撮像することができる。

ＲＦ照射部１６は、制御シーケンスに従って信号送信部２６から送信される高周波信号に基づいて、被検者２の関心領域にＲＦ（Radio Frequency）パルスを照射するものである。
なお、ＲＦ照射部１６は、図１において、磁場印加機構１１に内蔵されているが、寝台１８に設けられたり、あるいは、受信コイル２０と一体化されて、送受信コイルとして構成されていてもよい。

受信コイル２０は、被検者２からの応答波（ＮＭＲ信号）を検出するものであって、このＮＭＲ信号を高感度で検出するために、被検者２に近接して配置されている。
ここで、受信コイル２０には、ＮＭＲ信号の電磁波がそのコイル素線を切ると電磁誘導に基づき微弱電流が生じる。この微弱電流は、信号受信部２８において増幅され、さらにアナログ信号からデジタル信号に変換されデータ処理部３２に送られる。
（送）受信コイル２０については、上述の通り、ＳＮ比の向上のためにマルチアレイコイルが使用される。

すなわち、静磁界にＺ軸傾斜磁界を加えた状態にある被験者２に、共鳴周波数の高周波電磁界を、ＲＦ照射部１６を通じて印加すると、磁界の強さが共鳴条件になっている部分の所定の原子核、たとえば、水素原子核が、選択的に励起されて共鳴し始める。共鳴条件に合致した部分（たとえば、被験者２の所定の厚さの断層）にある所定の原子核が励起され、（古典的な描像では）スピンがいっせいに回転する。励起パルスを止めると、受信コイル２０には、今度は、回転しているスピンが放射する電磁波が信号を誘起し、しばらくの間、この信号が検出される。この信号によって、被験者２の体内の、所定の原子を含んだ組織を観察する。そして、信号の発信位置を知るために、ＸとＹの傾斜磁界を加えて信号を検知する、という構成になっている。

画像処理部４８は、記憶部３６に構築されているデータに基づき、励起信号を繰り返し与えつつ検出信号を測定し、１回目のフーリエ変換計算により、共鳴の周波数をＸ座標に還元し、２回目のフーリエ変換でＹ座標を復元して画像を得て、表示部３８に対応する画像を表示する。

たとえば、このようなＭＲＩシステムにより、上述したＢＯＬＤ信号をリアルタイムで撮像し、制御部４２により、時系列に撮像される画像について、後に説明するような解析処理を行うことで、安静時機能結合的ＭＲＩ（ｒｓ－ｆｃＭＲＩ）の撮像を行うことが可能となる。

図４においては、ＭＲＩ装置１００．１および他の計測サイトにおけるＭＲＩ装置１００．２～１００．Ｎｓからの計測データ、測定パラメータ、被験者属性データが、データセンター２００内の通信インタフェース２０２を介して記憶装置２１０に集積され、格納される。さらに、計算処理システム３００は、通信インタフェース２０４を介して、記憶装置２１０内のデータにアクセスする構成となっている。

図５は、データ処理部３２のハードウェアブロック図である。
データ処理部３２のハードウェアとしては、上述のとおり、特に限定されないが、汎用コンピュータを使用することが可能である。

図５において、データ処理部３２のコンピュータ本体２０１０は、メモリドライブ２０２０、ディスクドライブ２０３０に加えて、演算装置２０４０と、ディスクドライブ２０３０およびメモリドライブ２０２０に接続されたバス２０５０と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ２０６０とに接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのＲＡＭ２０７０と、アプリケーションプログラム、システムプログラム、およびデータを記憶するための不揮発性記憶装置２０８０と、通信インタフェース２０９０とを含む。通信インタフェース２０９０は、駆動部２１等と信号の授受を行うためのインタフェース部４４および図示しないネットワークを介して他のコンピュータと通信するためのネットワークインタフェース５０に相当する。なお、不揮発性記憶装置２０８０としては、ハードディスク（ＨＤＤ）やソリッドステートドライブ（ＳＳＤ：Solid State Drive）などを使用することが可能である。不揮発性記憶装置２０８０が、記憶部３６に相当する。

演算装置２０４０が、プログラムに基づいて実行する演算処理により、データ処理部３２の各機能、たとえば、制御部４２、データ収集部４６、画像処理部４８の各機能が実現される。

データ処理部３２に、上述した実施の形態の機能を実行させるプログラムは、ＤＶＤ－ＲＯＭ２２００、またはメモリ媒体２２１０に記憶されて、ディスクドライブ２０３０またはメモリドライブ２０２０に挿入され、さらに不揮発性記憶装置２０８０に転送されてもよい。プログラムは実行の際にＲＡＭ２０７０にロードされる。

データ処理部３２は、さらに、入力装置としてのキーボード２１００およびマウス２１１０と、出力装置としてのディスプレイ２１２０とを備える。キーボード２１００およびマウス２１１０が入力部４０に相当し、ディスプレイ２１２０が表示部３８に相当する。

上述したようなデータ処理部３２として機能するためのプログラムは、コンピュータ本体２０１０に、情報処理装置等の機能を実行させるオペレーティングシステム（ＯＳ）は、必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいればよい。データ処理部３２がどのように動作するかは周知であり、詳細な説明は省略する。

また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
また、計算処理システム３００内のハードウェアも、演算処理装置が並列化されていたり、ＧＰＧＰＵ（General-purpose computing on graphics processing units）が使用される等の構成上の相違はあり得るものの、基本的な構成は、図５に示したものと同様である。

（脳機能結合に基づく、疾患／健常の判別器の生成処理、およびクラスタリング処理）
図６は、図２Ｂで説明したような相関行列から、診断マーカーとなる判別器を生成する過程およびクラスタリング処理について説明する概念図である。

機械学習の処理としては、判別器の生成は、いわゆる「教師あり学習」の処理が実行され、クラスタリング処理については、「教師なし学習」の処理が実行される。

そして、クラスタリング処理自体は、「教師なし学習」であって、医師の診断などの情報を用いているわけではないので、この結果得られる個々のクラスタは、データ駆動によって得られる患者のグループであり、患者がサブタイプに分かれる場合は、脳機能結合を特徴量とする「患者の層別化」の基礎となる。

図６に示すように、まず、複数のＭＲＩ装置において、健常群、患者群において安静時のｆＭＲＩ画像データが撮像され、計算処理システム３００が、そのようなｆＭＲＩ画像データに対して、後述するような「前処理」を実施する。続いて、測定された安静時機能結合的ＭＲＩのデータから、それぞれの被験者について、計算処理システム３００が、脳領野のパーセレーション処理を行って、脳領域間（関心領域間）の活動度の相関行列を導出する。

続いて、相関行列の非対角成分について、後述するようにして、対応する測定バイアスを導出しておき、計算処理システム３００は、その測定バイアスを相関行例の要素の値から減算することで、ハーモナイゼーション処理を実行する。

さらに、ハーモナイゼーション処理が行われた相関行列の要素値と、各被験者についての疾患ラベル（疾患あるいは健常を示すラベル）との間で、計算処理システム３００は、過学習を抑制し、特徴選択を伴う識別器の生成を、後述するような「アンサンブル学習による識別器生成処理」として実施して、被験者の疾患あるいは健常を予測することが可能な疾患識別器（診断マーカ）として生成する。

一方で、計算処理システム３００は、アンサンブル学習中に、疾患ラベルに対する識別器の生成処理の中で特定された特徴量（脳機能結合）のうちから、後述するようなクラスタリングのための特徴量選択処理を実施した上で、「教師なし学習」により、多重共クラスタリング処理を実行する。

以下、図６中の各処理について、さらに詳しく説明する。
［前処理から疾患識別の生成およびクラスタリング処理までの概要］
（安静時機能的結合ＦＣマトリックスの前処理および計算）
計測されたｆＭＲＩデータの、たとえば最初の１０秒間はＴ１平衡を考慮に入れるために廃棄される。

前処理のステップでは、計算処理システム３００は、スライスタイミングの校正、頭部に見られる体動アーチファクトの補正を行うためのリアライン処理、脳機能画像（EPI画像）と形態画像の共登録（co-registration）、歪補正、Ｔ１強調構造画像の分割、モントリオール神経学研究所(ＭＮＩ)空間への正規化および、たとえば、６ｍｍの半値幅の等方性のガウスカーネルによる空間の平滑化などの処理を実施する。

このような前処理のパイプライン処理については、たとえば、以下のサイトに開示がある。
http://fmriprep.readthedocs.io/en/latest/workflows.html

（脳領域の区画化（パーセレーション：Parcellation））
脳領域のパーセレーションについては、特に限定されないが、上述した「方法３」に従って、「表面ベースの方法」で実施することが可能である。

（生理的な雑音回帰）
生理的な雑音回帰は、以下の文献１１に開示されるCompCorを適用して実行される。公知文献１１の記載の全体をここに参照により援用する。
公知文献１１：Behzadi, Y., Restom, K., Liau, J., and Liu, T.T. (2007). A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage 37(1), 90-101. doi: 10.1016/j.neuroimage.2007.04.042.
いくつかのスプーリアスの源（余計な信号源）を除去するために、６つの運動パラメーター、全脳、など回帰母数を備えた直線回帰が使用される。

（時間的なフィルタリング）
計算処理システム３００は、時間的なバンドパスフィルターとして、たとえば、０．０１Ｈｚと０．０８Ｈｚの間の通過帯域を備えたバターワース・フィルタを使用して、時系列データに適用し、ＢＯＬＤ活動の特性である低周波の変動に分析を限定する。

（頭部運動）
フレームの位置ずれ（ＦＤ：Frame-wise displacement）が、個々の機能的なセッションにおいて計算され、頭部運動による機能的結合ＦＣのスプーリアスの変化を低減するために、たとえば、ＦＤ＞０．５ｍｍであるボリュームが除去される。
ＦＤは、スカラ量（つまり並進運動と回転における絶対的な変位の加算）として時間的に連続する２つのボリュームの間の頭部運動を表わす。

たとえば、後述する具体例においては、上述したような具体的なデータセットでは、スクラビングの後に除去されたボリュームの比率が、（平均±３×標準偏差）を超過した場合、その参加者のデータは分析から除外される。その結果、データセット全体では、３５人の参加者が除去された。したがって、学習用データセットの中で６８３人の参加者（５４５人のＨＣ、１３８人のＭＤＤ)を使用し、独立した検証データセットの中で４４４人の参加者（２６３人のＨＣ、１８１人のＭＤＤを患者)のデータが、以下の分析に対して使用される。

（機能的結合（ＦＣ）マトリックスの計算）
本実施の形態の具体例においては、機能的結合ＦＣは、上述したようなパーセレーション法により領域の分割がされた後に、各参加者につき３７９個の関心領域（ＲＯＩ）にわたってＢＯＬＤ信号の時間的な相関として計算される。

機能的結合の計算では、特に限定されないが、ここでは、ピアソンの相関係数を使用する。
個々の可能な組のＲＯＩの前処理されたＢＯＬＤ信号の時間的経過間のフィッシャーのｚ変換されたピアソンの相関係数が計算され、要素がそれぞれ２つのＲＯＩの間の結合の強度を表わす、３７９行×３７９列の対称的な結合マトリックスが構築される。
さらに、分析のために、結合マトリックスの下三角行列の７１，６３１（＝（３７９×３７８）／２）の機能的結合ＦＣの値が使用される。

（脳活動バイオマーカーのためのハーモナイゼーション処理）
精神疾患と関連したビッグデータを収集するときには、上述した通り、１つのサイトが大規模な脳画像データ（人間の疾患と関係するコネクトーム）を集めることが、ほとんど不可能であるので、複数のサイトから画像データを収集することが必要である。
ＭＲＩ装置（スキャナー）の型、プロトコルおよび患者層を完全に制御するのは難しい。したがって、収集したデータを解析するには、異種条件下で撮像された脳画像データが用いられる。
特に、疾患要因はサイト要因と交絡する傾向があるので、そのような異種条件下のデータに機械学習技術を適用することで疾患要因を抽出する場合、サイト間差は最も大きな障壁となる。

１つのサイト（あるいは病院）は、ほんの少数のタイプの精神疾患（例えば、サイトＡからは主として統合失調症、サイトＢからは主として自閉症、サイトＣからは主として大うつ病など）をサンプリングしがちであるので、交絡が生じることになる。
適切にそのような異種条件下のデータを管理するために、サイト間のデータをハーモナイズ（調和）させることが必要である。

サイト間差は、本質的に２つのタイプのバイアスを含んでいる。
技術的なバイアス（つまり測定バイアス）および生物学的なバイアス（つまり標本バイアス）である。
測定バイアスは、撮像パラメータ、電界の強さ、ＭＲＩ装置メーカおよびスキャナの型のようなＭＲＩスキャナの特性の違いを含んでおり、標本バイアスは、サイト間の被験者のグループの相違に関連している。

そこで、このようなサイト間差を補償するための「ハーモナイゼーション処理」が必要となる。このハーモナイゼーション処理の詳細については、上述した非特許文献８（Ayumu Yamashita他）に記載されており、その内容については後述する。

（アンサンブル学習による疾患識別器）
本明細書においては、「アンサンブル学習」という用語は、原学習データから復元抽出してＫ組の学習データを作成し、それぞれの学習データについて、独立に、機械学習処理によりＫ個の識別器を生成し、これらのＫ個の識別器を統合して、判別器を生成する処理のことをいうものとする。

特に、ここでは、特定疾患について、ある被験者の脳機能結合パターンにより、その被験者が疾患あるいは健常であるかを判別することを目的とするので、個々の識別器は、２クラス識別問題に対する識別器ということになる。

そして、原学習データから復元抽出してＫ組の学習データを作成する際には、後述するように、「アンダーサンプリング」および「サブサンプリング」が実行される。

ここで、Ｌ１正則化による識別器（ＬＡＳＳＯ（Least Absolute Shrinkage and Selection Operator）法）のような「特徴量選択を伴う学習処理による識別器」や、リッジ正則化法（Ｌ２正則化）などを使用する正則化学習法を用いることも可能である。

ここでは、「正則化学習法」とは、元の訓練データの特徴量全体を学習の対象とするものの、学習アルゴリズムにおいて、モデルを学習する際に、複雑さが増すことに対するペナルティを設け、このペナルティを訓練誤差に加えた量が最も小さくなる学習モデルを求めるようにすることにより汎化性能を高めようとする学習方法のことを意味する。そして、Ｌ１正則化は、ペナルティとして学習モデルの（特徴量に対応する）パラメータの絶対値の総和を用いるものであり、Ｌ２正則化は、ペナルティとして学習モデルのパラメータの二乗の総和を用いるものである。なお、ペナルティとしては、モデルに使用される特徴量の個数そのものを用いるＬ０正則化もあり得る。

また、ＬＡＳＳＯ法（Ｌ１正則化）は、いわゆるスパース推定が可能な手法であり、その派生形として、Elastic Net法、Group Lasso法、Fused Lasso法、Adaptive Lasso法、Graphical Lasso法などがある。

一方で、「特徴量選択を伴う学習処理による識別器」としては、「ランダムフォレスト法」のように、識別器の生成において、特徴量の選択とともに、特徴量の重要度が、併せて、得られる手法を利用することも可能である。

なお、このような「アンサンブル学習による識別器の生成手法」としては、以後は、ＬＡＳＳＯ法を中心として説明するものの、その手法は、上述のようなものに、限定されるものではない。たとえば、パーセレーション法として、辞書学習（Dictionary Learning）の手法を用いてデータ依存的に解析対象脳領域を設定し、機能的脳結合の値としてタンジェント法（tangent-space covariance）を用い、データセット内の脳機能結合ＦＣを後述するComBat法を用いて施設間補正を行い、リッジ正則化で識別器を生成するアンサンブル学習法などであってもよい。パーセレーション法や、脳機能結合の計算方法、ハーモナイゼーション法、識別器の生成方法などは、別の組み合わせであってもよい。たとえば、脳機能結合の計算方法としては、距離相関を用いることもできる。ここで、「距離相関」とは、ピアソン相関法のように脳領域内の活動パターンを平均化せず、活動パターンの類似度を算出する方法であり、以下の公知文献１２および公知文献１３に開示がある。公知文献１２および公知文献１３の各々の記載の全体をここに参照により援用する。
公知文献１２： G.J. Szekely, M.L. Rizzo, N.K. Bakirov, Measuring and testing dependence by correlation of distances, Ann. Statist., 35 (6) (2007), pp. 2769-2794
公知文献１３： https://en.wikipedia.org/wiki/Distance_correlation
後述するように、本実施の形態では、このようなアンサンブル学習において、識別器を学習させる際に、各特徴量の識別という機能を達成するための「重要度」が特定される。

（クラスタリングのための特徴量選択とクラスタリング）
Ｋ組の学習データに対して、「アンサンブル学習」として、Ｋ個の識別器を生成する過程で、それぞれの識別器生成において使用された「第１の特徴量」の和集合から、さらに、「教師なし学習」によるクラスタリングを実行する際の第２の特徴量の集合を特定するための処理が実行される。

特に、限定されないが、たとえば、以下のようにして「重要度」が決定される。

ｉ）アンサンブル学習で、Ｋ個の識別器を生成する学習手法が、「特徴量選択を伴う学習処理」の場合は、識別器の生成において選択された「第１の特徴量」の和集合において、Ｋ個の識別器を生成する際に使用される頻度に応じて、特徴量をランキングする。

ｉｉ）アンサンブル学習で、Ｋ個の識別器を生成する学習手法が、「ランダムフォレスト法」のように、識別器の生成において、特徴量の重要度が得られる手法の場合は、特徴量のランキングリストは、このような重要度に従って生成される構成とすることができる。

ｉｉｉ）アンサンブル学習で、Ｋ個の識別器を生成する学習手法が、（必ずしも特徴量選択を伴わない）「リッジ正則化法（Ｌ２正則化）」であって、特徴量の重み付き和を引数とする識別器の生成手法である場合は、個々の識別器での各特徴量の重み係数の絶対値を、Ｋ個の識別器にわたって集計した中央値を重要度として特徴量のランキングリストが生成される。なお、重要度としては、このような「中央値」に必ずしも限定されるものではなく、たとえば、「Ｋ個の識別器にわたって積算した積算値」など他の統計的な代表値を用いることとしてもよい。

ｉ）～ｉｉｉ）のようにして生成されたランキングのリストの中で、上位の所定の個数の特徴量を、「第２の特徴量」とする、という構成とすることが可能である。

なお、「第２の特徴量」を特定する条件としては、ランキングリストの上位の所定の個数ということに限られず、たとえば、ランキングリストにおいて、所定の頻度以上であること（Ｋ個の識別器の生成で、ある割合以上選択されるという頻度そのものを用いる）を条件としてもよい。

以上のようにして、選択された特徴量に基づいて、後述するように、教師なし学習である「多重共クラスタリング法」により、クラスタリング処理（患者層別化）が実行される。

［２クラス分類のための識別器生成処理］
以下では、図６で説明した処理のうち、アンサンブル学習による識別器生成について、より詳しく説明する。

すなわち、２クラス分類のための分類器生成処理、より特定的には、疾患識別器（「健常」または「疾患」についての２クラス分類器）の訓練データとして、学習用データセットを使用して、ＭＤＤのためのバイオマーカーを構築する処理を例として説明する。

ここでは、精神疾患のうち、大うつ病性障害を例として、すなわち、従来の症候に基づく診断手法により医師によって大うつ病性障害と診断された患者群を例として、分類器を生成する処理を説明する。そして、図８に示した疾患識別器生成部３００８が、患者群と健常群とを判別するための診断の補助情報を出力するような分類器を生成するために実行する処理の例を説明する。
そこで、以下では、機能的結合ＦＣに基づいて、健常群（ＨＣ）およびＭＤＤ患者を識別するＭＤＤ識別器を構築する手続きについて説明する。
以下では、疾患識別器（ＭＤＤ識別器）を作成するための「特徴量選択を伴う学習処理による識別器」としては、例示として、Ｌ１正則化による識別器の学習方法（ＬＡＳＳＯ法）を用いるものとして説明する。
そして、後述するように、ＭＤＤ診断と関係する機能的結合ＦＣを特定するために、各機能的結合ＦＣの疾患識別器の構築に対する「重要度」に応じて、クラスタリングに使用する特徴量が選択される。

図７および図８は、データセンター２００の記憶装置２１０に格納されたデータに基づき、ハーモナイゼーション処理、疾患識別器生成処理、クラスタリング分類器生成処理、判別処理を実行する計算処理システム３００の構成を説明するための機能ブロック図である。

なお、ここでは、「判別処理」には、疾患の判別（疾患あるいは健常の判別）と、対象の被験者がいずれの「クラスタ」（サブタイプ）に属するのかを判断する分類処理を含むものとする。

図７を参照して、計算処理システム３００は、記憶装置２１０からのデータおよび計算途中で生成されるデータを格納するための記憶装置２０８０と、記憶装置２０８０内のデータに対して、演算処理を実行する演算装置２０４０とを含む。演算装置２０４０としては、たとえば、ＣＰＵが相当する。

演算装置２０４０は、プログラムの実行により、患者群および健常群のＭＲＩ測定データ３１０２について相関行列の要素を算出し、相関行列データ３１０６として、記憶装置２０８０に格納するための相関行列算出部３００２と、ハーモナイゼーション処理を実行するハーモナイゼーション算出部３０２０と、ハーモナイゼーション処理の結果に基づいて、疾患識別器の生成処理、クラスタリング分類器の生成処理および生成された疾患識別器またはクラスタリング分類器による判別処理を実行する学習及び判別処理部３０００とを含む。

図８は、図７の構成をより詳しく説明する機能ブロック図である。
また、図９は、アンサンブル学習による疾患識別器を生成するための機械学習の手続きを説明するためのフローチャートである。

そこで、まずは、図６に示したように、ハーモナイゼーション処理、アンサンブル学習による識別器生成（疾患識別器）の生成までの処理について、図８および図９を参照して説明する。

まず、データセンター２００の記憶装置２１０には、各計測サイトから被験者（健常者および患者）のｆＭＲＩ計測データ、被験者の属性データ、測定パラメータが、「学習用データセット」として収集されていることを前提とする。

図８および図９を参照して、疾患識別器（「健常」または「疾患」についての２クラス分類器）の訓練データとして、このような学習用データセットを使用して、ＭＤＤ判別のためのバイオマーカーを構築する。それは、つまり、７１，６３１の機能的結合ＦＣの値に基づいて、健常群（健常（ＨＣ）との診断ラベルの個体群）およびＭＤＤ患者群（大うつ病性障害との診断ラベルの個体群）を識別するものである。

以下に説明するように、ＭＤＤに対する識別器（以下、「ＭＤＤ識別器」と呼ぶ）の生成のための学習処理において、Ｌ１正則化（ＬＡＳＳＯ法）によるロジスティク回帰分析（スパースモデリング法の１種）を使用して、７１，６３１個の機能的結合ＦＣの中からの機能的結合ＦＣの最適な部分集合を選択する。
一般に、Ｌ１正則化を使用すると、いくつかのパラメータ（以下の説明では、ウェイトの要素）を０にすることができる。つまり、特徴選択を行っていることになり、スパースモデルになる。
ただし、スパースモデリングの手法としては、ＬＡＳＳＯ法には限定されず、後述するように、変分ベイズ法をロジスティック回帰に応用したスパースロジスティック回帰（ＳＬＲ：Sparse Logistic Regression）を用いるなど、他の手法を用いることも可能である。

図９を参照して、ＭＤＤ識別器に対する学習処理が開始されると（Ｓ１００）、予め準備され（記憶装置２０８０に格納され）た学習データセットを使用して（Ｓ１０２）、相関行列算出部３００２が、結合マトリックスの成分を算出する。

続いて、ハーモナイゼーション算出部３０２０が、算出した測定バイアスを使用して、ハーモナイゼーション処理を実施する（Ｓ１０４）。
後述するとおり、ハーモナイゼーション処理は、トラベリングサブジェクトを用いた方法が望ましいが、他の方法であってもよい。
たとえば、発見データセットおよび独立検証データセット間で、後述するようなcombat法を用いて、データセット間のハーモナイゼーションを行う構成とすることも可能である。

続いて、疾患識別器生成部３００８は、学習用データに対して、いわゆる「アンサンブル学習の手法」であって、「入れ子構造の交差検証（Nested Cross Validation）」の手法を修正した方法により、ＭＤＤ識別器を生成する。

まず、疾患識別器生成部３００８は、学習用データに対して、「Ｋ分割交差検証」（Ｋ：自然数）（外側の交差検証）を使用して学習処理を実行するために、たとえば、Ｋ＝１０として、学習用データを１０分割する（Ｓ１０６）。

すなわち、疾患識別器生成部３００８は、Ｋ分割（１０分割）されたうちの１つの部分データセットを検証用の「テストデータセット」とし、残りの（Ｋ－１）分割（９分割）分のデータを訓練データセット（トレーニングデータセット）に設定する（Ｓ１０８，Ｓ１１０）。

続いて、疾患識別器生成部３００８は、トレーニングデータセットに対して、「アンダーサンプリング処理」および「サブサンプリング処理」を実行する（Ｓ１１２）。

ここで、「アンダーサンプリング処理」とは、トレーニングデータセットにおいて、分類の対象となる特定の属性データ（２種以上）にそれぞれ対応するデータの数がそろっていない場合に、この数をそろえるために数が多い方の属性のデータを除いて同数となるように行う処理を意味する。

ここでは、トレーニングデータセットにおいて、ＭＤＤ患者群の被験者数と、健常群の被験者数が等しくないために、これをそろえるための処理を行うことに相当する。

さらに、「サブサンプリング処理」とは、トレーニングデータセットから所定数のサンプルをランダムに抽出する処理を意味する。

すなわち、ステップＳ１０８～Ｓ１１８、Ｓ１２２を経由して、Ｋ回繰り返される交差検証において、各交差検証では、トレーニングデータセットがＭＤＤ患者および健常者ＨＣの数に関して不均衡であるので、分類器を構築するためのアンダーサンプリング方法を行うこととし、かつ、サブサンプリング処理として、所定数、たとえば、１３０人のＭＤＤ患者および同数の１３０人の健常者が、トレーニングデータセットからランダムにサンプリングされる。

なお、１３０人という数値は、このような値に限定されるものではなく、学習データセット中のデータ数（後述する「データセット１」では６８３人）、分割数Ｋ（ここでは、たとえば、Ｋ＝１０）、ならびに、分類の対象となる特定の属性に含まれるデータ数の不均衡の程度に応じて、適宜、上述したようなアンダーサンプリングが可能となるように決定される。

このようなサブサンプリング処理を実行するのは、アンダーサンプリングは、除外されたデータを使用して識別器が学習することができなくなるという点で不利であるため、この不利を除去するために、無作為抽出手続き（つまりサブサンプリング）を、Ｍ回（Ｍ：自然数、たとえば、Ｍ＝１０）繰り返す処理とするためである。

なお、後述するように「層別化」のための「分類器」を生成する際の「特徴量選択」のためにも、このようなアンダーサンプリング処理およびサブサンプリング処理を実行することには、技術的な意義が存在するので、その点については、後述する。

続いて、疾患識別器生成部３００８は、サブサンプリングされたサブサンプル１～１０の各々に対して、ハイパーパラメータの調整処理を実行する（Ｓ１１４．１～Ｓ１１４．１０）。

ここで、各サブサンプルにおいて、以下のようなロジスティック関数を使用することで、識別器サブモデルが生成される。このようなロジスティック関数は、サブサンプル内において、ＭＤＤクラスに属する参加者の可能性を以下のように定義するために使用される。

ここで、ｙ_subは、参加者のクラス・ラベル(ＭＤＤ, ｙ＝１；ＨＣ, ｙ＝０)を表わし、ｃ_subは、与えられた参加者に対するＦＣベクトルを表わし、ｗはウェイトベクトルを表わす。

ウェイトベクトルｗは、以下の評価関数（コスト関数）を最小化するように決定される（ＬＡＳＳＯ計算）。

ＬＡＳＳＯ計算では、コスト関数において、ウェイトベクトルの各要素の絶対値（１次）の総和（Ｌ１ノルム）が第２項として存在する。

ここで、λはハイパーパラメータを表わし、評価に適用される収縮量をコントロールする。

疾患識別器生成部３００８は、各サブサンプルにおいて、これも特に限定されないが、所定数のデータをハイパ―パラメータ調整用データとし、残りのデータ（たとえば、ｎ＝２５０または２４８人のデータ）を用いてウェイトベクトルｗを決定する。このとき、疾患識別器生成部３００８は、特に限定されないが、たとえば、ハイパーパラメータλが０＜λ≦１．０であるものとして、この区間を、Ｐ等分（Ｐ：自然数）、たとえば２５等分した各値のλを使用して、上記のようなＬＡＳＳＯ計算により、ウェイトベクトルｗを決定する。

このとき、上述の通り、「入れ子構造の交差検証」として、ハイパーパラメータの調整については、「内側の交差検証」として実行する。内側の交差検証では、外側の交差検証の「テストデータセット」は、使用されない。

その上で、疾患識別器生成部３００８は、生成された各λの値に対応するロジスティック関数によりハイパ―パラメータ調整用データについて判別性能（たとえば、正確度）を比較し、最も判別性能の高いλに対応するロジスティック関数を決定する（ハイパーパラメータの調整処理）。

続いて、疾患識別器生成部３００８は、現在の交差検証のループの中で生成された各サブサンプルに対応するロジスティック関数の出力値の平均を出力するものとして「識別器サブモデル」を設定する（Ｓ１１６）。各サブサンプルにおいて算出された識別器の出力値の平均により識別性能を判定しているという点でも、これは、「アンサンブル学習」の一種ということができる。

疾患識別器生成部３００８は、ステップＳ１１０で準備されたテストデータセットを入力として、現在の交差検証のループの中で生成された識別器サブモデルの検証を実行する（Ｓ１１８）。

なお、アンダーサンプリングとサブサンプリングにより、サブサンプルを生成して、各サブサンプルにおいて、特徴選択を実行して識別器サブモデルを生成する方法としては、上記のようなＬＡＳＳＯ法とハイパーパラメータの調整とを実行する方法以外に、他のスパースモデリングの手法を利用してもよい。

疾患識別器生成部３００８は、交差検証のループをＫ回（ここでは１０回）分終了していないと判断すると（Ｓ１２２でＮＯ）、Ｋ分割されたデータにおいて、それまでのループで使用したのとは異なる別の部分データセットをテストデータセットに設定し、残りの部分データセットをトレーニングデータセットに設定して（Ｓ１０８，Ｓ１１０）、処理を繰り返す。

一方で、疾患識別器生成部３００８は、交差検証のループをＫ回（１０回）分終了している場合（Ｓ１２２でＹＥＳ）、入力データに対して、Ｋ×Ｍ個（この場合は、１０×１０＝１００個）のロジスティック関数（識別器）の出力の平均を出力するものとして、ＭＤＤに対する識別器モデル（ＭＤＤ識別器）を生成する（Ｓ１２０）。

結局のところ、ＭＤＤ識別器は、Ｋ×Ｍ個の識別器の出力の平均を、その識別出力とするという意味で、「アンサンブル学習」の結果得られる「識別器」であるといえる。

ＭＤＤ識別器の出力（診断の確率値）が、０．５を超えるときは、ＭＤＤ患者を示す指標と見なすことができる。

さらに、本実施の形態では、このようにして生成されたＭＤＤ識別器の性能の評価指標として、マシューズ相関係数（ＭＣＣ：Matthews correlation coefficients）、ＲＯＣ曲線(Receiver Operatorating Characteristic curve、受信者動作特性曲線)についてのＲＯＣ曲線下面積（ＡＵＣ: area under the curve)、正確度（Accuracy）、感度（Sensitivity）、特異度（Specificity）を使用する。

なお、それぞれのサブサンプルにおいて特徴選択された特徴量（この場合、測定バイアスに対するハーモナイゼーション処理後の相関行列の要素）を使用して、対象疾患（たとえば、ＭＤＤ）の識別器を生成する手法は、このような複数個の識別器サブモデルの出力の平均処理によるものに限られず、多数決による処理としてもよいし、あるいは、特徴選択された特徴量に対して、他のモデリング手法、特に、他のスパースモデリング手法を使用して識別器を生成する構成としてもよい。

（ＭＤＤ識別器に用いるデータの例と性能）
既述したように、機械学習アルゴリズムを使用する信頼できる分類器および回帰モデルの構築については、多数の撮像サイトから集められた大規模なサンプルサイズのデータを用いることが必要である。

そこで、以下では、４つの異なる撮像サイトから集められたＭＤＤ患者を含む、約７００人の参加者の学習用の安静時ｆＭＲＩデータセットを使用して検討する。

図１０は、このような学習用のデータセット（データセット１）のデモグラフィック特性を示す図である。

データセット１は、上述したＳＲＰＢＳ中のデータである。

図１１は、独立検証のデータセット（データセット２）のデモグラフィック特性を示す図である。

データセット２も、基本的には、上述したＳＲＰＢＳ中のデータである。

すなわち、以下の分析では、２つの安静時機能的ＭＲＩ（ｒｓ－ｆＭＲＩ）データセットを使用する
（１）図１０に示すように、データセット１は、７１３人の参加者に対するデータを含んでいる（４つのサイトから５６４人の健常群ＨＣ、３つのサイトからの１４９人のＭＤＤ患者群）。
（２）図１１に示すように、データセット２は、４４９人の参加者に対するデータを含んでいる（４つのサイトからの２６４人の健常群ＨＣ、４つのサイトからの１８５人のＭＤＤ患者群）。

また、併せて「抑うつの症候」は、各データセットの大部分の参加者から得られたベックうつ病自己評価尺度（ＢＤＩ：Beck Depression Inventory）ＩＩを使用して評価されている。
データセット１は、「学習用データセット」であり、ＭＤＤの識別器およびクラスタリングの分類器を構築するのに使用される。

参加者は、それぞれ１０分間の単一の安静時機能的ＭＲＩセッションにて計測が実行されている。
ここでも、統一された撮像プロトコルの下で安静時機能的ＭＲＩデータをが取得されている（http://www.cns.atr.jp/rs-fmri-protocol-2/）。
ただし、画像撮影が、すべてのサイトで同じパラメーターを使用して行われたことを保証することは実際には困難であり、２つの位相変調方式の方向（Ｐ→ＡとＡ→Ｐ）、２社のＭＲＩ装置メーカー（シーメンスとＧＥ）、コイルの３つの異なる数（１２，２４，３２）および３つの型番のスキャナが計測に使用されている。

安静時機能的ＭＲＩのスキャン中に、参加者は、原則として、以下のような教示がされている。
「リラックスしてください。眠らないでください。中央の十字線マーク上を注視して、特定のことに関して考えないでください。」
データセットにおける「デモグラフィック特性」は、いわゆる「人口統計学」において使用される特性であって、年齢、性別、などの他、診断名のような表中の属性を含む。
なお、図１０および図１１において、括弧中の人数は、ＢＤＩスコアのデータを有する参加者の数を示す。
デモグラフィック分布は、すべて学習用データセットの中でＭＤＤとＨＣの個体群間で統計的に有意な差はない（ｐ＞０．０５）。

データセット２は、「独立した検証データセット」であり、ＭＤＤの分類器およびクラスタリングの分類器をテストするために使用される。
データセット２に対する撮像を行ったサイトは、データセット１に含まれていない。

年齢のデモグラフィックな分布は独立した検証データセットの中でＭＤＤとＨＣの個体群間で一致しているものの（ｐ＞０．０５）、性比のデモグラフィックな分布は独立した検証データセットの中でＭＤＤとＨＣの個体群間で一致していない（ｐ＜０．０５）。

（サイト効果の制御）
また、以下では、機能的結合ＦＣの上のサイト効果を制御するために、後述するような、学習用データセットに対するトラベリングサブジェクトのハーモナイゼーション方法を使用するものとして説明する。
ただし、ハーモナイゼーション方法としては、この方法に限定されるものではなく、たとえば、これもＣｏｍＢａｔ法など他の方法を使用してもよい。

なお、ＣｏｍＢａｔ法については、たとえば、以下の公知文献１４に開示がある。公知文献１４の記載の全体をここに参照により援用する。

公知文献１４：Johnson WE, Li C, Rabinovic A. “Adjusting batch effects in microarray expression data using empirical Bayes methods.” Biostatistics 8, 118-127 (2007).
トラベリングサブジェクトのハーモナイゼーション方法を使用することで、純粋なサイト間差（測定バイアス）を除去することが可能となる。
なお、独立した検証データセットに含まれたサイトに対しては、トラベリングサブジェクトのデータセットが存在しなかったので、独立した検証データセットの中でサイト効果のコントロールのためには、ＣｏｍＢａｔ法によるハーモナイゼーション方法を使用している。

図１２は、全撮像サイトについて、学習用データセットに対するＭＤＤの予測性能（出力の確率分布）を示す図である。
学習用データセットに対して、０．５のしきい値によって、識別器モデルからの出力において、ＭＤＤ患者と健常者の個体群に対応する２つの診断の確率分布が、明白に、右（ＭＤＤ）と左（ＨＣ）へ分離されている。

識別器モデルは、ＭＤＤ患者を６６％の正確度でＨＣ個体群から分離する。
対応するＡＵＣは、０．７７であり、高い識別力を示した。
また、ＭＣＣは、約０．３３である。

図１３は、各撮像サイトについて、学習用データセットに対するＭＤＤの予測性能（識別器の出力の確率分布）を示す図である。
図１３から、全データセットだけでなく３つの撮像サイト（サイト１、サイト２、サイト４）の個々のデータセットに対して、ほとんど同じ程度に高い分類精度が達成されることがわかる。
なお、サイト３（ＳＷＡ）のデータセットには、単に健常者群があるだけであるが、その確率分布は、他のサイトの健常者群のものに相当している。

（識別器の汎化性能）
図１４は、独立した検証データセットにおけるＭＤＤの識別器の出力の確率分布を示す図である。
すなわち、独立した検証データセットを使用して、識別器モデルの汎化性能がテストされる。
ＭＤＤに対して、図１２の処理において１００個（１０分割×１０サブサンプリング）のロジスティック関数の識別器が機械学習により生成されており、独立した検証データセットを、生成された１００個の識別器のすべて（識別器の集合としての識別器モデル）に入力する。

そして、各参加者につき１００個の識別器の出力の平均（診断の確率）をとり、平均された診断の確率値が＞０．５だった場合、その参加者の診断ラベルとして、大うつ病性障害が該当するものとする。

独立した検証データセットにおいて、約７０％の精度で、生成された識別器モデルは、ＭＤＤ個体群をＨＣ個体群から分離している。
対応するＡＵＣは、０．７５となり、高い識別能力となる能力を示した（並べかえ検定ｐ＜０．０１）。
独立検証データセットに対して、０．５のしきい値によって、識別器モデルからの出力において、ＭＤＤ患者と健常者の個体群に対応する２つの診断の確率分布が、明白に、右（ＭＤＤ）と左（ＨＣ）へ分離されている。
感度は６８％であり、特異度は７１％である。これは、０．３８という高いＭＣＣ値となっている（並べかえ検定ｐ＜０．０１）。

図１５は、各撮像サイトについて、独立検証データセットに対するＭＤＤの識別器の出力の確率分布を示す図である。
４つの撮像サイトの全データセットだけでなく個々のデータセットに対して高い分類精度が達成されることがわかる。

［被験者データに対するクラスタリング処理］
以下では、図６で説明した処理のうち、クラスタリングのための特徴量選択と、選択された特徴量によるクラスタリング処理について、より詳しく説明する。
すなわち、図６において、「特徴量選択」と「クラスタリング処理」として説明した処理を、図８における疾患識別器生成部３００８とクラスタリング分類器生成部３０１０とが実行する処理として説明をする。

図１６は、特徴量を選択して、教師なし学習によりクラスタリングを行う処理を説明するためのフローチャートである。

以下では、図９において説明したような「２クラス識別器」の学習処理において、各識別器サブモデルの生成において使用された特徴量（脳機能結合）を、ランク付けして、上位の所定数の特徴量を使用して、教師なし学習によりクラスタリングを行う処理について説明する。

上述したように、脳機能結合は、脳に対するパーセレーションの手法に応じて、７００００次元を超えるような高次元であって、通常の方法では、教師なし学習でのクラスタリング処理を実行することは一般には困難である。本実施の形態の方法では、このようなクラスタリング問題に対して、「教師あり学習での識別器生成」において、特徴量の重要度に応じたランキングを実行し、このランキングに基づいて選択された特徴量による「教師なし学習によるクラスタリング」を組み合わせることで、このようなクラスタリング処理を可能としている。

なお、以下では、便宜上、疾患識別器の生成処理とは別のものとして、クラスタリング処理を説明するが、図１６において、ステップＳ２００～Ｓ２１０は、図９におけるステップＳ１００～Ｓ１２０と同等の処理であり、疾患識別器の生成処理とクラスタリング処理とは、一連の処理として実行することが可能である。

図１６を参照して、クラスタリングの学習処理が開始されると、疾患識別器生成部３００８は、健常群Ｎｈ人、うつ病群Ｎｍ人の被験者データを準備し（Ｓ２０２）、疾患識別器生成部３００８は、被験者の機能的脳活動データにおいて、脳領野の分割（パーセレーション）処理、脳機能結合値の算出、および、ハーモナイゼーション処理を実行する（Ｓ２０４）。

続いて、疾患識別器生成部３００８は、Ｎcv分割交差検証（Ｎcv：自然数、かつ、Ｎcv≧２）のためのデータ分割を実行して、各分割されたデータについて、トレーニングデータセットと、テストデータセットとを準備し、各トレーニングデータセットについては、Ｎs個のテストデータのサブセットを生成するように、アンダーサンプリングおよびサブサンプリングを実行する（Ｓ２０６）。

さらに、疾患識別器生成部３００８は、サブサンプルされた各サブサンプルに対して、特徴量選択を伴う学習方法により識別器を生成する（Ｓ２０８）。
なお、ここでは、図９と同様にして、Ｌ１正則化（ＬＡＳＳＯ）で特徴量選択が実行されるものとする。

このようなステップＳ２０６～Ｓ２０８の処理を、Ｎcv個に分割された学習データセットについて、トレーニングデータセット（分割されたデータセットの（Ｎcv－１）個分）とテストデータセット（分割されたデータセットの１個分）との組み換えを順次実行して、交差検証がＮcv回実施されるまで、繰り返す。

このようにして生成された(Ｎs×Ｎcv）個の識別器の平均を出力とする統合された識別器を、疾患識別器（診断マーカ）として生成する（Ｓ２１０）。

ここまでの処理は、上述したとおり、図９におけるステップＳ１００～Ｓ１２０と同等の処理である。

一方で、クラスタリング分類器生成部３０１０は、Ｎcv回繰り返されるステップＳ２０６～Ｓ２０８において、特徴量選択を伴う学習方法により識別器を生成する際に選択された特徴量（脳機能結合）の和集合において、特に限定されないが、この和集合中の特徴量について、選択された回数でランキングを実行する（Ｓ２２０）。

ここでは、「特徴量として選択された回数」を、このランキングにおける、その特徴量の重要度と呼ぶことにする。

言い換えると、ＬＡＳＳＯ法で、たとえば、図９に示した例によれば、１００（＝１０×１０）個の識別器が生成されるが、各識別器において重みがゼロでない脳機能結合について、＋１となるように、選択回数をカウントする。カウント回数の大きい順で、重要な結合としてランキングする。

続いて、クラスタリング分類器生成部３０１０は、うつ病患者群に対してクラスタリングを教師なし学習で実施するために、上記和集合中から重要度に応じて、たとえば所定数の特徴量を選択する（Ｓ２２２）。

さらに、クラスタリング分類器生成部３０１０は、教師なし学習の方法として、後述するような多重共クラスタリング法を用いて、クラスタリング処理を実施する（Ｓ２２４）。

以上の処理により、クラスタリング分類器生成部３０１０は、うつ病患者群に対するクラスタリング分類器を生成する（Ｓ２２６）。

すなわち、以上の処理により、クラスタリング分類器生成部３０１０により、各クラスタについて、観測データから、そのような観測データを生成するような確率分布のモデルがそれぞれ特定され、記憶装置２０８０に各モデルの情報が格納される。そして、判別値算出部３０１２は、クラスタリング分類器として、学習データ以外の入力データに対しては、このよう各確率分布のモデルに基づいて、それぞれのクラスタに当該入力データが属する事後確率を算出し、最大の事後確率となるクラスタに、当該入力データが属するとの分類結果を出力することになる（ＭＡＰ推定法（Maximum A posteriori Probability estimation method））。

なお、以上の説明では、クラスタリング処理は、疾患識別器生成部３００８により、健常群Ｎｈ人、うつ病群Ｎｍ人の被験者データに対して実施された識別器の生成処理に基づくものとして説明した。しかし、本実施の形態のクラスタリング方法は、このような場合に限定されることなく、「うつ病患者群」以外の疾患群、たとえば、「統合失調症の患者群」「自閉症の患者群」「強迫神経症の患者群」など、他の精神疾患の患者群のクラスタリングに使用することも可能である。

あるいは、より一般的に、人間が経験的に分類してきた属性ラベル（たとえば、その人の性格、その人の得意分野など）と、脳活動の時間変化の領野間の相関のパターンとが一定の関係性を有することが判明した属性については、当該属性ラベルにより分類される被験者について、データ駆動により「当該属性に属する被験者群のクラスタリング」（サブタイプへの分類）を実施することに使用することも可能である。

（アンダーサンプリングおよびサブサンプリング処理）
以上説明した処理では、「アンダーサンプリングおよびサブサンプリング処理」を実行しているので、その技術的意義について、簡単に説明する。
まず、「アンダーサンプリング処理」の効果については、識別器における識別の境界が適切に設定されることが挙げられる。
たとえば、２クラス分類のタスクを考えた場合、学習データ中で、各クラスに属するデータの数に偏りが少ないほど、処理フロー中で、識別器の性能（たとえば、正確度）の評価の精度が高くなる。

図９におけるステップＳ１１４．１～１１４．１０において、ハイパーパラメータの設定には、「最も判別性能の高いλに対応するロジスティック関数を決定する」との処理を実行するので、「判別性能」の評価が正確に行われることが必要になる。

極端な例として、クラス１に属するデータ数が１００個で、クラス２に属するデータ数が１個というような学習データについて識別器の学習をする場合、識別器として、すべてのデータがクラス１であると判別したとしても、正確度等には、大きな影響が出ない、という事態が発生してしまう。この点で、２つのクラスにそれぞれ属するデータ数を、ランダムなサンプリングによりそろえることには、意義がある。

また、サブサンプリングについては、そもそも、これは、以下の理由により、アンダーサンプリングと併せて、複数回実施することが前提である。
第１には、アンダーサンプリングおよびサブサンプリング処理は、仮に、ランダムサンプリングにより実施されるとしても、１回の処理だけでは、データに偏りが生じる可能性がある。
第２に、以下に説明するように、図９のステップＳ１０８～Ｓ１２２で繰り返される「識別器」の生成は、上述したように、「アンサンブル学習」により実施される。
このとき、各識別器の生成において、識別に対する特徴量の重要度が決定される。

重要度は、「特徴量選択を伴う学習処理」においては、当該特徴量が選択されること、あるいは、「特徴量選択を伴わない学習処理」においては、算出された識別に対する当該特徴量の重みというように、識別に対する各特徴量の貢献度に応じて、決定される。

以下では、このような重要度の決定における「アンダーサンプリングおよびサブサンプリング処理」の意義を、「特徴量選択を伴う学習処理」を例として説明する。なお、Ｌ２正則化のような「特徴量選択を伴わない学習処理」であっても、「当該特徴量の重みが大きくなる」という事象は、「特徴量として選択される」という事象と、基本的には、同様の技術的理由から生じるものと考えることができる。

ここで、「特徴量選択を伴う学習処理」としては、たとえば、上述したＬＡＳＳＯ法のような、いわゆる「スパースモデリング」の手法が挙げられる。
スパースモデリングでは、特徴量がスパースに選択される、すなわち、特定の特徴量については、重みが非ゼロとなるのに対して、他の特徴量については、重みが０となることで、特徴量が選択される。このような特徴量のスパースな選択が実現される理由の１としては、「判別（識別）処理」に対して「類似の寄与をする特徴量のグループ」が存在する場合に、当該グループのうちの１つの特徴が選択され、当該グループの他の特徴量の重みが０となるように、学習処理が進行させるための「特徴量の数に対応するペナルティ項」が存在していることが、挙げられる。ＬＡＳＳＯ法については、特に、この傾向が顕著である。

つまり、「判別処理」に対して、特徴量Ａと特徴量Ｂが、同様に関与する、たとえば、特徴量Ａと特徴量Ｂの相関が高い、という場合には、特徴量として、特徴量Ａのみを選択しても、判別性能を落とすことなく、判別処理を行うことができることになる。

しかしながら、たとえば、クラスタリングにおいては、特徴量Ａと特徴量Ｂとを双方考慮する必要がある場合が想定される。ところが、このような「スパース化」を実行して、１回の識別器の生成過程における判別処理に対する寄与度だけで、特徴量を選択してしまうと、クラスタリングに対する「特徴量選択」としては、十分でない可能性が生じてしまう。

図１７は、このような「特徴量選択を伴う学習処理」により、複数個（たとえば、Ｎch個）の特徴量が存在する場合に、特徴量の選択が実施される概念を示す図である。

図１７を参照して、被験者群には、健常群と患者群とが含まれるものとする。

健常群の被験者には、ラベルＨが対応付けられており、健常群には、サブタイプｈ１とサブタイプｈ２とが含まれるものとする。

患者群の被験者には、ラベルＭが対応付けられており、患者群には、サブタイプｍ１、サブタイプｍ２、サブタイプｍ３とが含まれるものとする。

ここで、観測量では、健常群、被験者群が、いくつのサブタイプに分かれるのかは未知であり、サブタイプの識別ラベルは明示的には被験者に対応付けられておらず潜在的なラベルであるものとする。

そして「クラスタリング」の目的は、観測量からこれらのサブタイプへのデータ駆動によるクラスタリングを実行することである。

上記のような健常群と患者群とから、被験者を「アンダーサンプリング」および「サブサンプリング」をランダムに実施するので、図１７の「被験者群」に示すように、健常群および患者群のそれぞれから、たとえば、点線で囲んだ部分の被験者が選択されるものとする。

さて、ラベルＭとラベルＨを識別するために使用し得る特徴量（脳機能結合の相関値）は、各被験者を特徴づける特徴量（全体でＮch個）の脳機能結合のうち、図１７において一点鎖線で示した範囲の特徴量（図１６のステップＳ２２０における「脳機能結合の和集合」）であるものとする。

そして、特徴量選択を伴う学習処理（ここでは、ＬＡＳＳＯ法による処理）により、ラベルＭとラベルＨを識別するための識別器を学習した結果として、図１７における一点鎖線内のさらに、黒丸で示した特徴量が選択されるものとする。

図１８は、アンダーサンプリングおよびサブサンプリングの処理後に、特徴量選択を伴う学習処理により１つの識別器を生成する際に最終的に選択される特徴量を示す概念図である。

図１８に示すように、一点鎖線内のラベルＭとラベルＨを識別するために使用し得る特徴量は、さらに点線の枠で示すように、互いに相関が強い特徴量のグループに分かれているものとする。

ＬＡＳＳＯ法では、このような点線の枠内のグループごとに、１つの特徴量が選択されることで、スパース化が実現していることになる。

図１９は、アンダーサンプリングおよびサブサンプリングの処理を複数回実施して、識別器を生成される際に、特徴量が選択される様子を示す概念図である。

図１９に示すように、サブサンプリングを、たとえば、Ｎs回実施するものとすると、各回において、健常群と患者群から、それぞれ異なる被験者がサブサンプルされる。

そして、各サブサンプルに対して、特徴量選択を伴う学習処理により、ラベルＭとラベルＨを識別するための識別器を学習した結果として、各サブサンプル対応して、上述した和集合である一点鎖線内のうちで、相関の強い各グループからは、黒丸で示したように、それぞれの識別器で異なる特徴量が選択されることになる。

結果として、アンダーサンプリングおよびサブサンプリングの処理を複数回実施することで、ラベルＭとラベルＨを識別するために使用し得る特徴量の和集合が選択されることになる。

本実施の形態では、各サブサンプルに対して、ＬＡＳＳＯ法により、特徴量選択を伴う識別器の学習処理により選択される特徴量を、選択される頻度に応じて、ランク付けする。

そして、ランク付けのうちの上位から所定の個数、たとえば、１００個の特徴量を用いて、図１６のステップＳ２２４において、後述するような「多重共クラスタリング」により、教師なし学習によるクラスタリングを実施する。

なお、以上の説明では、「特徴量選択を伴う識別器の学習処理」として、ＬＡＳＳＯ法を例として説明し、選択された頻度のランキングでクラスタリングのための特徴量の選択を実施する構成としている。

ただし、上述のとおり、本実施の形態のクラスタリング処理において、「特徴量選択を伴う識別器の学習処理」は、このような方法に限定されず、たとえば、ランダムフォレスト法のような方法でもよく、所定の重要度に従って、クラスタリングのための特徴量の選択を実施することも可能である。

たとえば、上述したように、ランダムフォレスト法では、識別器の学習処理において、ジニ不純度(Gini impurity)やpermutation importanceをもとに特徴量の重要度(importance)が併せて算出されるので、この重要度に基づいて、特徴量をランク付けし、ランク付けのうちの上位から所定の個数の特徴量を用いて、図１６のステップＳ２２４における「多重共クラスタリング」により、教師なし学習によるクラスタリングを実施するとの構成とすることも可能である。

なお、ここで、permutation importanceとは、「ランダム化された特徴量で作ったモデルの誤差ともともとのモデルの誤差の差分」であり、「どの特徴量がもっともモデルの精度に寄与しているか、もしくは精度に寄与していないのか」を算出するものである。たとえば、Permutation Importanceについては、以下の公知文献１５に開示がある。公知文献１５の記載の全体をここに参照により援用する。
公知文献１５：Breiman, Leo. "Random forests." Machine learning 45.1 (2001): 5-32. https://link.springer.com/article/10.1023/A:1010933404324

また、「アンサンブル学習として識別器を学習する処理」としては、リッジ正則化法などを使用し、上述したように、重み係数の絶対値を集計した中央値に対応する重要度に従って、特徴量をランキングして、ランク付けのうちの上位から所定の個数の特徴量を用いて、クラスタリングのための特徴量の選択を実施することも可能である。

［多重共クラスタリング処理］
以下では、図１６のステップＳ２２４における「多重共クラスタリング」について、その概念を説明し、「多重共クラスタリング」との用語の定義を行う。

前提として、「クラスタリング」とは、コンピュータにより実行される教師なし学習によるデータ分類手法をいい、より特定的には、与えられたデータを外的基準なしに自動的に分類する手法をいうものとする。これに対して、「クラス分類」は、一般に、「教師あり学習」による分類手法のことをいう。また、「クラスタ」とは、内的結合と外的分離の性質を持つデータの部分集合である、と定義される。ここで、外的分離とは、違うクラスタにある対象は類似していないという性質をいい、内的結合とは、同じクラスタ内の対象は互いに類似しているという性質をいう。さらに、「類似」の尺度として、集合の要素間の距離が定義される。一般には、いわゆる「距離の公理」を満たすように距離が定義され、距離としては、ユークリッド距離、マハラビノス距離、シティブロック距離、ミンコフスキー距離などが使われることがある。

また、一般に、教師なし学習によりクラスタリングを行う手法としては、非階層型手法である「k-means法」などのように、クラスタの良さを定義する目的関数を最適にする分割を探索する方法である「分割最適化クラスタリング」や、階層的クラスタリング手法である「凝集型階層的クラスタリング」や「分割型階層的クラスタリング」などと呼ばれる手法が知られている。

ただし、このような従来のクラスタリング手法は、特徴量をすべて使って対象をクラスタ（グループ）に分け、得られるクラスタの分け方は一通りとなる、という特徴がある。

したがって、特徴量に依存して複数のクラスタの分け方が存在する場合、うまく対処できないという問題が存在することになる。一般に、特徴量の数が多くなるほど、このような複数のクラスタ構造が存在する可能性が高いと考えられる。

なお、クラスタリングの手法としては、上記のような手法だけではなく、むしろ、クラスタリングを実施する複数の対象が、それぞれのクラスタにおいて、ある確率分布に従って生起したものと想定して、このような「確率分布」を推定するという方向で、クラスタリングを実行するというアルゴリズムもあり、このようなクラスタリング手法としては、たとえば、「混合ガウス分布によるクラスタリング法」などが知られており、より柔軟なクラスタリングを実行することが可能であることが知られている。

以下では、まずは、特徴量に依存して複数のクラスタの分け方が存在する場合を説明するために、クラスタリングを実行する複数の対象を含んだ対象群について、各対象は、複数の特徴量で特徴づけられるものとする。

図２０は、特徴量に依存して複数のクラスタの分け方が存在する場合を説明するための概念図である。
図２０に示すように、クラスタリングの対象となるデータ（以下、単に「対象」と呼ぶ）が、６つの文字「Ａ」「Ｂ」「Ｃ」「Ｄ」「Ｅ」「Ｆ」であるものとする。
そして、これらの文字は、異なる背景パターンや、異なる書体（文字スタイル）を有しているものとする。

そこで、これらの文字を特徴づける特徴量としては、「背景パターン」「文字スタイル」「文字に含まれる穴の数（線で完全に囲まれた領域の数）」を考えることができる。
したがって、同一の文字の集合を考えた場合でも、いずれの特徴量に基づいて、クラスタリングするかによって、異なるクラスタに分けられることになる。

図２０では、たとえば、「背景パターン」に基づく場合、｛Ａ，Ｄ｝｛Ｂ，Ｅ｝｛Ｃ，Ｆ｝の３つのクラスタに分けられ、「文字スタイル」に基づく場合、｛Ａ，Ｂ，Ｃ｝｛Ｄ，Ｅ，Ｆ｝の２つのクラスタに分けられ、「穴の数」に基づく場合、０個、１個、２個にそれぞれ対応して、｛Ｃ，Ｅ，Ｆ｝｛Ａ，Ｄ｝｛Ｂ｝の３つのクラスタに分けられる。
図２０では、１つのクラスタは、１つの特徴量により特徴づけられる場合を例としているが、一般には、１つのクラスタは、複数個の特徴量により特徴づけられる。

図２１Ａおよび図２１Ｂは、複数の特徴量により、複数の対象が特徴づけられる場合のクラスタリングの概念を説明するための概念図である。
まず、図２１Ａに示すように、クラスタリングする対象が行方向に配置され、これらの対象を特徴づける特徴量が列方向に配置されている「データ行列」を考える。
図２１Ｂに示すように、対象をクラスタリングする（対象を複数の対象クラスタに分割する）と同時に、各対象クラスタと関連するように特徴量もクラスタリングする手法は、「共クラスタリング」と呼ばれ、例えば、以下の公知文献１６にその手法が開示されている。公知文献１６の記載の全体をここに参照により援用する。

公知文献１６：Madeira SC, Oliveira AL. Biclustering algorithms for biological data analysis: a survey. IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB). 2004; 1(1):24±45. https://doi.org/10.1109/TCBB.2004.2
「共クラスタリング」では、図２１Ｂに示すように、データ行列の行、あるいは列を、入れ替えることで、すなわち、対象および特徴量を類似度によって再配列することで、たとえば、（ｉ，ｊ）（ｉ＝１，２：ｊ＝１，２，３）で示されるクラスタブロックに分割する。

このとき、各クラスタに含まれる対象ついての生成モデル（確率モデル）を仮定して、観測データに対して、その尤度が高くなるように各確率モデルのパラメータを決定する。

このように、各クラスタについて、確率モデルが推定されると、特定の観測データ（テストデータ）について、当該データが、いずれのクラスタに属するかが判別できる（分類できる）ことになる。

図２２Ａおよび図２２Ｂは、多重クラスタリングおよび多重共クラスタリングを説明するための概念図である。

図２１Ｂに示した「共クラスタリング」では、「データ行列」の行と列を入れ替えることにより、ブロック構造のクラスタを生成したので、特徴量が、複数の特徴量クラスタに分割された場合に、各対象は、この複数の特徴量クラスタに共通に並んでいるものとしてクラスタリングがされていることになる。

ただし、特徴量を複数の特徴量クラスタに分割し、特徴量クラスタごとに対象も対象クラスタに分割することを想定する場合、特徴量クラスタごとに、対象クラスタ中の対象の並び方（ある対象クラスタ中に含まれる対象の並び方）も異なるものとした方が、より尤度の高い確率モデルが推定できることが想定される。

このような場合、特徴量クラスタごとに、対象の分割方法（対象のクラスタリング）が異なることに対応して、特徴量クラスタを、特に、「ビュー（視点）」と呼ぶ。

以上のようにして、特徴量の視点ごとに、異なる対象のクラスタリングを実行することを、図２２Ａに示すように「多重クラスタリング」と呼ぶ。

さらに、各視点において、特徴量の列と対象の行を入れ替えてクラスタリングすることで、さらに、観測データに対する尤度の高い確率モデルが推定できる場合は、図２２Ｂに示すように、「多重共クラスタリング」と呼ぶ。

ここで、ビューが１つしかない場合や、複数のビューが存在する場合に、少なくとも１つのビューにおいて特徴量クラスタが１種類しかないときも含めて、「多重共クラスタリング」と呼ぶこととし、「共クラスタリング」および「多重クラスタリング」は、「多重共クラスタリング」の下位概念であるものとする。

なお、この実施の形態においては、単に「クラスタリング」という用語は、１つのビューにおけるクラスタの組を生成することを言うものとし、たとえば、図２２Ａのように特徴量のビューへの分割と対象のクラスタリングを行う場合は「多重クラスタリング」と呼び、図２２Ｂのようにビューへの分割と共クラスタリングとを同時に実行する場合は、「多重共クラスタリング」と呼ぶことで区別する。

図２３は、「多重共クラスタリング」において、１つのビューの中に、異なる種類の確率分布の確率モデルが想定されている場合を示す概念図である。

図２３においては、白のブロックと、網掛のブロックでは、異なる種類の確率分布の確率モデルに従っている場合を示している。

たとえば、白のブロックは、連続的な確率変数の確率分布であるのに対して、網掛の部分は、離散的な確率変数の確率分布を想定している場合などを示す。

以後に説明するように、本実施の形態の「多重共クラスタリングの学習方法」では、このように、異なる分布を含む分布族に対して、クラスタリング処理を実行することが可能である。

図２４は、多重共クラスタリングの学習方法の概要を説明するためのフローチャートである。

多重共クラスタリングの学習方法の処理が開始されると（Ｓ３００）、クラスタリング分類器生成部３０１０は、データ行列について、ランダムに特徴量を部分群に分割して、特徴量のビューとビュー内の特徴量クラスタを生成する（Ｓ３０２：後述するＹの生成（Ｙの初期化）に対応）。

続いて、クラスタリング分類器生成部３０１０は、ステップＳ３０２で生成された特徴量のビューおよび特徴量クラスタに対応して、対象クラスタの分割を生成し最適化する（ステップＳ３０４：後述するＺの生成に対応）。

さらに、クラスタリング分類器生成部３０１０は、得られた対象クラスタに対して、特徴量の分割を最適化する（Ｓ３０６：後述するＹの生成処理に対応し、生成されたＺを使ってＹを最適化する）。

続いて、クラスタリング分類器生成部３０１０は、目的関数が所定の条件を満たして収束しているかどうかを判定する（Ｓ３０８）。なお、この目的関数は、後述するような関数Ｌ(q(φ))に相当する。関数Ｌ(q(φ))は、これも後述するＹおよびＺをアップデートするに従って単調に増加する性質があり、その増え方が十分小さくなったと判断された時に収束したものと判定することになる。クラスタリング分類器生成部３０１０は、収束していなければ（Ｓ３０８でＮＯ）、処理をステップＳ３０４に復帰させ、収束していれば（Ｓ３０８でＹＥＳ）、処理を次のステップに進める。

そして、クラスタリング分類器生成部３０１０は、目的関数の大きさを記憶装置２０８０に格納する（Ｓ３１０）。

次に、クラスタリング分類器生成部３０１０は、ステップＳ３０２～Ｓ３１０までの処理を、所定回数実施したかを判定する。クラスタリング分類器生成部３０１０は、所定回数実施していない場合は（Ｓ３１２でＮＯ）、処理をステップＳ３０２に復帰させ、所定回数実施している場合は（Ｓ３１２でＹＥＳ）、処理を次のステップに進める。

クラスタリング分類器生成部３０１０は、目的関数が最大となる特徴量分割、クラスタの分け方を最終結果として（Ｓ３１４）、多重共クラスタリングについての学習の処理を終了させ、クラスタリング分類器を生成する。

図４８は、クラスタリング分類器の構成を示す概念図である。

図４８では、３つのビューに分割が実行され、ビュー１は特徴量群１に対応して分割されて、クラスター１～３に対象がクラスタリングされており、ビュー２は特徴量群２に対応して分割されて、クラスター１～４に対象がクラスタリングされており、ビュー３は特徴量群３に対応して分割されて、クラスター１～２に対象がクラスタリングされているものとする。

したがって、クラスタリング分類器生成部３０１０は、記憶装置２０８０の疾患識別器データ３１１２に、各ビューに対応する特徴量の情報、ビューごとの確率密度関数を特定するための情報（たとえば、確率密度関数が正規分布であるときは、分布の中心座標μと分散σ²）を格納する。

図４８に示すように、学習データ中に含まれていない「新規なデータ（新規被験者のデータ）」であって、特徴量群１～３を有する「新規なデータ」が、クラスタリング分類器に入力されると、判別値算出部３０１２は、ビュー１に対応する特徴量群１について、ビュー１内の確率密度関数に基づいて、各クラスタに属する事後確率を計算する。図４８では、ビュー１のクラスタ２に属する事後確率が最も高く、新規データは、ビュー１のクラスター２に属するとの分類の結果が出力される。同様にして、新規データは、ビュー２では、クラスター３に、ビュー３では、クラスター１に属する、という分類の結果が、クラスタリング分類器から出力されることになる。

（多重共クラスタリングの処理の詳細）
以下では、図２４で説明した多重共クラスタリングの学習方法について、さらに詳しく説明する。
なお、多重共クラスタリングの処理の詳細については、下記の公知文献１７に開示があるので、以下では、その概要について説明する。公知文献１７の記載の全体をここに参照により援用する。
公知文献１７：Tomoki Tokuda, Junichiro Yoshimoto, Yu Shimizu, Go Okada, Masahiro Takamura, Yasumasa Okamoto, Shigeto Yamawaki, Kenji Doya，“Multiple co-clustering based on nonparametric mixture models with heterogeneous marginal distributions”, PLOS ONE | https://doi.org/10.1371/journal.pone.0186566 October 19, 2017

図２５は、図２４の多重共クラスタリングの学習方法におけるベイズ推定のグラフ表現を示す図である。

多重共クラスタリングモデルは、図２５のグラフィカルモデルにまとめられており、このグラフィカルモデルは、関連するパラメータとデータ行列間の因果関係のリンクを明確にしている。

（多重共クラスタリングモデル）
特徴量（脳機能結合値）と被験者（ここでは、患者群の被験者）とは、図２１Ａに示したようなデータ行列として表現されるものとする。
そして、データ行列Ｘは、事前に知られているＭ個の分布からなる分布族で構成されていると仮定する。
分布族に属する確率分布としては、ガウス分布、ポアソン分布、およびカテゴリ分布/多項分布などが含まれ得るものとする。

クラスタリング分類器生成部３０１０は、Ｘ^(m)に対して、各データサイズがｎ× ｄ^(m)であるように、以下のように分割する。
Ｘ＝ {Ｘ⁽¹⁾，…，Ｘ^(m)，…，Ｘ^(M)}
ここで、ｍは、分布族（ｍ＝１，…，Ｍ）を示す指標である。さらに、ビュー（視点）の数をＶ(すべての分布族に共通)とし、ビューｖおよび分布族ｍの特徴量クラスタの数をＧ_ν ^(m)とし、ビューｖのオブジェクトクラスタの数をＫ_v(すべての分布族に共通) で示す。

さらに、表記の簡単化のために、特徴の数およびクラスタの数を示すために、空のクラスタが存在することを許して、Ｇ^(m)＝max_v Ｇ_v ^(m) およびＫ＝max_v Ｋ_vと表記する。
この表記では、分布族ｍに対する独立同分布（i.i.d.）のｄ^(m)次元ランダムベクトルＸ₁ ^(m), …, Ｘ_n ^(m)について、ｄ ^(m)× Ｖ × Ｇ^(m)の(３階の)特徴量分割テンソルＹ^(m)を考えることとし、分布族ｍの特徴量ｊがビューｖの特徴量クラスタｇに属するときは、Ｙ_j,v,g ^(m)＝１（それ以外は０）であるものとする。

異なる分布族に対してこれを組み合わせて、Ｙ= {Ｙ^(m)}_mとする。

同様にして、ｎ×Ｖ×Ｋのオブジェクト分割 (３階) テンソルＺを考え、オブジェクトｉがビューｖのオブジェクトクラスタｋに属する場合はＺ_i、v、k = 1であるとする。

特徴量ｊはビューのうちの１つ( Σ_v,gＹ_{j, v, g} ^(m) ＝１）に属し、オブジェクトｉは各ビューに属している (すなわち、Σ_kＺ_i,v,k ^(m)＝１)。さらに、Ｚはすべての分布族に共通しており、これは、推定される確率モデルは、すべての分布族の情報を使用して被験者クラスタリング解を推定していることを意味する。

まず、図２５に示されるように、Ｙの事前生成モデルに対しては、ビューと特徴量クラスタの階層構造を考え、最初にビューが生成され、次に、特徴量クラスタが生成されるものとする。したがって、特徴量はビューと特徴量クラスタのペアのメンバーシップで分割されることとなり、特徴量の分割の割り当ては、ビューと特徴量クラスタで共同的に決定される。

一方、図２５に示されるように、オブジェクト（対象）は各ビューのオブジェクトクラスタに分割されるため、Ｚに対して、オブジェクトクラスタの１つの構造だけを考えている。これらの生成モデルは、以下に説明するように、すべて「棒折り過程」（ＳＢＰ：Stick Breaking Process，スティックブレイキング過程）に基づいていると仮定する。

（特徴量クラスタＹの生成モデル）
Ｙ_j.. ^(m)が、階層棒折りプロセスによって生成される、分布族ｍの特徴量ｊのビュー／特徴量クラスタメンバーシップベクトルを示すものとすると、以下の式が成り立つ。

ここで、τ^(m) は、1 × ＧＶベクトル (τ_{1 ,1} ^(m) ,…,τ_G,V ^(m))^T を示す。
Ｍｕｌ(・|π)は確率パラメータπを持つ１つのサンプルサイズの多項分布である。
Beta(・|ａ，ｂ)は、事前サンプルサイズ（ａ，ｂ)を有するベータ分布である。
Ｙ_j.. ^(m)は１×ＧＶベクトル(Ｙ_j,1,1 ^(m), … ,Ｙ _j,V,G ^(m))^Tを示す。

ここでは、所定の条件に従って、十分な大きなＶのビューの数とＧの特徴量クラスタの数を切り捨てる。

公知文献１８：Blei DM, Jordan MI, et al. Variational inference for Dirichlet process mixtures. Bayesian analysis. 2006; 1(1) : 121－143，
https://doi.org/10.1214/06-BA104
Ｙ_j,v,g ^(m)＝１の場合、特徴量ｊはビューｖの特徴量クラスタｇに属する。デフォルトでは、ハイパーパラメータである集中度パラメータα₁とα₂は、１に設定される。

（オブジェクトクラスタＺの生成モデル）
ビューｖのオブジェクトｉの被験者クラスタメンバーシップベクトルであって、Ｚ_{i, v.}と表記されるベクトルが以下の式によって生成される。

ここでＺ_{i, v.}は、Ｚ_i,v= （Ｚ_{i, v, 1}，…，Ｚ_{i, v, K}）^Tによって与えられる１×Ｋ（Ｋは十分に大きい値を取る)ベクトルである。集中度パラメータβは、１に設定される。

（尤度と事前分布）
各インスタンスＸ_i,j ^(m) は、ＹとＺに条件付きで、独立に特定の分布に従うと仮定する。ビューｖ、特徴量クラスタｇ、およびオブジェクトクラスタｋのクラスタブロックにおける分布族ｍのパラメータを θ_v,g,k ^(m)と表す。

さらに、Θ={θ_v,g,k ^(m)}_v,g,k,mと表記して、Ｘの尤度の対数は以下の式に従う：

ここで、Ｉ（ｘ）は指標関数で、ｘが真の場合は１を返し、それ以外の場合は０を返す。尤度は、ｗ＝｛ｗ_v｝_v、ｗ´＝｛ｗ´_g,v ^(m)｝_g,vおよびｕ＝｛ｕ_k,v｝_k,vに直接関連付けられていない。

未知変数の結合事前分布φ＝｛Ｙ，Ｚ，ｗ，ｗ´，ｕ，Θ｝（すなわち、クラスのメンバーシップ変数とモデルパラメータ）は、次のように与えられる。

（変分推定）
ＹとＺのＭＡＰ（最大事後）推定には変分ベイズＥＭアルゴリズムを使用する。
このような変分ベイズＥＭアルゴリズムについては、以下の公知文献１９に開示がある。公知文献１９の記載の全体をここに参照により援用する。
公知文献１９：Guan Y, Dy JG, Niu D, Ghahramani Z. Variational inference for nonparametric multiple clustering. In:MultiClust Workshop, KDD-2010; 2010.

対数周辺尤度ｐ（Ｘ）はジェンセンの不等式を使用して以下のように近似される。

なお、ジェンセンの不等式については、以下の公知文献２０に開示がある。公知文献２０の記載の全体をここに参照により援用する。

公知文献２０：Jensen V. Sur les fonctions convexes et les inegalites entre les valeurs moyennes. Acta Mathematica. 1906; 30(1):175－193.
https://doi.org/10.1007/BF02418571
ここで、q(φ)はパラメータφの任意の分布である。q(φ)とp(φ)の間のカルバック・ライブラー・ダイバージェンス、すなわち、ＫＬ(q(φ), p(φ|Ｘ))によって左辺と右辺の差が与えられることが証明される。したがって、q(φ)を選択するアプローチは、ＫＬ(q(φ), p(φ|Ｘ))を最小化することであり、これは評価するのは、一般には難しい。

ここでは、異なるパラメータ(平均場近似)に対して因数分解されるq(φ)を選ぶ。

ここで、各 q(・)は、パラメータのサブセットｗ_v, ｗ´_g,v ^(m),Ｙ_j.. ^(m), ｕ_{k, v}, Ｚ_i,v.およびθ_v,g,k ^(m)に対して、さらに因数分解される。

一般に、ＫＬ（Π_l=1 ^L q_l(φ_l), p (φ|ｘ)）を最小化する分布 q_i(φ_i)は、以下の式によって与えられる。

この性質については、以下の公知文献２１に開示がある。公知文献２１の記載の全体をここに参照により援用する。
公知文献２１：Murphy K. Machine Learning: A Probabilistic Perspective. Cambridge, Massachusetts: MIT Press;2012.

この性質を現在検討しているモデルに適用すると、以下を示すことができる。

ここで、以下の形で表される関数を考える。

上記の式を除くハイパーパラメータは、以下の式で表される。

ψ(・)は、ガンマ関数の対数の１次導関数として定義されるダイガンマ関数を示す。
τ_j,g,v ^(m)は、各ペア（ｊ，ｍ）についての複数のペア（ｇ，ｖ）にわたって正規化される。一方、η_{i, v, k}は各(i, v) のペアでｋについて正規化される。
観測モデルとパラメータΘの事前分布は、後述する。

（観測モデル）
観測モデルでは、ガウス分布、ポアソン分布、およびカテゴリ分布/多項分布を考慮する。クラスタブロックごとに、クラスタブロック内の特徴が独立していることを前提に、これらのファミリーの１変量分布をフィッティングさせる。これらの分布族のパラメータは、共役事前分布を仮定する。

（最適化アルゴリズム）
ハイパーパラメータの更新方程式では、変分ベイズＥＭアルゴリズムを用いて、次のように計算が実行する。

まず、ランダムに{τ^(m)}_mと{η_v}_vを初期化し、式（１）の下限Ｌ(q(φ)) が収束するまでハイパーパラメータを更新する。これは、Ｌ(q(φ))の観点から局所的に最適な分布q(φ)を生成する。この手順を何度も繰り返し、近似事後分布q*(φ)として最大の下限を持つ最良の解を選択する。
ＹとＺのＭＡＰ推定値は、それぞれ argmax_Y q^* _Y(Y)およびargmax_Z q^* _Z(Z)として評価される。

下限Ｌ(q(φ))は以下の式で与えられる。

右辺の両方の項は、閉じた形式で導出できる。q(φ)が最適化されるにつれて、この値は、単調に増加することが示される。すなわち、上述した通り、関数Ｌ(q(φ))は、ＹおよびＺをアップデートするに従って単調に増加する性質があり、その増え方が十分小さくなったと判断された時（特に限定されないが、たとえば、増分が所定値以下となるなどの条件が成り立つ時）に収束したものと判定することになる。

まず、各特徴の分布族が特定され、対応する分布族のデータ行列を生成する。次に、データ行列のセットに対して、さらに、ＹとＺのＭＡＰ推定値を生成し、ＹおよびＺの推定値を使用して、各ビューにおけるオブジェクト／特徴量クラスタが分析される。

（モデル表現）
多重共クラスタモデルは、ビューの数と特徴/オブジェクトクラスタの数がデータ駆動のアプローチで導き出されるため、さまざまなクラスタリングモデルを表現するのに十分な柔軟性を備えている。たとえば、ビューの数が１の場合、モデルは、共クラスタモデルと一致する。特徴量クラスタの数がすべてのビューに対して１つである場合、多重クラスタリングモデルと一致する。さらに、ビューの数が1であり、特徴量クラスタの数が特徴の数と同じである場合、独立した特徴を持つ従来の混合モデルと一致する。さらに、このモデルは、オブジェクトクラスタを区別しない無情報特徴を検出することができる。このような場合、モデルはオブジェクトクラスタの数が１つのビューを生成する。モデルの利点は、このような基になるデータ構造を自動的に検出することである。

以上のような「多重共クラスタリング法」により、以下のことが可能となる。
１）データの背後にある複数のクラスターの分け方（対象の分け方だけでなく、特徴量の分け方も含める）とそれに対応する特徴量群を、データ駆動により同定することができる。
２）この手法によって他の手法では見出せなかったクラスターを同定することが可能となる。
３）さらに、それぞれのクラスターの分け方を特徴量によって意味づけることができ、個々のクラスターの解釈を容易にできる。

［データセットに対するクラスタリングの結果の評価］
以下では、上述したようなＳＲＰＢＳとして公開されている多数の被検者から収集された多施設大規模ｆＭＲＩデータを２つに分割して、それぞれを上述したような多重共クラスタリング法で使用することで、クラスタリングについての汎化性能を検証する。

図２６Ａおよび図２６Ｂは、それぞれ、このように２つに分割されたデータセット１とデータセット２とを示す図である。

図２６Ａに示すデータセット１は、施設１～４までで取得された健常者５４５人、うつ患者１３８人のデータからなり、図２６Ｂに示すデータセット２は、施設５～８までで取得された健常者２６３人、うつ患者１８１人のデータからなる。基本的には、図１０および図１１に示したデータセットに対応する。

図２７は、各データセットにおいてクラスタリングを実施する概念を説明する概念図である。
図２７に示すように、データセット１に対しては、それぞれ、独立に、図２４に示したフローに従って多重共クラスタリング法により、クラスタリングが実行される。

ここで、関心の対象となるのは、このようにして、データセット１と、データセット２に対して、それぞれ個別に実行されたデータ駆動によるクラスタリング法により得られた各クラスタが、相互にどの程度、類似しているか（一致の程度は、どの程度か）という問題である。

データセット１と、データセット２におけるクラスタリングが、同一または類似する特性を有するクラスタ（被験者グループ）へ分類（グループ分け）できているのであれば、このようなデータ駆動によるクラスタリングは、施設や計測装置等に依存することなく、汎化性の高い状態で実行されていることになる。そこで、問題となるのは、「同一または類似する特性を有するクラスタへ分類」されていることを、どのように定量的に評価するのかということになる。

図２８は、被験者データに対する多重共クラスタリングの例を示す概念図である。
図２８（ａ）に示すように、入力となるデータ行列として、行方向に被験者が並び、列方向に特徴量が並んでいるものとする。

この入力のデータ行列に対して、多重共クラスタリングを実行すると、たとえば、図２８（ｂ）に示すように、２つのビューに特徴量が分割され、各ビューにおいて、被験者がクラスタリングされることになる。

図２９は、データセット１およびデータセット２に対して、実際に多重共クラスタリング処理を実施した結果を示す図である。
図２９においては、クラスタリングのための特徴量としては、データセット１およびデータセット２とも、それぞれ９９個が選択されている。
その上で、データセット１については、１３８人のうつ患者に対して、データセット２については、１８１人のうつ患者に対して、多重共クラスタリング処理が実行されている。

データセット１については、ビュー１とビュー２の２つのビューに特徴量が分割される。ビュー１については、さらに２つの特徴量クラスタに共クラスタリングされ、被験者は５つの被験者クラスタに分割され、ビュー２についても、被験者は５つのクラスタに分割される。

データセット２についても、ビュー１とビュー２の２つのビューに特徴量が分割される。ビュー１については、さらに２つの特徴量クラスタに共クラスタリングされて、被験者は４つの被験者クラスタに分割され、ビュー２については、被験者は５つのクラスタに分割される。

図３０は、データセット１とデータセット２において、それぞれのビューに割り当てられた脳機能結合（ＦＣ）の個数を示す表である。
データセット１では、ビュー１に特徴量としてＦＣが９２個、ビュー２にＦＣが７個割り当てられる。
データセット２では、ビュー１に特徴量としてＦＣが９３個、ビュー２にＦＣが６個割り当てられる。

また、この表では、データセット１とデータセット２で、割り当てられた脳機能結合のうち、一致する個数を表の対角上に記載している。データセット１とデータセット２で、ビュー１とビュー２に割り当てられた脳機能結合は、ほぼ一致していることがわかる。
（クラスタリング（層別化）の汎化性（データセット間類似性）検証方法）
以下では、データセット１と、データセット２に対して、それぞれ個別に実行されたデータ駆動によるクラスタリングにより得られた各クラスタが、相互にどの程度、類似しているか（一致の程度は、どの程度か）を定量的に評価する。

図３１は、このようなクラスタリングの類似度（層別化の汎化性能）の評価手法を説明するための概念図である。

まず、図３１（Ａ）に示すように、データセット１とデータセット２を独立に、上述したような多重共クラスタリング法により、クラスタに分割した場合、被験者は、それぞれのデータセットのクラスタにおいて、相互に独立であるために、クラスタリングの類似度を比較する、ということは困難である。

ここで、データセット１で作成された分類器１で、データセット１の被験者の分類を実行した結果をクラスタリング１とする。一方で、データセット２で作成された分類器２で、データセット２の被験者の分類を実行した結果をクラスタリング２とする。

これに対して、図３１（Ｂ）に示すように、データセット１で作成された分類器１で、データセット２の被験者の分類を実行した結果をクラスタリング１´とする。一方で、データセット２で作成された分類器２で、データセット１の被験者の分類を実行した結果をクラスタリング２´とする。

この場合、クラスタリング１とクラスタリング１´の間、および、クラスタリング２とクラスタリング２´の間では、それぞれ、共通な被験者について分類を行っているので、それぞれの類似度を評価することができる。

（クラスタリング間の類似度（再現度）を測る評価尺度）
ここで、データ駆動によりクラスタリング処理を実行しており、図２９におけるクラスタインデックスの値自体（インデックスの値の順序）には、意義はないので、同じデータ集合に対して、異なるクラスタリングが実行された場合に、その類似度をどのように評価するかが問題になる。

たとえば、同じデータ集合Ｘに対する２つのクラスタリング結果πとρとがある場合において、この２つのクラスタリング結果の類似性（外的妥当性尺度）を評価する尺度として、ランド指数（Rand index）が知られている。

データ集合中のすべてのデータの対（ペア）｛ｘ₁，ｘ₂｝∈Ｘ（Ｍ＝Ｎ（Ｎ－１）／２個）については、対の種類として、以下のものがあり、それぞれの種類に属するペアの個数を以下のように定義する。

このとき、２つのクラスタリングにより分類されるクラスタが同じクラスタかどうかの判定の正解率として、ランド指数が以下の式により定義される。

ただし、たとえば、データ集合の各クラスタの要素数に偏りがある場合などでは、「ランダムにクラスタリングした場合でも、ランド指数が高い値となってしまう」場合があることが知られている。このため、より厳密には、以下のような調整ランド指数（ＡＲＩ：Adjusted Rand Index）が使用される。

図３２Ａおよび図３２Ｂは、ＡＲＩについて説明するための概念図である。
なお、ＡＲＩについては、たとえば、以下の公知文献２２に開示がある。公知文献２２の記載の全体をここに参照により援用する。
公知文献２２：Jorge M. Santos and Mark Embrechts，“On the Use of the Adjusted Rand Index as a Metric for Evaluating Supervised Classification”，ICANN 2009, Part II, LNCS 5769, pp. 175-184, 2009.

上述のとおり、同じデータ集合に対して、２つのクラスタリング結果をそれぞれ適用する場合、図３２Ａに示すように、２回とも同じクラスタに分類される場合と、２回とも別のクラスター分類される場合がある一方、図３２Ｂに示すように、１回は同じクラスタに分類されるが、もう１回は異なるクラスタに分類される場合も存在する。

ＡＲＩは、２つのクラスタリングが互いに独立だとした場合に、データ対がともに２つのクラスタリングで同一のクラスタに分類され、または、ともに異なるクラスタに分類される場合の期待値を計算して、その期待値をランド指数の分子と分母のそれぞれから減算して、計算される。したがって、ＡＲＩでは、クラスタリングに相関がない場合に、その値が０となるように調整されていることになる。

ここでは、Ａは、２つのクラスタリングに対して、［（２回とも同じクラスタに分類）＋（２回とも別のクラスタに分類）］された被検者ペアの数を表し、max（Ａ）は、全ペアの個数を表し、Ｅは、２回のクラスタリングが独立にもかかわらず、割り当て結果が一致する被検者ペアの数を表す。

図３３Ａおよび図３３Ｂは、それぞれ、クラスタリング１とクラスタリング１´の間、および、クラスタリング２とクラスタリング２´の間の類似度の評価結果を示す図である。

図３３Ａは、データセット１とデータセット２について、それぞれのビューに対するＡＲＩを算出した表である。
データセット１とデータセット２について、ビュー１については、ＡＲＩ＝０．４７、ビュー２については、ＡＲＩ＝０．５１となっており、有意な類似性があるといえる。

図３３Ｂは、図３３Ａに対応するパーミュテーションテストの結果を示す。ビュー１およびビュー２については、要素の交換を行ったとき（ヒストグラムで表示）と比較して、ビュー１とビュー２では、ＡＲＩ値（実線で示す）が統計的に有意に高い値となっていることがわかる。

なお、「パーミュテーションテスト」とは、被験者のクラスタ属性ラベルを被験者間でランダムに交換した場合においてＡＲＩ値を算出した結果であり、図中では、このような交換を所定回数実施したときのその分布がヒストグラムとして表示されている。もしも、クラスタリング間の類似が統計的に有意であるならば、要素をランダムに交換した場合に比べて、比較対象となるクラスタリング間のＡＲＩ値が有意に高い値となることになる。

以上より、データセット間でのクラスタリング（層別化）について、有意な類似性が認められた、すなわち、汎化したクラスタリングが実現されたと判断できる。

以上説明した通り、データセット１に対する多重共クラスタリングも、データセット２に対する多重共クラスタリングも、ともに、データ駆動により実現されたものであるから、脳機能結合を特徴量とする「患者の層別化」の基礎となるといえる。

図３４は、クラスタリング１とクラスタリング１´において、それぞれビュー１の各クラスタに割り当てられる被験者数の分布を示す表である。

被験者クラスタインデックスを適切に並べ替えることにより、ほとんどの被験者が、表の対角の近傍に分布するようにでき、２つのクラスタリングが、相互に類似していることが視覚的にも確認できる。

［ハーモナイゼーション処理］
以下では、以下の文献に開示され、図６において、ハーモナイゼーション処理と呼んでいる処理の内容について、説明する。

［トラベリングサブジェクト法のハーモナイゼーション］
以下では、以上説明してきた「疾患識別器」の生成や、層別化のための「クラスタリング処理」のために使用される手法であって、標本バイアスとは独立に、測定バイアスを評価して計測データをハーモナイズするための手法について説明する。

図３５は、本実施の形態のｒｓ－ｆｃＭＲＩ法において、サイト間を移動しながら計測を受ける移動被験者（以下、「トラベリングサブジェクト：旅行被験者」）によるサイト間差の評価手法を説明するための概念図である。

以下に説明するように、本実施の形態では、トラベリングサブジェクトのデータセットを使用して、測定バイアスだけを除くことができるようなハーモナイゼーション法について説明する。

図３５を参照して、計測サイトＭＳ．１～Ｍｓ．Ｎｓのサイトを横断して測定バイアスを評価するために、トラベリングサブジェクトＴＳ１（人数：Ｎts人）のデータセットを取得する。

健康なＮts人の参加者の安静時脳活動が、Ｎｓ個のサイトの各々で撮像されるものとして、Ｎｓ個のサイトは、患者データを撮像したすべてのサイトを含んでいるものとする。

取得されたトラベリングサブジェクトのデータセットは、データセンター２００における記憶装置２１０に移動被験者データとして格納される。

そして、後述するように、計算処理システム３００において、「脳活動バイオマーカーのハーモナイズ方法」のための処理が実行される。

トラベリングサブジェクトのデータセットは健常者群だけを含む。また、参加者はすべてのサイトにわたって同じとする。したがって、トラベリングサブジェクトについては、サイト間差は「測定バイアス」のみから成る。

以下に説明する本実施の形態のハーモナイゼーション法では、「脳活動バイオマーカーのハーモナイズ方法」として、各計測サイトにおける計測データについて、「測定バイアス」の影響を除いて補正する処理を行う。

つまり、以下では、「測定バイアス」および「標本バイアス」を、「統計モデリング」手法のなかの「一般化線形混合モデル（ＧＬＭＭ：Generalized Linear Mixed Model）」を用いて評価する。

ここで、通常は、ＧＬＭ（Generalized Linear Model:一般化線形モデル）は、「応答変数」の確率分布を説明する「説明変数」を組み込んだモデルである。ＧＬＭには、「確率分布」「リンク関数」「線形予測子」という主に３つの部品があり、この部品の組み合わせかたを指定することによって、さまざまなタイプのデータを表現できる。

さらに、ＧＬＭＭ（一般化線形混合モデル）は、ＧＬＭでは説明できない「人間が測定できない・測定しなかった個体差など」を組み込むことができる統計モデルである。ＧＬＭＭでは、たとえば、対象がいくつかの部分集合（たとえば、計測場所が異なる部分集合）からなるような場合も、その場所差をモデルに組み込むことができる。別の言い方をすれば、複数の確率分布を部品とする（混合させた）モデルということになる。

たとえば、ＧＬＭＭについては、以下の公知文献２３に開示がある。公知文献２３の記載の全体をここに参照により援用する。
公知文献２３：久保拓弥著、「データ解析のための統計モデリング入門」、岩波書店、2012年第1刷、2017年第14刷

ただし、以下に説明する本実施の形態の統計モデルでは、通常、「効果」と呼ばれるものについて、「バイアス」と「要因（ファクター）」という用語を、「バイアス」に関しては「測定バイアス」と「標本バイアス」に使用し、「要因（ファクター）」に関しては、その他の要因（被験者要因や疾患要因）に関して使用する。

そして、以下の解析は、単純なＧＬＭＭの流れとは異なり、「固定効果」と「ランダム効果」という分け方をせずに要因の解析を行う。これは、通常、ＧＬＭＭを使用すると、ランダム効果に関しては分散のみが推定され、各要因の効果の大きさがわからなくなるからである。そこで、以下では、各要因の効果の大きさを評価するために、各要因について、平均０の固定効果となるように、以下のように変数を変換して、推定を行う。

ｉ）各々の機能的なコネクティビティに対する相関値のすべてのサイトにわたる平均からの偏差として、各サイトの測定バイアスを定義する。
ｉｉ）健常者および精神疾患の患者の標本バイアスは、互いに異なると仮定する。したがって、健常者群および各疾患を有する患者群に対して個別に、各サイトの標本バイアスを計算する。
ｉｉｉ）疾患要因は、健常群の値からの偏差として定義する。

すなわち、以下では、患者が含まれるデータセットと、トラベリングサブジェクトのデータセットに対して、一般化線形混合効果モデルを以下のように適用する。
トラベリングサブジェクトが、Ｎts人であり、Ｎs個の計測サイトのうち、健常者の計測が行われたサイト数をＮsh個とし、ある疾患（ここでは、添え字“dis”で表現する）の患者の計測が行われたサイトの個数をＮsd個とする。
参加者要因（ｐ）、測定バイアス（ｍ）、標本バイアス（Ｓhc, Ｓdis）および精神疾患要因（ｄ）は、患者に対する計測結果のデータセットおよびトラベリングサブジェクトのデータセットからのすべての参加者の機能的なコネクティビティの相関値に、回帰モデルを適合させることにより評価される。

以下では、ベクトルは小文字（例えば、ｍ）によって表示され、ベクトルはすべて列ベクトルであると仮定される。
ベクトルの要素は、ｍ_kのように添字によって表示される。

脳領野間のｎ個の相関値からなる機能的コネクティビティベクトル（列ベクトルとする）の回帰モデルは、以下の式のように表される。

参加者の特性を表わすために、1-of-Kのバイナリーコード体系を使用し、サイトｋに属する測定バイアスｍに対するターゲットベクトル（たとえば、ｘm）は、１に等しい要素ｋ以外は、全てゼロに等しい。
参加者がどんなクラス（健常者、患者、トラベリングサブジェクト）にも属さなければ、ターゲットベクトルは、すべての要素が０に等しいベクトルである。
上付きの添字Tは、マトリックスかベクトルの置換を表示し、ｘ^Tは行ベクトルを表わす。

ここで、ｍは、測定バイアス（Ｎs×１の列ベクトル）を表し、ｓhcは、健常群の標本バイアス（Ｎsh×１の列ベクトル）を表し、ｓdisは、患者の標本バイアス（Ｎsd×１の列ベクトル）を表し、ｄは、疾患要因（２×１の列ベクトル、健常および疾患を要素とする）、ｐは、参加者要因（Ｎts×１の列ベクトル）を表し、ｃｏｎｓｔは、すべての計測サイトからの全ての参加者（健常者、患者、トラベリングサブジェクトを含む）にわたる機能的コネクティビティの平均を表し、ｅ～Ｎ（０，γ^-1）は、ノイズを表す。

なお、ここでは、説明の簡略化のための疾患の種類は１種類であるものとして説明している。疾患の種類が複数ある場合については、後述する。

各機能的なコネクティビティの相関値に対しては、回帰モデルの計画行列がランク不足であるので、Ｌ２正規化による最小二乗回帰を使用して、それぞれのパラメータを評価する。なお、Ｌ２正規化による最小二乗回帰法以外にも、たとえば、ベイズ推定法など別の評価方法を用いることも可能である。

上記のような回帰計算の後には、被験者ａのｂ番目のコネクティビティは、以下のように記述できる:

図３６は、被験者ａのｂ番目の機能的なコネクティビティの表現を説明するための概念図である。

図３６においては、１項目と２項目のターゲットベクトルの意義と、測定バイアスベクトルおよび健常者の標本バイアスベクトルについて示している。
第３項目以降についても同様である。

（ハーモナイゼーション処理のフロー）
図３７は、ハーモナイゼーションのために、測定バイアスを算出する処理を説明するためのフローチャートである。

まず、データセンター２００の記憶装置２１０には、各計測サイトから被験者（健常者および患者）のｆＭＲＩ計測データ、被験者の属性データ、測定パラメータが収集される（図３７Ｓ４０２）。

続いて、トラベリングサブジェクトＴＳ１の脳活動が、特に限定されないが、たとえば、所定の周期（たとえば、１年周期）で、各計測サイトを巡回して計測され、データセンター２００の記憶装置２１０には、各計測サイトからトラベリングサブジェクトのｆＭＲＩ計測データ、被験者の属性データ、測定パラメータが収集される（図３７Ｓ４０４）。

ハーモナイゼーション算出部３０２０は、上述したようなＧＬＭＭ（一般化線形混合モデル）を用いることで、機能的コネクティビティについて、各計測サイトの測定バイアスを評価する（図３７Ｓ４０６）。

このようにして算出された各計測サイトの測定バイアスを、ハーモナイゼーション算出部３０２０は、記憶装置２０８０に測定バイアスデータ３１０８として格納する（図３７Ｓ４０８）。

（判別器生成処理におけるハーモナイゼーション）
判別処理部３０００が、被験者に対する疾患あるいは健常ラベルについての疾患識別器を生成する処理における脳機能結合値に対するハーモナイゼーションについて簡単に説明する。

このような疾患識別器は、被験者を診断するための補助情報（支援情報）を提供する。

相関値補正処理部３００４は、記憶装置２０８０に格納されている各計測サイトごとの測定バイアスデータ３１０８を読出し、疾患識別器生成の機械学習のための訓練対象となる各被験者の相関行列の非対角成分について、以下の式のようにして、ハーモナイゼーション処理を実行する。

ここで、機能的コネクティビティConnectivityは、ハーモナイゼーション前の機能的コネクティビティベクトルを表し、Ｃsubは、ハーモナイゼーション後の機能的コネクティビティベクトルを表す。また、ｍ（ハット）（以下、文字ｘの頭部に＾が付されたものを“ｘ（ハット）”と表記する）は、上述したようなＬ２正規化による最小二乗回帰により評価された計測サイトにおける測定バイアスを表す。これにより、機能的コネクティビティConnectivityが計測された計測サイトに対応する測定バイアスが、機能的コネクティビティConnectivityから減算され、ハーモナイゼーション処理を受けることになる。
補正処理が実行された後のデータは、補正後相関値データ３１１０として、記憶装置２０８０に格納される。

なお、以下の公知文献２４にも記載されるように、疾患要因は、脳全体のコネクティビティに関係しているわけではなく、コネクティビティの特定の部分集合に関係している、と想定するのが妥当である。公知文献２４の記載の全体をここに参照により援用する。
公知文献２４：Yahata N, et al. A small number of abnormal brain connections predicts adult autism spectrum disorder. Nat Commun 7, 11254 (2016).

そこで、続いて、疾患識別器生成部３００８により、被験者についての疾患／健常ラベルを含む被験者の疾患ラベルと補正処理後の機能的なコネクティビティとについて、上述した通り、特徴選択を伴う学習処理により疾患識別器が生成される。

なお、過学習を抑制して、特徴選択とモデル化を実行するための手法としては、ＬＡＳＳＯによる正則化ロジスティック回帰に限定されるものではなく、たとえば、以下の公知文献２５に開示されるようなスパースロジスティック回帰や、他のスパースベイズ推定法など、他の手法を使用してもよい。以下の公知文献２５の記載の全体をここに参照により援用する。
公知文献２５：Okito Yamashita, Masaaki Sato, Taku Yoshioka, Frank Tong, and Yukiyasu Kamitani. “Sparse Estimation automatically selects voxels relevant for the decoding of fMRI activity patterns.” NeuroImage, Vol. 42, No. 4, pp. 1414-1429, 2008.

以上により、疾患識別器（脳活動バイオマーカによる診断のための判別器として機能する）が生成され、疾患識別器生成部３００８は、判別器を特定するための情報を記憶装置２０８０に、疾患識別器データ３１１２として格納する。

以上のような処理は、上述したように、特に限定されないが、所定の期間（たとえば、１年）ごとに実施されるものとしてよい。

そして、いずれかの計測サイトに、新たに１人の被験者についての機能的コネクティビティの計測が実行された場合、当該計測サイトの測定バイアスは、所定期間の間は一定と仮定できるので、判別値算出部３０１２は、入力される当該被験者についての入力データの相関行列の要素の値から、上述した手続きですでに算出されている当該入力データが計測された計測サイトに対応する「測定バイアス」を減算することにより、ハーモナイゼーション処理を実行する。そして、上述した手続きにより、すでに作成されている「疾患識別器」により、当該被験者についての判別ラベルを判別結果として出力する。

判別結果は、「疾患」「健常」の二者のいずれかを示す値であってもよいし、あるいは、「疾患」「健常」の少なくとも一方の確率を示す値であってもよい。

なお、判別値算出部３０１２が実施する判別処理に際して、入力される「入力データ」としては、被験者の脳機能活動を表すＭＲＩ計測データそのものであってもよいし、あるいは、各計測サイトにおいて被験者の脳機能活動を表すＭＲＩ計測データから相関行列の値を計算した後の相関行列の非対角要素である相関値そのもののデータであってもよい。

（新たな計測サイトが追加された場合のハーモナイゼーション算出処理）
図３８は、図３７において説明したような手続きにより、ハーモナイゼーション処理のための測定バイアスの算出が実行された後に、新たな計測サイトが追加された場合のハーモナイゼーション処理のための測定バイアスの算出処理を説明するための概念図である。

図３８を参照して、（Ｎｓ＋１）個目の計測サイトＭＳ．Ｎｓ＋１が新たに加わった場合、改めて、トラベリングサブジェクトＴＳ１を、また、この（Ｎｓ＋１）個の計測サイト全体について、巡回させ、上述したハーモナイゼーション算出処理と同様の手続きを行い、測定バイアスを改めて算出することが可能である。

以上説明したような処理を実行するシステムにより、複数の施設で計測される脳の活動の計測データに対して、各施設における測定バイアスを調整および補正することが可能である。これにより、複数の施設での計測データに基づく、脳機能結合相関値を調整することが可能となる。

また、本実施の形態のような処理を実行するシステムによれば、複数の施設で計測される脳の活動の計測データをハーモナイズして、健常または疾患の状態を客観的に判断するためのデータを提供することが可能な脳活動分類器のハーモナイズ方法、脳活動分類器のハーモナイズシステム、脳活動解析システムおよび脳活動解析方法を実現することができる。

あるいは、本実施の形態のような処理を実行するシステムによれば、複数の施設で計測される脳の活動の計測データに対して、神経・精神疾患について、脳機能画像法によるバイオマーカー装置、バイオマーカ装置のプログラムを実現することが可能となる。

［実施の形態１の変形例］
なお、以上の説明では、各計測サイトでは、以下の被験者のデータが取得されていることを前提として、ハーモナイゼーション算出処理について説明した。

ｉ）患者のデータ
ｉｉ）健常者のデータ
ｉｉｉ）トラベリングサブジェクトのデータ
ただし、上述した「測定バイアス」の評価のみを目的とするのであれば、トラベリングサブジェクトについて計測したデータに基づいて、ＧＬＭＭ（一般化線形混合モデル）により、「測定バイアス」を評価することとしてもよい。

すなわち、トラベリングサブジェクトが、Ｎts人であるとする。
参加者要因（ｐ）、測定バイアス（ｍ）は、トラベリングサブジェクトのデータセットからのすべての参加者の機能的なコネクティビティの相関値に、回帰モデルを適合させることにより評価される。
ここでも、ベクトルは小文字（例えば、ｍ）によって表示され、ベクトルはすべて列ベクトルであると仮定される。
ベクトルの要素は、ｍ_kのように添字によって表示される。

参加者の特性を表わすために、1-of-Kのバイナリーコード体系を使用し、サイトｋに属する測定バイアスｍに対するターゲットベクトル（たとえば、ｘm）は、１に等しい要素ｋ以外は、全てゼロに等しい。
上付きの添字Tは、マトリックスかベクトルの置換を表示し、ｘ^Tは行ベクトルを表わす。
以上のような構成によって、「測定バイアス」を算出して、ハーモナイゼーション処理を実行する構成とすることも可能である。

実施の形態１の手法と、実施の形態１の変形例の手法との相違点としては，実施の形態１の変形例の手法のようにトラベリングサブジェクトのみを使用する場合には、標本バイアスは考慮しなくてよくなるという利点がある反面、実施の形態１の手法のように患者のデータと健常者のデータとの双方を含んだ方が、推定に使用されるデータ数が多くなるという利点がある。

したがって、標本バイアスを考慮しないことによる測定バイアスの推定精度の向上と、データ数が多いことによる測定バイアスの推定精度の向上とは、互いにトレードオフの関係にあるといえる。

そこで、特に限定されないが、たとえば、現実の運用の中で、各計測サイトにおいて使用する参加者（患者および健常者）のデータ数については、実験的に最適化をすることも可能である。その場合、特に限定されないが、各計測サイトから抽出する参加者の人数を割り振っておき、各計測サイトでは、ランダムに当該人数分のデータを抽出することとしてもよい。

［実際の脳活動計測結果に対するデータおよびその解析］
以下では、実際に公開されている多疾患データベースのデータに対して、上述したようなハーモナイゼーション処理を実行した結果について説明する。

図３９Ａおよび図３９Ｂは、それぞれ、ハーモナイゼーション処理に使用する、多疾患データベースのデータとトラベリングサブジェクトのデータセットとを示す概念図である。

図３９Ａに示すように、脳科学に対する戦略的研究プログラム（ＳＲＰＢＳ：Strategic Research Program for Brain Sciences）における多疾患データベース(https://bicr-resource.atr.jp/decnefpro/)として、公開されているデータを患者群および健常群のデータセットを使用する。

多疾患データセットとして、複数の計測サイト（図３９Ａでは、サイト１、サイト２、サイト３を、典型的に示す）で計測されたデータは、サイト間差として、精神疾患の患者についてのサンプリングバイアスと、健常者についてのサンプリングバイアスと、測定バイアスとを含む。

一方、図３９Ｂに示すように、トラベリングサブジェクトについてのデータセットは、複数の計測サイト（図３９Ｂでは、サイト１、サイト２、サイト３を、典型的に示す）で計測されたデータは、サイト間差として、測定バイアスのみを含む。
データセットの同時分析によってサイト間差を測定バイアスおよび標本バイアスに分割することができ、安静時の機能的なコネクティビティ上のそれらの効果サイズを精神疾患のものと定量的に比較する。

測定バイアスについては、定量的に、異なる撮像パラメーター、ＭＲＩ装置メーカー、各ＭＲＩスキャナ中の受信コイルのおよび数について、効果の大きさを比較した。
サイト間差に関連したこれらの制限を克服するために、トラベリングサブジェクトのデータセットを使用して、測定バイアスだけを除くことができるようなハーモナイゼーション法を実行する。

新しく提案する方法および既存の方法によってハーモナイズされた複数サイトの安静時の機能結合ＭＲＩデータに基づいて、被験者の情報（属性）（例えば、年齢）の回帰モデルおよび精神疾患のバイオマーカーを構築した。

これらの予測モデルの性能が、ハーモナイゼーションの方法に基づいて、どのように変わるか調べた。

（使用したデータセット）
以下のような３つの安静時の機能ＭＲＩデータセットを使用した:
（１）ＳＲＰＢＳの複数疾患のデータセット、
（２）１９人のトラベリングサブジェクトデータセット、および
（３）独立した検証用データセット。

（ＳＲＰＢＳの多疾患のデータセット）
図４０（図４０Ａから図４０Ｄ）は、ＳＲＰＢＳの多疾患データセットの内容を示す図である。

９つのサイトで診察された、５つの異なる疾患患者および健常な対照群（ＨＣ）を含んだこのデータセットは、以下の合計８０５人の参加者が含まれていた:
９つのサイトからの４８２人の健常者、５つのサイトからの、１６１人の大うつ病（ＭＤＤ）患者、１つのサイトからの４９人の自閉症（ＡＳＤ）患者、１つのサイトからの６５人の強迫障害（ＯＣＤ）の患者、また３つのサイトからの４８人の統合失調症（ＳＣＺ）患者のデータが含まれる。

図４１（図４１Ａから図４１Ｄ）は、各計測サイトにおける撮像プロトコルを示す図である。
安静時の機能ＭＲＩデータは、３つのサイト以外のすべてで統一された撮像プロトコルを使用して得られた。（http://www.cns.atr.jp/rs-fmri-protocol-2/）
このデータセットのサイト間差は、測定および標本バイアスの両方を含んでいた。
バイアス評価については、統一プロトコルを使用して得られたデータのみを使用した。
なお、ＯＣＤ患者はこの統一プロトコルを使用して、スキャンされなかったので、ＯＣＤのために疾患要因を評価することは行っていない。

表中の略号は、以下のようである。
ATT：国際電気通信基礎技術研究所のシーメンスTimTrioスキャナ;
ATV：国際電気通信基礎技術研究所のシーメンスVerioスキャナ;
KUT：京都大学のシーメンスTimTrioスキャナ;
SWA：昭和大学;
HUH：広島大学病院;
HKH：ヒロシマKajikawa病院;
COI：COI(広島大学);
KPM：京都府立医科大学;
UTO：東京大学;

ASD：自閉症障害。
MDD：大うつ病性障害。
OCD：強迫障害。
SCZ：統合失調症。

SIE：シーメンスfMRI装置。
GE：GE fMRI装置。
PHI：フィリップスfMRI装置。

（トラベリングサブジェクトのデータセット）
ＳＲＰＢＳデータセット中の計測サイトを横断して測定バイアスを評価するために、トラベリングサブジェクトのデータセットを取得した。

健康な９人の参加者（すべての男性; 年齢層: ２４－３２歳; 平均年齢: ２７±２．６歳）が、１２のサイトの各々で撮像され、１２のサイトにはＳＲＰＢＳデータセットを撮像した９つのサイトを含んでおり、合計４１１セッションの操作が実施された。

ＳＲＰＢＳの複数疾患のデータセットと同じ撮像プロトコルを使用して、このデータセットを得ることを試みたものの、各サイトのパラメーター・セッティングあるいは慣習的なスキャン条件の制限のために、サイト間では、撮像プロトコルにいくつかの違いがあった。

たとえば、２つの位相エンコード方式の方向（P→AとA→P）、３社のＭＲＩメーカー（シーメンス、ＧＥおよびフィリップス）、４つのコイル・チャンネル数（8、12、24および32）および７つのスキャナ・タイプ（TimTrio、Verio、Skyra、スペクトル、MR750W、SignaHDxtおよびAchieva）というような相違が存在する。
同じ９人の参加者が１２のサイトにわたって撮像されたので、このデータセットのサイト間差は測定バイアスのみを含んでいた。

（独立した検証データセット）
安静時の機能結合ＭＲＩデータに基づいた、精神疾患の分類器および参加者の年齢を予測するモデルの汎化性能を検証するために、２つの疾患および７つのサイトをカバーする独立した検証コホートのデータを取得した。
このデータには、合計６２５人の参加者が含まれていた。
すなわち、６つのサイトからの４７６人の健常群（ＨＣｓ）、２つのサイトからのＭＤＤの９３人の患者、および１つのサイトからのＳＣＺの５６人の患者のデータが取得された。

（サイト間差と疾患効果の視覚化）
最初に、主成分分析（ＰＣＡ：principal component analysis）によって、ＳＲＰＢＳ複数疾患の安静時の機能結合ＭＲＩデータセットにおける、サイト間差および疾患効果を視覚化した。

図４２は、このような主成分分析によるサイト間差および疾患効果を視覚化を示すグラフである。
ここでの主成分分析は、教師なしの学習による次元低減手法に相当する。
被験者の機能的なコネクティビティは、上述したように、各参加者に対して２つの脳領域間の安静時の機能ＭＲＩの血中酸素濃度に依存する（ＢＯＬＤ）信号の時間的な相関性（ピアソンの相関係数を使用）として計算された。

機能的なコネクティビティは、全体の脳をカバーする２６８のノード（つまり脳領域）から成る機能的な脳地図に基づいて、定義される。
機能結合の相関を表すマトリックスの下三角行列の３５，７７８個（つまり、（２６８×２６７）／２）の結合強度（コネクティビティ）の値が使用される。

図４２に示すように、ＳＲＰＢＳ複数疾患のデータセット中の参加者のデータはすべて、最初の２つの主成分から成る２つの軸の上でプロットされた。
すなわち、ＳＲＰＢＳ複数疾患のデータセットのすべての参加者は、小さい薄い色のマーカーによって示されるように、最初の２つの主成分（ＰＣ）へ投影されている。
主成分１に対して、ＨＵＨサイトが明瞭に分離しており、それは、データの中の分散のほとんどを説明するほどであった。

ＡＳＤの患者はＳＷＡのサイトのみでスキャンされた。したがって、このサイトでスキャンされたＡＳＤ（▲）の患者および健常群ＨＣ（●）の平均は、ほとんど同一の位置へプロットされている。
図４２では、複数疾患のデータセットにおけるＰＣＡによる次元削減が示される。

（バイアス評価）
定量的に安静時の機能結合ＭＲＩデータのサイト間差を調べるために、測定バイアス、標本バイアスおよび診断の要因を特定した。
上述したように、各々の機能的なコネクティビティに対する相関値のすべてのサイトにわたる平均からの偏差として、各サイトの測定バイアスを定義した。健常者および精神疾患の患者の標本バイアスは、互いに異なると仮定する。

したがって、健常者群および各疾患を有する患者群に対して個別に、各サイトの標本バイアスを計算した。
疾患要因は、健常群の値からの偏差として定義された。
複数のサイトでサンプリングされた患者群であったという理由で、標本バイアスは、ＭＤＤとＳＣＺの患者に対して評価された。
対照的に、トラベリングサブジェクトのデータセットでは、参加者は固定されているので、測定バイアスのみを含んでいた。

トラベリングサブジェクトをＳＲＰＢＳ複数疾患のデータセットと組み合わせることによって、異なるサイトによって影響を受けた異なる要因として、同時に測定バイアスおよびサンプリングを評価した。機能的なコネクティビティに対する両方のタイプのバイアスと疾患の要因の効果を評価するために、以下のように、「線形混合効果モデル」を利用した。

（ＳＲＰＢＳ複数疾患のデータセットに対する線形混合効果モデル）
線形混合効果モデルでは、ＳＲＰＢＳ複数疾患のデータセットの各被験者のコネクティビティの相関値は、固定効果およびランダム効果からなる。
固定効果は、ベースラインとしてすべての参加者およびすべてのサイトにわたる平均相関値と、測定バイアス、標本バイアスおよび疾患要因の和を含む。
参加者の要因（つまり個人差）およびスキャン間の変化の結合した効果は、ランダム効果と見なされる。

（バイアスおよび要因の評価の詳細）
参加者要因（ｐ）、測定バイアス（ｍ）、標本バイアス（shc, smdd, sscz）および精神疾患要因（d）は、ＳＲＰＢＳ複数疾患のデータセットおよびトラベリングサブジェクトのデータセットからのすべての参加者の機能的なコネクティビティの相関値に回帰モデルを適合させることにより評価された。

この例では、ベクトルは小文字の太字（例えばm）によって表示され、ベクトルはすべて列ベクトルであると仮定される。
ベクトルの要素は、ｍ_kのように添字によって表示される。
上述したのと同様に、参加者の特性を表わすために、1-of-Kのバイナリーコード体系を使用し、サイトｋに属する測定バイアスｍに対するターゲットベクトル（たとえば、ｘm）は、１に等しい要素ｋ以外は、全てゼロに等しい。
したがって、参加者がどんなクラスにも属さなければ、ターゲットベクトルは、すべての要素が０に等しいベクトルである。
上付きの添字Tは、マトリックスかベクトルの置換を表示し、ｘ^Tは行ベクトルを表わす。

各コネクティビティについては、回帰モデルは以下のように書くことができる。

ここで、ｍは、測定バイアス（１２サイト×１）を表し、ｓhcは、健常群の標本バイアス（６サイト×１）を表し、ｓmddは、ＭＤＤ患者の標本バイアス（３サイト×１）を表し、ｓsczは、ＳＣＺ患者の標本バイアス（３サイト×１）を表し、ｄは、疾患要因（３疾患×１）、ｐは、参加者要因（９人のトラベリングサブジェクト×１）を表し、ｃｏｎｓｔは、すべてのサイトからの全ての参加者にわたる機能的コネクティビティの平均を表し、ｅ～Ｎ（０，γ^-1）は、ノイズを表す。

各機能的なコネクティビティの相関値に対しては、Ｌ２正規化による正規の通常の最小二乗回帰を使用して、それぞれのパラメーターを評価した。
正規化が適用されない場合、健常群に対する測定バイアスと標本バイアスの間のスプリアス反相関、および健常群に対する標本バイアスと精神疾患の患者に対する標本バイアスの間の擬似相関を観察した。これらの擬似相関の絶対平均を最小化するためにハイパーパラメータ・ラムダの調整をした。

（トラベリングサブジェクトのデータセットに対する線形混合効果モデル）
トラベリングサブジェクトのデータセットに対する線形混合効果モデルは、トラベリングサブジェクトのデータセットにおいて特定のスキャンに対する各参加者のコネクティビティの相関値は、固定効果およびランダム効果からなる。
固定効果は、すべての参加者およびすべてのサイトにわたる平均相関値と、参加者要因と、測定バイアスとの和を含む。
スキャン間の変化は、ランダム効果と見なされた。
各参加者に対して、すべての参加者にわたった平均からの脳機能コネクティビティの相関値の偏差として、参加者要因を定義した。
２つの異なるデータセットの機能的なコネクティビティの相関値に、同時に前述の２つの回帰モデルをフィッティングさせることにより、すべてのバイアスおよび要因を評価した。

以上を要約すると、バイアスあるいは各要因は、それぞれ、コネクティビティの相関値の数（つまり３５，７７８個）を反映する次元を含んだベクトルとして評価された。

（寄与サイズの分析）
定量的に要因間の大きさの関係を確認するために、線形混合効果モデルにおいて、寄与サイズを計算し比較して、各タイプのバイアスおよび要因が、どの程度、データの分散を説明するのかを決定した。

例えば、このモデル中で、測定バイアス（つまり最初の項）の寄与サイズは、以下の通り計算された。

ここで、Ｎmは，各要因の要素の数を表わし、Ｎは、コネクティビティの数を表し、Ｎsは、被験者の数を表し、Contribution size_mは、測定バイアスの寄与サイズの大きさを表す。

これらの式は、測定バイアス（例えば、位相エンコーディング方向、スキャナ、コイルおよびｆＭＲＩメーカー）と関係する個々の要因の寄与サイズを評価するために使用された。

特に、測定バイアスをこれらの要因へ分解し、その後、関連するパラメーターが評価された。
他のパラメーターは、それ以前に評価したものと同じ値に固定された。

寄与サイズの評価からは、以下のことが分かった。
１）機能的なコネクティビティに対する測定バイアスの効果サイズは、参加者要因より小さいことを示すものの、ほとんど疾患要因のものより大きく、測定バイアスが精神疾患に関する研究で重大な制限となることことを示唆する。
２）標本バイアスにおける最も大きな分散は、ＭＤＤ要因の分散よりも有意に大きく、標本バイアスにおける最も小さな分散は、疾患要因に対する分散の２分の１であった。
このことは、標本バイアスは、また、精神疾患に関する研究で重大な制限となることを示唆する。
３）参加者要因の標準偏差は、ＳＣＺ，ＭＤＤ，およびＡＳＤの疾患要因に対する標準偏差の約２倍であった。したがって、健常者の母集団内における個人の変動は、すべての機能結合を考慮した場合は、ＳＣＺ，ＭＤＤ，およびＡＳＤの患者間よりも大きかった。
４）さらに、測定バイアスの標準偏差は、疾患要因の標準偏差よりも、大抵は大きく、一方で、サンプリンバイアスの標準偏差は、疾患要因の標準偏差と同程度であった。

そのような関係性は、精神疾患あるいは発達障害に対する、安静時の機能結合ＭＲＩに基づく分類器を開発することを、大変、困難なものにしている。非常に少ない数の疾患に特定的な、あるいは、測定サイトに独立な異常な機能結合を、多数の結合から選択できる場合にのみ、複数のサイトについて、ロバストで、汎化可能な分類器を生成することを可能とする。

（測定バイアスに対する階層的クラスタリング分析）
各サイトｋに対して、測定バイアスｍk（Ｎ×１、ここで、Ｎは、機能的なコネクティビティの個数）間のピアソンの相関係数を計算し、測定バイアスにわたる相関係数に基づいた階層的クラスタリング分析を行った。

図４３は、階層的クラスタリング分析による樹状図である。
樹状図中の各リンケージの高さは、そのリンクによって連結されたクラスタ間で非類似性（１－ｒ）を表わす。

測定バイアスの特性を調べるにあたり、１２か所のサイトに対して、評価された測定バイアスベクトルの類似性が、位相エンコードの方向、ＭＲＩメーカー、コイル・タイプおよびスキャナ・タイプのようなＭＲＩスキャナの特定の特性を反映するものかどうかを調べた。
測定バイアスに対する同様のパターンのクラスタを発見するために、階層的クラスタ分析を使用した。

その結果、図４３に示すように、１２のサイトの測定バイアスは、最初のレベルで、位相エンコードの方向のクラスタに分割された。
測定バイアスは、第２のレベルで、ｆＭＲＩメーカーのクラスタに分割され、さらに、コイル・タイプ・クラスタへさらに分割され、そして、スキャナ・モデル・クラスタに分割された。

図４４は、各要因についての寄与サイズを示す図である。
図４４に示すように、各要因の寄与を評価するために、同じモデルを使用することにより、定量的に要因の大きさの関係を確認した。

寄与サイズは、位相エンコーディング方向(0.0391)に対して最大であり、続いて、ｆＭＲＩメーカー(0.0318)、コイル・タイプ(0.0239)およびスキャナ・モデル(0.0152)が、これに続くという関係であった。
これらにより、測定バイアスに影響を及ぼす主な要因が、位相エンコーディング方向、続いて、ｆＭＲＩメーカー、コイル・タイプおよびスキャナ・モデルの相違であることを示す。

（ハーモナイゼーション効果の視覚化）
次に、トラベリングサブジェクトのデータセットを使用して、測定バイアスだけを除くことを可能とするハーモナイゼーション法について説明する。
上述のとおり、線形混合効果モデルを使用して、測定バイアスを標本バイアスとは別に評価した。
この方法により、ＳＲＰＢＳの複数疾患のデータセットから測定バイアスだけを取り除き、生物学的な情報を含んでいる標本バイアスを維持することができる。

（トラベリングサブジェクトのハーモナイゼーション）
測定バイアスを評価するために、ＨＫＨサイトでの位相エンコーディング方向が、ＳＲＰＢＳ複数疾患のデータセットとトラベリングサブジェクトのデータセットの間で異なっていたので、以下では、測定バイアスとは別個に、式中に位相エンコード要因（ｐａ_q）を含めることとした。

測定バイアスおよび位相エンコード要因は、ＳＲＰＢＳ複数疾患のデータセットおよびトラベリングサブジェクトのデータセットを結合したものに、回帰モデルを適合させることにより評価された。

回帰モデルは、以下のようになる。

ここで、ｐａは、位相エンコード要因（２つの位相エンコード方向×１）を表わす。
各機能的なコネクティビティの相関値に正規のＬＳ正規化による最小二乗法回帰で正規化を実行して使用して、パラメーターをそれぞれ評価した。
サイト間差は、評価されたサイト間差および位相エンコード要因を引くことにより除去された。

したがって、標本バイアスがハーモナイズされた機能的なコネクティビティの相関値は以下のようになる。

ここで、ｍは、評価された測定バイアスを表し、ｐａ（ハット）は、評価された位相エンコード要因を表わす。

図４５は、ハーモナイゼーションプロセスの影響を視覚化した図であり、図４２と対比される図である。
図４５では、ＳＲＰＢＳの複数疾患のデータセットから測定バイアスだけを除いた後にデータをプロットした。

図４２ではハーモナイゼーションの前のデータを反映しており、図４２に示されたデータと比較して、図４５においては、ＨＵＨサイトは、原点側（つまり母集団の平均）により近くなるように、大きく移動しており、また他のサイトから大きく離れるということはなくなった。
図４５に示した結果は、図４２で見られたようなＨＵＨサイトの分離は、測定バイアスによりもたらされており、それは、ハーモナイゼーションによって取り除かれ得ることを示している。

さらに、ハーモナイゼーションは、同じサイトでスキャンされた患者と健常者群とを識別するのにも有効であった。
ＡＳＤの患者は、ＳＷＡのサイトのみでスキャンされたので、このサイトでスキャンされたＡＳＤ患者（▲）および健常者（○）の平均は、図４２では、ほとんど同一の位置へプロットされている。
しかしながら、２つのシンボルは、図４５では、お互いに明白に分離されている。

[実施の形態２]
実施の形態１では、脳活動計測装置（ｆＭＲＩ装置）で、複数の測定場所で計測された脳活動データを計測し、この脳活動データに基づいて、バイオマーカーの生成およびバイオマーカーによる診断ラベルの推定（予測）を行う構成として、分散処理により行う一例の構成を説明した。

ただし、ｉ）バイオマーカーを機械学習により訓練するための脳活動データの計測（データ収集）、ｉｉ）バイオマーカーの機械学習による生成処理および特定の被験者についてのバイオマーカーによる診断ラベルの推定（予測）の処理（推定処理）、ｉｉｉ）上記特定の被験者についての脳活動データの計測（対象者の脳活動計測）を、それぞれ、異なる施設で分散して実行する構成とすることも可能である。

図４６は、データ収集、推定処理および対象者の脳活動計測を、分散して処理する場合の一例を示す機能ブロック図である。
図４６を参照して、サイト１００．１～１００．Ｎは、患者群、健常者群のデータを脳活動計測装置により計測する施設であり、データセンター２００は、サイト１００．１～１００．Ｎsからの計測データを管理する。

計算処理システム３００は、データセンター２００に格納されたデータから識別器を生成する。
また、計算処理システム３００のハーモナイゼーション算出部３０２０は、サイト１００．１～１００．ＮsおよびＭＲＩ装置４１０のサイトを含めて、ハーモナイゼーション処理を実行するものとする。

ＭＲＩ装置４１０は、計算処理システム３００上の識別器の結果を利用する別サイトに設けられており、特定の被験者について脳活動のデータを計測する。
コンピュータ４００は、ＭＲＩ装置４１０が設けられる別サイトに設置され、ＭＲＩ装置４１０の測定データから特定被験者の脳の機能結合の相関データを算出し、機能結合の相関データを計算処理システム３００に送信して、返信されてくる識別器の結果を利用する。

データセンター２００は、サイト１００．１～１００．Ｎsから送信されてくる患者群および健常者群のＭＲＩ測定データ３１０２と、ＭＲＩ測定データ３１０２に関連付けられた被験者の人属性情報３１０４とを格納し、計算処理システム３００からのアクセスに従って、これらのデータを計算処理システム３００に送信する。
計算処理システム３００は、通信インタフェース２０９０を介して、データセンター２００からのＭＲＩ測定データ３１０２および被験者の人属性情報３１０４を受信する。

なお、データセンター２００、計算処理システム３００、コンピュータ４００のハードウェアの構成は、基本的に、図５で説明した「データ処理部３２」の構成と同様であるので、その説明は繰り返さない。

図４６に戻って、相関行列算出部３００２、相関値補正処理部３００４、疾患識別器生成部３００８、クラスタリング分類器生成部３０１０および判別値算出部３０１２、ならびに、機能結合の相関行列のデータ３１０６、測定バイアスデータ３１０８、補正後相関値データ３１１０および識別器データ３１１２については、実施の形態１で説明したのと同様であるので、その説明は、繰り返さない。

ＭＲＩ装置４１０は、診断ラベルの推定対象となる被検者の脳活動データを計測し、コンピュータ４００の処理装置４０４０は、計測されたＭＲＩ測定データ４１０２を不揮発性記憶装置４１００に格納する。

さらに、コンピュータ４００の処理装置４０４０は、ＭＲＩ測定データ４１０２に基づいて、相関行列算出部３００２と同様にして、機能結合の相関行列のデータ４１０６を算出し、不揮発性記憶装置４１００に格納する。

コンピュータ４００のユーザから診断の対象となる疾患が指定され、当該ユーザの送信の指示に従い、コンピュータ４００が、機能結合の相関行列のデータ４１０６を計算処理システム３００に送信する。これに応じて、計算処理システム３００は、ＭＲＩ装置４１０の設置されるサイトに対応したハーモナイゼーション処理を実行して、判別値算出部３０１２は、指定された診断ラベルについての判別結果やサブタイプについての評価結果を算出し、計算処理システム３００は、通信インタフェース２０９０を介して、コンピュータ４００に送信する。

コンピュータ４００では、図示しない表示装置などを介して、ユーザに対して、判別結果を知らせる。

このような構成とすることで、より多くの被験者について収集したデータに基づいて、識別器による診断ラベルの推定結果を提供することが可能となる。

また、データセンター２００と計算処理システム３００とを別個の管理者が管理する形態とすることも可能で、その場合、データセンター２００にアクセスできるコンピュータを制限することで、データセンター２００に格納される被験者の情報のセキュリティを向上させることも可能となる。

さらに、計算処理システム３００の運営主体からみると、「識別器による判別のサービスを受ける側（コンピュータ４００）」に対して、識別器についての情報や「測定バイアス」に関する情報は、一切提供しなくても、「判別結果を提供するサービス」を行うことが可能となる。

なお、以上の実施の形態１および実施の形態２の説明では、脳機能画像法により脳活動を時系列に計測するための脳活動検出装置としては、リアルタイムｆＭＲＩを用いるものとして説明した。ただし、脳活動検出装置として、上述したｆＭＲＩ、脳磁計、近赤外光計測装置（ＮＩＲＳ）、脳波計、またはこれらの組み合わせを使用することができる。たとえば、これらの組み合わせを用いる場合、ｆＭＲＩとＮＩＲＳとは、脳内の血流変化に関連する信号を検出するものであり、高空間分解能である。一方で、脳磁計や脳波計は、脳活動に伴う電磁場の変化を検出するための高時間分解能であるという特徴をもつ。したがって、たとえば、ｆＭＲＩと脳磁計とを組み合わせれば、空間的にも時間的にも高分解能で脳活動を計測することができる。あるいは、ＮＩＲＳと脳波計とを組み合わせても、同様に空間的にも時間的にも高分解能で脳活動を計測するシステムを小型で携帯可能な大きさで構成することも可能である。

以上のような構成により、神経・精神疾患に対して、脳機能画像法によるバイオマーカーとして機能する脳活動解析装置および脳活動解析方法を実現することが可能となる。

また、以上の説明では、被験者の属性として「診断ラベル」を含んだ場合について、機械学習による識別器の生成により、当該識別器をバイオマーカーとして機能させる例について説明したが、本発明は、必ずしもこのような場合に限定されず、事前に機械学習の対象となる測定結果を得る対象の被験者群が、客観的な方法により、複数のクラスに分けられており、被験者の脳領域間（関心領域間）の活動度の相関（結合）を測定し、測定結果に対する機械学習により、クラスに対する識別器が生成できるものであれば、他の判別のために使用されるものであってもよい。

また、上述のとおり、このような判別は、ある属性に属する可能性を確率として表示するものであってもよい。

したがって、たとえば、ある「訓練」や「行動パターン」をとることが、被験者にとっての健康増進に役立つかを客観的に評価することができる。また、実際には、疾患に至っていない状態（「未病」）であっても、ある「食物」「飲料」等のような摂取物や、ある活動などが、より健康状態に近づくために効果があるのか、ということを客観的に評価することも可能である。

また、未病の状態においても、上述したように、たとえば、「健常である確率は、○○％」といような表示が出力されれば、ユーザに対して、健康状態について客観的な数値としての表示を行うことができる。このとき、出力されるのは、必ずしも確率でなくとも、「健康度合いの連続値、たとえば、健常である確率」をスコアに変換したものを表示することとしてもよい。このような表示を行うことで、本実施の形態の装置を、診断の支援以外にも、ユーザの健康管理のための装置として使用することも可能である。

［実施の形態３］
以上の説明では、すべての計測サイトについて、平等に、トラベリングサブジェクトが移動して計測を行うことにより、測定バイアスを評価するとの構成であった。

図４７は、実施の形態３のトラベリングサブジェクトの巡回方式を示す図である。
図４７に示すように、拠点となる「拠点計測サイトＭＳ．１～ＭＳ．Ｎｓ」について、トラベリングサブジェクトＴＳ１が巡回していることを前提とする。

これに対して、たとえば、これらの拠点計測サイトＭＳ．１～ＭＳ．Ｎｓのうちの拠点計測サイトＭＳ．２については、予め配下となるＭＳ２．１～ＭＳ２．ｎが存在するとして、新たに「計測サイトＭＳ．２．ｎ＋１」が追加された場合は、この配下の範囲でトラベリングサブジェクトＴＳ２が巡回する。他の拠点計測サイトについても同様である。

すなわち、拠点計測サイトＭＳ．２の測定バイアスは、トラベリングサブジェクトＴＳ１で評価された値に固定して、計測サイトＭＳ．２、ＭＳ２．１、ＭＳ２．ｎ，ＭＳ２．ｎ＋１について、トラベリングサブジェクトＴＳ２が巡回した計測結果に基づいて、測定バイアスを決定するという構成とすることも可能である。

たとえば、「拠点計測サイトＭＳ．１～ＭＳ．Ｎｓ」は、予め定められた地域ごとに存在すとして、日本国内の北海道、東北、関東、…、関西、…九州というような各地域に、１か所ずつ設置されており、配下サイトというのは、たとえば、その地域の１つの関西に位置する計測サイト、というような構成とすることも考えられる。

あるいは、「拠点計測サイトＭＳ．１～ＭＳ．Ｎｓ」は、予め定められたＭＲＩ装置の機種ごとに定められていてもよい。この場合、配下の計測サイトは、拠点サイトと同一のＭＲＩ装置の機種が設置された計測サイトということになる。

あるいは、「拠点計測サイトＭＳ．１～ＭＳ．Ｎｓ」は、予め定められた地域ごとに、定められたＭＲＩ装置の機種ごとに定められていてもよい。この場合、配下の計測サイトは、拠点サイトと同一の地域内で、同一のＭＲＩ装置の機種が設置された計測サイトということになる。

このような構成であっても、実施の形態１と同様の効果を奏することが可能である。

今回開示された実施の形態は、本発明を具体的に実施するための構成の例示であって、本発明の技術的範囲を制限するものではない。本発明の技術的範囲は、実施の形態の説明ではなく、請求の範囲によって示されるものであり、請求の範囲の文言上の範囲および均等の意味の範囲内での変更が含まれることが意図される。

２被験者、６ディスプレイ、１０ＭＲＩ装置、１１磁場印加機構、１２静磁場発生コイル、１４傾斜磁場発生コイル、１６ＲＦ照射部、１８寝台、２０受信コイル、２１駆動部、２２静磁場電源、２４傾斜磁場電源、２６信号送信部、２８信号受信部、３０寝台駆動部、３２データ処理部、３６記憶部、３８表示部、４０入力部、４２制御部、４４インタフェース部、４６データ収集部、４８画像処理部、５０ネットワークインタフェース。

Claims

対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリングを実行するための脳機能結合相関値のクラスタリング装置であって、
前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者について、脳活動の計測値に基づいて、前記クラスタリングの処理を実行するための計算処理システムを備え、前記計算処理システムは、記憶装置と演算装置とを含み、前記演算装置は、
ｉ）各前記複数の被験者について、所定の複数の脳領野ペア間の脳活動の時間相関をそれぞれ表す複数の脳機能結合相関値に基づく特徴量を前記記憶装置に格納し、
ｉｉ）前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するように構成され、
前記演算装置は、
前記識別器モデルを生成する機械学習において、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成し、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択し、
前記演算装置は、さらに、
選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングしてクラスタの分類器を生成する、脳機能結合相関値のクラスタリング装置。
前記脳機能結合相関値のクラスタリング装置は、
複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各前記複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受け、
前記計算処理システムは、
各前記複数の被験者について前記複数の脳機能結合相関値について、前記計測サイトの測定バイアスを除去するように補正することで、補正された調整値を前記特徴量として前記記憶装置に格納するハーモナイゼーション算出手段を含む、請求項１記載の脳機能結合相関値のクラスタリング装置。
前記機械学習による識別器の生成の処理は、前記複数の学習用サブサンプルに対してそれぞれ複数の識別器サブモデルを生成し、前記複数の識別器サブモデルを統合して前記識別器モデルを生成するアンサンブル学習である、請求項１または請求項２記載の脳機能結合相関値のクラスタリング装置。
前記属性は、所定の精神疾患であるとの診断結果のラベルで表現され、
前記クラスタリングは、データ駆動による機械学習により、前記第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理である、請求項１～請求項３のいずれか１項に記載の脳機能結合相関値のクラスタリング装置。
前記演算装置は、前記機械学習による識別器の生成において、
ｉ）前記調整値を、機械学習用のトレーニングデータセットと検証用のテストデータセットとに分割し、
ｉｉ）前記トレーニングデータセットに対して、アンダーサンプリングおよびサブサンプリングを所定数だけ実行して、前記所定数の学習用サブサンプルを生成し、
ｉｉｉ）前記学習用サブサンプルごとに識別器サブモデルを生成し、
ｉｖ）前記識別器サブモデルの出力を統合して、前記属性の有無に対する識別器モデルを生成する、請求項２記載の脳機能結合相関値のクラスタリング装置。
前記機械学習による識別器の生成の処理は、外側の交差検証と内側の交差検証を有する入れ子構造の交差検証であり、
前記演算装置は、前記入れ子構造の交差検証の処理において、
ｉ）前記外側の交差検証をＫ分割交差検証として、前記調整値を、機械学習用のトレーニングデータセットと検証用のテストデータセットとに分割し、
ｉｉ）前記トレーニングデータセットに対して、アンダーサンプリングおよびサブサンプリングを所定数だけ実行して、前記所定数の学習用サブサンプルを生成し、
ｉｉｉ）前記Ｋ分割交差検証の各ループにおいて、前記内側の交差検証で、ハイパーパラメータを調整して、前記学習用サブサンプルごとに識別器サブモデルを生成し、
ｉｖ）前記識別器サブモデルに基づいて、前記属性の有無に対する識別器モデルを生成する、請求項２記載の脳機能結合相関値のクラスタリング装置。
前記機械学習による識別器の生成の処理は、特徴量選択を伴う機械学習法であり、
前記クラスタリングのための特徴量の選択においては、前記和集合に属する特徴量の重要度を、前記識別器サブモデルを生成する際に前記特徴量が選択される頻度のランキングにより決定する、請求項３に記載の脳機能結合相関値のクラスタリング装置。
前記機械学習による識別器の生成の処理は、ランダムフォレスト法であり、
前記クラスタリングのための特徴量の選択においては、前記和集合に属する特徴量の重要度は、各特徴量について前記ランダムフォレスト法においてジニ不純度に基づいて算出される重要度である、請求項３に記載の脳機能結合相関値のクラスタリング装置。
前記機械学習による識別器の生成の処理は、Ｌ２正則化による機械学習法であり、
前記クラスタリングのための特徴量の選択においては、前記和集合に属する特徴量の重要度を、Ｌ２正則化により算出される前記識別器サブモデルにおける特徴量の重みに基づくランキングにより決定する、請求項３に記載の脳機能結合相関値のクラスタリング装置。
前記記憶装置は、前記複数の計測サイトにおいて共通に計測対象となる複数の旅行被験者に対して、各前記旅行被験者の予め定められた複数の脳領域について脳活動を計測した結果を予め格納しており、
前記演算装置は、
前記複数の脳領域のペアについての脳活動の時間相関を表す脳機能結合行列の所定の要素を各前記旅行被験者について算出し、
一般化線形混合モデル法を用いることで、前記脳機能結合行列の所定の要素ごとに、前記複数の計測サイトおよび前記複数の旅行被験者にわたる当該要素の平均に対する各計測サイトにおける固定効果として、前記測定バイアスを算出する、請求項２記載の脳機能結合相関値のクラスタリング装置。
前記属性は、所定の精神疾患であるとの診断結果のラベルで表現され、
前記クラスタリングは、データ駆動による機械学習により、前記第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理であり、
前記演算装置は、前記複数の計測サイト以外の計測サイトで計測された対象者についての計測データに基づいて、前記サブタイプへの分類処理を実行する、請求項２記載の脳機能結合相関値のクラスタリング装置。
対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリングを実行するための脳機能結合相関値のクラスタリングシステムであって、
前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者についての脳活動を時系列で計測するために、複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置と、
前記複数の被験者について、脳活動の計測値に基づいて、前記クラスタリングの処理を実行するための計算処理システムとを備え、前記計算処理システムは、記憶装置と演算装置とを含み、前記演算装置は、
ｉ）各前記複数の被験者について、所定の複数の脳領野ペア間の脳活動の時間相関をそれぞれ表す複数の脳機能結合相関値に基づく特徴量を前記記憶装置に格納し、
ｉｉ）前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するように構成され、
前記演算装置は、
前記識別器モデルを生成する機械学習において、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成し、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択し、
前記演算装置は、さらに、
選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングしてクラスタの分類器を生成する、脳機能結合相関値のクラスタリングシステム。
前記計算処理システムは、
複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各前記複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受け、
各前記複数の被験者について前記複数の脳機能結合相関値について、前記計測サイトの測定バイアスを除去するように補正することで、補正された調整値を前記特徴量として前記記憶装置に格納するハーモナイゼーション算出手段を含む、請求項１２記載の脳機能結合相関値のクラスタリングシステム。
前記属性は、所定の精神疾患であるとの診断結果のラベルで表現され、
前記クラスタリングは、データ駆動による機械学習により、前記第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理である、請求項１２または請求項１３記載の脳機能結合相関値のクラスタリングシステム。
対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリング処理を計算処理システムが実行するための脳機能結合相関値のクラスタリング方法であって、
前記計算処理システムは、記憶装置と演算装置とを含み、
前記演算装置が、前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を前記記憶装置に格納するステップと、
前記演算装置が、前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを備え、
前記識別器モデルを生成する機械学習を実行するステップは、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、
さらに、前記演算装置が、選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを備える、脳機能結合相関値のクラスタリング方法。
対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳機能結合相関値の分類器プログラムであって、
前記分類器プログラムは、各前記クラスタの確率分布のモデルに基づき、コンピュータが、前記入力データが最大事後確率を有するクラスタへ分類する分類機能を有し、
前記計算処理システムは、記憶装置と演算装置とを含み、
前記計算処理システムは、前記クラスタリング処理に基づく前記分類器プログラムの生成処理において、
前記演算装置が、前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を前記記憶装置に格納するステップと、
前記演算装置が、前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、
前記識別器モデルを生成する機械学習を実行するステップは、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、
さらに、前記演算装置が、選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを実行する、脳機能結合相関値の分類器プログラム。
前記計算処理システムは、
複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各前記複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受信するステップと、
各前記複数の被験者についての複数の前記脳機能結合相関値について、前記計測サイトの測定バイアスを除去するように補正することで、補正された調整値を前記特徴量として前記記憶装置に格納するハーモナイゼーションを実行するステップとを実行する、請求項１６記載の脳機能結合相関値の分類器プログラム。
前記属性は、所定の精神疾患であるとの診断結果のラベルで表現され、
前記クラスタリングは、データ駆動による機械学習により、前記第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理である、請求項１６または請求項１７記載の脳機能結合相関値の分類器プログラム。
対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳活動マーカー分類システムであって、
前記脳活動マーカー分類システムは、各前記クラスタの確率分布のモデルに基づき、コンピュータが、前記入力データが最大事後確率を有するクラスタへ分類する分類機能を有し、
前記計算処理システムは、記憶装置と演算装置とを含み、
前記計算処理システムは、前記クラスタリング処理に基づく前記脳活動マーカー分類システムの生成処理において、
前記演算装置が、前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を前記記憶装置に格納するステップと、
前記演算装置が、前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、
前記識別器モデルを生成する機械学習を実行するステップは、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、
さらに、前記演算装置が、選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングしてクラスタの分類器を生成するステップを実行する、脳活動マーカー分類システム。
前記計算処理システムは、
複数の計測サイトにそれぞれ設けられた複数の脳活動計測装置から、各前記複数の被験者の所定の複数の脳領野ペア間の脳活動の時間相関を表現する情報を受信するステップと、
各前記複数の被験者について前記脳活動の時間相関を表す複数の脳機能結合相関値について、前記計測サイトの測定バイアスを除去するように補正することで、補正された調整値を前記特徴量として前記記憶装置に格納するハーモナイゼーションを実行するステップとを実行する、請求項１９記載の脳活動マーカー分類システム。
前記属性は、所定の精神疾患であるとの診断結果のラベルで表現され、
前記クラスタリングは、データ駆動による機械学習により、前記第１の被験者群を少なくとも１つのサブタイプのクラスタに分類する処理である、請求項１９または請求項２０記載の脳活動マーカー分類システム。
対象者の脳活動の計測結果に基づいて、前記対象者において、少なくとも１つの所定の属性を有する前記対象者に対するクラスタリング処理を計算処理システムが実行することにより生成され、入力データに対してクラスタリング処理の結果に対応する各クラスタへの分類をコンピュータが行うための脳機能結合相関値のクラスタリング分類器モデルであって、
前記クラスタリング分類器モデルは、学習データに含まれる前記対象者を特徴づける特徴量群を分割したビューごとに、各前記ビューに含まれる前記特徴量の情報と、各前記ビューの前記対象者のクラスタごとの確率密度関数を特定する情報とに基づいて、前記入力データに対して算出される前記確率密度関数の値に応じて、前記入力データが最大事後確率を有するクラスタへ分類する機能を有し、
前記計算処理システムは、記憶装置と演算装置とを含み、
前記計算処理システムは、前記クラスタリング処理に基づく前記クラスタリング分類器モデルの生成処理において、
前記演算装置が、前記所定の属性を有する第１の被験者群と前記所定の属性を有さない第２の被験者群とを含む複数の被験者の各々について、所定の複数の脳領野ペア間の脳活動の時間相関を表す脳機能結合相関値に基づく特徴量を前記記憶装置に格納するステップと、
前記演算装置が、前記記憶装置に格納された前記特徴量に基づいて、前記属性の有無を判別するための識別器モデルを生成する機械学習を教師あり学習で実行するステップとを実行し、
前記識別器モデルを生成する機械学習を実行するステップは、
前記第１の被験者群と前記第２の被験者群から、アンダーサンプリングおよびサブサンプリングを実行して、複数の学習用サブサンプルを生成するステップと、
前記学習用サブサンプルのそれぞれについて、機械学習による識別器の生成において使用される特徴量の和集合から、前記和集合に属する特徴量の重要度に応じて、クラスタリングのための特徴量を選択するステップとを含み、
さらに、前記演算装置が、選択された前記クラスタリングのための特徴量に基づいて、教師なし学習の多重共クラスタリング法により、前記第１の被験者群をクラスタリングして、前記特徴量を前記ビューへ分割し、各前記ビューにおける前記対象者のクラスタごとの確率密度関数を生成するステップを実行する、脳機能結合相関値のクラスタリング分類器モデル。