ぎょーむ日誌 2003-11-08

2003 年 11 月 08 日 (土)

• 0800 起床． 朝飯． コーヒー． 起床時間が比較的早くてもだらだらとすることに変わりなし．
• いまさらながら R 的には glmmML.stepAIC は stepAIC.glmmML とすべきだったなぁ， と気づいてみたり． まあ， そのうち修正しますか．
• 1255 自宅発北大構内走． 曇． そしてかなり寒い． 体重 71.6kg． 昼飯． 1430 自宅発． 北 9 東 4 のホーマックにいき， トレッキングシューズふう靴と衣類を購入． さいきん， ホーマックで購入できるものはなんでもホーマックで買ってしまう， という状況になってるよーな． 1530 研究室着．
• また一般化線形混合モデルにとりくむ． 結論として， 昨日かいたことはいろいろと間違ってる．
• その前に， 同じデータに R で立てつづけに glm() + stepAIC() するときの注意． 原因はまだ完全には理解できてないんだけど， glm(..., data = ...); stepAIC(...) というのを連続させるのはまずいようで， attach(what = data); glm(...); stepAIC(...); detach(...) としたほうがなぜかしら良いようだ． 少なくとも前者ではうまくいかない場合がある (けどその状況が把握できてない)．
• さて， 昨日の一般化線形モデル比較のハナシ，
  1. 個体差を区別しない glm + stepAIC
  2. 個体差を個体 ID (個体識別番号) の名義変数化で表現してみた glm + stepAIC
  3. 個体差をランダム変数として混合モデル化してみた glmmML + stepAIC.glmmML
  なんだが …… なお， ここでいう「個体差」ってのは 個体のサイズだのその場の明るさだの遺伝子型だの といった人間に観測された量ではなく， その場の栄養塩類の濃度だの地形のちょっとしたでこぼこだの といった人間には観測されなかった けど何か個体に影響を与えてる部分のことだ．
