本質的なランダム性とそうでないもの

本質的にランダムな現象というのは実際には少ないと思う。

コイントスでもサイコロでもいい。表か裏かは1/2の確率。1の目がでる確率は1/6。コイントスを20回続けて、全部表だったとする。次は裏だろう、と予測しても、次に表が出る確率も裏が出る確率も1/2。これはランダムな現象の本質である。

が、これをランダムたらしめているものは、観測不可能な情報というものだ。手の上のコインの初期配置、手の動きとコインに与えられる力。これらが厳密に測定された場合、コイントスはただの力学の問題になるはずだ。だから、表か裏かは高い確度で予測可能なはずである。

でも、「実際のコイントス」をするときに、これらの情報が入手できることはありえない。絶対にありえないのだ。それをやっていたら、コイントスではないのである。

言いたいことは、決定論的に予測可能であっても観測不可能な情報がある場合は確率論的にふるまっているということだ。こういうものは本質的なランダム性があるわけではない。

こう考えると、実は本質的にランダムなものというのはほとんどないのではないだろうか。知人との雑談の中で、量子力学的な不確実性よりも上のスケールでのランダム性ってどこまであるのかが話題になったことがあるが、自然科学で扱える現象の場合、それはないはずだ。

なので、手に入るデータの示すランダム性は、情報が不充分であるからこそのランダム性なのである。

これはPharmacometricsをやっている人なら容易に理解できることだし、既に理解しているだろう。我々はプレーンなモデルに対して、共変量という、結果を予測しやすくする因子を選択して含めていく。体重とか腎機能とか。結果に影響をもたらしているだろう因子を。これらは観測可能な情報であって、臨床試験でも当然評価可能である。観測された情報を含めると、ランダム性は減るのだ。

しかし、完全に説明しきれるものではない。PKであったって、ある被験者の肝や腎の血流量、酵素の発現量などPKに由来するクリティカルな情報は評価されない。だから、観測されれば当然よりデータを説明できるだろうものが観測されてない分、個体間変動というランダム性を完全に除去することなどできはしない。

そういう意味で、随所で繰り返し主張されているように、モデルはもちろん正しくはないのである。

より困難なケースでは、何が観測されたらランダム性が下げられるのか、見当がつかない状況というものもありえる。PKよりはPDのほうが、決定的な因子が何であるのかを基礎科学の知見から見抜くことが難しい。

また、同じ理由で経済に関するモデルにおいてはPharmacometricsのモデルよりは難しいことが多い。

ランダム性をある一定以上に下げることは難しいし、そしてそれは現実的でもない。臨床試験で評価不能な因子を、無理に評価しようとするとそれは臨床試験という枠組みでできなくなるわけだ。コイントスと同じである。だから、一定のランダム性を常に覚悟する必要がある。

それでもモデルが有用なのは、そのランダム性が我慢できる程度であるときだ。つまりこのくらいの予測のズレは判断に影響を及ぼさない、と言えるときなのである。その文脈で、シミュレーションが重要なのである。

充分なモデルができても、そこからの予測が役に立たなくなる原因は２つある。1つは短期試験から長期試験のような外挿である。厳密にはモデルが役に立たないのではなく、手元にあるデータそのものに病態の自然経過による情報が含まれていないだけだ。だからこれは観測情報の議論になる。ここをうめようとしているのがDisease model。

もう１つは、臨床試験にエントリーする被験者が、違ったふるまいを示す場合である。モデルから外挿できる範囲を越えている。次にやった試験では、モデルを構築した試験の被験者では影響力がなかった因子が突如として強い影響を示すようになったとか。もしくは新たな影響要因が含まれるようになったとか。

ただ、これも実際はモデルの問題ではない。こういうことが起きないように、試験は毎回同じような患者層からエントリーするものなのだ。患者のPopulationは常にオーバーラップしていて、その外側の人たちは入ってくるもののメジャーではないというやり方になっているはずである。

統計モデルを構築することは、こういう観測不可能なランダム性を覚悟しつつ、そうでないランダム性を削れれば削っていくことである。

説明がハイブローすぎるのもね

確率微分方程式を勉強しよう！と思って本を開いたり、webをあさるとする。そうすると、ブラウン運動（まぁこれはよし）、ルベーグ積分（統計やるならしゃあない）、伊藤積分（うほ）。。。と続く。高尚だよね、説明が。これがLipsitz条件などのもとで、一意な解を持つというのは確かにすごいと思うのだけど、そういう説明からがんがん行くというのも、説明が高尚すぎるような気がする。

ランジュバンとかブラック・ショールズとか、応用から行ったほうがいいね、

ここは(´∀｀)

確率過程自体が、自然なイメージを持っているものだから、そこから行きましょう。

Huber推定量

NONMEMでいうところの、Rinv*S*Rinv。サンドイッチ分散だが別名Huber推定量と呼ばれている。

モデルがincorrectな場合の、分散の推定値についての頑健化である。

原報はどれだろうと思っていたのだが、どうも以下らしい（？）

Huber, P. J. (1967). “The Behavior of Maximum Likelihood Estimates under Nonstandard Conditions,”Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability,　vol. I, pp. 221–33.

