ぎょーむ日誌 2005-02-12

2005 年 02 月 12 日 (土)

• 0920 起床． 朝飯． コーヒー．
• ふろりだの奥山さんから， 過去数日の GLMM 混乱ハナシに関するコメントいただく (生存解析その他でやるべきでは， という内容でした …… 昨日のぎょーむ日誌をアップロードする前だったもんで)．
• これに関連して， ``Negative binomial loglinear mixed models'' (負の二項分布の対数線形混合モデル) という「ひゃー」というような題名のふろりだ産論文おしえていただいた．

  Booth, J. G., Casella, G., Friedl, H. and Hobert, J. P. (2003). Negative binomial loglinear mixed models, Statistical Modelling, 3: 179-191.

  Poisson 分布をガンマ分布混合でばらつかせたのが 負の二項分布なんだけど， これにさらに正規乱数による random effects (個体差とかブロック差とか) を組み込んだもの． しかも Monte Carlo EM (MCEM) algorithm で計算します …… という何だかスゴいもので． 勉強になりそうだ．
• で， 奥山さんからのご指摘読みつつ返信書いてるうちに， 数日前 に「R の glmmPQL() では死亡率に関してブロック差の中に個体差がある， という状況ではうまく推定できん」 と騒ぎまくっていた問題が解決した …… これは 推定できなくて当然 だったのである． 嗚呼， どうしてこんなあたりまえのことに気づくのに えらく時間かかってしまったりするんだろう (答: 統計学をよく理解してないから)．
• まず， そもそもくだんの biometry 質問への回答は生存時間解析となるべきで， {生, 死} 二値パターンを説明する混合ロジスティック回帰 ってのはあまりよろしくない回答だろう． それはおいといて， glmmPQL() とかで「処理の差に加えて， ランダム効果としてブロック差があってその中に個体差ある」 というモデルで生成した乱数データの ブロック差・個体差をうまく推定できない理由は何だったのか， を考えてみよう．
• ブロックというのはこの場合だと生物をいれとく水槽かなんかで， その水槽のなかに個体が何匹かいる， という状況だ． ブロック差というのは， 実験処理を (同じ水準のなかでは) 均質にしてるはずなのに水槽ごとに何か違いが出てしまう， という差のことである． 個体差は， まあサイズとか遺伝子型の影響で， 同じ毒性にさらされても死にやすいとか死にやすくない， とか．
• 同じ水準の毒をあたえた (つもりの) ときに， ある水槽ではぜんぜん死ななくて， べつの水槽ではばたばた死んだ， という実験結果だったとしよう． このような「ブロックごとのばらつき」 は単純な二項分布モデルからは逸脱している， つまり overdispersion (過分散; 過大分散) が生じている状況である． 簡単に言えば， このときにこういう overdispersion の原因が水槽 (ブロック) 差なのか 個体差なのか原理的に見わけがつかないので， あわれな推定関数 glmmPQL() は錯乱した推定結果を出力していたのである． これはブロック数を増やしたり， あるいはブロック内の個体をいくら増やしても変わらない．
• えー， このあたりは数学の問題というより， まあ常識的なこととゆーか何とゆーか …… ある水槽では実験動物が他よりばたばた死んだ， という状況だったときに，
  • ブロック差説の例: これは水槽にまちがってちょっと多めに毒をいれたせいである (あるいは水循環装置の具合がいまいちだったとか， この水槽を置いてる場所がよくなかったとか)
  • 個体差説の例: この水槽には「毒に弱い」 遺伝子型の個体がたまたま集中してしまったせいである
  などともっともらしく並べてみたのは良いものの， 手もとにある {生, 死} データをいくらヒネくってもどちらが相対的により重要だったか， といったことは推定計算できない．
