遺伝子検査とライフスタイルの行動変容に関するSystematic Review

遺伝子検査サービスを行動変容理論と絡めた現状唯一？のSystematic Reviewだったのでサマリながら理解します。
（query："behaviour change" AND "personalization"）

論文を読む目的

  • 行動変容理論を組み合わせた介入に関して雰囲気掴みたい
  • 遺伝子検査を介した行動変容実現にはどのような課題が生じているか把握する

何をした論文か

  • 遺伝子検査の介入による行動変容 (栄養、身体活動、睡眠、および喫煙) に関する研究論文の詳細で包括的な評価を実施
  • 遺伝子検査を介した行動変容の介入に関する既存の研究では、理論的介入（TPB）の検討はほとんど行われておらず、介入の質としても低いことがほとんどであることを明らかにした
  • 質の高い（RCT）研究では、遺伝子検査を介して行われた実用的な推奨事項の提供は遺伝情報のみの提供よりも、行動の変化を促進する可能性を高めていることを明らかにした
  • 遺伝子検査を介した行動変容の介入においては、高品質の遺伝的介入が参加者に提供されることを保証する必要があり、研究デザインと分析においては計画行動理論（TPB）などの検証済みの理論を考慮するべきであることを提案した

先行研究と比較した優位点

  • 遺伝子検査を用いた行動変容研究では、研究結果に大きな影響を与える可能性のある交絡因子を制御するために、検証済みの行動変容理論を考慮することが重要だが、既存の研究では検証済みの行動変化理論と遺伝的介入の質の評価を使用した説明は為されていない
  • バイアスのリスクのみを評価する従来の体系的なレビュープロセスを使用するのではなく、バイアスのリスクを超えた要因（遺伝的介入の質、および理論的介入の検討（主に計画行動変容理論））が遺伝子検査と行動の変化に関連する研究結果に影響を与えることを実証した

追加の議論

  • 遺伝子検査を通した介入には明らかに偏りがあり、多くの論文が質の「悪い」介入が行われていると評価された。大部分の研究において遺伝子検査を通した介入効果が得られなかったと報告されたのは、研究において参加者に高品質の介入を提供していなかったことが要因である可能性が高い
  • 栄養、身体活動、喫煙習慣は複数の遺伝子介入研究で研究されてきたが、睡眠は遺伝学と行動変化の研究が不十分な領域のままである（本Systematic Reviewでhit数0件）。睡眠に関しては、血糖コントロール、過体重、肥満等との関連があることが既知の研究で明らかとなっている。健康転帰に影響を与える可能性のある遺伝子と睡眠の相互作用を決定していくようなアプローチが有効であると考えられる

その他メモ

  • 健康リスクや生活習慣に関する個別化された情報と推奨事項を比較的低コストで提供するため、遺伝子検査は臨床現場を中心に使用されてきている。一方で、遺伝子検査が生活習慣の変化を促進するかどうかを評価する研究では、相反する結果が得られている
  • 計画行動理論 (TPB) は、おそらく学界で最も広く受け入れられている行動変容理論。本理論では、行動への態度（attitudes）、主観的規範（subjective norms）、および知覚された行動制御（perceived behavioural control） が行動を予測するために使用できる重要な構成要素であると仮定している 引用：Theory of planned behavior - Wikipedia

読むべき参考文献

パーソナライズされた行動変容に関する研究を前進させるための重要な次のステップとして、TPB の観点からシステマティック レビューを完了することを推奨 Horne J, Madill J, Gilliland J: Incorporating the “Theory of Planned Behaviour” into personalized healthcare behaviour change research: a call to action. Pers Med 2017; 14: 521–529. （https://dx.doi.org/10.2217/pme-2017-0038）

計画的行動理論（TPB）に関して Ajzen I: The theory of planned behaviour: reactions and reflections. Psychol Health 2011; 26: 1113–1127. （https://pubmed.ncbi.nlm.nih.gov/21929476/）

分数多項式を使用した分析：重症患者における初期乳酸および正常化時間とハザードとの関連

生存時間分析 論文要約 統計

実務において、非線形な関係性を説明性の担保できる統計的手法でモデリングしたい状況となっております。David W. Hosmerの本において、分数多項式解析が使用されていたのでその実用例を調査します。

論文を読む目的

  • 分数多項式を使用したモデル構築フローを理解する
  • 分数多項式を使用することの優位点を理解する

何をした論文か

  • 十分に確立されていなかった乳酸と死亡率の関係性を明らかにした
  • 乳酸の正常化は、乳酸の正常化がない場合と比較して、死亡率を有意に低下させることを示した
  • 乳酸の正常化の遅延は、ICU入室後150時間以内に測定された場合、死亡のリスクが高くなることを示した

先行研究と比較した優位点

  • 既存研究のほとんどは観察研究であり、交絡因子の影響を受けているが、この影響を取り除くために多変数回帰分析が採用されていた。また、共変量の線形仮定がテストされていない。線形の仮定に違反した場合、モデルから得られたオッズ比は、乳酸値の全範囲には当てはまらない。本研究では分数多項式を使用しており、より柔軟なモデリングが可能となる。
  • 治療後の迅速な乳酸正常化は良好な転帰と関連していることが示唆されている。既存研究では、乳酸正常化の時点が事前定義されており、臨床医が乳酸をいつ再チェックする必要があるかを判断するのが困難であった。モデルに乳酸正常化の時間を含めることで正規化時間と死亡の危険性との関連を調査した。

追加の議論

分数多項式を使用することによって明らかになったこと

分数多項式を使用することで、ある定義域における説明変数とアウトカムとの関係性が明らかとなる。
本論文においては、

  • 初期乳酸レベルの増加（L0）はハザードに寄与し、その傾きは4〜10 mmol/Lの間で急勾配であり、その後、急勾配は徐々に減衰すること。これは組織灌流を改善し、乳酸クリアランスを伴う蘇生バンドルが、乳酸レベルが4〜10 mmol/Lの患者にとって最も有益である可能性があることを示唆。
  • 乳酸正常化の患者では、150時間以内において正常化の時間（T）がかかるほど、臨床転帰は悪化すること。すなわち、ICU入室後150時間以内においては、乳酸の正常化を迅速に行う必要があることを示唆。

非線形性を許容することにより生じる非単調性は説明できるものなのか

Tについて単調性が保証されておらず、250時間辺りを境目に低下傾向が見えるが、これは説明できるのか。
→ 高い正常化時間を持ったデータ数が少ないため、信頼区間幅が150h以降広くなっている。150h以降はデータから述べることはできないとみなし、説明は不要
→ 負の値をとりうることになるからといって非線形モデルをあきらめる必要はない