• まず， 1. と 2. を比較した場合， 個体差が無視できない観測データにおいては， だいたいにおいて 2. で個体 ID が説明変数として選択され， それは当然のごとくに 1. よりも AIC が高い． つまりこの問題に関しては stepAIC() によるモデル選択には矛盾はない と言っていいと思う．
• 蛇足ながら， R の version によって glm() の挙動が少しだけ異なる． R-1.7.1 では「(個体数が多いのに) 個体 ID を名義変数とする」 というような乱暴きわまりないことをやると， いくつかの ID では最尤推定値が NA になってた (原因不明)． R-1.8.0 にすると以前より計算時間が長くなるけれど， ことごとく推定値をだす． なお， 最大化対数尤度や AIC は version 間で違いはなかった．
• Vine Linux-2.6 用の R-1.8.0 が準備されてる …… 日本にある CRAN ミラーとしては 会津大学 とか．
• ともあれ R における一般化線形モデルにおける変数選択， すなわち glm() + stepAIC() (上でいう 1. と 2.) には問題ない， と判明． 問題は 3. のような一般化線形混合モデル との比較だ．
• これまた昨日の言明がひっくり返るわけだが …… 離散確率分布モデルの一般化線形モデルの結果 (glm() で family が binomial や poisson などの場合) と一般化線形混合モデルの結果 (glmmML() など) は (おそらく) 比較不可能， だと憶測する． 理由は混合モデルのほうには尤度方程式に random variable に関する連続確率分布をふくむから．
• 観測データへの確率論的モデルのあてはまりの良さ， つまり尤度ってのは便利なものではある． たとえばいささか乱暴な気もするけれど， 異なる確率分布を仮定したモデル間の定量的な比較でも使える． ある非負連続値データのばらつきに対して， 正規分布 or ガンマ分布のあてはまりがどちらがよいのか， とか (連続 vs 連続)． 同じように離散分布モデル間の比較， ポアソン分布と負の二項分布のどちらがどの程度マシなのか， といった比較 (離散 vs 離散) など．
• しかし， だ． 連続分布モデルと離散分布モデルの間 (連続 vs 離散) では比較できん． そもそもこんな比較には意味ないわけだが …… もし無理にやったとしてもダメである． 離散確率分布モデルに対して定義される尤度はホントの確率だ． すなわち 「ある確率論的モデルを仮定してる状況でその観測データが得られる確率 L」 という定義を完全に満たしている． しかしながら， 連続確率分布モデルで計算している尤度は 「(上の「……」内で定義された) L に比例する量」 つまり L × constant である． ふつー， この constant なる部分は考えない．
• 例を考えてみよう． ある無限集団における人間の身長が正規分布 N(μ, σ²) に従うとする (この確率分布関数は F(x | μ, σ²) で確率分布密度関数は f(x | ...) とでもしますか)． このモデルが正しいときに， その集団から身長 177cm なるサンプルが得られる確率は? 答はもちろんゼロである． 連続分布とはそういうものだ． さて， ここで Δx が小さい量としてみよう． これを使って連続分布上のある区間に サンプルが存在する確率を計算できるわけで …… 身長が 177cm から 177 + Δx cm に存在する確率， というのがよーやくゼロより大きい値 F(177 + Δx) - F(177) ≈ f(177)Δx となる． 上でいう constant とはこの Δx のことだ． 対数尤度になおすと log(f(177 | μ, σ²)) + log(Δx) となり， log(Δx) は最尤推定においてはパラメーター {μ, σ²} に依存しない項となるから無視して差しつかえないわけだ．
• ということで， 離散分布モデルと連続分布モデルでは尤度の意味が違うので， これらの値を比較しても意味がない．
• しかるに， 今回のシウリザクラあれこれ推定においては 一般化線形混合モデル (GLMM) の 二項分布版だのポアソン分布版だのを使ってるわけだが …… これが 一般化線形モデル (GLM) の二項分布版やポアソン分布版の推定結果と 比較できない， と言ってるのは 前者は離散確率分布+連続確率分布のハイブリッドモデル (離散確率分布モデルの中に正規分布モデルがたたみこまれたもの) であるのに対して， 後者は離散確率分布だけからなるモデルであるから (離散*連続 vs 離散)， ということ． どちらも尤度なる数値をだすけれど， 前者は確率に比例する量で後者は確率そのもの． ということで比較不可能， というのが私の現在の理解である． (後記: これって実はかなり間違ってるかも)
• もし無理矢理にでも比較する方法を考えるなら …… GLMM はこれ以上は簡単にはならないから， GLM をややこしくするしかない． うーん …… 何か Bayes 推定的な解釈をいれて， すべてのパラメーターが何かの確率分布にしたがっている， と仮定して …… すると GLMM では random variable が正規分布にしたがってる， としている部分を GLM では一様分布 (げーっ) とかにしたがってると仮定して …… うん， まともなデータ解析者であれば一様分布の推定なんかは 絶対にやりたくないはずだよね (嗚呼，しかし極値分布の推定とかで何かあるのかも)． ということで一様分布仮定の Bayes 的 GLM は推定不可能， この路線では GLMM と GLM の比較も不可能． めでたしめでたし …… とめでててはいかんな．
• random variable にしてる部分を 局外母数 (nuisance parameter) なるものに指定して， これを除いた profile likelihood を …… という作法もあるのかもしれんけど， これはいまいちよくわかってない．
• ともかく一般化線形混合モデルは避けられないという気がするんで， これを導入するにあたって， 何かもっと別の正当化が必要とされるわけだ． うだうだした楽しくないものになりそーだ．
• これまたとうめんの対策として， こういうのはどうだろう．
  1. 一般化線形モデルで推定を行う． もしここで (人間には観測されなかった) 個体差が無視できないばらつきをもつ量であるならば， 個体 ID が説明変数として選択されてしまうだろう．
  2. もし 1. で個体 ID が説明変数として選択されてしまったら， そのモデルは捨てる …… 理由は， 上の仮定では個体差という量が 何か一様分布のごとき確率分布にしたがうと 仮定していることになるから．
  3. あらためて個体差が正規分布などの確率分布に従う， と仮定しなおし， 一般化線形混合モデルを使った パラメーター推定・モデル選択を行う．
  すっきりせんけど， このへんで妥協するか …… 混合モデルシミュレイションとかやって， 個体差のばらつきと 1. におけるモデル選択の挙動を調べるべきだろな．
• 計算結果をみながら， 以上のようなことをうだうだと考えて終った． 今日は幸か不幸か， かかる問題を惹起しているシウリザクラ観測データ， それを 2 年間かけて集めちゃった大澤君は来ていない． そこで (データ隠匿の第一人者として著名なる) 小菅せんせーにこういう問題でひっかかっている， といったら 「そんな瑣末なことどもはさっさと解決して (隠匿せよってこと?)， データ解析全体を片づけてしまえ」 というような説教をされてしまった ……
• 2110 研究室発． 2130 帰宅． 晩飯．
• 2400 消灯．
• [今日の素読]
  • やってない
• [今日の運動]
  • 北大構内走 1255-1335． ストレッチング．
• [今日の食卓]
  • 朝 (0830): パン． キャベツ・ニンジン・タマネギ・卵の炒めもの．
  • 昼 (1350): キャベツ・ニンジン・タマネギの焼きそば．
  • 晩 (2150): キャベツ・ニンジン・タマネギの焼きそば．