ぎょーむ日誌 2003-11-(01-10)

2003 年 11 月 01 日 (土)

• 0930 起床． 朝飯． コーヒー．
• うーむ， 時間を有効活用できん ……
• 1430 自宅発． 手持ちの現金が 1000 円未満になってるので南方の ATM へ． 曇． で， いったん JR 札幌駅の南面に出てしまうと書店徘徊してしまうんだよねえ． 本屋が集中分布してるんで． 1630 研究室着．
• Loader 本勉強のつづき． ぢりぢりわかってくる． しかし仕事ススまず．
• 大澤君の解析結果を拝見する． R の glm() 関数 family=poisson で統計モデリング． 結論としては， 苫小牧のシウリザクラに被虐趣味あり． 痛めつけるほど芽の形成個数の期待値が高まる．
• 2300 研究室発． 2310 帰宅．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 18. Does Smoking Cause Cancer?
    • Material implication
  • With some redundancy, Koch was stating the conditions for material implication. This may be adequate to determine that a particular species of bacteria caused the infectious disease. When it comes to something and cancer, however, Koch's postulates are of little value. Let us consider how well the connection between lung cancer and cigarette smoking fits Koch's postulates (and hence Russell's material implication). The agent is a history of cigarette of the lung. There are cigarette smokers who do not get lung cancer. Koch's first postulate is not met. There are some people who get lung cancer who claim that they have not been smokers. If we are to believe their claims, Koch's second postulate is not met. If we restrict the type of cancer to small oat-cell carcinoma, the number of nonsmokers with this disease appears to be zero, so maybe the second postulate has been met. If we take the agent away, that is, if the patient stops smoking, the disease may still come, and Koch's third postulate is not met.
• [今日の運動]
  • 今日もうんどうしてないなあ
• [今日の食卓]
  • 朝 (1000): そば．
  • 昼 (1430): スパゲッティー． タマネギのシチュー．
  • 晩 (2200): 研究室お茶部屋． 特に匿名希望する料理人氏と小菅せんせーが釣ってきた ウサギアイナメ だの カジカ だのといった北海道な魚の味噌汁．

2003 年 11 月 02 日 (日)

• 1030 起床． 朝飯． コーヒー．
• どうも時間がうまく使えない． とはいえ， ぎょーむばかりやってるわけにもいかんが． 11 月は …… よろしくない予感として， 統計学授業準備・釧路大会 web site 準備・その他あれこれ， で終ってしまいそうな．
• 1620 北大構内走． すぐに日没． 一時間ほどゆっくりと走る． 体重 71.4kg．
• 1800 自宅発． 1810 研究室着．
• 粕谷さんからメイルいただく． やはり順位相関つかった多変量解析はいろいろあるようで …… ありうるだろうな．

さて、順位相関にもとづく多変量解析ですが、重回帰や分散分析なども含むなら、 Puri, ML & Sen, PK　の Nonparametric methods in general linear modelsという 本があります。

また、検定するなら、問題によっては MRPP (Multiresponse permutation procedures) とか MRBP (Multivariate rabdomized block permutation procedures) は使えないでしょうか。

• 教えていただいた Puri, ML & Sen, PK, 1985, Nonparametric methods in general linear models なる本は絶版か? 出版社 (Wiley) とかのサイトで見つから． Amazon.com いわく ``out of print''． 北大図書館には あるけど． 検定ねえ …… 三中さん， 無理して「統計学化」せずとも別の生きかたもありえるんでは， というご指摘．
• いつものごとく， まとまった進捗はなく 2220 研究室発． 2230 帰宅．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 18. Does Smoking Cause Cancer?
    • Cornfield's solution
  • An alternative approach is the prospective one. In such a study, a group of individuals is identified in advance. Their smoking histories are carefully recorded, and they are followed to see what will become of them. By 1958, three independent prospective studies had been run. The first (reported by the same Hill and Doll who had done the first retrospective study) involved 50,000 medical doctors in the United Kingdom. Actually, in the Hill and Doll study, the subjects were not followed over a long period of time. Instead, the 50,000 doctors were interviewed about their health habits, including their smoking habits, and followed for five years, as many of them began to come down with lung cancer. Now the evidence did more than suggest a relationship. They were able to divide the doctors into goups depending upon how much they smoked. The doctors who smoked more had greater probabilities of having lung cancer. This was a dose response, the key proof of an effect in pharmacology. In the United States, Hammond and Horn ran a prospective study (published in 1958) on 187,783 men, whom they followed for for months. They also found a dose response.
  • They addressed in two ways the question of whether the opportunity samples were nonrepresentative. They showed the range of patient populations involved, increasing the likelihood that the conclusions held across populations. They also pointed out that, if the cuase and effect relationship holds as a result of fundamental biology, then the patient's different socioeconomic and racial backgrounds would be irrelevant. They reviewed toxicology studies, which showed carcinogenic effects of tobacco smoke on lab animals and tissue clutures.
• [今日の運動]
  • 北大構内走 1620-1720． ストレッチング．
• [今日の食卓]
  • 朝 (1040): 「北欧」バゲット．
  • 昼 (1430): スパゲッティー． タマネギのシチュー．
  • 晩 (2300): マカロニ． タマネギのシチュー．

2003 年 11 月 03 日 (月)

• 1040 起床． 朝飯． コーヒー． 1140 自宅発． 曇． 1150 研究室着．
• Matt が MacOS X を Wake On Lan (WOL) する問題をもちだす …… ちょっと調べてみて， 最終的にはいつものごとくかとー先生に助けていただいて， 問題解決をみた．
  1. A 棟 8F 闇ネットの闇ルータに ssh
  2. wakeonlan -i 192.168.1.255 MACアドレス …… ここでブロードキャスト (192.168.1.255) に magic packet をばらまく， とゆーハタ迷惑なところがぽいんと (かとー先生の発見)
  3. 数秒まつと起きる
  4. terminal から MacOSX を寝かしつけたい場合は SleepNow をインストールして使う
  ただし， SleepNow に関しては， 発令してから避難するまでの時間をかせぐべく， てきとーなシェルスクリプトを書くべきだろう …… sleep 10; ./SleepNow/SleepNow てなかんぢで． おお， かの Apple 社 Macintosh で GNU bash が使えてしまう不思議さよ． いまさらながら．
• アカマツモデルを単純化する試みにとりくむ …… 実験データが例によってゑくせるファイルだったので， これを解毒して解読するしくみ， というのを作っていたんだが …… このレヴェルから手直しが必要なんだけど， なかなか複雑でよーわからん． 自分の 2 年ほど前に作ったプログラムを直すのに数時間を費す． 本日は， 暗中模索を脱した， というぐらいの段階で．
• M2 大澤君の苫小牧シウリザクラデータの解析結果が出はじめてる． なにしろ複雑なデータ構造に帰着されてしまう野外実験なんで， ふつーの検定のたぐいはおよそ役にたたづ， ひたすらに R で glm() & stepAIC() だ． 葉の回転は family = binomial， 芽の個数 family = poisson， 呪われ C-N データ (重量換算値) は family = Gamma． しかし， これはこれで首尾一貫していてよい．
• で， 出てきた結果がじつにじつに興味ぶかい …… stepAIC() すると， 「実験処理 A」「光量子量」 というような生態学的にそれっぽいファクターは 「あっ」というまに消し去られて， 残るのは季節変化だの「個体ごとに定められるパラメーター」 (つまり測定不可能な微環境など「個体差」) だののたぐいばかり …… 苫小牧のようにひどく平らな森林であっても， 野外実験における「処理の効果」などは 「個体差」の中に埋没して観測不可能になるのか? はたまたシウリザクラの「頑健さ」が 人間の傲慢に復讐しているのだろうか?
• あるいは応答が遅いだけ， とか． ともかく， ある程度おおきくなってしまった木本を使った野外実験とは， こういうモノかもしれませんな． 昨年度にもこーいう事例 あるわけでして．
• 2120 研究室発． 2130 帰宅． 体重 71.6kg． 晩飯．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 18. Does Smoking Cause Cancer?
    • Smoking and cancer versus agent orange
  • A contrast to this can be seen in the attempts to indict Agent Orange as a cause for health problems Vietnam War verterans have suffered in later life. The putative agents of cause are contaminants in the herbicide that was used. Almost all studies have dealt with the same small number of men exposed in different ways to the herbicide. Studies in other populations did not support these findings. In the 1970s, an accident in a chemical factory in northern Italy resulted in a large number of people being exposed to much higher levels of the contaminant, with no long-term effects. Studies of workers on New Zealand turf farms who were exposed to the herbicide suggested an increase in a specific type of birth defect, but the workers were mostly Maoris, who have a genetically related tendency toward that particular birth defect.
• [今日の運動]
  • 北大構内走 1650-1740． ストレッチング．
• [今日の食卓]
  • 朝 (1100): シリアル．
  • 昼 : なんとなく食ってない
  • 晩 (2150): 米麦 1.0 合． タマネギ・ニンジン・ニラ・エノキダケ・豆腐・煮干の味噌汁．

2003 年 11 月 04 日 (火)

• 0800 起床． 朝飯． コーヒー． アイロンがけ， とか． あ， 洗濯機うごかしとくの忘れてた． 0850 自宅発． 晴． 0900 研究室着．
• メイルよみかきとか． 釧路大会の要旨の文字数は 800 字とゆーことになったそうで．
• X11 上での日本語入力に関する質問いただいたんだけど， まぬけなことに即答できん． 自分の設定ファイルとかひっくりかえしてながめてみて， $HOME/.bash_profile 中の export XMODIFIERS="@im=kinput2" とか重要かなぁ (その場ですぐに変更したいときは source $HOME/.bash_profile とかが必要)， てなことを思い出していく．
• Matt の統計モデリングコンサルタント． 個体ごとに再測定期間が一定してない毎木データを 線形モデルとして解析したい， というもの． ある期間 Ti において死なない確率を exp(-Mi T_i) と仮定して， Mi を linear predictor にしたら， という素朴なモデルを提案してみる． 簡単に最尤推定可能．
• 生態環境専攻科 (学科のごときもの) の広報委員である 大原さんに送ったつもりのメイルが， 生態環境スタッフ全員にばらまかれてしまった． まぬけだ． しかも漢字まちがってるのを甲山さんに指摘されるし．
• アカマツりはびりに取り組む．
• 午前おわって北大構内走．
• 1300 より講座セミナー， 今日はあやしい M1 堀端君の林床植物のハナシ． あやしい M1 なのであやしい割算値ばかりが示される …… 生態学的現象というより生態学者のあやしい現象を観察する私としては， むしろこのあたりを「よく見られるパターン」と判断し， なぜ生態学者はすぐに割算したがるのか， を考える． まあ， ひとつには情報隠匿だろうな．
• 昼飯．
• 冬は院生たちが論文セミナーやるというので， メイリングリストと連絡用 web page 作る．

ezmlm-make -P -u -f ~foo/ml ~/.qmail-ml foo-ml hosho.ees.hokudai.ac.jp

• アカマツりはびり． ぢりぢりと． なんで時間がかかるのかな．
• 日本語入力問題， かとー先生にご教示いただいて少し理解がススむ． Gtk (Gtk2 のみ?) には GtkIMmodule なるものがあるそーで …… こいつの設定 (/etc/ 中) によっては， kinput2 が動かなかったり． 私のには Gtk2 が入ってないんでわからんかった．
• お茶部屋でアジア旅行雑談． みなさん， いろいろな所に行っておられるようで． と地図ばかりながめる私．
• 2030 研究室発． 2050 帰宅． 体重 71.8kg．
• 晩飯くったらいきなり寝る， 真夜中起床， なるまぬけ周期．
• [今日の素読]
  • 今日はやってない
• [今日の運動]
  • 北大構内走 1210-1245． ストレッチング．
• [今日の食卓]
  • 朝 (0810): 米麦 0.7 合． タマネギ・ニンジン・ニラ・エノキダケ・豆腐・煮干の味噌汁．
  • 昼 (1500): 弁当． 研究室お茶部屋． 米麦 0.7 合． タマネギ・ニンジン・ニラ・エノキダケ・豆腐・煮干の味噌汁．
  • 晩 (2100): 米麦 0.7 合． タマネギ・ニンジン・ニラ・エノキダケ・豆腐・煮干の味噌汁．

2003 年 11 月 05 日 (水)

• 0910 起床． なんと言うか， ですね …… 朝飯． コーヒー． 0950 自宅発． 曇． 1000 研究室着．
• 気がつくと， 一般化線形混合モデル (GLMM) 問題にトラップされている． またしても M2 大澤君シウリザクラ解析がらみで． データ解析ヘルパーから見た 大学院生にはふたとーりに分類できるのかもしれない ……
  1. ほっておくと後がたいへんオーラを発散してる
  2. しばらくほっといても大丈夫だろ感をにじませてる
  うむ， 問題なのは 1. だろうが 2. だろうが， 「最終到達たいへん量」は不変であることだ (この量は個体ごとに異なるけど …… すなわち後述するよーに非線形な混合モデルで表現される)． いましんどいか， 締め切り寸前にたいへんになるか， だ．
• 混合モデルってのはある確率モデルのパラメーターの ひとつもしくはいくつかがまた別の確率論的モデルから得られた， と仮定するモデル． 確率論的モデルのなかで 「無視できんけど関心がない量」 を表現するときに使う． 個体差とか人間の測定誤差とか． 粕谷さんに教えてもらったところによると， これは経験ベイズ的 (empirical baysian) と表現するむきもあるそーで (具体的には 岸野さん とか)． まあ， なんでもいいが．
• R ではこれを nlme ライブラリー中の， lme (確率分布はことごとく Gaussian にしてモデルはひたすら線形) だの nlme (その非線形版) だので扱うそーで．
• なにもかも Gaussian (正規分布) ってのはあまりにも役にたたん， と思ってたんだけど， MASS4 本みてると， 一般化線形モデルを混合モデルにした GLMM も計算できる関数も用意されてる， とわかった．
• ただし混合モデルっつーのはそもそも莫大な計算量を要求するので (尤度方程式を書いてみよ!)， R の中でも扱いは少々めんどーなのである．
  • GLMMGibbs library: 呪われ多重積分を回避するべく Markov Chain Monte Carlo (MCMC) 法を使う --- ただし現時点では CRAN にはない (作者が保守やめたんだろう)
  • glmmPQL 関数 (in MASS library): 罰則つき擬似尤度 (penalized quasi-likelihood; PQL) を使う推定計算法 --- これが良さそうだったんで最初使ってたんだけど， Ripley 御大自身が 「PQL でモデル選択とかやってんじゃねーよ」 と R-help mailing list 上でしつこく威嚇してまわってるのを発見したんで使用中止
  • gee library: McCullagh & Nelder (1989) GLM 本とかにのってる generalized estimating equations (GEE) を計算する方法 --- 使いかたわからずにほうりだす
  • glmmML library: R 上の一般化線形混合モデルぎょーかいの新顔， 定数項ひとつだけ数値積分やって尤度を計算します， という私が以前に Perl とかでやってたのに近い方式 --- なかなか制動むずかしいんだけど， これを使うことに (いわゆる random effect は個体差だけ，なんで)
• で， 使うことになった glmmML の現時点での問題点だが …… 数式オブジェクトの受け渡しがヘン． しょうがないんで， 文字列で渡して as.formula() する (そう，R の如何なる複雑怪奇な操作も文字列操作に必ずや帰着できてしまう， すなわちここが R の「どん底」だ)． それから start.sigma をうまく調節する必要あり． ついでに当然のごとくに stepAIC() とか使えん． これは自作するしかないなあ． いやはやー
• これまたあたりまえのことなンだけど， どの方法をとっても GLMM の推定計算にはえらく時間がかかる．
• 夕方に講座にきた苫小牧 M2 堀田さんに 「久保さん，顔色わるいですよ」と指摘される． 混合モデル呪われです．
• 2130 研究室発． いつものごとく自宅に隣接する北 12 生協で買い物． 野菜果物のたぐいの高価なる北海道であるが (輸送コストが莫大であるため)， 産地が近いのでリンゴだけは安い． 1 個 50 円以下． 2200 帰宅． 体重 71.4kg． 晩飯．
• [今日の素読]
  • 混合モデル呪われでなんとなくやってない
• [今日の運動]
  • 混合モデル呪われでなんとなく走ってない
• [今日の食卓]
  • 朝 (0920): 米麦 0.7 合． 納豆．
  • 昼 : 混合モデル呪われでなんとなく食ってない
  • 晩 (2220): 食わなかった弁当の残り． 米麦 0.9 合． シュンギク． 目玉焼き．

2003 年 11 月 06 日 (木)

• 0920 起床． ううう …… ひとつ明らかなのは夜中に起きてても仕事もススまなければ， 休憩にもなってない， というところだ． ということで「2300 消灯」なる習慣が定着するかどうか， しばらく試みてみるか． 朝飯． コーヒー． 1000 自宅発． 雨． 1010 研究室着．
• 朝から夜まで R プログラミング． 一般化線形混合モデルの数値積分関数 glmmML に stepAIC みたいなことをやらせる機能の構築． 先日の最尤折れ線選択より格段に難しい． ちょっと油断してると， すぐに予想してない方向に vector 展開しやがって …… ですね．
• 午後は農学部に 松木さん の博士論文発表聞きにいく． 農学部地図上では隠匿されてる共同研究棟なる新しい建物の 1F． 1530 開始． 1 時間ほどで無事に終了．

[D 論発表会会場]

どうもこの日は数人分まとめて発表会をやってたらしいとはいえ， すげーおおいりだった (松木さん曰く，内地からもわざわざ聞きにきたヒトありとのこと)．

• カバノキ科 (ケヤマハンノキ，シラカンバ，ダケカンバ，ウダイカンバ，アサダ，サワシバ) の葉の被食防御 (虫に食われない努力……毒物・とげ・trichome など) の種特性を野外調査と栽培実験で調べました， というもの． どーもハナシの圧縮度が高すぎるよーな， という感触だったんだけど …… 実際のところ， 北方食害マフィアにおける 松木直系子分の大澤君が教えてくれたところによれば， やはりこの内容の数倍に及ぶ成果が得られてるそうで． たしか， 博士課程に入ってからこの方面の研究に着手したはずなのに …… おそるべし．
• ところでその大澤君のさらに下働きとしてこき使われてる私などは 北方食害マフィアの準構成員ということになるのか?
• 研究室にもどって準構成員労働のつづき． 下請け， といえば農学部でお会いした苫小牧の鍋嶋さんに 直径成長解析のハナシを申しつけられたよーな ……
• プログラミングがうまくできなくて， 今日の仕事はススみませんでした， という挫折をひさびさに味わいつつ本日は撤退． 2030 研究室発． ちょい雨． 2050 帰宅． 体重 71.6kg． 晩飯．
• 2300 消灯．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 18. Does Smoking Cause Cancer?
    • Publication bias
  • This was one of Fisher's accusations. He claimed that the initial work by Hill and Doll had been censored. He tried for years to have the authors release detailed data to back up their conclusions. They had only published summaries, but Fisher suggested that these summaries had hidden inconsistencies that were actually in the data. He pointed out that, in the first Hill and Doll study, the authors had asked if the patients who smoked inhalted when they smoked. When the data are organized in terms of ``inhalers'' and ``noninhalters,'' the noninhalters are the ones with an excess of lung cancer. The inhalers appear to have less lung cancer. Hill and Doll Claimed that this was probably due to a failure on the part of the respondents to understand the question. Fisher scoffed at this and asked whi they didn't publicize the real conclusions of their study: that smoking was bad for you, but if you had to smoke, it was better to inhale than not to inhale.
• [今日の運動]
  • うんどうしてない三日目 …… えーい
• [今日の食卓]
  • 朝 (0910): 米麦 0.7 合． 大豆カレー．
  • 昼 : 今日もなんとなく食ってない． あ， お茶部屋リンゴ (余市産レッドゴールド) かじった
  • 晩 (2130): 米麦 0.7 合． 大豆カレー．

2003 年 11 月 07 日 (金)

• 0750 起床． 朝飯． コーヒー． 0845 自宅発． 晴． 0855 研究室着．
• また R プログラミング．
• Matt が先日のモデルのパラメーター推定計算がうまく収束せん， と言ってくる …… それは私には予想されてた困難だったんで， 線形予測子を exp(線形予測子) とすれば， と回答． これでうまく計算できるはずだ．
• で， そっちは片づいたんだけど， Mathematica 使いたいと言いだしたんで， それに対応せんといかん (前からも言ってたし)． 北大は Mathematica キャンパスライセンスっつーのを 取得してたはずで， それならばかなり安く使用可能なはず ……
• ところが， いつものごとく， というか北大計算機 & ネットワーク関連の できの悪いサイト (実際のところ死ぬほどめちゃくちゃで 学部学生が演習とかでいやいや作らされた web site のごとし …… あんたらそれでも 情報処理 の専門家か? と詰問したくなるしろもの) で右往左往な探索してよーやく 発見 したわけだが …… 専攻科単位で申請しろって? また意味不明なことを． とはいえ， 抗議するのも疲れるので岩熊専攻科長に 「Mathematica 利用申請だしたいんですが」 おうかがいメイルをだす．
• Matt 計算， optim() の Nelder-Mead と BFGS の計算結果が一致しないって? しかも simulated annealing (SANN) がコケる， と． 彼の計算命令 をちょい手直しして Nelder-Mead と SANN が一致することを示す (あたりまえだ……どちらも目的関数の微分に依存しないんだから)． っつーことで， あんたの導関数計算がまちがっとるんでしょーが， ということで落着．
• これで午前終了ですよ， ということで北大構内走． 天気はいいけど寒い． 本日は道内何か所かで「遅い」初雪． 札幌はまだだけど．
• 夕方ごろに， よーやくにして一般化線形混合モデルを glmmML で計算させた結果を stepAIC する関数， glmmML.stepAIC ができた (後記: これはむしろ stepAIC.glmmML と命名すべきだったなぁ) …… 苦闘した． とりあえず乱数を使った試験運転クリア．
• しかし， シウリザクラ観測データを投入する勇気がないので， お茶部屋で時間をつぶす．
• が， いつまでも先のばしにはできんので， 葉っぱ死亡速度を推定させる計算スタート． さすがに時間かかる (ThinkPad240Z Cerelon 500MHz (遅!) だと 20 分ちょい) …… 計算結果の 例． なかなかよろしい． ばてた．
• 苫小牧の鍋嶋さんから直径成長データ送っていただく． これもまた glmmML.stepAIC おくりになってしまうのか? なんでもかんでも混合モデルにすりゃあいいってもんじゃないけど．
• で， シウリザクラ当事者たる M2 大澤君がなかなかスルどい指摘を …… たしかに stepAIC 系のアルゴリズムは全面的に信用してはいかん， という可能性がある． というのも， 一般化線形モデル (たとえば二項分布 + ロジスティックモデルなどなど) において，
  1. 個体差を区別しない glm + stepAIC
  2. 個体差を個体 ID の名義変数化で表現してみた glm + stepAIC
  3. 個体差をランダム変数として混合モデル化してみた glmmML + glmmML.stepAIC
  これら， 1.-3. において得られた結果が常に必ずしも無矛盾とはかぎらん， というところ． たぶん 1. と 3. の比較は素朴にやればよろしい． 問題なのは 1. と 2. だ． モデル 1. は 2. に「含まれている」はずだから， 1. なんぞは省略してしまえ， という考えかたもありうるわけだが …… じつは計算結果をみると個体差無視したモデル (1.) のほうが AIC が低い (つまり良い) のに， 2. においては個体 ID が名義変数として「選択」されてしまった， という状況がありえた! ということ．
• このあたり， 詳しく調べてみる必要がある． とうめんの対応策としては， 愚直に 1. から 3. までを「予選」的に計算してみて， つぎに各計算における最良モデルたちのあいだで「決選」 すれば問題はない． (後記: 問題ありあり …… 翌日の記述参照)
• 2030 研究室発． 2040 帰宅． 体重 71.6kg．
• 2310 消灯．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 18. Does Smoking Cause Cancer?
    • Fihser's solution
  • The epidemiologists were using the tools Fisher developed for the analysis of designed experiments, such as his methods of estimation and tests of significance. They were applying these tools to opportunity samples, where the assignment of treatment was not from random mechanism external to the study bun an intricate part of the study itself. Suppose, he mused, that there was something genetic that caused some people to be smokers and others not to smoke. Suppose, that this same genetic dispostion involved the occurence of lung cancer. It was well known that many cancers have a familial component. Suppose, he said, this relationship between smoking and lung cancer was because each was due to the same event, the same genetic dispostion. To prove his case, he assembled data on identical twins and showed that there was a strong familial tendency for both twins to be either smokers or nonsmokers. He challenged the others to show that lung cancer was not similarly genetically influenced.
  • On one side there was R.A. Fisher, the irascible genius who put the whole theory of statistical distributions on a firm mathematical setting, fighting one final battle. On the other side was Jerry Cornfield, the man whose only formal education was a bachelor's degree in history, who had learned his statistics on his own, who was too busy creating important new statistics to pursue a higher degree. You cannot prove anything without a randomized experimental designs, said Fisher. Some things do not lend themselves to such designs, but the accumulation of evidence should prove the case, said Cornfield. Both men are now deceased, but their intellectual descendants are still with us. These arguments resound in the courts, where attempts are made to prove discrimination on the basis of outcomes. They play a role in attempts to identify the harmful results of human activity on the biosphere. They are there whenever great issues of life and death arise in medicine. Cause and effect are not so simple to prove, after all.
• [今日の運動]
  • 北大構内走 1220-1300． ストレッチング．
• [今日の食卓]
  • 朝 (0810): 米麦 0.7 合． 大豆カレー．
  • 昼 (1320): 研究室お茶部屋． リンゴかじっただけ 大豆カレーは高カロリーなのか， いったん食うとなかなか空腹にならん．
  • 晩 (2120): 米麦 0.7 合． 大豆カレー．

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): キャベツ・ニンジン・タマネギの焼きそば．