その他メモ

多変量解析のステップ

統計的に有意でも交絡因子でもない変数のみを除外した後、非線形を考慮できる分数多項式モデルを作成している。

  1. 単変量解析で P-value <0.2 のすべての変数を使用して初期の多変数モデルを構築
  2. Wald検定からP>0.1の場合、初期モデルの変数は削除
  3. 除去された共変量が乳酸係数に有意な変化（> 20％の変化）をもたらした場合、それは交絡因子であると考えられ、モデルに戻す
  4. 共変量を削除できなくなるまで繰り返し、予備的主効果モデルを構築
  5. –2、–1、–0.5、0、0.5、1、2、3のセットから構築された二次の分数多項式回帰モデルと一次の回帰モデルの逸脱度を検定し、P<0.05の場合採用（$\chi2$検定）

使用しているデータセット

データセットには、MIMIC-IIを使用。数万人の集中治療室（ICU）患者について、患者モニターから取得した生理学的信号とバイタルサインの時系列、および病院の医療情報システムから取得した包括的な臨床データが含まれている。

読むべき参考文献

分数多項式をSAS, STATA, Rで使用できる形で公開
Sauerbrei, W., C. Meier-Hirmer, A. Benner, and P. Royston. 2006. “Multivariable Regression Model Building by Using Fractional Polynomials: Description of SAS, STATA and R Programs.” Computational Statistics & Data Analysis 50 (12): 3464–85. (https://www.sciencedirect.com/science/article/abs/pii/S0167947305001623)

本論文の著者が分数多項式を解説
Zhang, Zhongheng. 2016. “Multivariable Fractional Polynomial Method for Regression Model.” Annals of Translational Medicine 4 (9): 174. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4876277/)

分析・開発環境構築（Windows編）

Windows Python

かなり時間が空いてしまった。ブログを継続するのはなんと難しいことか。 最近購入したデスクトップPCに環境構築する機会があったので、kimotonの分析・開発環境をここにまとめておきます。

WSL2のインストール

今やWindows使いの開発環境はWSL2一択です。Virtual Box経由で仮想マシンを使う時代は終わりました（まだ不安定な挙動があるのは拒めないですが…）

管理者権限でPower Shellを開いて以下のコマンドを一発たたいて再起動しましょう。

$ wsl --install

ディストリビューションを指定することもできるみたいです（デフォルトではUbuntu 20.04 LTSが入る）。

インストール可能なディストリビューションを表示

$ wsl --list --online
インストールできる有効なディストリビューションの一覧を次に示します。
'wsl --install -d <Distro>' を使用してインストールします。

NAME            FRIENDLY NAME
Ubuntu          Ubuntu
Debian          Debian GNU/Linux
kali-linux      Kali Linux Rolling
openSUSE-42     openSUSE Leap 42
SLES-12         SUSE Linux Enterprise Server v12
Ubuntu-16.04    Ubuntu 16.04 LTS
Ubuntu-18.04    Ubuntu 18.04 LTS
Ubuntu-20.04    Ubuntu 20.04 LTS

ディストリビューションを指定してインストール

$ wsl --install -d <Distribution Name>

参考：https://docs.microsoft.com/ja-jp/windows/wsl/install

ターミナル

長いことマルチプラットフォームに対応しているHyperを使ってきたわけですが、Electron製なためか、やや動作が重かったので、正式リリースとなったタイミングでWindows Terminalに乗り換えています。

Windows Terminalについて

インストールはMicrosoft Storeから。

見た目の設定

Ctrl + ,から開ける「設定」、もしくはC:\Users\owner\AppData\Local\Packages\Microsoft.WindowsTerminal_8wekyb3d8bbwe\LocalState以下に保存されているsetting.jsonをいじることで、諸々の設定を変更することができます。

まずはカラースキーマを変えましょう。個人的にちょっと薄目が好きなのでOne Half Darkにします。
参考：https://docs.microsoft.com/ja-jp/windows/terminal/customize-settings/color-schemes#included-color-schemes

次にtransparentを設定してスケスケにします。公式ドキュメントによると、"useAcrylic": falseのtransparentはWindows11でのみ使用可能です。

こんな感じでスケスケターミナルが使えます

起動時ディレクトリの変更

WSLの起動時ディレクトリはWSL側のホームディレクトリのほうが何かと便利。

            {
                "colorScheme": "Campbell",
                "commandline": "wsl.exe -d Ubuntu",
                "guid": "{c0a7939e-6455-4394-9958-f2f923b093f1}",
                "icon": "ms-appx:///ProfileIcons/{9acb9455-ca41-5af7-950f-6bca1bc9722f}.png",
                "name": "Ubuntu",
                "startingDirectory": "/home/kimoton"
            }

startingDirectoryをホームディレクトリに変えてあげます。

フォント

後続のfishテーマでPowerline記号とNerd Fontが必要になるため、対応したフォントをインストールします。今回は、大変美しく可読性の高い等幅フォント、白源 を使わせていただきます。

github.com

インストールはこちらから。

シェル

補完が優秀かつラグも少ないため、fishが好きです。

fish is a smart and user-friendly command line
shell for Linux, macOS, and the rest of the family.

fishはスマートでユーザーフレンドリーなコマンドラインシェルらしいです。

fishに関するパッケージマネージャとしてoh-my-fishを、機能拡張のためpeco、z、bobthefishを入れます。この辺は好みで選定しましょう。

$ curl -L https://get.oh-my.fish | fish
$ omf install peco
$ omf install z
$ omf install bobthefish

Dockerのインストール

昨今のシステム開発において、Dockerはもはや必須となってきています。
ちょっと前まで、Windows上でDockerコンテナを動作させるにはHyper-Vを無効化してVirtualBoxの仮想化ソフトウェアを使用するみたいなことが必要だったのですが、MSとDocker社の提携、WSL2の登場により随分と楽になりました（Docker on Windowsの遍歴についてはこちらが詳しそう）。

以下から、Docker Desktopをインストールして、起動するだけです。
Developers - Docker

なぜかインストールしても起動できない（dockerアイコンが出てこない）事象に苛まれましたが、どうやらWSLのバグのようです。Power Shellから再起動してあげると慌てて立ち上がってくれます。
参考：[SOLVED] Docker Failed to Start - Docker Desktop for Windows - Docker Desktop for Windows - Docker Community Forums

統合開発環境

私基本ターミナルで作業するのだけれど、特にファイルを行ったり来たりするフロントエンドの開発なんかではVSCodeを使います。

