線形代数とか線形とか

今日は線形ﾓﾃﾞﾙについて行列を用いた表記での講義をした。初めて線形代数の手法を用いて教えるときは慎重になるのだが、質問もでたし納得感もあったようでよかった。ちょっと意外だった。

固有値や標準化の話をしないだけで、行列の扱いについては線形ﾓﾃﾞﾙで相当でてくる。まずはこの辺りのレベルをこなすことが目標であろうか。クロネッカー積がでてくるところくらいまで統計の理解が進めばよいのだけれど。

最近の高校数学のカリキュラムがよくわからないので、教え方に困ったりもする。基本的に、大学の一般教養の数学についてﾏﾆｱｯｸに使いたおせれば足りないということはほとんどない。公理論的ｱﾌﾟﾛｰﾁは大切だとは思うが、実際は数理研究をするわけじゃないから実用的な計算力や式変形のﾃｸﾆｯｸを持っていた方が良いと思う。

線形代数の話をするたびに思うが、線形という言葉は自然科学においては重要である。基本的に比例関係とその和を扱うと思っておけばいいだろうか。

薬学卒業生にはなかなかの割合で、線形を1次関数か否かでとらえてしまう。つまり反応変数が直線的に変化するかである。これはなかなか困る。

物理学の場合、線形といえば線形微分方程式を連想する。線形微分方程式システムには解析解が存在するので非常に扱いやすい。PKでいうMammilary Modelも線形微分方程式システムである。線形微分方程式の解は、非線形関数である。

統計学の場合は、パラメータに対する線形性である。説明変数と反応変数の関係がどれほど曲がっていようとも関係ない。y=Xbというモデルであればよい。線形ﾓﾃﾞﾙの場合、パラメータ推定が解析的に可能であり数値的な反復推定は不要となる。これは線形混合効果ﾓﾃﾞﾙでさえ例外ではない（実際は分散成分の推定のために反復推定が必要なのだが）。

ことほどさようにいろいろな線形があるのだが、こういうことを真剣に考えるのもまたよし。

周辺という言葉

日本語としてイメージしやすいかどうかは別だが、統計解析においては”周辺”という言葉が用いられるときがある。英語の marginal の訳語である。

marginal は日常語でも時には使われる「マージン」の、まぁ形容詞なわけである。marginは余白という意味もあり、例えばWordで文書を印刷するときに、上下左右に文字が印刷されないスペースがある。余白である。あれをマージンと呼ぶ。

血中濃度が上昇しても安全性上懸念がないと考えられる、血中濃度の振れ幅のことをSafety Margin、安全域と呼ぶ。これも余白の意味であろう。さすがにマージンであるだけあり、「Safety Marginの範囲のうちは懸念がない」。Wordの文書で例えれば、紙をはみ出ることはないということである。紙をはみ出たら、別な世界に行ってしまう。そういうことだろう。何年か前まで、安全域の計算で毒性が発現する濃度域に対して計算している人がいて、驚いたことがあった。これはこれで意味を持つだろうが、安全域と呼称するのはちょっとまずい。

さて、話を元に戻すと、この余白という意味は統計においてはイメージしにくい。ずばり直撃するのは、分割表における周辺度数である。確かに分割表本体のはしっこに、集計の結果としてあるいはそもそものﾃﾞｻﾞｲﾝとして記述されていて、marginalである。ここから周辺確率を計算することができる。

そこで周辺確率、marginal probabilityである。これはmarginからやや意味を転じることになる。端的に言えばこれは、「他の要因は全くおかまいなしに考えて、ある事象が起こる確率」といってしまってよい。何でもよいのだが、有効／無効というﾃﾞｰﾀがあったときに、それは当然Placeboか、Activeかによって出方は異なるだろう。細かく言えば用量によっても異なるだろうし、被験者背景によっても異なるだろう。今着目している結果に影響を与える要因はたくさんあるわけだ。

これらを全部無視して、有効である確率を考える。こうするとこれはmarginal probabilityとなる。分割表の余白にあるという意味から、他の要因を無視したという意味合いに変わってくる。この言葉遣いは正しく、「本当の姿は全ての要因を考慮した状態の、同時確率にある」からであって、それからすると、単独の結果に着目している周辺確率というのは確かに周辺にすぎない。

おもしろいのは、その真の姿まで推定できないｹｰｽがあることである。確率という言葉で上の分掌を作成したが、ﾊﾟﾗﾒｰﾀ推定の立場に立つと確率という言葉がそのまま尤度に置き換わる。

そうすると、推定できないだろう要因については全て合計してしまえということになる。周辺確率でも同じなのだが、他の要因については合計を行なう。合計は積分と言い換えてもいい。

今日はある勉強会でχ2検定、分割表の独立性の検定についての基本的な話をする機会があり、こういうことも話した。ざっとしたイメージが大切だからであるが、もちろんあることを意識しているためでもある。

残業時間しか面白いことができない

最近、コア業務時間内の平均「科学」時間は10分以下で大量の雑用を処理している。今やっている雑用には終わりがないので、夜7時から7時半を過ぎるともうやめることにしている。

8時から9時過ぎまで最尤法についての教育を行なう。題材は自作ではないものを用いたので、notationがあれ？という感じになる。自分が作ったものではない資料で説明するのはいつも難しいのだが、それは自分にとっては一つの勉強にもなるということで。

論文をまだ読む

今日はある意味どうでもいい外出で、少ししか時間がなかった。

Liang, K.Y. and Zeger, S.L. Biometrika 1986 73:13-22, Longitudinal data analysis using generalized linear models

ﾊﾟﾗﾒｰﾀ推定についての続き。

論文

今日は急用が入ったりして忙しかった。論文を一つ読んでる。

Wei Pan, Biometrika 2001, 88(3):901-906 "On the Robust variance estimator in generalised estimating equations"

通常のSandwich varianceはPROC MIXEDではEMPIRICALと指定しないと計算してくれないが、推定値はあくまでも漸近正規性を示すのでﾃﾞｰﾀ数が大きくない場合はﾛﾊﾞｽﾄ分散を計算したほうがよいと考えられる。

また、非正規分布するﾃﾞｰﾀの場合は特にﾓﾃﾞﾙ分散よりもﾛﾊﾞｽﾄ分散の使用が勧められる。

このLiang と Zegerの提唱した推定値の分散共分散行列についての改良提案の論文。

カテゴリーデータのモデリング

カテゴリーデータや計数データのモデリングは難しい。一番の難点はモデル診断である。計量値に対するモデルは残差をもとにした診断が可能であり、診断から得られる情報が多い。一方カテゴリーデータはそういう情報を得ることができない。

説明変数もカテゴリーである場合、説明変数のカテゴリーごとに情報を要約することが可能である。それは有効率のようにモデリングする対象そのものとしてまとめることができる。これは非常に良いことで、データとモデルから予測される有効率の比較から、モデルの妥当性についての検討を行なうことができる。

説明変数が連続量である場合は特にその傾向が強いのは、単独の残差が意味を持たなくなってくることだ。包括的な適合度として尤度、目的関数を評価することができるが、ある特定の部分集団についてのあてはまりなど、細かい検討をすることが難しい。

連続型説明変数のロジスティックモデルについて言及されている書籍は決して少なくない。丹後、Agrestiなど入手しやすい書籍にも記述されている。Etteの本もほぼ同様のｱﾌﾟﾛｰﾁが例としてあげられている。診断としては説明変数のbinningを勧めているｹｰｽがほとんどであった。予測される有効率をモデルからSimulationする。それを実測と比べるということ。Uppsalaの発表にもそういう stack bar chart が出てくる。

なお、binaryﾃﾞｰﾀの場合はHosmer-Lemshow検定を適合度検定として用いることができる。実際の2値ﾃﾞｰﾀと、モデルから予測される確率を大小関係に並べたときのｸﾞﾙｰﾌﾟで分割表を作ると、χ2検定を実施できる。ただし順序ｶﾃｺﾞﾘｶﾙﾃﾞｰﾀには適用できない。

使い方としてはLack of fit。

PROC LOGISTIC data=XXX;

model Y = X1 X2 X3 / selection=stepwise slentry=0.05 slstay=0.01 lackfit;

RUN;

今日

A.Dobsonの本（3版、洋書の方）がどこにいったかわからなくなって困った。探す。

昼間はZeglerとLiangの論文が届いたのでとりあえずﾌｧｲﾙ。ﾛﾊﾞｽﾄ分散についてちょっと読みたくなって。ｻﾝﾄﾞｲｯﾁ推定量について調べたかった。Verbekeも読んだり。

多忙でやりたいこともやるべきことも遅れてるんだけど、しょうがない。

最近、わりとあっさり「自分はだめだなぁ」と思うがそう深刻でもない。

JMP

せっかく使える環境なのに手を出していない。ちょっともったいない。

過去の記事を読みかえしつつREML

過去の記事を読み返す。ブログではあまり饒舌ではなくなってきたことに気づく。

今日は勉強時間がとても少なかった。帰り道、REMLについて読む。意外にわかりやすい説明に出会って、ちょっとうれしい。

今年は2009年である。このブログは2007年の夏くらいから書いている。最近諸事情でややﾌﾞﾗﾝｸがあいたが、まあ適当ないたずら書きなので今後も続けていこうかと思う。

なんとなく遠慮がちな記事が増えてきたしｺｰﾄﾞもあまり載せなくなったと思う。なぜかといえば、ｺｰﾄﾞが重要というよりは自分がやっていることの意味を把握するほうが重要だと感じるようになったから。やりたいことがはっきりしていれば、ｿﾌﾄの使い方は必ずついてくるような気がするのだ。

統計ｿﾌﾄは使いにくいわけではない。ｿﾌﾄではなく統計への理解が深まれば使いやすい。理解が浅いとHELPに何が書いてあるのかがわからないから、使えない。結局、ｼﾝﾌﾟﾙな話だ。

ロジスティック回帰

いろいろと成書をあたるが、ﾛｼﾞｽﾃｨｯｸ回帰の変数選択やﾓﾃﾞﾙ診断については知っている以上の話はでてこない。標本割合と予測確率の差の診断以外には適合度統計量での判断となっている。通常の回帰分析における留意点はおなじように留意すべきであることは確かなのだが。

しかしこの環境、すぐに閲覧できる論文がほとんどない。ﾘｱﾙﾀｲﾑに取り寄せできない。なんというか一種の苦学生時代が続く。

今日の勉強

これをひたすら読んでまとめたりしていたおとなしい一日。他の書籍と比較すると、これは読みやすい方ではないかと思う。

NONMEMというか混合効果ﾓﾃﾞﾙのことが頭にありすぎて、最尤法においてσ＾2が推定対象ではないnuisance parameterであることに妙な違和感を感じてしまい、ずっと考え込んでしまった。誤差の構造が階層的でない場合、これを推定する必要はなくてあくまでも母数の推定ができればよいことになるということか。

逆に誤差の構造が階層的だったり相関があったりする場合、測定値の分散共分散行列も複雑ななりたちをするので分散成分も推定過程に入れる必要があることになる。

統計ｿﾌﾄがどうやっているかは別のお話であって、最尤法を用いているProcedureでは誤差分散を推定している。こういうこともあって、記述にひっかかってしまったのだろう。

やりすぎたかも

単回帰分析についての基礎的なｾﾐﾅｰ。時間かかりすぎ、しゃべりすぎ。

QT延長の解析

Time-matched control をおき、時点どうしの差の信頼区間の90％片側信頼区間を評価時点で求めて全点の上限を評価する。

ﾃﾞｰﾀとしてこれとそっくりのものを持っていて解析しているのだけど、けっこう大変だなぁと思う。幸か不幸か、解析の結果はそれほど不明確でもないDataだった。しかし、SAS Codeの量といったら、予備解析と描画を入れるとなかなかなﾎﾞﾘｭｰﾑである。なんとなくむずむずするCodingなのだけど、やっていると愛着がわくのも不思議だ。その後、SAS Codeの整理について考える。ﾌﾟﾛｸﾞﾗﾑ仕様書をどうつくっておくかとか。

今日はこの解析についてﾌﾟﾚｾﾞﾝしたあとで、Rの簡単なSimulationについての講義をした。ループってこうなんだよとか、簡単なもの。ちなみに、Rの条件分岐の構文を忘れてた。

PROC MIXED data=XX;

class subject period time drug;

model Delta = period time drug period*time drug*time / solution outp=Result;

random subject;

repeated time / type=AR(1) subject=subject * period;

lsmeans time * drug / corr cov;

RUN;

明日はこんな感じ？

今日の解析

PROC MIXED。線形混合効果モデルが必須なのでそれについての対応。非線形性はおそらくなさそうなので、Linear な PK-PDﾓﾃﾞﾙを試す。というか、この場合はE-Rﾓﾃﾞﾙか。

頭を悩ませるような内容に今回も直面するが、今日はどちらかというとDataset作成とGraphばかり。あらかじめ用意していたCodeでグラフ作成するのだが、探索解析なので準備してないことも多発。統計解析はあまり充分にできなかった。

背景因子の比較もしないといけないから、Dataset作成は続きそう。解析用Datasetを作る人が別にいるという環境はうらやましい。全部自前。

でも、ﾓﾃﾞﾙが作れそうな気配はある。

ラプラス近似

p(y|x)=p(x|y)p(y)/p(x)
但し，p(x)=∫p(x|y)p(y)dy

とあるときに、p(x)を求めることが困難である。p(x)はｙについて合計した周辺尤度になるのだが、これの計算は難しいと。