Levene検定みたいに、Proceedingsで報告されているということなのだろうか(この検定も、和文の本では1冊くらいしか言及されてなかったと思う)。

重み付き残差

今日の講義で重み付き残差について触れた。統計のベーシックな講義。

回帰分析では、weighted residualsというよりは、そもそもステューデント化した残差を用いる。それは標準化のための計算に、ハット行列の対角成分による補正が入ることと、誤差分散が未知であるため残差平均平方で代用するものだ。しかも残差平均平方については、該当データを抜いた回帰結果を用いる。

こうすると、残差の分布がt分布に従うので理論的に扱いやすくなる。標準的な回帰分析ではこういう細かい工夫がされている。こういう工夫が可能だということだ。回帰分析の診断については恐れ入るほどの検討がなされている。

もちろん、今日の初学者向けの講義ではこういうところは省く。残差はそのままでは診断に使えないことがあるよ、標準化がいるよ、という内容だ。

非線形混合効果モデルまでいけば状況は知ってのとおりである。近似が入らざるをえないから、t分布を用いるという感覚はなくなっていく。そもそも、自由度が計算できない。

自由度の計算が困難なのは、誤差構造を階層にした時点で既に発生する問題だし。線形混合効果モデルで既に生じる問題だ。だから、satterthwaiteとかkenward-rogersとか、怪しげな近似式が登場しているのである。

90分くらいで回帰分析から線形、非線形混合効果モデルまでを、それほど数理を用いず、エッセンスとモデリングにおけるPracticalな話として講義できれば、けっこう価値はあるだろうなぁと感じた。これは非常に難しい。「ほら、こうクリアしたでしょ」みたいな、スキップの妙みたいなものが必要かもしれない。

私は統計家ではないから、こういう話をしていいものか。出自や自分の被教育歴について迷うことが多いのだが、結局は学んだことを伝えるということを真面目にやるしかない。それ以外の答えは結局のところない。

でも、なんで私は統計が好きなんだろ？　面白いからだけど、言葉では説明できない。

統計の研究者を私はいつも尊敬している。

正規分布の歌

なんと世の中には正規分布の歌なるものがあるようだ。

試しにGoogleで正規分布の歌を検索してみると、初音ミクが歌っているものを目にするはずだ。

Ax=b

教育資料のための行列に関する資料づくりを少し。連立方程式についての見方について、まとめる。

係数行列A、未知なる解ベクトルx、定数ベクトルbとしたときに、Ax=bに解があるかどうかは、Aに逆行列が存在するか否かというところにかかっている。が、もう少し見方を変えれば、係数行列Aを構成する列ベクトルの張る部分空間にbが存在するかという見方ができる。

Aを構成する列ベクトルは決して空間全体をカバーできるだけのものでないことがありえる。Ax=0なる同時方程式系を考えた時に、Aのkernelを考えることができる。これがtrivialな解のみを持つ場合は、kernelは原点の1点のみである。一方そうでないときは、kernelは次元を持つ。つまり空間を持っている。このような場合、Aがnxn行列、つまり元の空間がRnであったとしても変換によって構成される空間の次元が落ちる。次元定理は有名だろう。よって、bがその部分空間に存在しない可能性がでてくるということだ。

これを簡単に判定するには、Aを構成する列ベクトルが線形独立であるかを確認すればよく、rankを調べればよい。rankは行列に対して基本変形を行えばよい。基本変形によってrankは変わらないからである。

どうも書き方がぎこちないな。

ちなみにrankをRで確かめるには、qr関数でqr分解し、$rankを確認すればよい。

統計を学ぶための。。。

英語教育研究に携わっているらしい人が、「卒論修論のために、これから統計に取り組む人のための五ヶ条」という記事をブログに書いている。生物系や経済系のように体系的な統計のトレーニングを受けなかった人のためとのこと。生物系や経済系に限らず、理系全般でも統計を体系的にトレーニングを受けられるところは少ないと思うが＾＾　書き方は鼻につくのだが、言ってることが間違ってるわけではないのでまとめてみた。

1. 先輩から教えてもらわない

2. 「論文」で学ばない

3. 「統計ができる」という保証がある人・本から学ぶ

4. （最低限の）数式の暗記・理解は、覚悟する

5. 最初は統計ソフトを使わない

この五ヶ条であるらしい。

1は、もしあなたが統計学初心者であるならば、その先輩が統計を理解しているかどうかを判別することはできないということだ。また、理解している可能性は極めて低いと。周囲の評判といっても、周囲も理解していないのだからあてにならない。だから学んではならないと。「統計を使って研究している人」が「統計を理解している人」とは必ずしも限らないということを言っている。ソフトウェアも発達し、クリック統計学派が存在するから。

2は、先行研究と類似の手法を用いたとはいえ統計がわかる人が査読していることはなかなかない。また、紙幅の関係で統計手法の使い方の吟味について議論されていることはなかなかないということである。

3は結局のところ教科書がよいということ。

4は、数式を一切使わない統計学的な本は読まないほうがいいということ。定理とその証明のレベルで統計学を理解することは、ユーザーにとってはさすがに必要ではないと。ただ、使い方とそこにおかれている仮定、前提の理解には数式が必須であるということだ。