以下からよしなにインストールしましょう。
Visual Studio Code – コード エディター | Microsoft Azure

VSCodeからWSL環境をいじる

WSL環境の操作にVSCodeを使うには、Remote - WSLを入れる必要があります。

細かい説明はこちらの記事を参考にしてください。

注意点として、VSCodeの拡張機能は、各WSLディストリビューションに個別にインストールする必要があります。

Python環境

minicondaのインストール

もう1年前の話ですが、Anacondaが「大規模な」商用利用では有償になりました。

個人で使う分にはAnacondaを入れても問題ないようですが、シンプルにサイズがでかいのでMinicondaを使っています。
こちらも以下からよしなにインストールしてください。
Miniconda — Conda documentation

自動整形ツール（black）

以前の記事でPythonの自動整形ツールを比較しました。

今やPythonの自動整形ツールとしては、blackが一般的になってきていると感じます。blackにはvimのpluginが存在するので、Python使い＆Vimmerの皆さんはぜひ入れておきましょう。

vim-plugでインストールする場合、~/.vimrcに以下を追加して:PlugInstallしましょう。

call plug#begin()
" Python black lint
Plug 'psf/black', { 'branch': 'main' }
call plug#end()

さらに以下を追加してあげると、ファイル保存の度にBlackを実行することができます。

autocmd BufWritePre *.py execute ':Black'

PEP8に準拠したimport文

import文にまで気を使いたいあなたはisortを使いましょう。 isortに関してもvimのpluginが存在するので入れておきましょう。

Plug 'fisadev/vim-isort'

:Isortと叩くだけで並び替えてくれます。便利～

血液由来DNAと唾液由来DNAのWGSサンプルの品質評価

巷ではやり始めている遺伝子検査キット。このサービスの大半は唾液を送って解析しますよね。唾液からDNAが取れるのはなぜかというと、白血球が大量に含まれているためだそうです（下記参照）。

遺伝子検査の試料として、なぜ唾液が使えるのか。一般的な試料として血液や、口腔の粘膜、髪の毛などが挙げられますが、試料の条件として体を大きく傷つけることなく、遺伝情報を担うDNAを豊富に含む細胞を採取できるところがポイントになります。 核を持つ白血球を豊富に含む血液は検査の精度を考えると好ましいのですが、注射などで血液を取り出す必要があるため、大量に手軽に採取できません。 この血液と同じくらい検査に向いていると考えられているのが唾液なのです。国際的な遺伝情報の研究でも標準的に使われることが多くなっています。一見透明な唾液ですが、多くの白血球を含んでいます。この白血球の細胞にDNAが含まれていて、遺伝子検査の解析対象になります。流動性があって、液体で均一であるところも試料として好ましいのです。 （引用：https://mycode.jp/howitworks.html）

変異検出の観点からみて、血液と唾液の差がどの程度なのか。
今回はそんな論文を要約してみました。

bmcmedgenomics.biomedcentral.com

PS. 最近、論文要約用のフォーマットを色々お試し中です。ちょろちょろ変わるかと思いますが、迷走中なんだな、と温かい目で見守って頂ければと思います。

この論文を読む目的は？

  • WGSにおいて、血液と唾液の品質にはどの程度の差があるのか知る
  • サンプルの品質評価方法について知る

どんな論文か？

厳格な QC 基準を満たす高品質の唾液サンプルは、血液由来の DNA が利用できない場合、または適切でない場合に WGS に適用可能であることを示した。唾液 DNA は、SNV においては血液DNAと大差ない検出が可能だが、CNVの検出量は低くなる。

使用されている技術や手法は？

Illumina HiSeq X（30x, 150bp）
↓
Isaac Alignerを用いてhg19に対してマッピング
↓
Isaac Variant Callerを用いて変異検出（SNV・short INDEL）、Control-FREEC/Canvas を用いて変異検出（CNV）
↓
American College of Medical Genetics and Genomics (ACMG)を用いて病原性をアノテーション

口腔内マイクロバイオームへのマッピングにはHuman Oral Microbiome Database (HOMD) を利用

リサーチギャップは何か？

血液と唾液のWGSサンプルの品質について、比較検証した研究はまだ存在しなかった。
唾液由来 DNA が血液由来 DNA の代替案として機能するか否かについて、唾液由来サンプルの有用性を検証した点が新しい。

どのように検証しているか？

研究で使用した5サンプルのシーケンスカバレッジ（引用：https://bmcmedgenomics.biomedcentral.com/articles/10.1186/s12920-020-0664-7/tables/2）

ヒトリファレンスゲノムおよびヒト口腔マイクロバイオームへのマッピング結果（引用：https://bmcmedgenomics.biomedcentral.com/articles/10.1186/s12920-020-0664-7）

  • アガロースゲル状のバンド、DNA濃度（20 ng/μl）、吸光度比260/280（1.3）を指標としてDNAの品質を検証
    • 46%の唾液由来のサンプル、6%血液由来のサンプルがWGS の品質管理 (QC) 要件に適合しなかった。
  • 血液サンプルと唾液サンプルについて、ヒト以外のリファレンスゲノムにマッピングされるリードの割合を比較
    • 唾液由来のサンプルでは、血液由来のサンプルと比較して口腔内マイクロバイオームの配列が取得される傾向があることを確認
    • 人リファレンスゲノムを用いて非ヒト配列を除去することが可能なことを検証
  • 血液サンプルと唾液サンプルをついて、SNVとCNVの一致率を調査
    • 心血管疾患関連遺伝子コーディング領域及びすべてのコーディング領域内の変異において、唾液由来のサンプルから特定されたSNVの95%以上が血液由来のサンプルで特定されたSNVと一致していた。
    • MAF < 0.01 で定義される希少なSNVについては、血液・唾液間の一致率が90％低下した。CNV は、血液と唾液のサンプル間で 76% しか一致しなかった。

追加の議論は？

  • 唾液は長期保存が可能で、常温で輸送できるため、血液に比べて輸送中のサンプルの保存が容易なため、遠隔の参加者からのサンプル収集において有用
  • 唾液サンプルからの DNA 品質のばらつきが大きいことを考えると、唾液採取のための適切な技術を確保し、DNA 品質に厳格な QC 基準を適用することが重要
  • CNVの検出において、リードデプスが使用されているため、唾液サンプルの CNV 収量が低いのは、唾液サンプルのカバレッジが減少した結果である可能性がある。
  • CNV検出におけるベストプラクティスなパイプラインが欠如していることもこの一因と考えられる。

PackageCompiler.jl を使った初回実行高速化