このときラプラス近似を用いて計算するやり方がある。ラプラス近似というのは、その積分の値を計算するときに、yならyの最頻値が最もその積分に影響を与えていることを利用しており、よってその近傍について近似できればよいとするものだ。

ｱｲﾃﾞｨｱ自体は実はそれほど難しくはなくて、極値近辺で近似するのだから、yでの1階微分は0となる。だから、展開は2次の項まででかつ1次の項を含んでいない。

そうするとexp(-y2)の形の積分計算へ帰着できるから、積分が求まるということだ。

さて、こういうことは論文を読んだほうがよいのだが、これあたりだろう。

Pinheiro, J.C. Bates, D.M. Approximations to the Log-likelihood Function in the Nonlinear Mixed-effect Model. Journal of Computational and Graphical Statistics 1995; 4: 12-35.

2次近似までを用いている。そして2階微分の計算は数値的にタフになることが多いのだが、2階微分自体についても1階微分係数の積で表現するといった手段がとれるということがわかる。2階微係数の寄与が小さいのでこれでよいということだ。

この論文では手法の比較も行なっている。FOCE（PinheiroとBatesの）、Laplace、Adaptive Gaussian Quadratureである。

A.C. Davison Approximate predictive likelihood 1986; 73:323-332 でラプラス近似について書いてある。この論文はﾌﾘｰではない。

特徴

統計の特徴は、数理科学としては非常に例外的で一番やりたいことが、演繹的に物事を説明しつくすことではなく、帰納的に真実を推定しようとしているという点かもしれない。

本について迷う

  線形代数の（何度目かの）復習を知人が始めようとしており、自分も一緒にできないかと考える。考えて手元にある線形代数の本をぱらぱらとめくり、どれをやろうかとか迷っている。こういう数学の勉強というものは、自分のようなレベルの人間にとっては「見栄をはってはいけない」のが鉄則で、古今東西の名著というよりは、最近出版されていてわかりやすく、数学のユーザーフレンドリーであるということを重視する。そういう意味で、左の本はAmazonで最も売れている線形代数の本の一つだから、やはりこれがいいと思った。自分は毎年ゴールデンウィーク辺りにこれを復習することがなぜか多い。

この本の著者は今度「ﾌﾟﾛｸﾞﾗﾐﾝｸﾞのための確率統計」を出版するようであり、そちらもなかなか優れている。こういうわかりやすい本を書けるということは、それだけ自分で咀嚼していることになるのだが、ほんとうに尊敬する。

最近購入したのはこれ。この著者の教科書は微積分で読んでいたが、こちらもかなりわかりやすくするための工夫がある。同著者は教養課程のための本とはいえ、現在となってはかなり難しい本も執筆されているが、そちらとはまるで違う。

数学の本は手広くやる意味はあまりなく、相性のいい書籍を見つけて心中するのがよいようだ。他に線形代数だと川久保先生の本があったり、統計むきの本もあったりする。2009年の題材として迷い中。

終わりに

ちょっと疲れたので一日の終わりにこの本を勉強して、PPTにまとめ始めた。この本は原著では3版が出ているから買おうかどうか迷っている。章末問題の答えを作らないといけない。

この本のよいところは、簡潔というかﾎﾞﾘｭｰﾑが少ないため読みきることが可能なことだろう。大部な統計の本に、特に英語でｱﾀｯｸするのは大変だから、こういう書籍は非常に大切である。最低限の理論をｺﾝﾊﾟｸﾄにまとめるというのは簡単なことではないから、存在価値が高い。

今日の勉強

ほんの少しだった。MCMC。階層型ﾓﾃﾞﾙの場合に、hyperparameterが大切！とか、そういうﾚﾍﾞﾙを出ていない範囲で終わり。

optim

Rのoptim関数の使い方をまとめている。最適化手法は自分でﾌﾟﾛｸﾞﾗﾐﾝｸﾞする意義は総じて少なく練習にしかならない。よって計算ｴﾝｼﾞﾝがあればよい。統計解析も同様なのだが、最適化ｴﾝｼﾞﾝ以外を作って試すことはけっこう大切な経験のように思うので、これで最尤法くらいやってみればよいかと思う。

なお、私の持っているS-PLUSにはoptim関数はなく一変数版のoptimize関数のみが実装されている。最新versionではどうなのかはわからない。S-PLUSは好きなんだけど。

Condition Number

明確な基準が設定されないことが多いが行列の条件数が大きい場合、わずかなﾃﾞｰﾀの変化が大きな推定値の誤差をもたらすために数値解析がしにくい。条件数は行列の固有の性質であり、丸めなどそういう誤差とは関係ない。条件数が大きい問題を、ill-conditionedであるという。Bonateは10^6にﾗｲﾝを引いていた。絶対値として明確な基準はないが、計算機イプシロンの逆数以上に大きいと、まともに解いていることにはならないだろう。そういう意味で考えると、計算機イプシロンは2^-20近辺であるから（IEEE754では-23乗）、だいたいあうと思う。行列とその逆行列のノルムから計算される。

調べるかどうかについては要注意。逆行列を求めること自体は、方程式を解くことよりも計算量が多いからだ。

行列のノルムは何でもよかったように記憶しているけど、どうだっただろうか。それぞれについて条件数が計算されうるということだ。

Googleさんに聞いてみたところ、ノルムが定義されていればよいようだ。

今日の勉強

この本の中の例題をRの簡単なﾌﾟﾛｸﾞﾗﾑにしたりする。

R のループは、 for (i in 1:10){}のように記述するが、厳密には 1:10 のところにはﾘｽﾄも使えるらしい。それは知らなかった。

反復推定があるから if (判定条件) break で抜けるようにする。

例題を扱うﾌﾟﾛｸﾞﾗﾑは簡単であるが、汎用にしようとすると難しい。当たり前だけど。

準Newton法

非線形最小二乗法がおそらくたいていの人にとって最初に数値計算ｱﾙｺﾞﾘｽﾞﾑに触れる機会だろう。この辺りはなかなか難しい。

最急降下法というのが伝統的で、目的関数曲面の低いほう低いほう、Gradientベクトルのマイナス方向へひたすら降りていく。最も早く収束しそうなものだが実は曲面が谷状になっていたりすると、収束には大きな時間を要するのでイメージほどの使い勝手はなく、歴史上の存在である。

ニュートン法はヘッセ行列を計算することによって最小方向を見つけるｱﾙｺﾞﾘｽﾞﾑである。最急降下法が評価関数の1次近似を利用するとみなせば、これは2次近似である。計算が難しいが局所収束性に優れている。実用上の問題はヘッセ行列の計算が難しい問題では適用しづらいことや、ヘッセ行列の正定値が保証されるわけではないのでsaddle pointにのりあげることがあること。

準ニュートン法は、ヘッセ行列とみなせるような疑似行列を更新していくやり方である。初回推定についてはこの行列は単位行列であり、最急降下法と同じになる。2回目から、反復推定の結果を利用して行列を更新し、最終的にはヘッセ行列の近似行列が得られているというものである。最急降下法とニュートン法の中間的性質を持っていて、BFGS法が有名である。収束に困ると、というかニュートン法っぽく停滞をはじめると、ヘッセ行列をリセットして才急降下方向へのﾊﾟﾗﾒｰﾀ探索を1回実行したりするやり方もある。

最後にシンプレックス法。NelderとMeadが発展させた方法である。ﾊﾟﾗﾒｰﾀ空間で三角形を構成し目的関数を計算、最大の目的関数を示す点を残りの点の中点に対して反転させることで、目的関数の最小化を考える。この三角形をシンプレックスといい、シンプレックスの大きさが変わらないといつまでも精度が上がらないので、拡大や収縮、縮小を行なっていく。導関数計算を行なうニュートン系列の方法と比べて、単純作業でありﾌﾟﾛｸﾞﾗﾑﾘｿｰｽが小さい利点がある。このﾀﾞｳﾝﾋﾙｼﾝﾌﾟﾚｯｸｽの有用性は混合効果ﾓﾃﾞﾙにおいてはあまり検討されていない。

統計の話じゃないが、ｶﾃｺﾞﾘｰは一緒にした。

推定法とか

非線形混合効果モデルの推定法だったり、各種ソフトウェアの違いはどこが大きいかというと、ほとんどが尤度の扱い方に帰着されると考えられる。ある意味当然といえば当然なのだ。200例のPKﾃﾞｰﾀがあるとし、i.v.1-ｺﾝﾊﾟｰﾄﾒﾝﾄﾓﾃﾞﾙであるとする。そうすると素直に尤度を考えるとPKﾊﾟﾗﾒｰﾀ数×200もの推定をしなければならないことになる。

L(θ, Ω, Σ|ηi, Yobs)を考えるということになるが、各被験者における変量効果ηiは未知だから、この尤度を理論的に考えることはできても現実的ではない。そもそもそれを求めることができないからモデル解析していると考えることができるし、母集団を興味の対象とするのだから、θ、Ωが主役である。個別値を求める必要はなく、分布形がわかればよいのだ。

だからこれを、L(θ,ηi, Σ|Yobs)*L(Ω|ηi)とする。そしてηについて積分してしまう。そうすれば得られた尤度はηに依存しない。全員分、合計してしまうというイメージでもいいかもしれない。これをmarginal likelihood、周辺尤度といっている。多変量正規分布の勉強をしたことがあれば、周辺分布という言葉にはであっているから、その尤度版を想像することは容易だろう。

この積分をいかに行なうか、近似の仕方や積分の仕方が異なっている。それが非線形混合効果モデルのｵﾌﾟｼｮﾝの違いと考えればいい。

線形近似はここに関係している。PKモデルは非線形であり変量効果はその構造上、おもいっきり非線形項として存在する。正規分布を仮定したとすると、尤度はexpのべき乗の部分にさらにPK式が入っている形になり、これは相当にとんでもない。そこで対数尤度とするが、それでもなおてごわい。微分と異なり、積分は必ずしも解析的にできないからだ。学生の頃、簡単そうな積分に苦労した経験は誰しもある。微分は機械的作業なのにだ。

よってTaylor展開することになる。特にη＝0の周りで、Maclaurin展開すると積分は容易である。微分係数を求めるのは簡単なので、式はηの1次関数の積分になるから、これは単純だ。η＝0でなく、POSTHOC推定値を求めてそれを代入し、η＝ηiの周りのTaylor展開としてもよい。

さて、この近似はなかなかナイスなｱｲﾃﾞｨｱなのだが、波紋を呼ぶことになる。ηは確率変数である。確率変数の2乗というものは、分散を意味する。正確にはその期待値を分散として代表させるのだが、要はそれを無視していることと1次近似は同じであった。

そこで2次まで近似しようという話もでてくる。これは正しいことではあるが、一般に数値計算の世界では2階微分というのは扱いづらい存在であり計算時間がかかる。なので正しいとはいえ、なかなか使いがたい。

まぁ、こんなかんじでつかんでおけばよいのではないかと思う。ソフトウェアの難しいところは、この推定法というもの以外に、数値計算的なｱﾙｺﾞﾘｽﾞﾑの話もまざってくるところだ。ﾊﾟﾗﾒｰﾀ推定を反復計算で実行するが、それ自体はここで示した推定法の話とは違う。反復計算のｱﾙｺﾞﾘｽﾞﾑを理解することは、目的関数について理解しておくこととは違う。後者は解析のｺﾞｰﾙを事実上決める、統計解析の前提であるので重要なのだ。だからﾃｷｽﾄなどでは後者を扱っており、前者については記載さえないことがある。一方で線形最小二乗法や非線形最小二乗法では、推定値の求め方やｱﾙｺﾞﾘｽﾞﾑが書いてあることもある。

このあたりの記述のバランスの崩れも、混合効果モデルを理解しづらくしているのかもしれない。

今週はかなり仕事で疲れたが、やっと一息ついてこんなことを書いてしまった。自分も決して深く理解しているわけではないから、これも解説ではなくて随筆程度のものだけど。

誤差について

指数誤差モデルがあるとする。

θind=θpop×exp(η)　として、η~N(0,ω^2)　という典型的なものだ。

このとき、指数誤差モデルであれ相対誤差モデルであれ、母集団平均値が大きくなればなるほど誤差分散が大きくなるから、CV（％）が誤差の指標として使われることがある。非常に通りのよい指標だからだ。

exp(x) = 1 + x + ・・・ というように、Maclaurin展開するとexp(η)は1+ηとしてあらわすことができる。よって、NONMEMで指数誤差モデルを用いてETAを推定したとして、そのルートをとって100かければ、ηが小さいときはMaclaurin展開の妥当性が成立するからCVとみなすことができる。θindの分布の形も裾を引くというよりは左右対称な分布になると考えやすい状態だからだと。

指数誤差の場合、得られたﾃﾞｰﾀからCVを考えると有名な関係が導かれる。CV=SQRT(EXP(ω^2)-1)×100 (%)という関係である。

よってω^2が大きい場合、CVとして変動の大きさを報告するには、この式を用いればよい。

これはみんなどうしてるのかと思っていたときもあったのだが、2次情報を見つけた。Bealはこんなこと意味ないから、とにかくSQRT(ω^2)×100 (%)をCVとして表記してよいぞと言ったという。自分もそうしてきてるから困ることはないが、原典にあたるものが今ない。

こんなことでもAuthorityが言ってくれるとうれしいのだが、どういう形だったんだろ。

なお、理由はω^2、要はNONMEMのETAは推定法の影響であれ解析対象ﾃﾞｰﾀの被験者数が少ないせいであれ、母集団平均θと比べれば精度がいつだって悪い。よってこういうところは好みでこだわってもいいが、正確さが増すというよりはむしろいらんことをするということだったようだ。

