elsur.jpn.org >

« 2018年4月 | メイン | 2018年6月 »

2018年5月28日 (月)

Dotson, J.P., Howell, J.R., Brazell, J.D., Otter, T, Lenk, P., MacEachern, S., Allenby, G.M. (2018) A Probit Model with Structured Covariance for Similarity Effects and Source of Volume Calculations. Journal of Marketing Research, 55(1), 35-47.
 新製品のSOV(Source of volume)測定について考えていて、なにかの足しになるかと思って読んでみた。先日出たばかりの論文である。

 数ページ読んで、コンジョイント分析に特化した分析手法提案であることに気づき、ちょっとやる気が萎んだ。前半の理論編を読み終えるあたりでようやく、ベイズモデリング界の有名人・Allenby兄貴が共著者であることに気づき、急速にやる気が衰えた。でもまあ、いちおう最後までめくりました。

 いわく。
 SOV(source of volume)、つまり需要変化の源について知ることは大事だ。わが社の新製品が自社の現行製品とカニバリゼーションを起こさないかしらん、とか。
 製品間の代替性のパターンを完全にモデル化するとなると、製品数を$J$個として$J^2-J$個の交差効果が必要になる。コンジョイント分析なんかだと架空の製品がたくさんあるわけで、とてもじゃないがモデル化できない。
 ロジットモデルの世界では、代替性のパターンを柔軟にモデル化するよりも倹約性のほうが重視される。というわけで、ふつうIIAが仮定される。IIA仮定をおいたロジットモデルでは、個人レベルでのSOVは計算できないけど、消費者には異質性があるから、累積レベルではSOVが計算できるんだよ、という理屈である。我々はこういう路線は採らない。ほんとうの問題は、ロジットモデルは選択肢間での誤差の相互依存性を説明していないという点にある。
 そこで、多項プロビット(MNP)モデルに、構造化した共分散行列を導入する。この共分散行列は選択肢間の知覚距離に基づいて決める。知覚的に類似した選択肢の間では選択確率が相関すると考えるわけである。
 多項モデルに誤差の相関を導入する提案はすでにある。Hausman&Wise(1978)は正規誤差を、McFadden(1978)は極値誤差を考えた。しかしそういうモデルでは、誤差項は提示順序とかブランドとかと結びつけて考えられており、それぞれの選択肢の具体的な属性・水準と関連付けられることはなかった。[←そうなんですか?]
 
 モデル。
 マーケティングにおける選択モデルは、ふつう誤差項がiidだと仮定する。そのため、どんな選択肢でも誤差項の実現値によっては勝ちうるということになり、弱い製品への需要が過大評価され、予測シェアは0にならない。
 このようにモデルの倹約性に価値が置かれていたのは、マーケティング・データはふつう、対象者数は大きいけど対象者あたりの観察数が小さいからである。誤差項に相関を許した先行研究もあるけれど、異質性の表現だけが目的だったり、nested logitモデルによって相関構造を制約したりしていた。
 さて、我々のモデルは次の通り。
 製品の属性のうち観察可能なのを$x_o$、観察不能なのを$x_u$とする。世帯特性のうち観察可能なのを$h_o$, 観察不能なのを$h_u$とする。$x= [x_o, x_u], h=[h_o, h_u]$とする。
 効用関数を
 $u(x, h) = V(x_o, h_o) + \epsilon(x, h)$
と仮定する(誤差項は加算されるのだと仮定する)。意思決定者は選択を通じて効用を最大化すると仮定する。
 決定論的効用について、選択課題を$c$、選択肢を$k$として
 $V_h(x_{hck}) = x'_{hck} \beta_h$
と考える。
 さらに、誤差項も$x, h$の関数と考える。誤差$\epsilon_{hck}$は平均ゼロの正規分布に従うと考え、その共分散行列を世帯x選択課題別に$\Sigma_{hc}$とする。これはサイズ$K$の行列で(つまり行と列は選択肢)、対角要素は1、非対角要素は$\sigma_{ij}$である。なお、対角要素を1にしないというのも試したが、モデルの性能は上がらなかった。
 当然ながら、世帯x選択課題別に共分散行列を推定できるだけのデータはない。そこで$\sigma_{ij}$についてこうモデル化する。
 $\sigma_{ij} = \exp(-d_{ij}/\theta_h)$
$d_{ij}$はスカラーで、選択肢$i, j$の距離を表す。$\theta_h$は距離をリスケールするパラメータで、未知の世帯特性によって決まると考える。ゼロに近いときにIIAに近づく。

 では、選択肢$i$と$j$のあいだの距離$d_{ij}$をどうやって求めるか。4案考えました。

 なお、選好と知覚とは別だろうと考え、距離の式の中に出てくる$\beta$を別の自由パラメータに置き換えるというのも試したが、モデルの性能は上がらなかった。

 こうしてみると、従来の選択モデルと比べて新たに導入されたパラメータは、結局$\theta_h$だけである。
 [以下、$\theta_h$の識別条件とか、もうちょっとだけモデルを拡張して負の共分散を許しました、とかといった話題。パス]

 理論編の最後になりましたが、先行研究紹介。[ここは勉強のために詳しくメモ]
 iidでない誤差項をいれることでIIA仮定を避けるという文献は山ほどあるが、大きく3つの次元を考えることができる。

本モデルは、3つの次元すべてにあてはまる。

 実験やりまーす。
 コーヒー・クリーミング・パウダーについてのコンジョイント実験。属性は、ブランド(4水準)、乳製品か(2)、無脂肪か(2)、香りつきか(2)、サイズ(2)、価格(5)。951名。14試行。完全プロファイルで4選択肢並べてひとつ選ばせる。1試行をホールドアウト。
 すべてのモデルパラメータをランダム効果とする。対象者を$h$として
 $\beta_h \sim N(\bar{\beta}, V_\beta)$
でもって、次のモデルを比較する:

 結果。インサンプルの適合度でも[Gelfand & Dey(1994, JRSS)の推定量のlogmarginal densityというのを使っている。そんなん知らんがな]、ホールドアウトの予測でも、距離3による提案モデルが優れていた。パラメータ推定値を比べると...[略]。誤差項の相関をこまかくみていくと...[略]。
 別のデータセットでやっても、やっぱし提案モデル(距離3)が優れてました。
 
 SOVの推定について。
 新規参入のSOVをどう推定するか。まずは多項ロジットモデルの場合について考えてみる。IIAが仮定されているから、SOVのレートは参入前の選択確率と同じである。既存品を$i$, 新規品を$k$, すべての製品を$j$とする。$i$の決定論的効用を$V_i = x'_i \beta$とする。参入前は
 $Pr(i)_{Before} = \exp(V_i) / \sum_j \exp(V_j)$
参入後は
 $Pr(i)_{After} = \exp(V_i) / ( \sum_j \exp(V_j) + \exp(V_k) )$
 $Pr(k)_{After} = \exp(V_k) / ( \sum_j \exp(V_j) + \exp(V_k) )$
製品$k$のSOVにおける$i$の割合は
 $\%SOV = (Pr(i)_{Before}-Pr(i)_{After})/Pr(k)_{After}$
これは結局$Pr(i)_{Before}$となることが示せる。代替性を調べているはずが、実際には選好を調べているわけである。たとえば低価格のために選択確率が高くなっている選択肢があったら、そこからのSOVも高くなるわけで、これはおかしい。
 ロジットモデルのもつこの制約的な性質は、異質な消費者集合について累積することで緩和されるだろうか。
 クリームパウダーの実験に基づき選択肢シミュレータをつくった。パラメータとしてランダム効果の事後分布の平均$\bar{\beta}$を使う版と、事後分布$\pi(\{\beta|Data\})$を使う版をつくった。なお後者は$\{\beta_h \}$を使うのと似た結果になる。このことは実務家ならみんな知っている。Huber, Orme, Miller(1999, Sawtooth Software Conf.)をみよ。
 まず3選択肢でシェアを求める。選択肢2,3は属性がひとつしか違わないものにする。次に選択肢1と属性がひとつしか違わない選択肢4を追加し、再度シェアを求める。これを繰り返す。
 選択肢4のSOVは1で大きくなるはずである。しかし多項ロジットモデルの場合はそうならない。これは$\{\beta_h \}$版でも対して改善しない。つまり、個人レベルでIIAを仮定しちゃったことによって生じる問題は、異質性を盛り込んでも取り返せない。問題の深刻さは異質性の程度によって決まるわけだけど、それは事前にはわからない。
 [さらに選択肢1の全体効用を操作した実験を行っている。主旨は同様なので省略]
 [別のデータでも試しているけど、もう疲れたのでパス。やりすぎでしょう...]

 というわけで、誤差相関が選択肢の価値と構造的に関連している選択モデルを御提案しました。相関構造は製品属性から勝手に決まります。
 距離3が優れていたということは、選択肢間の相関は、どの属性のせいで効用が高い・低いかとは無関係に、補償的に決まる模様です。マーケティングや心理学の研究では、製品の類似性ってのは製品属性で決まると考えるのがふつうで、効用で決まるとは考えない。いっぽう本研究の結果は、製品が異なっているときだけでなく、製品の効用が異なっているときにも、対象者はそれを弁別しやすい、ということを示唆している。
 実務家のみなさん、ロジットモデルをよく使っておられますが、製品間の代替性を調べるときにはIIA仮定の影響が強すぎて怖い、ということを知っていただきたい。
 今後の課題は、異なる代替関係をもたらすような選択理論を構築すること。Tversky&Simonsonの文脈依存的選好の研究とかがこれに近いが、そういうのは遷移律からの逸脱に焦点を当てているのに対して、代替性が製品によって違うというのは遷移律の違反ではない、という点が異なる。
 云々。

 ... やれやれ、いちおうめくったぞ。疲れたぞ。
 私の仕事との関係で言うと、「コンジョイント分析におけるIIA仮定は非現実的だけど個人効用をHB推定してるから大丈夫」というのは半ば幻想である、という点が、ちゃんとしたシミュレーションでわかったというのが収穫であった。そうなんすかー。参っちゃいますねー。
 いっぽう、あんまり真面目に読んでないので読み落としてるのかもしれないけど、提案手法じゃなくてHBプロビット($\Sigma_h$)でも別にいいんじゃね?とも思った。モデルのパフォーマンスは遜色ないみたいだし。SOV算出時のシミュレーションは示されてないのだが、実は提案手法と対して変わらんのではないかしらん。えーっと、HBプロビット($\Sigma_h$)なら、Allenby先生のbayesmパッケージでできちゃいません?

 実務的問題とはちょっと離れてしまうけど、距離3が優れているのはなぜか、って話が興味深かった。大きく言うと、実験者は選択肢を属性の束として捉えているけど、被験者はその属性の束から潜在属性(ここでは全体効用)を生成したうえでやおら課題を遂行しているかも、ということであろう。Barsalouの「アドホック・カテゴリ」みたいな話ですね。ああ栄光の80年代認知心理学。

 ベイズ・モデリング界のスターであるAllenby兄貴の論文には、アクション映画界のスターであるドニー・イェン兄貴のアクションが似合うので(前者は読者に、後者は共演者に、無理難題を押し付けるところが似ている)、今回も動画を貼っておくことにする。ピーター・チャン監督「武侠」(2011, 日本公開タイトル「捜査官X」)より。気弱な農民のはずのドニーさんが、クララ・ウェイ姉さんの横暴を前にし、ついにその本性を明らかにするという名場面である。
 屋根の上を疾走するシーンも忘れがたいですが、この映画の最大の見どころは、リアリティがありそうでなさそうな、さまざまな謎のキメポーズだと思う次第である。牛小屋の対決でクララ姉さんが刃物振り回しているときの、ドニーさんのあのポーズはなんなの? 両手のこぶしを顔の前に10cmくらい縦に離して構える奴。全然意味がわかんないんだけど、やだもう、超かっこよくないですか? 私もいまディスプレイの前で試してます。


論文:データ解析(2018-) - 読了:Dotson, et al. (2018) 選択肢の全体効用が似ているときに選択肢間の誤差相関を高くする階層プロビット選択モデルを使ったコンジョイント分析でSOVを推定しよう (feat. 武侠)

2018年5月25日 (金)

 こないだ学会(消費者行動研究学会)のワークショップを拝聴していて気が付いたんだけど、反応時間についての私の知識は古い。全然勉強してないんだから当然である。
 まあいまさら私が勉強してもしょうがないのかもしれないが、いずれ反応時間を使ってやってみたいこともあるし、少しでも勉強しておくか... と思って、適当そうな論文を読んでみた。当座の仕事とは関係ないわけで、まあ半分趣味みたいなもんである。

Krajbich, I., Bartling, B., Hare, T., Fehr, E. (2015) Rethinking fast and slow based on a critique of reaction-time reverse inference. Nature Communications, 6, 7455.
 第一著者はオハイオ州立大の心理学者。ワークショップで某先生が口走っていた名前をカタカナでメモし、しばしの試行錯誤の末に突き止めた。反応時間の研究で有名な研究者らしいのに、恥ずかしながら初耳だ...と思ったが、調べたところ2011年に学位取得とのこと。そりゃあ知らないわ。

 いわく。
 意思決定プロセスの研究は、二重過程説の観点に立ち[直観的・自動的過程と熟慮的過程、Type IとType IIってやつですね]、ある行動はどっちの過程の結果か、と問う。それを区別する手段のひとつが反応時間(RT)の検討である。直観的過程は熟慮的過程より速かろうという理屈である。
 というわけで、意思決定の研究ではRTに基づく主張が溢れている。しかし、ある自動的過程が熟慮的過程よりも速いだろうと予測するということと、ある選択についてこれは速いから自動的だと分類するということとは、別の問題である。
 RTに寄与する認知過程にはいろいろあることがわかっている。RTに基づく推論はそれらの認知過程についても説明しないといけない。たとえば、RTは弁別性と関連する。これは記憶・知覚から経済的選択にいたる広い範囲で示されている。90年代、視覚探索過程が系列的か並列的かを区別するためにRTが使えるかという点についての論争があったが、そのとき批判者は、RTの主要な規定因がむしろ弁別性であることを示したものである。また経済的選択の分野では、選択肢が類似している決定課題は遅いことが示されている。
 つまりこういうことだ。選択は実は単一の熟慮的過程によってなされていて、RTの変動は複数種類の過程の競合のせいで生じているのではなく、むしろ選択肢の類似性知覚のせいで生じている... こういう可能性について、ちゃんと検討しないといけない。
 本論文は社会的選好と対人的選択の文脈において、上記の点を詳細に示す。選択肢をうまいこと作り込めばお望みのRTが得られる、その様子をとくとご覧いただきたい。RTは二重過程理論の証拠にはならないのである。

 まず、RT逆向き推論問題(RT reverse-inference problem)について紹介しよう。
 「人が直観的に好むのはエールかバーボンか」を調べる実験について考える。エールの選択肢セットからひとつ、バーボンの選択肢セットからひとつを提示し、どちらかを選ばせる。これを繰り返し、選択とその反応時間を調べる。
 当然ながら、選択の結果は選択肢セット次第である。いま、実験者1がつくった選択肢セットではエールの勝率が高く、実験者2がつくった選択肢セットではバーボンの勝率が高かったとしよう。
 選択課題のRTはどうなるか。RTは選択肢間の差異の知覚に依存する。いま選択肢AとBの間で選択するとき、横軸に選好の差(A-B)、縦軸にRTをとると、0を中心にした山形になると期待できる。実験者1の場合、選択肢のなかに良いバーボンが少なかった、つまり左裾に位置する選択肢ペアが少なかったわけだ。だから、Bが選ばれる試行(真中から左側)のRTの平均は、Aが選ばれる試行(右側)のRTの平均より高くなる。こうして、実験者1は「人はエールを直観的に選ぶ」とうっかり結論してしまうことになる。実験者2はその逆になる。
 [こういう順序で説明されると、いやあそんな勘違いはしませんよ...とつい思ってしまうけど、実際にはこういうの結構多いですよね。「XXの選択率は高くRTは短かった、つまりXXは消費者に直観的に選ばれているのだ」的なの]

 ここからは実際の研究に即して説明しよう。3つ例を挙げる。

 その1、社会的選好。
 独裁者ゲームというのがある。独裁者役の被験者が、独裁者と受け手に金をどう配分するかについて、自己中心的選択肢と向社会的選択肢のどちらかを選ぶ。
 実験やりました。独裁者役は25名、ひとり70試行。
 まずは選択タイプとRTの関係をみてみよう。自己中心的選択のほうがRTが短かかった。ってことは、自己中心的意思決定は直観的で、向社会的意思決定は熟慮的だ...と思っちゃいますよね?
 こんどは各被験者について、自己中心的選択肢の選択率、自己中心的選択のRT、向社会的選択のRTを調べる。すると、自己中心的な選択をしやすい被験者は自己中心的決定が速く、向社会的な選択をしやすい被験者は向社会的選択が速かった。また、自己中心的選択肢の選択率が5割近辺の人は、RTの差も0に近かった。
 さらに、全体を通して、自己中心的選択肢の選択率が高かった。これは選択肢の作り方の問題であって、各試行において「向社会的選択肢を選ぶことで独裁者が損する金額」がもっと低くなるようにしておけば、向社会的選択肢の選択率はもう少し上がったはずである。
 というわけで、選択肢の作り方のせいで自己中心的選択肢の選択率が高く、よって自己中心的選択をしやすい被験者が多く、よって自己中心的選択のほうがRTが短くなった、と解釈できる。
 分析してみよう。自己中心的選択肢$i$と向社会的選択肢$j$のペイオフを$x_i, x_j$とする。選択肢ペアの効用関数を次のようにあらわす。
 $U(x_i, x_j) = (1-\beta r - \alpha s)x_i + (\beta r + \alpha s) x_j$
$r, s$は独裁者のペイオフが受け手のペイオフより高い/低いことを示すダミー。$\beta,\alpha$は個人パラメータ。
 それぞれの試行について、推定した効用を横軸、RTを縦軸にとると、どちらを選んだ場合でもRTの分布は変わらない。混合効果回帰でも確認できる。このように、「自己中心的選択は速い」というのは実験パラメータのアーティファクトである。

 その2、時間選好。
 被験者41名が選択課題216試行を行った。選択肢は「25ドルをいまもらう」「xドルをt日後にもらう」の2択(xは25より大)。全試行の約半分(53%)で即時選択肢(25ドル)が選ばれた。RTに有意差はなかった。
 以下では、一部の選択肢ペアだけを取り出し、即時選択肢のほうが魅力的であるデータを作って分析する。するとRTは即時選択のほうで速くなる。即時ペイオフの選択は直観的だと結論しかねない。
 遅延ペイオフについて次の時間割引関数を考える。
 $U(x, t) = x / (1+ kt)$
$k$は個人パラメータ。このモデルで推定した効用差を横軸、RTを縦軸にとると、どちらを選んだ場合でもRTの分布は変わらない。混合効果回帰でも確認できる。

 その3、Rand, Greene, & Nowak (2012 Nature)。[←読んでてだんだん眠くなってきてうとうとしはじめてたんですが、急にシャキッとしました! 性格の悪い論文大好き!!]
 彼らは人間の利他的協同が二重過程に支配されていると主張している。彼らによれば、直感的システムは協同を支持し、熟慮的システムは自己中心性を支持する。
 彼らが行った公共財ゲームについてみてみよう。被験者を4人ずつのグループにする。各被験者に40単位の金銭を与え、いくらを手元に残し、いくらを公共財のために寄付するかを決めさせる。実験者は寄付金を倍にしてメンバーに均等に配る(1単位の寄付につき0.5単位が戻ってくる)。
 彼らの結果では、は多く寄付するという選択でRTが短かった。なお、彼らはほかにもいろんな手続きをやっているんだけど、ここでは脇においておく。
 この実験を次のように改変した場合について考えよう。
 A1. 1単位の寄付につき0.9単位が戻ってくる(全額寄付すると36単位戻ってくる)。
 A2. 1単位の寄付につき0.3単位が戻ってくる(全額寄付すると12単位戻ってくる)。
 A1において向社会的な被験者はRTが短くなり、A2において自己中心的な被験者はRTが短くなるだろう。
 再現実験やりました。{オリジナルと同じ設定, A1, A2}のそれぞれについて被験者175名。二重過程説の説明によれば、どのバージョンでも、寄付金額とRTは負の相関を持つはずである。
 結果。A1では正の相関、オリジナルでは0に近く、A2では負の相関になった。オリジナル条件での元論文とのずれは被験者の寄付額のちがいのせいだろう。

 このように、RTの解釈に際しては選好の強度を考慮しなければならない。RTによる二重過程理論の支持には落とし穴がある。証拠蓄積モデルやベイジアン更新メカニズムのような単一過程モデルだって、選択バイアスとRTを同時にうまく説明できる。
 要約しよう。ある行動が二重過程理論でいう直観的コンポーネントによって支配されているかどうかを推論する際に、決定タイミングを根拠とするのは問題がある。選択タイプによるRTのちがいは、背後にある過程のちがいのせいではなくて、選好強度や弁別性のちがいのせいかもしれない。意思決定のモデルを比較しようとするみなさん、反応時間データはもっと注意して扱いなさい。

 ... いやー、面白かった!
 メモを読み返しながらあれこれ考えるわけだけど、この論文はあれだろうか、反応時間にはいろんな要因があるから、そこからあるモデル(例, 二重過程モデル)に基づく示唆(例, 向社会的行動は直観的だ)を引き出すのはそもそも難しい、という悲観的な教訓として捉えるべきであろうか。それとも、ちょうど著者らが個人ごとの主観効用の差を推定して回帰に放り込んでいるように、うまい実験計画なりうまい分析計画なりを組む必要があるのだ、という楽観的な教訓として捉えるべきだろうか。
 早い話、もしRand, Greene, & Nowakがいろんなインセンティブ設計を通じて寄付金額とRTとの負の相関を示していたとしたら、選好強度説は排除され、二重過程説に基づく「向社会的行動は直観的だ」という示唆が得られていたのではないかしらん? それとも、排除しきれない第三の説がありえたかしらん?

 これは私のような立場の素人にとってもちょっと切実な話で...
 マーケティングのための消費者調査で回答の反応時間を分析するとしたら、それは判断の背後にある認知過程のタイプを知るためではないだろう。おそらく、(1)回答態度が反応時間と関連するというモデル(「長すぎる/短すぎる回答は真剣でない」)に基づき回答態度について推論しようとするか、(2)態度の強度が反応時間と関連するというモデル(「選好判断が速いブランドは選好強度が強い」)に基づき態度の強度について推論しようとするか、のいずれかであろう。でもよく考えてみると、反応時間に影響するであろう要因は他にいっぱいあり(「自動的認知過程に支配されていると速い」「知覚的流暢性が高い刺激に対する判断は速い」)、いつだって代替的説明が可能なのである。
 そういうのを潰すうまい手続きはありうるのか? それとも、反応時間というのは結局はよくわからんもんだと割り切り、あるドメインにおけるある事柄に対して反応時間が予測的妥当性を示すという知見を、そのドメインのなかで地道に積み上げていくべきなのか? ...という問題である。

論文:心理 - 読了:Krajbich et al.(2015) 直観的判断は速く熟慮的判断は遅いと考えられているが、では速い判断は直観的で遅い判断は熟慮的だといえるか

水野将樹 (2003) 心理学研究における「信頼」概念についての展望. 東京大学教育学研究科紀要, 43, 185-195.

 仕事の都合で目を通した。えーっと、信頼ってのは測定の信頼性とか推定の信頼区間とかの話じゃなくて、日常用語でいうところの信頼、trustのことね。
 修論のまとめだと思うんだけど、とにかくレビューというのは助かります。

 いくつかメモ:

...職場における信頼関係についての研究っていっぱいありそうな気がするけど、ああいうのはみんなBに分類されるのかしらん。そうしてみるとBもずいぶん幅がありそうだなあ。

論文:心理 - 読了:水野(2003) 心理学における「信頼」概念レビュー

2018年5月14日 (月)

 仕事の都合で不意に対応分析を使ったりすると、いつも途中で混乱しちゃうので、思い余ってメモをとった。
 でもあれなんだよな、メモを取れば覚えるかというと、これがそうでもないんだよな... 思えば、こういうのはじめて勉強したのは20代前半であった。それから時代は流れ、幕府は倒れ、江戸は東京となり、日本は急速な近代化を遂げたが(←ちょっと記憶が混乱してます)、私は全然進歩してないような気がする。

 以下、Greenacre(1984) "Theory and Applications of Correspondence Analysis", Chap.4 の前半のメモ。この本だと、固有値分解ではなくて一般化SVDでスーッと説明しちゃうのである。

準備
 $I \times J$行列を$\mathbf{N}=\{n_{ij}\}$とする。値は非負、行和と列和は非ゼロとする。そうでなくてもかまわないけれど、$\mathbf{N}$は二元クロス表だと考えるとわかりやすい。
 $\mathbf{N}$で総和で割って$\mathbf{P}=\{p_{ij}\}$とし、これを「対応行列」と呼ぶ。$\mathbf{P}$の行和ベクトルと列和ベクトルをそれぞれ$\mathbf{r}=\{r_{i}\}, \mathbf{c}=\{c_{j}\}$とし、これを持たせた対角行列を$\mathbf{D}_r, \mathbf{D}_c$とする。
 
行プロファイルと列プロファイル
 $\mathbf{P}$のある行$\mathbf{r}_i$をその行の和で割ったもの$\tilde{\mathbf{r}}_i$を「行プロファイル」、ある列$\mathbf{c}_j$をその列の和で割ったもの$\tilde{\mathbf{c}}_j$を「列プロファイル」と呼ぶ。プロファイル行列は
 $\mathbf{R} \equiv \mathbf{D}^{-1}_r \mathbf{P}$
 $\mathbf{C} \equiv \mathbf{D}^{-1}_c \mathbf{P}^T$

行プロファイル空間と列プロファイル空間
 行プロファイル$\tilde{\mathbf{r}}_i$を、$J$次元空間におけるその行の座標とみる。その行の行和、つまり行和ベクトル$\mathbf{r}$の$i$番目の要素$r_i$を「質量」と呼ぶことにする。
 二点間の距離を、重み$\mathbf{D}^{-1}_c$による加重ユークリッド距離(=カイ二乗距離)とする。つまり、次元$j$の重みは$1/c_j$となるわけである。
 そもそも行$i$の行プロファイルの$j$番目の要素は$p_{ij}/r_i$であった。すべての行プロファイルについてその重心を求めると、$\sum_i r_i (p_{ij}/r_i) / \sum_i r_i = \sum_i p_{ij}= c_j$であるから、つまり行プロファイルの重心は列和ベクトル$\mathbf{c}$となる。
 以上は行と列をひっくり返しても成立する。つまり$\tilde{\mathbf{c}}_j$を$I$次元空間における座標とみて、二点間の距離を重み$\mathbf{D}^{-1}_r$による加重ユークリッド距離としたとき、その重心は行和ベクトル$\mathbf{r}$である。

全慣性
 ちょっと話が逸れるけど... 上の行プロファイル空間なり列プロファイル空間について、各点の重心からの二乗距離の、質量を重みにした加重和を「全慣性」と呼ぶ。
 行$i$の重心からの二乗距離は、もし次元に重みがなかりせば
 $(\tilde{\mathbf{r}}_i - \mathbf{c})^T (\tilde{\mathbf{r}}_i - \mathbf{c})$
だが、軸に重み$\mathbf{D}^{-1}_c$がつくので
 $(\tilde{\mathbf{r}}_i - \mathbf{c})^T \mathbf{D}^{-1}_c (\tilde{\mathbf{r}}_i - c)$
で、この子の質量が$r_i$である。質量で重みづけて合計すると
 $in(I) = \sum_i r_i (\tilde{\mathbf{r}}_i - \mathbf{c})^T \mathbf{D}^{-1}_c (\tilde{\mathbf{r}}_i - \mathbf{c})$
同様に
 $in(J) = \sum_j c_j (\tilde{\mathbf{c}}_j - \mathbf{r})^T \mathbf{D}^{-1}_r (\tilde{\mathbf{c}}_j - \mathbf{r})$
これが、行プロファイルの空間の全慣性$in(I)$, 列プロファイルの空間の全慣性$in(J)$である。
 導出は省略するけど、結局
 $in(I)= in(J) = \sum_i \sum_j \{(p_{ij} - r_i c_j)^2/r_i c_j\}$
であることが示せる。この式をよーく見ると、$\mathbf{N}$を二元クロス表とみて、その独立性検定をしたときのカイ二乗値を、総和で割った値になっている。

一般化SVDの登場
 さあ、ここからが本番だ!

 上の行プロファイル空間・列プロファイル空間から、少数次元の下位空間を取り出したい。その際、オリジナルの空間からの各点の二乗距離の、質量を重みにした加重和を最小化したい。
 こういうのは、ご存じ一般化SVD
 $\mathbf{A} = \mathbf{N} \mathbf{D}_\mu \mathbf{M}^T$ ただし $\mathbf{N}^T \mathbf{\Omega} \mathbf{N} = \mathbf{M}^T \mathbf{\Phi} \mathbf{M} = \mathbf{I}$
の得意技である。
 主はいわれた。あなたがたはこれまで、$\mathbf{A}$の行を点、列を次元とみなしていた。はっきり言っておくが、$\mathbf{N}$の左$K^*$列, $\mathbf{D}_\mu$の左上$K^* \times K^*$, $\mathbf{M}$の左$K^*$列を取ってきて
 $\mathbf{A}_{[K^*]} \equiv \mathbf{N}_{[K^*]} \mathbf{D}_{\mu [K^*]} \mathbf{M}^T_{[K^*]}$
とし、$\mathbf{N}_{[K^*]} \mathbf{D}_\mu [K^*]$の行を点、列を次元とみなしなさい。なぜなら、ランク$K^*$以下のあらゆる行列$\mathbf{X}$のなかで、二乗距離
 $trace\{\mathbf{\Omega}(\mathbf{A}-\mathbf{X}) \mathbf{\Phi}(\mathbf{A}-\mathbf{X})^T\}$
 を最も小さくする$\mathbf{X}$こそ、$\mathbf{A}_{[K^*]}$だからである。