Juliaは初回ロード時にJIT（Just In Time）コンパイルを行っているため、巨大なパッケージをロードする場合、実行に時間がかかってしまいます。2回目以降の実行ではコンパイル済みなのでC並の速度が出ますが、ファイルを実行してデバッグするようなスタイルを取ると実行の度にJuliaセッションが切り替わってしまうため、毎回毎回コンパイルに時間をかけてしまうことになります。

sysimageと呼ばれるセッション保存機能を使うと、JITコンパイル済みの状態を保存することが可能となります。このsysimageを利用することにより、sysimageに含まれるパッケージや関数についてコンパイルの必要がなくなるため、パッケージをリロードして再コンパイルするよりも高速に動作させることが可能となります。

Juliaには起動時にデフォルトで使用されるsysimageが付属しており、Baseという標準ライブラリ等最低限の要素はコンパイル済みで用意されています。（参考：Building the Julia system image）

sysimageをカスタムで作成したい！そんな願いを叶えてくれるのが、今回ご紹介するPackageCompiler.jl です。 github.com

sysimageの作成コマンド

sysimageの作成には、PackageCompiler.create_sysimage()を使います。

function create_sysimage(packages::Union{Symbol, Vector{Symbol}}=Symbol[];
                         sysimage_path::Union{String,Nothing}=nothing,
                         project::String=dirname(active_project()),
                         precompile_execution_file::Union{String, Vector{String}}=String[],
                         precompile_statements_file::Union{String, Vector{String}}=String[],
                         incremental::Bool=true,
                         filter_stdlibs=false,
                         replace_default::Bool=false,
                         cpu_target::String=NATIVE_CPU_TARGET,
                         script::Union{Nothing, String}=nothing,
                         base_sysimage::Union{Nothing, String}=nothing,
                         isapp::Bool=false)

覚えておきたい引数リスト

引数 値
packages パッケージ名のSymbolをリストで指定
sysimage_path 出力するsysimageのパス
precompile_execution_file プリコンパイル対象の関数実行を記載したファイルを指定（後述）
precompile_statements_file コンパイル命令が記載されたファイルを指定（後述）
incremental 追記するか、新規作成するか。デフォルトのsysimageを上書いてしまうと、REPLやPkgのプリコンパイルがされないため、特別な場合を除いてデフォルトのtrueでよい
replace_default デフォルトのsysimageを置き換えるか否か。デフォルトのsysimageに依存しているパッケージが存在するため、基本的には非推奨（後述）

デフォルトのsysimage置き換え

replace_default=true引数により、デフォルトのsysimageを置き換えることも可能ですが、Julia-VSCode等のパッケージがデフォルトのsysimageに依存しているため、非推奨とされています。

なお、初回のデフォルトsysimage置き換え時にバックアップが作成されるため、上書いてしまった場合でもrestore_default_sysimage（）コマンドで元のsysimageに戻すことができます。

julia> create_sysimage([:Debugger, :OhMyREPL]; replace_default=true)

※デフォルトのsysimageパスは以下で確認できます。

julia> unsafe_string(Base.JLOptions().image_file)
"/usr/local/src/julia-1.5.2/lib/julia/sys.so"

パッケージのロードを高速化

Exampleパッケージを例にとってコンパイルの行い方を見ていきます。 単純にパッケージをコンパイルしたい場合は、パッケージ名とsysimageの出力先（sysimage_path）だけ指定してあげれば良いです。

julia> using PackageCompiler

