NONMEM

Uppsala way

$ERROR
IPRED = F
W = IPRED * THETA(.)
Y = IPRED + W * EPS(1)

$SIGMA
1 FIX

個体内変動モデルのUppsala流の記述。ヨーロッパ帰りのコントロールファイルというと、通常これをさす(笑)

数値計算としてのメリットはなく、THETA()として標準偏差あるいはCVのレベルで個体内変動が出力できるというわかりやすさがあるくらい。

が、この記載方法は iWRESを正確に求めようとすると必須になる。日頃やっているとあまり気づかない。

| | コメント (0) | トラックバック (0)
|

わかったこと

時期をみることというのは大切だ。

NONMEM のversion 7 は2009年の10月にリリースされた。Distributionは順次行われ、おそらく日本の多くの会社が手にしたのは2010年のはずだ。version は、7.1.0 というリリースとなった。

NONMEM7 のアピールポイントが何だったか。データアイテムが50まで可能になったりと、いろいろ細かい点はあったが、基本的にはさまざまな新規推定アルゴリズムが利用可能になったところがもっとも大きかったと思う。

私は思った。NONMEMはそんなに新しいことを追わないでもっと枯れたソフトであってもよかったのではないかと。そしてライセンスとかお金にこだわらずにずっと生きていけばよかったんじゃないかなぁと。

なぜか。新規アルゴリズムは別に、Sheinerが考案したわけではなく、90年代の後半に盛んになった非線形混合効果モデル解析のアルゴリズム実装の中で出てきたものをimplementしただけにみえたから。

だって、SAEMなんてMONOLIXが先じゃないかと。MCMCだってWinBUGSとかみたいに先手を取ったソフトがあると。

それを後追いで実装してどうする? しないよりはましだけど。と思ったのだ。

今となってみると、私はこれは正解だったと思う。

今回、Phoenix NLMEを使ってみて感じたのは、想像以上にプロ向けのソフトウェアだったことだ。これはNONMEMにとっては無視しがたいライバルである。

おそらくだ。日本の企業で行うような解析であれば、NLMEでカバーできることがほとんどである。NONMEMのような、一種のベアボーンなプログラムを要するような解析は、欧米ではされているが、日本ではどうだろう。

NONMEMでやりやすいのは、相互フィードバックがかかってしまうような複雑なモデルを、部分分割してモデリングしていくときに、かなりうまくコーディングできるところだ。

そうでないようなモデリングは、ほぼNLMEが匹敵している。

こういう状態になっている2013年を考えたときに、仮にNONMEMが線形近似をベースにした推定法のみを実装しているのだとすれば、NONMEMにとってはかなり危険な状態だったはずだ。

なぜなら、今。SAEMをはじめとする新規推定法は選択肢の一つだからだ。経験が少ないだけであって、モデリングをするときには候補推定法の一つに入っている。当局も、経験が少ないからSAEMがだめとは言わないわけだ。なぜなら、統計学的に優れた特性を持っているからである。

そして、だからといって線形近似が全否定されたわけではない。というのは、非線形混合効果モデリングは難しいモデリングであって、データによってどういう推定法が効率的にパフォーマンスするかが異なっている。FOCEのような条件付き推定法は収束困難性はあっても、驚くほどのパフォーマンスをみせてきたのがこれまでの歴史であり、inconsistencyはモデルそのものの複雑さを考えると、ある程度やむなしというところがあるわけである。

よって、手持ちのデータをどう解析するかという点については、推定法におけるヒェーラルキーがあるわけではない。そのデータに対してうまく結果をもたらしてくれるアルゴリズムがよい、と考えられることが多いのである。

こういう状況を考えたときに、あの2009年でNONMEMが新しい推定法を実装しておき、改善を続けるという路線をとったことは、ソフトウェアの寿命という意味ではやはり得策だったのだろう。2012年の今から始めたのでは遅すぎるからだ。

健全な競争があることは望ましい。非線形混合効果モデルの世界がより豊かになってくれればそれが一番だ。そのためには、NONMEMはこれからも重大な役割を果たしえるのだろう。

| | コメント (0) | トラックバック (0)
|

MAP推定

EBE推定のアルゴリズムはBFGS法であることが、ICONから初めてコメントされた。ようやく、やっていることのイメージがつくようになってきたと思う。

| | コメント (0) | トラックバック (0)
|

7.3 beta

ベータ版のインストール。ちょっとコントロールファイルの書き方で変わったところがあるというのでそこを試すのと、線形近似でないアルゴリズムがどうなっているか。