細かい話（？）というか、大切な話というか。ともかく会議ばかりの今日、ささいな発見。

未読

2ヶ月くらい前に購入したのだが、まだ読んでいない。。。 もういい加減読まないと。

混合効果ﾓﾃﾞﾙの勉強

今日はこれをやっている。ﾗﾝﾀﾞﾑ係数ﾓﾃﾞﾙ。面白い冗談がある本に書いてあった。”OK. I've understood. By the way, I have one question. What's the G matrix?”。勉強はちゃんとしないと、これを笑うことはできない。

非線形混合効果ﾓﾃﾞﾙと比べると、線形混合効果ﾓﾃﾞﾙは誤差構造は複雑に指定できる。これは線形ﾓﾃﾞﾙゆえなのだが、NONMEMから入っていくとここにつまづきやすい。というか、自分もきっちりつまづいているのだが。

しかし、PROC GLM と結果が一致するあたりも実にすばらしい。構造が複雑なﾃﾞｰﾀの検定は、PROC MIXEDでやらないとかなり繁雑な手作業になることからも、格段に進歩したﾌﾟﾛｼｰｼﾞｬｰであることがわかる。

なお、慣れないうちは、GやR行列を出力しておくと何をやっているのかわかりやすいと思う。

新刊

書店に立ち寄ったときに新刊を見つけて反射的に購入する。衝動買いであるが、こういう本が出るようになったとは、時代も変わったものだと思う。

R & WinBUGS とまでタイトルに記載されている和書は、初めてであろう。まだ中身は眺めてもいない。さすがに衝動買い過ぎるという気がするが。

ついでに入門書も購入してみる。この本はAmazonでも上位に引っかかるもので、けっこう売れているようである。

中学生レベルの数学で読めるというふれこみである。統計教育というところに苦心している現状を考えると、買ってみようかと考えた。中身はほとんどみていない。こちらは新刊ではない。NONMEM を日ごろ動かして解析しているような人たちには不要な本だと思うけれど。

今日明日くらいでやること

ｸﾛｽｵｰﾊﾞｰ試験の例数設計について、SASﾌﾟﾛｸﾞﾗﾑを作ろうとしている。というかほとんどできているといえるから、説明というか解説用のPPTが一番時間がかかる。PROC MIXEDで差の信頼区間をSASﾃﾞｰﾀｾｯﾄに出力するには、ods output Diffs = out;で出力できる。これを集計すればよい。これについてはいかなるﾃﾞｰﾀをｼﾐｭﾚｰﾄするかのほうが、大切である。

M&Sの重要点は2つある。ﾓﾃﾞﾘﾝｸﾞにおいては、いかに現実のﾃﾞｰﾀから目的に応じたﾓﾃﾞﾙを作れるかというところ。これが大切かつ難しいのは誰も否定しない。ｼﾐｭﾚｰｼｮﾝはいかなるﾃﾞｰﾀを発生させるかというところ。これはあまり意識されないが、かなりの重要性を持っている。ここはﾌﾟﾛｸﾞﾗﾐﾝｸﾞﾃｸﾆｯｸが必要なのだが、それに意識が引っ張られるあまり、重要なｼﾐｭﾚｰｼｮﾝの前提への考察が弱くなってしまうからだ。

生物学的同等性

例数設計が必要になった。おそらく2×3のクロスオーバーﾃﾞｻﾞｲﾝになる（してほしいらしい）。例数設計のためのマクロの作成を始める。

オッズ比

今日やった演習で、オッズ比についての説明があまりよくなかったかもしれない。ということは自分でもまだよくわかっていないということか。

x をｶﾃｺﾞﾘｰ変数（というかダミー変数）としたロジスティック回帰について ln(p/1-p) = b0 + b1 * x + ・・・　というところで、exp(b1) がオッズ比になるという話。それほど難しくはないのだが、p/1-pがオッズであるということと、ln(p'/1-p') - ln(p/1-p) = b1 となること、よってオッズ比はexp(b1)となる。数学的に変形してそうでしょ、と話すのではなく、もう少し生き生きとイメージを伝えないといけない。

オッズ比がいささか初学者にわかりづらいのは、「比」であるからだろう。オッズそのものの絶対値に興味があるわけではなく、ある原因もしくはある介入によってどのくらいオッズを変化させたのかに興味があるというところだ。このときに、考えている対象が思ったよりも多いことがﾎﾟｲﾝﾄだろう。

オッズは、”有効というﾃﾞｰﾀが無効というﾃﾞｰﾀと比べてどのくらい得られやすいか”と言い換えることができる。この絶対値には興味がなく、これが実薬群でどのくらい上昇したかになる。よって、オッズ比を日本語でいうと、”有効というﾃﾞｰﾀが無効というﾃﾞｰﾀと比べてどのくらい得られやすいかが、ﾌﾟﾗｾﾎﾞから実薬にしたら何倍上昇したかという比”となる。オッズ比という言葉に「比」という言葉が一つしか入っていないにもかかわらず、オッズ自体が○／△という値になっているために、意外と言葉にしにくいのだ。だからこその専門用語でオッズ比として呼称する。

いろいろと調べてみたのだが、文章として簡単な必殺技はなかなか見つからない。丁寧に説明すべき概念なんだろう。

演習自体では用いなかったが、準備に対して参考にしたのはこの2冊。

Compound Symmetry

Compound Symmetry がどうしてああいう構造になっているのかがわかった。反復測定ﾃﾞｰﾀのときに、R行列の構造だけで説明しようとするとああなる。変量効果として被験者をいれてやると、ZGZ' + R = V を計算するとCompound Symmetryになっているのだが、このときはRはDiagonalでよい。わかってみるとあたりまえのように思える。これはSASで確認することもできる。R行列はrepeatedステートメントで、G行列はrandomステートメントで構造を指定する。同じ結果が得られる。

AIC

赤池の情報量基準、Akaike's Information Criteria なる AIC は非常に有名な統計量であり、多くの統計ソフトがその算出機能を持っている。算出としては、-2×Ln(尤度) + 2 ×ﾊﾟﾗﾒｰﾀ数である。

モデルは複雑になればなるほどデータへの適合性が増加するが、それは小さいノイズを無限に拾うことになりかねない。極端な話、100人から得られたある指標を説明するのにおいて、100のﾊﾟﾗﾒｰﾀがあればモデルはﾃﾞｰﾀに完全に適合するがそのモデルにはなんら普遍性がないこととなる。

このﾊﾞﾗﾝｽをとるために情報量基準が用いられていて、両者のﾊﾞﾗﾝｽをとるための指標として用いられる。情報量基準は、TIC（竹内の情報量基準）をはじめとしてほぼアルファベットの数だけあるといわれている。

全てのモデルは間違っているが、役に立つモデルがあるという Box の言葉はよく引き合いにだされるのだが、これは PPK-PDを行なうときにはよく考えておかないといけない。というのは、統計的にモデル選択ができることは大切であるが、それを意識すればするほど、モデル選択が機械的というか、オートマティックなものとなりかねない。結果を公表するというよりは社内での意思決定に用いるというときが多いこの技術では、統計的な判断以外の要素もモデル選択にとっては重要視すべきであろう。客観と主観の意味が問われるわけである。

多変数の微分

復習する。多変数の微分についてだが、けっこう穴があることがわかった。簡単な内容でももう一度丁寧に勉強しないといけないかもしれない。

多変数の微分は統計学を学ぶときの必須アイテムである。行列表現に隠されて、とっつきにくさが表面化しにくいのだが、やっていることはそういうことだ。非線形を解析的には扱わなず数値的に扱うから、あまりやっているようにみえないかもしれない。

極値という問題が重要であり、saddle point ではだめという状況が多い。よって固有値問題になったりするのだが、多変数の微分について慣れてないと理解が浅くなる。今日はちょっとへこむ。

最尤法

朝予想していたよりも意外にも一日の業務がすんなりと終了し、ちょっと時間があいた。早めに退社し、今日はベイズ統計の本はちら見にとどめて、最尤法とその講義のためのPPT資料を作った。

最尤法をNONMEMで実行するには、$EST に LAPLACE LIKELIHOODのｵﾌﾟｼｮﾝ指定が必要である。これはもう有名なことで、尤度計算の近似法としてLaplace近似といわれる方法が用いられているからだ。2階微分までを用いたTaylor展開である。理屈は難しいからとても把握できない。できている人は決して多くはないと思う。Wang の PDFﾌｧｲﾙとにらめっこしているし。

さて、最尤法の場合、$PREDに Y = ○　として記述するのは尤度になる。個々の被験者のﾃﾞｰﾀからの計算方法であるから、例えばLogistic回帰だと以下のようになる。

IF (DV .EQ. 1) THEN Y = P ELSE Y = 1-P ENDIF

NONMEMの最尤法ｵﾌﾟｼｮﾝは discrete ﾃﾞｰﾀについてのｵﾌﾟｼｮﾝでそれほど extensiveに作られてはいないとHELPにも言い訳がかいてあるが、PKについても尤度を書き込むことで推定が可能になるようだ。

ということは、Cox比例ﾊｻﾞｰﾄﾞﾓﾃﾞﾙとかもできるということになる。すごいなぁ。NONMEMのある意味での幅の広さには驚かされる。

こういうことは今日のＰＰＴにはまだほとんど書けない。それよりも自分も含めて基礎が大切という資料である。

最尤推定量は漸近的に正規分布する。そしてその漸近分散はクラメール・ラオの下限に等しいので、漸近有効推定量となる。フルモデルの尤度とﾃﾞｰﾀから計算される尤度の差は乖離度として定式化され、それはﾕﾆｰｸﾌﾟﾛﾌｧｲﾙ数とﾊﾟﾗﾒｰﾀ数から自由度が決まるχ二乗分布からﾓﾃﾞﾙのOverallの有意性が決まる。そしてﾊﾟﾗﾒｰﾀやﾓﾃﾞﾙ間に尤度比検定をﾊﾟﾗﾒｰﾀの数の差を自由度とするχ二乗分布を用いて構成できる。

これだけある意味で美しい統計解析手法である最尤法を考えたときに、NONMEMは最尤法を用いているのか？　そうみなすことができるのか？と気になる。　この質問を誰かがすると、いきなり意見がいろいろとわかれてしまう。やはりﾌﾞﾗｯｸﾎﾞｯｸｽなのだろうか。

ベイズ統計の続き

昨日に引き続き、ベイズ統計の本を読む。

導入0－3章の構成が良い本で、基礎統計を学んだ後の学習者が確率分布について復習しながらベイズ統計に導入されていくというやり方をとっている。数学的な道具立てはほとんど不要であり、ある意味平易なのだが、やはり概念を理解することは難しいものであることを改めて実感した。

ベイズ統計は最尤法という概念にさらに事前分布が加わったものとみなすことができる。事前分布という過程と、ﾃﾞｰﾀが取得された後に計算される状態になる尤度関数から、推定のために必要な事後分布が得られるという流れである。事前分布の設定や尤度の決め方に自由度があり、同じデータであってもとうぜんのことながら事後分布は設定の仕方によって異なる。ここがわかりにくいという人は多いようだ。

3章までは推論の本論に入っていない。よってベイズ推定の基礎と概念になるのだが、適宜戻ってくることになる。

統計の書籍を読むのは楽しい。しかし何度読んでも思うのは、これはほんとうは検証のためのﾂｰﾙとして存在するというよりは、決して不可知なものをどれだけ論理的に推定するかという科学であって、フィロソフィーに近いと思う。

ベイズ統計の本

初学者むけのベイズ統計の本として、Amazonでもかなり評判の良いこの本を読んでいる。確かに初学者が読むにはわかりやすく、きちんと読み込めば理解が進む印象を受ける。MCMCについての記述はまったくないが、つまりそれだけ初学者の導入のための役割に徹していると考えられよう。しばらくこれを読んでいくことにする。

Bootstrap 回帰分析

回帰分析における係数の Bootstrap 信頼区間の算出をやってみようかと考えていて、プログラムを書いてみている。回帰分析の場合、例えば10点の回帰分析だったとき、残差が10個ある。この残差について反復抽出したのちに、推定された回帰線から目的変数の値を計算する。新しく作られた目的変数の値と説明変数に対して回帰分析を行ない、回帰係数を求める。

これを1000回繰り返して回帰係数の信頼区間を求める。理屈は簡単である。母回帰直線が推定されてるのだから、新しい残差から目的変数を再計算させるときの回帰直線は、当然最初にｵﾘｼﾞﾅﾙﾃﾞｰﾀから推定されたものである。残差は独立に発生する確率変数なのだから、これをBootstrap反復抽出する。そして目的変数を計算する。母回帰直線が何かはほんとうは知らないのだから、それに対して回帰分析を行ない回帰係数を求める。これが、偶然得られるかもしれなかった回帰係数の組み合わせの一つとなる。あとはCIを求めればよい。

NONMEMのことを思い出しているのだが、BootstrapはValidationだろうか？　解析の安定性と推定の精度の算出であり、それこそ信頼区間算出が重要であるように思う。

Giltinanの本

なかなか難易度が高い本なのだが、ひさびさにこれを読んでいる。体調不良が続いたせいで今日は出社しないといけないのだが、気分転換に。他の人に聞いてみてもあまり簡単な本はやっぱりこの分野にはないから、この辺りになってしまう。

用量比例性の回帰分析評価