• これは， 処理の水準とブロックを決めてしまえば， そのブロック内個体の {生, 死} パターンは (個体差があろーがなかろーが) 必ず二項分布にしたがうので (オモテでる確率が 0.5 でないコインであっても， コイン投げ独立反復試行においては オモテでる回数は必ず二項分布にしたがう， というのと同じ)， 個体差の有無やその大小をいえないためだ． 個体差の有無大小がいえないので， ブロックごとに見られる死亡率の相違が， ブロックの影響によるものなのか， その中の個体差によるものなのか観察者にはわからない …… というだけの， ぢつにぢつにあたりまえなハナシで．
• 推定計算してる関数 glmmPQL() の視点にたってみると， 「個体差ぜんぜんナシ」 だろーが 「いや個体差とても大きい」 だろーが， ともかく何であってもデータへのあてはまりかたはほとんど変わらない (実際には乱数生成で作ったデータセットのちょっとした 偏りに左右されて「個体差ゼロ」⇔「個体差最大」 の間のどこかにさまよいこむ)， という状況なんで， ああいう無茶苦茶な計算結果を出していたわけだ．
• いやはや …… {生, 死} データはある期間においては同じ個体から一回しか得られない， というところが特徴的なんだよなぁ． これが複数回観測できる事象なら， ブロック差・個体差を同時に推定できる場合もあるだろう．
• たとえば， その実験動物が元気なときには 一日に何度も X 行動をとる (毒でカラダが弱ると回数が減る)， とする． このときはブロック内の個体数と観測日数が十分であるなら， その個体ごとの X 行動のカウントデータを利用することで， nested なランダム効果である ブロック差と個体差 (in ブロック差) を区別して推定できるはず．
• 1640 自宅発． 曇． 1650 研究室着． 途中で奇妙なモノをみてしまった．

[雪壁上のらくがき]

これは無頼な若者の無軌道な表現活動の帰結ではなく， 北海道警察交通課の 「ここに路駐してた車はレッカー移動しました」 表示なのである． うーむ．

• あ， またズボンのひざの部分が破れた． これは先日ぼろぼろになったのと同じ時期・同じ場所で買ったものである …… このような「壊れる時期がそろう」のは overdispersion の逆の underdispersion なのかもしれない (現象論的には……メカニズムのほうは不明)． まあ， たぶんこちらも破れるだろうと予感してたので， 新しいのは二本買ってあるわけだが． 私は Sony 製品もってないんでよく知らないんだけど， 世に言う ソニータイマー なんかも underdispersion 現象の例なのか?
• いやいや， このあたり Cox の比例ハザードモデルなら， 時間依存な破損率の増大， といったモデルで説明されるにちがいない． 勉強しなくては ……
• その前に， うだうだと上述の 「glmmPQL() 騒動」 の顛末を書く …… なンか最近， 過去のぎょーむ日誌に 「後記」をつける機会が多いような気がする． うかつな日々．
• いやはや． ともあれそもそもの原因は 「個体差だろーがブロック差だろーが， とにかく random effects はどんどん nest して片っ端から混合モデルにつっこめばよいのだ」 という 完全にまちがった 久保ドクトリンにあった， と結論するのが妥当なところでありましょう． いやはや．
• 月曜日の修論発表会にむけて， 修論生をさぽーとする利他的な活動がおこなわれてるよーで， 私も修論生用晩飯のおすそわけにあずかる．
• 2320 研究室発． 2340 帰宅．
• [今日の運動]
  • いかん， 今日も運動休養日になってしまった ……
• [今日の食卓]
  • 朝 (1010): 米麦 0.7 合． ブナシメジ・海藻・豆腐・煮干の味噌汁．
  • 昼 (1430): 米麦 0.7 合． モヤシ・ニラ・エノキダケの炒めもの． ブナシメジ・海藻・豆腐・煮干の味噌汁．
  • 晩 (2030): 研究室お茶部屋． 修論生用晩飯として他の院生がつくったもの． 米麦 0.5 合． ジャガイモ・ナガイモ・ニンジン・鶏肉の煮物． ヤマクラゲ・豆腐の味噌汁．