2003 年 11 月 09 日 (日)

• 0810 起床． 朝飯． コーヒー．
• じつは昨日のうだうだたるぎょーむ日誌は 本日になってうだうだと書いたものである． うだうだ．
• 1315 自宅発北大構内走． 曇天． 帰路の自宅ちかくでちょい雨． うわー， 白いものがまじってる． みぞれぎみだ． 体重 71.6kg． 減らん． 昼飯． アイロンがけなど．
• 1530 自宅発． 歩いて 300m ほどの投票所に行ってみる． 開校 103 年目の 北九条小学校 で公民権とやらの行使． 私と同時刻に投票したヒトは 3 人ぐらい． 空いてる． NHK の出口調査に捕まる． ヨドバシカメラで電池買う． 単4×8本 620 円． 1620 研究室着．
• かとーさんとのすごい会話．

  くぼ

  あ，かとーさんが日曜日まで……お疲れさまです

  かとー

  依田さんの測定関係は片づきましたよ

  くぼ

  なるほど， で， いまは MacOS X の hacking を

  かとー

  そうなんですよ …… MacOS X 用 OpenOffice で使える僕が作った日本語化 patch を今の開発 version に取りこみたい， という依頼があってですね

  くぼ

  ああ， かとーさんペイジで 公開されてるのを 見つけたんでしょうね …… MacOS X 用は開発者が少ないとはいえ， すごいことですね

  かとー

  そう， 検索して見つけたらしいんですよ …… 先方は日本語読めないのに

  くぼ

  え?

  かとー

  その MacOS X 用 OpenOffice 開発者は 外国人で日本語読めないヒトなんですよね

  くぼ

  ということは， C のコードの断片を見て ……

  かとー

  そう， そのようです

  くぼ

  実に驚くべきことですね …… hacker と hacker の間には自然言語は不要， プログラミング言語のコードを見せあえば， おのずと意志は疎通できる， と ……!!

