「調査観察データの統計科学」3.1章 傾向スコアの数式メモ（前半）

【追記】
社内勉強会資料を整えて公開しました。

• 「調査観察データの統計科学」読書会資料を公開しました（数式周りをフォロー） - 木曜不足

【／追記】

みどりぼん(「データ解析のための統計モデリング入門」)を読み終わったから、というわけではないが、同じ岩波・確率と情報の科学シリーズの「調査観察データの統計科学」（星野崇宏）を読んでいる。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

社内で週一開催している勉強会の自分の担当回でもこの「調査観察データの統計科学」を紹介。今は3章の傾向スコアの途中。
勉強会向けに資料も作っているが、closed ということでいろいろ遠慮なくやらかしており、そのままだとちょっと公開できない（苦笑）。どうせ傾向スコアまでいかないと内容無いので、3章終わったくらいに再構成して公開しようと思う。

最近は、社内のマーケティングなデータ解析っぽい仕事もしてて、先日開催された WebDB Forum 2014 の技術セッションでは「クラウドサービスにおけるユーザ動向予測について」といった発表をさせてもらった。資料は残念ながら公開できないのだが、概略を簡単に説明すると、BtoB クラウドサービスのように、利用できるデータに制限があってもデータ解析をするには、という感じの話。
そうした「いわゆるデータサイエンス」なお仕事をちょびちょびしていると、やりたいことの半分は「意思決定」なのだろうけど、もう半分は「施策の効果を測定したい」で占められるのでは、と思えてくる。ただしその「施策」もコストがかかるとなると、できるだけ効果の大きい（と希望的に推測している）ところを狙うことになるわけで。
そんな「無作為抽出した実験群」うんぬんが夢物語な状況でも公正な効果測定を行いたい。そういう欲張りな要件に、この本の言う調査観察研究が役立っちゃったりなんかするのかも？　と、これが動機。

前置きが長くなったが、今日の本題。
そういうノリで読み始めた「調査観察データの統計科学」だが、この本の数式周りには「むむむ？」というツッコミどころがちょいちょいある。特に傾向スコアが導入される 3.1章は、わずか2ページちょっとながらこの本の中でも重要度が高いだろうと思われるので、ここの数式は疑問をなくしておきたい。
そのあたりをフォローするメモをまとめてみた。

まずバランシングスコアを定義する。
b(x) がバランシングスコアであるとは、 が成立することとする(条件付き独立の記号が出ないので、縦棒一本で勘弁)。
なお Notaion は本の通りなので、新たに自前の記号を導入でもしない限りここでは一切説明しない。潜在的な結果変数が〜因果効果が〜ってやってたらキリないし。あくまで、「調査観察データの統計科学」を読んでて困ってる人向けのメモである。

続けて、「バランシングスコアは、関数 g を使って と表現できることが必要」とあって、式 (3.2) によってそれを証明しようとしているが、これは残念ながら証明になっていない。
A=B を証明せよという問題で、 A=B を使って「この等式が成り立つには A=B が必要である」と答え、「それは証明ではありません」と×される、という流れに身に覚えのある人もいるかと思う。それである。A=B を証明するのに A=B を使ってはいけない。
式 (3.2) を逆に見ると、 を仮定して つまり が成り立つことを示すことになら使える。よって、 の十分性なら 式 (3.2) によって証明できる。

それを踏まえつつ、そもそもこの命題は何のためにあるのかというと、実は次ページで登場する「傾向スコア」が を満たしていることから、傾向スコアはバランシングスコアであることを示すのに用いる。つまりここで示すべきは必要性ではなく十分性の方なのだ。
というわけで、ここは式 (3.2) はそのままに、その説明を「バランシングスコアは、関数 g を使って と表現できることが十分条件である。というのも であるとき 式 (3.2) が成立するからである」とし、次のページに「先の証明から傾向スコアがバランシングスコアであることがわかる」という一文を追加すれば論理がつながる（後半は傾向スコアを定義したあとにもう一度説明するのでご心配なく）。

ちなみに、実はこの命題は必要十分で、本の脚注にて「十分条件については Rosenbaum and Rubin(1983)の定理 2 を参照」と付記されている。まあ上での指摘を入れれば、参照すべきは必要条件の方であることは言わずもがなだろう。
せっかく勧めてくれたので Rosenbaum and Rubin(1983) "The Central Role of the Propensity Score in Observational Studies for Causal Effects" の定理 2 を見てみると、必要性を示すのに背理法を用いていた。めんどくさい。
ここは「 x を止めたとき b(x) は一意に決まるので、p(z=1|x,b(x))=p(z=1|x) が常に成り立つ」ことを認め、一方で b(x) がバランシングスコアであることから p(z=1|x,b(x))=p(z=1|b(x)) でもあるので p(z=1|b(x))=p(z=1|x) が言える。よって g(b(x)):=p(z=1|b(x)) とおけばよい。これで必要性の証明になる。背理法とかいらない。かんたん。

この関係式 p(z=1|x,b(x))=p(z=1|x) は本には明記されていないが頻繁に使用されており、例えば先の式 (3.2) の中でもバッチリ使われている。
でも今「 x を止めたとき b(x) は一意に決まるので、p(z=1|x,b(x))=p(z=1|x) が常に成り立つ」とさらっと流そうとしたけど、これもし数学科の演習だったら見逃してもらえない。「それ本当に成り立つの？」「はい、明らかに……*1」「明らかなら証明できるよね？（にっこり）」と真正面から刺されて見事な窮地に陥る。
まったく細かいことであり、傾向スコアを使いたいだけなら無条件に認めても困らないのだろう。が、ここは「3.1章の数式から疑問をなくすメモ」なので、道筋くらいはきっちりつけておこう。