切片の信頼区間が0を含めばよいではないか、という感覚があるだろう。これはうまくいきそうに思えるのだが、一群例数が少なかったり例数はあっても個体間変動が大きかったときには切片の信頼区間は広がるから、曲線的な関係であっても線形という判定がなされることになってしまう。切片は外挿によって求められるﾊﾟﾗﾒｰﾀだから、信頼区間は広がりやすい。よってこのアプローチは避けている。

この説明をするのに適切な図を作っておくべきかもしれない。あと、誤って線形と判定する、False Negative がどのくらいあるかという率を見直すべきかもしれない。

回帰分析にはもう一つ誤用があって、CmaxやAUCが相対誤差もしくは指数誤差に従うため各群の等分散性が成立しない。こういうﾃﾞｰﾀのときに回帰分析を行なうのは典型的誤用である。

ﾊﾟﾜｰﾓﾃﾞﾙの優れている点は対数化することで誤差の等分散性を達成させた状態で回帰しているところである。また、用量比例性評価は切片ではなく傾きからなされるので、外挿ﾊﾟﾗﾒｰﾀではないところもよい。

では用量比例性評価におけるﾊﾟﾜｰﾓﾃﾞﾙの回帰係数をどのように評価するか、それが今後の課題である。

notation

統計の勉強で敷居が高いのは、notationへの慣れが形成しにくいこと。というか、元々けっこうめんどうくさいものを解析するのだから、表記が細かくなるのは当然である。suffixと行列表現に慣れることができるかどうかが、勉強の勢いを決めてしまう。

自分も理論的なところが決して得意なわけではない。しかしやらねばならない時間が増えてきた。立場がそうさせているともいえる。

R新刊

これも新刊。Rは出版ラッシュが続いているようで、どれもなかなか装丁のきれいな本であることが多い。スタイルが R の本にはどれもあって、プログラミング重視、統計重視などあるのだが、これはどちらかといえば後者であり人工ﾃﾞｰﾀ発生法やｻﾝﾌﾟﾙｻｲｽﾞ計算などが他書よりはよく書かれている。

SAS新刊

こちらはSASの新刊。珍しくSASの一般書籍が刊行になった。内容はver9.1を用いて執筆されているので新しく、古い本にはないods outputなどの記述もみられる。

内容はベーシックなもので、統計解析の基礎知識についても記述されていてｺｰﾃﾞｨﾝｸﾞのみの本ではない。買おうと思ったのだが、今日は買わなかった。

BEのコード 続き

生物学的同等性について、SASコードを調べたりしていた。ようやくある程度調べることができた。生物学的同等性は2×2クロスオーバー試験で評価できる。近年、より拡張したﾃﾞｻﾞｲﾝも多数考えられているが、まずこれを基本的と考えるべきであろう。デザイン自体は古く、自分自身を対照に持つという明解さから1950年代から臨床適用例がみられる。2×2クロスオーバー試験の解析方法はGrizzleが1965年に議論しており、Grizzleﾓﾃﾞﾙと呼ばれることがある。

PROC GLM data=Demo2;
class Seq ID Form Period;
model LnCmax = Seq ID(Seq) Form Period /clparm alpha=0.1;
TEST  H = Seq E = ID(Seq);
lsmeans Seq / E = ID(Seq) STDERR PDIFF;
lsmeans Form Period / STDERR PDIFF;
estimate 'T vs. R' Form 1 -1 ;
RUN;
QUIT;

これはPROC GLMによる解析である。IDが被験者であるが、これはSequenceにnestしている。B(A) は、B nested within A ということだが、ある本では上のBEコードにおいて Seq Seq(ID) となっていた。出力は変わらないのだが、これはID(Seq)の誤りじゃないかと思う。

Sequenceは事実上、群のことである。BE試験は被験者をどちらかのSequenceに割り付けるので、枝分かれ実験になっている。1次誤差で Sequence 間の変動を検定させている。製剤や時期については、2次誤差で検定している。PDIFFｵﾌﾟｼｮﾝで差の検定をしている。

検定は基本的には今のBE試験においては補助的な意味しかもたない。基本は信頼区間法であるので、estimate ｽﾃｰﾄﾒﾝﾄで計算させた。clparmｵﾌﾟｼｮﾝがないと区間推定はしてくれない。

PROC GLMを用いると分散分析表が出力されるから多くの人にとってはわかりやすい。PROC MIXED でも実行することが可能で、こちらが現在はむしろ主流である。どちらを用いてもいいのだが、中途脱落があったりするときは、PROC MIXEDの方が適している。PROC GLMの方が計算速度は速いのだが、現在はそれも気にならないだろう。

PROC MIXED と PROC GLM はComputationが違うので、あらゆる出力が一致するわけではない。PROC MIXED はREML法で計算している。また本質的ではないが、細かいｵﾌﾟｼｮﾝ指定が教科書や人によっても異なっているようだ。PROC MIXED のほうについては、まだ検討中である。

BEは意外に難しい。これはなぜかというと、一つはいまだに概念が追加されたりして発達中であること。もう一つは解析方法が時代によって変化してきていることである。国内は1980年代に、分散分析での有意差検定で同等性が主張されることを、最小検出差20％で検出力80%を保証することで担保するというやり方から、同等性を許容する範囲を決めた信頼区間法へ1997年から変わった。これは1992年のFDAのｶﾞｲﾄﾞﾗｲﾝを受けている。そして2001年にFDAがPopulation BE、Individual BEの概念を含めたBEへの統計学的ｱﾌﾟﾛｰﾁなるｶﾞｲﾄﾞﾗｲﾝを公表してきている。こういう変遷をふまえてかかれる教科書や論文があるため、学ぶことが多く敷居が高いのだろう。

分割実験

split-plot について調べる。分散分析についてもう少し深めたい。

分散分析の復習

分散分析を復習している。今日は分割実験とラテン方格なのだが、これはクロスオーバーの基本となるからである。特に分割実験は重要だと思う。テキストは昔の本だが、東大出版会の実験データの解析なるSAS本。このシリーズは読み返すとやはり優れたテキストだと思う。SASのコーディングは大きく変わらないが、ODSアウトプットの使用法など含めて、改訂してみればよいのにと思った。逆に、そう思うならば自分用のものを自分で作れば良いのだが。
最近は共分散分析なども含めて、復習を余儀なくされる事態が多い。医薬品開発において古くさい問題とはあまりなく、常に物事は古くて新しい問題となるのであろう。

非心度

対立仮説の Non-centrality parameter は標準偏差にもｻﾝﾌﾟﾙｻｲｽﾞにも依存して変わるから、検出力が変化する。不慣れだとわかりにくい概念。良い説明法を探す。

ncp = (mu - mu0) / (sigma / sqrt(n)) なのだが、けっこうあっさりと書いてある。mu0は帰無仮説が正しいときの母集団平均。muは対立仮説の母集団平均である。sigmaは標本標準偏差、nはｻﾝﾌﾟﾙｻｲｽﾞである。

くどく考えると、通常の t 分布というものは検定論では、”帰無仮説が正しいとするとこの統計量は t 分布するはずですよ”という文脈で使われている。検定論の場合、差がないことを帰無仮説として設定するから通常の t 分布を用いる。つまり mu は mu0 を母平均にもつ母集団から得られていることを疑わない。たまたまずれているだけなのだ。mu0 = 0 のように仮説をおいていれば、得られた t 値がどれだけ稀なものかがわかるので検定に用いることができる。逆に区間推定の場合、mu0 に仮定をおかずこれを推定しているわけで、非心度を考えるような問題ではない。

非心 t 分布を考えるときはどんなときか。これは2群もしくは対応のあるﾃﾞｰﾀであってもいいのだが、母平均に真に差があるとしたときに、その分布を考えねばならないときである。こういう場合、分布の形状がゆがんでしまう。この分布は通常の t 分布として扱えない。

よって対立仮説としては非心 t 分布を考える。これはあくまでも”差の”ﾃﾞｰﾀの分布であるため、単純に t 分布を平行移動したものにはならないのである。このとき、非心度は2群の母平均の差が標準誤差の何倍にあたるかで定義する。あとは帰無仮説が正しいとしたときの棄却限界が、対立仮説の分布においてどのくらいの片側確率を持つか確かめればよい。

ここまで書いたが、やはり理解が浅い部分がある。次の機会にもう少しうめておこう。

検出力と例数設計

検出力と例数設計のための基本的講義資料の作成をする。これはなかなか難しい。例数設計には2通りのｱﾌﾟﾛｰﾁがあり、一つは信頼区間がある一定の自分たちが許す幅におさまるように設計する、精度ﾍﾞｰｽのｱﾌﾟﾛｰﾁである。もう一つが、検定で有意差が得られることを目標にしたものでこれが検出力ｱﾌﾟﾛｰﾁとなる。後者は検証試験で行なわれるやり方なので、PPK/PD解析者にとっては実は Critical ではない。多くの場合、Phase II であっても前者が重要視される。そういう意味では、検出力計算についてどこまで教えるかというのもあるのだが、やはり正統派でいくというのが目標である。

簡単なケースでは PROC POWER で一撃で計算自体はできる。これで答えが出るまでに課題を落とし込むことが大切で、複雑な検出力計算はそれをやってる自体OUTということもありえるのだが、今の目的は PROC POWER の答えをより素過程の計算からトレースできる人をちゃんと育成することだ。

PK の線形性解析

訳あっていろいろと考えた。やはり興味が深くなってきた。もっと勉強しないといけない。

Power Model

線形性解析にこだわっている。Power Model のふるまいに調べたいことが多い。用量規格化曝露量の分散分析よりも非線形への検出力は高いことはそうであろうと思う。しかし評価用量が少ないときのメリット・デメリットは何か？　まだ自分は何かをつかめていないらしい。

本

入手した。notation も含めて、かなり難解であるが読み応えもある。

ほしい本

ほしいのだがなかなか本屋にもいけず、入手できない。早く読みたいのだが。

LME

そろそろ復習しないといけない。NONMEM に最近ずいぶんと入れ込んでいるし。

赤本

Bayes 統計を R で学ぶ。Use R! という右上のヘッダーがおもしろい。Robert Gentleman が editorだからしょうがないのかも。

積命題と和命題

どれか一つでも有意（とか何とか、とにかく自分の主張したい仮説が正しいと言える）であれば、対立仮説が正しいとすることを和命題を設定するという。これとこれとこれが正しいときにはじめて対立仮説が正しいと主張することを積命題という。ORとANDであり、当然ANDが厳しい。

多重性の問題は、和命題において生じる。積命題については生じない。

この話は、”閉検定手順”というところにつながる。Marcusが1976年に提唱したものだが、定式化はその10年くらい後である。Type I Familywise Error の抑制であり、これによって仮説族が明確であれば多重比較を構成できるようになった。

容易にわかるように、厳しすぎる基準をおくことは簡単である。それは強い抑制といわれているが、自分たちが検出したいｼｸﾞﾅﾙ強度を変えすぎるのは良くない。ﾃﾞｻﾞｲﾝによって有意水準を変えるということは甘くするにせよ厳しくするにせよ、一貫性がないことになる。統計解析の場合、同じ土俵で議論することは大切なことである。

母集団薬物ﾃﾞｰﾀの解析（本）

復習する。Population Pharmacokinetics の和書としては数少ない一冊である。Bayes推定で1点の測定値から曲線を推定することへのｱﾗｰﾄは有名であるが、ちょっと読み返す。

対数正規分布

ｾﾐﾅｰ用に少し丁寧な教材を準備。対数正規分布は Web 上でわかりやすい記事が転がっているので、参考にしやすい。使う割にはﾃｷｽﾄの記述が薄いからなのかもしれない。Web では幾何平均や幾何標準偏差の検索例も多いようだ。

こつがある。対数正規分布の確率密度関数を丁寧に観察することだ。そしてｸﾞﾗﾌに描く。R を使えばこれは一発だ。そして、算術平均、標準偏差の関係を確認する。分位点を書き加えてもいい。母集団ﾓﾃﾞﾙでも重要な分布だから、よく時間を使うべきだろう。ただ、NONMEM の場合、相対誤差ﾓﾃﾞﾙとして表記しても結果が変わらない。Taylor展開してみればそれも理解できる。

動態の人が混乱しやすいのは、2つの要因が混ざるからだ。つまりﾃﾞｰﾀの分布が対数正規分布しているということと、独立変数に対しての関数関係が指数的で対数をとると直線回帰が適用できるということの2つが。両者は全く無関係と考えることができるのだが、ここは非常に混乱をもたらしやすい。

そして対数をとった直線回帰は、幾何平均の回帰線が推定されるという点だろう。これも理解しておかないと、間違いかねない。これは Power Model の理解には欠かせない。

動態の場合、様々な変換や複雑な計算を要するﾊﾟﾗﾒｰﾀがある割には、統計的なﾄﾚｰﾆﾝｸﾞが弱くなってしまう。自分も含めて常に謙虚であれかしと思う。

BEのコード

BE 解析の SAS Code。ありふれているためか(?)、Webでこの辺りを調べてもあまり見つからない。ﾓﾁﾍﾞｰｼｮﾝの問題か？

PROC MIXED data=XX;

CLASS DRUG PERIOD ID;

model LnAUC = DRUG PERIOD;

random ID(PERIOD);

lsmeans DRUG / CL alpha = 0.1 diff = control('1');

estimate 'Average BE' DRUG 1 -1;

RUN;

今日

佐和隆光先生の”回帰分析”を出張の帰りに購入。先日立ち読みしてから、気になっていたので。

