ぎょーむ日誌 2005-02-02

2005 年 02 月 02 日 (水)

• 0800 起床． 朝飯． コーヒー． 0910 自宅発． 晴． 0920 研究室着．
• いろいろなヒトたちのスケジュールのかねあいから， ずいぶんと早い時期に開催される「追い出しコンパ」情報めも． ここに提示する意味はないんだけど， もしかしたらこういうところに置いといたほうが便利かも， という実験 …… というほどでもないか．

　　[日にち]　    2月16日（水）　
　　[会場と時間]  一次会　「ひだまり庭」　　　　19:00-21:00
　　　　　　　　            　住所：北11西3
                            　 tel：011-709-5850
　　　　　　      二次会　「がらり」
　　　　　　　　　　　　　　　住所：北18西4
                            　tel：011-707-7155
　　　　　　　　※地図をお茶部屋に掲示します。
　　[会費]　4000円（一次会）
　　　　　　2300円（二次会）

• ここしばらくの修論・D 論発表など見ていると， R の使用とともに， いよいよ一般化線形モデル (glm) も普及してきた． その「使いたおし」ぐあいはヒトそれぞれで， 単に「なんでも『分散分析』」が「なんでも glm()」 になっただけのもいるし， その一方で 「統計モデルは現象を表現する手段である」 という理解にもとづいて 適切なモデルを構築できる大学院生たちもいる．
• そのなかで特に賢明なる大学院生たちは， 単純な一般化線形モデル (GLM)では 十分に説明しきれない現象があることに気づき (個体差問題など)， そのあたりを一般化線形混合モデル (GLMM) で何とかできるかも， と模索している (ただし混合モデル化が常にこのあたりの解決策というわけではない)． 一年ほど前に R で一般化線形混合モデル のパラメーターを推定計算する関数として glmmML() を紹介 してみたんだけど， このあたりの現状について簡単にまとめてみよう．
• まずは最尤法にもとづく推定用の関数:
  • glmmML() (in glmmML package): Göran Broström 作． CRAN に入ってる． 推定結果に stepAIC() を適用できない (で， 私は自作のモデル選択関数でそれに対応している)．
  • glmm() (in repeated package): Jim Lindsey 作． CRAN に入ってない． 推定結果に stepAIC() を適用できる． AIC なども自動的に計算してくれる．
  後者は苫小牧の平尾君に教えてもらった． 最尤法を使う 利点は本モノの尤度を使っているのでモデル選択などができる， ということ． 欠点は， 現時点ではランダム効果をひとつしかいれることができない (しかも定数のみ)， 推定計算の収束に失敗する可能性が PQL 法より高いかも， というあたりか．
• 次に罰則つき擬似尤度 (penalized quasi-likelihood; PQL) 法にもとづく推定関数 (どちらも CRAN に入っている):
  • GLMM() (in lme4 pacakge): Douglas Bates and Deepayan Sarkar 作．
  • glmmPQL() (in MASS pacakge, nlme package も利用)
  これらの相違はふたつある: 一つはオブジェクトの属するクラスで S4 class (GLMM) か S3 class (glmmPQL) か …… まあ， これは一般ユーザーにはあまり関係ない; 第二点は計算結果の一部が両者で異なる (計算方法が異なるため)． PQL 法を使う利点は簡単に multi-level モデルの推定計算ができる， というもの． 欠点は PQL は尤度ではないので， 尤度にもとづくいろいろな方法 (モデル選択など) が使えない， というところである． GLMM() や glmmPQL() を使った計算結果オブジェクトには最大化対数尤度や AIC (後者のみ) も含まれているけれど， これらは最尤法で得られる尤度とは別モノである．
