ぎょーむ日誌 2005-02-09

2005 年 02 月 09 日 (水)

• 0830 起床． 朝飯． コーヒー． 0930 自宅発． 曇． 0940 研究室着．
• 昨晩はまりこんでしまった混合モデル推定問題のつづき． どうも個体 in ブロックという二階層 random effects はなかなかうまく推定できんようなので， とりあえずブロックの効果だけのモデルに単純化してみる．
• random effects をひとつにすれば， 最尤法を使った混合モデルの推定ができる． 乱数を使ってブロック差・個体差のあるデータを生成する． ブロック (plot) のランダム効果は rnorm(1, mean = 0, sd = 0.5) で， 個体のそれは rnorm(1, mean = 0, sd = 0.1) という設定にしている． いつものことながら， この同じデータを使っていても glmm() (repeated package) の推定結果と

  Coefficients:
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)    0.586      0.258    2.27   0.0231
  toxicity.L     3.646      0.653    5.58  2.4e-08
  toxicity.Q    -0.840      0.612   -1.37   0.1701
  toxicity.C    -1.482      0.635   -2.33   0.0196
  toxicity^4     1.176      0.642    1.83   0.0669
  toxicity^5     0.306      0.635    0.48   0.6301
  sd             0.764      0.266    2.87   0.0041
  

  と glmmML() (glmmML package) の推定結果

                coef se(coef)      z Pr(>|z|)
  (Intercept)  0.580    0.291  1.992  4.6e-02
  toxicity.L   3.614    0.822  4.396  1.1e-05
  toxicity.Q  -0.834    0.679 -1.230  2.2e-01
  toxicity.C  -1.471    0.710 -2.071  3.8e-02
  toxicity^4   1.168    0.709  1.648  9.9e-02
  toxicity^5   0.305    0.695  0.439  6.6e-01
  Standard deviation in mixing distribution:  0.726 
  

  はびみょーに異なる (AIC はだいたい同じ)． glmmPQL() (MASS package) だとこんなかんぢで．

  Random effects:
   Formula: ~1 | plot
          (Intercept) Residual
  StdDev:      1.0456  0.81762
  Fixed effects: dead ~ 1 + toxicity 
                Value Std.Error DF t-value p-value
  (Intercept)  0.5811   0.27448 84  2.1173  0.0372
  toxicity.L   3.5539   0.66403 36  5.3520  0.0000
  toxicity.Q  -0.7927   0.66284 36 -1.1959  0.2395
  toxicity.C  -1.4142   0.67645 36 -2.0907  0.0437
  toxicity^4   1.1011   0.68171 36  1.6152  0.1150
  toxicity^5   0.2699   0.67646 36  0.3989  0.6923
  

• 同じデータをもちいて， こんどはブロック (plot) 差だけでなく個体差 (id) も考慮して glmmPQL() で推定させると (推定のさせかたは昨日のぎょーむ日誌 参照)，

  Random effects:
   Formula: ~1 | plot
          (Intercept)
  StdDev:  5.9131e-06
   Formula: ~1 | id %in% plot
          (Intercept)   Residual
  StdDev:      6.7657 4.4077e-06
  Fixed effects: dead ~ 1 + toxicity 
                Value Std.Error DF t-value p-value
  (Intercept)  1.5665   0.61762 84  2.5364  0.0131
  toxicity.L  12.9624   1.51286 36  8.5681  0.0000
  toxicity.Q  -2.0533   1.51286 36 -1.3572  0.1832
  toxicity.C  -4.1781   1.51286 36 -2.7618  0.0090
  toxicity^4   3.0496   1.51286 36  2.0158  0.0513
  toxicity^5  -0.8537   1.51286 36 -0.5643  0.5761
  

  このように， むちゃくちゃになる (とくに random effects)．