ぱらぱら読んだだけでも、すごく丁寧な本であることがわかる。Econometrics の高名な碩学でいらっしゃることは知っていたが。回帰分析を実用する人の立場をよく配慮されて執筆されていると。

非線形

言わずと知れた名著。NONMEM も言及されているところがなんとなく好きだったりする。内容はきっちり難しいが、必読書だと思う。この本はまだ読みきれていない。

ほんの少し

ほんの少しだけ復習して就寝。

非線形まで到達してないなぁ。今月中に進めるだけ進もう。

今日は

少し気分転換して BE。すぐに出張が入るのでそのための解析も必要だが、そちらはほとんどやってある。資料作成にもう少しかかりそうだ。

しかし欧米の書籍のﾚﾍﾞﾙは高い。和書よりも2－3段高いところに、まだまだたくさんの出版がある。この差が生じているのは、何故なのだろう。

分位点の信頼区間

PROC UNIVARIATE data = ○○ CIPCTLNORMAL; で分位点の信頼区間をﾃﾞｰﾀが正規分布すると仮定したときに算出可能である。CIPCTLDEFでﾉﾝﾊﾟﾗﾒﾄﾘｯｸ(分布ﾌﾘｰ)の算出ができる。

少数例だと分位点は信頼性が低下するが、それはこの区間の幅が示していると考えてよいのだろうが、算出するという発想がなかった。もう少し研究が必要である。意外と血中濃度というのは、こういう考えが大切になってくる局面があるので。

混合効果モデルの書籍2冊

統計解析ﾊﾟｯｹｰｼﾞｿﾌﾄや特に R のﾗｲﾌﾞﾗﾘなどはそうだが、解説書を読むのはほとんど必須である。基本解析作業のようなものは、Helpや最近ではWeb上で調査することも可能であるが、やはり細かいところは書籍で読まないとわからない。

nlmeﾗｲﾌﾞﾗﾘは今までずっと使ってきてはいたが、この本を読むのと読まないのでは理解度が相当変わってしまう。今、半分読んだがずいぶんと使いやすくなってきた。groupedDataってこんなに便利なのかと感心したりして。ﾃﾞｰﾀにｸﾗｽﾀｰがあったり、階層構造があるときは、非常に使いやすく特に描画能力が高い。Excelでやってられない作業が一瞬である。

こっちはいかにもSAS本である。どちらかといえば、連続型変数に対する回帰というよりはｶﾃｺﾞﾘｶﾙな説明変数に重点を置いているかなというような気がする。乱塊法やSplit-plotなどの話から入っていく辺りがそれっぽい。

MEM in S については理論の章と実践の章がかなりはっきりわかれている。実践の章だけ読んでも分析ができるようになるような工夫がされているというわけだ。SAS for Mixed Modelsのほうは、もちろん実践主体であるがその前に SAS の本である。細かいｽﾃｰﾄﾒﾝﾄの話にどっぷり入っている。構成としては、自分がやりたいことを探すのが大変だったりするのだが、検出力計算などにも章が割かれているのは良いと思う。

MEM in S については前半部分は輪講できるくらいの資料の準備ができてきた。今日、線形混合効果モデルについてもう少しやったら、明日からは SAS のコース準備に入る。

lmer

lme4パッケージのlmer関数を初めて使ってみる。やはり売りは変量効果の記述が楽になって、より複雑なﾓﾃﾞﾙにも対応できるようになったということなのだろうか。random=○○と記述しなくとも、()で書いておけばいいというのは楽である。かつ複数組み込むときに、random=list() でくくらなくてもいいというのは、なかなか良い。

しかしaugPredとかintervalsとか便利関数がまだ使えないのが玉に瑕である。使わないかといったら、それは使うでしょうというような機能なので。

あと、同じ計算法を指定してもlme関数と微妙に結果が異なる。lmer関数のどの辺りに修正が入っているのだろう。しばらくはlmer関数を使う予定はないのだが、勉強は必要そうだ。

資料作成中

教育用資料の作成中。自分のこれまでを振り返ると、あるカリキュラムをこなしてきたということではなく無手勝流に必要なものを研究してきただけだから、カリキュラムというようなかっちりしたものを作りづらい。

PK解析、薬力学解析と統計学を融合させた領域にのぞめる人材を作るというのが来年の目標である。もちろん、自分自身も充分ではないのだが、その状態でも育成を始めないといけないのが教育というもの。

2008年は混合効果モデルについての理解をもっと深めたい。しかし線形混合効果モデルというものは、線形モデルや特に実験計画法の理解が深くないと理解しづらいものだ。後者は苦手なので、今も困っている。

BE

SAS コードも書かれており意外にわかりやすい本だった。内容についてはいろいろと述べておきたいところもあるが、BEや最近多くなっている2期以上のCrossover BE、他に用量比例性、有効性、QT、母集団解析についても章が設けられている。

各章の冒頭に Introduction として一種の小話（？）が書いてあるが、けっこう楽しく読めた。

等差3水準の回帰

10、30、50 のように等差で3つの説明変数をとった回帰分析をするとする。このとき、中央の30 のﾃﾞｰﾀは回帰分析の傾きに影響を与えない。つまりどんな数値であっても、回帰直線が平行移動するだけである。最小二乗法を使う限りはそうなる。

繰り返しがあってもそうなる。ただし、傾きの信頼区間は低水準、高水準の両方から得られる直線から中央の測定値がずれるほど拡大する。

混合効果モデル

今日もこの本。

混合効果モデルを R で

今日は休日だが仕事を3時間。その後で、Pinheiro の混合効果モデルの本に戻る。

Chapter 1 は線形混合効果モデルの基本。混合効果モデルへの動機付けとも言える。

Chapter 2 はlmeの理論部分。ここは読むのに骨が折れる。初読の際は飛ばしたほうが無難な部分。

Chapter 3 は、R の groupedData の作成法。lme にも nlme にも共通するから重要な部分だ。そして Chapter 4 は lme の細かい実行に関する部分だ。

この本はかなり丁寧だと思う。nlme ライブラリを用いるのに、これを読まないという手はないだろう。線形モデルを卒業するかしないかの人にも、混合効果モデルにエントリーさせるための書籍であるといえる。

もちろん、非線形混合効果モデルにも対応していて、そのために読んでいる。しかし、基本は常に線形モデルにある。

こっちもやっている。プラットフォームをまたがって、2足のわらじをはかないといけないのは、母集団薬物動態解析の宿命かもしれない。

PKの線形性解析

用量漸増で試験している場合は線形性解析に、Population Linearity よりは Individual Linearity をまず考えるべきだから、線形混合効果ﾓﾃﾞﾙを適用して解析する。意図せぬタイミングで非臨床TKに関して（！）今日そういう解析をせざるを得なくなった。今月の予想外ｲﾍﾞﾝﾄの発生率は異常に高い。

ｲﾍﾞﾝﾄが起こるたびに成長できるチャンスを得ているのだが、ちょっと食傷ぎみである。

PROC MIXED data = ○○;

class ID;

model LnCmax = LnDose / cl;

random INT / sub = ID;

RUN;

なお、今回の Data では、PROC GLM を用いて個体間の対応を考慮しない解析をしたときとほとんど差はなかった。

対比

また対比について考えている。分散分析・実験計画法の知識が弱いからだろうと、自分に見当をつけている。線形モデルの中で勉強することでもあるのだが。

F分布を用いて検定するとき、包括的帰無仮説の残差平均平方に対して対比平均平方を用いて検定する。これはいいんだっけ。たぶん勉強したことないかもしれない。さらに最大対比法もやり方は知っているのだが、なぜこれでいいんだっけ。

どうも CONTRAST は弱い。久々に考える機会ありだった。

多重比較

今日は統計講義だったが、帰り道にふと考えた。よく多重比較が生物統計では顔を出すのだが、受け取られ方がちょっと違うような気がする。今日のようなコースをレギュラーに聞いている統計に関心の深い研究者たちは特に問題がないのだが、一般の人たちについてだ。

多重比較は family-wise に危険率を制御できるということは、たいていの場合納得してもらえる。逆に制御しないと Type I error が大きくなることもである。しかし頭で納得してもらえても、依然として意識が「どれかに差があるだろう」というような差がある対比較をまさぐるような意識が強く、その探索が甘くならないような防波堤に多重比較を使おうとしているような気がする。

そもそも多重比較はそういうものではなく、解析計画自体は事前に立案されるものである。そしてそれは、比較する対象を事前に明確にしてErrorｺﾝﾄﾛｰﾙするためのものだ。多重比較は、Exploratoryな解析にはあわないものだと思う。Dataは結局のところある幅をもって分布するのだから、広く離れた値のもの同士をとれば差がある結果を得るのは当たり前だかｒだ。そういうことを、わざわざ統計解析で飾る意味はあまり感じない。DUNNETTは好例であろうが、最初から比較したい構造を決めているのだ。そしてそのために比較対の数を絞っている。各論に入る意味はないのだが、あまり後づけ解析にいそしむのもどうかと思う。

DUNNETTも誤解されているのではと感じる。対照群がある対比較にはDUNNETTと決まり文句のように使われている。教育の中でそういう言い方をするときはある。ただ、DUNNETTがむくのは対照薬とそれ以外の薬剤A、B、Cというような構造のDataである。よく誤用されているのだが用量反応性のDataにむいているわけではない。対照群以外の群内に構造が存在するときは、それにむいている解析を使うべきである。説明変数が量的あるいは順序的なときに、それを無視した手法を選択するというのは情報量を落としやすいと思うのだが。

なかなか難しいなと思う。私もそれほど自信があるわけではないが。

最近

教育ということに関しては大スランプといえる。大も大なりという感じで、本来なら教えることで自分も学んでいったはずなのだが、今は少し内向きになっているようだ。自分のためのことは、それなりに進歩しているのだが。あせっているようだ。

教育法としては、Contrast や SAS の Estimate ｽﾃｰﾄﾒﾝﾄについてはけっこううまい話し方ができるような自信がでてきた。この辺りは線形モデルへの理解度が深まるほど、使いようがでてくる部分である。もう少し粘着質に理解しようとすべきところだ。M & S という局面よりももっとオーソドックスな解析として汎用している部分だけに軽視できない。しかし教育面ではスランプだ。周囲からの必要性があまり感じられないからかもしれない。

医薬統計

かなり著名な本であるのでほとんどの人が読んだことがあると思うのだが、医薬統計Q&Aの改訂本が出版された。これはなかなかよい本だと思う。実務寄りの記述が多いのと、本の焦点がプログラミングでもなくかといってがちがちに数理であるわけでもなく、医薬業界で統計に携わる人にとっては優れたバランスであると思う。

NONMEM とか Levenberg-Marquardt法といった記述もあり、動態ﾓﾃﾞﾙについても臨床で用いるPKについては基本的なところが書かれている。今日買ったばかりだから細かい感想はまだ先になるのだが、前版のことも考えるとよい本であることは間違いないだろう。

最近は統計の重要性がかなり認知されてきていて、日本語で勉強する環境は加速度的によくなってきている。その分、不勉強の言い訳はどんどんきかなくなってきている。

今日

癖がある本という評判も多いが、Dalgaardのこの本は和訳もあるしなかなか使いやすい。

今日はこれを読む時間があった。通勤時間にこれを読みながら、英語のﾋｱﾘﾝｸﾞ練習をする。

会社では英文レポートの業務が今多い。

ちなみにこの本は、以外にﾛｼﾞｽﾃｨｯｸ回帰がよく書いてある。テキスト向きだと思う。

ﾛｼﾞｽﾃｨｯｸ回帰から最尤法、その先

PK/PD でﾛｼﾞｽﾃｨｯｸ回帰が利用できることはかなり知られてきている。正確にはE-Rモデルでの適用となる。まだ試したことはないが、有効／無効の2値判定が24時間の間に頻繁に判定されることは臨床試験ではほぼないといえるので、24時間という単位で薬物の血中濃度推移と有効率を対応つける必要性はほとんどない。むしろ、ある治療期間中の有効率が判明していて、各被験者についてその期間での平均曝露量を推定できるような、一日AUCやCmaxが測定されているかPPK/PDモデルから推定できるという状況が多いと思う。よって、代表的曝露量と有効率の関係を解析するから、ﾛｼﾞｽﾃｨｯｸ回帰が応用される。

細かいことを言えば被験者ごとのﾌﾟﾗｾﾎﾞ効果が異なることが想定されるし、場合によっては曝露に対する反応性さえ異なるかもしれない。両方考えられるが、臨床経験上はﾌﾟﾗｾﾎﾞ効果のほうが大きい。裏を返すと病態重篤度がﾍﾞｰｽﾗｲﾝ反応性を支配していると考えることもできて、これも薬剤非曝露状態での病態改善を左右するから同じことになる。クロスオーバー試験以外は同一被験者でﾌﾟﾗｾﾎﾞと処理の両方を体験することはないから、ﾍﾞｰｽﾗｲﾝの差か曝露に対する反応性の差かは切り分けることが難しいことは確かなのだが、こう考えることは無理は少ないのではないかと思われる。ただ、こういう個体差があるので、単純なﾛｼﾞｽﾃｨｯｸ回帰というよりは、そこに変量効果を含める必要がある。GLMMとなるわけだ。

しかし、変量効果まで教えるのはなかなか大変ではあるから、通常の一般化線形ﾓﾃﾞﾙの枠組みでまずは教えることになる。