| | コメント (0) | トラックバック (0)
|

すごい

Robert Bauerはすごい。彼の記述はNONMEMに親しんでいる者にとってはかなりわかりやすい。Yaning Wang とともに、NONMEMの透明性を上げたといっていいだろう。もちろんBauerは今は開発者なのだから当然といえばそうだけれど、記述のわかりやすさはそれだけが原因じゃないな。

| | コメント (0) | トラックバック (0)
|

NONMEMの難しいところ

非線形混合効果モデルの難しいところとも一致するのだけど、NONMEMが何をしているのかを理解するのは大変である。

まずクラシカルな方法における推定について。NONMEMは公表当時から厳密には今にいたるまで、拡張最小二乗法 Extended Least Squares Method(ELS)を用いているとされており、User's Guideにもそう書かれている。だが現在、拡張最小二乗法という呼称で呼ばれることは海外ではなく、日本の研究会などで時折使用されるだけになっている。ELSの目的関数は最尤法の原理から求めることが可能であり、現在NONMEMの推定法は最尤法を用いているとみなされる。事実、NONMEMの目的関数は、定数項を除いた-2×対数尤度となっている。定数項というのは、データの数に応じた定数 nLog(2pi) の部分であり、規格化されていない尤度を用いていると考えればよい。

最尤法を用いるにも非線形モデルでは尤度が厳密には計算できないところが泣き所である。そのために近似が必要となる。NONMEMの推定については、この近似法についての議論がながらく焦点にあったといってもいい。現実問題として計算機性能の観点から、FO法が20年程度用いられてきたが、当然のことながら不適切であった。この線形近似を用いた推定法とをいかに理解するかというところは、NONMEMや母集団解析の初学者にとっては高いハードルであった。FO法の欠点は早くから統計家たちによって指摘されてきていたが、2000年代の初頭にUlrika Wahlbyによって出された2報の論文によってプラクティカルな意味での欠点が明らかとなり、エンドユーザーにとってもデフォルトの推定法の変更を余儀なくされたところがきっかけであっただろう。NONMEMの推定法についての理解はその後もいくばくかの謎を残していたが、2007年のYaning Wangの論文によって数学的な側面が明らかにされている。また、同年に発表されたAndrew Hookerの残差に関する論文も、NONMEMにおける線形近似に対する理解を促進した。おおまかにいって、2000年から2007年までのあいだで、推定法の持っている特徴がプラクティカルな意味でも明らかになったといえよう。

数値計算的な側面については尤度の近似の話ほどにはdocumentationがされていない。これについては限定的な情報しか手に入らないが、Bab Bauerらの論文の中にはNONMEMが原理的にはQuasi-Newton法による推定を用いていることが記されている。また、これについてはPharsightの技術者らもコメントしている。Quasi-Newton法、準ニュートン法とも呼ばれるが、これは基本的には推定において目的関数の1階微分しか利用しない。が、パラメータの変化方向の計算には「擬似的な」ヘッセ行列が用いられるという点でわかりづらい。これは推定の当初は単位行列からスタートするが、逐次的に更新され、徐々に目的関数のパラメータによる2階微分行列に「収束する」という特徴を持っている。これを採用している理由は、最急降下方向への探索が必ずしも計算時間を最短にしないことが理由である。これは2階微分行列を用いる単純なNewton法と、1階微分飲みを用いる最急降下法の折衷的な方法である。NONMEMはBFGS法を用いている。

パラメータ推定がうまくいかなくなると、つまりパラメータを変化させることによって目的関数を低下させることができにくい状態になると、NONMEMはその疑似ヘッセ行列を一度リセットする。これが、RESET HESSIANという推定中に表示されるエラーメッセージであり、このときに、NONMEMは最急降下方向を選択する。

このヘッセ行列は対数尤度についての2階微分を用いていないことに注意したい。推定中のHESSIANは、あくまでも逐次更新される補助行列である。

FOCE法においてはこの1回1回のITERATIONの中で、POSTHOC推定値を必要とする。Intermediateモデルは当然未熟なので、POSTHOC推定値を求める計算は失敗したり不適切な値が求まることがある。これはFOCE法の持っている収束不安定性につながっている。この作業は母集団パラメータを既知とした上で、それについての罰則をかけた状態の最小二乗法が用いられているが、Documentationが少ない部分である。知られている限りでは、NONMEMのUser's Guide以外には、Sheinerが1982年と最初期に出した論文で言及されているくらいである。これについては1999年に、当時GlobomaxであったMarc Gastonduayが行ったプレゼンテーションに他のBayesアプローチなどとともに紹介されている。