• これはいくら何でもひどい …… ということで glmmPQL() の使いかたを間違ってる可能性を検討してみる …… うーむ， ヘンなことはやってないハズ， なんだが． しかし， ブロック内個体数を多くする設定のもとでも， 事態は改善されんのだよねえ． (後記: これは glmmPQL() の使いかたが間違ってるとゆーより， 私のモデルのたてかたが良くなかった， とわかった …… これでは推定できなくて当然 …… 詳しくはこの週末のぎょーむ日誌)
• さらに， だ …… 上のいずれの推定においても fixed effects の推定値もよくわからないものになっている． これは toxicity なる処理を順序因子 (ordered) として表現してる影響のよーで …… たとえば (Intercept) をここに書いてあるとーり， と理解してはいかん． うーむ， 選点な直交多項式 (orthogonal polynomial) め．
• ナゾに満ちた選点直交多項式系， John Fox ``An R and S-PLUS Companion to Applied Regression'' によると， これが使える必要条件は
  1. 水準間が「等間隔」
  2. 各水準における標本数が同じ
  だそーだ． 1. は知っていただけど， 2. は知らなかった． 嗚呼， 苫小牧某修論とかどうなることやら …… ともあれ， R における順序因子の取りあつかいは， いっそう慎重でなきゃならん， とわかった． このあたり， 単純ではないよーで， 慶應大の統計学講義ノート PDF ファイルには 「あえて非順序因子として」 うんぬんと書かれていたり， あるいは John Fox は ``Helmert contrasts may also be of interest for an ordered factor.'' などと注記している --- R では Helmert 対比は非順序因子 (factor() な連中) の default の coding として使われてるアレ．
• とはいえ， いま乱数例を作って調べてる問題は， 上記の 1. も 2. も完ペキに満たしてるんだけどねえ …… だれかさっさと biometry に模範解答を投稿してくれ．
• たしかに Helmert 対比のほうがわかりやすい …… しかし， 矛盾する結果が出やすいのもこちらだろう．
• 三浦さんの ThinkPad に Norton AntiVirus インストールする作業． 二回も再起動を要求された． あいかわらず， ゐんどーづはヘボい．
• 研究費の残り 2800 円， か． いつのまにか消失してますなぁ． 昼飯．
• 昨日の私の軌道おんちぶりを見かねて， プラネテス でも鑑賞して勉強しなさい， という助言メイルいただく …… とある女性大学院生から． これが理系わーるどだ．
• インドネシアなデータ解析の下請けぎょーむ． 発注者側， なんかだんだん意味不明になってきた． とりあえず， 計算プログラムまわりの作業は終了． あとは引き渡して終りにしますか． とりあえず， 1740 当方の作業は終了． ここからは撤退．

• …… いや， 待て． この計算に関する文書かき作業が残されていた， か …… 作業終了 1810． オワった．
• 1850 研究室発． 1900 帰宅． 運動． 晩飯． 明日の輪読会の予習．
• お． 粕谷さんが biometry に投稿． しかし現時点では指摘 3 点のみ， か．
  1. GLM こそが不等分散な状況をあつかう
  2. 「不等分散ならのんぱら」信仰はまちがい (cf. 先月末のぎょーむ日誌)
  3. 生存時間のモデル化したいんだったら負の指数分布とか?
  …… 負の指数分布ってあったのか …… どうも調べてみると， ふつーの指数分布をそう呼んでる場合もあれば， 生存時間解析で使ってるあの分布 (私が勝手に「二重な指数分布」と読んでるやつ) のことのよーで． しかし， おそらくこの問題は時間分布ではあるまい． (後記: またしても， まちがい …… 時間分布のハナシである)
• [今日の運動]
  • エアロバイク 50 分間．
• [今日の食卓]
  • 朝 (0850): 米麦 0.5 合． 卵焼．
  • 昼 (1350): 研究室お茶部屋． 米麦 0.6 合． ジャガイモ・ニンジン・タマネギ・大豆のシチュー．
  • 晩 (2130): 米麦 0.7 合． ジャガイモ・ニンジン・タマネギ・大豆のシチュー．