# Example パッケージについてExampleSysimage.soという名前のsysimageを作成
julia> PackageCompiler.create_sysimage([:Example], sysimage_path="ExampleSysimage.so")
[ Info: PackageCompiler: creating system image object file, this might take a while...

作成したsysimageを使用してJuliaを起動する際は -J, --sysimage で指定します。

$ julia -JExampleSysimage.so

# ロードが早くなる！
julia> using Example

関数のロードを高速化

PackageCompilerでは、パッケージのロードに加え、内部の関数についても事前にコンパイルすることが可能です。関数のコンパイルには、以下の2つの方法を取ることができます。

  1. 関数実行のスクリプトをprecompile_execution_fileとして指定する
  2. コンパイル命令をprecompile_statements_fileとして指定する

1. 関数実行のスクリプトを使用

例として、precompile_example.jlという以下のスクリプトを作成します。

$ cat precompile_example.jl
using Example
Example.hello("friend")

precompile_execution_fileに作成したスクリプトを指定してsysimageを作成することで、関数を含めてコンパイルすることができます。

julia> using PackageCompiler

# Example パッケージ及びprecompile_example.jlに含まれる関数についてExampleSysimagePrecompile.soという名前のsysimageを作成
julia> PackageCompiler.create_sysimage(:Example; sysimage_path="ExampleSysimagePrecompile.so",
                                         precompile_execution_file="precompile_example.jl")
[ Info: PackageCompiler: creating system image object file, this might take a while...

2. コンパイル命令を使用

trace-compileの出力を利用

Juliaを--trace-compile = file.jl引数を付けて実行すると、Juliaプロセスの実行中にコンパイルされた関数のコンパイル命令をfile.jlに出力することができます。 コンパイル対象のコードが記述しづらい場合（インタラクティブシェルの起動をトリガーに呼ばれる関数等）に、この引数を用いた方法は便利です。
OhMyREPLの例が公式にあるので参考にどうぞ。

なお、--trace-compileの出力は標準出力にも出せます。どのようなコンパイルが行われているかサクッと確認したい場合は標準出力で確認するのが良いでしょう。

$ julia --trace-compile=stderr -e 'import Example' --startup-file=no                                  (base)
precompile(Tuple{typeof(fzf_jll.__init__)})
precompile(Tuple{getfield(Core, Symbol("#Type##kw")), NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.FreeBSD, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{typeof(JLLWrappers.get_julia_libpaths)})
precompile(Tuple{typeof(OhMyREPL.__init__)})

テストスイートを利用

以下のようなテストスイートをコンパイル追跡対象とすることで、使用される関数について網羅的なコンパイル文を生成することができます。

import Example
include(joinpath(pkgdir(Example), "test", "runtests.jl"))

テストスイートのコンパイル文をcompile_example.jlに出力

$ julia --trace-compile=compile_example.jl -e 'import Example; include(joinpath(pkgdir(Example), "test", "run
tests.jl"))'

コンパイル文からsysimageを生成

julia> using PackageCompiler

# Example パッケージ及びcompile_example.jlに含まれるコンパイル文についてExampleSysimagePrecompile.soという名前のsysimageを作成
julia> PackageCompiler.create_sysimage(:Example; sysimage_path="ExampleSysimagePrecompile.so",
                                       precompile_statements_file="compile_example.jl")
[ Info: PackageCompiler: creating system image object file, this might take a while...

Example パッケージを用いたベンチマーク

sysimageを作成

# ~/.julia/config/startup.jl を読まないようにする
$ julia -q --startup-file=no

julia> using PackageCompiler

# 仮想環境を作り、Exampleパッケージをインストール
(@v1.5) pkg> activate .
 Activating new environment at `~/test/Project.toml`
 
(test) pkg> add Example
   Updating registry at `~/.julia/registries/General`
  Resolving package versions...
  Installed Example ─ v0.5.3
Updating `~/test/Project.toml`
  [7876af07] + Example v0.5.3
Updating `~/test/Manifest.toml`
  [7876af07] + Example v0.5.3

# "ExampleSysimage.so"というファイル名でsysimageファイルを作成
julia> create_sysimage(:Example; sysimage_path="ExampleSysimage.so")
[ Info: PackageCompiler: creating system image object file, this might take a while...

julia> exit()

sysimageを利用しない場合

# sysimageファイルを指定せずREPL起動
$ julia -q --startup-file=no

# 普通に読み込んでみると1sec以上かかる
julia> @time using Example
  1.203322 seconds (652.89 k allocations: 37.429 MiB, 5.15% gc time)

sysimageを利用した場合

# -Jで作成したsysimageファイルを指定してREPL起動
$ julia -q --startup-file=no -JExampleSysimage.so

# 起動直後からload済みとなっている
julia> Base.loaded_modules
Dict{Base.PkgId,Module} with 34 entries:
...
  Example [7876af07-990d-54b4-ab0e-23690620f79a]          => Example
...

# めちゃ早い
julia> @time using Example
  0.000259 seconds (203 allocations: 12.609 KiB)

その他

毎回決まったsysimageを読み込みたい！

juliaコマンドのエイリアスを作りましょう。

$ alias julia = 'julia -J /path/to/sysimage.so'

Rosalindを解く - フィボナッチ数列

rosalind Julia Python bioinformatics

Rosalindの紹介はこちらから。

www.kimoton.com

本日は自然界でも見られる数列、フィボナッチ数列についての実装を見ていきます。

rosalind.info

生物学的知識のおさらい

参考：相補性 (分子生物学) - Wikipedia

中世で最も才能があったと評価されるイタリアの数学者、レオナルド＝フィボナッチは、ウサギの個体群の繁殖に関する問題を考案しました。

  1. 個体数は最初の月に生まれたばかりのウサギのペアが存在する時点からカウントする。
  2. ウサギは1ヶ月後に生殖年齢に達するとする。
  3. 生殖年齢のすべてのウサギは、生殖年齢の別のウサギと交尾する。
  4. 2匹のウサギが交尾してからちょうど1か月後に、オス1匹とメス1匹のウサギが生まれる。
  5. ウサギが死んだり、繁殖を停止したりすることは考慮しない。

このとき、1年間に何組のウサギが残るか。 上記の問題の答えとしては、1年後、ウサギの個体数は144ペアとなります。

ウサギの不死には異議を立てられてしまいますが、捕食者のいない環境での繁殖ではこのモデルは当てはまるとされており、現に花びらの数や花や実に現れる螺旋の数ははフィボナッチ数であることとなることが知られています。

どの月にも、前月に生きていたウサギと新しい子孫が含まれます。すなわち、任意の月における子孫の数が2か月前に生きていたウサギの数に等しくなります。

結果として、 $F_{n} = F_{n-1} + F_{n-2}$（$F_1 = F_2 = 1$）という漸化式によって定義されるフィボナッチ数列を得ることになります。 なお、この数列にはフィボナッチの名前が付いていますが、2千年以上前にインドの数学者に知られていたようです。

問題

「$n\leq40$ の $n$、$k\leq5$ の $k$ が与えられる。1ペアから始めて、各世代において繁殖年齢のペアのウサギが $k$ ペアのウサギの子孫を残すとしたとき、$n$ か月後に存在するウサギのペアの総数を求めよ」

Given: Positive integers n≤40 and k≤5.
Return: The total number of rabbit pairs that will be present after n months, if we begin with 1 pair and in each generation, every pair of reproduction-age rabbits produces a litter of k rabbit pairs (instead of only 1 pair).

具体例と方針

具体例として、 $n=5$, $k=3$ の場合のウサギの個体数増加を下記に示します。

f:id:kimoppy126:20201127090952j:plain

任意の月における子孫の数が2か月前に生きていたウサギの数の $k$ 倍に等しくなります。

これは、漸化式で書くと下記のようになります。

$$ F_{n} = F_{n-1} + k \times F_{n-2} $$

この漸化式における $n$ 番目の値を求めればよいわけです。

フィボナッチ数列の算出方法には計算量の異なる様々な方法があります。今回は動的計画法を用いた実装を見ていきます。

解 (Julia)

Juliaでフィボナッチ数列を実装している良い例があったのでこちらを参考にさせて頂きました。

qiita.com

Juliaでは多重ディスパッチを定義することができます。

多重ディスパッチ（英: Multiple dispatch）またはマルチメソッド（英: Multimethods）は、多重定義された関数やメソッドからそこで呼び出されるべき1つの定義を選出し実行する（ディスパッチする）際に、2個以上の複数の引数が関与してどれかひとつを選ぶこと（特殊化）がおこなわれるものである。

ポイントは引数の判定が呼び出し時に行われる点です。
オーバーロード（継承）では、メソッド引数の定義型を基に静的にコンパイル時にメソッドの呼び出しが決定されますが、多重ディスパッチではこれが呼び出し時に行われるわけです。

この特性を活用したフィボナッチ配列の関数を作成します。

0. ディレクトリ構成

ディレクトリ構成は下記とします。

$ tree 
.
├── main.jl # 実行ファイル
├── data 
│   └── rosalind_fib.txt # 入力データ
└── utils
    └── file_handler.jl # 実行ファイルで使用するモジュール

1. ファイルを読むためのutil関数を作る

DeliminatedFiles パッケージを使用します。readdlm関数では、ファイル名、デリミタ、型の順で指定してあげます。

using DelimitedFiles

function read_file(file)
    contents = readdlm(file, ' ', Int)
    return contents
end

2. メイン関数を作成する

動的計画法で求めます。

動的計画法とは、メモ化を活用しながら計算することで処理時間を短縮するアルゴリズムを指します。

動的計画法は、計算機科学の分野において、アルゴリズムの分類の1つである。対象となる問題を複数の部分問題に分割し、部分問題の計算結果を記録しながら解いていく手法を総称してこう呼ぶ。

毎回計算するとフィボナッチ数列の演算は $O(2^N)$ かかってしまいますが、この動的計画法を用いることで計算量は$O(N)$となります。 

#!/usr/bin/env julia

include("utils/file_handler.jl")

# 多層ディスパッチで使用するInteger型を引数にとるメソッドを定義
function fib(n::Integer, k::Integer)
    # メモ初期化
    d = Dict(1=>big"1", 2=>big"1")
    fib(n, d, k)
end

# 引数が3つの場合は下記の関数に渡される
function fib(n, d, k)
    # メモ内に値があれば返す
    if haskey(d, n)
        return d[n]
    end
    # n-1, n-2項目の値からn項目の値を求める
    result = fib(n-1, d, k) + k * fib(n-2, d, k)
    # メモ化
    d[n] = result
    return result
end

n, k = read_file("data/rosalind_fib.txt")
println(fib(n, k))

Tips
Juliaではmethodsという便利関数があり、これを使うことで定義されているすべてのメソッドとその引数タイプを表示してくれます。

> methods(fib)
# 2 methods for generic function "fib":
[1] fib(n::Integer) in Main at main.jl:4
[2] fib(n, d) in Main at main.jl:11

解 (Python)

いつも通り、Pythonの実装に関しても載せておきます。

0. ディレクトリ構成

ディレクトリ構成は下記とします。

$ tree 
.
├── main.py # 実行ファイル
└── data 
    └── rosalind_fib.txt # 入力データ

1. ファイルを読むためのutil関数を作る

with句を使ってファイルを読み書きすると安全にファイルを開けるため便利です。

#!/usr/bin/env python

import os
from typing import List

def read_file(file_name: str, dir_name: str="data") -> List:
    file_path = os.path.join(dir_name, file_name)
    with open(file_path, 'r') as rf:
        contents = rf.read()
    return contents

2. メイン関数を作成する

Juliaの場合と大枠同様です。Pythonでは外部ライブラリを使用しない限りは多層ディスパッチを用いることはできないため、 n=1, n=2のときの値は条件分岐で入れてあげる必要があります。

#!/usr/bin/env python

from utils.file_handler import read_file

# メモの初期化
memo = {}
def fib(n, k):
    # メモ内に値があれば返す
    if n in memo:
        return memo[n]
    if n == 1:
        # n=1, n=2 の場合は1を返す
        result = 1
    elif n == 2:
        # n=1, n=2 の場合は1を返す
        result = 1
    else:
        # それ以外の場合は漸化式に従い、n-1, n-2項目の値からn項目の値を求める
        result = fib(n-1, k) + k * fib(n-2, k)
    # メモ化
    memo[n] = result
    return result

n, k = read_file("rosalind_fib.txt").split()
print(fib(5, 3))

実行

$ cat data/rosalind_revc.txt
5 3

$ python main.py
## 19

コード例

下記で公開しているので参考にしてみて下さい。

repl.it

（おまけ）バイナリ公式を使用するver

通常のフィボナッチ数列を求める場合下記のような公式を使うことで計算量を短くすることができます。

$$ \begin{eqnarray*} F_{2n} &=& (2F_{n-1} + F_n) F_n \\ F_{2n+1} &=& F_{n}^2 + F_{n+1}^2 \end{eqnarray*} $$

こちらを使うと、計算量は$O(logN)$となります。

#!/usr/bin/env julia

# 多層ディスパッチで使用するInteger型のメソッドを定義
function fib(n::Integer)
    # メモ初期化
    d = Dict(zero(n)=>big"0", one(n)=>big"1")
    fib(n, d)
end

# 引数が2つの場合は下記の関数に渡される
function fib(n, d)
    # メモ内に値があれば返す
    if haskey(d, n)
        return d[n]
    end
    # バイナリ公式を利用
    m = n ÷ 2
    result = if iseven(n)
        (2 * fib(m - 1, d) + fib(m, d)) * fib(m, d)
    else
        fib(m, d) ^ 2 + fib(m + 1, d) ^ 2
    end
    # メモ化
    d[n] = result
    return result
end

println(fib(20))

生存時間データへのDeep Learningの適用 - DeepSurv

論文要約 生存時間分析 Deep Learning Python

生存時間データの分析に関してちょこちょこ取り上げていますが、今回はそんな生存時間データにDeep Learningを適用してみた論文、DeepSurv論文を読んでまとめてみました。

bmcmedresmethodol.biomedcentral.com

1分で理解するDeepSurv

  • Cox比例ハザードモデルのような標準的な生存モデルでは、個人レベルでの治療の相互作用をモデル化するために、広範な事前の医学的知識が必要となる。
  • 通常のCox比例ハザードモデルでは、共変数の線形性を仮定していることから、個人レベルでの治療の相互作用をモデル化することが困難であった。Neural NetworkやSurvival Forestといった非線形手法は、原理上高次元の相互作用項をモデル化できるが、効果的な治療推奨システムとしての有用性はまだ示されていない。
  • DeepSurvは古典的なCox比例ハザードモデルの拡張として、生存時間データに適用可能なDeep Learning手法。DeepSurvを用いると患者の共変量と治療効果の相互作用をモデル化することが可能。さらに、患者の特徴と様々な治療の有効性との関係をモデル化することで、治療推奨に有用であることを示した。
  • 実際の臨床研究データを使用して、個別化された治療の推奨によって各データセットにおける患者の生存期間がどれだけ延びるかを検証した。治療を受けるように推奨した群は、治療を受けないように推奨した群に対して有意に生存時間が延長することが示された。

マテメソ

Cox比例ハザードモデルの復習

Cox比例ハザードモデルにおいて、ハザード関数はベースラインハザード関数 $\lambda_{0}(t)$とリスク関数 $r(x)=e^h(x)$ の積で表されるのでした。 $$ \lambda(t | x) = \lambda_{0}(t) \cdot e^{h(x)} $$ Cox回帰を実行するには、重み $\beta$ を調整してCox比例ハザードモデルの部分尤度関数を最適化（最大化）します。 部分尤度関数 $L_{c}(\beta)$ は、$T_i$: イベント発生時間、 $E_i$: イベントの発生有無、$x_i$: ベースライン時の共変量としたとき下記の式で表されます。 $$ L_{c}(\beta) = \underset{i : E_{i} = 1}{\prod} \frac{\hat{r}_{\beta}(x_{i})} { \underset{j \in \Re(T_{i})}{\sum} \hat{r}_{\beta}(x_{j})} = \underset{i : E_{i} = 1}{\prod} \frac{\exp (\hat{h}_{\beta}(x_{i}))} { \underset{j \in \Re(T_{i})}{\sum} \exp (\hat{h}_{\beta}(x_{j}))}, $$

モデルの構成

モデル自体は至ってシンプルな順伝播型ニューラルネットワークとなっています。
隠れ層は全結合層（Fully-Connected Layer）とそれに続くドロップアウト層（Dropout Layer）で構成されています。最終層は線形結合層となっており、出力は対数リスク関数$\hat{h}_{\theta(x)}$となっています。

なお、ドロップアウト層は、全結合層のノード群と出力層の間の接続の一部をランダムに切断することで過学習を防ぐ働きを担っています。
ネットワークに関するハイパーパラメータ（隠れ層の数、各層におけるノード数、ドロップアウト率等）はRandom Searchによって決定しています。

ここで、目的関数は負の対数部分尤度を平均化した値にL2ノルムを追加した値を使用しています（下記）。

$$ {}l(\theta) \!:=\! - \frac{1}{N_{E=1}} \sum_{i: E_{i} = 1} \left(\hat{h}_{\theta}(x_{i}) - \log \sum_{j \in \Re(T_{i})} e^{\hat{h}_{\theta}(x_{j})} \right) + \lambda \cdot ||\theta||^{2}_{2}, $$

治療推奨システムへの応用

すべての患者をn個の治療グループ $\tau \in \{0,1、…、n-1\}$ の1つに割り当てます。 各治療 $\tau = i$ を行った患者は独立したリスク関数 $e^{h_{i}(x)}$ を持つと仮定すると、$\tau = i$ の患者のハザード関数は下記で表されます。 $$ \lambda(t; x | \tau = i) = \lambda_{0}(t) \cdot e^{h_{i}(x)}. $$ 対数ハザードの差を、推奨関数として定義します。 $$ \begin{aligned} {}\text{rec}_{ij}(x) &= \log \left(\frac{\lambda(t;x | \tau = i)} {\lambda(t; x | \tau = j)} \right) = \log \left(\frac{\lambda_{0}(t) \cdot e^{h_{i}(x)}}{\lambda_{0}(t) \cdot e^{h_{j}(x)}} \right) \\ &= h_{i}(x) - h_{j}(x). \end{aligned} $$ 推奨関数 $\text{rec}_{ij}(x)$ が正となる場合、治療 $i$ は治療 $j$ よりも高い死亡リスクにつながります。したがって、患者は治療 $j$ を処方されるべきであると判断することができます。 一方、推奨関数 $\text{rec}_{ij}(x)$ が負となる場合、治療 $i$ は治療 $j$ よりも効果的であるとみなすことができます。したがって、患者は治療 $i$ を処方されるべきであると判断することができます。 なお、Cox比例ハザードモデルにおいて推奨関数を考えると、対数リスク関数の線形制約からモデルが患者の特徴に関係なく定数を返すことが分かります。 $$ \begin{aligned} \text{rec}_{ij}(x) &= \log \left(\frac{\lambda(t;x | \tau = i)} {\lambda(t; x | \tau = j)} \right) \\ &= \log \left(\frac{\lambda_{0}(t) \cdot e^{\beta_{0} i + \beta_{1} x_{1} + \ldots + \beta_{n} x_{n}}}{\lambda_{0}(t) \cdot e^{\beta_{0} j + \beta_{1} x_{1} + \ldots + \beta_{n} x_{n}}} \right) \\ &= \log \left(e^{\beta_{0} i + \beta_{1} x_{1} + \ldots + \beta_{n} x_{n} - (\beta_{0} j + \beta_{1} x_{1} + \ldots + \beta_{n} x_{n})} \right) \\ &= \beta_{0} i - \beta_{0} j \\ &= \beta_{0} (i-j). \end{aligned} $$

これは解釈として、モデルが重み $β_0$ を正または負と推定するかどうかに基づいて、すべての患者に同じ治療オプションの選択を推奨することを意味します。

治療（患者）ごとに異なるリスク関数を推定する特性から、個別化された治療推奨システムへの応用はDeep Learningベースの手法に分があるといえます。

結果

本論文では、大きく下記の2点を主張しています。

  1. （非線形性を持った対数リスク関数に対する）予測精度の高さ
  2. 治療推奨システムとしての有用性

予測精度の高さ

精度指標にはC-Indexを使用し、CPH（Cox比例ハザード）、RSF（Random Survival Forest）と比較しています。
シミュレーションにより生成した線形データへの当てはまりは極わずかにCox比例ハザードモデルに負けているものの、それ以外のデータセットではCox比例ハザードモデルの精度を上回っており、多くのデータセットでRSFと同等、またはそれ以上の精度が出ていることが示されています。

予測精度をシミュレーションにより検証

共変量を2変数としたときの対数リスク関数をシミュレーションすることでその予測精度の高さを示しています。

  1. リスク関数に線形性を仮定
    共変量 $x_0$, $x_1$、実際のリスク関数が共変量の線形関数（$h(x) = x_0 + 2x_1$）で表されると仮定したときの各手法で予測されるリスク値、誤差を表示しています。 CPH（Cox比例ハザード）、DeepSurv共に精度よく推定できていることが分かります。

  2. リスク関数に非線形性を仮定
    共変量 $x_0$, $x_1$、実際のリスク関数が共変量の非線形関数（$h(x) = \log (\lambda_{\max}) : \exp \left(-{\frac{x_{0}^{2} + x_{1}^{2}}{2 r^{2}}} \right)$）で表されると仮定したときの各手法で予測されるリスク値、誤差を表示しています。 当たり前ですが、CPH（Cox比例ハザード）では全く推定できていないことが分かります。
    一方、DeepSurvでは精度よく推定できていることが分かります。

治療推奨システムとしての有用性

シミュレーションデータ（非線形）、臨床データセットの2セットにおいて検証されています。
有用性の指標としては、推奨グループと非推奨グループにおける生存関数の中央値、およびログランク検定におけるp-valueを使用しています。RSF（Random Survival Forest）と比較しています。

治療推奨のグループ間差異を検証

  1. リスク関数に非線形性を仮定したデータセットを使用
    データセット内のシミュレートされた各患者に対して、治療グループ $\tau\in{0,1}$ を均一に割り当てた後、治療グループ $\tau = 1$ にはガウス分布に従う非線形の治療効果を与えています。
    その後、推奨グループと非推奨グループの両方について、各モデルによって推定されたカプラン・マイヤー生存曲線をプロットしています。ログランク検定の結果から、DeepSurvでは有意な差が得られた一方、RSFでは有意性が得られなかったことが分かります。

  2. 臨床データセットを使用
    続いて、治療の推奨、非推奨に関するデータを含む実際の臨床データセット（Rotterdam & German Breast Cancer Study Group (GBSG) ）を使用します。推奨グループと非推奨グループの両方について、各モデルによって推定されたカプラン・マイヤー生存曲線をプロットしています。
    DeepSurv、RSFのいずれにおいても有意な差が得られています。一方で、RSFのp-valueの方が高く、治療推奨システムとしての有用性はDeepSurvの方が高いことが予想されます。

DeepSurvを試してみる

DeepSurvはPythonのパッケージとして、GitHub上で提供されています。

※インストールの際には、h5pyのインストール及び最新版のLasagneのインストールが必要でした。

$ pip install --upgrade https://github.com/Lasagne/Lasagne/archive/master.zip
$ pip install h5py

上記のインストールが完了したらDeepSurvのインストールに移ります。

$ git clone https://github.com/jaredleekatzman/DeepSurv.git
$ cd DeepSurv
$ pip install .

./example_data.csvにテスト用のデータセットが公開されています。
こちらを使用してテストランを行います（テストランのノートブックはGitHub上に公開されています）。

train_dataset_fp = './example_data.csv'
train_df = pd.read_csv(train_dataset_fp)
train_df.head()

##    Variable_1   Variable_2  Variable_3  Variable_4  Event  Time
## 0            0           3           2         4.6      1    43
## 1            0           2           0         1.6      0    52
## 2            0           3           0         3.5      1    73
## 3            0           3           1         5.1      0    51
## 4            0           2           0         1.7      0    51

{
    'x': numpy array of float32
    'e': numpy array of int32
    't': numpy array of float32
    'hr': (optional) numpy array of float32
}

pandasのDFをDeepSurvの入力フォーマットに変換するdataframe_to_deepsurv_ds関数を作成・実行します。

# event_col にはイベントの発生有無を示す列を指定 
# time_col にはイベント発生時間を示す列を指定
def dataframe_to_deepsurv_ds(df, event_col = 'Event', time_col = 'Time'):
    # 各データをnumpy arrayに変換
    e = df[event_col].values.astype(np.int32)
    t = df[time_col].values.astype(np.float32)
    x_df = df.drop([event_col, time_col], axis = 1)
    x = x_df.values.astype(np.float32)
    
    # DeepSurvの入力フォーマットとして返す
    return {
        'x' : x,
        'e' : e,
        't' : t
    }

train_data = dataframe_to_deepsurv_ds(train_df, event_col = 'Event', time_col= 'Time')

各ハイパーパラメータは辞書型で定義します。
本来であればデータセットごとにチューニングする必要がありますが、今回はテストランのため適当に与えます。

# ハイパーパラメータを定義
hyperparams = {
    'L2_reg': 10.0,
    'batch_norm': True,
    'dropout': 0.4,
    'hidden_layers_sizes': [25, 25],
    'learning_rate': 1e-05,
    'lr_decay': 0.001,
    'momentum': 0.9,
    'n_in': train_data['x'].shape[1],
    'standardize': True
}

DeepSurvインスタンスを作成し、.trainメソッドにより学習を行います。
学習の過程では、損失関数 (loss) の値と C-Index (ci)の値が表示されます。

# ハイパーパラメータを与えてDeepSurvインスタンスを作成
model = deepsurv.DeepSurv(**hyperparams)

# DeepSurv は学習データと検証データの追跡にTensorBoardを使用できます。

# ログを出力したくない場合は下記の3行をコメントアウトした上でlogger = Noneを指定します。
# logger = None

experiment_name = 'test_experiment_sebastian'
logdir = './logs/tensorboard/'
logger = TensorboardLogger(experiment_name, logdir=logdir)

# 学習
update_fn=lasagne.updates.nesterov_momentum # 使用する最適化手法を選択 \
                                            # 参照： http://lasagne.readthedocs.io/en/latest/modules/updates.html \
n_epochs = 2000

# If you have validation data, you can add it as the second parameter to the function
metrics = model.train(train_data, n_epochs=n_epochs, logger=logger, update_fn=update_fn)

## 2020-11-22 12:36:03,215 - Training step 0/2000    |                         | - loss: 26.2085 - ci: 0.3693
## 2020-11-22 12:36:32,527 - Training step 250/2000  |***                      | - loss: 13.5646 - ci: 0.3617
## 2020-11-22 12:37:02,550 - Training step 500/2000  |******                   | - loss: 8.7107 - ci: 0.3730
## 2020-11-22 12:37:32,982 - Training step 750/2000  |*********                | - loss: 7.0760 - ci: 0.3919
## 2020-11-22 12:38:03,217 - Training step 1000/2000 |************             | - loss: 6.4698 - ci: 0.4375
## 2020-11-22 12:38:33,590 - Training step 1250/2000 |***************          | - loss: 6.2745 - ci: 0.5634
## 2020-11-22 12:39:04,227 - Training step 1500/2000 |******************       | - loss: 6.1954 - ci: 0.6801
## 2020-11-22 12:39:36,684 - Training step 1750/2000 |*********************    | - loss: 6.1768 - ci: 0.7048
## 2020-11-22 12:40:07,556 - Finished Training with 2000 iterations in 245.75s

学習の過程（学習曲線、検証曲線）は .plotメソッドにより描画することが可能です。
学習中に表示される各種精度指標は.trainメソッドの返り値（metrics）に含まれています。

# 最終的なC-Index指標を出力
print('Train C-Index:', metrics['c-index'][-1])
# print('Valid C-Index: ',metrics['valid_c-index'][-1])

# 学習曲線、検証曲線を描画
viz.plot_log(metrics)

## Train C-Index: (1999, 0.695762840452166)

終わりに

METABRICのデータセットでは、患者の臨床的特徴（ホルモン治療指標、放射線療法指標、化学療法指標、ER陽性指標、診断時年齢）に加え、遺伝子発現量（MKI67, EGFR, PGR, ERBB2）のデータを使用していました。
こういう分析はワクワクしますね。

本論文によりDeep Learningの生存時間データへの適用可能性が示されたことにより、 時系列の医療画像を入力データとするようなことも可能となるのではないでしょうか（もうあるかも）。

一方で、Deep Learningとなるとその説明性・解釈性は皆無となっており、あくまで「結果としてこれだけの精度が出せたので臨床応用が期待できる」といった書きぶりになってしまっています。実際にこのようなリコメンドシステムが医療機関で利用されるようになるには、精度のみならずその精度の説明性が担保される必要があると感じました。