ここで最尤法が登場する。正規分布でない誤差分布を考えるとすると、最小二乗法はとうとう出番がなくなってしまう。統計教育コースにおいて最尤法はとうとう出てきた、というタイプの概念であって、これが出てくると統計学がある意味自然に思えてくる。ﾊﾟﾗﾒｰﾀｰ推定はこうやりたいという姿というかありようを、素直に示してしまうわけだ。けっこうチャレンジングではあって、最尤法のありがたみは自分が知っている確率分布の種類の多さに依存するところがある。Gaussianではありがたみはわからず、BinomialやPoissonみたいな基礎分布にどれくらいすんなり入れるかが最尤法については大切である。こういう意味でﾛｼﾞｽﾃｨｯｸ回帰の重要性は高いのだ。

そしてこれは徐々にではあるが確実に、Bayes推定というもの後姿をかいま見させる入り口になってしまう。これは望むと望まざるとに関わらずである。もちろんやってることは違う。ただし、そろそろBayesがちらつき始めるタイミングなのだ。経験分布を考慮しないとだめだと感じた瞬間あたりから。

このあたりの教育は楽しいだろうと想像する。教える方だって絶対面白いからだ。ある意味、やっと統計解析で実用的な部分を俯瞰できるところまでくる場面だからである。私には楽しむ余裕がないのだが、ここは充実感を覚える部分だろう。

そこでBayesまでいって、モンテカルロの基礎がつまれたくらいで一度教育を振り返る局面に来ると思う。マルコフ連鎖まで休まずにいくのは難しいと思う。MCMCまでいくと、頭がごっちゃごっちゃになってくる学生は多くなるんではないだろうか。このあたりの実装となると、GibbsとかMetropolis-Hastingsの話になってしまうし。統計教育をどこかで区切るなら、きっとここの前あたりであろう。

「PPK-PD と M&Sのための統計」について考え続けてきていて、こんなところにいる。混合効果ﾓﾃﾞﾙはある意味当たり前だから死に物狂いでがんばればいいのだが、この辺りを「切り出して」勉強するのは狭いと思う。やはりその基盤や周辺に敬意を払って学んでこそ、いいことがあるんじゃないかと思う。

自分はちゃんとできているか？　いや全く無理だと思う。無理だけど努力するという世代にいるし、そこにとどまるのではないかと感じる。でも基盤を作っておけば、次、その次の世代は、もっと高く飛べるはずなのだ。教育とか伝承の意味はまさにそこにあるのだし。

線形モデルの先

この本を読んでいる。Extending the linear model with R という Faraway の有名な本なのだが、Extend する以上、一般線形ﾓﾃﾞﾙはほとんどでてこなず、最初から binomial dataでありそこから混合効果や一般化線形ﾓﾃﾞﾙやTree analysis に進んでいく。NonparametricやNeural Networkも含まれているので、要は”Linear Modelの次に読む本”ということだと思う。

M & S というのはやるだけならそれほど高度な統計学が必要であるとは思わないのだが、解析対象であるﾃﾞｰﾀの種類や方法論が多いほうが融通が利くことも確かである。一般化線形ﾓﾃﾞﾙなんてもう必修科目みたいになっている。

M & S というと少し熱くなるところもあるから、とりあえずクールダウンして統計リテラシーとしてやってるくらいのテンションではいるつもりだが。

一般化線形混合効果ﾓﾃﾞﾙ

R2.6.0をｲﾝｽﾄｰﾙしたことがきっかけで、今までやってなかったことをやろうとしてみている。とりあえず、glmmMLを試してみる。

glmmMLは通常のﾛｼﾞｽﾃｨｯｸ回帰と異なって、変量効果を一つ入れることができる。薬力学で考えられる場合というのは、各被験者から複数回評価されている場合であって、各被験者ごとに有効率が異なるような状況である。これには、Baselineが異なっていて反応性は同一であるﾗﾝﾀﾞﾑ切片ﾓﾃﾞﾙと反応性にも個体差があるﾓﾃﾞﾙが考えられるが、glmmMLの場合はﾗﾝﾀﾞﾑ切片ﾓﾃﾞﾙに対応している。

使い方はけっこう楽な関数だと思う。解析の自由度が決して高い関数ではないからだろう。

欲しい本

これ良さそうだなぁ。

今日、本屋で見かけたけど、買わなかった。積読は怠惰であるが、今までの人生振り返ってそういうことは多かった。でも、いざというとき積読が効くこともあったので。

とりあえず、心のウィッシュリストに入れておく。

本を購入

今日、この本を購入。教育用であるが、SVMやSOMの章もあったりしてなかなか面白い。

統計ソフトの普及と大量のデータ取得・管理の世界から、データマイニングがかなり流行している。今日昨日始まったことではないが、おそらく今、統計の初学者が加速度的に増大しているのは、データマイニングがらみであろう。

R の書籍を執筆されている先生方も必ずしもマイニングを志向しているというわけではないと思うが、ニーズに答えるということではマイニングを前面に出していくことになっているのではと想像する。

あと2冊、R の和書では欲しい本があるが、今日はこれで。

質問を受けるのと軸足と

多重比較について質問を受けたり、PK-PD解析の質問を受けたり、戦略を話したりと。。。

ただ一般線形ﾓﾃﾞﾙや線形混合効果ﾓﾃﾞﾙとか言われると少し安心する。いや、安心するほど理解してはいないのだけど。NONMEM, PK-PD, SAS, R とか言われると、やはり拙くてもそこが軸足だと思う。

数学をもっと鍛えるべきだとは思う。ベクトルや行列表記には大分慣れてきてはいるのだけど、まだまだどっぷり感が足りない。もっと埋まらないと。SAS ﾌﾟﾛｸﾞﾗﾐﾝｸﾞもだいぶ好きになってきた。一方、Rで先駆的な解析のパッケージを追うことはなかなかできてないし、WinBUGSもやらないといけないのだが。。。

この軸足、ある程度周囲はわかってくれているようで、今はだいぶ助けられている。

幾何学的なｲﾒｰｼﾞ

勉強していて思ったのだが、統計における幾何学的なイメージは理解を容易にするのだが、応用、計算にはあまり適用できない。

計算をBlackboxにするのも良くないから、やはり実Dataで行列計算をトレースさせることになる。幾何はイメージ構成のガイドということで。

線形代数

線形代数を少しかじったくらいでも、PCAがやっていることがすごく感動できると思う。線形代数というのもだいぶ苦労させられている分野だ。その理由は、統計解析にとってはなくてはならないが、線形代数の初学者が学ぶことと、統計解析を行列表現で統一的にあらわすことは決してイコールではないからだ。

線形代数を身にしみて理解しようとすると、統計も含められるとは思うが特に工学的な立場だと計算してなんぼであるということがある。数字の入った行列を扱うことが大切で、その中の数値に物理的意味があるわけだ。一方線形代数を数学側から理解しようとするとそういう位置づけはなくなり、抽象的な扱いが重要になる。線形代数は、数学的な”構造”を理解する第一歩になるからである。もっとも扱いやすい空間を扱うのだ。そこを具体的な演算もやらせて行列を身近に感じさせようとすると、手計算で済むような一桁の数値が入った行列の計算問題にまで下りていってしまうし、それはそれで仕方がない。

ただこれも、R や SAS が助けてはいる。行列演算のﾚﾍﾞﾙで統計解析を理解させてくれる手段があるからだ。

最近、線形代数がさびついてきた気がする。もともと研ぎ澄まされてはいないのだが、不安になってきた。

難しい本を読んでかっこつけたいという気分もあるが、この本をしっかり理解するのももちろん大切だ。思い出すときはここからやるようにしている。本というものは、背伸びして得るものもある。でも自分に正直になることも必要である。この本は読んでいて本当におもしろい。

t 検定と分散分析

2標本の t 検定と2群しかないDataでの一元配置分散分析のP値が一致することは有名な話で、それは F 統計量が t 統計量の2乗にあたるからである。

計算機統計の基本として、P値、例数設計と検出力、信頼区間の算出が基本ﾃｸﾆｯｸであるが、 t 検定の P 値を R で1000回シミュレーションするコードと、SAS で 1000 回シミュレーションするコードをつくったのだが、SAS での結果がおかしい。SAS では PROC GLM を用いて PROC TTEST を用いなかったのだが、まさかこれが？と一瞬疑ってしまった。

調べてみたら、Data生成のときの単純なコーディングミスで一安心した。output を 2回書いていたので、群内でつねに重複するDataが発生していた。こうなると群内平均平方が小さくなってしまうので有意差がつきやすくなってしまう。1000回のシミュレーションで、150回くらい有意になってしまっていた。

計算機シミュレーションでは小さいミスでも結論が異なってしまうことを、いまさらながら実感する。

R のコードはこれ。

Sim.ttest <- function(rep) {
alpha <- 0
for (i in 1:rep) {
   A <- rnorm(6, mean=1.5, sd = 0.15)
   B <- rnorm(6, mean=1.5, sd = 0.15)
   D <- data.frame(x = A, y = B)
   out <- t.test(x = D$x, y=D$y, var.equal=TRUE)
   if (out$p.value < 0.05) {
   alpha <- alpha + 1
  }
}
return(alpha)
}

SAS のコードはこれ。

Data d1;
array mu(2) mu1-mu2;
mu1 = 1.5; mu2 = 1.5;
do i = 1 to 1000;
  do group = 1 to 2;
  do j = 1 to 6;
   if group = 1 then
    y = RAND('NORMAL', 0, 0.15) + mu1;
  if group = 2 then
   y = RAND('NORMAL', 0, 0.15) + mu2;
   output;
  end;
end;
end;
RUN;

ods NORESULTS;
ods listing close;
ods output OverallANOVA = out;
PROC GLM data=d1;
by i;
class group;
model y = group;
RUN;

DATA result;
set Out;
if (ProbF < 0.05) then flag = 1;
else                 flag  =0;
if (ProbF > 0) then output;
RUN;

ods listing;
PROC FREQ data=result;
tables flag / norow nopercent;
RUN;
QUIT;

今日の統計

今日はこれに取り組んでいる。読んでいて、基礎的なところで誤解していたことに気づいたり、χ2分布がΓ分布の特別な場合であることを思い出さされたりして、ちょっとびびる。日ごろメインテナンスしないと、いろいろと基礎がこぼれていってしまうと感じる。

久しぶりに統計ソフトから離れた一日を送っているが、これはこれでよいものだと思う。やはり計算機と理論のバランスをその人なりにとっていくのがいいんだろう。私の場合はどちらも弱いのだが、その人なりでよいと思う。

夜はこれで。SAS マクロもだいぶ多用するようになってきたが、基礎が弱い。こちらは業務の中では重要度高だから、PPT資料とかも作らないといけない。

脱落の統計

COPDとか喘息でよくあることとして、”急性増悪”がある。あるとき、病態のNatural Courseとは思えぬくらいに一過的に病態が悪化する現象である。多くは可逆的でありまた回復することになる。要因は疲労だったり別のアレルゲンに誘引されたり様々ある。糖尿病のシックディのようなものだが、患者のつらさは呼吸器疾患のほうがつらいだろう。

これが臨床試験中に発生すると、多くは脱落を引き起こす。こういうときに、LOCFを用いるとバイアスがかかる。経時的な観察をしていると、有効性を小さく見積もってしまう。呼吸器疾患の場合、たいていは経時的測定とみなして解析することもあるから、これはあまりよろしくない。では有効性解析対象から外すかといっても、情報量の減少を招くことになる。

当然だが”急性増悪”は患者においてのみ発生する現象である。そしておそらくそれはランダムではない。FEVが悪化している患者であればあるほど、急性増悪をおこす可能性は高い。FEVの絶対値が低い患者でもそうだろうし、悪化速度が速い患者でもそうだろう。MCARとは考えにくいわけだ。

こういうとき、Dropの統計ﾓﾃﾞﾙが必要になる。これを M&S に組み込むのが重要である。脱落のロジットに線形ﾓﾃﾞﾙを仮定することから始まるのだが・・・　これは経験が浅い企業ではなかなか使えまい。こういうところのﾓﾃﾞﾙ化というのは、その治験薬に限らず企業の当該疾患領域の開発経験がダイレクトにでてしまうからだ。

困ったことに糖尿病でもまたあるのだ。効きすぎたりして低血糖までいってしまうと原則投薬中止されてしまい脱落することが多い。こういう被験者を抜いた解析は、有効性についてうまく説明しているとはいえないだろう。SU剤では良くあることではあって、血糖ｺﾝﾄﾛｰﾙを新規治験薬でもちゃんとしてくれるというのは治験実施施設に高い力量を必要とするだろうし、薬物動態に変動要因が多いとそれだけ難しい。さらに長期試験を行なうと、薬理作用に基づく初期脱落と、多くはないが二次無効のようにｺﾝﾄﾛｰﾙを逸する被験者も現れてきて、脱落のメカニズムに複数の異なる要因が寄与する。

脱落を統計的に考え始めて私はまだ日が浅いのだが、疾患領域や臨床からの知識でもかなりいろんなパターンを考えることができる。これをそれなりの仮定をおきつつ定量化するというのは、勇気がいることでもある。

環境

これはｶﾃｺﾞﾘｰ的には「統計」にしてよいだろう。仕事で韓国のことを調べる機会があったのだが、あの国は70年代後半から統計学教育を整備しはじめていて、現在は20学科以上を国内に抱えているようだ。もう一息で日本の私立薬学部の数を越えるのではないかと思う。

なぜ日本の大学には統計学科が皆無なのだろう。実際に業務で必要とする人は多数いるのだが企業内や民間のｾﾐﾅｰへの依存度をあげざるを得ないし、教育者も不足している。

