ぎょーむ日誌 2008-07-11

2008 年 07 月 11 日 (金)

• 0810 起床． またよく眠れなくてねむい． 朝飯． コーヒー． 0920 自宅発． 雨． 0935 研究室着．
• データ解析こんさるメイルがたまっていたので， いろいろとお返事かき． 衛生管理の微生物学まわりで， なかなかおもしろいポアソン分布の使われかたをされている， と勉強になった．

> （１）添付ファイルのp. 170 の例6-2の問題
> この中で，計数値として＞100といったデータを含めてポアソ 
> ン分布を仮定して解析しております。
> どう解釈して良いのか分からず，悩んでおります。

たとえばデータが N_1 = a, N_150 = b, N_170 = c だったとすると，ふつ
うの尤度方程式では

 (尤度) = (1 個だった結果が a 個存在する確率)
          × (150 個だった結果が b 個存在する確率)
          × (170 個だった結果が c 個存在する確率)

と書くわけですが，この教科書の例題では 101 個以上はカウントしないとしてい
るので，尤度方程式は

 (尤度) = (1 個だった結果が a 個存在する確率)
          × (100 個より大だった結果が b + c 個存在する確率)

を計算しているわけです．上記の確率はすべてポアソン分布から計算されてい
ます．「100 個より大 …… の確率」はポアソン分布で 101 から上を全部積
分しているわけです．

• えーと， 微生物学なのでシャーレ上の寒天培地とかに液体状になった標本をぬって， いくつのコロニーが出現するのかをカウントしているわけだが …… コロニー数が多すぎる場合には「> 100」とか記録されるよーで． で， こういう 「コロニー数 100 まではちゃんと数える， それ以上は (たぶんそういう例が少ないので) いちいち数えない」 という場合のポアソン分布を使った統計モデルで最尤推定するには どうしたらよいか? という問題だ． 上で私が解説してるように， データシートに「> 100」と記録されている場合は 「コロニー数が 101 以上である確率」を尤度とすればよいのである． ふーむ， なるほど．
• さらにこのハナシのつづきとして，

> （２）添付ファイルのp. 171-172，量分析による推定
> 細菌が「増殖する／増殖しない」といった実験結果から
> 増殖の場合の試料中の菌数を推定するという話です。

> 増殖する→1個以上は菌が存在する→菌が0個でない確率
> 
> q = 1- exp(-λ)
> 
> を上記のp = の式に入れております。
> このような方法はどうなのでしょうか？

よく使われる方法です．

> また，もしこのような解析をGLMで扱うにはどうしたらよいで 
> しょうか？

R の glm() 関数の場合，family = binomial(link = "cloglog") を使った場
合に該当します (ただし推定されるのはλではなく log(λ))．

• R で help(family) (あるいは ?family) するとわかるけど， family = binomial(link = ...) の場合， link はおなじみ "logit" (ロジスティック回帰) となるわけだが， それ以外にも"probit", "cauchit", "cloglog" が選択可能である．
• で， この "cloglog" こと complementary log-log なる link 関数はいったい何なのだろうか? たとえば link = "logit" の場合には確率 p と線形予測子 z の関係は
  
  つまりこうなっている．
  
  ところが link = "cloglog" では
  
  つまりこうなる．
  
• 生態学まわりでは発症者の少なくない， 「わからない数式をみても図もつくらずに， そのままふりーづしてますますわからなくなる」 症の大学院生たちのために図を作ってみるとこんなかんぢで．

• で， 使いどころとしては上の質問にあるように， ポアソン分布との連係わざ， たとえば
  • ポアソン分布にしたがっているカウントデータ (のはずなんだけど)
  • 『いない』『いる』データしかない
  • しかしポアソン分布の平均値を推定したい
  といった状況で使える． うーむ， きわめて個体群密度のひくい動物の目撃情報 (「ここらへんでは見た」「見なかった」) データから密度を推定するとか， かな?
• 地衣類繁殖論文 revision 作業のつづき．
• 午後からデータ解析こんさるなので， いつもより早めの昼飯．
• 1300 より 地球圏科学専攻の重光さんのデータ解析こんさる． 海洋沈降粒子の同位体比のハナシとかだったらイヤだなぁ (割り算値の統計モデリングは難しいので) …… と思っていたのだけど， わりと単純なモデル選択の質問だったので， 助かった．
• そしてこのあとも別件のデータ解析こんさるメイルかきが …… さらにネット雑用も 3 件ほどわりこんでくる． さらにプリンタートナーカートリッジの交換． A802 室の CANON LBP5400． 研究室雑用 Wiki に記録 …… ふーむ， 二ヶ月で K トナーがきれたか． このプリンターは場所をとらないんだけど， カートリッジもちっこいんだよね．
• 窓の外はひたすら雨．
• よーやく地衣類繁殖原稿 revision 作業にもどる． えーい， 何がどうなっていたのやら ……
• けっきょく cover letter ひねくるだけの段階は終了し， 原稿本文を変更していかねばならない状況にいる， とよーやく気づいた． 原稿にとりくむ …… こりゃー， なおすべきところがたくさんあるわ ……
• まあ， 文章のなおし自体はそれほど難しくないので， ぢりぢりと進捗していくのだが …… OpenOffice の「わーぷろ」上で作業しなければならぬので， ちょっといらいらする． 私が主著者ならテキストエディター (vim) + LaTeX なんだけど， まあこれはしょうがない． MS わーどじゃないだけマシだと考えることにしよう．
• まだまだ作業は終わらず． 1930 研究室発． 雨． 1945 帰宅． ちょっとうんどう． 晩飯の準備． 晩飯．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 10 × 3 回．
• [今日の食卓]
  • 朝 (0850): 米麦 0.7 合． ワカメ・豆腐の味噌汁． キャベツ・ニラ・ショウガ・豚肉・卵の炒めもの．
  • 昼 (1210): 研究室お茶部屋． 米麦 0.7 合． キャベツ・ニラ・ショウガ・豚肉・卵の炒めもの．
  • 晩 (2100): 米麦 0.8 合． ワカメ・豆腐の味噌汁． キャベツ・ニラ・ショウガ・豚肉・卵の炒めもの．