• しばらく放置しておいたメイル雑用など片づける．
• 昨日分のうだうだした記述はところどころ不正確でもあるな． 一般化線形モデルは離散確率分布モデル (二項分布・Poisson 分布など) でも計算できるし， 連続確率分布 (Gamma 分布・Gaussian 分布など) も OK (ついでに言えば擬似尤度も使える)． ということで …… 離散確率分布モデルの一般化線形モデルの結果 (glm() で family が binomial や poisson などの場合) と一般化線形混合モデルの結果 (glmmML() など) は (おそらく) 比較不可能 …… と書くべきだろうな．
• 今日は大澤君が来てるんで混合モデル議論． 結局， ほとんどすべての統計モデリングは混合モデルになってしまいそうだ． prunusdb に説明を追加．
• 2100 研究室発． 2120 帰宅．
• 2330 消灯．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 19. If You Want The Best Person ...
  • Gertude Cox came to study at Iowa State during this time. She had dreams of becoming a missionary and saving souls in far-off lands. For almost seven years after high school, she had dedicated her life to social service in the Methodist Church. She needed a university education to apply to the missionary service that interested her. Snedecor convinced her that statistics was more interesting, and she stayed on after she graduated to work with him at the Statistical Laboratory. In 1931, she received the first master's degree in statistics granted by Iowa State, and Snedecor hiered her to teach in his department. She became particularly interested in Fisher's theories of experimental design and taught the first courses on experimental design at Ames. Snedecor found a place for her in the graduate program in psychology at the University of California, where she continued her studies for two more years. She returned to Ames with her doctorate, and Snedecor put her in charge of the Statistical Laboratory.