といっても、この「簡単で難しい関係式」が「難しげ」なのは記号 b(x) のせい、ただそれだけである。
b(x) は確率変数でもあり x の関数でもあるかのように都合よく使いまわされている*2。要は手抜き。普段はそうしてズボラしているが、何がどうなっているのか説明しろと言われ、手抜きのしっぺ返しを食らってる、って寸法だ。
これは確率変数と関数をちゃんと分離した記号を導入するだけで解消できる。

今、b は関数 b(x) を表す記号とし、それとは別に確率変数 B を p(B|x)=1 (B=b(x) のとき)、p(B|x)=0 (B≠b(x) のとき) と定義する*3。この記号の下で、今まで p(z=1|x,b(x)) と書いていたのは実は p(z=1|x,B=b(x)) のことだったのだ。これが p(z=1|x) に一致することは確率の加法乗法定理でさくっと出るので以下略。記法重要。

ここまででやっと1ページ進んだが、それなりの分量になってしまい。
3.1 の残りの1ページのためのメモは更に長いので、残りは別記事に切り分けることにする。
社内勉強会もこんな調子でやってるので、なかなか進まない(笑)。1回のセミナーで1行しか進まないとかまあよくある話なんである。数学屋さんは。

*1:ここで「本にそう書いてあったから」と答えたら退場

*2:確率論の公理で「確率変数とは関数である」的な話とは全く別なので混同しないように

*3:B の動く範囲(= b の値域)が離散ならこれでいいが、連続なら本来はディラックデルタ使って指示分布？　を構成しないといけない。実際、傾向スコアは連続値をとる。が、計算が煩雑になるだけで本質的には全く変わらないので、理解を得る分にはこの定義で十分だろう

とりあえず plot だけでもしてみるのススメ #みどりぼん

10/21 に開催された「データ解析のための統計モデリング入門」（以下みどりぼん）読書会の最終回にのこのこ参加。主催の yamakatsu さん、参加者＆発表者のみなさん、会場を提供してくださったドワンゴさん、ありがとうございました。
懇親会はちょっと断念した。無念。

• 「データ解析のための統計モデリング入門」読書会 最終回 あの日作ったモデルの名前を僕達はまだ知らない - connpass

最後なのでちょっと口はばったいことを言ってみる。
WinBUGS がインストールできなくて試せなかった的な話もあったが、参加された60人ほどの方たちでサンプルデータをとりあえず plot だけでもしてみたって人はどれくらいいただろうか。
みどりぼんにもちろん plot 図はもともと載っている。ただ（学習という面で考えると特に）残念なことに、全ての図に「正解」の点線や、「正解モデル」で推定した分布などが重ねられており、生のデータのものがない。

たいしたことではない。久保先生のサイトで配布されているデータを R で load して plot するだけだ。慣れてれば1分、慣れてなくてもまあ5分くらいの作業。
個人的には RData で配布されるのはちょっとめんどくさい。中を気軽に見れないし。50個の整数データくらいなら、テキストの方が扱いやすくて嬉しい。
というわけで、11章のデータをテキストにしたものを貼り付けておこう。これならコピペして5秒だ。

Y <- c(0,3,2,5,6,16,8,14,11,10,17,19,14,19,19,18,15,13,13,9,11,15,18,12,
       11,17,14,16,15,9,6,15,10,11,14,7,14,14,13,17,8,7,10,4,5,5,7,4,3,1)
plot(Y)

みどりぼんで見た「正解入り」 plot 図と比べて、ずいぶん印象が違う。本当は「正解」なんて知っているはずがないので、実際に目にすることができるのはこちらの方だ。
「正解」を知らないでこの図を見たとしたら、と考えると、正解以外の解釈の可能性がちらちらよぎらないだろうか。例えば「両端は外れ値だな！」とか。
ちなみに、このデータの平均値は 10.9 であるのに対し分散は 27.4 もある。みどりぼんが言うところの「過分散」が起きており、単純なポアソン分布ではモデリングできない → 階層ベイズで空間構造や！　というのが 11章のストーリーなわけだが、両端の5点を外れ値とみなして捨てれば、平均 12.7、分散 16.3 と過分散はかなり抑えられる。図もこうなる。

おっとこれって右下のもう2点捨てれば……みたいな作為的な後出しジャンケンは統計の嫌うところではあるが、試してみるのも面白いと思うし、すぐに試せる（から、ここではこれ以上やらない）。