行プロファイル行列と列プロファイル行列の一般化SVD
 行プロファイル行列$\mathbf{R} \equiv \mathbf{D}^{-1}_r \mathbf{P}$について考えよう。あらかじめ$\mathbf{R}$の各列$j$からその列の重心$c_j$を引き($\mathbf{R} - \mathbf{1} \mathbf{c}^T $)、各列の重心を0にしておく。
 おさらいしておくと、点$i$の質量は$r_i$で、次元$j$の重みは$1/c_j$であった。ってことは、主よ、行$i$に重み$r_i$, 列$j$に重み$1/c_j$をつけた一般化SVDで
 $\mathbf{D}^{-1}_r \mathbf{P} - 1\mathbf{c}^T = \mathbf{L} \mathbf{D}_\phi \mathbf{M}^T $ ただし$\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{M}^T \mathbf{D}_c^{-1} \mathbf{M} = \mathbf{I}$ [式1]
と分解し、$\mathbf{F} \equiv \mathbf{L} \mathbf{D}_\phi$を座標にすればよろしいのですね。
 同様に、列プロファイル行列$\mathbf{C} \equiv \mathbf{D}^{-1}_c \mathbf{P}$については、
 $\mathbf{D}^{-1}_c \mathbf{P}^T - 1\mathbf{r}^T = \mathbf{Y} \mathbf{D}_\psi \mathbf{Z}^T $ ただし$\mathbf{Y}^T \mathbf{D}_c \mathbf{Y} = \mathbf{Z}^T \mathbf{D}_r^{-1} \mathbf{Z} = \mathbf{I}$ [式2]
と分解し、$\mathbf{G} \equiv \mathbf{Y} \mathbf{D}_\psi$を座標にすればいいわけだ。

式1と式2は結局同じ行列の一般化SVDである (その1)
 さて、式1に左から$\mathbf{D}_r$を掛ける。
 $\mathbf{P} - \mathbf{rc}^T = (\mathbf{D}_r \mathbf{L}) \mathbf{D}_\phi \mathbf{M}^T$
 左辺の分解は、もはや重み$\mathbf{D}_r, \mathbf{D}^{-1}_c$の一般化SVDではない。ええ、はい、わかってます、確かにそうなんです。しかしですね、$\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{I}$の$\mathbf{L}$に左から$\mathbf{D}_r$を掛けまくるという極悪非道な真似をしますとですね、
 $(\mathbf{D}_r \mathbf{L})^T \mathbf{D}_r^{-1} (\mathbf{D}_r \mathbf{L}) = \mathbf{I}$
である。$\mathbf{D}_r$を二回も掛けた分、真ん中を逆数にして取り返していることに注意。また、
 $\mathbf{M}^T \mathbf{D}^{-1}_c \mathbf{M} = \mathbf{I}$
であることには変わりない。結局、重み$\mathbf{D}^{-1}_r, \mathbf{D}^{-1}_c$の一般化SVDになっているわけだ。[←この理屈に気が付かなかったせいで、かなり時間を食いました...]

 同様に、式2に左から$\mathbf{D}_c$を掛ける。
 $\mathbf{P}^T - \mathbf{cr}^T = (\mathbf{D}_c \mathbf{Y}) \mathbf{D}_\psi \mathbf{Z}^T$
 $\mathbf{P} - \mathbf{rc}^T = \mathbf{Z} \mathbf{D}_\psi (\mathbf{D}_c \mathbf{Y})^T$
 ここで
 $\mathbf{Z}^T \mathbf{D}_r^{-1} \mathbf{Z} = \mathbf{I}$
 $(\mathbf{D}_c \mathbf{Y})^T \mathbf{D}^{-1}_c (\mathbf{D}_c \mathbf{Y}) = \mathbf{I}$
こちらも重み$\mathbf{D}^{-1}_r, \mathbf{D}^{-1}_c$の一般化SVDになっている。

 あらやだわ、二人とも同じ行列$\mathbf{P} - \mathbf{rc}^T$を、同じ重み$\mathbf{D}_r^{-1}, \mathbf{D}_c^{-1}$で一般化SVDしてたのね。おほほほほ。
 というわけで、一般化SVDの解が一意に決まるとして、
 $\mathbf{D}_\phi = \mathbf{D}_\psi=\mathbf{D}_\mu$
 $\mathbf{D}_r \mathbf{L} = \mathbf{Z} = \mathbf{A}$
 $\mathbf{M} = \mathbf{D}_c \mathbf{Y} = \mathbf{B}$
と置き、
 $\mathbf{P} - \mathbf{rc}^T = \mathbf{A} \mathbf{D}_\mu \mathbf{B}^T$ ただし $\mathbf{A}^T \mathbf{D}_r^{-1} \mathbf{A} = \mathbf{B}^T \mathbf{D}_c^{-1} \mathbf{B} = \mathbf{I}$
 厳密にいうと、一般化SVDの解が一意に決まらない場合(特異値のなかに同じ値がある場合)もありうるけど、要するに上のように書ける。

 座標は次のようになる。
 行プロファイルの座標: $\mathbf{F} = \mathbf{L} \mathbf{D}_\phi = \mathbf{D}_r^{-1} \mathbf{A} \mathbf{D}_\mu$
 列プロファイルの座標: $\mathbf{G} = \mathbf{Y} \mathbf{D}_\psi = \mathbf{D}_c^{-1} \mathbf{B} \mathbf{D}_\mu$

式1と式2は結局同じ行列の一般化SVDである (その2)
 ここでChap.4 からちょっと離れるけど...

 式1に右から$\mathbf{D}^{-1}_c$を掛けたら何が起きるか。
 $\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T = \mathbf{L} \mathbf{D}_\phi (\mathbf{M}^T \mathbf{D}^{-1}_c) = \mathbf{L} \mathbf{D}_\phi (\mathbf{D}^{-1}_c \mathbf{M})^T $
ここで
 $\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{I}$
 $(\mathbf{D}^{-1}_c \mathbf{M})^T \mathbf{D}_c (\mathbf{D}^{-1}_c \mathbf{M}) = \mathbf{I}$
というわけで、こんどは重み$\mathbf{D}_r, \mathbf{D}_c$の一般化SVDになる。

 同様に、式2に右から$\mathbf{D}^{-1}_r$を掛ける。
 $\mathbf{D}^{-1}_c \mathbf{P}^T \mathbf{D}^{-1}_r - 11^T = \mathbf{Y} \mathbf{D}_\psi (\mathbf{Z}^T \mathbf{D}^{-1}_r) = \mathbf{Y} \mathbf{D}_\psi (\mathbf{D}^{-1}_r \mathbf{Z})^T$
 $\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T = (\mathbf{D}^{-1}_r \mathbf{Z}) \mathbf{D}_\psi \mathbf{Y}^T$
ここで
 $(\mathbf{D}^{-1}_r \mathbf{Z})^T \mathbf{D}_r (\mathbf{D}^{-1}_r \mathbf{Z}) = \mathbf{I}$
 $\mathbf{Y}^T \mathbf{D}_c \mathbf{Y} = \mathbf{I}$