最近実感したのだが、大学という組織が多数の専門家を世に送り出すことは大切なのだ。それもとてつもなく重要な役割があるのだ。彼らのほんの一部が大学の研究者になり、残りの大多数は社会で、ある意味普通の人として暮らすことになる。もちろん職能人になるものもいる。大学のステータスは一流研究者の輩出だろうが、それ以上に社会で活躍する常識的な知識レベルの人間を多数輩出する意味合いはあなどれないということだった。

有機化学の重要性が論をまたず、分子生物学の知識も当たり前であって、これら知識によって議論をすることはありふれているが、統計的な議論はそもそもそれをしようとする者がいないようだ。しようとすると「よそへ行け」という風潮が強い。私は統計の専門家ではないが、私でさえ正直そういうプレッシャーを受けることがある。何でこんなことになっているのか、私はわからない。

今日はひとりごとが多かった。SAS で NONMEM を実行して集計するプログラムについても進捗したのだが、少し気分が沈んでいるのかもしれない。

基礎本

いろいろ統計の本を使ってきたが、基礎本はこれが優れていると思う。丁寧な本なのだが、決してゆっくりとしたﾃﾝﾎﾟではなく、重要な概念はかなり早いスピードで触れられる。先を急ぐのがちょっともったいなくなるような本である。分厚いわけではないのだが、薄めの本でこういうスピード感を出しながら大切なことを落っことさないという本はあまり多くない。ある意味、細かい定式化をうまくはしょることに成功している。

もちろん、これは言わずと知れた名著である。これを中心に据えて勉強することに何の異論もない。このシリーズも「統計でData解析したい人」のために照準が合わされている。丁寧さがより増すので、きちんと勉強したい人むけである。確率分布の部分なんかはとても読みやすくできていると思う。ものすごく売れている本だが、著者のWebサイトもお勧めである。

でもこちらも読むべきである。というか、このシリーズはここまで読まないと薬学の統計はできないだろう。ここで一般線形ﾓﾃﾞﾙや最尤法、ﾛｼﾞｽﾃｨｯｸ解析やマルコフ連鎖の話がでてくる。アラカルト的な部分もあるが、他の本を考えるとこういう内容の本は国内の統計学習者には重要である。数理的でありすぎないという意味でもすばらしい。このﾚﾍﾞﾙの本を充足することが、国内教育環境には大切ではないだろうか。ニーズが多いﾚﾍﾞﾙだと思う。

この本も評判が高くて事実売れているようだ。装丁からはドリルみたいな印象を受けるのだが、重要なことが良く書いてあって、特に統計解析が実技としての面を持つことが強調されている。M & S とかを知る前に読んでいた本だが、導入する力は高いのではないかと思う。藤澤先生の本はすばらしいが数理的面が強い。この本は実技としての面が強い。

この本も良い。この本はアラカルトであり、これで体系的に学ぶための構成にはなっていない（学べるのだが）。他の本に書いてないことが書いてある。ある程度学習したらこの本を読むと良い。

ちなみにこれは昔からあって買うかどうするか迷ったのだが、最近買った。まだあまり読んでない。

これは定番の名著だ。原著へチャレンジした方が良いのだが、今日は和書ということで日本語版をのせておく。もちろん趣はかなり異なる。これは疫学の本である。ｶﾃｺﾞﾘｰデータ解析について学ぶにも良いと感じさせられることが多い本だ。もちろんその辺りをきっちり学ぶには、とても良い和書がある。

基礎本についてあげてみた。中級本については別の機会にする。中級本は解析ソフトの使い方を意識している本の重要性が高く、必然的に洋書が多くなってしまう。最近では、SAS と R の双方をタイトルにあげた本が増えている。なお、これらの本の書評はこの記事では不充分であって、今日は顔見せに過ぎない。後日もう少し丁寧に紹介しようかと思う。

薬物動態の書評も載せるかもしれない。あるところのカスタマーレビューに別々のニックネームでだいぶ書いてしまっていて、あれはあれで良いかと思う。ブログならもう少し細かく書けるから、2007年の立場から再度書きなおしてみようかと思う。

統計ばかり？

M & S へ勧誘するときに最大の障壁となるのは、統計の重要性の高さである。

M & Sで臨床薬力学を考えるときに統計がどのくらい必要かと考える。NONMEM でのﾓﾃﾞﾘﾝｸﾞからシミュレーションにいたるまで、統計の深い知識がなくとも実行するだけならやれる環境になってきたのが昨今の状況だと思う。

しかし「統計の深い知識」が必要でなくても「統計の基本的な知識」をﾌﾙに使うということは M & S でも重要だと思う。そしてそれは、実はかなりの統計経験への熟達がいるということを覚悟しないといけないと思う。PK-PD ﾓﾃﾞﾘﾝｸﾞ は臨床統計家が持っているバックグラウンドの中で整理して扱われてこなかったから、結局はちょっと毛色の違うメンバーで実行せざるをえない。そして臨床試験においては統計的知識はエッセンシャルである。

だから統計ばかりやってるという批判は的をはずれている。統計以外のことも大量にやっているのだ、ただ、統計知識をおろそかにしていると地に足がつかなくなってくる。少なくとも統計の非プロパーであるバックグラウンドを持つ場合は、やりすぎくらいでちょうどいい。

ちなみに、統計教育は日本の場合はそれほど充実していない。ﾚﾍﾞﾙの高い講座が少なく、誰でもが知っている。だからそれ以外の出身だったり、普通の薬学や理学の出身者は自分を専門家であるとは思っていない。

よって質問するときに「自分は統計のずぶの素人」と名乗ったり断りを入れることが多いようだ。Webの掲示板でも決まり文句である。しかし質問自体は「○○効果ﾓﾃﾞﾙで△△を計算し、□□を調整したところ・・・」とか、やりたいことはかなり込み入っている。正直、アマチュアのやる解析ではない。

仕事として解析をやっている以上、「自分は素人です」という断りを入れることを禁じた方がよい。私も素人だが、そこに甘えたくはない。「そんなこともわからないとは低レベルだ」と思われたくないからこその予防線を張る意味は何なのか？　そういわれたらプロとして傷つく方が結果として成長するには良いのではないかと思う。

夜になるとこんなことも考えてしまう。まぁ、下手にへりくだっても逆効果だと私は思う。混合効果を業務で扱うのは決してアマチュアではないから、立場から逃げてはいけないんじゃないかと思う。日本でこの分野をとりまく環境は特殊であるが、それはある意味みんな知っていることである。

輪講初回

演者としてまずANOVAの基本を話す。SAS の出力で１つだけわからないものがあって、そこをつっこまれて見事に「わかりません」と言ったが、こういうのは楽しい。どうしても標準誤差の出力で意味がわからないものがある。もちろん、LSMEANSｽﾃｰﾄﾒﾝﾄというのは行列計算した上での出力であり、表計算や出力どうしの加減乗除から簡単にトレースできない出力があるのは当然なのだが、意味はわかっていてもどういうことなのかわからないことがある。

PROC GLM の出力は全てを常に使うわけではない。ただなかなか味があるというか、意味は理解しておくべき出力が多いから、やはり勉強量は多いという実感を受ける。一般線形ﾓﾃﾞﾙは統計学の基礎だから、どうしてもそうなるんだろう。t 値などは使わないことが多いとは思うが、改めて説明したり頭に入れることは大切である。

総じて好評だったと思う。話を少しふくらませてしまいすぎたかもしれない。同時信頼区間とかは今日のイントロ的なChapterからはふれずとも良い話なのだが、テクノロジーの進歩がこういう計算も簡単にできるようにしてしまっているから結局話すことになった。ANOVAは検定のニュアンスが強く教育されることが多いが、それだけではないとふれるためである。

初級の R コースがあと2～3回で終了できるんじゃないかと思う。そっちは次に、R ﾌﾟﾛｸﾞﾗﾐﾝｸﾞとNONMEMによるSimulationコースに入る。

同時信頼区間

回帰分析にも入っていかないといけないのだが、この辺りをどう学ぶかを考えている。R でシミュレーションするのが一番よいと思うのだが、SAS でもできないだろうか。できれば常に両方並べていたいのだが。

概念を直感的にわかればよいと思うから、数理的につっこむつもりはあまりない。単一の母数の信頼区間の構成は楽であるが、それを組み合わせると長方形や立方体のような領域ができあがる。そしてそれは狭くなっている。信頼領域を頭になんとなく描ければよく、それぞれの単一の母数については保守的であるということを実感できれば良いのだと思う。

しかし輪講準備は思ったよりも大変である。やりやすい本を選んだつもりだったのだが、結果的に濃厚に準備せざるをえない。

もちろん、とてもいいことだ。

検定の方向

改めて”検定の方向”を考えようとすると結構難しい。分散分析は両側検定であるという話にふれざるを得ない機会があるのだが、これもわかりにくい話かもしれない。

わかりにくいのは、対立仮説の取り方だろう。帰無仮説の場合、各群の母平均が全て等しいという仮説になりこれはわかりやすい。独立標本 t 検定のときと変わらないからだ。一方で2群比較で話がおさまる独立標本 t 検定では、帰無仮説は差の母平均＝0が仮説となるが、対立仮説については差の母平均が正であるか負であるか母平均が異なるかという仮説の立て方がある。

一方で分散分析の多群比較の場合は、C>A>Bであったり、B>C>Aであったり、A>B>Cであったりと群間の大小関係は多数存在する。その関係をF値は”区別することができない”。平均平方の比が1より大きいことは、Dataが異なることでのばらつきかたよりも所属する群が異なることによるばらつきの方が大きいことを意味しているが、これは群間の関係とは直接関係がないのだ。だからある群とある群の大小を固定した状態で比較する対立仮説は置けない。よって意味的に両側検定になる。

誤解されやすいのはこういうことを考えて検定するときに、F分布の使い方が片側確率を計算するために利用されることであることだろう。だから分散分析は片側検定であるように思われやすいのだろう。分布の利用の仕方と、対立仮説のおける置き方には違いがあるのだ。

検定法の場合、通常ことわりがなければ両側である。片側検定が可能な場合はたいていそれと示されている。そして両側、片側は対立仮説の置き方を示しているのだと考えればよい。難しいことは多々あるが、こういう理解からまず学ぶことになる。

輪講準備

行列表記を使っていないが、確かに優れてわかりやすい本というものはある。こういうのはいい本なのだろう。ただやはり、R で行列計算を追って欲しい。IMLも学びたいのだが、こちらは時間が足りない。

他人に講義すると話しながらたくさんのことが自分でも整理されてくる。経験が浅いから、もっともっと講義も必要だ。何気ない質問に、いつも学ばされる。

逆推定続き

思い出した。2003年くらいの頃に発表が多かったようだ。佐久間昭先生の「生物検定法」の内容からヒントを得たものだったと思う。高橋ｾﾐﾅｰでも検討されていたはずじゃなかっただろうか。

答えを得るのに

Simulation することが多い。今日も R と SAS でやっていた。ただし、PK ﾊﾟﾗﾒｰﾀｰ同士に共分散をおいたときの Simulation が難しい。R では非常にやりやすいが、SAS の場合が大変だ。IML でやれるほどの技量はまだない。

実行した Simulation を社内教育資料の基礎につかうことにし、試しに話してみた。ことのほか好評だったので、もう少し作りこむことにした。

Davidian & Giltinan

結構有名な本で読みたいのだが、なかなか時間がない。というか、まだ買っていない。

プログラミングばかりしていないで、理論を学ぶことも大切なのだが。どうやるのかばかりを鍛えるのではなく、何をやるのかとか、あの人は何を言っているのかわかるということが最終的には重要になってくる。

明日から

Ette の Pharmacometrics の本のアタックを再開。SASﾏｸﾛﾌﾟﾛｸﾞﾗﾐﾝｸﾞについてはまだ本が届かないから、おあずけ。線形混合効果ﾓﾃﾞﾙとBEについては、がんばらないと。

2値Dataに思うこと

計数 Data についてあまり手を出したくない手法があって、実際試したことはない。Arcsin関数で変数変換するやり方だ。

変数変換というものはあまり好きじゃない。好きじゃない人が多いと思う。統計を専門とする人たちも、変数変換を推奨する人たちにはあまり会ったことがない。

私にとって最大の理由は非常にわかりにくいからである。要は、変数変換して○○分布になるということはその前のDataは、かなり特殊な△△分布をしているということになることが多い。○○として狙われるのは当然正規分布だからおのずと無理が出て、△△分布というものは数学的になんというか特殊で想像しがたい分布になる。

それなら元のDataの分布をもっと直感的にわかりやすい分布だと考えて、最尤推定した方がいいのではないかと思う。逆正弦変換や平方根変換というのはなんとも難しい。もっとも、そういう手法を使ったりしていると、「勉強しているなー」と思われるのかもしれない。

というわけで今日も Logistic 回帰という感じである。もちろん、わかりやすいわかりにくいは主観の問題かもしれない。尤度やオッズ比やハザードなんてものだって、わかりやすいものじゃないかもしれない。単に私が見慣れてきたということに過ぎないのかもしれない。

一般線形ﾓﾃﾞﾙ

Amazonからこの本を買った。邦訳が出版されていたので、日本語で読む誘惑に負ける。

Oxford大のWebから扱われているData Setを、SAS、SPSS、Minitab、ExcelでそれぞれDownloadでき、演習問題やそれに対する解析 Code も得られるようだ。

こういう本は欧米の場合、学部教育でよく使われている。Undergraduate のための書籍である。こんな本で勉強できるという学生は恵まれていると思う。私の学生時代の場合、そもそも PC がまだほとんど普及していなかった。「大学の計算機センターに通わねばならなかった」と書くともっともらしいのだが、残念なことに私はコンピューターに触れずに大学3年生までを過ごしている。4年生で初めて Mac を触ったことが PC 事始めである。