• [今日の運動]
  • 北大構内走 1315-1410． ストレッチング．
  • 腹筋運動 30 × 3 回．
• [今日の食卓]
  • 朝 (0840): 蕎麦． 納豆 & コンブ．
  • 昼 (1450): スパゲッティー． キャベツ・タマネギ・ニンジン・ジャガイモ・豆腐のシチュー．
  • 晩 (2200): キャベツ・タマネギ・ニンジン・ジャガイモ・豆腐のシチュー．

2003 年 11 月 10 日 (月)

• 0730 起床． 朝飯． コーヒー． 0820 自宅発． 晴． 0830 研究室着．
• stepAIC.glmmML() を改造して family = binomial(logit) だけでなく family = poisson() も指定できるようにする．
• 呪われ成分分析は nlme するしかないのかなぁ …… ということで nlme を調べる．
• 計算機センターより連絡あり， いまさら (つまり年度の途中では) Mathematica のライセンス発行してやるつもりナシ， なる連絡うける． すばらしー硬直ぶりで． まあ， 私は使わんからいいけど．
• PC よせあつめ計算機ばなし． 宣教師殿は何か祭壇のたぐいを必要とするものらしい． そういえば東京の某集合計算機会社はまだ存続してるんだろーか， と調べてみると …… あった． あいかわらず， いかにも理系な見苦しい web page で．
• といったことをうだうだとハナしてるうちに， 下野さんが関東から無事帰還．
• これで午前終了． 北大構内走．
• シウリザクラ観測データ処理， 最後の面倒 (であってほしい)． 2002 年の個葉自動解析過程はできてるんで， 2003 年データをこれに適合させるもの． 将来 H 社ソフトウェアぎょーむ者たる M2 大澤君がススめてきたデータ転換作業が まだまだハンパというか …… とゆーか， H 社人事部 (いや，リスクマネジメント部か?) におたくの新入社員予定者が こんなコード書いてますぜ， と転送したくなったりするわけだが．
• で， 転換もできたんで， これをデータハンドラー + 統計モデリング方面に流していく． その過程のエラーチェック機構で「観測データのバグ」 が見つかっていく …… 入力ミスのたぐいはじつは前後左右をみれば 自動的に間違いを特定できるできるものが多い． ただし， たいはんの生態学者はこういうことすらやらずに 「解析」しておられるわけだが． 今年はデータシートの構造を改善しましたってことで (まあ，そのせいで上記データ変換なるよけーな作業は必要になるんだが)， さすがに 2002 年度データよりは間違い数が 1/5 ぐらいに減った． 合計 20 箇所ぐらいだ．
• 途中でばてたので， 晩飯は院生たちと北大ちかくで外食．
• 2002 年と 2003 年の観測データをまとめて パラメーター推定したりモデル選択できるように， Makefile その他を書き換えていく．
• まだこまごまとしたことは残ってるとわかりつつも， 今年度修論における「たいへんなこと」は終った， と信じたい． これでもう終ってるはずだ． きっとそうに違いない． 2230 研究室発． 2240 帰宅．
• 明日以降の懸案あれこれ ……
  • 個体ごとの葉の生き死に数変化
  • 統計学授業の準備
  • カエデ論文原稿
  • 苫小牧直径成長
  • parametric vs nonparametric
  • 某土建コンサル無料下請け
  • アカマツりはびり
  • 直達光と散乱光
  書き出してみると， かえってスケジュール破綻， という実感がつのってくるなぁ．
• 2400 消灯．
• [今日の素読]
  • うー
• [今日の運動]
  • 北大構内走 1310-1340． ストレッチング．
• [今日の食卓]
  • 朝 (0740): 米麦 0.7 合． キャベツ・タマネギ・ニンジン・ジャガイモ・豆腐のシチュー．
  • 昼 (1400): 研究室お茶部屋． なんとなく食う気がしなかったので， お茶部屋リンゴ (本日は紅玉) をかじる．
  • 晩 (2030): 地下鉄北 12 条駅すぐそばの CHELSEA (review 1 & 2) で札幌ふうスープカレー． 豚角煮カレー 800 円．