というわけで、こんどは重み$\mathbf{D}_r, \mathbf{D}_c$の一般化SVDになる。くどいですね。

 あらやだわ二人とも同じ行列$\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T$を同じ重み$\mathbf{D}_r, \mathbf{D}_c$で一般化SVDしてたのねおほほほほ、というわけで、
 $\mathbf{D}_\phi = \mathbf{D}_\psi=\mathbf{D}_\mu$
 $\mathbf{L} = \mathbf{D}^{-1}_r \mathbf{Z} = \mathbf{A'}$
 $\mathbf{D}^{-1}_c \mathbf{M} = \mathbf{Y} = \mathbf{B'}$
と置き、
 $\mathbf{D}^{-1}_c \mathbf{P} \mathbf{D}^{-1}_r - 11^T = \mathbf{A'} \mathbf{D}_\mu \mathbf{B'}^T$ ただし $\mathbf{A'}^T \mathbf{D}_r \mathbf{A'} = \mathbf{B'}^T \mathbf{D}_c \mathbf{B'} = \mathbf{I}$
という一般化SVDで表現できる。

 座標は次のようになる。
 行プロファイルの座標: $\mathbf{F} = \mathbf{L} \mathbf{D}_\phi = \mathbf{A'} \mathbf{D}_\mu$
 列プロファイルの座標: $\mathbf{G} = \mathbf{Y} \mathbf{D}_\psi = \mathbf{B'} \mathbf{D}_\mu$
つまり、一般化SVDで得られた左特異行列と右特異行列に特異値を掛けるだけで、座標が得られるわけである。その点ではこの表現のほうがスマートである。実際、この本のAppendix Aではこういう表現になっている。

 Chap.4 に戻り、最後に2つだけメモしておこう。

遷移方程式
 導出は省略するけど、
 $\mathbf{G} = \mathbf{C} \mathbf{F} \mathbf{D}^{-1}_\mu$
 $\mathbf{F} = \mathbf{R} \mathbf{G} \mathbf{D}^{-1}_\mu$
という関係が成り立つ。つまり、列$j$の列プロファイル(=$\mathbf{G}$の$j$行目)は、行プロファイル行列$\mathbf{F}$の各行の重心$\tilde{\mathbf{c}}_j F$を、次元$k$について$1/\mu_k$したものなのだ。うわーすごいわねー[←出来レース]。これを遷移方程式という。

中心化せずに一般化SVDしたら?
 さきほどは$\mathbf{P}-\mathbf{rc}^T$を一般化SVDしたが、うっかり$\mathbf{P}$を一般化SVDしちゃったらどうなるか。
 $\mathbf{P} = \tilde{\mathbf{A}} \mathbf{D}_\tilde{\mu} \tilde{\mathbf{B}}^T$ ただし$\tilde{\mathbf{A}}^T \mathbf{D}_r^{-1} \tilde{\mathbf{A}} = \tilde{\mathbf{B}}^T \mathbf{D}_c^{-1} \tilde{\mathbf{B}} = I$
 めんどくさいので証明は省略するけど、$\tilde{\mathbf{A}}$は「$\mathbf{A}$の左にベクトル$\mathbf{r}$をくっつけた行列」、$\tilde{\mathbf{B}}$は「$\mathbf{B}$の左にベクトル$\mathbf{r}$をくっつけた行列」となり、第一特異値は1となり、元の特異値は一個づつ後ろにずれるのである。ああ、なんかわかるような気がしてきた。中心化してないので「質量因子」みたいのが出てきちゃうわけだ。

雑記:データ解析 - 覚え書き:対応分析を一般化特異値分解で説明する

2018年5月11日 (金)

Nishisato, S., Clavel, J.G. (2003) A note on between-set distances in dual scaling and correspondence analysis. Behaviormetrika, 30(1), 87-98.
 調べ物の都合で読んだ。朝野先生がSASのサイトに載せていた、対応分析についてのコラムで引用されていたもの。
 第一著者は世界的に知られる偉大な計量心理学者。ここだけの話(ここってどこだ)、一般向けセミナーの聴講者のなかになぜかこの先生がおられて、穏やかな笑みを浮かべておられて...という背筋も凍る怪談を、私は複数の方から別々に聞いたことがある。起業家向けセミナーの聴衆に孫正義がいたというような話である。正直、もう怖くてセミナーなんてできない!と思いました。

 いわく。
 カテゴリカルデータを数量化する際、行変数の空間と列変数の空間はふつう異なる。行間の距離を座標から求めること、列間の距離を座標から求めることはできるけど、行と列の距離を座標から求めることはできない。
 いっぽう、行変数を列変数に射影することはできるし、逆もできる。これがいわゆる非対称的スケーリングであり、真の意味での同時グラフである。しかし、射影された変数のノルムは元の変数よりふつう小さくなる。つまり射影された変数は原点付近に集まる傾向がある。対応する特異値が小さいほどこのずれは大きくなる。言い換えると、最初のふたつの特異値がどちらも1に近いとき(上位2成分をプロットしたときとか), normed weightsとprojected weightsの exact plot はreasonableだが、第3成分と第5成分のexact plotを描いた日なんぞは、normed weightsの分散とprojected weightsの分散とのギャップが大きくなり、距離の比較はもはやreasonableでなくなる。[←このくだり、理解できたかどうか自信がない... 二元クロス表$Y$を総和で割って$P$とし、$P$の行和と列和を持たせた対角行列を$D_r$, $D_c$とし、$P$をうまいこと中心化したうえで、$D_r, D_c$を重みとして$NDM'$と一般化SVDしたとしよう。特異行列$N, M$がnormed weightsで, それに特異値を掛けた$ND, MD$がprojected weightsで、たとえば行のprojected weightsと列のnormed weightsを同時布置するのは正しいが、行の布置の分散は絶対に小さくなるよ。ということだと思うのだけれど、あっているだろうか]

 もっとも一般的な妥協案は、対称的スケーリングによって同時プロットをつくることである。その場合、2つの異なる空間が重ね書きされているということに気をつけないといけない。それはexactな記述ではなくて、practicalな近似だ。
 行と列を同じ空間に布置しようという試みの一つが、ご存じCGSスケーリングだが、明確な理論的基盤がなかった。Nishisato(1980 双対尺度法本), Greenacre(1989, JMRに載った批判論文)をみよ。実はNishisato(1990, 1997 いずれも手に入りにくそうな論文集)でその方法を提案している。
 本論文の主旨は:行と列との距離を求める方法はあるんだよ。あまりに単純なんだけど。

 $F$を$n \times m$分割表とする。成分$k$について、$y_k$を規準化された行ウェイト(長さ$n$)、$x_k$を規準化された列ウェイト(長さ$m$)とする。$y_k$で重みづけた反応と、$x_k$で重みづけた反応との相関$\rho$を最大化するように$x_k, y_k$を決めたい。
 良く知られた方法はこうである。$F$の行和をいれた対角行列を$D_r$, 列和をいれた対角行列を$D_c$とし、$D_r^{-1/2} F D_c^{-1/2}$の第$k$特異値を$\rho_k$として、
 $y_k = (1/\rho_k) D_r^{-1} F x_k$
 $x_k = (1/\rho_k) D_c^{-1} F' y_k$
Nishisato(1980)はこれをdual relationsと呼んでいる。

 $y_k, x_k$はnormed weights(ないし標準座標)と呼ばれている。なぜそう呼ぶかというと、$y'_k D_r y_k$と$x'_k D_c x_k$が事前に決めた定数と等しくなるからである。これに特異値を掛けた奴($\rho_k y_k, \rho_k x_k$)をprojected weights(ないし主座標)と呼ぶ。
 $y_k$と$\rho_k x_k$は同じ空間を張るし、$\rho_k y_k$と$x_k$も同じ空間を張る。でも、ある成分について行の軸と列の軸は$cos^{-1}\rho_k$だけ離れる。Nishisato(1996, Psychometrika)をみよ。
 だから行変数と列変数の同時布置は、ひとつの成分について2次元になり、ふたつの成分について3次元になる。いいかえると、行空間と列空間は違うもので、ある成分セットについて完全に記述するためにはもうひとつ次元を追加する必要がある。

 行と行の距離、列と列の距離についてはふつうカイ二乗距離を使う。行$y_i$と$y_j$があり、和の割合が$p_i, p_j$だとして、
 $\displaystyle d_{ij}^2 = \sum_k \rho^2_k (\frac{y_{ik}}{\sqrt{p_i}} - \frac{y_{jk}}{\sqrt{p_j}})^2$
 ところで、いま二次元空間上に点AとBがあり、原点からそれぞれまでの距離が$a, b$, AOBのなす角度が$\theta$だとしたら、AB間の距離は
 $d_{AB}^2 = a^2+b^2-2ab \cos \theta$
 さて、2成分について行空間と列空間を張ったとして、行$y_i$と列$x_j$とのカイ二乗距離を求めると、次のようになる。それぞれの和の割合を$p_{i.}, p_{.j}$として...[導出プロセスで落ちこぼれちゃったんだけど、先生を信じます]
 $\displaystyle d_{ij}^2 = \sum_k \rho^2_k \left( \frac{x^2_{jk}}{p_{.j}} + \frac{y^2_{jk}}{p_{i.}} - 2 \rho_k \frac{ x_{jk} y_{ik} }{ \sqrt{ p_{.j} p_{i.} } } \right)^2$

 [数値例。Rで途中まで試してみたけど省略]

 というわけで、行と行、列と列、行と列について距離行列をつくれる。これを非計量MDSで観察してみよう...[略]

 このように、行空間と列空間の乖離を表現するためには次元をひとつ追加する必要があり、その次元の寄与は、特異値が小さいときほど大きい。思うに、対称的スケーリングで同時プロットを描いていいのは特異値がコレコレ以上の時だ、というようなガイドラインがあったほうがいいのではないか。
 双対尺度法・対応分析は行と列の距離を無視している。現時点で、我々は行と行、列と列、行と列のカイ二乗距離を同時に分析する方法を持っていない。今後の発展が望まれる。云々。

 ...なんとか読み終えたけれど、やはりもっときちんと勉強しておかないとだめだなあ...

論文:データ解析(2018-) - 読了:Nishisato & Clavel (2003) クロス表の対応分析の結果から「行と列の距離」を求めるには

西里静彦(2014) 行動科学への数理の応用:探索的データ解析と測度の関係の理解. 行動計量学, 41(2), 89-102.
 前に読んだ気もするんだけど、調べものの都合で読んだ。前後の文脈を調べてみると、2014年に行動計量学会が開いたワークショップを踏まえた寄稿らしく、題名はそのワークショップのタイトルに由来したもの。論文の内容としては、データを名義尺度で表現しなさい、十分な多次元空間を用いなさい、行と列を同時空間で見なさい、Nishisato & Clavel(2004)の全情報解析へようこそ、という感じ。

いくつかメモ:

論文:データ解析(2018-) - 読了:西里(2014) 全情報解析への招待

2018年5月 8日 (火)

Baillon, A. (2017) Bayesian markets to elicit private information. PNAS, 114(30), 7958-7962.
 アブストラクトに目を通して青くなった。ベイジアン自白剤と予測市場の合いの子という、私の心のどまんなかを撃ち抜く論文。これ去年の6月じゃん。なぜこれに気が付かなかったんだ...

 ベイジアン市場を提案します。二値の私秘的情報を引き出すための市場です。そんじゅそこらの予測市場とは異なり、結果についての客観的検証ができない場合も大丈夫です。
 ベイジアン市場の基盤にあるのは、私秘的情報は他者についての信念に影響する、というベイジアン推論の想定です[ここでDawes(1989 JESP)を引用]。いまある事柄にYesと答える人は、その事柄に対する他者のYes割合についての期待を更新する際に自分の答えを使います。[←おおお、ベイジアン自白剤と全く同じ話だ]
 予測市場では、あるイベントの賭けがその人の信念を表します。ベイジアン市場では、他の人の回答への賭けが、他者についての信念を表し、ひいては当該の問いへのその人の真の答えを表します。

 私秘情報を引き出す手法としてはすでにベイジアン自白剤やピア予測法があります[Prelec(2004 Sci.), Miller, Resnick, & Zeckhauser(2005, MgmtSci), Parkes & Witkowski (2012 Proc.AAAI)(←たぶん Witkowski & Parkes(2012)の間違い), Radanovic & Faltings(2013 Proc.AAAI; 2014 Proc.AAAI)]。でも確率推定やメタ信念推定をしているぶん複雑です。いっぽう提案手法はただの賭けなので単純。ただし二値の質問限定です。

 エージェントの数を$n$とする。私秘情報についての二値設問を$Q$、値を{0,1}とする。$i$が持つ真の情報(=$i$のタイプ)を$t_i$とし、$\omega = \sum^n t_i/n$とする。
 先行研究と同じく、すべてのエージェントは「自分のタイプを知らない場合の事前信念」$f(\omega)$を共有していると仮定する。なおこの事前信念の共有という仮定は Harsanyi(1968 MgmtSci)が支持しているぞ。
 [次の段落は大事なので全訳]

Prelec(2004)と同様に、タイプが非個人的に情報的だということ、すなわち$f(\omega|t_i)=f(\omega|t_j)$と$t_i=t_j$が等価だということ、を共通知識とする。この特性は2つの側面を含んでいる。
 第一に、タイプは非個人的である。$t_i=0$であるすべてのエージェント$i$は共通の更新後信念$f(\omega|t_i=0)$を持ち(その期待値を$\bar{\omega}_0$と書く)、$t_j=1$であるすべてのエージェント$j$は共通の更新後信念$f(\omega|t_j=1)$を持つ(その期待値を$\bar{\omega}_1$と書く)。このように、エージェントのタイプはすべての非共有情報を含んでいる。
 第二に、タイプは情報的である(ないし「確率的に関連性を持つ」)。エージェント$i$のタイプが1ならば、このシグナルのせいで彼は、$\omega$は彼が事前に想定していたよりも大きな割合だと考えるようになる。いっぽうタイプ0のエージェントは小さな割合だと考えるようになる。よって$\bar{\omega}_0 < \bar{\omega}_1$である。

 話を単純にするため、$n$は無限大であり、$f$は「すべて0」や「すべて1」でないと仮定する。

 提案手法。
 $Q$についての市場をつくる。全員が同時に参加するワンショット市場である。
 参加者はあるアセットを取引できる。そのアセットとは、価値$v$が「1と報告する人の割合」であるアセットである。
 この市場では、参加者は主観的な期待ペイオフを最大化する、参加者は主観的期待ペイオフが正の時しか市場に参加しない、というのが共有知識になっている。

  1. 参加者はまず回答$r_i$を報告する。
  2. 次に、$p$が一様分布からランダムにドローされる[←すごく混乱したんだけど、この$p$がマーケットメーカの提案価格であり、全参加者に対して共通なのだ]。
  3. $r_i$が1だったら「価格$p$でアセットを買うか」、0だったら「価格$p$でアセットを売るか」を問われる。[←これはいわば注文であって、成立するとは限らない]
  4. すべての取引はマーケット・メーカとの間でなされる。取引が成立するかどうかはあるルールで決まる。
  5. アセットを清算する。清算価格は$r_i$における1の割合とする。つまり、アセットを買った人に$v$を配り、売った人から$v$を徴収する。買い手の手元には$v-p$, 売り手の手元には$p-v$が残る。

さて、取引の成立・不成立を決めるルールとは...

話の先取りになるけど、ここで「多数派」というのを「全員一致」に置き換えても、「三分の一以上」に置き換えても、実はこの論文の結果は変わらない。実装の上で「多数派」としておくのが自然なだけで。

 結論からいうと、この市場では真実報告がベイジアン・ナッシュ均衡になる。以下、その説明。[毎度のことながら頭がこんがらがってくるので全訳する]

 まず、すべてのエージェントが市場に参加すると仮定する。後述するように、すべての期待ペイオフは0より大なので、エージェントは実際に参加することは保証されている。
 エージェント$i$について考える。他のすべてのエージェントは真実を報告する、すなわち$v=\omega$と仮定する。
 タイプ1のエージェントは買い手側となり、市場価格$p$が、アセットの価値についての彼らの期待値$\bar{\omega}_1$を下回るときに買い注文を出すだろう。同様に、タイプ0のエージェントは売り手側となり、市場価格が$\bar{\omega}_0$を上回るときに売り注文を出すだろう。
 仮定により、両方のタイプのエージェントがいることは確実であり、それぞれの側の多数派が取引を求めた時、マーケット・メーカが[逆側の]取引希望者の全員と取引することも確実である。従って、取引が生じるのは$\bar{\omega}_0 < p < \bar{\omega}_1$のとき、そのときに限られる。
 エージェント$i$はどうすべきだろうか? もし彼のタイプが1ならば、彼はアセットの価値を$\bar{\omega}_1$と期待する。売り手として利益を出すためには、彼は市場価格[$p$]が$\bar{\omega}_1$を上回ったときだけ売り注文を出すことになるが、そんな高値での取引は起こらないだろう。しかし、市場価格が$\bar{\omega}_1$までであれば彼は買い注文を出したい。従って彼は、市場価格が$\bar{\omega}_0$ と $\bar{\omega}_1$の間である時に取引で利益が得られると期待する。
 このペイオフを獲得するためには、彼はまず1と報告しなければならない。すなわち、真実を報告しなければならない。そうすれば、彼は
 $\int_{\bar{\omega}_0}^{\bar{\omega}_1} (\bar{\omega}_1 - p)dp = (\bar{\omega}_1-\bar{\omega}_0)^2)/2 > 0$
を受け取ると期待できる。つまり、$\bar{\omega}_0$と$\bar{\omega}_1$との間にあるすべての市場価格$p$について、ペイオフは$\bar{\omega}_1 - p$となる。
 いわゆる混合方略(ある確率で1と報告し、そうでないときに0と報告する)には意味が何。なぜならそれはペイオフを獲得する確率を低くするだけだからだ。
 タイプ0のエージェントにとっても真実申告が最良である。その証明は上と対称であり、期待ペイオフも同一である。

 [よーし、具体例で考えてみよう。
 人々に「幸せですか」と訊ねる。全員から回答を集めたのち、マーケット・メーカの提案価格(1円~99円)をルーレットか何かで決め、全員に提示する。

取引後に証券は清算される。清算価格は提案価格や取引とは無関係に、最初に「はい」と答えた人の割合で決まる。ここまでがルールね。
 架空例として、参加者のなかには幸せな人が60%、そうでない人が40%いるとしよう($\omega=0.6$)。このことを参加者は知らない。幸せな人は幸せ率を70%と見積もり($\bar{\omega}_1=0.7$)、不幸な人は幸せ率を50%と見積もる($\bar{\omega}_0=0.5$)、としよう。話を簡単にするために、このことを参加者全員が知っているとする(本当は$\bar{\omega}_0 < \bar{\omega}_1$であるということさえ知っていれば良い)。
 太郎は幸せである。ゆえに幸せ率は0.7、清算価格は70円となるとみている。太郎は考える:

ということは、正直に回答したほうが得だ。
 次郎は不幸である。ゆえに幸せ率は0.5, 清算価格は50円となるとみている。次郎は考える:

ということは、正直に回答したほうが得だ。
 というわけで、全員が正直に回答することになる。]

 なお、

 他の研究との比較。

 前提となる仮定について。

 最後に、応用領域について。

 。。。うっわー。。。 これ、面白い。。。
 著者はフランス出身の経済学者で、すごく若そうな人。実をいうとしばらく前に、エラスムス大ロッテルダム校の院生によるベイジアン自白剤の修論というのがネット検索でひっかかり、まあ修論なら読まなくてもいっかとブラウザのタブを閉じてから、おいちょっと待て、それを指導してる研究者がいるってことじゃん、それ誰? とひっかかっていたのである。ああ、いたよ、ここに張本人が。
 あまりに面白いので、あれこれ考え込んでしまい、まだ感想がまとまらない。とりあえずいま気になっているのは、実証実験はあるのか、やったらどうなるのかということだ。この手法はBDMメカニズムに似ていると思うんだけど、BDMメカニズムは実証的には必ずしも機能しないと聞いたことがある。
 ともあれ、論文本文を読んだメモとして記録しておく。次は付録を読もう。いやしかし、これ、面白いなあ...

論文:予測市場 - 読了:Baillon (2017) ベイジアン・マーケット

2018年5月 5日 (土)

 SNSを検索していて、「マスメディアの世論調査はあてにならない」という強烈な感情表現を数多く目の当たりにし(その論拠はピンキリだが、ほとんどはキリのほう)、サーベイ調査に対する人々の信頼とは... と、深い感傷に浸ったのであった。

Nir, L. (2011) Motivated reasoning and public opinion perception. Public Opinion Quarterly, 75(3), 504-532.
 とはいえ、ぼけーっとしていてもしょうがないので、これを機にちょっと勉強してみるか、と関連しそうな論文を探して読んでみた。予備知識ゼロだけど、タイトルからして面白そうなので。

 ハバーマスいわく[おっと、カッコ良い出だしだね]、publicとはただの個人の集まりでない。publicはコミュニケーション行為を通じて意見を形成する。それを可能にする主要な情報源がマスメディアである。
 しかし、個々人が他人の政治的選好についての情報を求めるかどうかは別の問題である。人々が自分の意見を支持する意見ばかりに耳を傾けるようでは、熟議は成り立たない。
 motivated reasoning(目標志向的な認知処理方略)の研究では、accuracy目標とdirectional目標のちがいが注目されている。前者は正しい信念を維持しようという目標で、推論に投じられる努力が大きくなり、情報処理が深くなる。後者は望ましい結論を維持しようという目標で、それを支持する証拠が重視される。Kunda(1990, Psych.Bull.)をみよ。
 Lodge & Taber (2000)はmotivated reasoningと政治的評価のタイポロジーを提案している。横軸にaccuracy目標の強さをとり、縦軸にdirectional目標の強さをとる。

 先行研究が主に注目してきたのは党派的推論者の情報処理だった。意見にバイアスを与えるメカニズムとして、選択的接触、選択的処理、選択的知覚が指摘されてきた。しかし、Lodge & Taberのタイポロジーは、以下の3つの方向の研究を示唆している。

 仮説です。

 推論目標をどうやって操作的に定義するか。accuracy目標はCacioppo & Petty(1982 JPSP)のNeed for Cognition尺度で測る。direction目標はJavis, Blair, & Petty(1996 JPSP)のNeed to Evaluate尺度で測る。[←すげえ... 社会心理学版の打ち出の小槌だね...]
 
 方法。
 Electronic Dialogue 2000 (ED2K)プロジェクトのデータを使う。脚注もみながらちょっと詳しくメモしておくと...
 対象者はKnowledge Networks社のパネルから抽出した。これはまずSSI社が世帯をRDDで抽出し、そのうちWebTVのデータベースにない世帯をリクルートして、WebTVの装置をただで配るかわりに毎月40分程度の調査に参加してもらいますという全米パネル。年齢や人種が代表性を持つようにリクルートしている。
 ここから2000年2月に3967人抽出し、2014人の同意をとった。ベースライン調査は長いので2回にわけてやった(2~3月)。回答者は1684人。これを「議論」群(900人)、「コントロール」群(139人)、「セットアサイド」群に割り当てた。議論群とコントロール群は毎月調査参加する羽目になる。議論群は、10-15名ずつの60チームにわかれ、モデレータつきのオンライン議論に毎月参加する(調査はその前後)。実は60チームはうまく組んであって、20チームはリベラルor中道、20チームは保守or中道、20チームは混在である。最終の調査が2001年1月(これも2ウェーブに分けた)。実に一年がかりのプロジェクトであった。お疲れ様でした。
 Need for CognitionとNeed to Evaluateの項目は、4回目の議論のあとの調査(2001年1月)、ないし最終調査に含まれていた(うまいこと設計されていて、ある対象者にはこのうちどこかで一回だけ訊いている)。
 分析に使う指標は以下の通り。

 お待たせしました!結果です!

 考察。
 本研究の貢献は以下の3点。

  1. motivated reasoningの研究の拡張。推論タイプと意見知覚の関係を示した。
  2. 観察研究における、推論タイプの操作的定義を示した。
  3. 人々の推論過程の個人差、それが近くにに及ぼす影響についての我々の理解を示した。

 本研究の限界:

 今後の研究への示唆。

 。。。いやー、さっすがはPOQ。私のような素人の心も掴んで離さない、すごく面白い論文であった。
 Lodge & Taberのタイポロジーというのがとにかく面白い。引用されている論文はLupia, et al. (eds.) "Element of Reason"に所収。第一著者のMilton Lodgeについて調べてみると、世界のあちこちに同名の宿泊施設があって困ったが、どうやら政治心理学者。70年代にJPSPに書いているほかは、政治心理学の専門誌での論文が多いようだが、Sageの緑本でマグニチュード尺度について書いていたりする。67年に学位取得というから結構なお年であり、今年お弟子さんたちかなにかによる記念論文集が出ているようだから、さぞや偉い学者なのだろう。2013年にTaberと共著で"The rationalizing voter"という本を書いているが、残念ながら翻訳はないようだ。

 書き方から推測するに、著者自身はED2Kプロジェクトにはタッチしておらず、あとでデータを借りてきたような雰囲気だ。こんなんあれですよ、凡百の心理学者がやったら、「回帰分析してみたら、認知欲求が高い人はxxxが高くて、評価欲求が高い人はxxxが高かったですー」的な、吹けば飛ぶよな相関研究ですよ。つくづく、大事なのはデータじゃなくて理論的枠組みなのである。

 自分の仕事に当てはめると、消費者に対して他者の選好を推測させるcitizen forecastを考えた時、うっかりしていると「カテゴリ関与が高いほうが推測が正確」と素朴に仮定してしまいそうになるが、正確性を目標にした推論をもたらすようなカテゴリ関与もありうるし、自分の選好の支持を目標にした推論をもたらすようなカテゴリ関与もありうるわけだ。その違いは、もちろん環境との相互作用の履歴が生んだ産物ではあるんだけど、カテゴリを超えた個人的傾向の尺度で、ある程度は予測できるかもしれない。

論文:心理 - 読了:Nir (2011) 推論の動機づけ的個人差が、世論についての知覚のバイアスと関連している

先日ふとした機会に、人々がサーベイ調査に対して持つ信頼感ってどうやって形成されているんだろう?と不思議に思った。茫漠とはしておりますが、結構切実な疑問であります。メディア研究の文脈とかで、きっと理論とか実証研究とかあるんだろうけれど、探し方が悪いのか、ぴったりなのが見つからず...

稲増一憲(2016) メディア・世論調査への不信の多面性:社会調査データの分析から. 放送メディア研究, 13, 177-193.
2015年5月、NHK放送文化研究所の調査でメディア・世論調査への信頼について調べたという報告。著者は社会心理の若い研究者の方。Webで拾って読んでみた。

 調査は住基台帳ベース、留置法。さっすがー。
 報道媒体への信頼度を5件法で訊いたんだけど、TVニュースと新聞記事に対する回答は意外に高く、T2B(「とても信頼している」「まあ信頼している」で7割越え。
 新聞社への信頼とテレビ局への信頼(4件法、こちらはT2Bが3割台)を目的変数として重回帰した(順序ロジットとかではなくて、単純なOLS)。テレビニュースや新聞記事への接触がポジティブに効き、ネットでの他者の意見への接触頻度がネガティブに効いた(ネットの政治ニュースへの接触は効いてない)。テレビ局への信頼には民主党支持がポジティブに効き、政治関心の高さがネガティブに効いた(性・年齢・年齢二乗・学歴・年収をコントロールしてもなお。へー。でもこれ、性x年代カテゴリで入れたほうが良かったんじゃないかな、きっと交互作用があるから...)
 途中省略するけど、世論調査の中立性の知覚について調べているところがあるのでメモ。世論調査を「マスメディアが操作している」という項目(5件法)を目的変数とした重回帰で、政治への関心の高さ(自己報告)がポジティブに効き、実際に持っている政治についての知識レベル(簡単なテストで調べている)がネガティブに効く。なるほどねえ。
 なお、上記に支持政党は効いていない。ただし、世論調査が「公平に意見を反映している」という項目には、自民党・公明党支持がポジティブに効いている(上と同じく、デモグラはコントロール済なのに)。へー、そうなのか。
 考察。Ladd(2012) "Why Americans Hate the Media and How it Matters"という本で、アメリカでは党派的に分極化した政治家たちのメディア批判がメディア不信の一因になっているという指摘があるんだけど、この調査結果をみるかぎり、日本ではそうはなっていない模様。云々。

 いまの日本の首相は国会での質疑で、別に聞かれてもいないのに特定の新聞社を口汚く罵るという、私にはよく理解できない振る舞いをなさっているのだが、ああいうのって、メディアへの信頼と内閣支持率にどう効くんでしょうね。反応の異質性がすごく大きいので、簡単には答えの出ない話だろうけど。

論文:その他 - 読了:稲増(2016) メディア・世論調査への不信感が高いのはどんな人か

Liu, I., Agresti, A. (2005) The analysis of ordered categorical data: An overview and a survey of recent development. Test, 14(1), 2005.
 順序カテゴリカルデータの分析手法についてのレビュー。
 かつてAgrestiの教科書とかで読んだ内容とも重複しているんだけど、こういう知識は突然に必要になるので、時々おさらいしておこうと思って... というのは建前で、いろいろ疲れることが重なったあげく、不意に出先で時間が空いたので、リハビリのつもりで読んだ奴。こういう話題はややこしいことを考えなくていいので助かる。
 掲載誌はSpanish Society of Statistics and Operations Researchの発行。Google様的な被引用回数は241件、結構少ない。

1. イントロダクション (略)

2. 順序カテゴリカル反応のモデル
 順序カテゴリカル反応への関心を引き起こした主な研究がふたつある。McCullagh(1980)による累積確率のロジットモデルと、Goodman(1979)によるオッズ比の対数線型モデルである。

2.1 比例オッズモデル(累積ロジットモデル)
 順序反応のモデルとして現在もっとも一般的なのは累積確率のロジットを使うモデルであろう。
 McCullaghはこう考えた。$c$カテゴリの順序反応を$Y$、予測子のセットを$x$として、
 $logit[P(Y \leq j | x)] = \alpha_j - \beta' x$
$j$は$1$から$c-1$まで動く。
 このモデルは$c-1$個の累積確率のすべてに対して予測子が同一の効果を持つと仮定している。$c$カテゴリの反応をどこかで切って二値に落とすとして、どこで切ろうが説明変数の効果を表すオッズ比は同じだということになる。このタイプのモデルは比例オッズモデルと呼ばれることも多い。

2.2 累積リンクモデル
 McCullagh(1980)は、プロビット、loglog, cloglogといった二値データでおなじみのリンク関数を使うという手も考えた。そういうのを取り込んで一般化したのが累積リンクモデルである。
 $G^{-1}[P(Y \leq j | x)] = \alpha_j - \beta' x$
$G$は連続累積分布関数(cdf)。

2.3 多項ロジットモデルによる代替案

 いま、カテゴリカル反応を長さ$c-1$のダミーベクトルで表現し、対象者$i$におけるその平均を$\mu_i$として
 $g(\mu_i) = X_i \beta$
という多変量一般化線型モデルを考えると、累積リンクモデルも隣接カテゴリロジットモデルもcontinuation-ratioロジットモデルも、このモデルの特殊ケースとして捉えられる。

2.4 その他の多項反応モデル

2.5 順序反応の連関のモデリング
 $r \times c$分割表の行を$X$, 列を$Y$とし、セル頻度を$\{n_{ij}\}$, その期待値を$\{\mu_{ij}\}$とする。
 Goodman(1985)はこう考えた。
 $\log \mu_{ij} = \lambda + \lambda_i^{X} + \lambda_j^{Y} + \sum_k^M \beta_k u_{ik} v_{jk}$
ただし$M \leq min(r-1, c-1)$。$M = min(r-1, c-1)$とすると飽和する。ふつうは$M=1$とする。

関連する研究は対応分析と正準相関モデルの文脈でも行われている。一般正準相関モデルは
 $\pi_{ij} = \pi_{i+} \pi_{+j} \left(1+\sum_k^M \lambda_k u_{ik} v_{jk}\right)$
Goodmanの定式化と似ているけど、連関の項を使って$\mu_{ij}$と独立なときの値とのずれをモデル化している点に注意。この路線の限界は、多元表への一般化が大変だという点。

3. 順序反応のクラスタ化データ・反復データのモデリング
 クラスタ化データ・反復データのモデルには大きく分けて2つある。

 以下、あるクラスタにおいて$Y_1,\ldots,Y_T$の反復反応が得られているとする。話を簡単にするため$T$は全クラスタで共通だということにするが、クラスタ$i$について$T_i$だと一般化してもよろしい。

3.1 GEEアプローチによる周辺モデル
 累積ロジットリンクを持つ周辺モデルについて考えよう。
 $logit[P(Y_t \leq j|x_t)] = \alpha_j - \beta' x_t$
 これだと$T$個の反応の周辺分布の説明変数への依存性だけを考えていて、$T$回の反復反応の間の多変量的な依存性はモデル化していない。
 さて、この周辺モデルを推定する際、対数尤度関数を最大化するというのはおかしな話になる。対数尤度関数は、予測子の多様な水準から得られた多項分布の積になる。それぞれの多項分布は、$T$個の反応のクロス分類における$c^T$個のセルから定義される。尤度関数はこの完全な同時分布を指しているのだから、周辺モデルの公式を対数尤度に直接置き換えることができない。[←ううむ...よくわからない...]
 ここでは一般化推定方程式(GEE)を使うのが簡単である。GEEは擬似尤度の多変量への一般化に基づいており、周辺回帰モデルだけを指定すれば、$T$個の反応の相関構造については適当に推測してくれる。
 GEEはもともと、二値とかポワソンとかの単変量分布の周辺モデルとして開発されたが、反復順序反応の累積ロジットモデルや累積プロビットモデルに拡張されている。
 [ここからGEEによるモデリングの説明に入るんだけど、残念ながら良く理解できなかった。短い説明だし、写経してもいいんだけど、面倒なので省略する。ううむ、GEEっていつも勉強しかけて挫折するんだよなあ]

3.2 MLアプローチによる周辺モデル
 上記のように、周辺モデルを最尤推定するのはおかしい。同時セル確率と、同時分布の高次パラメータならびに周辺モデルのパラメータとを一対一対応させる多変量ロジスティックモデルという提案もあるが、次元数が増えると無理になる。周辺モデルを制約方程式として扱い、ポワソン尤度と多項尤度を最大化するというアプローチもあるが、計算が大変。
 [だめだ、3.1-3.2はさっぱり理解できない。Agrestiの教科書も見たけどやっぱりむずかしい。まあふだんはGLMMでやっておりますので、今回はあきらめるけど、なにかやさしい教科書はないかしらん]

3.3 一般化線形混合モデル
 こんどは周辺分布じゃなくて、同時分布をクラスタのランダム効果を使ってモデル化する。反応を指数分布族として、一般化線形モデル(GLM)にランダム効果を加えたのを一般化線形混合モデル(GLMM)という。同様に、多変量GLMも多変量GLMMに拡張できる。
 ランダム効果が表すものは状況によって異なる。ヒトの異質性とか、未知の共変量とか、なんらかの過分散(overdispersion)とか。
 基本的なモデルは
 $logit[P(Y_{it} \leq j | x_{it}, z_{it})] = \alpha_j - \beta^{'} x_{it} - u^{'}_i z_{it}$
ここで$z_{it}$はランダム効果の説明変数ベクトルで、$u_i$はiidに$MVN(0, \Sigma)$に従う。リンク関数を変えてもいいし、continuation-ratioロジットモデルに変えてもいい。[嗚呼... 周辺モデルよか全然わかりやすいよ... 胸なごむ...]

3.4 マルチレベルモデル
 クラスタのレベルが2つ以上ある場合もある。地域-学校-学生、とか。
 アプローチとしては、GEE, 多変量GLMM, 階層モデルでそれぞれのレベルについて事前分布を与えるベイズモデル、が考えられる。

3.5 その他のモデル(遷移モデルと時系列)
 反復測定データを遷移モデルで分析するという手もある。反応を説明変数と過去の反応で条件づけて説明するわけだ。二値データについてマルコフ連鎖構造を使うことが多い。
 順序データについてはEkholm et al.(2003)というのがあって、反応間の連関をdependence ratioというので表す。これは所与のセルの確率を独立性の下での期待値で割った値のこと。このモデルは、連関メカニズムが純粋に交換可能であればヒトのランダム項をいれたGLMMと等しくなり、連関メカニズムが純粋にマルコフであればマルコフ連鎖構造の遷移モデルと等しくなる。[←悪いけどなにをゆうておるのかようわからん]
 連続時間の遷移モデルのGEE推定とML推定というのもある。

 周辺モデルが良いか、クラスタ特定的なモデルが良いか、遷移モデルが良いか。これは、母集団レベルで解釈したいか個人レベルで解釈したいか、説明変数の効果を以前の反応で条件づけて記述したいか、によって決まる問題である。
 モデルのタイプがちがえばパラメータのサイズも変わる。たとえば、周辺モデルよりもクラスタ特定的モデルのほうが効果の推定値が大きくなる。

4. 順序反応へのベイジアン分析
4.1 多項パラメータの推定
 とりあえず、$c$カテゴリの多項変数、説明変数なし、ということにしよう。セルの頻度$(n_1, \ldots, n_c)$がサイズ$n=\sum n_i$, パラメータ$\pi=(\pi_1, \ldots, \pi_c)'$の多項分布に従うとする。多項確率質量分布の共役密度関数はディリクレで...[分布とその特性値の説明。省略]。事前のパラメータを$(\alpha_1, \ldots, \alpha_c)$とすると、事後のパラメータは$\{n_i+\alpha_i\}$となり...[特性値の説明。省略]。
 さて、多項確率の推定の際に$\pi_1 \leq \ldots \leq \pi_k \geq \pi_{k+1} \geq \ldots \geq \pi_c$という制約をいれるという考え方があり、これは順序カテゴリに適することが多い。[←えっ、なぜ??]
 多項ロジットのパラメータ事前分布としてディリクレじゃなくてMVNを仮定するという路線もある。相関行列を自己回帰の形式にすると、隣接カテゴリの確率が似てくるので、順序カテゴリに向いている。
 ディリクレ・パラメータの分布を階層的に指定するという手もある。二階の事前分布をいれるわけね。共役分布を使わないと決めてしまえば、ロジットモデルでMVN事前分布の階層モデルのほうが計算しやすい。

4.2 クロス表における確率の推定
 Good(1965)はディリクレ事前分布と階層モデルを使ってクロス表のセル確率を推定した。その後、データ依存な事前分布を使うというのも出てきて... [面倒になってきたのでこの節はスキップ]

4.3 順序反応のモデリング
 [いよいよ、順序反応そのもののベイジアン・モデリングの話、なんだけど... 1996年ごろ以降の先行研究が、一本につき数行づつ紹介されている感じなので、省略する。多変量反応にも拡張されている由]

5. 順序反応についてのモデルベースでない手法
5.1 層別クロス表のCMH法
 カテゴリカル変数$X, Y$の連関を、第三の変数$X$をコントロールして分析したいとき、ふつうは三元クロス表を書く。では、$Z$をコントロールした下での$X, Y$の条件つき連関は?
 モデルに基づかないアプローチとして一番有名なのはCochran-Mantel-Haenszel(CMH)検定である。$X$と$Y$がともに順序なら、層を結合した連関が自由度1のカイ二乗分布に従うことを利用する。これは行と列になんらかスコアを振って線形トレンドをみていることになる。実はこれ、多項ロジットと深い関係がある。Agrestiの本を読め。
 検定統計量だけでなく、後述する順序オッズ比も手に入る(それが層を通じて一定だという仮定の下で)。各層が十分大きければ多項ロジットモデルのML推定量と類似するし、層が多くなってデータがスパースになってくると、実はML推定量より優れている。

5.2 順序オッズ比
 順序反応の二元表のオッズ比にはいろんなタイプがある。

5.3 ランク・ベースのアプローチ
 先ほど述べたように、CMH統計量は連関の要約にカテゴリのスコアを使っていると捉えられるが、スコアを使わず厳密に順序情報だけで推論するのもある。Kendallの$\tau$統計量, Goodman-Kruskalの$\gamma$, Jonckheere-Terpstra検定、など。[←最後の奴、初めて聞いた...]

5.4 不等性制約の使用
 順序カテゴリを活用するその他の方法として、連関構造を記述するパラメータに不等性制約を置くというのがある。たとえば、二元クロス表での連関を記述する際、累積オッズ比、局所オッズ比、大域オッズ比、連続オッズ比が使えるが、その際に$(r-1)(c-1)$個の対数オッズ比のすべてが非負だという制約を掛ける。するとセル確率のML推定値が手に入ったり、尤度比検定で独立性を検定したりできる。Bartolucci et al(2001, JASA)をみよ。なお、局所オッズ比を制約するのがもっともきつい制約になる(局所対数オッズ比が一様に非負なら大域対数オッズ比も非負)。
 2行の場合、非負対数オッズ制約のもとで尤度比検定の漸近分布が手に入って...でも3行以上には一般化できなくて... [とかなんとか... 略]
 それとは別の話として、対数線形モデルとかのパラメータに非負制約を掛けるというのもある。

5.5 一致の測定
重みつき$\kappa$とか、ROC曲線を使うとか、いろんな方法がある。

6. その他の問題
6.1 正確推論
[ごく短い説明だけど、力尽きてきたのでパス]

6.2 欠損データ
 多変量GLMMの場合、MCARやMARがあっても大丈夫。いっぽうGEEのように尤度ベースでない手法では、ランダム欠損でも無視できるとは限らない。
 他にも欠損を扱うモデルの提案がいろいろある。反応じゃなくて共変量の欠損を扱うモデルの提案もある。

6.3 標本サイズと検定力
 $c$カテゴリの順序反応を2群間で比較する場合について考える。比例オッズモデルの場合、周辺確率が変わらないとして[←?]、ある検定力に達するために必要な標本サイズ$N(c)$には次の性質がある:
 $N(c) / N(2) = 0.75/(1-1/c^2)$
ここからわかるのは、カテゴリを2つに潰すと情報の損失は大きい、しかしカテゴリが4~5個あればそれ以上あるときと変わらない。
 ほかに、Wilcoxonのrank sum検定の場合のガイドラインとか、個人について処理前と処理後に測定すするsubjec-specificモデルの場合のガイドラインが提案されている。

6.4 順序データ分析のためのソフトウェア
[略]

6.5 結語
[略]

 ... いやあ、長かった... 疲れた... めんどくさかった...

 この論文が面白そうだと思ったのは、73頁のうち本文は29頁で、残りは識者たちのコメントと返答になっているところ。せっかくなので対話調でメモしておく。

Tutz: 順序回帰モデルの利用について伺います。
 順序モデルとして良く知られているのは累積型のモデル、つまり
 $P(Y \leq j | x) = F(\eta_j(x))$
 $\eta_j(x) = \alpha_j - x'\beta$
というものです。
 でも、逐次型のモデル
 $P(Y=j | Y \geq j, x) = F(\eta_j(x))$
も有用だと思うんですよね。ここで$F$は狭義単調分布関数であればよくて、たとえばロジスティック分布関数ならcontinuation-ratioロジットモデルになります。
 このモデルはカテゴリ$r$から$r+1$への遷移をモデリングしているとみることもできます。だから、カテゴリ$r$に到達しているという条件のもとで、どんな二値回帰モデルでも使えるんです。
 逐次型モデルは拡張が容易です。たとえば累積モデルによる近似がうまくいかない場合は、$\beta$をカテゴリごとに$\beta_j$として推定してもいいです。[以下、累積モデルに比べて逐次型モデルがいかに柔軟かという話。略]
 というわけで、累積型モデルばかりが注目されるのって納得いかないんですけど、どう思います?

Liu & Agresti: 逐次型モデルがあまり使われてないのは、カテゴリを低いほうから並べるか高いほうから並べるかで結果が変わってきちゃうからじゃないですかね。

Tutz: 予測子をどう組み合わせるかという話なんですが。
 反復測定データの周辺モデルや混合モデルの場合、一番単純なモデルだと、クラスタ$i$の観察$t$のカテゴリ$j$について
 $\eta_{tj}(x_{it}) = \alpha_j - x'_{it}\beta$
というふうに、観察$t$を通じて$\alpha_j$と$\beta$は一定にすると思います。カテゴリ数なり観察数なりが小さい場合はともかく、これはちょっと制約し過ぎと思います。先々の研究を考えたら、やっぱり
 $\eta_{tj}(x_{it}) = \alpha_{tj} - x'_{it}\beta_j$
とか
 $\eta_{tj}(x_{it}) = \alpha_{tj} - x'_{it}\beta_t$
というふうに、カテゴリ別・観察別のパラメータを入れたほうがいいと思います。でも全パラメータを自由推定しちゃうとノイズに敏感になっちゃうし解釈しにくいので、なんらか制約を入れたほうがいいでしょうね、たとえば罰則付きML法で、対数尤度$\ell$を
 $\ell_p = \ell - \sum_{j,s} \lambda_s (\beta_{j+1,s} - \beta_{j,s})^2$
に置き換えるとか($\lambda_s$は平滑化パラメータ)。
 累積型モデルで追加制約を考慮する場合、罰則付き尤度モデルは比例モデルと非比例モデルの間のどこかに位置します。$\lambda_s$がすごく大きい場合、$x_{it}$の$s$番目の要素に比例オッズモデルを当てはめることになります。もちろん$\alpha_{tj}$は$t,j$を通じて一定にしないといけません。[←??? 力不足で理解できない...]
 予測子を平滑化するという方向の拡張もあります。線形予測子の代わりに加算構造を使うというのは興味深い拡張ですね。[GEEのこと???] 横断データだと加算順序回帰モデルとしてすでに使われていますが、反復測定の場合、$t$別にスムーズな構造を見つけるというのはチャレンジングな課題です。

Liu & Agresti: 予測子のより一般的な構造を用いるという解説、ありがとうございます。おっしゃるとおり、単純・解釈可能でかつスムーズな構造を見つけるのが難しいですね。でもこれからは応用事例も増え、もっと一般的になると思いますよ。

Tutz: 平均反応モデルについて一言申し上げたい。
 実務家にとって魅力的だというのはわかりますけど、あれは順序モデルじゃないですよ。カテゴリにスコアを与えているということは、反応をmetricallyにスケールした離散反応モデルにすぎないです。

Liu & Agresti: おっしゃる通り、順序反応データに平均反応モデルを使うのはよろしくないです。でも道具箱のなかにはいれておきたいと思います。非統計家にとってはロジスティック回帰よりわかりやすいですし、反応カテゴリ数が多ければ、潜在構造を仮定することなしに順序回帰への道を開くという利点があります。メトリックを選ばないといけないというのは確かに難点ですが、順序説明変数を扱うときにだって同じ問題は起きているわけですし。

Tutz: 順序データ分析の研究では予測が無視されがちじゃないですか? 二値データとかカテゴリカルデータの場合には、分類の枠組みであんなに予測の研究があるのに。なんで判別分析とか機械学習で順序データが注目されてないんでしょうかね。

Liu & Agresti: そうですね。

Simonoff: あんまし触れてくださいませんでしたが、ノンパラメトリック・セミパラメトリックなアプローチについて補足させてください。私の専門なもので。
 まずノンパラメトリックなアプローチについて。大域的にパラメトリックな仮定をおくのをやめて、局所的な多項関係を仮定します。たとえば... [以下、説明省略。難しいよお]
 少なくとも局所的にはパラメトリックなモデルがあてはまると思う場合には、セミパラメトリックなアプローチがよいでしょう。たとえば... [説明省略]

Liu & Agresti: いやー、良く知らなかった話で、勉強になりました。[←とは書いてないけど、まあそういうニュアンス]

Kateri: コメントの前に、触れてくださらなかった話題について紹介してもいいですか? 正方順序行列の分析という話題です。
 分類変数の正方行列の場合、対称モデル、準対称モデル、周辺等質性モデルがあることは広く知られています。これが順序変数の場合だと、適切なモデルのクラスがさらに広がって、条件付き対称モデル、対角非対称モデルというのが出てきます。
 これらのモデルは、一致の測定とか、移動表の分析とかと関連していますし、行列がマルコフ連鎖の遷移確率行列である場合には「ランダム・ウォーク」「均衡状態」「リバーシビリティ」という概念と関わってきます。

Liu & Agresti: なるほど、ありがとうございます。順序変数の正方行列だったらほかにもこんなモデルがありますね...[とひとしきり紹介があるが、話についていけなくなったのでパス。このへんはAgrestiの教科書をみたほうがよさそう]

Kateri: クロス表の連関のモデル化についてコメントします。
 連関はKL距離の観点からの指標、相関はピアソン距離の観点からの指標です。Goodman(1996) はこの二つを含むより一般的なクラスのモデルを考えました... [以下、自分の研究紹介。省略]
 順序クロス表の分析における興味深い問題として、行ないし列の併合という問題があります... [自分の研究紹介]

Liu & Agresti: (特に返答無し)

Lesaffre: 私、自分が相談を受ける際は、順序ロジスティック回帰じゃなくて二値ロジスティック回帰を勧めることが多いんですよね。相手はたいてい医者で、二値ロジスティック回帰の出力ならよく理解できますし、二値ロジスティックでも順序ロジスティックでも結果はたいていほとんど変わらないから。おそらく他の統計家の方もそうなんじゃないかと。
 お伺いしたいんですけど、二値じゃなくて順序ロジスティックを使うことの利点を明確に示したシミュレーションとかないでしょうか? それにソフトはあるんでしょうか?

Liu & Agresti: Lesaffre先生のこのご指摘にはちょっとガックリなんですけど、ええとですね、順序モデルだとそれぞれの効果についてひとつのパラメータが出ます。二値に潰して分析するのとでは結果は結構違います。Whitehead(1993)をみてください。

Lesaffre: Anderson(1984)のモデルについて触れて下さらなかったのにはちょっとびっくりです。このモデルは評価の順序カテゴリ変数を分析するもので、
 $P(Y=s) = \exp(\beta_{0s} - \phi_s \beta^T x) / $(分子の合計)
ただし$1 =\phi_1 > \ldots > \phi_k = 0$とします。これは多項ロジット回帰の特殊ケースですが、説明変数との関係が一元的だという点が古典的な順序データ分析とは違います。
 この尺度上での評価者の誤分類をどう修正するのかというのは重要な研究分野だと思うんですけど。

Liu & Agresti: ご指摘の通り、見落としていた研究です。$\{\phi_t\}$が固定なら対数線形モデル、等間隔なら隣接カテゴリロジットモデルとして捉えられますね。このモデルがあまり注目されていない理由のひとつはパラメータ数が多いからだと思いますが、Andersonさんが早く亡くなってしまったせいもあるかもしれません。[←へー]
 誤分類の話ももっと注目されていいと思います。先生がレビュー論文書いてください![←おまえら実は仲良しだな]

Loughin: [ものすごく大雑把に要約すると] たとえば我々の研究だと、GLMMモデルでカウントデータを生成してGLMMモデルを当てはめると、完全に正しいモデルなのにパラメータが保守的に歪んだりします。いい加減だけど簡単なモデルじゃなくて、複雑で正しいモデルを使うべしという根拠とかありますかね? 分析者も忙しいわけで、証拠がないと移行しないと思うんです。

Liu & Agresti: 大事なご指摘と思います。たとえば、データがスパースなときは名義尺度モデルより順序モデルのほうが有利です。こういう風に、実務的な状況での手法の良さを調べていくことが大事ですね。

Svensson: 私は対応のある順序データをカテゴリ数と無関係に評価する方法について研究しています。視覚的アナログ尺度のデータなどに使える方法です。基本的なアイデアは、評価の間の依存性についての情報を得ることで説明できるaugmented rankingというのを考えるということです。変数ペアにランクを振るわけです。[...以下説明が続く。面白そうなんだけど、この説明だけではちょっと理解できそうにないので、メモ省略]

Liu & Agresti: [頑張ってください的コメント]

Svensson: [統計教育についてのコメント。パス]

Liu & Agresti: [Agrestiさんのフロリダ大での学部教育の話。パス]

Aguilera: 論文に出てきませんでしたが、PLS順序ロジスティック回帰というのがあります。これはPLS一般化回帰の特殊ケースで、変数の数が事例数より多かったり多重共線性があるときにも使えます。主成分ロジスティック回帰というのもあります。
 関数的データ分析(functional data analysis)というのもあります。これは観察がベクトルじゃなくて関数であるときに、関数的共変量の観点から反応変数を説明しようというものです[←??? よくわからん]。

Liu & Agresti: ご指摘ありがとうございます。いずれも興味深い問題ですね。

 ... はい!終了!疲れた!
 すぐに忘れてしまいそうなので、最後に「ここがわからん」「ここをもっと知りたい」という箇所をメモしておこう。

論文:データ解析(2018-) - 読了: Liu & Agresti (2005) 順序カテゴリカルデータの分析方法レビュー (全体的に遠慮がちな質疑応答つき)

Rosipal, R., Kramer, N. (2006) Overview and Recent Advances in Partial Least Squares. in Saunders C. et al. (eds) Subspace, Latent Structure and Feature Selection. Lecture Notes in Computer Science, vol 3940. 34-51.
 PLS回帰にはいろいろバージョンがあるようで、イライラするので読んでみた。先日読んだAbdiさんの解説にも引用されていた。

 いわく。
 PLSというのはPLS回帰だけでなく、もっと幅の広い手法である。それらの共通点は、観察データとの背後に少数の潜在変数が存在すると考える点である。ふつうは2つの変数セットについて(回帰の場合だと予測子セットと反応セット)、その間の共分散を最大化するようにして、直交する潜在変数スコアを求める。
 統計学者はもともとPLS回帰を無視する傾向があり、いまでもPLS回帰を統計モデルというよりアルゴリズムとして捉えていることが多い。もっとも最近では、PCRやリッジ回帰などとあわせ、continuum回帰という統一的なアプローチとして捉えられるようになった。
 PLSを分類手法としてみるとFDAと密接な関連がある。またPCAと関連付けることもできる。潜在変数をつかってSVMするという提案もある[へー]。カーネル法と併用することもできる。

 標本サイズを$n$とし、$N$変数データ$X$と$M$変数データ$Y$があるとしよう(どちらも変数の平均は0にしてある)。PLSではこれらを
 $X = TP^T + E$
 $Y = UQ^T + F$
と分解する。$T$と$U$は$p$個の潜在変数のスコア行列でどちらもサイズ$n \times p$。$P$と$Q$は負荷行列で、サイズは$N \times p$, $M \times p$。$E$と$F$は残差行列。
 オリジナルのNIPALSアルゴリズムでは、$Xw$と$Yc$の標本共分散を最大化するような、長さ1のウェイト・ベクトル$w, c$を求める。手順としては、まずランダムなスコアベクトル$u$を用意しておき、

  1. $w = X^T u/(u^T u)$
  2. $||w|| \to 1$ [←これがXウェイトね]
  3. $t = X w$ [←X得点]
  4. $c = Y^T t / (t^T t)$
  5. $||c|| \to 1$ [←Yウェイト]
  6. $u = Yc$ [←Y得点]

これを収束するまで繰り返す。なお、$Y$が1変数($y$)の場合には$u$が$y$そのものになるから、繰り返しはいらない。また、実はウェイトベクトル$w$は、
 $X^T Y Y^T X w = \lambda w$
という固有値問題の第一固有ベクトルになっている。
 以上が第一因子。$X, Y$をデフレーションし、今度は第二因子を求める。

 さて、PLSにはいろいろな形式がある。それらのちがいは、スコア$t, u$を求めたあとで、$X, Y$をどうやってデフレーションさせるかという点の違いである。4つの形式について述べる。
 以下、$p = X^T t / (t ^T t)$, $q = Y^T u / (u^T u)$を負荷ベクトルと呼ぶ。

 PCA, CCAとの関係について。
 PCAは$Xw$の分散を最大化するような$w$を求める。
 CCAは$Xw$と$Yc$の相関を最大化するような$w, c$を求める。これは正準リッジ分析という概念からみることもできる。いま、
 $\displaystyle max_{|r|=|s|=1} \frac{cov(Xr, Xs)}{([1-\gamma_x] var(X r) + \gamma_x)([1-\gamma_y] var(Y s) + \gamma_x)}$
という最適化問題があるとしよう。$\gamma_x, \gamma_y$は正則化項で0以上1以下。両方を0にするとCCAとなり、両方を1にするとPLSとなる。つまり、その中間(正則化CCA)もあるわけだ。また、$Y$が1次元であればリッジ回帰もこのなかにはいる[へー!]。
 なお、CCAは固有値分解で一発で求めるという点ではPLS-SBに近い。

 PLS1, PLS2についてもう少し詳しく解説すると... [OLS, PCR, リッジ回帰, PLS回帰の回帰係数を統一的に説明する。途中でついていけなくなったので略]
 回帰係数$b$の推定量$\hat{b}$のMSEをバイアス-バリアンス分解すると、OLS推定量はバイアス0で、そのかわり($X$の共分散行列の固有値が小さい時に)バリアンスが大きくなっている。そこで推定量を縮小して
 $\hat{b}_shr = \sum_i f(\lambda_i) \hat{b}_i$
とすることを考える。$\lambda_i$は$X$の第$i$特異値の二乗、$f(\lambda_i)$はシュリンケージ・ファクター。$i$は共分散行列のランクだけ回る。この枠組みで考えると、PCRとは「$i$が$p$以下の時$f(\lambda_i)=1$, そうでないとき0」に相当し、リッジ回帰は$f(\lambda_i) = \lambda_i / (\lambda_i + \gamma)$に相当する。
 実はこの縮小推定の枠組みでPLS回帰も説明できてしまう。PLSではなんと
 $f(\lambda_i) = 1 \prod_j (1-(\lambda_i/\mu_j))$
なのである。[$\mu_j$というのは... 上で読み飛ばしたくだりで説明されているんだけど、なにがなんだかさっぱりわからない。なんでもRitzペアというのがあって、それを理解するためにはKrylov空間というのとLanczosアルゴリズムというのを理解する必要があるらしい... 知らんがな!!!]

 PLSを判別・分類に使う場合は... [関心はあるんだけど、疲れちゃったのでパス]
 非線形PLSとは... [カーネルPLSとか。疲れちゃったのでパス]

... というわけで、中盤でPLS回帰をいろんな線形回帰の推定量のひとつとして位置付けるあたりで頭がパンクしてしまい、肝心の後半(判別・分類や非線形PLS)は断念。また体力のあるときに。
 PLS回帰の解説ってすぐに反復アルゴリズムの話になっちゃうので、そういうものかと思っていたんだけど、ちゃんとこういう統計学的な議論もあるのね。正準相関分析の正則化として捉えるとか、線形回帰の縮小推定量として捉えるとか。知りませんでした。
 昔からPLS回帰とPLS-SEMとの関係がいまいちわからなくて困っているんだけど、PLS-SEMについては言及がなかった。ま、言及されても理解できないかもしれないけれどさ。

 先日読んだAbdiさんの解説と見比べると、説明がちょっと違っていて混乱する... この論文で言うPLS-SBとはAbdiさんのいうBPLSのことだろうか? Woldのオリジナルの提案は、結局$X$と$Y$を対称的に扱っていたのか(この論文)、そうでないのか(Abdiさん)?

論文:データ解析(2018-) - 読了: Rosipal & Kramer (2006) PLSレビュー in 2006

Abdi, H. (2010) Partial least squares regression and projection on latent structure regression (PLS regression). WIREs Computational Statistics. 2(1), 97-106.
 「PLS回帰というのはこういうもので...」と人に説明する必要が生じたのだが、めんどくさいのでよさそうなPDFを送り付けて済ませてしまった。ごめんなさい。渡しておいて読んでないというのもなんなので、自分でも読んでみた次第。
 著者AbdiさんはPLS回帰の専門家として私のなかで評価が高い方である(超・上から目線。要はわかりやすく書いてくれるのでスキよということです)。いま検索してみたら、なかなか渋い感じのおじさんであった。

 途中で気が付いたんだけど、これ、前半は以前読んだAbdi(2007)とほぼ同じ。どうせ三歩歩けば忘れてしまうので、気にせずメモを取りましたけど。

 いわく。
 PLS回帰は、予測子がたくさんあるときの予測に向いている手法。もともとPLS回帰は社会科学から生まれたのだが(Woldは経済学者だしね)、ケモメトリクスとか官能評価とかで人気がでて、社会科学に逆輸入された。

 事例数を$I$、従属変数の$I \times K$行列を$Y$, 予測子の$I \times J$行列を$X$とする。
 $X$で$Y$を予測したい。$Y$が単変量ならふつうはOLSを使うところだが、$J$が$I$より大きかったりすると$X$が正則でなくなる。さあどうするか。変数を削るか? リッジ回帰に走るか?
 ひとつのアイデアは主成分回帰(PCR)である。すなわち、$X$を主成分分析し、主成分を予測子とする。テクニカルには、$X$をSVDにかけて
 $X = R \Delta V^T$
右特異行列$V$を負荷と呼ぶ。$G = R\Delta$を因子得点と呼ぶ。で、$G$を予測子にしてOLSをやる。
 このアイデアの短所は、予測子から最適なサブセットを選ぶという問題が未解決だという点である。上位主成分を選ぶのは、$X$の説明という意味では最適だが、$Y$の説明という意味では最適でない。

 そこでPLS回帰でございます。
 $X$と$Y$の両方を、共通の直交因子の積として表現することを考える。つまり、まず
 $X = TP^T$
と分解する。$T$の列を「潜在ベクトル」と呼ぶ。ふつうは$T^T = I$とする(実のところ、ここの規準化の仕方によって微妙な差が出てくる。ソフト間で結果が違う原因がここにある)。PCAとのアナロジーで言うと、$T$が得点行列、$P$が負荷行列である。
 $Y$は
 $\hat{Y} = T B C^T$
と推定する。$C$は「ウェイト行列」, $B$は対角行列で「回帰ウェイト行列」という。

 潜在ベクトルの行列$T$をどう決めるか。PLS回帰ではこれを、$X$の線形和と$Y$の線形和を最大化する問題として捉える。
 オリジナルのアルゴリズム(NIPALS)だとこうである。$X$, $Y$はあらかじめ標準化しておく。これを$E$, $F$に代入する。ランダムベクトル$u$を用意しておく。
 ステップ1. $E^T u$を求め、標準化して「Xウェイト」$w$とする。
 ステップ2. $Ew$を求め、標準化して「X得点」$t$とする。
 ステップ3. $F^T t$を求め、標準化して「Yウェイト」$c$とする。
 ステップ4. $u=Fc$を「Y得点」とする。
これを収束するまで繰り返す。で、係数$b = t^T u$と、「X負荷」$p=E^T t$を求める。$E$の予測値$tp^T$, $F$の予測値$b t c^T$を求め、$E, F$から引く。
 以上が第一因子。これを繰り返していく。$b$は$B$の対角にいれていく。
 $X$の二乗和のうち潜在ベクトルによって説明された部分は$p^T p$、$Y$の二乗和のうち潜在ベクトルによって説明された部分は$b^2$となる。

 ところで、行列の固有値を求める方法にpower methodというのがあるんだけど、NIPALSはこれと似ている。
 NIPALSはSVDの観点からも捉えることができる。上のをつなげると、「求めて標準化する」のを$\propto$と書くとして、
 $w \propto E^T u \propto E^T Fx \propto E^T F F^T t \propto E^T F F^T E w$
なので、
 $S = E^T F$
をSVDしたときの第1右特異ベクトルが$w$, 第1左特異ベクトルが$c$であることがわかる。[←ううううむ... ごめんなさい、私にはよくわかんないです... ]

 PLS回帰モデルをどう評価するか。
 手元のデータについての記述が目的なのであれば、モデルの評価にもっともよく使われているのはRESS(残差平方和)である。$L$因子PLS回帰による予測値の行列を$\hat{Y}^{[L]}$として、
 $RESS = || Y - \hat{Y}^{[L]} || ^2$
が小さければ良い予測だということになる($L$を大きくすれば小さくなる)。潜在変数については、潜在変数を実質的に解釈することに意味がある文脈ならば、解釈すればよろしい。
 そうではなくて、手元の標本に基づく新データの予測が目的だとすると、それぞれの潜在変数による寄与と、潜在変数の数の決定が問題になる。これはジャックナイフ法やCVで調べる必要がある。相関係数の二乗とか、PRESS(予測誤差平方和)とかが良く使われている。予測値の行列を$\tilde{Y}^{[L]}$として、
 $PRESS = || Y - \tilde{Y}^{[L]} || ^2$
である。
 $L$をどう決めるか。記述の場合とちがい、PRESSは$L$を大きくし過ぎると逆に上がってしまう。ひとつの方法は、PRESS最小の$L$を選ぶというもの。もうひとつは、$L=l$のときのPRESSとRESSを$PRESS_l$, $RESS_l$として、
 $Q^2_l = 1 - (PRESS_l) / (RESS_{l-1})$
が閾値を上回っている最大の$l$を選ぶというもの($RESS_0=K(I-1)$とする)。閾値の案としては、$(1-0.95^2)=0.0975$とか、「事例数が100以下だったら$.05$、そうでなかったら$0$」とか。[このくだり、理解できない。挙げられている事例で計算してみたんだけど、$Q^2$の値が全然合わない。なにか読み間違えているのか、誤植があるのだと思う]

 ところで、ここまでで述べてきたPLS回帰はWold, Martensから発する系統。もうひとつ、Booksteinという人の系統もある。BPLSと呼ぶことにしよう。
 BPLSでは、$S = X^T Y$のSVDが
 $S = W \Theta C^T$
だとして、
 $T = XW$
 $U = YC$
とする。BPLSの場合、$X$と$Y$を入れ替えても同じ結果になる。正準相関分析とすごく似ている(ちがいは、正準相関分析では$X$, $Y$の線形和の相関を最大化するのに対して、BPLSでは共分散を最大化するという点)。
 BPLSのバリエーションとして、

 云々。

 。。。一番勉強になったのは、PLSとBPLSのちがい。そういうことだったのか、とようやく腑に落ちた。以前、PLS回帰についての説明を複数種類読み比べて、微妙なちがいに頭を抱えたことがあるんだけど、以前取ったメモと見比べて思うに、Rのplsパッケージの解説書Mevik & Wehren(2007)で説明されているのは、PLSではなくてBPLSなのではないかと思う。

論文:データ解析(2018-) - 読了:Abdi (2010) PLS回帰への招待 in 2010

Indahl, U.G, Liland, K.H., Naes, T. (2008) Canonical partial least squares: A unified PLS approach to classification and regression problems. J. Chemometrics, 23, 495-504.
 PLS回帰と正準相関分析ってどう違うのだろうか... と検索していて見つけた論文。Rのplsパッケージにはcanonical powered PLSという手法が載っていて、その参照文献のひとつになっている。あ、第二著者のLilindさんという人はplsパッケージの作者のひとりだ。
 掲載誌が掲載誌だけに(ケモメトリクスはPLS回帰の本場)、PLS回帰そのものについての知識を前提として、新手法を提案している論文。

 いわく。
 本論文ではPLS2 (多重反応PLSのことらしい) よりももっと自然なやりかたで、PLSによる判別分析を提供します。その際、すでに提案されているpowered PLS(PPLS)という方法を使います。
 それではご紹介しましょう。じゃじゃーん。その名も正準PLS(CPLS)。PLSと正準相関分析(CCA)をあわせた手法です!

 背景。
 予測子$n \times p$行列をX, 反応の$n \times q$行列を$Y$とする。
 $Y$が連続変数である場合に、$p$次元単位ベクトル$u$, $q$次元単位ベクトル$v$について
 $f_1 (u, v) = u^t X^t Y v$
を最大化する手法をPLS2という。
 これはSVDで解ける。$X, Y$を中心化しておいて、$W = X^t Y$をの一階SVD近似を$W_{(1)} = s a b^t$とする。このとき、$cov(Xu, Yv) = (1/n) f_1(u, v)$となる。[←あれれ...どこかに誤植があるはずだ。$cov(Xa, Yb) = (1/n) f_1 (u,v)$の間違い??]

 $Y$が$n \times g$のダミー行列である場合、ちょっとやり方を変えたほうがよい。$W = X^t Y$のSVDから特異ベクトルを得るのは、
 $B_Q = n \bar{X}_g^t Q \bar{X}_g$
のSVDをするのと同じことになる。ここで$\bar{X}_g = (Y^t Y)^{-1} Y^t X$、すなわち$g \times p$の群平均行列。$Q$は重みを表す対角行列で、対角要素は群サイズ$n_k^2$に比例し合計1となる。それよか、
 $f_2(u, v) = u^t X^t Y (Y^t Y)^{(1/2)} = u^t W_\Delta v$
を最大化することを考えたほうがいい。なぜなら、$W_\Delta$のSVDから特異ベクトルを得るのは
 $B_P = n \bar{X}_g^t $P$ \bar{X}_g$
のSVDと同じことになり、ここでPの対角要素は$n_k$に比例し合計1だ。これをPLS-DAという。[はああ? どういうこと? 狐につままれた気分]

 今度はCCA。
 PLS2では$cov(Xu, Yv)$の最大化を行ったわけだが、今度は$corr(Xu, Yv)$の最大化を行う。これは結局... [以下、説明が続くけどパス]

 で、ここからが新手法CPLSの説明なんだけど、哀しいかな、理解できなかった。なんでも、$corr(Xu, Yv)$じゃなくて$corr(X X^t Y c, Yd)$を最大化するのだそうだ。えっ、なぜ???
 説明が書いてあるんだけど、これがもう全然理解できない... 涙を呑んで諦めます。今回はご縁がなかったということで。
 えーと、とにかく、PLSやPLS-DAと比べて、よりアグレッシブに$Y$を予測し、PLS因子数は少なくなるんだそうだ。

 。。。いやいや。深刻に考えるのはよそう。ある日ふと読み返したら、意外にスルスルと理解できるかもしれない。あるいは加齢によってこの分析手法への関心ごと消え失せ、理解できなくてもなんら痛痒を感じなくなるかもしれない。いいよ、もう、どうだって。

論文:データ解析(2018-) - 読了:Indahl, Liland, Naes(2008) 正準PLS回帰分析のご提案 (おまえらにわかるかどうかは別にして)

 どういうわけか、マーケティングのための消費者調査で糊口を凌いでいるのだが、いろんな調査課題があるなかのひとつに、価格調査と呼ばれるものがある。最初は市場における価格についての実態調査のことかと思ったんだけど、そうではなくて、企業が価格について意思決定する際にそれを支援するための消費者調査、という意味である。いくら値下げしたらなんぼ売れまっせ、というような。
 市場調査についての解説も、価格戦略についての解説もあまたあるけれど、価格調査についての詳しい解説は意外に少ない。入手可能な解説書として思いつくのは、朝野・山中「新製品開発」など朝野先生の一連の書籍か、杉田・上田・守口「プライシング・サイエンス」くらいである(他にあったらぜひ教えてほしいです)。
 仕方がないので英語の本に頼るんだけど、読んでいるとどうしても眠くなってしまう。仕方がないのでメモをとったりして... せっかくなので整理・追記して載せておくことにする。と、このように、私の人生には「どういうわけか」「仕方がないので」「どうしても」「せっかくなので」といった形容語が多い。

 以下、Monroe(2003) "Pricing: Making Profitable Decisions", Third Edition, Chapter 9. Research Methods for Pricing Decisions からのメモ。大学の教科書だと思う。こういう教科書があるなんて、アメリカの高等教育はさすがだと感心すべきなのか、資本主義の支配ここに極まれりと呆れるべきなのか、よくわからない。私が持っているのは2005年のInternational Editionで、表紙がぺらっぺらで持ちにくい。

1. 価格決定における問い。以下が挙げられる。

2. 価格調査のポイント。次の3つの点を決める必要がある。

3. 価格調査のアプローチ。大きく分けて4つある。

4. 具体的な調査手法について。とにかくポイントは、調査対象者に適切な参照枠を提示することである(顧客が持っている参照枠の違いが反応に及ぼす効果を調べるというのなら別だけど)。

 その1、受容価格の上下限を調べる方法。

 その2、購入意向を知る方法。
 80年代のパーカー・ペンの実例が挙げられている(元ネタはTomkovick&Dobie(1995, JPIM))。低価格帯への参入にあたり、次の5つのステップを踏んだ。

  1. ターゲット顧客を特定。
  2. 重要属性を特定。まずConsumer Reportsから属性を洗い出した。消費者調査で差別化ポイントは "feel" と "control"だと突き止めた。で、前者は重さ、太さ、形などで決まり、後者は紙の上を滑るスムーズさと疲労感で決まるのだろうということになった。[←どうやって決めたんだろう]
  3. 市場にある製品について、小売価格を目的変数、上記属性を説明変数にとった回帰分析。全属性を平均にすると価格は1.07ドルだと推定。[←そうか、ヘドニック価格モデルってことか]
  4. 再び消費者調査。全属性を平均にしたペンを実際につくり、参照ペンとして渡し、「専門家パネルによればこれは1.07ドルの価値があるそうです」と教示。で、「これより10パーセント増しでスムーズに書けるペンがあったらいくら払う?」と聞いていく。[←まじか... そう訊かれても困るよな]
  5. このデータを回帰分析。

というわけでパーカーは、参照ペンよりちょっと重くてちょっと書きやすいペンが最適なんだけど、いずれにせよ支払意思額はすごく低いということを発見し、安っすいペンを売り出して大儲けしましたとさ。[←嗚呼... わざわざメモして損した...]

 その3、価格感受性を推定する方法。既存製品ラインに対する製品追加や価格変更の場合、受容価格帯やある価格に対する需要が分かっただけでは困るわけで、価格弾力性や交差価格弾力性を知る必要がある。

 ... というわけで、頭の整理になりましたです。
 コンジョイント分析に関しては首をひねる記述が多かった。ちょっと内容が古いような気がする。

雑記 - 覚え書き:価格調査の方法

赤穂昭太郎 (2013) 正準相関分析入門:複数種類の観測からの共通情報抽出法. 日本神経回路学会誌, 20(2), 62-72.

 仕事のなかで正準相関みたいな問題が出てきちゃったので、勉強の足しになるかと思って読んだ。こういうの、昔から苦手分野である。まあ得意分野なんてどこにもないけどな。

 あれこれ「へーっ」という話がありましたが、ひとつだけメモ。
 ある測定で$m$次元ベクトル$x$と$n$次元ベクトル$y$のペアが観測されるとする(いずれも標本平均を0とする)。それぞれを線形変換して得られる
 $u(x) = a^{T}x$
 $v(x) = b^{T}y$
について、その相関を最大化したい。これは一般化固有値問題に帰着する(要は解けるってことですね先生)。最大で$min(m,n)$個の成分がとれる。これが線形正準相関分析。

 さて、$u(x), v(x)$を非線形化したい。
 カーネル関数を$k_x(x, x')$とする。なんでもいいけど、ガウスカーネルなら
 $k_x(x, x') = \exp(-||x-x'||^2/2\sigma^2)$
である。ある観測ベクトル$x_i$について
 $\phi_x(x_i) = (k_x(x_i, x_1), k_x(x_i, x_2), \ldots, k_x(x_i, x_n))^T$
を求める。同様に$\phi_y(y_i)$も求める。それぞれから平均を引き、それを$x, y$だと思ってやおら線形正準相関分析をやる。実際にはオーバーフィットしちゃうので正則化項をいれる(残念ながら、ここで人手によるパラメータ調整が必要になる)。
 $u(x) = a^{T} \phi_x(x)$
 $v(x) = b^{T} \phi_y(y)$
というわけで、非線形な特徴$u(x), v(y)$が手に入ったことになる。

 ところが。これとは別の非線形化アプローチがある。
 実は、$(x, y)$の同時確率密度関数$p(x, y)$が既知ならば、$u(x), v(x)$は連立方程式
 $\int \{p(x|y)-p(x)\} u(x) dx = \lambda v(y)$
 $\int \{p(y|x)-p(y)\} v(x) dx = \lambda u(y)$
の解になるのだそうだ。実際には$p(x, y)$は未知なので、実用性は低いんだけど、著者の先生いわく、「xやyが実数ベクトルではなくアンケートデータのような質的データの場合の解析法である数量化III類は、この定式化から理解することが容易となる」とのこと。へええ、そうなんですか。

 。。。いやー、もうね、嫌になっちゃいますね、勉強しないといけないことが多すぎて。
 そういえば、正準相関とPLS回帰ってどう違うんだろう。相関を最大化するか、共分散を最大化するかという違いに過ぎないのだろうか? 疲れちゃったのでいまは考えたくないけど。

論文:データ解析(2018-) - 読了:赤穂(2013) 正準相関分析入門

Jain, A.K., Murty, M.N., Flynn, P.J. (1999) Data Clustering: A Review. ACM Computing Surveys, 31(3), 264-323.
クラスタリングについてのレビュー論文。Google scholar上の被引用件数が14000件超というすごい奴である。勉強のために読んだ。20年近く前の論文を、なにをいまさらという気もするけれど...

1. イントロダクション
1.1 モチベーション
 探索的データ分析であれ検証的データ分析であれ、その手続きにおいて鍵となる要素のひとつにグルーピングがある。パターンの集まりを類似性に基づくクラスタへと組織化することをクラスタ分析という。判別分析(教師つき分類)とは異なり、与えられたパターンにはラベルがついていない(教師なし分類)。いいかえれば、カテゴリのラベルはデータ駆動的である。
 本論文ではクラスタ分析におけるコア概念と技法をサーベイする。統計学・決定理論にルーツを持つ分野を中心とするが、機械学習などのコミュニティに関しても言及する。

1.2 クラスタリング課題の構成要素
 ふつう、クラスタ分析は次の段階からなる:

  1. パターンの表現。3章で扱う。
  2. 近接性の定義。4章で扱う。
  3. クラスタリング。5章で扱う。
  4. (必要ならば)データ抽象化。たとえば、クラスタについての記述とか。
  5. (必要ならば)出力の評価。ふたつにわかれる。
    • cluster tendencyの分析。そもそもクラスタ分析していいことがあるのかどうかを調べる。研究分野としてあまり活発でないので本論文では割愛。
    • cluster validityの分析。得られたクラスタが偶然の産物だったりアルゴリズムのアーティファクトだったりしないかを調べる。(1)外的評価(アプリオリに与えられた構造を復元できるか)、(2)内的評価(クラスタがデータに対して内在的に適切か)、(3)相対的テスト(2つの構造を比べてどちらが良いか)、に分けられる。

1.3 ユーザのジレンマと熟達の役割
 手元の問題に適したアルゴリズムをどうやって選ぶか。Dubes & Jain (1976)はFisher & Van Ness (1971)が定義した許容性基準によってアルゴリズムを比較している。そこでの許容性基準のベースにあるのは、(1)クラスタの形成方法、(2)データの構造、データの構造に影響しない変化への敏感性、である。でも、「データは規準化すべきでしょうか」「どんな類似性指標がどんな状況で適切でしょうか」といった問いには答えてくれない。[←ああなるほど、アルゴリズムのadmissibilityの比較だけでは役に立たんということか]
 多次元データセットにおける構造をいつでも明らかにしてくれるという便利な方法はない。ユーザは個々の手法について十分に理解するだけでなく、データ収集過程の詳細についての知識と、ある程度の領域知識を持たねばならない。
 データソースに対する適切な制約もクラスタリング手続きの一部である。たとえばmixture resolvingがそうだ。データが未知の数の確率密度の混合から抽出されていると仮定するせいで、混合要素の数と各要素のパラメータを推測することがクラスタリング問題となる。

1.4 歴史(略)
1.5 本論文の構成(略)

2. 定義と表記
 あるパターン(特徴ベクトル、観察事例、データ)を$\mathbf{x}$とする。その要素$x_i$を特徴(属性)と呼ぶ。パターンの次元数を$d$とする。
 パターン集合を$\mathcal{X} = \{\mathbf{x}_1, \ldots, \mathbf{x}_n\}$とする。通常は$n \times d$行列である。
 パターン生成過程を支配する状態ををクラスと呼ぶ。クラスタリングとは、パターンをグルーピングして、異なるパターン生成過程を反映するクラスを手に入れようとすることである。
 パターン集合$\mathcal{X}$にラベル$\{l_1, \ldots, l_n\}$を付与することをハード・クラスタリングという。ラベルの種類数(クラスタの数)を$k$とする。
 パターン$\mathbf{x}_i$に、それがクラスタ$j$に属する程度$f_{ij}$を付与することをファジー・クラスタリングという。
 パターンの類似性を定量化するのに用いる、特徴空間におけるメトリック(ないし準メトリック)を距離指標という。

3. パターン表現、特徴選択、特徴抽出
 [良いパターン表現がいかに大事か、という説明が続いて...]
 特徴は、量的(連続, 離散, 間隔)、質的(名義, 順序)にわけられる。特徴が木として構造化されている場合もある。
 特徴選択は大事。本質的にアドホックで試行錯誤を要する。
 PCAなどで特徴抽出することもある。
 [この節、期待してたんだけど、ごく短い。がっくり。まあそういう主旨のレビュー論文ではないのだろう]

4. 類似性指標
 連続的特徴の場合、一番一般的なメトリックはユークリッド距離である。
 $d_2(\mathbf{x}_i, \mathbf{x}_j) = (\sum_k (x_{ik}-x_{jk})^2)^{1/2} = || \mathbf{x}_i - \mathbf{x}_j||_2$
 一般化するとミンコフスキー距離。
 $d_p(\mathbf{x}_i, \mathbf{x}_j) = (\sum_k (x_{ik}-x_{jk})^p)^{1/p} = || \mathbf{x}_i - \mathbf{x}_j||_p$
 良い点:(1)直感的にわかりやすい。(2)「コンパクト」ないし「よく分離されている」クラスならうまく機能する。
 悪い点:(1)スケールが大きい特徴が支配的になる(なので、ふつうは事前に規準化する)。(2)特徴間に線形相関があると歪む。そういうときは正規化マハラノビス距離
 $d_M(\mathbf{x}_i, \mathbf{x}_j) = (\mathbf{x}_i - \mathbf{x}_j) \Sigma^{-1} (\mathbf{x}_i - \mathbf{x}_j)^T$
を使うこともある。$\Sigma$は既知の共分散行列または標本共分散行列。これは、各クラスの条件付き密度が単峰で、多次元空間における広がりで特徴づけられるという仮定(つまり多変量正規だという仮定)に基づいている。Mao & Jain (1996 IEEE Trans.Neur.Net.)をみよ。ほかにHausdorff距離というのを使うという提案もある。 Huttenlocker, et al. (1993 IEEE Trans.PatternAnal.Mach.Intell.)をみよ。

 クラスタリングアルゴリズムはふつうパターン集合を使い、パターン間の距離を求めるのだが、最初から距離行列を使うアルゴリズムもある。
 連続量でない特徴から距離を求めるというのは厄介な問題である。Wilson & Martinez (1997 J.Artif.Intell.Res.), Ichino & Yaguchi (1994 IEEE Trans.PatternAnal.Mach.Intell.)をみよ。
 パターンを文字列で表現するというのもある(syntactic clusteringという)。統計的手法に比べるとすべての面で劣っているといわれているので、本稿では割愛。パターンを木構造で表すというのもある。
 近隣のデータ点(文脈)の効果を考慮した距離指標もある。たとえば相互近隣距離(MND)。データ点$i$に関して$j$が持っている近隣数を$NN( \mathbf{x}_i, \mathbf{x}_j)$とする。[←よくわからないのだが、$j$からみて$i$より近い点の数、ってことかな?] で、
 $MND(\mathbf{x}_i, \mathbf{x}_j)  NN( \mathbf{x}_i, \mathbf{x}_j ) + NN(\mathbf{x}_j, \mathbf{x}_i)$
この距離は三角不等性を満たさないのでメトリックでないんだけど、クラスタリングでは良く使われている。このように、距離・類似性は必ずしもメトリックでない。

 特徴が十分にたくさんあれば、いかなるふたつのパターンであっても類似性は同程度である。類似性に差をつけるためにはなんらかの付加的な領域情報が必要になる。[ここでconceptual clusteringの例を紹介。3つの事例があって、ユークリッド距離では同程度なんだけど、うちふたつは同一のカテゴリに属しているようにみえて...という話。Michalski & Stepp (1983) というのが引用されているのだけれど、おそらくMichalski, Stepp, & Diday (1983 IEEE Trans. Pattern Anal. Mach. Intell.) の誤りだろう]

5. クラスタリング技法
 クラスタリングはとりあえず次のように分類できる。

といっても、ほかにもたくさんの分類がありうる。

5.1 階層クラスタリング・アルゴリズム
 [基本的な考え方と、デンドログラムについての説明があって...]
 たいていの階層アルゴリズムは、{単一リンク、完全リンク、最小分散}の変形である。[←書いてないけど、Ward法は最小分散アルゴリズムに分類されるのであろう]
 一番有名なのは、クラスタ間の距離を、それぞれ1個ずつ事例を出してその距離を調べた時の最小値(単一リンク)ないし最大値(完全リンク)とするやりかた。単一リンクだとchainingしやすい。完全リンクはコンパクトなクラスタになりやすい。実務的にもふつう完全リンクのほうが優れている。
 階層アルゴリズムは分割(Partitional)アルゴリズムに比べて使い勝手がよい面がある。たとえばk-meansはデータセットが等方的(isotropic)なクラスタを含んでいないとうまくいかない[←なるほどー]。いっぽう、アルゴリズムは階層アルゴリズムのほうが複雑。両者のハイブリッド案もある。
 なお、ここまでの話は全部hierarchicalかつagglomerativeだったけど、hierarchicalかつdivisiveなアルゴリズムもある。

5.2 分割アルゴリズム
 この方法の問題点は、クラスタ数を決めなきゃいけないという点。Dubes(1987, Pattern Recogn.)を参考にせよ。
 この方法はふつう、なんらかの関数を最適化するクラスタをつくるわけなんだけど、網羅的に探すのは大変なので、初期値を変えながら複数回繰り返し、最良の奴を選ぶことが多い。

5.2.1 二乗誤差アルゴリズム
 一番よく使われていて、直感的にもわかりやすい。二乗誤差とは、クラスタ$j$に属する$i$番目のパターンを$x_i^{(j)}$, クラスタ$j$の重心を$c_j$として、
 $\sum_j \sum_i || x_i^{(j)} - c_j ||^2$
のこと。
 もっとも良く使われているのがk-meansアルゴリズム。問題点は、初期値に対して敏感で、局所最小解に落ちちゃうかもしれない点。[図で説明している。親切...]
k-meansのバリエーションはたくさんある。(1)初期値選択を工夫する。(2)重心の距離が閾値を下回る2クラスタを併合する。例としてISODATAアルゴリズムがある。(3)二乗誤差以外の基準関数を併用する。例としてdynamic clusteringというのがある。

5.2.2 グラフ理論的クラスタリング
 いちばんよく知られているのは、データの最小全域木(minimal spanning tree; MST)をつくり、一番長いエッジを消すという手法。
 階層クラスタリングはグラフ理論の観点からみることもできて...[略]
 非階層、overlapありという手法でドロネーグラフ(Delaunay graph; DG)というのもある。ボロノイ近隣の点をつなぐ。

5.3 mixture-resolvingアルゴリズムとmode-seekingアルゴリズム
 mixture-resolvingアルゴリズムは、パターンが複数の分布から抽出されたものだと仮定し、それぞれの分布のパラメータを推定するアルゴリズム。ふつう分布はガウシアンだと仮定される。伝統的にはパラメータを最尤推定する。最近ではEMアルゴリズムを用いることが多い。[←とこのように、有限混合モデルについては意外にそっけない説明である。機械学習系の論文だからかしらん?]
 密度ベースのクラスタリングのノンパラメトリックな手法も提案されている。ノンパラ密度推定の方法にParzenウィンドウというのがあるけど、あれと同じようにして、入力パターンセットの多次元ヒストグラムのなかで、カウント数が大きくなるビンを探す。[←ああそうか、これってmodeを探索しているわけか...]
 ほかに、ノンパラ密度推定に基づく距離を使った階層手法・分割手法もある。

5.4 最近隣クラスタリング
 最近隣距離を使ったクラスタリングもある(Lu & Fu, 1978 IEEE Trans. Syst. Man Cybern.)。距離が閾値を下回る最近隣ペアを探しては、ラベルを一方から他方にコピーするというのを繰り返す。[←ちょ、ちょっと待って、最初のラベルはどうすんの?]

5.5 fuzzyクラスタリング
 対象を$N$個、クラスタを$K$個とする。まず$N \times K$の初期メンバーシップ行列$U$をつくる。要素$u_{ij}$は対象$x_i$のクラスタ$c_j$へのメンバーシップを表していて、ふつう$u_{ij} \in [0,1]$。
 で、$U$からファジー基準関数の値を求める。たとえば
 $\sum_i \sum_k u_{ij} ||x_i - c_k||^2$
というような関数。で、この基準関数が小さくなるように$U$を少しずつ変えていく。
 もっともよく用いられているのはfuzzy c-meanアルゴリズム(FCM)。他にもいろいろある。

5.6 クラスタの表現
 クラスタリングの最終目的が、データの分割ではなく、むしろデータの抽象化だということもある。クラスタの表現という問題はこれまであまり研究されていない。表現のスキーマとして以下が挙げられる。

 意思決定におけるデータ抽象化は重要である。なぜなら:

5.7 人工ニューラルネットワークによるクラスタリング
以下の特徴がある:

良く使われているのは、Kohonenの学習ベクトル量子化(LVQ)、自己組織化マップ(SOM)、適応的レゾナンス理論モデル(ART)。これらはアーキテクチャとしては単層。入力ノードと出力ノードのあいだのウェイトを更新する手続きは、実は古典的なクラスタリングとよく似ている。たとえば、LVQはk-meansと深いつながりがある。
 SOMは2次元のマップを使うので直感的にわかりやすいけど、初期ウェイトをしくじると局所解に落ちるし、収束に関わるパラメータがいっぱいあって、同じ入力なのに反復回数によっては異なる出力となりうる(つまり安定性がない)。学習率を下げていけば安定性が生じるけどこんどは可塑性がなくなる。
 ARTは安定的かつ可塑的といわれているんだけど、パターン入力の順序に依存する。また、クラスタの数とサイズがビジランス閾値と呼ばれる定数に依存する。
 SOMもARTも超球になっているクラスタの検出にしか向いていない。これに対して、2層のネットワークで楕円もみつけようという提案がある。

5.8 進化アプローチによるクラスタリング
 まず、クラスタリングの解(データの分割)のランダムな母集団をつくる。で、それぞれの解について適合度を求める(通常は二乗誤差に反比例する値)。で、それぞれの解を染色体と見立てて、selection, recombination, mutation操作を行い、第二世代の母集団をつくる。これを繰り返す。
 有名な手法として、遺伝的アルゴリズム(GA)、進化方略(ES)、進化プログラミング(EP)がある。GAが一番良く使われている。GAでは解は二値文字列になる。recombination操作にはいろいろあるけど、交差が一番良く使われている。selection, recombination, mutationのいずれも、適合度に基づいて確率的に行われる。
 ES, EPとGAとのちがいは解の表現と操作である(EPはrecombination操作を行わない)。要するにいずれも二乗誤差を最小化しているわけである。
 たとえばk-meansの場合、「次の反復」で得られる解は今ある解の近傍である。fuzzyクラスタリングもANNもそうだ。こういうのを局所化探索と呼ぶことができる。これに対して、GAでは交差や変異の操作によってとんでもない解が得られるわけで、大域化探索であるといえる。
 進化的アプローチによるクラスタリングにはいろいろ拡張案があって...[省略]
 GAの問題点としては、母集団サイズや交差・変異確率といったいろいろなパラメータに対する敏感性の高さが挙げられる。問題特有的なヒューリスティクスを組み込んだほうがよいという意見もある。
 クラスタリングを最適な分割を見つけるという課題ではなくて、最適な重心の位置を決めるという最適化問題として捉えるなら、ESやEPが使いやすい(GAとちがって、重心を実数ベクトルとして直接符号化できるから)。

5.9 検索ベースのアプローチ
 検索テクニックを使って最適な分割を検索するというアプローチ。ここでいう検索は決定論的検索と確率的検索に分かれる。
 決定論的検索としては、Branch-and-bound法で最適な分割を探すというアプローチがある。計算が大変。また、決定論的アニーリングを使うという手もある。大域最適解が保証されないというのが欠点。[←なるほどね... データ生成メカニズムは一切問わず、事例にクラスタ番号を付与する整数計画として解いちゃうということだろうか]
 確率的アプローチとしては、シミュレーテッド・アニーリング(SA)とか、Tabuサーチとかがある。[説明略]

5.10 手法の比較
各手法の特徴を概観すると、

 実証的な比較研究はいろいろあって...[めんどくさいのでパス]...
 要約すると、大データ向きなのはk-meansとKohonen net、あとは小データ用。どの手法でも領域知識をうまく使うことが大事。基準関数を使うアプローチだとクラスタは超球になりやすいので、それが困る場合は、階層アルゴリズムのほうがいいかも。

5.11 領域に基づく制約をクラスタリングに統合する
 そもそもクラスタリングは本質的に主観的な課題である。だから主観性をうまく組み込みたい。
 どんなクラスタリング・アルゴリズムであっても、なんらかのタイプの知識をつかっている。ただそれが暗黙的か明示的かがちがうだけだ。[←おおお、名文句だ...」
 パターンを表現するスキーマの選択、類似性指標の選択、グルーピングのスキーマの選択、ANN, GA, TS, SAにおけるパラメータの選択、いずれにおいても、知識が暗黙的に影響している。もっと明示的に知識を使うこともできる。知識を活かした特徴選択とか。
 領域知識の組み込みはアドホックなものなので、事例と先行研究を挙げることしかできない。こういうのに関心がある人は機械学習系の学術誌を読みなさい。

 実務的には、次の3つの問題が重要。(1)領域知識の表現、利用可能性、完全性。(2)知識を用いた推論の構築。(3)知識の変化への適応。

5.12 巨大なデータセットのクラスタリング
 大データに対して使われているのはk-meansとKohonen netである。超球じゃないクラスタが作れるという点では階層アルゴリズムが便利なんだけど、時間とメモリが大変になる。
 大データ向き手法もいろいろ開発されている。CLARANS, BIRCHが有名。[いろいろ説明が書いてあるけどパス]
 データが大きいとメモリにデータが載りきらなくなる。この問題に対する対処策として以下が挙げられる。

6. 応用
[画像分割、物体と文字の認知、書籍を検索するシステムというような状況検索、原油を掘るというようなデータマイニングでの活用例、を消化。内容も古いと思うし、気力も尽きたのでパス]

7. 要約 [略]

... やれやれ、疲れた。正直、途中から「これ20年前の話だよな...」という思いにかき乱されて、モチベーションが落ちてしまったんだけど、勉強になりましたです。
 グルーピングの手続きについて、この論文ではまず大きく階層的手法とそうでない手法という観点から分けていて、前者のほうが使い勝手が良い(versatile)だと何箇所かで述べている。クラスタの形状についての仮定がない分柔軟だ、という意味なのだろうと思う。うーん、それって良し悪しだと思うんですけど、どうなんすかね。
 admissibilityについての議論が短いところも、ちょっと意外であった。クラスタリング手法のレビューなんだから、当然Fisherのadmissibility基準で比較するんだろうと思ったのである。著者らが説明していたように、アルゴリズムのadmissibilityだけで形式的に語れる範囲はそう広くない、ってことなんでしょうね。

論文:データ解析(2018-) - 読了:Jain, Murty, Flynn (1999) クラスタリング手法レビュー

Hanley, J.A., Lippman-Hand, A. (1983) If nothing goes wrong, is everything all right? Interpreting zero numerator. Journal of the American Medical Association. 249(13), 1743-1745.
 「いずれ読む」フォルダにPDFが入っていたんだけど、なぜ入れたんだか思い出せない。時々こういうことがあって、結局読まないままになっちゃうんだけど、これはたったの3pなのでざっと目を通した。
 タイトルの通り、割合の推定において分子が0だったらどうするのかという啓蒙論文であった。誌名に見覚えがないな...と思ったのだが、よく見たら現在でいうJAMAである。

 いわく。医師は患者に特定のアウトカムが生じる確率についての推定値を示さなければならない。そのアウトカムが一般的なものでない場合、正確な推定は難しくなる。特に、ある研究においてそのアウトカムの事例が一件もない場合には、一件以上ある場合と比べて質的にも量的にも異なる事態となる...と考えられている。
 本論文では以下を強調したい。(1)分子が0であることは「リスクなし」を意味しない。(2)分子が0だからといってリスクのサイズについて推測できないわけじゃない。(3)分子が非ゼロのとき適用できる推測統計の諸原理は、分子がゼロの時にも適用できる。

 いま、標準的な造影剤があって、深刻なリアクションを示す人が10000人中15人であることがわかっているとする。新しい造影剤を167人に試したところ、深刻なリアクションを示した人はいなかった。リスクはどちらが高いか。
 まずはふたつの造影剤のあいだでリスクが同一(10000人中15人)だとしよう。167人試して出現数0となる確率は、(1-0.0015)^167 = 0.78。ごくありうる事柄である。
 では、新造影剤のリスクが100人中1人だったら? 167人試して出現数0となる確率は、(1-0.01)^167 = 0.19。まだまだ、驚くような値ではない。
 今度は、新造影剤のリスクが25人中1人だったら? 167人試して出現数0となる確率は、(1-0.04)^167 = 0.01。この場合は驚く。真のリスクはもっと低いに違いないという話になる。
 ここで慣例に従い、我々の「驚きたい程度」を5%とし、これを下回ったら「真のリスクはもっと低いに違いない」と考えることにしよう。[←というわけで、途中でやっと気が付いたんだけど、著者らは95%信頼区間についてわかりやすく説明しているわけである]

 さて、観察された割合が$0/n$であるとき、真のリスクは0からある値までの間にある。それと整合するような真のリスクの最大値(=95%信頼区間の上限)をMaximum Riskと呼ぼう。すると
 $(1-Maximum Risk)^n = 0.05$
 $Maximum Risk = 1-0.05^{1/n}$
これは$3/n$で近似できる。なお、ポワソン分布から直接算出しても同じことがいえる。

 [数学の勉強をしてこなかったせいで、こういう話には途中でついていけなくなる。自分向けにゆっくり書くと、
 $0.05^{1/n} = \exp(\log(0.05^{1/n})) = \exp(\log(0.05)/n)$
さて指数関数というものは
 $\exp(x) = \sum_i \frac{x^i}{i!} = 1+x+\frac{x^2}{2}+\frac{x^3}{6}+\cdots$
と展開できる。$x=\log(0.05)/n$と代入して
 $\exp(\log(0.05)/n) = 1 + \frac{\log(0.05)}{n} + \frac{\log(0.05)^2}{2n^2} + \frac{\log(0.05)^3}{6n^3} \cdots$
第二項の$\log(0.05)$はおよそ-3。第三項以降は分母が$2n^2, 6n^3, \ldots$と大きくなっていくので無視できる。というわけで
 $0.05^{1/n} \approx 1 - 3/n$
 $1- 0.05^{1/n} \approx 3/n$
ということであろうかと思う]

 手っ取り早くいうとこういうことだ、観察された割合が$0/n$だったら、真の割合は$3/n$以下である。100人観察して0人だったら、真の割合は3%より以下である。これが"The rule of three"だ、諸君、覚えておきなさい。

 。。。というわけで、なぜ読もうと思ったんだか結局最後まで思い出せなかったんだけど、面白い文章であった。The rule of threeね。こんどセミナーとかやることがあったらネタにしよう。

 途中で余談として出てくるんだけど、かつてラプラスは、過去5000年間にわたって日が昇らなかった日はない、この日数を$n$として、今後日が昇るオッズは$n:1$だ、と述べているそうである。へええ。

論文:データ解析(2018-) - 読了:Hanley & Lippman-Hand (1983) 調べたい事柄の発生率が0で困っているあなたのための「ルール・オブ・スリー」

中田秀基 (2018) 人工知能・機械学習の応用・研究に適したパブリッククラウドサービス. 人工知能, 33(1), 15-22.
タイトルの通りの内容。世の中どんどん変わってキャッチアップできないので、こういう話をまとめて紹介して下さるの、すっごく助かりますです。
音声認識・発話のSaaSとして、AWS, Watson, Google Cloud Platform, MS Azureが挙げられていたが、Azure Custom Speech Serviceというのはアプリケーション固有の言語モデルを与えることができるのだそうだ。へえー。

論文:データ解析(2018-) - 読了:中田(2018) 機械学習のクラウドサービス in 2018

Qi, Z., Wang, H., Li, J., Gao, H. (2018) Impacts of Dirty Data: an Experimental Evaluation. arXiv:1803.06071 [cs.DB]
 たまたま見つけて、面白そうなのでめくってみたもの。著者らの所属はハルビン工科大。
 タイトルの通り、データが汚れていたら機械学習の結果にどう響くか、という実証研究である。面白いことを考えるなあ。実験結果をどこまで一般化していいものかどうかわからないが、テーマとしてはなかなか切実な話である。

 データの品質として次の3つに注目する。

 アルゴリズムとして以下に注目。

 実験です。UCIデータセットから、irisなど13個のデータセットを借りてきて、わざと汚してはアルゴリズムを適用し、汚さない場合と比較した。
 手順の説明ののち、5頁にわたって結果の報告があったが、めんどくさいのでパス。すいません、そこまでの気力はないんです...

 最後に示された、ユーザ向けガイドラインだけメモしておく。
 分類の場合のガイドライン。なお、F指標ってのは(2 x Precision x Recall) / (Precision + Recall)のこと。

  1. まずは所与のデータのエラー率(欠損率、不整合率、矛盾率)を求めましょう。
  2. アルゴリズムを選びましょう。まずprecision, recall, F指標のどれに注目するかを決める。それが70%を超えるアルゴリズムを候補とする。データが小さい場合は、ロジスティック回帰も候補にいれる。
  3. 一番エラー率が高いタイプのエラーに注目し、上で決めた指標について、sensitivityが一番低いアルゴリズムを選ぶ。[sensitivityとは、データが汚れているときにそれに影響される度合いのことで、実験結果が表になっている。実際のところ、どのエラータイプでもどの指標でも、KNNが最低みたいだ]
  4. そのアルゴリズムのkeeping pointを調べ、それに達するまでデータをクリーニングする。[keeping pointとは許容できるエラー率のことで、実験から求めた結果が表になっている。表をみると、注目する指標がF指標の場合、欠損のkeeping pointが一番高いのは決定木だが、不整合のkeeping pontが一番高いのはベイジアン・ネットワーク。ランダム・フォレストはrecallに対する欠損のkeeping pointがやたらに高いが、F指標に対する欠損のkeeping pointはやたらに低い。などなど、いろいろ眺めているといろいろ面白い]

 クラスタリングの場合も上とほぼ同じで、

  1. 所与のデータのエラー率(欠損率、不整合率、矛盾率)を求める。
  2. precision, recall, F指標のどれに注目するかを決め、それが70%を超えるアルゴリズムを候補とする。データが大きい場合はDBSCANも候補にいれる。
  3. 一番エラー率が高いタイプのエラーに注目し、上で決めた指標について、sensitivityが一番低いアルゴリズムを選ぶ。[全体にLVQが低めだ]
  4. そのアルゴリズムのkeeping pointを調べ、それに達するまでデータをクリーニングする。[LVQ, DBSCANが高め。CLARANSは欠損のkeeping pointがめっちゃ低い。へー]

 回帰の場合のガイドライン。なお、NRMSDってのはRMSDを(予測値の最大値)-(最小値)で割ったもの、CVRMSDってのはRMSDを予測値の平均値で割ったもの。

  1. 所与のデータのエラー率(欠損率、不整合率、矛盾率)を求める。
  2. RMSD, NRMSD, CVRMSDのどれに注目するかを決め、それが1(NRMSDなら0.5)%を超えるアルゴリズムを候補とする。
  3. 一番エラー率が高いタイプのエラーに注目し、上で決めた指標について、sensitivityが一番低いアルゴリズムを選ぶ。[全体に最小二乗法が低めだわ...]
  4. そのアルゴリズムのkeeping pointを調べ、それに達するまでデータをクリーニングする。[最小二乗法が高めだ... これが伝統の強さというものか...]

 最後に、今後の課題。
 データ分析・データマイニング関係者のみなさん向けには、

 データ品質・クリーニング関係者のみなさんにはこういいたい。アルゴリズムがエラーに対して持つ頑健性は、エラーのタイプによってもアルゴリズムによっても異なる。だから、いきなりデータを100%クリーンにする必要はないわけだ。むしろ分析課題が決まってから、それにとって適切なkeeping pointまでクリーニングするのがお勧め。

 というわけで、あんまり真面目に読んでないけど、面白い論文であった。
 自分の仕事に引き付けて考えると、私の仕事上のデータ分析はほぼすべてがアドホックなもので、初期の段階では分析課題さえきちんと定義できていないことが多い。分析の進行とともに課題が徐々に明確になり、複数の下位課題にブレイクダウンされていくが、その下位課題によってアルゴリズムも使う変数も変わってくるし、クリーニングが必要な程度も変わってくる。そのせいで、初期段階で時間をかけたクリーニングが結果的には不要であったり、逆にあとでクリーニングをやり直す必要に迫られることも少なくない。どうやって改善すればいいのか、悩ましいところである。
 この論文のアプローチは、「目的とする指標とアルゴリズムによって、クリーニングの目標水準を決める」というもので、分析課題がはっきりしている状況では有用だと思う。でも私の場合、データを一通りクリーニングしないことには身動きが取れないが、しかしクリーニングしてる段階ではまだ先が読めないんですよね...クリーニングをいくつかのフェイズに分け、本格的なクリーニングは個別のモデリングごとにオンデマンドにやるものなのだと覚悟を決める、というのがよいのだろうか。ううむ...

論文:データ解析(2018-) - 読了:Qi, Wang, Li, & Gao (2018) 機械学習におけるデータの「汚れ」の影響

 これは私による私のためのメモでございます。まあ、そう云ってしまえば全部そうなんですけど。

 Rのdplyrパッケージにはscoped variantsと呼ばれる関数(dplyr風にいうと、verb)がある。使い方がいまいちわかっておらず、毎回イライラするので、きちんと調べてみることにした。

 scoped variantsを持つverbは、mutate(), transmute(), summarize(), filter(), group_by(), rename(), select(), arrange()の8つ。
 verbのscoped variantsには all, at, ifの三種類がある。mutate()の場合には、mutate_all(), mutate_at(), mutate_if()があるわけだ。

mutate(), transmute(), summarize()
 以下、変数変換のverbであるmutate()を例に、3つのvariantsについて使い方をメモしておく。transmute(), summarize()も同様である。

_all() variant
 mutate_all(.tbl, .funs, ...)
 オブジェクト.tblに含まれているすべての変数について、.funsに指定した変換を行う。

 さて、このfuns()の書式がちょっと不思議である(そうそう、ここでいつもとまどうのだ)。
 funs(..., .args=list())
 ...に関数を指定する。たとえば変数の平均を求めるとして、...のとこには単にmeanと書いてもいいし、関数名の文字列"mean"と書いてもいいし、mean(., na.rm=T)という風に書いてもいい。しかし無名関数はだめ。function(x) mean(., na.rm=T)というのは許してもらえない。

 funsの中に複数の関数をならべてもよい。ただし、文字列ベクトルは受け取れない。つまり、funs("min", "max")はありだがfuns(c("min", "max"))はなし。後者の場合にはfuns_(c("min", "max"))と書く。
 複数の関数を並べるとなにが起きるかというと...
 data.frame(v1=1:10, v2=1:10) %>% mutate_all(funs(min, max))
 返ってくるデータフレームには、v1, v2, v1_min, v1_max, v2_min, v2_maxが含まれる。接尾辞("min", "max")は勝手につけてくれるけど、明示的にfuns(a =min, b=max)とすれば"a", "b"になる。

_at() variant
 mutate_at(.tbl, .vars, .funs, ...)
 mutate_all()との違いは、.varsに指定した変数だけについて変換が行われるという点である。

 vars()の書式は
 vars(...)
 ...で変数名を指定する。select()と同じように書ける。よってvars(starts_with("hoge"))という風に書ける。

_if() variant
 mutate_if(.tbl, .predicate, .funs, ...)
 mutate_all()との違いは、.predicateに該当する変数だけについて変換が行われるという点である。(そうか... 引数.varsはないのか...)

 試してみたところ、どうやら引数.predicateを無名関数にするのはありらしい。mutate_if(df, function(x) is.numeric(x), as.factor)とか。

 さて、他のverbだとどうなるか。特にわかりにくいのは、引数.funsの意味である。mutate(), transmute(), summarize()は変数変換のverbなので、引数.funsの役割ははっきりしている。他のverbではなにを意味するのか。

filter()
 filter()のscoped variantsは引数.funsを持たない。かわりに引数.vars_predicateを持つ。この引数はall_vars(expr)かany_vars(expr)で生成する。exprは条件式。
 わけわかんなくなって参りましたが、たとえば
 filter_all(df, all_vars(. > 100))
とすると、すべての変数の値が100より大である行が選択される、ということである。

group_by()
 引数.funsを持つ。グループ化変数が変数変換の引数を持つというのも変な感じだが、これはgroup_by()のあとにmutate()するのと同じことになる。たとえば
 group_by_all(df, as.factor)
とすると、すべての変数でのグループ化が行われ、かつそのグループがfactorになる。

select(), rename()
 引数.funsを持つが、これは変数に適用する関数ではなく、変数名に適用する関数となる。ここが直観に反するんだけど、rename()のscoped variantsだけでなく、select()のscoped variantsも引数.funsを持つのである。たとえば
 select_all(df, toupper)
とすると、全ての変数が大文字になるんじゃなくて、全ての変数名が大文字になる。(あれ? じゃあ結局、select_all()とrename_all()って全く同じなの?)

 試してみたところ、無名関数が使える。たとえば
 select_all(df, function(x) paste0("hoge", x))
とすると、全変数の変数名の頭に"hoge"がつく。
 ということはだ。たとえば質問紙調査のデータで、Q1S1, Q2S2, ... に5件法の回答(1,2,3,4,5)がはいっている。これを反転した変数 nImage1, nImage2, ...と、T2Bのとき1になるダミー変数bImage1, bImage2, ..をつくりたい。この場合、
 df %>%
  mutate_at(vars(starts_with("Q1S"), funs(n = 6 - ., b = if_else(. %in% 1:2, 1, 0))) %>%
  rename_at(vars(starts_with("Q1S"), sub("Q1S(.+)_(.+)", "\\2Image_\\1", .))
とすればいいわけ? それとも
  rename_at(vars(starts_with("Q1S"), function(x) sub("Q1S(.+)_(.+)", "\\2Image_\\1", x))
と無名関数にするわけ? 試してみたところ、後者のみOKであった。

arrange()
 引数.funsを持つ。たとえば
 arrange_all(df, desc)
とすると、左の変数から順に降順になる。

雑記:データ解析 - 覚え書き: Rのdplyrパッケージのhogehoge_{all, if, at} 動詞の使い方

Richardson, D.B., Kinlaw, A.C., Keil, A.P., Naimi, A.I., Kaufman, J.S., Cole, S.R. (2018) Inverse Probability Weights for the Analysis of Polytomous Outcomes. American Journal of Epidemiology.

 多項回帰モデルは共変量が多いと推定が大変なので、ウェイティングでなんとかしましょう、という話。
 いずれ読もうと思い、「いずれ読む」という名のフォルダ(実質的には墓場)にPDFを放り込みかけたんだけど、行間の広いdraftなのに10pしかないことに気づき、そんならさっさと目を通しておこう、と仕事を中断して目を通した。要するに現実逃避である。

 アウトカムを$D$, そのカテゴリを$0, \ldots, G$とする。曝露変数を$E$とし、話を簡単にするために二値とする。共変量を$\mathbf{Z}=\{Z_1, \ldots, Z_k\}$とする。
 研究者は、$Z$を調整しつつ、曝露ありとなしの間で$D$の分布を比べたい。多項回帰なら、参照レベルを$D=0$として
$\displaystyle \log\frac{P(D=g|E=e, \mathbf{Z}=z)}{P(D=0|E=e, \mathbf{Z}=z)} = \alpha_g + \beta_g E + \mathbf{Z} \gamma_g$
とするところだ。するってえと、$\alpha$を$G$個、$\beta$を$G$個、長さ$k$の$\gamma$を$G$個、推定する羽目になるわね。そんなん推定してられっか、と思いません? そんなあなたのための新手法!名付けて「inverse-probability-of-exposure weighted多項ロジスティック回帰モデル」です!
 [正直言って、名前をみただけでなにすんのか見当が付いちゃうわけで、すごい出落ち感がある...まあ最後まで読みますけど...]

 ここからは実例でご覧いただきましょう。
 子宮内膜がん患者の横断調査(n=288)。知りたいのはがんの組織学的下位タイプと年齢の関連性。説明変数は年齢のみ、二値とする(64歳まで, 65歳以上)。アウトカムは3タイプある。共変量は喫煙有無など3つ、すべて二値。
 まず共変量で年齢の確率を予測するモデルを推定します。次に、この確率の逆数(に、年齢の周辺確率を掛けた値)を重みにし、年齢だけを説明変数にして、重みつきの多項回帰を推定します。
 シミュレーションしてみると、ふつうの多項モデルが収束しない時にもこの方法なら収束する。パラメータ数が多すぎる多項モデルをどうにかする方法としては、ほかに縮小推定なんかもあるけれど、多くの場合こっちのほうが簡単なのではないでしょうか。
 なお、この手法は要するに曝露の傾向スコアと関連しておりましてですね、本命の多項モデルが一致性を持つためには、その前の曝露予測モデルの指定が正しいことが必要であります。
 云々。

 。。。という、ごく短い報告であった。ま、気分転換っていうことで、ひとつ。
 詳細は付録をみよとのことだが、入手できない。なんだかなあ。

論文:データ解析(2018-) - 読了:Richarson, et al. (2018) 共変量のある多項回帰モデルを、共変量で予測した曝露確率の逆数でウェイティングして共変量なしで済ませる

 対応のある二値データについて考える。たとえば、n人の人に製品PとQをみせ、それぞれについて買ってみたいかどうかを判断してもらった、というような場合ですね。両方の製品について「はい」と答えた人数を$a$, 製品Qについてのみ「はい」と答えた人数を$b$, 製品Pについてのみ「はい」と答えた人数を$c$, 両方の製品について「いいえ」と答えた人数を$d$とする(当然ながら$a+b+c+d=n$)。
 製品PとQのあいだで「はい」率は異なるか。これを調べる古典的な方法が、ご存じMcNemar検定である。母集団における「はい」率が製品P,Q間で同一であるという帰無仮説の下で、$(b-c)^2/(b+c)$が自由度1のカイ二乗分布に従うことを利用する。
 分野の違いというのは面白いもので、市場調査の参考書ではこれを「対応のあるZ検定」と呼び、検定統計量$(b-c)/\sqrt{b+c}$を$N(0,1)$と比べよ、と書いてあることが多い。まあどっちでも同じことである。

 ところがこのMcNemar検定、自分で手計算してみるとわかるけど、ちょっと不思議な面がある。
 たとえば、サンプルサイズは$n=100$、うちQについてのみ「はい」と答えた人数が$b=15$, Pについてのみ「はい」と答えた人数が$c=10$だったとしよう。このときの検定統計量の値は$(b-c)^2/(b+c)=25/25=1$。
 今度は、サンプルサイズは$n=100000$, うち$b=15, c=10$だったとしよう。このときも$(b-c)^2/(b+c)=25/25=1$。
 つまり、サンプルサイズがどんなに変わろうと、2x2クロス表の対角セル$a, d$がどんなに変わろうと、非対角セル$b$と$c$さえ同じなら、検定統計量の値は同じなのだ。
 でも、そもそもここで比較したいのは、Pに「はい」と答えた人の割合$(a+c)/n$と、Qに「はい」と答えた人の割合$(a+b)/n$でしょう? サンプルサイズ$n$が大きいとき、これらの割合の信頼性は高くなるはずでしょう?

 ...という素朴な疑問は、私自身がかつて感じたことでもあり、前職の社内研修のときに訊ねられたことでもある。このたび仕事のなかで、ふとこの話を思い出す機会があり、試しにAgrestiの分厚い教科書をめくってみたら、その話はAgresti & Min (2003)に書いたから読んでね、とのことであった。

Agresti, A., & Min, Y. (2003) Effects and non-effects of paired identical observations in comparing proportions with binary matched-pairs data. Statistics in Medicine, 22. 65-75.

 というわけで、読んでみました。難しくてよくわかんない箇所が多かったんだけど...

 あるペア(たとえばあるヒト)$i=1, \ldots, n$が提供する$t=1,2$個目の値を$y_{it}$とする。値$1$を成功、$0$を失敗と呼ぶ。
 このデータを表す形式を2種類考えよう。

 subject-specific形式のデータに当てはめられる標準的モデルは、
 $link[P(y_{i1})=1] = \alpha_i$
 $link[P(y_{i2})=1] = \alpha_i+\beta_c$
であろう。これをconditional modelという(効果$\beta_c$の定義がペアに条件づけられているという意味)。linkというのはたとえばlogit。モデル・フィッティングにおいては、$(y_{i1}, y_{i2})$が$\{\alpha_i\}$の下で独立だと仮定するのが普通である。

 いっぽう、ランダムに取り出したあるヒトの最初の観察を$y_1$, ランダムに取り出した別のヒトの二番目の観察を$y_2$として
 $link[P(y_1)=1] = \alpha$
 $link[P(y_2)=1] = \alpha+\beta_m$
とモデル化するという手もある。これをmarginal modelという。ペアの母集団を通じた平均効果のモデル化である。population-averaged形式の2x2クロス表の周辺分布のモデル化だといえる。

 ここでちょっとこの論文のメモから離れて、Agrestiの"An Introduction to Categorical Data Analysis"からメモしておくと;
 一般化線形モデルのパラメータ$\beta$について$H_0: \beta = 0$を検定する方法には3つある。

 論文に戻って。。。
 $\pi_1=P(y_1=1)$と$\pi_2=P(y_2=1)$が等しいという帰無仮説の検定について考えよう。この帰無仮説は$\beta_m=0$と等価であり、$\pi_{12}=\pi_{21}$と等価である。
 marginal modelに基づいて考えると、$\beta_m=0$という条件の下で最大化された尤度と、条件なしで最大化された尤度の比は、$b,c$だけに依存している。尤度比統計量(=-2*対数尤度比)は$2b \log(2b/(b+c)) + 2c \log(2c/(b+c))$となる。
 $b+c$が小さい場合、$c$を試行数$b+c$, パラメータ$1/2$の二項分布と比べる正確検定が可能である。大標本の場合、これは正規近似できて、
 $\displaystyle z=\frac{c-(1/2)(b+c)}{\sqrt{(b+c)(1/2)(1/2)}} = \frac{c-b}{\sqrt{b+c}}$
これを二乗したのがMcNemar統計量である。これは周辺分布が等しいという検定のためのスコア統計量になっている。
 ワルド統計量はどうか。ワルド統計量は、推定された効果をSEで割ったものである。だからリンク関数の定義によって変わってくる。でも後述するように、$a,c$への依存性は無視できる程度である。
 このように、この検定において情報を提供するのは$b$と$c$だけである。

 ここからは、割合の差、オッズ比、相対リスクの推定、という3つの観点から考える。

 まず割合の差の推定について。ここ、短いんだけど途中からいきなりわかりにくくなるので、思い余って全訳した。

 identity linkの場合、subject-specificな効果とpopulation-averagedな効果は同一である。たとえばconditional modelでは、すべての$i$について$\beta_c = P(Y_{i2}=1)-P(Y_{i1}=1)$である。これを母集団内のペアを通じて平均すると、marginal modelにおける周辺確率の差$\beta_m=\pi_2-\pi_1$と等しくなる。
 別の見方として次の結果に着目することもできる。indentity linkのconditional modelを$2 \times 2 \times n$表に適用したとき、層別変数が説明変数と周辺独立であるときにcollapsibilityが成立する。このとき、これら2変数の標本クロス表は全セルが1となる。[←先生すいません、この段落、なにいってんだかさっぱりわかんないです]
 標本における割合の差は
 $p_2 = p_1 = p_{21} - p_{12} = (c-b)/n$
となる。これはindentity linkのmarginal modelにおける$\beta_m$の最尤推定値である。conditional modelの観点からいえば、これはヒト別の表のMantel-Haenszel型加重平均からも得られる。またこれは、Chen(1996)が示しているように、$\{\alpha_i\}$がなんらかのparametric familyからランダムに抽出されたものだと仮定した、ランダム効果版conditional modelでの最尤推定値になっている。
 周辺表の多項抽出において
 $Var(p_2-p_1) = [(\pi_{12}+\pi{21}) - (\pi_{21}-\pi_{12})^2]/n$
その標本推定値は
 $\hat{Var}(p_2-p_1) = \frac{(b+c)-(c-b)^2/n}{n^2}$
周辺等質性($\pi_1=\pi_2$)という帰無仮説の下で、分散は$(\pi_{12}+\pi_{21})/n$に、その推定値は$(b+c)/n$に、それぞれ帰着する。ここからMcNemar検定が導かれる。$n=a+b+c+d$であるから、帰無仮説が真でない時、効果のサイズの推定値とその標準誤差は、$(a+d)$が増大するにつれて減少する。
 ここでも、また以下のモデルでも、$(a+d)$の効果について考える際には、$b$と$c$を固定して$(a+d)$だけを大きくすると考える。$n$自体も大きくなるので、帰無仮説のもとでの漸近性において$a$, $d$が果たす役割をチェックすることができる。
 $(p_2-p_1)$の分散の推定値は、$n^{-2}$のオーダーで、$(b+c)/n^2$すなわち帰無仮説の下での分散の推定値になる。同様に、$n$が大きければ、$(p_2-p_1)$とその標準誤差の推定値との比は近似的に$(c-b)/\sqrt{b+c}$となる(これはMcNemar統計量の標準正規形式である)。このように、とりわけ$n$が大きいとき、割合の差のワルド検定における$(a+d)$の寄与はごくわずかである。

 では、リンクlogitである場合(パラメータは対数オッズ比になる)は...[パス]
 リンクがlogである場合(対数相対リスク)は...[パス]

 ベイズ推論の観点から見るとどうか。
 marginal modelの場合。$\{\pi_{jk}\}$の事前分布としてパラメータ$\{\mu_{jk}\}$のディリクレ分布を与えるとしよう。$\pi_2 > \pi_1$の事後確率は、$\pi_{21} / (\pi_{12}+\pi_{21}) > 1/2$の事後確率に等しい。左辺の事後分布はパラメータ$(c+\mu_{21}, b+\mu_{12})$のベータ分布になる。だから$a, d$は効かない。ところが、conditional modelになると話が変わってきて...[めんどくさいので省略]

 云々、云々。。。というわけで、$a, d$は効果のサイズやそのSEに全然効かないわけじゃないけど、たいして効かない、ということなのだそうである。すいませんAgresti先生、途中で心が折れました。

論文:データ解析(2018-) - 読了:Agresti & Min (2003) 対応のある二条件間で割合を比べるとき、「両方1」「両方0」の事例を無視してよいのか

2018年5月 3日 (木)

Bookcover 新装版 隠し剣孤影抄 (文春文庫) [a]
藤沢 周平 / 文藝春秋 / 2004-06-01
Bookcover 新装版 隠し剣秋風抄 (文春文庫) [a]
藤沢 周平 / 文藝春秋 / 2004-06-01
藤沢周平に手を出したら人生も終わりだと思ってたんだけど(すいません。そのくらい面白そうだってことです)、ついつい読んでしまった。たまたま映画「必死剣鳥刺し」(2010)を観て、トヨエツ演じるオッサンと池脇千鶴演じる姪がついついヤッちゃうところがどうにも納得いかなくて、いったい原作はどうなっているのかいなと読んでみた次第。原作だと... これはまあしょうがないかな、人間だもの、という感じでした。

Bookcover マルタの鷹〔改訳決定版〕 (ハヤカワ・ミステリ文庫) [a]
ダシール ハメット / 早川書房 / 2012-09-07

Bookcover 沿岸―頼むから静かに死んでくれ (コレクション現代フランス語圏演劇) [a]
ワジディ ムアワッド / れんが書房新社 / 2010-06
ケベックの劇作家Wajdi Mouawadの初期作品。日本では2010年に静岡でMouawad演出による来日公演があり(「頼むから静かに死んでくれ」)、今年になって世田谷で上村聡史演出で上演された(「岸 リトラル」)。

Bookcover プレイバック [a]
レイモンド チャンドラー / 早川書房 / 2016-12-08
Bookcover 高い窓 [a]
レイモンド チャンドラー / 早川書房 / 2014-12-05
村上春樹によるチャンドラー7作品の翻訳、先日完結した模様。いずれもすでに読んでいるんだけど(たしかすべて清水訳で)、村上訳という魅力にはあらがえず、ちびちび読み直している次第。
 「高い窓」はとにかく細部が面白く、特に終盤で娘を実家に送り届けるくだりが素晴らしい、という記憶があった。読み直してみたところ、その印象は変わらないんだけど、二人が車で実家に向かい、途中で湖畔を通りかかる、という情景描写があったような気がしてたんですよね。実際には、あっというまに実家についちゃいました。いつのまにか記憶のなかででっち上げていたようだ。
 えーっと、残っているのは「大いなる眠り」と「湖中の女」か。いつ読めるのかわからないけれど。

フィクション - 読了:「高い窓」「プレイバック」「マルタの鷹」「隠し剣孤影抄」「隠し剣秋風抄」「沿岸 頼むから静かに死んでくれ」

Bookcover JKハルは異世界で娼婦になった [a]
平鳥コウ / 早川書房 / 2017-12-06
 ライト・ノベルは全く読まないんだけど、評判を聞いて好奇心に駆られ、年末年始の移動中に読んだ。投稿サイト連載の書籍化、早川書房としては初のケースだそうだ。
 娼館で生きる娘の一人称小説。彼女はもとは現代の女子高生で、車にはねられて異世界に転生したのだという設定だが、そこはどうでもよくて、これは一種の青春小説というか成長小説なんでしょうね。正直なところ特に感想がないが、それはこの小説が扱っている問題といまの私の関心がずれているからであって、面白い小説ではあった。

Bookcover ラブラバ〔新訳版〕 (ハヤカワ・ポケット・ミステリ) [a]
エルモア レナード / 早川書房 / 2017-12-06
 このたびポケミスで出た新訳で再読。レナード一流の小粋な犯罪小説である。えーと、83年の作品だそうだ。
 調べてみたところ、あれだけたくさん翻訳されていたレナードの小説は、この本を除き現在新刊書としてはいずれも入手できない模様。なんと残念な...と思う反面、なんだかレナードらしいなという気もする。
 追記:上記は確か年明けにメモしたんだけど、その後なんと村上春樹の訳で短編集が出た。本屋さんに訊いてみたところ、売れている由であった。これを機に他の本も復刊されるといいんですけど。

Bookcover コリーニ事件 (創元推理文庫) [a]
フェルディナント・フォン・シーラッハ / 東京創元社 / 2017-12-11

Bookcover シェイクスピア全集29 アテネのタイモン (ちくま文庫) [a]
シェイクスピア / 筑摩書房 / 2017-10-06

Bookcover シンパサイザー [a]
ヴィエト タン ウェン,Viet Thanh Nguyen / 早川書房 / 2017-08-24

フィクション - 読了:「JKハルは異世界で娼婦になった」「ラブラバ」「コリーニ事件」「アテネのタイモン」「シンパサイザー」

Bookcover 応援される会社 熱いファンがつく仕組みづくり (光文社新書) [a]
新井範子,山川悟 / 光文社 / 2018-01-17

マーケティング - 読了:「応援される会社」

Bookcover アレント入門 (ちくま新書1229) [a]
中山 元 / 筑摩書房 / 2017-01-05

哲学・思想(2011-) - 読了:「アレント入門」

Bookcover 戦前日本のポピュリズム - 日米戦争への道 (中公新書) [a]
筒井 清忠 / 中央公論新社 / 2018-01-19

Bookcover 高度成長―シリーズ日本近現代史〈8〉 (岩波新書) [a]
武田 晴人 / 岩波書店 / 2008-04-22

Bookcover 戦争調査会 幻の政府文書を読み解く (講談社現代新書) [a]
井上 寿一 / 講談社 / 2017-11-15

日本近現代史 - 読了:「高度成長」「戦争調査会」「戦前日本のポピュリズム」

Bookcover 私はガス室の「特殊任務」をしていた (河出文庫) [a]
シュロモ ヴェネツィア / 河出書房新社 / 2018-04-06

ノンフィクション(2018-) - 読了:「私はガス室の特殊任務をしていた」

Bookcover 睡眠負債 『ちょっと寝不足』が命を縮める (朝日新書) [a]
NHKスペシャル取材班 / 朝日新聞出版 / 2018-04-13
せっかくNHKの頭のいい人たちがお金をかけて取材して、朝日新聞出版から大部数の本を出すんだから、もうすこしきちんと校閲すればいいのに... 「母数」とか「サンプル数」とか「母集団の偏り」とかさ... まあ細かいことですけど...

Bookcover 内村鑑三 悲しみの使徒 (岩波新書) [a]
若松 英輔 / 岩波書店 / 2018-01-20

Bookcover クリムト 官能の世界へ (角川新書) [a]
平松洋 / KADOKAWA / 2018-01-10

Bookcover 公文書問題 日本の「闇」の核心 (集英社新書) [a]
瀬畑 源 / 集英社 / 2018-02-16

Bookcover 官僚たちのアベノミクス――異形の経済政策はいかに作られたか (岩波新書) [a]
軽部 謙介 / 岩波書店 / 2018-02-21

Bookcover 矢内原忠雄――戦争と知識人の使命 (岩波新書) [a]
赤江 達也 / 岩波書店 / 2017-06-21

ノンフィクション(2018-) - 読了:「睡眠負債」「内村鑑三」「クリムト」「公文書問題」「官僚たちのアベノミクス」「矢内原忠雄」

Bookcover 千の顔をもつ英雄〔新訳版〕上 (ハヤカワ・ノンフィクション文庫) [a]
ジョーゼフ・キャンベル / 早川書房 / 2015-12-18

Bookcover ブッダたちの仏教 (ちくま新書) [a]
並川 孝儀 / 筑摩書房 / 2017-12-06

Bookcover 労基署は見ている。 (日経プレミアシリーズ) [a]
原 論 / 日本経済新聞出版社 / 2017-03-09

Bookcover 兼好法師 - 徒然草に記されなかった真実 (中公新書) [a]
小川 剛生 / 中央公論新社 / 2017-11-18

Bookcover 北朝鮮核危機! 全内幕 (朝日新書) [a]
牧野 愛博 / 朝日新聞出版 / 2018-02-13

Bookcover 和本のすすめ――江戸を読み解くために (岩波新書) [a]
中野 三敏 / 岩波書店 / 2011-10-21

ノンフィクション(2018-) - 読了:「和本のすすめ」「北朝鮮核危機全内幕」「兼好法師」「労基署はみている」「ブッダたちの仏教」「千の顔を持つ英雄」

Bookcover 釜ケ崎と福音――神は貧しく小さくされた者と共に (岩波現代文庫) [a]
本田 哲郎 / 岩波書店 / 2015-02-18
著者はカトリックの神父さんで、釜ヶ崎の支援活動で知られる人。感想は省略するけど、とても興味深く、印象に残る本であった。

Bookcover 剣と清貧のヨーロッパ - 中世の騎士修道会と托鉢修道会 (中公新書) [a]
佐藤 彰一 / 中央公論新社 / 2017-12-20

Bookcover ポピュリズムとは何か - 民主主義の敵か、改革の希望か (中公新書) [a]
水島 治郎 / 中央公論新社 / 2016-12-19

Bookcover 大量生産品のデザイン論 経済と文化を分けない思考 (PHP新書) [a]
佐藤 卓 / PHP研究所 / 2017-12-17

Bookcover 朱子学と陽明学 (ちくま学芸文庫) [a]
小島 毅 / 筑摩書房 / 2013-09-10

ノンフィクション(2018-) - 読了:「剣と清貧のヨーロッパ」「ポピュリズムとはなにか」「大量生産品のデザイン論」「釜ヶ崎と福音」「朱子学と陽明学」

三浦しをん「広辞苑をつくるひと」岩波書店.
 今年発売された広辞苑第七版、ついつい買ってしまった。実のところ、普段はPCにいれた第五版のEPWING版を使っており(とても便利です)、紙を使う用事はあまりないんだけど、ちょっとしたノスタルジーに駆られた次第である。仕方がないので、机の脇に置いてたまに意味もなく眺めている。
 これは予約特典でついてきた文庫判の小冊子。辞書編集部を舞台にした人気作「舟を編む」の作家による、広辞苑の製作に関わった裏方さんたちの取材記である。
 買ってきた広辞苑の紙の函、場所を取るから潰して捨てるつもりだったのに、これを読んだせいでなんだか捨てにくくなってしまい、ちょっと困っている。

Bookcover 図説 和菓子の歴史 (ちくま学芸文庫) [a]
青木 直己 / 筑摩書房 / 2017-08-07

Bookcover 奇想の図譜 (ちくま学芸文庫) [a]
辻 惟雄 / 筑摩書房 / 2005-04-01

Bookcover ルポ タックスヘイブン 秘密文書が暴く、税逃れのリアル (朝日新書) [a]
朝日新聞ICIJ取材班 / 朝日新聞出版 / 2018-04-13
内容もさることながら、朝日新聞・共同通信・NHKが国際報道プロジェクトに関わるまでのそれぞれの経緯が面白かった。やはり個人的つながりが重要であった模様。
 朝日新聞ではもともと調査報道NPOとの連携を模索していた由。その背後には、調査報道を生き延びさせていくために、いつか特別報道チームを朝日新聞からスピンアウトさせ、NPOとして独立させなければならないような事態が来るのではないか...という問題意識があったのだそうである。へえええ。

ノンフィクション(2018-) - 読了:「広辞苑をつくるひと」「図説 和菓子の歴史」「奇想の図譜」「ルポ タックスヘイブン」

Bookcover ε-δ論法とその形成 [a]
中根 美知代 / 共立出版 / 2010-07-23
  たまにふらふらと、絶対に理解できそうにない本を買ってしまうことがある。これもその一つで、微分・積分の歴史についての本。そもそも「ε-δ論法」というのがなんなのかもよくわかっていないのに、なぜこんな本を買おうと思ったのか、自分でもさっぱりわからない。家に帰って袋から出して、ちょっとあっけにとられた。
 というわけで、もちろん9割9分9厘までちんぷんかんぷんなのだが、せっかくなので最後のページまでめくった。いくつか印象に残ったところをメモしておく。

データ解析 - 読了:「ε-δ論法とその形成」

Bookcover 喜劇としての国際ビジネス: 私が出会った〈一流〉という名の怪優たち [a]
ダニエル・レヴィン / 創元社 / 2017-12-20
帯にD.カーネマンの推薦文が載ってたのでついうっかり買っちゃった本。法律家の体験談というか自慢話というか、まあそういう読み物であった。カーネマン先生ってば、もう。

Bookcover 辺境の怪書、歴史の驚書、ハードボイルド読書合戦 [a]
高野 秀行,清水 克行 / 集英社インターナショナル / 2018-04-05

Bookcover 闘いを記憶する百姓たち: 江戸時代の裁判学習帳 (歴史文化ライブラリー) [a]
八鍬 友広 / 吉川弘文館 / 2017-09-15

Bookcover 習近平帝国の暗号 2035 [a]
中澤 克二 / 日本経済新聞出版社 / 2018-03-09

Bookcover 都市の舞台俳優たち:アーバニズムの下位文化理論の検証に向かって (リベラ・シリーズ11) [a]
田村公人 / ハーベスト社 / 2015-05-28
小劇場演劇に関わる青年たちのエスノグラフィー。都市社会学の研究としての価値はよくわからないんだけど(修辞表現ではなく、本当に判断がつかない)、ノンフィクションとしてとても面白い本であった。

ノンフィクション(2018-) - 読了:「喜劇としての国際ビジネス」「都市の舞台俳優たち」「習近平帝国の暗号2035」「闘いを記憶する百姓たち」「辺境の怪書、歴史の驚書、ハードボイルド読書合戦」

Bookcover 鎌倉幕府の滅亡 (歴史文化ライブラリー) [a]
細川 重男 / 吉川弘文館 / 2011-02-01

Bookcover 例外時代 [a]
マルク・レヴィンソン / みすず書房 / 2017-11-16

Bookcover 和本への招待 日本人と書物の歴史 (角川選書) [a]
橋口 侯之介 / KADOKAWA/角川学芸出版 / 2011-06-25

Bookcover 誰がアパレルを殺すのか [a]
杉原 淳一,染原 睦美 / 日経BP社 / 2017-05-25
帯にいわゆる"有名ブロガー"の推薦文が載っていて、レジに持っていきかけて「え、俺こういう人が推薦するような本買うの...まじ...」としばし躊躇ったのだが、面白い内容の本であった。

Bookcover 暴君誕生――私たちの民主主義が壊れるまでに起こったことのすべて [a]
マット・タイービ / ダイヤモンド社 / 2017-12-21
毒舌コラムニストの時評集。これは別に読まなくてもよかったな...

ノンフィクション(2018-) - 読了:「鎌倉幕府の滅亡」「例外時代」「和本への招待」「誰がアパレルを殺すのか」「暴君誕生」

Bookcover ポピュリズムとは何か [a]
ヤン=ヴェルナー・ミュラー / 岩波書店 / 2017-04-19

Bookcover AI vs. 教科書が読めない子どもたち [a]
新井 紀子 / 東洋経済新報社 / 2018-02-02

Bookcover 炎と怒り――トランプ政権の内幕 [a]
マイケル ウォルフ,Michael Wolff / 早川書房 / 2018-02-23

Bookcover 軌道 福知山線脱線事故 JR西日本を変えた闘い [a]
松本 創 / 東洋経済新報社 / 2018-04-06

Bookcover エスカレーション――北朝鮮VS.安保理 四半世紀の攻防 [a]
藤田 直央 / 岩波書店 / 2017-12-21

ノンフィクション(2018-) - 読了:「炎と怒り」「軌道」「エスカレーション」「AI vs 教科書が読めない子どもたち」「ポピュリズムとは何か」

Bookcover 3億人の中国農民工 食いつめものブルース [a]
山田 泰司 / 日経BP社 / 2017-11-09

Bookcover はじめての世界音楽―諸民族の伝統音楽からポップスまで [a]
柘植 元一,塚田 健一 / 音楽之友社 / 1999-05-01
しばらく前に浜松の楽器博物館に行って、世界の楽器のあまりの多様さに腰を抜かし、帰ってきてから思わず買っちゃった本。大学の教科書みたいだ。

Bookcover 仮想通貨 [a]
岡田 仁志,高橋 郁夫,山﨑 重一郎 / 東洋経済新報社 / 2015-05-29
2015年刊。世の中の動きについていけないという不安に駆られ、ついついこんな本を買っちゃったりして。

Bookcover テヘランからきた男 西田厚聰と東芝壊滅 [a]
児玉 博 / 小学館 / 2017-11-15

Bookcover 東芝の悲劇 [a]
大鹿 靖明 / 幻冬舎 / 2017-09-21

ノンフィクション(2018-) - 読了:「食いつめものブルース」「はじめての世界音楽」「仮想通貨 技術・法律・制度」「テヘランから来た男」

Bookcover 実録 泣くまでボコられてはじめて恋に落ちました。(1) (BUNCH COMICS) [a]
ペス山ポピー / 新潮社 / 2018-04-09
こ、これは... 読んだ人は誰もがなにかしら言いたくなるであろう、強烈なエッセイマンガである。うかつに語れない内容なので、感想を書くのは控えておくけれど、この企画を通した新潮社はたいしたものだ...

Bookcover 恋は雨上がりのように 10 (ビッグコミックス) [a]
眉月 じゅん / 小学館 / 2018-04-27
美しい女子高生がなぜか冴えない中年男に恋をして、という人気作の最終巻。とてもさわやかな結末であった。

Bookcover 北北西に曇と往け 1巻 (ハルタコミックス) [a]
入江 亜季 / KADOKAWA / 2017-10-13

Bookcover 銃座のウルナ 5 (ビームコミックス) [a]
伊図透 / KADOKAWA / 2018-03-28
いまもっとも続きが気になる連載作。なぜこれがもっと評判にならないのか、打ち切りになっちゃったらどうすんだ、と心配していたのだが、文化庁芸術祭の優秀賞となったそうだ。よかったよかった。

コミックス(2015-) - 読了:「実録 泣くまでボコられてはじめて恋に落ちました」「恋は雨上がりのように」「銃座のウルナ」「北北西に曇と往け」

Bookcover 鉄工所にも花が咲く [a]
野村 宗弘 / 太田出版 / 2017-12-07
私はこの作家さんの作品をデビュー作「とろける鉄工所」以来ずっと追いかけていて、どんどん腕を上げていく様を目の当たりにしているんだけど、この連作の第五話「元工場長(70)」が現時点での到達点、最高傑作だと思う。隅から隅までもう完璧で、溜息が出る。

Bookcover ボンクラボンボンハウス 2 (フィールコミックス) [a]
ねむようこ / 祥伝社 / 2018-03-08

Bookcover インド夫婦茶碗 (24) (ぶんか社コミックス) [a]
流水りんこ / ぶんか社 / 2018-03-10
もはや面白いかどうかではなく、季節のお便りのような気分で読み続けていたんだけど、ここでいったん休載とのこと。お疲れ様でした。

Bookcover よつばと!(14) (電撃コミックス) [a]
あずま きよひこ / KADOKAWA / 2018-04-28
現代日本マンガを代表するお化け作品の最新刊。帯によれば、翻訳版が13言語・300万部を超えるそうだ(目を疑うが、翻訳版だけで300万部である。国内版は累計1370万部)。

コミックス(2015-) - 読了:「よつばと!」「鉄工所にも花が咲く」「ボンクラボンボンハウス」「インド夫婦茶碗」

Bookcover テンジュの国(1) (KCデラックス 週刊少年マガジン) [a]
泉 一聞 / 講談社 / 2018-03-09
講談社ってば、よくもまあ、ここまで「乙嫁語り」に似ているコンセプトの連載を... と呆れたけれど、でもまあ別のマンガだし、面白ければそれでいいですよね。

Bookcover 蟇の血 (ビームコミックス) [a]
近藤 ようこ / KADOKAWA / 2018-02-10
原作は田中貢太郎という大正・昭和期の作家だそうだ。へえー。

Bookcover ランド(6) (モーニング KC) [a]
山下 和美 / 講談社 / 2018-03-23

Bookcover ど根性ガエルの娘 4 (ヤングアニマルコミックス) [a]
大月悠祐子 / 白泉社 / 2018-03-29

Bookcover カフェでカフィを 2 (集英社クリエイティブコミックス) [a]
ヨコイ エミ / 集英社クリエイティブ / 2018-04-25

山本直樹「分校の人たち 3」, 太田出版, 2018.
なんと、amazonでは取り扱っておらず、APIで書影を得ることもできない。おそらく性描写のせいだとおもうんだけど... うーん、これって良いことなのだろうか...

コミックス(2015-) - 読了:「テンジュの国」「蟇の血」「ランド」「ど根性ガエルの娘」「カフェでカフィを」「そしてボクは外道マンになる」「分校の人たち」

Bookcover ものするひと 1 (ビームコミックス) [a]
オカヤ イヅミ / KADOKAWA / 2018-03-12
純文学作家の青年の静かな生活を描いた高野文子風のマンガ。正直いって好きなタイプのマンガではないんだけど(勝手にやってろという気がする)、好きな人の気持ちもわからないでもない、という感じ。

Bookcover ゴールデンカムイ 1 (ヤングジャンプコミックス) [a]
野田 サトル / 集英社 / 2015-01-19
Bookcover ゴールデンカムイ 2 (ヤングジャンプコミックス) [a]
野田 サトル / 集英社 / 2015-02-19
Bookcover ゴールデンカムイ 3 (ヤングジャンプコミックス) [a]
野田 サトル / 集英社 / 2015-05-19
2015年から連載している、明治期の北海道を舞台にした大冒険活劇。主人公の兵士とアイヌの少女が、互いを尊重する対等な存在として共闘する、というところがとても現代的だと思う。なるほど、こりゃあ評判になるわけだ。

Bookcover ダンジョン飯 6巻 (ハルタコミックス) [a]
九井 諒子 / KADOKAWA / 2018-04-13

Bookcover タイムスリップオタガール(3) (ポラリスCOMICS) [a]
佐々木 陽子 / フレックスコミックス / 2018-04-12

Bookcover 7人のシェイクスピア NON SANZ DROICT(4) (ヤンマガKCスペシャル) [a]
ハロルド 作石 / 講談社 / 2018-03-06
おかげさまで、「ヘンリー六世」の上演は成功を収めました。次巻はいよいよ「リチャード三世」であろうか。

コミックス(2015-) - 読了:「ものするひと」「ゴールデンカムイ」「ダンジョン飯」「タイムスリップオタガール」「7人のシェイクスピア」

Bookcover 乙嫁語り 10巻 (ハルタコミックス) [a]
森 薫 / KADOKAWA / 2018-02-15

Bookcover ワカコ酒 10 (ゼノンコミックス) [a]
新久千映 / 徳間書店 / 2018-02-20

Bookcover あんたさぁ、 (ビッグコミックススペシャル) [a]
ルネッサンス吉田 / 小学館 / 2018-02-09

Bookcover 健康で文化的な最低限度の生活 6 (ビッグコミックス) [a]
柏木 ハルコ / 小学館 / 2018-01-30

Bookcover あさひなぐ 25 (ビッグコミックス) [a]
こざき 亜衣 / 小学館 / 2018-01-30

Bookcover ディザインズ(3) (アフタヌーンKC) [a]
五十嵐 大介 / 講談社 / 2018-01-23

Bookcover 二月の勝者 ー絶対合格の教室ー 1 (ビッグコミックス) [a]
高瀬 志帆 / 小学館 / 2018-02-09
中学進学塾を題材にした、いかにも小学館的な情報マンガなんだけど、それはそれで興味深い。

Bookcover あれよ星屑 7 (ビームコミックス) [a]
山田 参助 / KADOKAWA / 2018-02-10
最終巻。ここ数年で読んだ漫画のなかで、「たそがれたかこ」「銃座のウルナ」に並ぶ文句なしの大傑作であった。

コミックス(2015-) - 読了:「乙嫁語り」「ワカコ酒」「あんたさぁ、」「健康で文化的な最低限度の生活」「あさひなぐ」「ディザインズ」「あれよ星屑」「二月の勝者」

Bookcover 大砲とスタンプ(7) (モーニング KC) [a]
速水 螺旋人 / 講談社 / 2017-12-21

Bookcover 甘々と稲妻(10) (アフタヌーンKC) [a]
雨隠 ギド / 講談社 / 2018-01-05

Bookcover 大人スキップ 2 (ビームコミックス) [a]
松田 洋子 / KADOKAWA / 2017-12-11

Bookcover BLUE GIANT SUPREME 4 (ビッグコミックススペシャル) [a]
石塚 真一 / 小学館 / 2018-02-23

Bookcover ペリリュー ─楽園のゲルニカ─ 4 (ヤングアニマルコミックス) [a]
武田一義,平塚柾緒(太平洋戦争研究会) / 白泉社 / 2018-02-28

Bookcover アルテ 7 (ゼノンコミックス) [a]
大久保圭 / 徳間書店 / 2017-07-20
Bookcover アルテ 8 (ゼノンコミックス) [a]
大久保圭 / 徳間書店 / 2018-01-20

コミックス(2015-) - 読了:「大砲とスタンプ」「甘々と稲妻」「大人スキップ」「BLUE GIANT SUPREME」「ペリリュー 楽園のゲルニカ」「アルテ」

Bookcover 血の轍 2 (ビッグコミックス) [a]
押見 修造 / 小学館 / 2017-12-27

Bookcover 私の少年(4) (アクションコミックス(月刊アクション)) [a]
高野 ひと深 / 双葉社 / 2017-12-12

Bookcover イムリ 22 (ビームコミックス) [a]
三宅 乱丈 / KADOKAWA / 2017-12-11

Bookcover 赤狩り THE RED RAT IN HOLLYWOOD 1 (ビッグコミックス) [a]
山本 おさむ / 小学館 / 2017-11-30

Bookcover からかい上手の高木さん 1 (ゲッサン少年サンデーコミックススペシャル) [a]
山本 崇一朗 / 小学館 / 2014-06-12
Bookcover からかい上手の高木さん 2 (ゲッサン少年サンデーコミックススペシャル) [a]
山本 崇一朗 / 小学館 / 2014-11-12
Bookcover からかい上手の高木さん 3 (ゲッサン少年サンデーコミックススペシャル) [a]
山本 崇一朗 / 小学館 / 2015-12-11
Bookcover からかい上手の高木さん 4 (ゲッサン少年サンデーコミックススペシャル) [a]
山本 崇一朗 / 小学館 / 2016-10-12
面白いマンガだけど、なんというかこの、ちょっと腹が立ってきた...

コミックス(2015-) - 読了:「血の轍」「私の少年」「イムリ」「赤狩り」「からかい上手の高木さん」

Bookcover めしばな刑事タチバナ 28 (トクマコミックス) [a]
坂戸佐兵衛,旅井とり / 徳間書店 / 2018-01-31
Bookcover めしばな刑事タチバナ 29 (トクマコミックス) [a]
坂戸佐兵衛,旅井とり / 徳間書店 / 2018-03-31

Bookcover イノサン Rouge ルージュ 7 (ヤングジャンプコミックス) [a]
坂本 眞一 / 集英社 / 2017-12-19

Bookcover はじめアルゴリズム(1) (モーニング KC) [a]
三原 和人 / 講談社 / 2017-11-22

Bookcover 中間管理録トネガワ(6) (ヤンマガKCスペシャル) [a]
福本 伸行,三好 智樹,橋本 智広 / 講談社 / 2017-11-06

Bookcover アイとアオの境界 (2) (バーズコミックス) [a]
天堂 きりん / 幻冬舎コミックス / 2017-11-24

Bookcover どこか遠くの話をしよう 下 (ビームコミックス) [a]
須藤 真澄 / KADOKAWA / 2018-01-12
正直なところちょっと体質に合わないんだけど、でも面白い漫画である。

コミックス(2015-) - 読了:「めしばな刑事タチバナ」「イノサン Rouge」「はじめアルゴリズム」「中間管理録トネガワ」「アイとアオの境界」「どこか遠くの話をしよう」

Bookcover 辺獄のシュヴェスタ 6 (ビッグコミックス) [a]
竹良 実 / 小学館 / 2017-12-12
中世ドイツの修道院を舞台にした復讐劇、最終巻。面白かった。

Bookcover サトコとナダ 2 (星海社COMICS) [a]
ユペチカ / 講談社 /

Bookcover チェイサー 5 (ビッグコミックス) [a]
コージィ 城倉 / 小学館 / 2018-01-04

Bookcover いちげき 1 [a]
松本次郎,永井義男 / リイド社 / 2016-11-30
Bookcover いちげき 2 (SPコミックス) [a]
松本次郎,永井義男 / リイド社 / 2017-07-27
Bookcover いちげき 3 [a]
松本次郎,永井義男 / リイド社 / 2018-03-27

Bookcover 大奥 15 (ヤングアニマルコミックス) [a]
よしながふみ / 白泉社 / 2017-12-28

Bookcover レベレーション(啓示)(3) (モーニング KC) [a]
山岸 凉子 / 講談社 / 2017-12-21

コミックス(2015-) - 読了:「辺獄のシュヴェスタ」「サトコとナダ」「チェイサー」「いちげき」「大奥」「レベレーション(啓示)」

長いこと記録をさぼっていたので、読了の本が溜まってしまった。

Bookcover ひなつば [a]
昌原 光一 / リイド社 / 2018-02-27

Bookcover 機械仕掛けの愛 5 (ビッグコミックス) [a]
業田 良家 / 小学館 / 2018-02-23

Bookcover めしにしましょう(3) (イブニングKC) [a]
小林 銅蟲 / 講談社 / 2017-07-21
Bookcover めしにしましょう(4) (イブニングKC) [a]
小林 銅蟲 / 講談社 / 2017-12-22
著者は前衛的ギャグマンガ家として知られていた人。この作品は一般誌連載の料理マンガなんだけど、ところどころに顔を出す異様な言語感覚が楽しい。

Bookcover 闇金ウシジマくん 40 (ビッグコミックス) [a]
真鍋 昌平 / 小学館 / 2017-07-28
Bookcover 闇金ウシジマくん 42 (ビッグコミックス) [a]
真鍋 昌平 / 小学館 / 2018-03-12
うっかり40巻を飛ばして読んでいた。

Bookcover 少女終末旅行 6 (BUNCH COMICS) [a]
つくみず / 新潮社 / 2018-03-09

Bookcover プリンセスメゾン 5 (ビッグコミックス) [a]
池辺 葵 / 小学館 / 2018-03-12

コミックス(2015-) - 読了:「ひなつば」「機械仕掛けの愛」「めしにしましょう」「闇金ウシジマくん」「少女終末旅行」「プリンセスメゾン」

相馬黒光「碌山のことなど」. 碌山美術館, 2008。
 休日に新宿三丁目付近を通りかかり、少し時間があったので、カレーの中村屋の上にあるという美術館に立ち寄って(中村つねの作品があると聞き覚えていた)、荻原守衛の裸婦のブロンズ像などをぼけーっと眺めて帰ってきたんだけど、受付前に置いてあった薄い小冊子が妙に面白そうで、つい買ってしまった。安曇野市にある碌山美術館の発行とある(碌山とは荻原守衛のこと)。中村屋創業者夫妻の妻・相馬黒光が、最晩年に荻原について語った追想の聞き書きだという(初出は昭和31年。黒光は30年に亡くなっている)。正直、ちょっと下世話な好奇心もあって。

 萩原の絶作となった上記の裸婦像について黒光はこう述べている。「足が地について立ち上がれないあの姿を見て私は正視してゐられませんでした。あの像の顔をみて山本安曇さん[工芸家、この像の鋳造者]は私に「誰の顔を作ったかわかるぢゃありませんか」と云ひました。去年も展覧会であの像を見まして、この年になっても猶息のつまる思ひをしました」
 荻原は死の病床で、黒光と戸張狐雁(相馬夫妻と荻原の友人、のちに彫刻家・画家となる)にアトリエの机の鍵を渡す。荻原の死後、二人は引き出しから日記を見つけ、読まずに焼却する。「狐雁が男泣きに泣き乍ら一枚一枚むしってはくべ、私を「ヘッダガブラーだ」と罵り乍ら灰にしたのです。狐雁はやさしい人でした」と、ここから黒光に対する萩原の秘めた愛の話になって...

 いやあ、こういっちゃなんだが、大変面白かった。相馬黒光という人は只者ではない。
 それにしても、これは映画になるんじゃないかしら、なんて思ったのだが、いま調べてみたら、2007年に信越放送が「碌山の恋」というTVドラマを制作しており、黒光は水野美紀さんが演じたそうだ。

ノンフィクション(2018-) - 読了:「碌山のことなど」

 仕事の関係でデータ解析について御相談をいただくことがあるのだけれど(嗚呼、流されていく人生...)、これこれこういうことをSPSSでやるためにはどうしたらいいでしょうか? なあんて問われることがあって、そのたびにオタオタしてしまう。SPSSなんて長らく使っていないもので。
 いやいや、SPSSの悪口を言うつもりはないです、全然。神は地上にさまざまな命をもたらし、人類にさまざまなソフトウェアをお与えになりました。いずれも等しく素晴らしい(両手を広げるポーズ)。素敵ですよね、SPSS!
 そんなこんなで、自分では普段使っていないソフトウェアについて考えるというのも、畳の上の水練みたいで、それはそれで興味深い。先日はSPSSでクラスタリングする際の事例間の近接性指標の選択について御相談いただき、かなり焦ったので、今後のためにメモを作っておく。

 手順は次の通り。

  1. SPSS Statistics Base 25の英語版マニュアルを検索し、事例間の近接性を求めていそうなプロシジャについて片っ端からメモする。BaseというのはSPSS Statisticsの基本パッケージで、基本的なプロシジャしか入っていない。オプションパッケージにも近接性指標を使うのがあるかもしれないけど、ええい、省略だ。
  2. SPSSのなにが一番困るかというと... 訂正、SPSSのどこが一番面白いかというと、統計用語の独創的な日本語訳である。というわけで、日本語版マニュアルから該当ページを探し、日本語訳を[]内に付記する。
  3. さらにさらに、アルゴリズム解説書から、各指標の定義と思われるものをメモする。なにしろ、IBMの膨大なドキュメント群からこのPDFファイルを見つけるまでが大変なのである。このPDFには日本語版がないというのがポイント。さらにいうと、このPDFには目次もノンブルもなくて... いや、そういうところも楽しいですよね!ええ楽しいです!

 さあ行きましょう! (と自分で自分を励まして)

距離
 PROXIMITIESコマンドでは以下の指標が利用できる(英語版p.43-; 日本語版p.45-; アルゴリズム解説書p.864-)。

 二値データの指標は山ほどある。アルゴリズム解説書ではいくつかのタイプにわけて説明しているので、このタイプ分けに従ってメモしておく。
 以下、2x2のクロス表について、両方presentの頻度を$a$, 両方absentの頻度を$d$, 片方だけpresentの頻度を$c,d$とする。またメモを楽にするために、$n=a+b+c+d$とする。

タイプ1. 二値データの指標 [←いや、そういってしまえばみんなそうなんですけどね]。いずれも類似性の指標。

  1. Russell and Rao[RusselとRao]: $RR(x,y) = a / n$
  2. simple matching[単純整合]: $SM(x,y) = (a+d) / n$
  3. Jaccard[Jaccard]: similarity ratioともいう。$JACCARD(x,y) = a / (a+b+c)$
  4. Dice[Dice]: Czekanowski, Sorensonともいう。$DICE(x,y) = 2a / (2a+b+d)$
  5. Sokal and Sneath 1[SokalとSneath 1]: $SS1(x,y) = 2(a+d) / (2(a+d)+b+c)$
  6. Rogers and Tanimoto[RogersとTanimoto]: $RT(x,y) = (a+d) / (a+d+2(b+c))$
  7. Sokal and Sneath 2[SokalとSneath 2]: $SS2(x,y) = a / (a+2(b+c)))$
  8. Kulczynski 1[Kulczynski 1]: $K1(x,y) = a/(b+c)$。分母が0なら9999.999とする。
  9. Sokal and Sneath 3[SokalとSneath 3]: $SS3(x,y) = (a+d)/(b+c)$。分母が0なら9999.999とする。

いま気が付いたんだけど、普段頼りにしている齋藤・宿久(2006)「関連性データの解析法」(共立出版)では、JaccardとRussell & Raoの定義が入れ替わっている(p,18-19)。あらら。おそらく単純な誤植だろう。

タイプ2. 条件付き確率。いずれも類似性の指標である。

  1. Kulczynski 2[Kulczynski 2]: $K2(x,y) = (a/(a+b) + a/(a+c))/2$. 一方がpresentであるという条件の下での、他方がpresentである条件付き確率の平均。
  2. Sokal and Sneath 4[SokalとSneath 4]: $SS4(x,y) = (a/(a+b)+a/(a+c)+d/(b+d)+d/(c+d))/4$。一方がなにかであるという条件の下での、他方がそれと同じである条件付き確率の平均。
  3. Hamann[Hamann]: $HAMANN(x,y) = ((a+d)-(b+c))/n$. 値が同じである確率から、値が異なる確率を引いたもの。

うーん、Hamannの指標を条件付き確率と呼ぶことができるのだろうか? ま、条件付き確率に基づく指標ではある。

タイプ3. 予測可能性の指標。いずれも類似性の指標である。あらかじめ次のように定義しておく。
 $t_1 = max(a,b)+max(c,d)+max(a,c)+max(b,d)$
 $t_2 = max(a+c, b+d) + max(a+b, c+d)$

  1. Lambda[ラムダ]: $LAMBDA(x,y) = (t_1 - t_2)/(2n-t_2)$。Goodman-KruskalのLambda。一方の値が分かったことで、他方の値についての予測誤差が何割下がるかという指標。
  2. Anderberg's D[AnderbergのD]: $D(x,y) = (t_1 - t_2)/2n$。Lambdaと似ているが、実際の予測誤差の減少を測っている。
  3. Yule's Y[YuleのY]: $Y(x,y) = (\sqrt{ad}-\sqrt{bc})/(\sqrt{ad}+\sqrt{bc})$
  4. Yule's Q[YuleのQ]: $Y(x,y) = (ad-bc)/(ad+bc)$. Goodman-Kruskalのgammaの2x2バージョン。

タイプ4. その他(類似性)

  1. Ochiai[落合]: $OCHIAI(x,y) = \sqrt{(a/(a+b)) (a/(a+c))}$
  2. Sokal and Sneath 5[SokalとSneath 5]: $SS5(x,y) = ad/\sqrt{(a+b)(a+c)(b+d)(c+d)}$
  3. phi 4-point correlation[ファイ4分点相関]: $PHI(x,y) = (ad-bc)/\sqrt{(a+b)(a+c)(b+d)(c+d)}$
  4. dispersion[散らばり]: $DISPER(x,y) = (ad-bc)/n^2$

タイプ5. その他(非類似性)

  1. Euclidean distance[ユークリッド距離]: $BEUCLID(x,y) = \sqrt{b+c}$
  2. squared Euclidean distance[平方ユークリッド距離]: $BSEUCLID(x,y) = b+c$
  3. size difference[サイズの差]: $SIZE(x,y) = (b-c)^2/n^2$
  4. pattern difference[パターンの違い]: $PATTERN(x,y) = bc/n^2$
  5. shape[形]: $BSHAPE(x,y) = (n(b+c)-(b-c)^2)/n^2$
  6. variance[分散]: $VARIANCE(x,y) = (b+c)/4n$
  7. Lance and Williams[LanceとWilliams]: Bray-Curtisともいう。$BLWMN(x,y) = (b+c)/(2a+b+c)$

階層クラスタリング
 Hierarchical Cluster[階層クラスター]プロシジャの場合 (CLUSTERコマンドのこと。英語版p.87-;日本語版p.93)... とメモしようと思ったら、順序は違えど上記のPROXIMITIESコマンドと同じであった。
 ということは結局、SPSS Statisticsで二値データを階層クラスタリングしようとするとと、豊富な距離指標を選べてしまうわけだ。そのいっぽう、IBMのサポートページには「二値データで階層クラスタリングするのはやめろ」というアドバイスが載っているわけで... 落とし穴掘ったから落ちないようにね、と親切に教えてもらったような気分である。

最近隣法
 Nearnest Neighbor Analysis[最近傍分析]プロシジャでは、以下の指標が利用できる(KNNコマンド。英語版p.66; 日本語版p.70; アルゴリズム解説書p.578)。

  1. Euclidean metric[ユークリッド計量]: $Euclidean_{ih} = \sqrt{\sum_p w_p (x_{pi}-x_{ph})^2}$
  2. City block metric[都市ブロック計量]: $CityBlock_{ih} = \sum_p w_p |x_{pi} - x_{ph}|$

なお、$w_p$は通常は$1$、特徴の重要度を使って重みづけする場合には、特徴の重要度を合計1に規準化した値になる。

TwoStepクラスタリング
 TwoStep Cluster[TwoStepクラスタ]プロシジャでは、クラスタ間距離として以下が選べる(TWOSTEP CLUSTERコマンド。英語版p.79; 日本語版p.84。定義がアルゴリズム解説書p.1078-にあるけど、そこだけメモしてもしょうがないので割愛する)。

  1. Log-likelifood[対数尤度]
  2. Euclidean[ユークリッド]

K-Means法
 K-Means Cluster[大規模ファイルのクラスター分析]プロシジャでは、距離はユークリッド距離に固定されている。(QUICK CLUSTERコマンド)

MDS
 Multidimensional Scaling[多次元スケーリング法]プロシジャの場合、以下の指標が選択できる。(ALSCALコマンド。英語版p.124; 日本語版p.133。さすがに訳語はPROXIMITIESコマンドでの訳語と揃えているようなので、略記する。アルゴリズム解説書p.14-を見ると、重みつきユークリッド距離の定義しかないんだけど、なぜだろう? なにか見落としているのかもしれない。)

雑記:データ解析 - 覚え書き:SPSS Statisticsにおける近接性の指標

Carvalho, A. (2016) A note on Sandroni-Shmaya belief elicitation mechanism. The Journal of Prediction Markets, 10(2), 14-21.
 ほんの出来心で目を通したSandroni & Shmaya(2013)がよく分からなかったので、毒を食らわば皿までという気分で(すいません)、こちらもめくってみた。短いし。
 Google Scholar上での被引用件数は... 1件だ。いやあ、風情があるなあ。

 以下にメモを取るけど、原文にある$x_{max}, x_{min}$は書きにくいので、Sandroni-Shmayaにあわせて$x, y$に書き換える。要するに、クジがあたったときのペイオフ金額と、はずれたときのペイオフ金額のことである。

 網羅的かつ相互排他なアウトカムを$\theta_1, \ldots, \theta_n$とする。専門家はアウトカムについての真の信念$p=(p_1, \ldots, p_n)$を持っている。彼が報告する信念$q=(q_1, \ldots, q_n)$を、$p=q$となるようにしたい。
 そのための伝統的テクニックとしてプロパー・スコアリング・ルールがある。アウトカム$\theta_x$が観察されたら専門家をスコア$R(q, \theta_x)$で評価し、これに応じてなんらかの報酬を渡す。$p=q$のときそのときに限りスコアが最大化されるとき、これをプロパーという。もっともポピュラーなプロパー・スコアリング・ルールとして、対数ルール$R(q, \theta_x) = log(q_x)$, 二次ルール$R(q, \theta_x) = 2q_x - \sum_k^n q_k^2$がある。
 プロパー・スコアリング・ルールは専門家がリスク中立だという仮定に基づいている。しかし、たとえばリスク志向的な専門家は信念をシャープに報告しがちになるだろう。

 リスク中立性が維持できない場合の方法として、効用関数を$U(\cdot)$としてスコアを$U^{-1}(R(q, \theta_x))$とする方法がある(Winkler, 1969 JASA)。[つまり効用がR(q, \theta_x)となるようにあらかじめ変換しておくということね。Winklerってひょっとして、おととしベイジアン合意について調べていたときに出てきた、あのWinklerさん? 世間狭いなあ...]
 このアプローチは次の2つの条件に依存する。(1)専門家の振る舞いが、期待効用理論に従って既知。(2)専門家の効用関数$U(\cdot)$が既知。これらの条件には無理がある。

 そこで登場する回避策が、あらかじめ専門家のリスク態度を規定する要素を調べておこうというもの。そうすれば、それらの要素の影響を取り除くことで、専門家の報告を事後的にキャリブレートできる。
 この路線においても、専門家がなんらかの決定モデルに従って振る舞うという仮定が必要になる。決定モデルが誤っていたらおじゃんである。

 そこでいよいよ登場するのが、支払いをクジで決めるという路線である。
 Allen(1987 MgmtSci)は、専門家の効用関数が未知である場合に、条件つきクジをつかった効用の線形化によって誠実な報告を引き出すという手法を提案した。
 またKarni(2009, Econometrica)は、金額を2つに固定し、専門家が報告した確率を[0,1]の一様乱数と比較することで支払関数を決めるという方法を提案した。この方法だと、専門家がprobablistic sophistication and dominanceを示すなら、リスク態度と関わりなく、誠実な申告が専門家にとって最適になる。
 AllenとKarniのアプローチは、考え方として古典的なBDMメカニズムと似ている。実験場面ではBDMメカニズムにうまく対処できない被験者がいることが知られている(Cason & Plott 2014 J.PoliticalEcon.; Plott & Zeiler 2005 Am.Econ.Rev.; Rutstrom 1998 IntJ.GameTheory)。

 お待たせしました。本論文の主役、Sandroni & Shmaya (2013)の登場です。
 信念報告というのはクジの選択みたいなものである。$n=2$の場合について考えよう。$q=(q_1, q_2)$と報告するということは、
 $[R(q, \theta_1):p_1, R(q, \theta_2):p_2 ]$
というクジを選んだのとおなじことである。
 報酬$x > y$, 確率$0 \leq \rho, \rho' \leq 1$として、
 クジA: $[y:\rho, x:1-\rho]$
 クジB: $[y:\rho', x:1-\rho' ]$
があるとき、$\rho < \rho'$のときそのときに限りBよりAが選好される、ということをprobabilitic dominanceという。
 Sandroni & Shmaya (2013)の主張は次の通り。proabilistic sophisticationは成り立っているとする[←これについては説明がない...]。誠実な信念報告を引き出すためには、probabilistic dominanceさえ成り立っていればよい。
 彼らが提案した支払スキーマはこうである。[0,1]に規準化されたプロパー・スコアを$S(q, \theta_x)$とする。(1)アウトカム$\theta_1$が起きたら、確率$S(q, \theta_1)$で$x$を払い、確率$1-S(q, \theta_1)$で$y$を払う。(2)アウトカム$\theta_2$が起きたら、確率$S(q, \theta_2)$で$x$を払い、確率$1-S(q, \theta_2)$で$y$を払う。

 この提案のポイントは、専門家のリスク態度についてなにも仮定していないという点である。なお、BDMメカニズムとの違いは、支払い決定にあたって外部のランダム化装置が不要であるという点である[原文: "This mechanism differes from Becker-DeGroot-Marschak based mechanisms in that no external randomization device other tha nature is required to determine an expert's payment." よくわからない。BDMメカニズムでもSandroni & Shmayaでも、ペイオフ決定にあたってはなんらかの確率乱数の生成が必要じゃないの?]

 さて、Sandroni & Shmaya が見落としている点がある。
 彼らは、専門家からみて、次の2種類のクジが等しいと考えている。
 クジその1、上述のクジ。(1)アウトカム$\theta_1$が起きたら、確率$S(q, \theta_1)$で$x$を払い、確率$1-S(q, \theta_1)$で$y$を払う。(2)アウトカム$\theta_2$が起きたら、確率$S(q, \theta_2)$で$x$を払い、確率$1-S(q, \theta_2)$で$y$を払う。
 クジその2。アウトカム$\theta_1, \theta_2$の真の主観確率を$p_1, p_2$とする。確率$p_1 S(q, \theta_1)+ p_2 S(q, \theta_1)$で$x$を払い、確率$p_1 (1-S(q, \theta_1))+ p_2 (1-S(q, \theta_1))$で$y$を払う。
 この2つのクジが等しいというのは、自明ではない。Sandroni & Shmayaは暗黙のうちに、合成くじの分解(reduction of compound lotteries; ROCL)の定理を仮定しているのである。

 残念ながら、現実にはROCL定理は必ずしも維持されない。Harrison, et al.(2014 J.Econ.Behav.Org)をみよ。彼らによれば、選択が二値の場合はROCLは維持されるが、多値の場合には維持されない。
 では、アウトカムが二値であればSandroni & Shmayaのメカニズムは真実報告を引き出すか。そうかもしれないし、そうでないかもしれない。信念報告をクジの選択と捉えた時、$q=(q_1, q_2)$と報告するということは、(2個ではなくて)無限個のクジのなかから
 $[R(q, \theta_1):p_1, R(q, \theta_2):p_2 ]$
を選ぶということだからである。
 さらにいえば、Harrisonらの実験はすべての専門家の信念が互いに等しいような場面でのものであって、一般的な不確実性についていえるのかどうかはオープン・クエスチョンである。

 ...ふうん。
 知識が足りなくて、このモヤモヤをうまく表現できないんだけど... 「人にはリスク選好ってのがある」という指摘と、「人は必ずしもROCL定理に従わない」という指摘は、同じレイヤの話なんだろうか?
 というわけで、この論文の意義については判断できないけど、先行研究概観は勉強になりました。

論文:予測市場 - 読了:Carvalho (2016) 当たり外れがプロパー・スコアリング・ルールで決まるクジを報酬にすれば参加者のリスク選好がどうであれ真実申告メカニズムが得られるというのは本当か

Sandroni, A., Shmaya, E. (2013) Eliciting beliefs by paying in chance. Economic Theory Bulletin, 1, 33-37.
 昨年のJ. Prediction Marketsにこの論文についてのコメントが載っていて、どうやらベイジアン自白剤のことが引き合いに出されているらしいので、読んでみた。
 雑誌名からして素人が読むべきものではないのかもしれないけれど、たったの5pだし、数式も少ないので、試しに目を通してみた次第。どうせ何についても専門家じゃないんだから、いいじゃないですか、何を読んだって。
 google scholarによる被引用回数は... 6件。渋い。

 いわく。
 専門家に自分の主観確率を誠実に報告させるインセンティブを求める方法としてプロパー・スコアリング・ルールがある。そういうのは多くの場合、専門家がリスク中立であること(ないしリスク選好が既知であること)を仮定している。本論文では非常に単純な原理を述べる。この原理を使えば、専門家の選好が既知であるという想定をdisposeすることができる。

 これから出来事$E$が起きるかもしれないし起きないかもしれない。ボブは$E$についての主観確率を持っている。Bobの主観確率をどうやって引き出すか。
 ボブが確率$\hat{p}$を申告したとして、$E$が起きたら金銭報酬 $S(\hat{p}, 1)$, 起きなかったら$S(\hat{p}, 0)$を支払うとしよう。ここで
 $S(\hat{p}, 1) = 2 - (1-\hat{p})^2$
 $S(\hat{p}, 0) = 2 - (\hat{p})^2$
とする支払スキーマ$S$をBrierスコアという。これはプロパー・スコアリング・ルールの例である。ペイオフの期待値
 $p S(\hat{p}, 1) + (1-p) S(\hat{p}, 0)$
を最大化するのは$\hat{p} = p$なので、ボブがリスク中立的なら、ボブは自分の主観確率を申告する。

 問題は、ボブのリスク選好がわからない場合、つまり$S(\hat{p}, 1), S(\hat{p}, 0)$がボブにとっての効用なのかどうかわからない場合である。
 ひとつの路線は、別の実験をやってボブの選好を調べるというものである。いっぽう、Karni(2009)は別の路線を考えた。以下で説明しよう。なお、より包括的な定式化としてLambert(2011)がある。
 金銭報酬$x, y$($x > y$)を使った次の2つのクジがあるとしよう。
 A: 確率$\mu$で$x$ドルもらえ、確率$1-\mu$で$y$ドルもらえるクジ。
 B: 確率$\mu'$で$x$ドルもらえ、確率$1-\mu'$で$y$ドルもらえるクジ。
ボブがBよりAを選好するのは、$\mu > \mu'$のとき、そのときに限ると仮定する。この仮定をprobabilistic sophistication and dominanceという。
 さて、次のランダム・スコアリング・ルールを考える。$P(\hat{p},1) = S(\hat{p},1)/2, P(\hat{p},0) = S(\hat{p},0)/2$とし、$E$が起きたら「確率$P(\hat{p},1)$で$x$ドルもらえ、確率$1-P(\hat{p},1)$で$y$ドルもらえるクジ」、$E$が起きなかったら「確率$P(\hat{p},0)$で$x$ドルもらえ、確率$1-P(\hat{p},0)$で$y$ドルもらえるクジ」を渡すのである。2で割っているのは単に確率を0~1の範囲に収めたいから。
 このとき、ボブが$x$を得る確率は
 $\{p S(\hat{p}, 1) + (1-p) S(\hat{p}, 0)\}/2$
なので、$S$がプロパーであれば、$p$を申告するのが最適となる。これはボブのリスクへの態度に関わらず成り立つ。
 [うううううう... わからないいいい... なぜそういえるの...??? リスク中立でないってことは、$S(\hat{p}, 1), S(\hat{p}, 0)$が効用でなくて、$U(S,p)$というような形の効用関数が別にあるってことだよね? それがどういう形であれ、sophisticationとdominanceという条件を満たしていれば、$\hat{p}=p$が効用を最大化するといえる、ってこと??? それって自明なの? どうも話の肝になるところが理解できていないみたいだ...]

 上の例は、支払額を偶然で決めることによって信念を引き出すという原理を示している。基本手続きは以下の通り。(1)なんらかのプロパー・スコアを基準化して、スコアが0~1に入るようにする。(2)この基準化されたプロパー・スコアを確率とみなし、この確率で高いほうの金銭報酬$x$を渡す。要するに、スコアが高いとき、高い報酬がもらえる確率が高くなるわけである。

 この原理は、多エージェントのゲーム理論的セッティングにも使える。たとえばPrelec(2004)のベイジアン自白剤について考えよう。[以下、ややこしいので全訳する]

Prelecは、それぞれの専門家の意見は共通の分布を持つある確率変数の実現値だと想定した。彼は次のようなベイジアン・ゲームを設計した。そのゲームにおいて、専門家$i$の行為空間$A_i$は、yesかnoかを述べること、そしてyesと答えた専門家の割合を予測することである。プレイヤー$i$の純戦略は、彼の実際の意見(yesないしno)を行為空間$A_i$にマップするものである。プレイヤー$i$のペイオフは、彼が構築する具体的な効用関数
 $U_i: \prod_k A_k \rightarrow R$
で与えられる。彼の論文のキーポイントは、全ての専門家たちが自分の意見を誠実に申告するナッシュ均衡が存在するというものである。
 Brierスコアの場合のように、$U_i$の下でのペイオフが金銭を単位として与えられている場合、そこではリスク中立性が仮定されている。ペイオフが効用であると仮定されている場合、そこでは専門家のリスク態度が既知だと仮定されている。しかし、いま線形の狭義単調増加関数$\tau$があり、全てのプレイヤー$i$、すべての行為プロフィール$a \in \prod_l A_k$について、$\tau(U_i(a))$が0と1の間だとしよう。行為プロフィールが$a \in \prod_l A_k$であるゲームにおいて、プレイヤー$i$が確率$\tau(U_i(a))$で金銭報酬$x$を承けとり、確率$1-\tau(U_i(a))$で金銭報酬$y$を受け取るとする(ただし$x > y$)。$x$と$y$のどちらになるかを決めるランダム化は、本質的に、それぞれのエージェントからもそれぞれの行為プロフィールからも独立である。Prelecのゲームにおける真実申告ナッシュ均衡は、この修正されたゲームにおいてもやはりナッシュ均衡である。このことは、probablistic sophistication and dominanceの下で、専門家のリスク態度と無関係に成立する。

 ...忘れちゃったんだけど、ベイジアン自白剤ってプレイヤーのリスク中立性を仮定しているんだったっけか。あとで調べておこう。
 えーと、要するに、報酬をプロパー・スコアリング・ルールで与えたときは、参加者がリスク中立でないと真実申告メカニズムにはならないんだけど、報酬を「プロパー・スコアリング・ルールに基づく確率」で決めれば、リスク選好がどうであれ真実申告メカニズムが作れるんだよ、という話なんだと思う。で、それはベイジアン自白剤にもあてはまるんだよ、ということなんだと思う。そうなんすか。
 ときに、報酬を確率的に決めるというのはBDMメカニズムもそうなんだけど、どういう関係にあるんだろうか。

論文:予測市場 - 読了:Sandroni & Shmaya (2013) 当たり外れがプロパー・スコアリング・ルールで決まるクジを報酬とせよ、さすれば参加者のリスク選好がどうであれ君は真実申告メカニズムを得るだろう

2018年5月 2日 (水)

Jennings, W., Wlezien, C. (2018) Election polling errors across time and space. Nature Human Behavior.
 ネットで話題になってたので目を通してみた。NatureとかScienceとかの記事って、構成が普通の論文と違ってて、ちょっと読みにくいんですよね...

 いわく。
 2015年UK総選挙、2016年US大統領選挙をみるにつけ、選挙調査は誤差が増大しているのではないか、危機にあるのではないか...という声が高まっている[←ネイト・シルバーのブログ記事やNYTの記事がreferされている]。ほんとうだろうか。
 1942年から2017年までの45ヶ国の国政選挙351件について分析する(地方選や国民投票は対象外。うち日本が5件)。分析対象となる調査は実に30916件。そういうすごいデータベースがあるのです(Harvard Dataverse)。
 分析に使用するのは、各調査における候補ないし政党への投票意向の数値である(DKと回答拒否を除いたシェア)。各調査における実査日を特定しておく(複数日にまたがる場合には中央とする)。各選挙について、候補・政党への投票シェアを調べておく。
 目的変数は投票-調査間の絶対誤差。必要に応じ、投票-調査間の対数オッズ比とか、首位2つの合計の間での絶対誤差とかも調べる。

 考察するに、回収率は下がり、調査モード・抽出・ウェイティングはますます多様になったけど、そのせいで選挙前調査のパフォーマンスは下がってはいない。おそらく、異なる調査が結合されているせいでキャンセルアウトされたのだろう(実際、誤差の調査間分散は大きくなっている模様)。
 本研究の限界:誤差に影響する要因は他にもあろう(調査結果公表に禁止期間を設けている国があるとか、制度変更とか、有権者それ自体のボラタリティとか、調査モードとか)。
 
 というわけで、調査機関はさまざまな難題に直面しているけれど、選挙調査の正確性が危機にあるとはいえない。
 選挙調査は大きくはずれることもあるだろう、注目を集めることもあるだろう。しかあし!それは方法論的反省、イノベーション、そして改善を導きうるのだっ。
 云々。

論文:調査方法論 - 読了:Jennings & Wlezien. (2018) 世界の選挙調査は危機にあるか?→意外にそうでもない

Trajman, A., & Luiz, R.R. (2008) McNemar chi^2 test revisited: comparing sensitivity and specificity of diagnostic examinations. The Scandinavian Journal of Clinical & Laboratory Investigation. 68(1), 77-80.
 McNemar検定(対応のある2条件間の割合の差の検定)についてちょっと知りたいことがあり、タイトルに惹かれて目を通した。なんだかマイナーそうな誌名だが、Google Scholar上は被引用100件。
 新しい診断基準ができたとき、フォローアップで患者と分かった人だけに絞って、通常の基準での診断を表側、新基準での診断を表頭にとった2x2クロス表をMcNemar検定すれば、敏感度の差を検定できますね、という論文だった。
 メモを取りながら読んだんだけど、そりゃそうだろうなという感じで、正直言って話のポイントが全然つかめなかった...。もっとも、読みながら何度も意識を失ったので、なにか読み落としているのかも。ま、少なくとも知りたい内容についての論文ではなかった、いいや、次にいこう。

論文:データ解析(2018-) - 読了:Trajman & Luiz (2008) 二つの診断基準の敏感度と特異度をMcNemar検定で比較する

 大きなことから小さなことまで、悩みの種は尽きないが、そのなかには勉強不足に起因するとしかいいようがない問題もある。もっとも、勉強すりゃいいんだから、ある面では気楽である。
 予測モデルの性能評価というのもそうした問題のひとつである。というわけで、数か月前に仕事を中断してメモをつくった。
 以下はHastie, Tibshirani, Friedman「統計的学習の基礎」第7章「モデルの選択と評価」からのメモ。英語版なら無料なのに高価な訳書を買っちゃったので、せいぜい頑張ってモトをとらねばならぬ。

 まずは損失関数loss functionの定義から。
 目標変数を$Y$, 入力ベクトルを$X$、予測モデルを$\hat{f}(X)$、訓練データを$T$とする。訓練データにおける損失関数としては、$Y$が量的ないし間隔の場合は、
 二乗誤差 $L(Y, \hat{f}(X)) = (Y - \hat{f}(X))^2$
 絶対値誤差 $L(Y, \hat{f}(X)) = |Y-\hat{f}(X)|$
が用いられることが多い。$Y$が質的な場合は、水準数を$K$とし、モデルを$p_k(X) = Pr(Y=k|X)$、予測を$\hat{Y}(X) = argmax_k \hat{p}_k(X)$として、
 0/1損失 $L(Y, \hat{Y}(X)) = I(Y = \hat{Y}(X))$
 逸脱度 $L(G, \hat{p}(X)) = -2 \sum_k^K I(Y=k) \log \hat{p}_k(X)$
が使われることが多い。確率密度に関する損失関数としては
 対数尤度 $L(Y, \theta(X)) = -2 \log Pr_{\theta(X)}(Y)$
が用いられる。
 標本における損失の平均
 $\bar{err} = (1/N) \sum_i^N L(y_i, \hat{F}(x_i))$
訓練誤差training errorという。

 いっぽうテスト誤差test error(=汎化誤差 generalization error)は
 $Err_T = E[L(Y, \hat{f}(X)) | T]$
として定義される。その期待値、すなわち訓練データのランダム性を含むすべてのランダム性についての平均
 $Err = E[L(Y, \hat{f}(X))]$
期待テスト誤差 (期待予測誤差) という。ほんとにやりたいことは、手元の$T$についての$Err_T$の推定なんだけど、実際には$Err$を推定するわけである。
 モデルに対する期待テスト誤差は、モデルが含んでいるパラメータ$\alpha$まで含めて $Err = E[L(Y, \hat{f}_\alpha(X))]$
と書くべきところだが、以下では略記する。

 さて、(期待)テスト誤差を最小にするモデルを選ぶにはどうしたらよいか(モデル選択)。また、最終的なモデルの(期待)テスト誤差を推定するにはどうしたらよいか(モデル評価)。というのがここからのテーマである。

 もしもデータが十分にあるならば、データを訓練集合、確認集合、テスト集合に分割するのがよい。確認集合はモデル選択に使い、テスト集合は最終的な(期待)テスト誤差の推定だけに使う。理想のサイズは一概にいえないが、たとえば50:25:25という感じ。
 もっとも、データはふつう十分ではないので、なにかと大変なのである。

 ここでちょっと脇道にそれて、バイアス-分散分解について。
 $Y=f(X)+e, E(e)=0, Var(e)=\sigma_e^2$と仮定しよう。$f(X)$は回帰だとして、2乗誤差損失を採用しよう。入力点$X=x_0$における$f(X)$の期待テスト誤差は
 $Err(x_0)$
 $= E[(Y-\hat{f}(x_0))^2|X=x_0]$
 $= [E\hat{f}(x_0)-f(x_o)]^2 + E[\hat{f}(x_0)-E\hat{f}(x_0)]^2 + \sigma_e^2$
と分解できる。第1項がバイアスの二乗、第2項が分散を表す。
 もっともこれは二乗誤差損失の場合の話。バイアス-分散トレードオフの挙動は、損失をどう定義するかで変わってくる。

 テスト誤差をもう一度きちんと定義しよう。
 データの同時分布から生成された新しいテストデータ点$(X^0 ,Y^0)$におけるモデル$\hat{f}$のテスト誤差を
 $Err_T = E_{X^0, Y^0}[L(Y^0, \hat{f}(X^0)) | T]$
とする。期待テスト誤差は、これを訓練集合$T$について平均した値
 $Err = E_T E_{X^0, Y^0}[L(Y^0, \hat{f}(X^0)) | T]$
である。
 ご存知のように、訓練誤差
 $\bar{err} = (1/N) \sum_i^N L(y_i, \hat{F}(x_i))$
はテスト誤差の推定値としては楽観的である。
 いま、訓練入力点$x_1, \ldots, x_N$に対して新しい応答変数$Y^0$が与えられたとしよう。
 $Err_{in} = (1/N) E_{Y^0}[L(Y^0, \hat{f}(X^0)) | T]$
訓練標本内誤差in-sample errorと呼ぶことにする。これ自体に関心はないんだけど(だってテスト入力点は訓練入力点とは異なるから)、話の都合上これを考えるわけである。
 $\bar{err}$は$Err_{in}$に対してもなお楽観的である。そこで最善度optimismを
 $op \equiv Err_{in} - \bar{err}$
と定義する。で、訓練集合上での予測変数のほうは固定して、出力変数にたいして期待値をとり[←ここの発想がややこしい...]
 $\omega \equiv E_y(op)$
と定義する。
 すると、幅広い損失関数について一般に
 $\omega = (2/N) \sum_i^N Cov(\hat{y}_i, y_i)$
が成り立つ。つまり、$y_i$が予測値に強く効いているとき、$\bar{err}$はより楽観的なわけだ。

 さて、話を戻して... (期待)テスト誤差をどうやって推定するか。路線がふたつある。
 ひとつの路線は、どうにかして最善度$\omega$を推定するというもの。それを訓練誤差$\bar{err}$に乗せてやれば、すくなくとも訓練標本内誤差$Err_{in}$は推定できるわけで、モデル間の相対的な比較はできるわけである。
 では、$\omega$を具体的にどうやって推定するか。
 単純なケースとして、$\hat{y}_i$が$d$次元の線形フィッティングで得られている場合を考えよう。$Y=f(X)+e$において$\sum_i^N Cov(\hat{y}_i, y_i) = d \sigma_e^2$だから、
 $\omega = 2 \cdot (d/N)\sigma_e^2$
ここから得られるのが$C_p$統計量
 $C_p = \bar{err} + 2 \cdot (d/N) \sigma_e^2$
である。
 これをもっと一般的にしたのがAICで...[ここから、AIC, 有効パラメータ数、BIC, MDL, VC次元の話が続くんだけど、ちょっと時間がないのでパス]

 もうひとつの路線は、$Err$そのものを直接に推定しようというものである。
 この路線としては、まずは交差確認(CV)が挙げられる。
 推定対象について誤解しないように。leave-one-out(LOO)とかだとついつい$Err_{T}$を推定しているような気がしちゃうけど、あくまで$Err$を推定してるんだからね。
 K-fold CVの$K$をどう決めるか。$K=N$の場合(つまりLOOの場合)、$N$個の訓練集合は互いに似ているから、気が付かないけど分散が大きくなる。いっぽう$K$が小さいと、訓練曲線によっては$Err$を過大評価するかもしれない。折衷案として5とか10とかを推奨する。
 CVでは「1SE」ルールが用いられることが多い。これは、最良モデルの誤差とSEを求め、最良モデルより1SE以内のモデルのなかでもっとも簡素なモデルを選ぶというルール。
 [ここで一般化CVの話。パス]
 予測変数がたくさんある分類問題で、(1)予測変数を相関とかで絞り込み、(2)それらで分類器を組み、(3)CVでパラメータを推定して最終モデルにする、というやりかたをとる場合がある。これはよくある間違い。最初の選別が教師無しである場合を例外として、とにかく最初から標本を抜くのがポイント。(1)標本をKグループにわけ、(2)各群について、それを抜いた全群で予測変数簿を絞り込み分類器を組み立て、その群でテストする、というのが正しいやりかた。
 分類誤差を最小にするような1変数を選ぶという場面を考える(こういうのをstumpという)。「すべての訓練集合に対してあてはめを行うと、データをとてもよく分ける予測変数をひとつみつけることになる。よってCVでは誤差を正確に推定できない」...これも間違い。CVの繰り返しのたびにモデルを完全に再学習しないといけないということを忘れている。

 もうひとつはブートストラップ法。訓練集合から同じサイズの標本集合を復元抽出で取り出す。これをB回繰り返す。B個の標本のそれぞれから任意の量を求める。B個の値の分布を調べる。
 予測誤差の推定という文脈ですぐに思いつくのが、各ブートストラップ標本にモデルをあてはめて、予測値$\hat{f}^{*b}(x_i)$を得て、
 $\hat{Err}_{boot} = (1/B)(1/N)\sum_b^B \sum_i^N L(y_i, \hat{f}^{*b}(x_i))$
を得るというのがひとつのアプローチ。しかしこれは$Err$の良い推定値ではない。訓練集合とテスト集合が分かれていないから。
 そうではなくて、個々の観察値について、それを含まないブートストラップ標本だけを使って予測値を得る、というのが正しい。観測値$i$を含まないブートストラップ標本の集合を$C^{-i}$, そのサイズを$|C^{-i}|$として
 $\hat{Err}^{(1)} = (1/N) \sum_i (1/|C^{-i}|)\sum_{b \in C^{-1}} L(y_i, \hat{f}^{*b}(x_i))$
とするわけである。これがLOOブートストラップ誤差である。
 ところが、CVの場合と同じく、訓練集合が小さくなったことに寄るバイアスが生じる。
 バイアスの大きさはどのくらいか。ある観測値$i$があるブートストラップ標本$b$に含まれる確率は$1-(1-(1/N))^N$。およそ0.632である。つまり、それぞれのブートストラップ標本における「異なり観測値」の数は、平均しておよそ$0.632N$。これはCVでいえば$K=2$の場合に相当する。だから、訓練曲線の傾きが$N/2$のあたりでまだ大きいようであれば、LOOブートストラップ標本による推定は過大になる。
 
 バイアスを修正するにはどうすればよいか。ふたつのアイデアがある。
 その1,0.632推定量。
 $\hat{Err}^{(0.632)} = 0.368 \bar{err} + 0.632 \hat{Err}^{(1)}$
となる[この導出は面倒らしい]。この推定量は過学習のことを全然考えていない。
 その2,0.632+推定量。まず、$y_i$と$x_{i'}$の全ての組み合わせを通じて
 $\hat{\gamma} = (1/N^2) \sum_i \sum_{i'} L(y_i, \hat{f}(x_{i'}))$
を求め、
 $\hat{R} = (\hat{Err}^{(1)} - \bar{err}) / (\hat{\gamma} - \bar{err})$
を求める。これを相対過学習率と呼ぶ。で、
 $\hat{\omega} = 0.632 / (1-0.368 \hat{R})$
 $\hat{Err}^{(0.632+)} = (1-\hat{\omega}) \bar{err} + \hat{\omega} \hat{Err}^{(1)}$
を求める。なぜなら...[ちゃんと丁寧に説明されているのだが、ここで力尽き、頭が停止した。いずれまた気力のあるときに...]

 最後に、CVで得た期待テスト誤差$Err$の推定は、我々が本当に知りたい対象であるテスト誤差$Err_T$の推定になっているか、という問題。
 実験してみるとわかるんだけど[...中略...]、意外にうまくいかない。訓練集合$T$に対するテスト誤差の推定は、同じ訓練集合からのデータだけではわからない。

 嗚呼... いろいろ勉強不足を痛感する次第である。LOOブートストラップってそういう意味だったのか。
 いや、勉強不足なのはまだよい。本当の問題は、このトシになって勉強してても無駄なんじゃないかという絶望感がひたひたと押し寄せてくる、という点である。辛いなあ。

雑記:データ解析 - 覚え書き:教師あり学習モデルの評価方法

2018年5月 1日 (火)

Choi, S.S., Cha, S.H., Tappert, C.C. (2010) A survey of binary similarity and distance measures. Journal of Systemics, Cybernetics and Informatics, 8(1), 43-48.
 6頁の短い論文。掲載誌についてはよくわからない。
 何の気になしに印刷し、帰宅する電車で眺めていたらこれが妙に面白く、図を「へぇぇぇ...」と感心して眺めていたら乗り過ごしてしまった。仕方がないので次の駅で降り、駅前の深夜スーパーでうろうろして安売りのワインを買った。

 どういう話かというと、長さが等しいふたつの二値ベクトルの類似性を測ることがありますわね。値が一致する個数を調べるとか、値が異なる個数を調べるとか。そういうときに使う指標が山のようにある。そこで、どの指標とどの指標が似ているかというデンドログラムを作りましたという話である。類似性指標の類似性を調べているわけです。

 要するに、なんらかのランダムデータを使って、二値ベクトル間の類似性をいろんなやり方で測ったところ、いつも似た値になる指標とそうでない指標があった、これをデンドログラムで表現しました、と。
 面白いのは、検討する類似性指標のコレクションである。集めも集め、実に76種類もの指標について調べている。各指標の発表年表もついていたりして、眺めているだけで飽きない。いやあ、いろんな指標があるものだ...

 著者らはデンドログラムに基づき、いくつかの指標のグループを指摘しているので、メモしておく。以下、2x2のクロス表について、両方presenceを$a$, 両方absenceを$d$, 片方だけpresenceを$b, c$とする。合計を$n=a+b+c+d$とする。

ほかにも一匹狼的な指標がある。Yuleのw $(\sqrt{ad}-\sqrt{bc})/(\sqrt{ad}+\sqrt{bc})$ とか。

 ... いやー、週末にできちゃうような実験だけど(すいません)、面白かった。
 一番面白かったのは指標の年表である。Jaccard類似性は1901年発表だが、これより古い(つまり19世紀の)指標として、二値ユークリッド距離、Peirceの類似性$(ab+bc)/(ab+2bc+cd)$、Yuleの類似性(YuleのQ$(ad-bc)/(ad+bc)$のことかな?)が挙げられるのだそうだ。ハミング距離なんていったらえらく古めかしい印象があるけど(情報理論っていうんでしょうか?)、発表は1950年、統計学の長い歴史の中では中堅どころに過ぎないのであった。恐れ入りました。
 仕事に役立つ学びは、というと...正直あんまりないけれど、ネガティブ・マッチを抜くか抜かないかというのが重要な対立軸なのね、というのが勉強になった点である。

 いちおうメモしておくけど、この論文、細かい点は良くわからないことが多くて...
 76個の指標をよく見ると、定義が全く同じ奴が散見されるんだけど、これはどういうことなんだろうか。
 実験手続きもよくわからない。該当箇所を逐語訳すると「ランダムな二値データセットをデータセットとして用いる。参照セットは30個の二値事例からなり、それぞれが100の二値特徴を持つ。この参照セットを使った一回のテスト・クエリで、距離ないし類似性の値が100個、それぞれの指標について生成される。2つの指標の間の相関係数を使ってデンドログラムを作る」「30回の独立な試行を平均してデンドログラムを作る」とのこと。どういうこと? 長さ100のランダム二値特徴ベクトルを1本作り、それとは別に100本作り、1本と100本の類似性をk種類の指標で測り、100行k列の行列からk行k列の相関行列をつくり、以上を30回繰り返して平均し、できあがったk行k列の相関行列を使ってk種類の指標を階層クラスタリングした、ってこと? というか、ランダム二値データって母比率はなんなの? ま、細けぇこたぁいいんだよ、って気もするけど。

論文:データ解析(2018-) - 読了:Choi, Cha, & Tappert (2010) 二値類似性指標はたくさんあるので、その類似性を調べてみた

Wright, M.N., Ziegler, A. (2017) ranger: A fast implementation of random forests for high dimensional data in C++ and R. Journal of Statistical Software.

 Rの高速ランダム・フォレストパッケージrangerの紹介。うっかりarXivのバージョン(2015)を読んじゃったけど、中身は同じだと思う。

 rangerについてご紹介します。分類木、回帰木、生存木に対応してます。速いです。高次元データに最適化してます。C++のスタンドアロン版もあるけど、Rパッケージ版をお勧めします(速さは同じです)。Rにおける標準的実装であるrandomForestパッケージに近い結果が得られます。云々。

 rangerは手法が新しいというより、高速・低メモリ消費を謳うパッケージ。紹介されている実験の結果をみると、

論文:データ解析(2018-) - 読了:Wright & Ziegler (2017) ランダム・フォレストのRパッケージranger

 たまにSNSとかを見るたびに、なんだか知らんがすごそうな若い人たちが、なんだか知らんがすごそうなことを語っていて(データサイエンティストっていうんですか?すごいですねえ)、彼我の違いにクラクラしてしまう。
 私の仕事のかなりの部分をデータ解析が占めるんだけど、データ解析というものは、私の知る限り実に地味なものでして、9割がたは前処理でございます。あんまし憧れるようなものではないと思うんだけど、きっと世の中にはそうでない世界もあるのだろう。

 Rには前処理用の便利機能を提供するパッケージが結構あるが、私はそういうのをあんまり使っていない。自力で書いちゃったほうが早いし、なんとなく安心できるからだ。
 その一方で、そういう機能を活用した方が効率的なのかな、と後ろめたい気分になることもある。でもああいうの、ほんとに忙しいときには、いちいち調べるのも面倒なんですよね。

 という事情により、ちょっと調べておくことにした。
 caretパッケージは機械学習系の予測モデル構築で定番のパッケージだが、前処理用の機能もいろいろ用意している。パッケージ解説の3章 Pre-Processing で紹介されている前処理用の関数についてメモする。
 (あとで気づいたのだが、逐語訳に近いものを作っておられる方がいた。ありがたいことであります)

 ほかにnearZeroVar(), findCorrelation(), spatialSign()が紹介されているが、preProcess()でも実行できるようなので省略する。

 。。。というわけなのだが、preProcess()の機能がとても多くて面食らう。
 ここからはpreProcess()について、マニュアルを参照してメモしておく。

 preProcess()は変数変換や削除そのものを行う関数ではなく。変換のためのパラメータを推定したり、どの変数を削除べきかを探してくれたりするだけ。実際の処理はこの関数が戻すオブジェクトとデータを引数にとったpredict()で行う。
 基本的な使用方法は:
 preProcess(x, method, ...)
 xは予測子がはいっている行列かデータフレーム。numericでない変数は単に無視される。
 methodは以下の文字列のうちひとつ、ないし複数個のベクトルをとる。
 ...でとれる引数(以下では便宜的に「追加引数」と呼ぶ)には3種類ある。

 さて、methodがとる文字列は以下の通り。複数個指定した場合の処理順序も決まっているので、その順にメモする。

  1. "zv": 分散がゼロの変数を削る。
  2. "nzv": 分散がゼロに近い変数を削る。caret::nearZeroVar()をコールするのであろう。追加引数としてfreqCut, uniqueCutをとる。
  3. "corr": caret::findCorrelation()をコールし、相関の高い変数を削る。追加引数としてcutoffをとる。
  4. 以下の変換。うーん、目的変数じゃなくて予測子をBox-Cox変換するというのもなんだか変な気分だが、そういうこともあるんでしょうね。
    • "BoxCox": ご存じBox-Cox変換 $(x^\lambda-1)/\lambda$。xに含まれている各変数について$\lambda$を決めてくれる。追加引数としてfudge(tolerance値), outcome(numericないしfactorのベクトル。どう使われるんだろう、よくわからない), numUnique(xのユニークな値の数がこの値より小さかったら変換しない)がある。内部でcaret::BoxCoxTrans()を呼ぶようだが、結局それはMASS::boxcox()のラッパーらしい。
    • "YeoJohnson": Yeo-Johnson変換という、Box-Cox変換みたいなものがあって、負の値もゼロもとれるらしい。試してみたら$\lambda$なるものを決めてくれた。webで拾ったPDFによれば、$x$が0以上のときは$((x+1)^\lambda-1)/\lambda$つまり$x+1$のBox-Cox変換で、$x$が負の時は$-[(-x+1)^{2-\lambda}-1]/(2-\lambda)$となる由。いったいなにがなんだか。上記PDFにも「解釈は困難」と書いてあって、笑ってしまった。
    • "expoTrans": exponential transformation. caretのマニュアルでreferされているManly(1976)の最初のページによれば$(\exp(\gamma x)-1)/\gamma$だそうである。これもBox-Cox変換みたいなもので、負の値がとれる。濱崎ら(2000)は「指数型べき変換」と呼んでいる。$\gamma$を決めてくれるのかなあ? 試してみたけどよく分からなかった。内部でcaret::expoTrans()を呼ぶらしいのだが、マニュアルに説明がない。
  5. "center": 平均を引く。
  6. "scale": SDで割る。
  7. "range": 指定した範囲に収める。範囲は追加引数rangeBounds(長さ2のnumericベクトル)で指定する。
  8. 以下の欠損値補完。
    • "knnImpute": K最近隣法による欠損値補完。追加引数としてk, knnSummary(近隣の値を平均するための関数。デフォルトではmean)がある。
    • "bag-Impute": 他の変数を全部使ったbagged treeモデルによる欠損値補完。
    • "medianImpute": 単に中央値で欠損を埋める模様。
  9. "pca": PCAで新しい変数をつくる。追加引数としてthresh(分散の累積割合のカットオフ),pcaComp(主成分の数。threshをoverrideする)がある。
  10. "ica": ICAで新しい変数をつくる。fastICA::fastICA()をコールする。追加変数はfastICA()にパスされる。
  11. "spatialSign": caret::spatialSign()をコールする。spatial sign変換というのをやるらしいのだが、意味がよくわからん、今度調べよう。

このほかに、マニュアルでは実行順序がはっきりしないが、次のmethodがある(きっと最初に実行されるんだろうな)。

 。。。ううむ。
 自分の仕事にあてはめて考えると、dummyVars()には代替がいっぱいある。その場になってみないとわからないけど、私はdummiesパッケージを使うか、model.matrix()で計画行列を作ることを考えると思う。classDist()はちょっと使い道が思いつかない。preProcess()の機能のうち、center, scale, rangeは自前で書くだろうし、変数の単変量分布の観察にかなり時間をかけることが多いので、zv, nzvは手作業でやる。Box-Cox系の変換は、目的変数ならともかく、予測子のほうに掛けたくなったという記憶が無い。欠損値補完や変数直交化をやるとしたらかなり気合いを入れるので、caretではやらないと思う。
 こうしてみると、使うことがありそうなのは、ひとつはfindLinearCombos()。こういうの、以前SASで苦労して自作したことがある。もうひとつは、意外にもpreProcess()のconditionalXだ。確かにね、こういうことが問題になるケース、ありそうですね。naive bayesのCPTにゼロを出したくない、とか。
 ... いやいや、そういう問題じゃないか。あれもこれも全部統一的な枠組みでできます、ってところがポイントなんでしょうね。業務改革への道は遠いぞ。

雑記:データ解析 - 覚え書き:caretパッケージが提供する前処理用の関数たち

« 2018年4月 | メイン | 2018年6月 »

rebuilt: 2018年6月20日 17:45
validate this page