• これら PQL な GLMM 計算関数は， 「一般化」ではない線形混合モデル (LMM?) を解く関数を利用して擬似尤度を最大化している． 線形混合モデルというのは， 例によって何もかもが等分散正規分布の世界である． しかしこの仮定をうまく利用することで 多重階層なモデルの計算が可能になっている．
• 最尤法を使って multi-level なモデルの推定やりたければ …… 自分で関数をかくしかないと思う． この場合， 重積分はしんどいので， MCMC 法などモンテカルロ法を援用したほうが良いだろう． あるいは， どうせ MCMC 使うならもはやベイズ推定で， というふうになってくるのかもしれない．
• 1000 フィールド科学センターの 前川さん の院生の渡辺君の統計学こんさる． 私が 先日 の修論発表会であれこれと統計学的難クセをつけたため． データの基本構造について教えてもらい， 魚セット (メス 1, 回遊型オス 1, 残留型オス 1) の「個性」を考慮した一般化線形混合モデル (GLMM) で何とかなりそう， という方針を検討する． 今日は生データなく (生態学の悪しき伝統のひとつというべき) 「わりざん値」 な情報欠落ずみデータだけだったんで， また後日に解析をススめることに．
• 北大生協から 大阪大会 の航空券と内地入国の査証が届いたので， カードで支払う． 領収書は三浦さんにお渡しして処理をお願いする． 昼飯．
• なんとなく統計学勉強モードに． GLMM など混合モデルでのモデル選択規準は， 単純にランダム効果のパラメーター数くわえて良さそう， というあたりの再確認 …… まあ， Akaike's Bayesian Information Criterion (ABIC) なんかがそういうカタチになっているから， というのだけでは正当化もへちまもないわけだが．
• さらに統計学勉強が続いてしまう． 昨日書いたことにまたたまた一部まちがいあったわけで …… バハさんの simultaneous autoregressive (SAR) モデルの推定計算， あれで最尤推定になってると思ってたんだけど， そうではなく pseudo likelihood estimator にすぎなかった． 難しい．
• ある限定された SAR で最尤推定を行なう Ord (1975) の方法というのを原論文

  Ord, K. 1975. Estimation methods for models of spatial interaction. Journal of the American Statistical Association 70: 120-126.

  を読みなおして (以前はいいかげんに眺めていただけだった) 勉強してみる． 空間相関のパラメーターを推定するのは， それほど簡単ではない． ひとことで言えば， 距離荷重行列の固有値をすべて計算して …… となるわけで， 複素固有値がまじってきたらどうすんだと思うわけだが， Ord 先生いわく絶対値が大きいところだけ取ればよい， とのこと． まあ， 絶対値の大きい固有値が複素数になったらこの方法はそれでもうダメ になるわけだが．
• 小林さんがアカマツデータ解析をしつつ， 論文の準備にとりくんでおられるようなので， とりあえず依頼された LaTeX ファイル変換して先方に送ってみる． 簡単な Perl スクリプトでフィルタリングしただけのものなんだけど．
• で， けっきょく Ord のアルゴリズムとは
  1. SAR モデルを空間依存項 (推定すべき空間相関なパラメーター ρ を含む) とそうでないところに分割する
  2. 空間に依存してないところ Y = X β + u について ordinary least square (OLS) を適用して 「残差」 u を求める
  3. その u を使って u = ρ W u + ε から ρ の推定計算するわけだが， このときに固有値計算が必要になる
  4. 3. の結果にもとづいて 2. に必要な新しい Y を計算する
  5. 値が変化しなくなるまで 2. -- 4. をくりかえす
  というもののよーだ．
• まあ， ふつーのユーザーはこういう問題に出あったら (以前に苫小牧の平尾君にご教示していただいたとーり) R の spdep package の lagsarlm() 関数でもって計算してしまえばよいわけだが． これはまさに上述のアルゴリズムを実装した関数である．
• 1920 研究室発． 1930 帰宅． 晩飯． 明日の輪読会の予習． いやはや， なかなか手ごわい内容なんで．
• 上記「追いコン」幹事長として東奔西走する尊敬すべき M1 女性大学院生殿から教官どもあての 「嘆願書」と血書したメイルいただき， 「慢性的財政難の学生だけでは、賄い切れません」 とのこととて 「二人分の料金を頂きませんことを、 此処に、切にお願い申し上げる次第でございます」 なる清廉潔白というほかない直訴である． 空腹のあまり勉学すらままならぬ 大学院生たちのまさにまさに窮乏のきわみというべき現状に袖を濡らしつつ， 御心配には及びますまいとの返信したためるうちにも， 以下のような見事なるしめくくりの 感銘のあまりの引用についてのお許しいただかずにはいられないのであった．

なお、当日は出発前に、集金に上がらせて頂く所存でございます。
快く承諾して下さる皆様の懐の深さを信じて疑いません。

• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0830): 米麦 0.5 合． ピーマン卵炒飯．
  • 昼 (1320): 研究室お茶部屋． 米麦 0.7 合． レトルトパウチドカレー．
  • 晩 (2130): 米麦 1.0 合． タマネギ・マイタケ・豆腐・カツオブシの味噌汁． モヤシ・ニラの炒めもの．