みどりぼんはいくつものモデルを紹介してくれているし、この本で紹介されていないモデルももちろんまだまだたくさんある。その数多あるモデルの中から、実際の場面ではどのモデルを使うべきか決めてくれる論理的な根拠というものは、残念ながら存在しない（せいぜい消去法。例：過分散だから生ポアソンは×）。だから、そこの判断は人間が適切にやるしかない。
データを見て、データに関する事前知識とすり合わせ、「ふむふむ、どうやら空間構造があるかも？(ドヤ」とか推測し、11章のモデルを使ってみるところにたどりつき、実際に試してみて、空間構造を入れた場合と入れてない場合とでナントカ IC を比べたりして、ビミョウな結果に「やっぱ外れ値かも……」とか凹んじゃうわけだ。
でもそれってモデルの上っ面の知識だけでできることだろうか。データを愛で、解釈やモデルを取っ替え引っ替えし、ハマった場合とハマらなかった場合のモデルの挙動に一喜一憂したことがなくてできることだろうか。

LT で berobero11 さんが「みんなもっと plot しよう！　WinBUGS しよう！」(意訳)とおっしゃっていたとおり、ホントもっと plot しよう。
WinBUGS は確かにセットアップがいろいろめんどくさい（特に環境によっては）が、JUGS だって Stan だってある。
みどりぼんをせっかく最後まで読んだのだから、意義のあるものにして欲しいと期待。

「データ解析のための統計モデリング入門」6.6章 割算値はなぜダメなのか？ #みどりぼん

昨日 7/29 の「第6回「データ解析のための統計モデリング入門」読書会」は参加しなかったのだが、ニコ生で中継してくださっていたので後半を聞くことができた。
6.6章「割算値の統計モデリングはやめよう」では、タイトルの通り観測データ同士を割り算するなと話しているわけだが、読んでいていろいろ疑問に思うところがあり。
読書会中継でちょうど 6.6 章以降を担当された 0kayu さんの発表を聞くことができたのだが、気になっていたあたりは特に質疑でも話題にならず残念。
というわけで、誰かがツッコミを入れてくれることを期待して自分の疑問をここに書いておく。

「データ解析のための統計モデリング入門」6.6章では統計モデルに観測データ同士の割り算値を持ち込むことを批判している。その理由として、

• 比率にすることで元のスカラー値の情報が失われる
• 値それぞれが分布を持っている場合、それらの割り算値の分布がよくわからない

があげられている。
が、前者は元の値をモデルから除かずに割り算項を加えるだけならいいとも言えるし、後者の理由なら直前の 6.5章で導入した交互作用項（つまり掛け算）もアウトなはずである。
つまり、上の2つではその理由付けには足りていないのではないか（それとも何か誤解してる？）のが疑問の1つ目。

6.6.1 では「ついつい割り算したくなるようなデータ」でも、オフセット項という「わざ」を使えば割り算を回避できるという話を具体例を交えて紹介している（どうしてここで「わざ」という言葉使いをするんだろうという疑問もあるのだが、さすがに枝葉末節かなとも思うので置いておく）。
しかし、紹介されている具体例を見ると、

という式で表現されている。これは要するに をモデリングしているわけだ。
しかしそれと を比べると、割り算がモデルに入ってしまう後者は「絶対ダメ」と言われてしまうほど本質的に違うものだろうか。

最も大きな違いは、y は整数値(植物個体数)なので、 ならポワソン分布でモデル化できるけど、調査地の面積 A で割った はそういう訳にはいかない、という点にある。
0以上の実数値を扱う分布と言えば真っ先に思い浮かぶ1つはガンマ分布だが、図ったように みどりぼん 6.8 章はそのガンマ分布を使った GLM の話なわけで、これは試しにやってみろということだろう。

を、y/A がガンマ分布に従うとしてパラメータ推定したのがこちら。

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/binomial/data4b.csv")
(model.poisson <- glm(y~x, offset=log(A),family=poisson,data=d))
(model.gamma <- glm(y/A~x, family=Gamma(link="log"),data=d))

> (model.poisson <- glm(y~x, offset=log(A),family=poisson,data=d))

Call:  glm(formula = y ~ x, family = poisson, data = d, offset = log(A))

Coefficients:
(Intercept)            x
     0.9731       1.0383

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    261.5
Residual Deviance: 81.61 	AIC: 650.3

> (model.gamma <- glm(y/A~x, family=Gamma(link="log"),data=d))

Call:  glm(formula = y/A ~ x, family = Gamma(link = "log"), data = d)

Coefficients:
(Intercept)            x
     0.9699       1.0516

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    5.753
Residual Deviance: 1.863 	AIC: 192

上がオフセット項付きポワソン分布、下がガンマ分布。
AIC が出てるとついつい比べたくなるが、それはいくらなんでもダメに決まってる（ヒント：片方だけ連続分布）のでグッと我慢しつつ、推定された係数を見ると非常に近い値が得られている。
線形予測子をプロットしてみても、ほとんど差がないことがわかる（分布は描くのが面倒なので略）。

plot(d$x,log(d$y/d$A))
abline(model.poisson,col="red")
abline(model.gamma,col="blue")

赤がポワソン、青がガンマ。
「結果だけ」見ると、どちらかがどちらかより確実に良い悪いとは言えなさそう。

みどりぼん p133 には「このようにオフセット項を使うと、個体数平均は調査地の面積 A_i に比例する、といった仮定を反映させつつ明るさ x_i の効果を推定できます。個体数を調査地面積で割って密度にする、といった観測値同士の割り算は全く不要です」とあるが、割り算しないことがまるで目的化しているようにも読めて、個人的にはちょっと気持ち悪い。
ここでポワソンのほうがいいんだろうなと判断する根拠としては、ガンマ分布の方は y が整数値しか取らないことをモデリングできていないからであって、割り算が直接要因ではない、はず。

つまり割り算を避けるべき理由がまだいまいちよくわかってない。
おそらくは割り算値を導入した「全然ダメダメなモデル」をたくさん見てきたことから「割り算ダメ」という経験則が生まれたのだと思うのだが、もうちょっと納得できる理由付けはできるんだろうか、というのが2つ目の疑問。
「割り算を使った特徴を組み込む場合はその妥当性に注意をはらいましょう」くらいの言い方なら特に異論ないのだけど、任意の特徴についても「妥当性に注意をはらいましょう」は言えちゃうので、何も言ってないのに等しいか（苦笑）。

そもそも、「比率にすることで元のスカラー値の情報が失われる」という理由を持ち出すのなら、「割り算したくなる値」が目的変数にあたるこの例は適していない気がする。
そして「割り算したくなる値」が説明変数に来る場合は「オフセットわざ」が使えないのだが、その時はどうするのだろう（それでも割り算はダメ？）、というのが最後の疑問。

【追記】

「データ解析のための統計モデリング入門」6.6章 割算値はなぜダメなのか？ #みどりぼん - Mi manca qualche giovedi`? に書いた疑問にコメント・フォローをいただいた。ありがとうございます。

「データ解析のための統計モデリング入門」6.6章 割算値はなぜダメなのか？ - URL @shuyoさんから#みどりぼん §6.6 は「なぜ割算値がダメなのか」の説明・例示になってないとのご指摘．たしかに，もっとうまい例示が必要ですね…

2014-07-31 16:08:01 via Twitter Web Client

• 統計モデルに観測値と観測値の割り算値を入れても問題ない: ニュースの社会科学的な裏側
• 統計モデルに観測値と観測値の割り算値を入れたときのバイアス - 餡子付゛録゛

ぶっちゃけると「割り算はダメ」も「割り算してもいい」も額面通りに信じるつもりはもともとなくて、「○○○○な場合は割り算は筋が悪い」くらいの話なのだろうな、と思っていた。みどりぼんはそこの「限定」がうまく表現できていないのかな、と。

現時点で確実に言えることは、割り算固有の話としては「0やその前後を含む分布が除数に来るような割り算はダメ」という算数レベルの常識くらい。
あとは「ヒストグラムや散布図(R で言う pairs)を描いてみて、これはヤバそうという場合は使わないでおきましょう」という、全ての特徴について言えることを守っておけば、割り算をあえて避ける理由はない（少なくともまだ言語化されていない）という理解で大丈夫そうな雰囲気。

ブートストラップの適切なサンプル数 -「データ解析のための統計モデリング入門」第5章 #みどりぼん

7/8 に開催された「データ解析のための統計モデリング入門」、通称「みどりぼん」の第 5 回読書会にのこのこ参加。主宰のやまかつさん、発表者＆参加者の皆さん、会場を提供してくださったドワンゴさん、ありがとうございました。

• 第5回「データ解析のための統計モデリング入門」読書会 - connpass

今回は第５章、尤度比検定。
みどりぼんは、この前の回の AIC も含めたモデル選択について比較的ゆるふわな説明で、これはこれでありなんだけど、しっかりしゃっきりしたものも読んでみたいなあ（でも数式少なめでよろしく！）という人は、伊庭さんの「モデル選択とその周辺」がおすすめ度高し。

• 伊庭幸人, モデル選択とその周辺

検定は（有用だろうけど）あんまり好きになれなくて、気分的な盛り上がりは今ひとつなのだが、知らなかったら好きも嫌いもないよね！　ということでまじまじ手を動かす。
勉強会の質疑応答も、前回参加した第３回よりおとなしめだったかな。モデル選択って自明な例見せられても実感わかないよねー。

そんな中、パラメトリック・ブートストラップについて、p104 の脚注には「実はサンプルの個数は 10^3 ぐらいでは十分なサイズではありません」とあるが、適切なサンプルサイズはいくつ？　そりゃあ増やすほどいいんだろうけどどこまで増やす？　とかいう質問が上がっていた。
サンプルを増やすほど推定量の精度が上がる（＝分散が下がる）ことは、直感的にわかるだろうし、事実その通りなのだが、残念なことに分散の下がり方は遅い。
実際、一番シンプルな問題として「サンプルの平均」を平均の推定量とする問題を考えるとき、簡単な計算ですぐわかるのだが、サンプル数を N 倍するとその推定量の分散は 1/N に、つまり標準偏差は 1/√N となる（分布の形が漸近的に正規分布に一致することを言うには中心極限定理が必要だけど、ここではそこまでの話は求めていない）。
つまり頑張ってサンプルを 100倍にしても精度は 10倍にしか上がらず、どこかでコストが見合わなくなるわけで、そこが「どこまで増やすか」の答えとなる。

と、いうのが机上のお話なわけだが、「1000 じゃあ足りないよね〜」と訳知り顔でうなずいている人たちは常にそういう事を意識してるかというと必ずしもそういうわけではなく、もっと単純に 1000 では足りない事を経験として知っているのだ。

みどりぼんの 5.4.1 では、パラメトリック・ブートストラップで 1000 個の逸脱度の差を得て、その分布の 95% 点を 3.953957 と推定するお話が出てくる。
ランダムサンプリングしている以上、実行するたびに得られる値は異なる（つまりこの推定値は確率変数）わけだが、さてはて、実際どれくらい結果に幅があるのだろうか。
コードが全て掲載されているのだから、実行して確かめてみればいい。

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv") 
get.dd <- function(d) {
  n.sample <- nrow(d);
  y.mean <- mean(d$y);
  d$y.rnd <- rpois(n.sample, lambda=y.mean)
  fit1 <- glm(y.rnd~1, data=d, family=poisson);
  fit2 <- glm(y.rnd~x, data=d, family=poisson);
  fit1$deviance - fit2$deviance;
}
pb <- function(d, n.bootstrap) {
  replicate(n.bootstrap, get.dd(d));
}

ここまで準備をしておけば、あとは quantile(pb(d, 1000), 0.95) と叩くだけで毎回異なる推定値が得られる。
ぜひ、答えを見る前に、みどりぼんの中の値 3.953957 と比べて、どれくらい推定値に幅が出るか予想してみよう。
テキトーでいいよ〜。

---

OK？　では実際に5回実行して得られた5つの推定値がこちら。

> quantile(pb(d, 1000), 0.95)
     95% 
4.165724 

> quantile(pb(d, 1000), 0.95)
     95% 
4.312189 

> quantile(pb(d, 1000), 0.95)
     95% 
3.513216 

> quantile(pb(d, 1000), 0.95)
     95% 
4.064914 

> quantile(pb(d, 1000), 0.95)
     95% 
3.751253 

おそらく多くの人が「そんなに変わるんかー……」と目を丸くしてくれたのではと期待しているのだが、どうだろう。
これでもう次からは「1000 じゃあ足りないよね〜」といっしょにうなずく仲間入り。

ただし、みどりぼん 5.4.1 の目的は「逸脱度の差 4.5 の p 値が 0.05 より小さいか」を確認することだったのを思い出すと、上の推定値は最大でも 4.3 ちょっとと、p<0.05 という結論には影響ない。
つまり、今回の問題にとって必要な精度は 1000 サンプルでも確保できてそう（まあ、もう少しくらい精度高いほうが安心できるだろうけど）。
しかし逆に、もうちょっときわどい値だったら、例えば「逸脱度の差 4.0 の p 値が 0.05 より小さいか」という問題だったときに、1000サンプルで推定した値がたまたま 3.953957 で、「あーよかった、ギリギリ超えたから棄却できた！」なぁんて言っちゃうたら……。
そういったことを適切に判断することが一番大事で大変だったりする。

「データ解析のための統計モデリング入門」第3章メモ #みどりぼん

6/10 に開催された「データ解析のための統計モデリング入門」、通称「みどりぼん」の第3回読書会にのこのこ参加＆ LT してきた。主宰のやまかつさん、発表者＆参加者の皆さん、会場を提供してくださったドワンゴさん、ありがとうございました。

• 第3回「データ解析のための統計モデリング入門」読書会 - connpass

LT は……正直、ネタを思いついた時は、ええやんおもしろいやん、とかなぜか思いこんでた。当日の朝くらいに「測度論とかないわー」ということにやっと気づき、一週間前の中谷の後ろ頭をしばき倒したくなったが、もはや後戻りはできず。この LT をちゃんと最後までやり通したという一点で褒めてあげて欲しい……

猫に教えてもらうルベーグ可測 from Shuyo Nakatani

さて気を取り直して。
質疑応答の時にあれこれ好き勝手なことを突っ込ませていただいてたのだが、読書会終了後にそのあたりをブログとかにまとめてあると嬉しい、という声をもらったので、かんたんにメモしてみる。

3.6 説明変数が数量型＋因子型の統計モデル

この節では、体サイズ x と施肥 f の両方を組み込んでポアソン回帰をしてみたら、施肥の係数 β3 が負になっちゃった！　という話をする。
しかし、「みどりぼん」はこの教科書的おあつらえ向きな展開にもかかわらず、「前の節(=施肥のみを説明変数とするモデル)では肥料の効果 fT がプラスであったのに、このモデルではマイナスだと推定されています。肥料の効果についてはいよいよわからなくなってきました」とコメントしたのみで終わってしまう。
案の定、読書会ではこの点が議論の的となり、みんなでこの「直感と異なる結果」をどのように解釈するのかを熱く語り合うことになった。このように議論してもらう事を期待して、あえて説明しなかった……？

先に結論めいた事を一言でいうと「説明変数である体サイズと施肥効果の間に、いわゆる弱い相関があるから」となるだろうか。

一般化線形モデルでは説明変数間の独立性をどこか期待している部分があって(「説明変数」「目的変数」ではなく「独立変数」「従属変数」と呼ぶ流儀もあり、このことと関係している……のかな？　由来は調べてない)、そうでない場合はパラメータが「期待する値」から外れたり、精度が悪くなる、などと言われたりする。
といっても、読書会でも指摘のあったことなのだが、「期待する値」とは「説明変数にその変数だけを使って、他の変数の影響を除外した場合の係数」だと解釈することができる。すると、その変数が他の変数と独立なら「期待する値」に一致するだろうし、逆に独立でなければ、その「期待する値」に一致しないのは不思議でもなんでもない。

つまり、パラメータが「期待する値」から外れているように見えて、その「理由」を考察する必要をもし感じたなら、まずは説明変数間の相関を確認するのが一つの手。
そこで体サイズと施肥効果の関係を見てみると、

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv") 
plot(d$f,d$y) 
plot(d$f,d$x) 

左図が施肥の有無ごとの種子数の関係(みどりぼん図3.3)、右が施肥の有無ごとの体サイズの関係。種子数は施肥の有無と関係なさそうだが、右の体サイズはどうも関係がありそうに見える。
実際 C=0, T=1 として各変数間の相関係数を見てみると(TRUE を 1 、FALSE を 0 とするために +0 している)、

> cor(d$y,(d$f=="T")+0)
[1] 0.01914445          # 種子数 と 施肥の有無

> cor(d$x,(d$f=="T")+0)
[1] 0.2806593           # 体サイズ と 施肥の有無

> cor(d$x, d$y)
[1] 0.2276909           # 体サイズ と 種子数

となり、施肥と種子数は相関がなく、施肥と体サイズ、体サイズと種子数には弱い相関がある、と一般に言われるような相関係数が得られる。
つまり、施肥と種子数の間にも（かろうじて）正の相関があるが、体サイズとの相関で十分説明できてしまっている。したがって施肥による効果は微調整レベル、負になっても不思議はない。そういう解釈をすることもできる。

また、勉強会の質疑で出ていた意見として、「同じ体サイズなら、施肥してそうなったものより、施肥をせずにそうなったもののほうがポテンシャルが高いと考えれば、施肥の係数が負になることは『直感に反しない』」という解釈が納得感が高くておもしろかった。そういう解釈をすることもできる。

一般化線型モデルに限らず統計モデルを現場で使ってみたときの「あるある」ネタとして、得られたパラメータの値、特に正負が「直感と異なる」場合に（上司などから）理由の説明を求められて困る、というのは鉄板の１つだと思うので、「みどりぼん」のような教科書でここらへんフォローしてくれててもいいのに、という気もちょっとする。
そういった時に上のような作文が求められるわけだが、逆にそういう主観的で作為的な作文術に結びついてしまうことを嫌って、この話題を深追いしたくなかったのかもしれない。

ただそういう作文は置いといたとしても、「データ解析で最も重要なのは、まず何はともあれ、そのデータをさまざまな方法で図示してみることです」(p16)とか、「統計モデリングにとりくむときに、そのデータを「いろいろな図にしてよく見る」点は何度でも強調しておきたいところです」(p44)とか、図示の重要性をたびたび強調する「みどりぼん」なら、上で描いた体サイズと施肥の関係を図示しといてほしかったかなあとは思うかな。

あと、p51 で z 値や Pr(>|z|) といった指標の説明をしておきつつ、施肥を説明変数に加えてパラメータの効果が「いよいよわからなくなってきました」というとき、なぜこのモデルでの z 値などを示さなかったのかも疑問。まあモデルの選択は次の章の AIC に全て託した、ということなのかなあ。

ちなみに z 値類を出してみるとこんな感じ。

> fit <- glm(y~x+f, data=d, family=poisson)
> summary(fit)

Call:
glm(formula = y ~ x + f, family = poisson, data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3977  -0.7337  -0.2023   0.6795   2.4317  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.26311    0.36963   3.417 0.000633 ***
x            0.08007    0.03704   2.162 0.030620 *  
fT          -0.03200    0.07438  -0.430 0.667035    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 89.507  on 99  degrees of freedom
Residual deviance: 84.808  on 97  degrees of freedom
AIC: 476.59

Number of Fisher Scoring iterations: 4

fT の z 値は -0.430 ( 0 から十分離れていれば有効なパラメータ)、Pr(>|z|) は 0.667035 ( 0 に近いほど有効なパラメータ)ということなので、fT は有効なパラメータとはなさそうだ、ということが３章の知識からだけでも推測できる、と言っていいのかな？

ちなみに AIC は、上の summary にも出てるが 476.59。
施肥を加えなかった場合の AIC は、

> fit <- glm(y~x, data=d, family=poisson)
> summary(fit)
(中略)
AIC: 474.77

なので、こちらからも施肥を加えない方がいい、という結論に。

PRML ガール 〜とある文芸部の統計女子〜

これは「PRMLガール 〜 文芸部のマネージャーが「パターン認識と機械学習」を読んだら 〜」の幕間的なお話です。
未読の方は先にそちらをどうぞ。

• PRMLガール 〜 文芸部のマネージャーが「パターン認識と機械学習」を読んだら 〜 - 木曜不足

---

「先輩、先輩。我らが PRML ガールの評判聞きました？」
「聞いてない」
「やっぱり気になりますよねえ。ねえ？」
「僕は別に」
「……」
「……わかったわかった。どんな評判なの？」
「評判というか苦情ですかね。『 PRML ガール言うくせに PRML 関係ないやん！』みたいな」
「まあ、情報量とエントロピーの話しかしてないし」
「『あとがきがわりの AC ガール』（暗黒通信団刊行の書籍『PRML ガール』に収録）に至っては選択公理ですからねえ。そんな評判を跳ね返すためには、紛れもなく PRML ってテーマをここらで一発取り扱っておくべきなんじゃないか！　と愚考するわけですよ。このあたくし」
「そういうことを言い出すと、黒幕がウォーミングアップを始めるからやめて欲しいのだけれど……それで本音は？」
「あの〜、PRML 8 章のグラフィカルモデルでどーしてもスッキリしないとこがあって〜」
「はいはい。余計な前振りいらないから。普通に聞けばいいのに」
「矢印のある方のグラフィカルモデルなんですけど」
「ベイジアンネットワークだね」
「はい。あれの矢印の向きに、本質的な意味なんか無いような気がしてきちゃってて……なんのためにあるのかなあ、と」
「なるほど。どうしてそう思ったか……の前に、ベイジアンネットワークの矢印の定義は？」
「むー、さすがにそれくらいはわかってますよー」
「わかっているなら答えられるよね」
「はうっ。……えーとぉ、改めて聞かれると自信なくなるんですが、変数 a から変数 b に矢印が引かれているのは、b が a に依存している、みたいな……」
「依存とは？」
「待って待って、依存、うーん、条件付き分布 p(b|a) がある？」
「……今のところはそれでいいか。後でもう一度確認しよう」
「つまりなんか間違えてるんですね……」
「間違いというほどではないけどね。とりあえずその定義で考えても、a→b なら p(b|a) が、b→a なら p(a|b) があるということだから、矢印の向きに意味はある。どう？」
「あれれ？　いや、そうじゃなくて、生成モデルって言うんですか、つまり同時分布 p(a, b) があるときに、これを乗法定理で p(a)p(b|a) と分解するとグラフは a→b になるけど、p(b)p(a|b) と分解すると b→a になるわけですよね。まあ p(a|b) と p(b|a) のどちらを考えるのが自然かとかはあるんでしょうけど、それに目をつむれば分解の選択は恣意的なので、矢印の向きに本質的な意味はないのかなあと」
「そうだね。その例なら」
「やっぱり！　……って、その例なら？」
「変数2個で例を作るのは無理があるから3個にして、設定も具体的にしようか。変数 a と b を『サイコロを振って出る目』、c を『その和の偶奇』とする。いい？」
「鉄板例ですね……んーでもそれってグラフは元からこの形に決まってるんじゃないですか？」

「自分から『分解の選択は恣意的』と言っておきながら、都合でそれを忘れるのはさすがにどうかと」
「ぐぬぬ」
「……というわけで、同時分布の分解を『恣意的に選んで』このモデルを解釈してみて」
「はーい……ではまず p(a, b, c)=p(a)p(b|a)p(c|a, b) と分解してみます。グラフはこんな形。

んで……それから……？」
「今、各変数は具体的にわかっている」
「え？　ああ、p(b|a) とかが書けるってことですね。えと、a はサイコロの目そのままだから p(a) は毎度おなじみ 1/6 ずつ。p(b|a) は……もう一つのサイコロの目に影響受けないから、やっぱり 1/6 ずつ。p(c|a, b) は…… a と b が与えられたら c は確定だから、a+b が偶数なら p(c=偶|a, b)=1、以下省略」
「気づくことは？」
「気づく、というか、最初からわかってたことですが a と b は独立ですよね」
「独立とは？」
「独立…… b は a によらない……？」
「今何を求めていたの？」
「えっ？　あ、そうか、p(b|a) は a によらない、p(b|a)=p(b) だ」
「そう。独立の定義そのまま。『b は a によらない』とかあいまいなことを言う必要はない」
「ぐはっ……もうあたしの HP はゼロです……」
「……このときグラフィカルモデルは？」
「ん〜、p(b|a)=p(b) だから a から b への矢印が消えてこの形に。なるほどな〜。

要領わかってきましたよ。次は p(a, b, c)=p(c)p(b|c)p(a|b, c) でやってみますね。まずグラフ描いといて。

c は偶数奇数だから p(c=偶)=p(c=奇)=1/2 ですね。p(b|c) は……ビミョウにややこしい……」
「表を書いたら？」
「そっか。a と b とで表にすると、

a＼b	1	2	3	4	5	6
1	偶	奇	偶	奇	偶	奇
2	奇	偶	奇	偶	奇	偶
3	偶	奇	偶	奇	偶	奇
4	奇	偶	奇	偶	奇	偶
5	偶	奇	偶	奇	偶	奇
6	奇	偶	奇	偶	奇	偶

ふむふむ、これなら簡単ですね。c が偶数の時、b のそれぞれの確率は

p(b=1|c=偶) = 3/18 = 1/6
p(b=2|c=偶) = 3/18 = 1/6

あ、あれ？　変だな〜」
「どうかした？」
「全部 1/6 になります」
「そうだね」
「p(b|c)=p(b) になっちゃいます」
「なるね」
「ダメでしょう!?」
「なぜ？」
「だって、c から b への矢印が消えて

ってグラフになっちゃいます。同じモデルが2種類の全く違うグラフを持つとか、どう見てもまずいじゃないですか」
「つまりどこかが間違っている？」
「それがさっきから見直しているんですが、どこも間違えてなくて……」
「間違ってないからねえ」
「えーっ！　いやでも。え〜……。そんなんアリなんですか？」
「そもそも、モデルの情報をわずかな丸と線でどれくらい表せると思う？」
「そりゃまあ全部ってことはないでしょうけど」
「一部しか表せないなら、『どの一部を表すか』という選択が生じる。つまり、グラフィカルモデルによる表現は一意ではない」
「はぁ〜なるほど。分解に恣意性があるようにグラフにも恣意性があるわけですね」
「普通は構成から自然に決まるグラフしか見ないから意識しないけれどね。これがまず一つわかること」
「もう一つあるんですか？」
「最初のあなたの疑問の答え。今得た2つのグラフィカルモデルを見比べてみて」
「最初の疑問って矢印の向きの話ですよね。でも、向き以前にグラフの形というか、切れているところがぜんぜん違いますよ」
「いいね。それが答え」
「む〜〜〜。禅問答じゃないんですから、もうちょっとわかりやすくプリーズ」
「では、切る前と切った後を並べてみよう」

「これ見せられてもやっぱりわかりませんよう。これ見てわかることなんて、向きが違うと切れるところが違うってことくらい……あ」
「気づいた？」
「左は a→b, b→c, a→c 『だから』 a→b が切れる。右は c→b, b→a, c→a 『だから』 c→b が切れる。向きが逆なら切れる場所が変わる！」
「そう。つまり、ベイジアンネットワークがこの構造を採るためには、各辺はこの向きである必要がある」
「この向きがあって初めてこの形のグラフになるんですね……ん？　ちょっと待ってください」
「何？」
「こっちの後から描いた方のグラフでも、a と b って独立じゃあないですか。だからもともとそこは切れるんじゃないですか？」
「そして先ほど b と c が独立とわかったのと同様に a と c も独立と示せるので、すべての辺が切れる、と」
「あれれ？」
「ちょうどいい。矢印の定義の話に戻ろう」
「あたしのが間違えてたやつですね」
「だから間違っていたわけではないよ。矢印が条件付き分布に対応するのはいい。ただ、親の変数は複数を許す。つまり、p(x|y_1, …, y_n) があることをグラフを用いて次のように表す」

「そうか！　指している先が同じ矢印は全部セットなんですね」
「だから a と b が独立であっても、p(a|b, c) で引いた辺の片方を勝手に消すことは出来ない」
「a と b、a と c が両方とも独立でも、 a と b, c の組が独立とは限らないですもんねー」
「どう？　疑問はかなり解消されたと思うけど」
「はい、ストンときました！　でも考えたんですけど、『この構造のためにこの向き』ということは、もちろん一つ一つの矢印の向きも大事なんでしょうけど、『どこが切れているか』のほうがもっと重要なのかなあ、なんて思えてきました」
「そうだね。PRML にもそう書いてあるし」
「うんうん……えっ!?　書いてあるんですか？」
「確かこのあたり……あった、下巻 73 ページ」
「『グラフはリンクが存在しないことをもって分布のクラスの性質に関する情報を表現する』うわーほんとだ、まんま書いてある……なんかもっと目立つように書いて、大事なんですよ！　って主張してくれないと見逃しちゃいますよう」
「一応太字で記されてはいるのだけれどね。四倍角で書いてあってもいいくらいではあるだろうね」
「よんばいかく？」
「知らないの？　文芸部なのに。ほら、ワープロで文字を縦横二倍の大きさに……」
「わーぷろ？」
「……文章を打つときに。会誌用に原稿を清書とか」
「ああ、ワードのことですか」
「ワードはパソコンでしょう。そうじゃあなくて、専用の。インクリボンで印刷する……」
「いんくりぼん？」
「……」
「……そういえば先輩の原稿って、文字がギザギザしてムラがあって、なんかちょっと昔風で味があるなあと思ってたんですが、あれはわざとではなく？」
「おとうさ……父が貸してくれて……最新の機械だって……まだ普及していないからインクリボンも手に入りにくいので、大事に使うようにと……」
「ひょっとしてですけど、その機械、印刷うるさかったりしません？」
「……する……」
「あの〜、大っ変申し上げにくいんですが、それ、たぶん絶対からかわれてますね。きっと」
「……また、だまされ……いや、なんでもない」
「……許せませんね」
「いや、別にあなたが怒ることじゃあ……」
「PRML では鋭い切れ味を見せる才媛ながら、普段は迂闊なお人好しキャラとか、ずるいです！」
「え……えええー！　ぼ、僕！？」
「あたしの立ち位置のことも考えてください！！　読者が許してもあたしが許しません！！！」
「そんなこと言われても……いえ、ごめんなさい……」
「今度から気をつけてくださいね！」

（おしまい）

---

新春お年玉駄文……のはずが、仕上げに時間がかかりすぎて危うく２月になるところだった。ふう。

今までの分：

• 独断と偏見によるノンパラ入門 - 木曜不足
• PRMLガール 〜 文芸部のマネージャーが「パターン認識と機械学習」を読んだら 〜 - 木曜不足

「機械学習 はじめよう」最終回

gihyo.jp で続けていた連載「機械学習 はじめよう」がとうとう最終回を迎えることができた。

• 第21回（最終回） 機械学習 はじめよう：機械学習 はじめよう｜gihyo.jp … 技術評論社

初回の 2010年 7月から 3年以上にわたる連載がなんとか完結できたのも、読んで下さった方々、担当編集者の高橋さんのおかげ。大感謝。
後半は2〜4ヶ月に1回のスローペースで全21回、年内にすっきり終わりたいね、12月25日公開ならクリスマスプレゼントって言えるねｗ　と、難産ながらそれなりにまとまったかと。

言いたいことはだいたい連載の中に置いてきたので、付け加える事はそんなにないのだけど。
PRML で機械学習を勉強していたときに、機械学習の国ではいろんなことが「あたりまえ」になっていて、特に説明もなしに使われるということに非常に戸惑った（まあどんな分野でも多かれ少なかれあることだけど）。
そういう「あたりまえ」は中の人になってしまうと「あたりまえ」すぎてわざわざ説明しないといけないということに気づきにくくなっていくので、その気持ちが新しいうちに後に続く人たちのためになんとか残しておきたい、という意識で書き続けてきた。

「機械学習」でぐぐったら Wikipedia の次にこの連載が並ぶので、今では初学者の入り口の一つとなっている、と言ってしまったら自意識過剰だろうか。
そういうものを書いて、終わらせることができたのはかけねなしに嬉しい。

内容への指摘・お叱りは今後も随時受け付けているので、何かお気づきの事があればブログのコメントや twitter (@shuyo) などへ。