5は、結局のところ統計ソフトはかなりブラックボックスであると。Rといえでもクリックしないだけだである。最近のソフトは融通が利きすぎていて、データの一部に間違いがあってもとにかく数値ははじき出す。基本的な統計はエクセル表計算でプリセット関数の使用なしでやれということだ。

1については会社内でも部長クラスの上司から致命的な誤りを含んだ統計解析の指示がなされたときが幾度となくあったため、身につまされることしきりだ。

2はPharmacometricsにおいては、あまりあてはまらないだろう。

3は正しいが、「正統な教科書」が何であるか見抜けるということは初学者では既にないことを意味する。初学者が適切な教科書を選択できるかというと、微妙かもしれない。

4はその通り。混合効果モデルも基礎理論を紙とエンピツで追うことは必須である。

5もまぁその通りではあるが、Pharmacometricsをエクセルで実行することはできない。VPCをエクセルで実行したケースがあり、それはまた驚嘆であるが＾＾　理論の勉強とソフトの勉強を並行させるのが最も良い。

私は統計を理解しているか？　いやそれはあるまい。だから私に聞いてはいけない＾＾；　といいたいところだが、社内的な立場としては「わかっていないといけない位置取りにいるのだから、その責任を果たせるように全力で勉強しなさい」というところだろう。いくつになってもなかなか大変なのである。4

線形混合効果モデル

Pinheiro & Bates の MEM in S の日本語版 Chapter 2を読む。パラメータ推定のための最適化式である eq 2.10 を求めるところ、§2.2がまずは大切なところ。このあたり久しぶりなので、線形代数というか、少し忘れていて苦労した。

数理の部分は式の展開を追うことができることが重要であるが、気をつけたところは。

式(2.5)のところ。Ψ-1 は変量効果の分散共分散行列の逆行列。

これは定義からΔT・Δ/σ2　であることがわかるから、そうするとbiT*ΔT*Δ*biとなり、Δbiのノルムの2乗となる。

式(2.6)の指数項の中身。最小二乗解が存在する。それはyi-Xiβをひとまとまりとし、Ziを計画行列、biをパラメータと考えれば、変量効果の推定値についての式が得られる。要はy=Xβの最小二乗解が、b=(XT*X)-1*XT*y であることと同じことをしている。これにより変量効果の条件付きモードを推定することができる。

式(2.9)は、多変量正規分布の密度関数が全域で積分すれば1になることを利用し、スケーリング項となることを利用。これはラプラス近似にもでてくる考え方。

2時間くらいかけて3ページしか進めないような場所だ、この部分は。

検出力の講義

仮説検定と検出力、例数設計についての勉強会。きわめて初歩的なところ。
しかし、検出力の計算は特に非心分布を使うわけではない、正規分布をベースにした解説であっても初学者には大変である。

基本的に、数式を用いた展開にはついてきてほしい。R を使えば、イメージや感覚だけで検出力を計算することが可能である。が、教科書にあるような対立仮説を構成するある一つの分布において、検定においてあらかじめ設定される棄却限界値に対して右側確率を計算する。両側の場合は寄与は小さいかもしれないが一応、両側の分を計算するということだ。これについて、累積分布関数で確率を計算しているなんてよくあることだが、これもマスターしてほしい。

CDFの特徴は、積分で定義される関数ということである。それが解析的であるかないかは関係がなく、数値積分でもかまわないわけだ。この積分関数というものを難しく考えてほしくないのだが、抵抗感は確かにあるだろう。

さて、検出力の理解において何が難しいんだろう。まず分布を複数考えるところが難しいのかもしれない。検定は考え方をつかめば理解はたやすい。帰無仮説が正しいとして、検定統計量の従う確率分布から p値を算出する。確かにこれは難しくない。

一方、検出力を考えるのは難しい。帰無仮説ではないときの分布を考えないといけない。でも検定の力を考えているから、棄却限界値はそのままだ。そして、帰無仮説ではないときの分布は無数にあり、それら一つ一つに検出力が計算されうる。

決して難しくはないのだ。しかし、作業に落とし込むと意外に手間取るのがわかる。ましてや数学的に示されると、実体以上に理解が困難になる。

講義のゴールは、検出力曲線を、検出したい平均の差、ばらつき、サンプルサイズの3つに対して作成できるようになることにある。今日は座学であったが、なんとか手と理論の両方が身についてくれないかと思う。

あっていた

セミナーで自由度1の t 分布について、これをCauchy分布と呼ぶはずだといった。あっていた。

この分布、Pharmacometricsではおそらくほとんどでてこないはずである。期待値も分散も計算できない分布だ。ただ、分光学など物理学の世界で確かRolentz分布という名称でよばれていて重視されてたはず。

ということはモデルによればPharmacometricsでも重視されるはず?　たぶんそんなことはない。たしか強制振動か共鳴か何かの現象の解か何かだったはず。そういうモデルがでてくることもないでしょう。

。。。 Wikipediaにも載っているのか。東大出版会の基礎本で確認してたんだけど。