一度推定がなされると、COV STEPと呼ばれる段階でパラメータ推定の標準誤差が求められる。ここの考えは最尤法の考え方であり、目的関数の2階微分をベースにしている。が、これについてもサンドイッチ推定量という概念が導入されている。これは1991年にWhiteが出した文献にあるが、変量効果に正規性が仮定できない場合にも頑健である。元来、ここで必要とされるFisher情報行列には、漸近収束する2つの計算方法が存在し、それがHessianを用いるものと、1階微分係数ベクトルをかけ算するものがある。前者はNONMEMではR行列と呼ばれ、後者はS行列と呼ばれる。NONMEMはR-1SR-1を用いる。また、推定中とは異なり、ここでは数値的に2階微分を計算するので、計算時間が必要になる。

FIMについては教科書における定義と異なるのもポイントであろう。厳密にはFIMの計算においては、平均操作が必要となる。が、そのプロセスを行わない、observedなFIMというものが存在する。NONMEMはこのIntegrationは行わず、observed FIMを利用する。また、これは推定値近傍で計算されるがその推定はデータから行われているので、Empirical FIMとも呼ばれる。$COVではR行列などを出力可能だが、NONMEMは目的関数の2階微分を用いているため、R行列を出力するときには0.5倍される。S行列の場合は0.25倍である。

NONMEMのいくつかのエラーはこのあたりから生じてくるので、イメージが重要である。が、これらは決してwell-documentedではないため情報の収集は難しい。

もちろん、何も知らずともNONMEMを使うことはできるのだが。だがこれらについて理解することはNONMEM学習のときの初期の大きなハードルの一つでもあるので、一応まとめ書きしてみた。

| | コメント (0) | トラックバック (0)
|

NONMEM 7.3.0

β版の使用ができるようになった。どの辺が改良されたのか?

| | コメント (0) | トラックバック (0)
|

コントロールファイルの流儀

コントロールファイルの書き方は実際のところは目的に従って動けばいいのだから、どういう書き方でもいい。

が、確かに組織によって流儀がある。セミナーとかで教えるときは、できるだけ一般的な書き方にするように留意している。

私のいるところでは一番目立つのは個体内変動の誤差の部分の書き方であるが、あの書き方は存外難しい。IWRESの部分も。よって、決して初学者向けではない。

(が、思いっきりここの書き方を出してしまったことも一度ある)

ほかにも、出自がわかるような流儀はけっこうある。

なお、南半球の大御所は、「こういう書き方をしている人たちもいるけど、ほとんど意味ないね。昔あったコントロールファイルが、まだVirusのように(!)生きてるだけだよ」とか、NMUserに書いたりするので、侮れない(笑)。

まぁ、書き方で出自がわかるというのも面白いといえば面白いのだろう。

| | コメント (0) | トラックバック (0)
|

NONPARA

NONMEMで$NONPARAを使ったことがない。関連論文も読んだことがないので、この2-3日、勉強している。別に使い方が難しい類のものではない。また、パフォーマンススタディの論文もよくわかる。

が。。。何をNONMEMがやっているかが今ひとつわからない(´Д` )

自分の研究に使う予定はないから、もう少し寝かせる(=放っておく)か?

NONMEMのノンパラメトリック解析は、ETAの分布に関して正規性を仮定しないというもの。これは2段階で実施されていて、最初に従来法で解析しEBEも求める。その後、求められたEBEをsupport pointとして。。。 ノンパラの解析に入るのだけれども、ここが今ひとつ。

これが学会で盛んなころ、私はむしろカテゴリカルデータ解析に関心を持っていて、そういう本ばかり読んでいた。だから、NONPARAはあまり知らないのだ。だいたい2007年くらいの話だ。

よし、寝かせよう(´∀`)

| | コメント (0) | トラックバック (0)
|

気づく

午前中ミーティングで、どうしてもうまくいかないコントロールファイル上の問題についてのアドバイスを受ける。

最初理解できなかったのだけれど、それを理解しようとしてマニュアルやNMUsersをひっくり返したら、わかるようになった。つまり、NONMEMがどんなふうに動いているのかについて具体的なイメージを持つことが重要。

NONMEM依存の問題という部分もあるけど、結局同じような課題はPhoenix だろうがMONOLIX だろうがかならず発生する。というか、おそらくやろうとしていることは、NONMEMでないと不可能なのかもしれない。

午後の4時間くらいを使い、異常に疲れた。。。

| | コメント (0) | トラックバック (0)
|

より以前の記事一覧