だからというわけではないのだが、一般線形ﾓﾃﾞﾙが分散分析や回帰分析を包含する一般化であることを学ぶまでには存外の時間がかかった。統計手法の場合、その使い分けが重要であることは間違いない。よって手法の紹介を行なうとき、Data のタイプと手法の対応関係や決定樹を頭に入れることに時間が用いられる。もちろんこれは大切なことであって、いや、極めて大切といってもよい。

一方で異なって見える手法を統一化する概念を学ぶことは、特に企業現場で統計を実行する場合にはチャンスが少ない。時間が取れないと言ったほうが正確かもしれない。生物系の人の場合は特に普遍原理への傾倒は少ないのではないかと思う。

一般線形ﾓﾃﾞﾙの枠組みは、思考のメモリを節約してくれる。陳腐な言い方をすれば、「わかった」と思わせてくれる。私がそう感じた最初の本がどれだったか、記憶が定かではない。読んでみるとたいていの本に書いてあるからで、おそらく名著と呼ばれるものではない本だと思う。ただ、そういう目線を持って読まないと意外にこぼれおちてしまう内容でもある。最初に学んだときはかなりうれしかった覚えがある。そのときの Power Point 資料が残っている。

今度はこの本を使って、社内セミナーの一つもできればと考えている。計算機統計的な側面は少ないので、SAS セミナーにしようかと思う。

2値Dataの PD model

今日は1日、2値変数を解析対象としたﾓﾃﾞﾘﾝｸﾞを行なっていたが、かなり疲れた。今ひとつ残差診断やﾓﾃﾞﾙのあてはまり具合を診断して共変量の有無を診断することに難しさがある。解析の効率化のための小さいｺｰﾄﾞもいくつか書かないといけない。やった解析については PPT に経緯をまとめておいた。

夜は小さいﾌﾟﾛｸﾞﾗﾐﾝｸﾞをいくつかやっておこうと思う。頭を切りかえて、通常のﾛｼﾞｽﾃｨｯｸ回帰の解析も少しやってみるつもりである。

一般化線形混合効果モデル

ロジスティック回帰を試していてふと、被験者のようなクラスターが存在する場合に通常の一般化線形モデルが適用できないことに気づいた。

一般化線形混合効果モデル(GLMM) というものが存在しているが、まさかこれを研究しないといけない状況になるとは思っていなかった。PK-PDは多くが連続量を対象とした微分方程式モデルだからなのだが、現在は離散量のモデリングが大きな関心を集めており、ロジスティック回帰や生存時間分析の重要度が高いからだ。そこで個体間変動と個体内変動のように複数の誤差構造をモデリングできるGLMMに注目することになった。

GLIMM の場合は R でやってみるのが一番楽で情報量が多い。S-PLUS では ver 7 から対応しているようだが、現状では S-PLUS を購入するには研究費が少なすぎる。前にいた部署では何とか買えそうになったのだが、結果的にはどうでもいい機器の予算にまわされてしまった。

GLMM は NONMEM で実行することは可能だから、$PRED でもやっている。ただこういう解析は、統計解析パッケージで行なうほうが適しているのではないだろうかとも思うので、R をとっかかりとして研究していこうと思う。

ﾃｷｽﾄ

SAS for Mixed Models, Second Edition、Introductory Statistics with R、Pharmacometrics　の3冊をアタックしている。レベルはまちまちなのだが、欲しい内容をきれいに網羅できてはいると思う。1冊目や3冊目は「アタック」という言葉がふさわしい。

意外であるが、何であれ基礎の本というものはよく読み返すものだ。例えば線形代数の基礎的な本とか、ホーエルの数理統計学とか。本は読み返してこそのものなんだろうとは、最近よく思う。

けっこう大変だが楽しい。自分が学生だったらなと、最近思うことが多い。一番足りないのは時間だからだ。1ヶ月フリータイムがあったらどれだけ進歩できるか。

これも「雑」でよかったかもしれないが、「統計」にしておく。考えていることについての話題が多くなったのは、今日は少し疲れたからだろう。コードをのせることが目的のブログなのだが、日記や随筆になっている。Tips はけっこうあるのだけど。

PKBugs

MCMCの場合、いつモデル推定が収束したかということについての考えが最小二乗法や最尤法と異なる。診断については厳密ではなく感覚的なものも含まれる。

PK解析をしてみた。最初の1000回の反復は初期値の影響を強く受ける burn-in phaseととらえる。

これはクリアランスの推定履歴をプロットしたものだが、最初の250回を越えた辺りから真値に向かっていることがわかる（上段）。下段はMetropolis-acceptance rateである。

burn-in phase を終えたあとで10000回のサンプリングを実施したところである。もはや反復によっても推定値はほとんど変わらない。

原理的に言って収束性の良さはNONMEM よりは格段に優れている。アルゴリズムの安定性である。一回一回がある意味シミュレーションであるから、計算速度がほとんど低下しないとイメージしておけばよいのかもしれない。WinBUGsの画面スナップショットを示すが、これがフリーソフトであるということは驚きである。

確かに定番となるだけの使い勝手と性能であると思う。もう少しDoodleを使えればかなり解析に自由度を出せるだろう。しかし階層ベイズを強く意識させられるソフトであり、勉強はけっこう必要だ。微分方程式モデルまでいけるかどうかは、まだわからない。Pharmacoをまだ動かせない。

WinBUGs (3)

R2WinBUGs　を用いてRからWinBUGsをリモートで実行する方法がわかった。ExampleのRatデータの例をあげておく。

model.txtとして、Exampleからmodelの部分をコピーしたテキストファイルを、作業フォルダに作ってから以下を実行する。

# Rプログラム
x <- c(8.0, 15.0, 22.0, 29.0, 36.0)
xbar <- 22
N <- 30
T <- 5
Y <- matrix(c(
151, 199, 246, 283, 320,  145, 199, 249, 293, 354,  147, 214, 263, 312, 328,
155, 200, 237, 272, 297,  135, 188, 230, 280, 323,  159, 210, 252, 298, 331,
141, 189, 231, 275, 305,  159, 201, 248, 297, 338,  177, 236, 285, 350, 376,
134, 182, 220, 260, 296,  160, 208, 261, 313, 352,  143, 188, 220, 273, 314,
154, 200, 244, 289, 325,  171, 221, 270, 326, 358,  163, 216, 242, 281, 312,
160, 207, 248, 288, 324,  142, 187, 234, 280, 316,  156, 203, 243, 283, 317,
157, 212, 259, 307, 336,  152, 203, 246, 286, 321,  154, 205, 253, 298, 334,
139, 190, 225, 267, 302,  146, 191, 229, 272, 302,  157, 211, 250, 285, 323,
132, 185, 237, 286, 331,  160, 207, 257, 303, 345,  169, 216, 261, 295, 333,
157, 205, 248, 289, 316,  137, 180, 219, 258, 291,  153, 200, 244, 286, 324),
nr = 30, nc = 5, byrow = T)
matplot(x, t(Y), type = "l", col = "black", xlab = "day", ylab = "weight")

library(R2WinBUGs)

dat <- list("x", "xbar", "N", "T", "Y")
inits <- list(alpha = rep(250, 30), beta  = rep(6, 30),
+     alpha.c = 150, beta.c = 10, tau.c = 1,
+     sigma.alpha = 1, sigma.beta = 1)
params <- list("alpha0", "beta.c", "sigma")
setwd("モデル.txtを作成した作業フォルダのパスを入れておく")
inits <- list(list(alpha = rep(250, 30), beta  = rep(6, 30),
+     alpha.c = 150, beta.c = 10, tau.c = 1,
+     sigma.alpha = 1, sigma.beta = 1))
rats.sim <- bugs(dat, inits, params, "rats.txt", n.chains = 1,
+     n.thin = 1, n.iter = 2000, n.burnin = 1000, digits = 5)

WinBUGs (2)

PKBugs のHELP機能が使えるようになった。というかPKBugsはWinBUGs ver 1.3にあわせてつくられており、WinBUGs 1.3をインストールしてPKBugsを組み込んだらHELPがちゃんと読めるようになった。Webページをよく読むと、判断することができた。どうやらPKBugsをメニューから順々に実行することで、PPK modelingができるようだ。

Example Data でようやく実行することができた。NONMEM 慣れしているとまったく異なる使い勝手に驚かされるが、おそらくこれでよい。PKBugsはかなりWinBUGsから独立して使っている感じだ。収束判定に明確な基準がなく、ﾊﾟﾗﾒｰﾀｰ推定の履歴をPlotで確認して、もう変化しないと判断したところでやめるという感覚である。heulisticだが、まずはこれでよいのかと思う。

モデル診断など含めて考えると、まだできないことが多いソフトのように思える。プロット機能の多くはPKBugsからさらにUpdate予定であるようだ。かなり探索解析的なソフトであるが、計算速度が非常に速い。Pharmaco、WBDiffでの微分方程式モデルはまだ試していないが、技術レベルを上げていく価値はあるだろう。

となると、とりあえずモデルや初期値を与えるコードを自分で書けるようになる必要がある。やることが多いが、これはまぁしょうがないだろう。

WinBUGs

WinBUGsに取り組んでいる。基本的な使い方はわかってきたものの、いくつかのことが致命的にわからない。

通常の統計解析の部分はそれなりにHELPとマニュアルから理解できてくるのだが、このモデル記述は大変である。Doodleで描くこともできるしそちらのほうがどうやらお勧めのようなのだが。

PKBUGsについてはまだ動かし方がわからない。NONMEM と同じDataﾌｫｰﾏｯﾄになっているのは極めて好感が持てるのだが。何か惜しい。NONMEM を扱うレベルまで自分を引き上げるのは、まだまだ先になりそうだ。

WinBUGs

ギブスサンプリングによるベイズ推定のスタンダードパッケージであるが、使い方がまだわからない。Doodleで簡単にモデル設定できるという話もあるが、とりあえず Model をコーディングするところから。

Web 検索すると、生態系や経済学の研究者で利用者は多い。医薬は少なそうだが、これから大きく入ってくる可能性が高い。既に欧米では PKBugs として Pharmacometrics への応用も増えてきている。

線形混合効果モデルと非線形混合効果モデル

非線形混合効果モデルの前に、線形混合効果モデルを学んでおくことは大切だと思う。特に NONMEM をやるとして、日本ではこの2年注目されているが CWRES の話のように残差計算のことを考えようとすると、どうしても計算原理に踏み込まないといけなくなることがでてくる。数値計算アルゴリズムとまでは決していかないが、数式レベルでの概念の話として。

思うに混合効果モデルのこの辺りの話がわかりづらくなるのは、各被験者の変量効果を固定効果と同時推定しないということではないかと思う。対数尤度を最大化するために反復計算が必要なのだが、その対数尤度関数を求めるときに、変量効果については「考えないようにする」。確率の問題なのだから要は合計してしまうわけだ。これを周辺尤度と呼ぶ。

そうして、後は対数尤度関数を各固定効果で微分したものが0となるような条件を目指して、SAS だったら Newton-Raphson法を用いて反復計算する。各被験者ごとの変量効果は、固定効果が求められたあとに Bayes 推定して求める。固定効果を求めてから、つまり事後分布がわかってから求める量だから、必然的に母集団平均に推定値は偏りがちになる。こういう量を縮小推定量と呼ぶ。

このように線形混合効果モデルには既に NONMEM がやっていることのほとんどが現れている。NONMEM の場合の違いは当然非線形性なのだが、それによって変量効果をまとめて周辺尤度を構成する部分が簡単でなくなってしまう。モデル関数が複雑な分、積分しないといけないのだが、それとて難しい。

そこで関数を簡略化して周辺尤度計算を楽な形にする。変量効果について Taylor 展開してやるわけだ。ここで変量効果を 0 としてその周りに展開すればFO法に、POSTHOC推定値の周りに展開すればFOCE法になる。SAS の場合は級数を用いた数値積分を実装している。Bugs の場合はMCMC法であるが、これがどんなものかはちょっとイメージを持っていない。

後はほとんど同じである。上記の内容を非線形混合効果モデル解析として動態出身の研究者が NONMEM を入り口にして理解することは難しい。ほぼ不可能とさえ思える。非線形最小二乗法の理解のために線形最小二乗法が重要であることと同じであり、統計的素養をゆっくりと醸成して積み上げていくしかない。

しかしこの積み上げを日々しながら思うことは、やはり NONMEM は「早すぎた」ソフトである。動態解析の分野にいきなり持ち込まれたときのショックはどれほどのものだったかと思う。そして、おそらくであるが数値計算には本当はもっと工夫の余地があったんだろうと思う。現在、NONMEM 以外のより高速で安定性の高い、先駆的アルゴリズムを取り入れたソフトの開発がなされつつあり、この流れに NONMEM は今ひとつ対応できていない。

ただ NONMEM は投与スケジュール設定に対する自由さ、微分方程式モデルの作成の容易さのように、臨床PK-PDへ応用するための強力な利点を持っている。

この混合効果モデルの統計理論を考えるたびに、NONMEM とそれをもたらした Sheiner と Beal の先駆性に驚かされる。いまだにこれら統計理論は自分でも咀嚼できていないのだが、ソフトまで解析してリリースした彼らは、臨床PK-PD解析の発展を 20年は加速したと思う。薬物動態の世界でこういうケースは Sheiner だけではないかと思う。