elsur.jpn.org >

メイン > 雑記:データ解析

2017年1月22日 (日)

 ここんとこ山のようなデータ整形を抱え込んでいて、これはほんとに終わるのかという不安と、こんなことをしていてよいのかという不安に両側から押し潰されそうだが... とにかく仕事があるだけありがたいと考えるべきか。
 ともあれ、ちょっと驚いたことがあったのでメモしておく。なんだか恥をさらしているような気も致しますが。

 Rで、ユニークなキーを持つとても大きなテーブルを左外部結合(LEFT JOIN)したいとき、どうするのが速いか。
 思いつくのは、

  1. 素直にmerge()する
  2. merge()より高速だといわれている, Hadley神つくりたもうた dplyr::left_join() を使う
  3. match()で行位置を探し、テーブルの行を並び替えてくっつける
  4. match()で行位置を探しておき、テーブルの各変数を ひとつづ並び替えてくっつける

 4番目のは一種の冗談のつもりであった。そんなん遅いにきまってるじゃん? たぶん2,1,3,4の順だろう。いやひょっとして2,3,1,4かな。私はそう思いました。

 衝撃の結果は下図。クリックで拡大。Rpubsにものっけてみました

xxx.png

 4がぶっちぎりで優勝。ぐぬぬぬ。これまでの苦労は何だったんだ。

RでどうやってLEFT JOINするか (ああなんて最先端な話題だろうか)

2016年6月10日 (金)

 仕事の都合で必要になって書いたメモなんだけど、ブログに載せておこう。

 二値変数を従属変数とする回帰モデルについて考える。よく使われているのはロジスティック回帰モデル
 $\log(\frac{\pi(x)}{1-\pi(x)}) = \beta_0 + \beta_1 x$
だけど($\pi(x)$は生起確率ね)、リンク関数はほかにもある。諸君、視野を広く持ちたまえ。

 どういうときにどのモデルを使うといいのか?
 まず、分析の主な目的が確率の推定で、共変量の効果の推定はその次だ、という場合。こういうときは、ロジスティック回帰のかわりにプロビットやcloglogやloglogを使うのも悪くない。
 分析の目的は共変量の効果の推定なんだけど、ロジスティック回帰だとオッズ比で表現されちゃうのがいやだ、という場合には、logリンクか線形モデル。

 以上にあげた6つのモデルを比べる際には、とりあえずロジスティック回帰の確率推定値でケースを10群にわけ、適合を比べるのがよろしかろう。
 左右対称なロジットやプロビットを使うのがよいか、非対称なcloglogやloglogを使うのがよいかを調べる方法がある。Stukel検定という、共変量を追加すべきかどうかの検定手法の応用である。残念ながらこの手法が載っているソフトはないので、手でやること。
 まずロジスティック回帰をやって、ロジット推定値$\hat{g}(x)$と確率推定値$\hat{\pi}(x)$を得る。次に、次のふたつの人工的な共変量をつくる。
 $z_1 = 0.5 \times [\hat{g}(x)^2] \times I[\hat{\pi}(x) \geq 0.5]$
 $z_2 = -0.5 \times [\hat{g}(x)^2] \times I[\hat{\pi}(x) < 0.5]$
すべての確率推定値が0.5よりも右なり左なりだったら片方だけつくればよろしい。
 で、これを追加してもう一度ロジスティック回帰モデルを推定し、係数の信頼区間を求める。もし両方の信頼区間に0.165がはいったらプロビットがよろしい。もし$z_1$の信頼区間に-0.037, $z_2$の信頼区間に0.620がはいったらloglog、逆になったらcloglogがよろしい。[←へええええええ!]

 以上、Hosmer, Lemeshow, Sturdivant (2013) "Applied Logistic Regression", Third Edition, Section 10.3 より。

 このメモは、amazonから届いた箱をあけて本をパラパラ捲っているときにこの章に引き込まれ、なるほどー!cloglogやloglogという選択肢も頭に置いておこう!と感心して書いたんだけど、読み返してみて思うに、自分にとっては、cloglogやloglogを使う場面ってやはりそんなにはなさそうだ。もし係数の解釈に関心がなく単に予測したいだけで、かつロジスティック回帰の予測精度が悪かったら、きっと機械学習系の手法を試すと思う。

 いっぽう、ロジスティック回帰はもうイヤ!説明変数の効果はリスク差で示したいの!ああもう加法モデルにしちゃいたい!誰かボクをここから連れ出して!海辺の町に連れてって!と切実に思うことなら、それは頻繁にある。そういう哀れな分析者のためのガイドラインが欲しいんだけどなあ。加法モデルを使う際の注意点でもいいし、誤魔化しながらロジスティック回帰を使い続けるコツでもいいから。

ロジスティック回帰に飽きたときのための個人的な覚え書き

2016年4月20日 (水)

 新製品・サービスの普及プロセスを説明する古典的モデルである、Bassモデルの推定方法について、ちょっと混乱しちゃったので、覚え書きをつくっておく。個人的メモでごさいます。

 そもそもBassさんが考えた理屈はこうだ。原典であるBass(1969)に沿ってメモする。
 ある製品の初回購入だけに注目する。総購入個数を$m$としよう。
 連続時間上の時点 $t$ における「購入の尤度」を $f(t)$ とする[わかりにくい表現だけど、ある潜在的購入者に注目したときにこの瞬間に購入が発生する確率密度、という意味合いであろう]。
 $F(t)$ を以下のように定義する。
 $F(t) = \int f(t) dt, \ \ F(0) = 0$
 時点 $t$ における売上を $S(t)$、累積売上を $Y(t)$ とする。

 時点$t$において、まだ購入が発生していないという条件の下での「購入の尤度」を$P(t)$とすると、
 $P(t)=f(t)/(1-F(t))$
である。また、
 $S(t) = m f(t) = P(t)(m-Y(t))$
 $Y(t) = \int S(t) dt = m \int f(t) dt = m F(t)$
である。

 さて、Bassモデルとは次のモデルである。
 $P(t) = p + (q/m) Y(t)$
 ここから以下が導かれる。
 $f(t) = p + (q-p) F(t) - q F(t)^2$
 $S(t) = pm + (q-p) Y(t) - (q/m) Y(t)^2$

 さらに、$f(t), F(t)$を$t$について解くと、
 $f(t) = \frac{(p+q)^2}{p} \frac{\exp(-(p+q)t)}{(q/p \exp(-(p+q)t)+1)^2}$
 $F(t) = \frac{1-\exp(-(p+q)t)}{q/p \exp(-(p+q)t)+1}$

 ここまではいいっすね。ここからはBass's Basement Research Instituteによる説明のメモ。

 実際には$t$は離散変数だ。すると、次のような困った事態が発生する。
 区間$[t-1, t]$について考えよう。$F(t)$は、この区間の右端までの$f(t)$の合計? では$f(t)$は? 時点$t$における瞬間の速度? それとも、$F(t-1)$から$F(t)$までの増分?
 アルゴリズムとして実装するには、どうにかつじつまを合わせないといけない。ここで3つのアプローチが登場する。

 第一のアプローチ。Sinivasan & Mason (1986)のアプローチである。$F(t)$のほうはさきほどの解をそのまま用いる。すなわち
 $F(t) = \frac{1-\exp(-(p+q)t)}{q/p \exp(-(p+q)t)+1}$
いっぽう$f(t)$のほうを変えてつじつまを合わせる。
 $f(t) = F(t) - F(t-1)$
 ただし、$t=1$においては$f(t) = F(t)$。
 このやり方の良い点:$F(t)$の定義はそのまま。$f(t)$は小さな値なので、これを直接求めずに引き算で求めるのは、数値的に正確。

 第二のアプローチ。$f(t)$のほうの解をそのまま用いる。すなわち
 $f(t) = \frac{(p+q)^2}{p} \frac{\exp(-(p+q)t)}{(q/p \exp(-(p+q)t)+1)^2}$
いっぽう$F(t)$のほうを変えてつじつまを合わせる。
 $F(t) = \sum_{i=0}^{t} f(i)$
これはおかしい。$f(t)$はその瞬間の速度であって、ある期間の平均ではないからだ。

 第三のアプローチ。Bass先生が1969年の段階で推していた方法。$f(t), F(t)$の差分方程式に戻る。
 $f(t) = p + (q-p) F(t-1) - q F(t-1)^2$
 $F(t) = F(t-1) + f(t)$
これはおかしい、すごくおかしい。Bass先生はなぜこんな方法を推したのか? コンピュータじゃなくて手回し計算機を使っていたからに他ならない。なんとかしてOLS回帰の形に持ち込む必要があったのだ。この$f(t)$の式なら、$f(t)$を当期購入数、$F(t-1)$を前期までの累積購入数として、重回帰で解けるからね。
 教訓。Bass(1969)の推定方法を用いるなかれ。コンピュータを使いなさい。第一のアプローチを使って、非線形回帰で推定しなさい。

 ... なるほどね。
 最後に、手元にあるいくつかの教科書を比較しておこう。

 日頃より頼りにしております、朝野・山中(2000), p.152。
 $m$をN, $F(t)$をn_T, $p$をp, $q$をr, $Y(t-1)$をZ_t (添え字に注意!), $S(t)$をX_tと呼んでいる。
 本文中でZ_tを「t期における購入者」、X_tを「t期の購入者」と呼んでいるので、前後の文脈を読み取ってその意味をよく考えないといけない。
 推定方法はOLS。推定する式はこうなる(8.5式):
 $X_t= pN + (r-p) Z_t - \frac{r}{N} Z_t^2$

 照井・ダハナ・伴(2009), pp.155-157。この本も頼りになります。
 $m$をN, $F(t)$をF(t), $p$をp, $q$をq, $Y(t)$をN(t), $S(t)$をS(t), と呼んでいる。
 推定方法はOLS。本文中の記号とRのコード例での変数名が違うので要注意。コード例ではなぜか$N(t-1)$にFt.1, $n(t)$にNtという名前をつけ、lm(Nt ~ Ft.1 + Ft.1^2)を推定している。(どうもよくわからないんだけど、formulaの指定はこれでいいのだろうか? Nt ~ Ft.1 + I(Ft.1^2) ならわかるんだけど...)

 里村(2010)、pp.64-67。
 $m$をm, $F(t)$をF(t)、$f(t)$をf(t)、$Y(t)$をN(t)、$S(t)$をX(t)と呼んでいる。$m-Y(t)$をY(t)と呼んでいるので注意。
 この本は、$F(t)$を$t$の関数として解いた式を非線形回帰で推定している。推定にはRのnls()を使っている。シリーズの主旨に沿って、Rのコード例がとても充実している。よく見ると日本語の誤字が結構多いのだが(「統計的に有為」とか)、そういうことをいちいち気にしていてはいかんのである。

 Lilien & Rangaswamy (2004), pp.255-258。ちょっと古い版の教科書だが、もったいないので使い続けている。前職のときに私費で買ったのだ、高かったのだ。
 $m$を$\bar{N}$、$F(t)$をF(t), $f(t)$をf(t)、$p$をp, $q$をq, $Y(t)$をN(t), $S(t)$をs(t), と呼んでいる。推定方法として、OLSと非線形回帰の両方を紹介している。
 分量も違うし、そもそも日本語の教科書と英語の教科書を比較してはいけないのだが(なにしろ市場規模が違う)、とにかくこの本の説明は最高にわかりやすい。もちろんBass(1969)なんかよりはるかにわかりやすい。

メモ:Bassモデルをどうやって推定するか (いろんな教科書を比べてみた)

2016年1月14日 (木)

 たとえば、プロのサッカー選手が来ているシャツにはスポンサーのロゴなんかが入っている。小さなロゴの裏では莫大な金が動いている。
 スポンサー契約を通じて、スポンサーは知名度を向上させたりなんだり、なんらかの価値を得るだろう。ひょっとするとクラブの側も、契約金とは別になんらかの価値を得るかもしれない(ないし、変な企業とスポンサー契約したせいで損をするかもしれない)。
 スポンサー契約はどれだけの価値を生むのか。

 あるスポンサー契約がクラブとスポンサーにもたらす、(契約金以外の)価値の合計について考える。その価値を手に入れたのがクラブ側かスポンサー側かはいったん脇に置いておく。この価値の合計のことを、経済学者っぽく「生産価値」と呼ぶことにする。
 スポンサー契約の生産価値はなにによって決まるか。クラブにもいろいろあって、スポンサー契約が価値をもたらすようなクラブもあれば、そうでもないクラブがあるだろう。スポンサー企業にもいろいろあって、スポンサー契約が価値をもたらすような企業もあれば、そうでもない企業もあるだろう。さらに、クラブとスポンサーの相性というものもあるだろう。
 次のように考えよう。市場$t$におけるクラブ$a$の特徴をベクトル$X_{at}$、スポンサー$i$の特徴をベクトル$Y_{it}$で表す。生産価値は
 $f(a, i, t) = \alpha X_{at} + \beta [X_{at} Y_{it}] + \gamma Y_{it} + \epsilon_{ait}$
 第1項は、いうなればクラブの力。第2項は相性の力で、ブラケットの中身をどうするかはあとで考える。第3項はスポンサーの力。最後は誤差だ。

 ここで知りたいのは、係数$\alpha, \beta, \gamma$だ。これらが推定できれば、スポンサー契約の価値を推定する仕組みが手に入る。クラブとスポンサーの特徴をインプットすれば、スポンサー契約がもたらす価値の推定値がアウトプットされる仕組みだ。素晴らしい。
 そこで、スポンサー契約の事例を片っ端から集めてくる。さらに、クラブとスポンサーの特徴についてのデータを片っ端から集めてくる。クラブとスポンサーの相性に影響しそうなデータも集める。
 たとえば、大チームと大企業を組み合わせると相乗効果が生まれたりするかもしれない。 クラブとスポンサーが地理的に近いところにある方が相性は良いかもしれない。強くて人気がある大チーム、大企業、国際的企業、サッカー向きの業種の企業は力を持ち、スポンサー契約の価値を高めるかもしれない。でもその力も距離次第では損なわれてしまうかもしれない。なんであれ、スポンサー契約を続けていると価値が増すかもしれない。。。という具合に、思いつく仮説を、上の数式にどんどん入れていく。

 さあ、準備はできた。データを上の式に放り込み、係数$\alpha, \beta, \gamma$を推定しよう... と思うところですよね。私もそう思いました。しかし、話はそのようには進まない。これがこのメモを書き始めた理由である。
 なぜか。スポンサー契約がもたらした価値についてのデータがないからだ。クラブの観客数や企業の売上からスポンサー契約の価値を割り出すのは難しい。契約金からどうにか推定できるとしても、そもそもスポンサー契約の金額は企業秘密だ。式の右辺についてはデータがある、しかし左辺についてのデータがないのである。
 さあ、ここからが本題。スポンサー契約がもたらした価値についてのデータなしで、スポンサー契約の価値を推定するモデルをどうやってつくるか。

 サッカークラブは金をくれるスポンサーを求め、スポンサーはロゴをつけてくれるサッカークラブを求めている。クラブたちとスポンサーたちはひとつの市場を形成している。たとえばあるシーズンにおけるある国のクラブたちとスポンサーたち、これがひとつの市場だ。
 いま、市場 $t$ においてクラブ $a$ がスポンサー $i$ と契約したとしよう。スポンサーがクラブに渡す金額を $r_{ait}$、スポンサーが得る価値を$\Delta V (a,i,t)$, クラブが得る価値を$\Delta U(a,i,t)$としよう。スポンサーの利得は
 $\pi^S (a,i,t) = \Delta V(a,i,t) - r_{ait}$
クラブの利得は
 $\pi^C (a,i,t) = \Delta U(a,i,t) + r_{ait}$
この契約によって生まれる価値、すなわち生産価値は
 $f(a,i,t) = \Delta V(a,i,t) + \Delta U(a,i,t)$
である。両者の間でどれだけのカネが動いたか($r_{ait}$)は、もはやどうでもよくなっていることに注意。

 同じ市場$t$において、別のクラブ $b$ が別のスポンサー $j$ と契約したとしよう。このとき、世の中がうまく回っているならば、
 $f(a,i,t) + f(b,j,t) \geq f(a,j,t) + f(b,i,t)$
であるはずである。つまり、2つの契約から生まれる生産価値の合計は、仮にクラブとスポンサーのマッチングを入れ替えたときに生まれる生産価値の合計と同じ、ないしそれよりもマシであるはずだ。マッチング理論ではこれを「局所生産最大化条件」と呼ぶのだそうである。

 市場$t$においてなされたスポンサー契約の集合から、2ペア$\{a,b,i,j\}$を取り出すすべての取り出し方について考える。契約が3つあったら3x2=6通り、k個あったらk(k-1)通りあるんでしょうね。で、すべての取り出し方について、局所生産最大化条件が満たされていたら1点、そうでなかったら0点とカウントする。このカウントを合計しよう。
 さらに、市場の数が$H$個あるとして、それらの市場を通じて、カウント合計の平均を求めよう。
$Q_H(f) = \frac{1}{H} \sum_{t \in H} \sum_{\{a,b,i,j\} \in A_t} 1[f(a,i,t) + f(b,j,t) \geq f(a,j,t) + f(b,i,t)]$
 $1[\cdot]$は、カッコ内の不等式が成立しているときに1, そうでないときに0を返す関数である。上記数式、原文にはどうもミスプリがありそうなので、勝手に表記を変えている。誤解していないといいんだけど。

 この関数を$f$について最大化した解を、スポンサー契約市場における均衡状態と捉えることができるのだそうだ。
 つまり、もし世の中がうまく回って回って回り続けていれば、いずれはスポンサー契約市場 がそうなるであろう姿。一旦世の中がそうなってしまった暁には、どの(合理的な)クラブもスポンサーもそこから抜け出すことができない、そんな姿。それがわかるというわけだ。まじですか。
 ということは、もし世の中がこれまでうまく回って回って回り続けているならば、スポンサー契約市場は均衡状態に陥っているはずだ。ということは、スポンサー契約から生まれる価値の関数 $f$ は、上の$Q_H(f)$を最大化するような関数であるはずだ。ということは、実際のスポンサー契約、ならびにクラブとスポンサーの特徴についてのデータの下で、$Q_H(f)$を最大化する$f$を求めれば、それがこの世の中における、スポンサー契約の価値を求める関数となるはずだ。
 ... というロジックなのではないかと思う。えーと、こういう理解であっているんでしょうか。こういう考え方そのものに、私、 いささか戸惑ってしまうんですが...

 このスコア関数$Q_H(f)$にさきほどの$f(a, i, t)$の式を放り込もう。ここで面白いのは、$Q_H(f)$の中身の不等式をよく見ると、$\alpha X_{at}, \alpha X_{bt}, \gamma Y_{it}, \gamma Y_{jt}$がひとつづつ出てくる、という点だ。つまりこれらは無視してしまってよい。問題はクラブとスポンサーの相性の力 $\beta [X_{at} Y_{it}]$ だけなのだ。おおお、なるほど。
 ともあれこのようにして、スポンサー契約、クラブ、スポンサーについての十分なデータがあれば、スポンサー契約がもたらす価値を推定するモデルを手に入れることができる。
 イギリスのサッカークラブのスポンサー契約のデータをつかって推定したところ、クラブとスポンサー企業は規模が釣り合っているときに相性がいいとか、クラブ本拠地と企業の本社所在地が地図上で近いほうが相性がいいとか、そういったことがわかった由。

 以上、Yang & Goldfarb (2015, J. Marketing Research) からメモ。
 実は上記は、私には途方もなく難しいこの論文のごく最初のほうに出てくる話で、ここから論文は「酒やギャンブルに関連する企業がサッカーのスポンサーになるのを禁止したら何が起きるか」という分析へと進んでいくのだが、分析の仕組みというか建付けの部分で途方に暮れてしまった。頭を冷やすために、論文のロジックを組み立て直し、私にも理解できるくらいに平易な筋立てに落として、メモにしてみた次第である。

 均衡という概念に基づいたこういう分析が、ビジネス・ リサーチやデータ解析でもこれから重要になってくるのか、そうでもないのか、そういう大きな話は私にはよくわからない。だから、手探りで勉強することにどれだけの投資対効果があるのかはよくわからないんだけど、 とにかくその、読むたびに途方に暮れてしまうのである。なんというかその... 現象を理解するための分析に、いつのまにか規範的言明が入ってくる感じ、というか...
 もっとも、ふつうの統計的分析だって、インプットは常に仮定とデータだ。上の場合でいうと、仮にスポンサー契約の価値のデータが手に入っていたら回帰分析でモデルを推定できただろうけど、でもそのときだって、きっと誤差項の分布についてなんらかの仮定を置くだろう。上の分析ではその確率的な仮定の代わりに「スポンサー契約市場が均衡状態にある」というゲーム理論的な仮定を置いただけだ、ということなのかもしれない。ううううむ。。。

メモ:スポーツ・スポンサー契約がもたらす価値をどうやって推定するか:マッチング理論の巻

2015年12月17日 (木)

 なんというかその、判別モデルや回帰モデルを組んだ時、「で、どの説明変数が重要でしたか?」って聞かれること、あるじゃないですか。いやいや、ひとことで「重要」っていってもいろいろありましてですね... と交互作用や非線形性の話をしても、ある意味、煙に巻くようなもんじゃないですか。素朴な気持ちとしては、こっちも内心では 「いやあこの変数は効くなあ」なんて思っているわけだし。なんだかんだいいながら、どうにかして変数を選ばないといけないわけだし。
 魚心あれば水心で、世の中には実にいろんな変数重要度指標があり、誠に困ったことである。さらに困ったことに、こういう話は仕事の中で不意に顔を出すので、いつかきちんと調べます、では間に合わない。
 とりあえず、いろんなタイプの予測モデルを統合したパッケージである、R の caret パッケージが実装している重要性指標についてメモしておこう。出典はこちら。純粋に自分用の覚え書です。

 まずは、モデルに依存しないタイプの指標。

モデルに依存するタイプの指標。

 ふうん...
 ランダムフォレストのような協調学習の分野では、変数重要性をpermutationベースで求めるとき、ふつうにその変数だけかきまぜちゃっていいのか (いわばmarginalな重要度になる)、それとも共変量で層別して層の中だけでかきまぜるべきか (conditionalな重要度になる)、という議論があるらしい。これは他の手法でもいえることであって、その辺のところが知りたいんだけど、caretパッケージではさすがにそこまで面倒はみてくれないようだ。

説明変数の重要度指標のいろいろ (by caretの中の人)

2015年6月16日 (火)

 先日から仕事の都合で生存時間(イベント発生時間)データの分析について考えていたのだけど、資料をみていると、サンプリング・デザインの用語がいろいろあって混乱してしまう。同じことを違う名前で呼んでいることもあって、ちょっと困る。
 というわけで、クライン&メシュベルガー「生存時間解析」より関連用語をメモ。この本、原題を"Survival Analysis: Techniques for Censored and Truncated Data"というだけあって、サンプリング・デザインの問題についてやたらに詳しい。

1. 打ち切り(censoring)。

2. 切断(truncation)。ある期間中にイベントが発生した人だけが観察される。

やれやれ...

生存時間データの打ち切り・切断に関する用語総ざらえ

2015年6月11日 (木)

調査におけるウェイティング(いわゆるウェイトバック集計)について、このブログがきっかけでお問い合わせを受け、えええ? いったいこのブログにどんなことを書いてんだっけ... と読み返してみたところ、全然忘れているエントリもあってびっくりした。三歩歩けば全て忘れる、ニワトリなみの記憶力である。
 せっかくなので、これまでに書いたウェイティングに関するエントリを時間順に並べてみた。すいません、完全に自分用のメモです。

ウェイティング回顧録 (私は結構ヒマな男なのではないか)

2014年10月17日 (金)

 以下は純粋に個人的なメモであります。

 ここにn変量時系列 $y_t$ がございます。
 これを以下のように分解する。まずは
 $y_t = B + Lf_t + u_t, \ \ \ u_t \sim N(0, \Sigma_u)$
$B$は切片(サイズ n x 1), $L$は係数行列(n x p), $f_t$はp変量の潜在動的因子(p x 1), $u_t$はn変量の不規則要素(n x 1)である。
 潜在動的因子をトレンドと季節に分解する。季節要素を観測値のモデルにではなく因子の時系列モデルのほうに入れている点に注意。
 $f_t = \alpha_t + \gamma_t$
トレンドのほうは
 $\alpha_t = \alpha_{t-1} + \beta_{t-1} + \epsilon_t, \ \ \ \epsilon_t \sim N(0, \Sigma_\epsilon)$
トレンドの傾きも時間変動させて
 $\beta_t = \beta_{t-1} + \delta_{t-1} + \eta_t, \ \ \ \eta_t \sim N(0, \Sigma_\eta)$
傾きの変化率をランダムウォークさせる。
 $\delta_t = \delta_{t-1} + \zeta_t, \ \ \ \zeta_t \sim N(0, \Sigma_\zeta)$
季節のほうも確率的に変動させる。以下、添え字が煩雑になって混乱するので、季節周期を4に決め打ちする。
 $\gamma_t = -\sum_{j=1}^{3} \gamma_{t-j} + \xi_t, \ \ \ \xi_t \sim N(0, \Sigma_\xi)$
 以上、すごく複雑に見えるけど、因子の時系列構造を工夫しているだけで、要するに因子分析である。びびってはいけない。

 さて、状態空間表現ではどうなるんだよ、という話である。さあ深呼吸!
 観察方程式は
 $y_t = A z_t + B + u_t, \ \ \ u_t \sim N(0, \Sigma_u)$
ここで状態変数ベクトルは
 $z_t \equiv [\alpha_t \ \beta_t \ \delta_t \ \gamma_t \ \gamma_{t-1} \ \gamma_{t-2}]'$
実際には$\alpha_t, \beta_t, \ldots$は縦ベクトル(p x 1)なので、こういう書き方でいいのかどうかよくわかんないんだけど(自慢じゃないけど高校で数学の時間はずっと寝ていた)、原文でもこうなってるし、まあいいことにしよう。なお、原文で最後の項は$\gamma_{t-s-2}$となっているが、$\gamma_{t-(s-2)}$の誤りであろう。
 係数行列は
 $A \equiv [L_{n \times p} \ 0_{n \times p} \ 0_{n \times p} \ L_{n \times p} \ 0_{n \times p} \ 0_{n \times p}]$
 順に、$\alpha_t$の係数, $\beta_t$の係数、$\delta_t$の係数, 季節ダミー1の係数($\alpha_t$の係数と同じ。なぜなら季節要素が状態方程式に突っ込んであるから), 季節ダミー2の係数、季節ダミー3の係数、である。結局、$z_t$の長さは$m = 6p$。季節周期を$s$として、$m=(2+s)p$である。

 さあ状態方程式だ!
 $z_t = C z_{t-1} + v_t, \ \ \ v_t \sim N(0, \Sigma_v)$
 状態撹乱要素は
 $v_t \equiv [\epsilon_t \ \eta_t \ \zeta_t \ \xi_t \ 0 \ 0]'$
 問題は遷移行列$C$。原文では次のようになっている。
 $C \equiv \left( \begin{array}{llllll} I_p & I_p & 0 & 0 & 0 & 0 \\ 0 & I_p & I_p & 0 & 0 & 0 \\ 0 & 0 & I_p & 0 & 0 & 0 \\ 0 & 0 & 0 & -I_p & \cdots & -I_p \\ 0 & 0 & 0 & I_p & 0 & 0 \\ 0 & 0 & 0 & 0 & I_{(s-3)p} & 0 \end{array} \right)$
 本当かなあ? ここにどうも納得がいかなかったのである。
 いちから考えてみたい。$C$の右側にあるのは、たとえば$[\alpha_4 \ \beta_4 \ \delta_4 \ \gamma_4 \ \gamma_3 \ \gamma_2]'$である。これに$C$を掛け、つくりたいのは$[\alpha_5 \ \beta_5 \ \delta_5 \ \gamma_5 \ \gamma_4 \ \gamma_3]'$である。以下、大きな行列を書くのが面倒なので、$C$をp行p列づつ区切り、上の段から順に考える。
 最初のp行は$\alpha_5$をつくる奴だから、
 $(I \ I_ \ 0 \ 0 \ 0 \ 0)$
 次のp行は$\beta_5$をつくる奴だから
 $(0 \ I \ I \ 0 \ 0 \ 0)$
 次のp行は$\delta_5$をつくる奴だから
 $(0 \ 0 \ I \ 0 \ 0 \ 0)$
 次のp行は$\gamma_5$をつくる奴だから
 $(0 \ 0 \ 0 \ -I \ -I \ -I)$
 次のp行は$\gamma_4$をつくる奴だから、
 $(0 \ 0 \ 0 \ I \ 0 \ 0)$
 最後のp行は$\gamma_3$をつくる奴だから、
 $(0 \ 0 \ 0 \ 0 \ I \ 0)$
 ほんとだ...! 合っている!
 やれやれ。分散行列 $\Sigma_v$はブロック単位行列で、左上から順に$\Sigma_\epsilon, \Sigma_\eta, \Sigma_\zeta, \Sigma_\xi, 0_p, 0_p$である。

 なお、実際にはすべての分散行列を対角行列に制約した由。$\Sigma_v$のみならず、観察方程式の分散行列$\Sigma_u$も対角行列にしちゃうのである。不規則要素に共分散はないと考えたわけだ。ま、いいか。
 パラメータ推定はEMアルゴリズムで行った由。詳細は略するが、ステップのなかでバリマクス回転している。いやーん。そんなんようせんわ。
 なお、データはあらかじめ対数変換した由(歪度が大きかったから)。さらに平均0, 分散1に標準化した由(分析の主旨は動的因子にあるので、別にかまわない)。

 以上、Du & Kamakura (2012, JMR)のWeb Appendix Bより。とても歯が立たないだろうと敬遠していたが、実際に読んでみたら、もちろん難しいんだけど、思ったほどではなかった。
 本文を読んだときも思ったことだが、季節要素をこうして動的因子系列に入れるべきかどうかはデータの実質的な性質によるなあ、と思った。著者らは自動車ブランド名の検索量の時系列を分析していて、因子として「欧州車」というようなものを想定しているので、それなりに筋が通っている(欧州車に共通の季節要素があると考えているわけだ)。でもこれをマルチカテゴリに拡張して、著者らが主張するようにtrendspottingという観点から捉えたときには、観察方程式に入れたほうが自然ではないかと思う。
 推定のところにオリジナルな工夫があるところがちょっと嫌だ。このモデルって、あらかじめ予備的分析で因子数と負荷行列にあたりをつけ(変量ごとに季節を除去しEFAにかけるとか)、$L$を何箇所か適当に固定すれば、最尤推定できたりしないのかなあ。

Du & Kamakura (2012) メモ

2014年10月 2日 (木)

 仕事の都合で状態空間モデルのソフトを使っているのだけれど、参考書によってもソフトによっても、記法がバラバラで大変に混乱する。メモをつくったので載せておく。純粋に自分のための覚え書きであります。
 話を簡単にするために、単変量のガウシアン状態空間モデルについてのみ考える。状態変数の数をmとする。

Commandeur & Koopman(2007): とりあえず参考書の例として。Section 8.1に基づく。

Koopman, Shepard & Doornik (1998): SsfPackの基になっている論文なのだが、記法がすごくややこしい...

SsfPack: Commandeur&Koopman(2007)のSection 11.2 より。

Rのdlmパッケージ: vignetteによれば、

RのKFAS パッケージ: reference manual によれば、

よく使うdlmとKFASを比較すると、

と呼ばれているわけである。やれやれ。

追記: はじめてMathJaxを使ってみた。うわあ、これは便利だ...

統計学者のみなさん、記法は統一してもらえないかしら (状態空間モデルの巻)

2014年7月 3日 (木)

多変量時系列データの分析手法を指す、「動的因子分析」とか「ダイナミック・ファクター・モデル」とかいう言葉があるけど、この言葉で指しているモデルの形式が人によってバラバラであることに気が付いた。ちょっと混乱しちゃったので、メモしておく。どうも恥をさらしているような気がしないでもないんだけど...

 なぜ混乱しちゃったかというと、次の2つのタイプのモデルが、両方とも動的因子分析と呼ばれているからである。

 Bを「動的因子分析」と呼ぶのは何となく違和感があるんだけど、そう呼んでいる例も少なくないのである。以下、順不同でメモしておくと...

ほか、気になるけど、どっちだかわかんないやつ: Du Toit & Browne (2001, in Cudeck et al.(eds.)), Wood & Brown (1994, Psych.Bull.), Browne & Zhang (2007, in Cudeck & MacCallum(eds.)), Molennar, De goooijer & Schmitz(1992, Psychometrika)

動的因子分析ってなんですか

2014年6月22日 (日)

先日リリースされたMplus 7.2の改訂点を、こちらで公開されているPDFからメモしておく。すいません、純粋に自分のためのメモです。

 一番大きな追加機能は、混合モデルで非正規分布を指定できるようになったこと。ANALYSISコマンドでTYPE=MIXTUREのときDISTRIBUTIONオプションで指定する。NORMAL(正規分布), SKEWNORMAL(歪正規分布), TDISTRIBUTION(t分布), SKEWT(歪 t 分布)から選べる。数値積分時は不可。うーむ、まだ使い方が良くわからない... Web Note #19を読めばいいのだろうか。
 なお、ESTIMATOR=BAYES;POINT=MODEのときにモード推定のための反復回数の上限をDISTRIBUTIONオプションで与えることができるけど、それと同名だがちがうオプションである。あれれ、じゃ「TYPE=MIXTURE; ESTIMATOR=BAYES;POINT=MODE;」としたらDISTRIBUTIONオプションはどちらの意味になるのだろうか。今度試してみよう。
 実はこの非正規分布の指定、TYPE=GENERALのモデルでも指定できるのだが、そちらは実験版である由。
 この新機能と関係していると思うのだが、OUTPUTコマンドにH1MODELオプション、ANALYSISコマンドにH1STARTSオプションが追加された。まだよく理解できていない。

 MODEL INDIRECTコマンドの変更。因果推論の観点から定義された直接効果・間接効果が出力されるようになった。こないだ読んだ奴に書いてあった話だ。
 INDコマンドで、一番右の変数が連続変数であるとき、その後ろにかっこをつけて比較する2つの水準を指定できるようになった。
 媒介変数(メディエータ)があるときの特別な間接効果を出力するMODオプションが追加された。メディエータ、メディエータの水準指定、外生モデレータ、外生モデレータのとる値の範囲、原因変数とメディエータの交互作用変数、外生モデレータとメディエータの交互作用変数、原因変数がとる値の2つの水準の指定、をカバーしている。なにもそこまでしなくても...

 潜在クラスモデル・潜在遷移モデルで、二値ないしカテゴリカルな指標の残差共分散を指定できるようになった。ANALYSISコマンドのPARAMETERIZATIONオプションでRESCOVARIANCESを指定する。たとえばMODELコマンドで「%OVERALL% u1 WITH u3;」とだけ指定すると、u1とu3の残差共分散がクラス間等値制約のもとで推定される。さらにクラス別に追記して、あるクラスでだけ残差共分散を自由推定したりできる。そうそう、これ、これまでできなかったんだよな...

 連続時間生存モデルの見直し。
 まず、VARIABLEコマンドのSURVIVALオプションが変わった。生存時間変数を t として、「SURVIVAL = t;」とすると、従来は定数ハザードが推定されたのだが、改訂後はCox回帰みたいなノンパラメトリック・ベースライン・ハザード関数が推定される。ただし、もし t から連続潜在変数にパスが延びていたり、マルチレベルモデルだったり、モンテカルロ数値積分をするモデルだったりすると、セミパラメトリック・ベースライン・ハザード関数となる。このとき、10個の時間間隔が勝手に決定されるが、「SURVIVAL = t(10);」とか「SURVIVAL = t(4*5 1*10);」という風に時間間隔を明示的に指定することもできるし、「SURVIVAL = t (ALL);」とするとデータ中の時間間隔がすべて用いられる。「SURVIVAL = t (CONSTANT);」とすると定数ハザードになる。
 これに伴い、ANALYSISコマンドのBASEHAZARDオプションも変わった。従来は、ONでパラメトリックなベースライン・ハザード、OFFでノンパラメトリックなベースライン・ハザードで、OFFのときのみ、その後ろの(EQUAL)ないし(UNEQUAL)でクラス間等値制約の有無を指定できた。改訂後は、パラかノンパラかではなくて、モデル・パラメータかどうかを指定するためのオプションとなった。デフォルトはOFFで、このときベースライン・ハザード関数のパラメータはモデルのパラメータではなく補助的なパラメータになる。ONにするとモデル・パラメータになり、MODELコマンドで t#1, t#2, ..., [t]が指定できる。またONでもOFFでも、後ろに(EQUAL)ないし(UNEQUAL)を指定できる。
 とこのように、同名のオプションの使用方法が変わっちゃったのでややこしい。たとえば単純なCox回帰は、これまでVARIABLEコマンドで「SURVIVAL = t (ALL);TIMECENSORED = tc (0 = NOT 1 = RIGHT);」、ANALYSISコマンドで「BASEHAZARD = OFF;」であったが、改訂後はSURVIVALオプションの(ALL)は不要となり(事情がない限りノンパラになるから)、BASEHAZARDオプションは不要となる。

 DEFINEコマンドが微妙に変わったらしい。 ちゃんと読んでないんだけど、たとえば、CENTERオプションのあとで交互作用変数をつくると、これまでは中心化される前の変数が使われたけど、改訂後は中心化された後の変数が使われるようになった、とかなんとか。

 DEFINEコマンドなどで使うDOオプションがネストできるようになった。ふーん。
 以上!

追記: いやいや、まだ細かい続きがあったのを見落としていた。

Mplus 7.2の改訂点

2014年3月30日 (日)

 このブログの趣旨に反して、これから書くのは私の覚え書きではなく、他の方に読んでもらうための文章である。どうか私が事柄を正しく説明していますように。そしてこのエントリが、検索結果の上の方に出てきますように。

 市場調査に関連する仕事で糊口を凌いでいるのだけれど、折々に受ける業務上のご質問のなかでかなりの領域を占めているのが、いわゆる「ウェイトバック集計」、すなわち確率ウェイティングに関するあれこれである。
 もっとも基本的な質問はこういうのだ。「ウェイトバックしたほうがいい?しないほうがいい?」 これはですね、ほんとに答えに困ります。

 いま手元に標本抽出確率が不均一な標本があり、集計・分析の際にウェイティングすべきかどうかが問われたとき、考慮すべきポイントは、主に4つあると思う。

 (1) については誤解がないと思う。(2)についてはいろいろ誤解があり、以前義憤に駆られて思うところを書いた。(4)は話がすごく長くなる。
 問題は(3)である。実のところこの論点こそが、調査実務に関わる人にとっての道しるべ、ウェイティングの是非という難しい問題に光を差してくれる灯火となると思うのだけれど、不思議なことに、この点についてきちんと説明した日本語の資料をほとんどみたことがない。同業種の方々と話していても、誠に失礼ながら、この点についてきちんと理解しておられる方は非常に少ないように思う。あまりに少ないので、ひょっとして私が全てを勘違いしているのではないかと不安になるほどだ。
 (3)について誰かに説明しなければならない羽目に陥ることもある。これが結構大変なのである。できれば避けて通りたい。いまはなんでも検索する時代だから、いちど文章にしてブログかなにかに書いておけば、そのぶん説明の回数も減らせるかもしれない。もし私が全てを勘違いしているのならば、それを指摘してくださる奇特な方が現れないとも限らない。
 以上が、この文章を書く主たる理由である。実はもうひとつ、ついうっかりwebアプリをつくってしまったというささやかな理由もあるのだけれど、それはあとで。

問題
 以下では、いわゆる「ウェイト・バック集計」、すなわちデータの集計・解析における確率ウェイティングについて考える。上記の4つのポイントのうち(1)(2)(4)はすでにクリアされていることにする。
 調査データの解析で確率ウェイティングが登場する局面は意外なほど多岐にわたるのだけれど(Kish(1990)は7種類挙げている)、説明に際しては、もっともわかりやすくて一般的な、非比例の層別抽出を想定する。集計対象は二値変数で、母集団全体における割合に関心があることにする。母集団は十分に大きいことにする。

 あまりにシンプルすぎて非現実的ではあるが、次の場面について考えよう。ある新製品コンセプトに対する消費者の購入意向を調べたい。そこで潜在的顧客から抽出した調査対象者に、その製品を買ってみたいどうかを「はい/いいえ」で聴取した。以下、「はい」と答える人の割合のことを購入意向率と呼ぶことにする。
 性別で層別した標本抽出を行った。つまり、男と女のそれぞれについて、決まった数の調査対象者を得た。潜在的顧客における男女比は5:5であることがわかっている。しかしなにかの事情があって、調査対象者は男60人、女40人とした(男女比は6:4)。これらの対象者は、その性別のなかでは無作為に抽出されたものとみなすことにする。

 いま知りたいのは、潜在的顧客(母集団)の全体における購入意向率である。その推測の手段として、調査対象者(標本)から得た回答を集計し、購入意向率を算出する。
 ここでunweightedの集計値とは、「はい」と答えた人の人数を100で割った値である。
 いっぽうweightedの集計値は次のように説明できる。男60人における購入意向率、女40人における購入意向率をそれぞれ求め、得られた2つの値を、母集団における割合(ここではそれぞれ0.5)で重みづけて合計する。これがweightedの集計値である。
 あるいは次のようにいいかえてもよい。まず各対象者に適切な「ウェイト」を与える。たとえば、男性の対象者には5/6=0.833, 女性の対象者には5/4=1.250というウェイトを与える。次に、「はい」を1点、「いいえ」を0点とし、各対象者の持ち点にウェイトを掛けて合計し、最後にウェイトの合計で割る。これがweightedの集計値である。

 さて、unweightedの集計値とweightedの集計値、どちらを使うべきか。ウェイティングすべきか、すべきでないか。

 たいていの人はこう答える。ウェイティングすべきです。なぜなら、母集団における男女比が標本における男女比と異なっているので、標本から得たunweightedの集計値は、男性の回答の方向に偏ってしまうからです。
 もう少し経験のある方はこう答える。確かにウェイティングすべきでない場合もあります。この例では男女だけが問題になっていますが、男女、年代、地域、などなどと複数の層が登場すると、ウェイトを算出するためのサンプルサイズが小さくなり、極端なウェイトが得られてしまい、集計値はかえって不適切になってしまうことがあります。気をつけましょう。

 僭越ながら、この説明は間違っていると思う。少なくとも、ストーリーの半分しか語っていない。
 私はこう答えたい。ウェイティングすべきだとも、すべきでないともいえません。

シナリオA. ウェイティングすべき場合
 いくつかの架空のシナリオを考えよう。
 まず最初のシナリオ。私たちは知らないのだが、実は母集団における購入意向率は、男性で0.4, 女性で0.6であるとしよう(女性の評判のほうがよいわけだ)。関心があるのは全体の購入意向率だから、つまり0.5が「正解」である。
 このとき、標本における購入意向率は、男性で0.4のまわりの値、女性で0.6のまわりの値となる。標本は男性を多く含んでいるので、unweightedの集計値は男性の方向に向かって偏る。つまり、「正解」よりも少し低めになりやすい。
 下の図はその様子を示したものである。
case 1 - distribution
赤の山は、この架空の調査をコンピュータ上で10000回繰り返し、unweightedの集計値がどのような値をとるかを調べた実験結果を示している(微妙にデコボコしているけど、あまり気にしないでください)。赤の山の位置は、「正解」よりも少し左側にずれている。
 いっぽう青の山はweightedの集計値の分布を示している。この山の真ん中(平均)は「正解」に一致する。つまり、ウェイティングによって、標本の男女比が母集団の男女比と異なることによる偏りを取り除くことができるわけである。
 よかった、よかった。これがウェイティングにとってハッピーなシナリオである。私たちはウェイティングすべきだ。

シナリオB. ウェイティングすべきでないシナリオ
 別のシナリオ。私たちは知らないのだが、実は母集団における購入意向率は、男性で0.5, 女性も0.5であるとしよう(購入意向は性別とは無関係なわけだ)。ここでも「正解」は0.5である。
 標本における購入意向率は、男性で0.5のまわりの値、女性で0.5のまわりの値となるだろう。標本は男性を多く含んでいるけれど、男女のあいだに差がないので、unweightedの集計値は偏らない。でもweightedの集計値も偏りはしない。だから「ウェイティングしてもしなくてもよい」。と考える人が多い。
 ここに落とし穴がある。このシナリオでの集計値の分布を示す。
case 2 - distribution
確かに、赤の山も青の山も、「正解」である0.5のまわりに分布している。でもよく見ると、青の山のほうがほんの少し、下方向につぶれて、裾の広い形になっている。つまり、unweightedの集計値に比べてweightedの集計値は、「正解」からより離れた値をとりやすい。このシナリオでは、「ウェイティングしてもしなくてもよい」のではない。ウェイティングすべきでないのである。
 このように、weightedの集計値の分布は、unweightedの集計値の分布よりもばらつきが大きくなる。ウェイティングとは、集計値の分布のばらつきを犠牲にして、集計値の分布の偏りを取り除く手続きであるといえる。

シナリオC. 手に汗握るシナリオ
 シナリオAは、集計値の偏りをウェイティングによって取り除くのが望ましいシナリオであった。シナリオBは、集計値に偏りが存在せず、ウェイティングによって集計値がばらついてしまう害が顕在化するシナリオであった。
 では、その中間のシナリオ。母集団における購入意向率は性別によってわずかに異なり、男性で0.48, 女性は0.52であるとしよう。「正解」は0.5である。
case 3 - distribution
unweightedの集計値の分布は少し左にずれている。いっぽうweightedの集計値の分布は、その平均が「正解」に一致している。その代償として、分布のばらつきは少し大きい。
 さて、どちらの集計値がよいだろうか?

 そもそも、私たちが求める集計値とはどんな集計値だろうか。
 weightedの集計値は、「正解」からみて偏りがない、すなわち、長い目で見て「正解」より大きくなりやすくも小さくなりやすくもない、という美点を持っている。いっぽうunweightedの集計値は、その分布の平均からみてばらつきが小さい、すなわち、長い目で見てその値が安定しているという美点を持っている。
 でも私たちは、長い目で見て性質の良い集計値を得るために調査を行っているのではない。たった1回の調査を通じて、母集団(潜在的顧客)の性質について少しでも正しい知識を得たいと思っているのである。だから本当に重要なのは、調査によって得られた集計値が「正解」により近いことである。
 上のシナリオで、集計値と「正解」0.5との差がどうなるかを計算してみよう。図でいえば、赤の山と青の山のそれぞれについて、山を形作る無数の点のそれぞれが、(その山の平均でなくて)「正解」0.5からどのくらい離れているかを調べてみる。統計学の伝統に従い、差を二乗した値を平均して表すことにする。これを平均二乗誤差という。
table3.png
unweightedの集計値の分布は、平均0.496(男性寄りに偏っている)。分散0.002。平均二乗誤差0.0025。
 weightedの集計値の分布は、平均0.5(偏りが取り除かれている)。分散0.003(すこし大きくなっている)。平均二乗誤差は0.0026。
 僅差ではあるが、unweightedのほうが誤差が小さい。つまり、このシナリオでは、私たちはウェイティングすべきでない。

シナリオのマップ
 これらの3つのシナリオでは、母集団の購入意向率が男女それぞれについてわかっていた。実際の調査は、母集団の性質がわからないからこそ行うわけだから、上のような図は手に入らない。
 しかし、こういう図を描くことはできる。
map3.png
この図は、母集団の男女比は5:5, 母集団全体での購入意向率は0.5, 標本サイズは100, という3点だけを決めて描いている。
 横軸は母集団における男性の購入意向率を示す。3つのシナリオでは、それぞれ横軸0.40, 0.50, 0.48であった。縦軸は男性の標本サイズを示す。3つのシナリオではすべて60であった。図中の×印はシナリオCを示している。
 図中の色は、unweightedの集計値の平均二乗誤差と、weightedの集計値の平均二乗誤差との差を表している。青のエリアは、ウェイティングした方が誤差が小さい(ウェイティングすべきである)。赤のエリアは、ウェイティングしないほうが誤差が小さい(ウェイティングすべきでない)。

 この図からいろいろなことがわかる。

「極端なウェイト」
 ウェイティングについての説明をみていると、ウェイティングの弊害として「ウェイトが極端なときは集計値の誤差が大きくなる」という点が挙げられていることが多い。この説明は当たっている面もある。上の図でいえば、ウェイト値は極端な値になるのは天井付近や床付近である。たしかに赤くなっている。
 しかし、赤のエリアは天井や床だけでなく、もっと内側にも広がっている。たとえば、さきほどの3つのシナリオはすべて、ウェイトは男0.833, 女1.250であった。これを「極端なウェイト」と呼ぶのは難しいだろう。それでも、ウェイティングによって誤差が大きくなる事態は生じた。
 このように、ウェイトが極端な値になるかどうかは、ウェイティングの是非にとって本質的問題ではない。

標本サイズとの関係 
 この図は状況によって変わる。全体の標本サイズを500に増やしてみよう。赤のエリアは小さくなってしまう。
map4.png
このように、標本サイズが大きいとき、weightedの集計値はunweightedの集計値に比べて有利になる。
 これは次のように考えるとわかりやすいだろう。ウェイティングとは、集計値の偏りを取り除く代わりにばらつきを拡大する手続きである。集計値のばらつきは標本サイズと関係している。標本サイズが小さいときは、集計値のばらつきがもともと大きいので、ウェイティングによる拡大の影響が深刻でなる。それに対し、標本サイズが大きいときは、集計値のばらつきがもともと小さいので、ウェイティングによる拡大の影響はさほど問題にならない。
 このように、ウェイティングの是非を判断する際のひとつの重要なポイントは、全体の標本サイズである。

母集団特性との関係 
 もっと評判の悪い製品だったらどうなるか。標本サイズを100に戻し、今度は母集団全体の購入意向率を0.2に下げてみよう。
map5.png
興味深いことに、赤のエリアは左右非対称になる。よくみると、男性の標本サイズのほうが大きい場合(図の上半分)、男性のほうが購入意向が高い場合にはウェイティングすべきエリアが広く(図の右上)、男性のほうが購入意向が低い場合にはウェイティングすべきでないエリアが広い(図の左上)。
 これは割合という集計値の持っている性質に由来している。購入意向率とはすなわち購入意向者の割合である。割合という集計値は、真の値が0.5に近いときにばらつきが大きくなる性質を持っている。層別抽出においては、ばらつきが大きい層の標本サイズを増やしておけば、全体の集計値のばらつきを小さく押さえることができる。図の右上と左下がそれである。集計値のばらつきがもともと小さいぶん、ウェイティングにとって有利になる。

この実験からわかること 
 この簡単な実験から、次の点がわかる。ウェイティングすべきかどうかという問いに簡単な答えはない
 標本抽出デザインだけからは、ウェイティングすべきだともすべきでないともいえない。たとえばこんな批判を耳にすることがある。「この調査、性別と年代で均等割り付けしている調査なのに、なんでウェイトバック集計していないんだ?」 こういう非難はいささか早計に過ぎると思う。非比例の層別抽出でも、あえてウェイティングしないという判断があってしかるべきである。ウェイティングは集計値の誤差を減らすこともあれば増やすこともあるのだ。

ではどうすればいいのか
 ウェイティングすべきかどうかという問いに簡単な答えはない。でも、あれこれ考えた上で、それなりの答えを出すことはできる。
 上の実験から、集計値の誤差という問題について考えるための簡単なガイドラインを手に入れることができる。その調査で測定した変数のうち重要な変数(たとえば購入意向)について、次の3点に注目すれば良い。

この3つの注意点に、冒頭に挙げた3つのポイント(適切なウェイトを構成できるか、能力とやる気があるか、分析の諸前提と整合するか)を加えれば、ウェイティングの是非を判断するための材料が揃うのではないかと思う。

シミュレータつくっちゃいました
 ところで、この文章を書いたもうひとつの理由、ささやかなほうの理由は...この文章で説明したシミュレーションを行うためのwebアプリを作ってしまったためである。図表はこのアプリでつくりました。
capture.png URLは今後変わると思いますのでご注意ください。
 Shinyというシステムの使い方を勉強するつもりで作り始めて、たった半日でできてしまった。おそろしい時代になったものだ。
 ご関心あらばどなたでもお試しくださいまし。なかなか楽しいのではないかと思います。

2015/03/11 追記: アプリは停止しました。お粗末様でございました。

「ウェイトバック集計」すべき場合とすべきでない場合がある(ということを説明するためのwebアプリを作ってしまった)

2014年2月16日 (日)

 Mplus 7.1に搭載された多群因子分析のための新機能が、魅力的なんだけどちょっとわかりにくいので、自分のための備忘録としてメモをとった。誰かの役に立たないとも限らないからブログに載せておく。きっとどこか間違えていると思うので、ほんとにお使いになる方はマニュアルをごらんください。

 この新機能 「MODEL=」は、「GROUPING=」か「KNOWNCLASS=」をつかった多群モデルで、測定不変性についてのカイ二乗検定を自動的にやってくれるというもの。従来は「DIFFTEST=」を使ったりなんだり、非常に面倒だった。

 「MODEL=」が使えるのはCFAかESEMのモデル。BY文は一回しか使えない(一次因子しか定義できない)。部分測定不変性は扱えない。変数のタイプは以下に限られる。

つまり、inflatedな変数(日本語ではなんて訳すのかしらん... 「ゼロ過剰モデル」っていうときの「過剰」ですね)、名義変数、連続時間生存変数、負の二項分布の変数、は扱えない。

 「MODEL=」は値として、CONFIGURAL, METRIC, SCALARの並びをとる。たとえば「MODEL = CONFIGURAL METRIC SCALAR;」という風に指定する。

 さて、「MODEL=」を指定すると正確にはどんなモデルが構築されるのか、という点だが、これが結構ややこしい...
 以下、因子分散の指定を、どこかの項目の因子負荷を固定するやり方で行う場合を[A]、どこかの群の因子分散を固定するやり方で行う場合を[B]と略記する。なお、Mplusが勝手に [B]を適用するとき、「GROUPING=」を使っている場合は最初の群、「KNOWNCLASS=」を使っている場合は最後の群が選ばれる。

あれ? 二値や順序カテゴリカル変数でベイズ推定したらどうなるんだろうか。今度試しにやってみよう。

Mplus 7.1 で多群因子分析をやるときの魅惑の新機能 MODELオプションについて

2013年5月19日 (日)

 このブログを書き始めてからずいぶん経つが、所詮は自分のための備忘録に過ぎない。いちばん熱心な読者はたぶん私だ。書店の棚に平積みになっている連載マンガのこの巻を、俺はもう読んだだろうか? と疑問に思った際、このブログはとても重宝する。店頭でスマフォと格闘している姿は、あまり格好のよいものではないけれど。
 そもそも私には世の中に訴えたいことが特にない。とりたてて語るべき知恵もない。アクセスログを調べると、過去のいくつかの記事が予想外に多くの方に読まれていることに驚くが(片側検定について書いた記事とか)、もっときちんとした専門家がお書きになったものをご覧いただいたほうが良いわけで、少し申し訳なく思っている。

 これから書くのは、このブログを書き始めて以来ほとんどはじめて、純粋に自分以外の誰かのために載せる記事だ。特に、検索エンジン経由でたどり着く方に向けて載せる記事だ。
 振り返ってみると、私は何年か前から、だいたい年に数回のペースで、以下の内容をどなたかに説明している。聞き手にとっては面倒で退屈な内容だ。説明に成功することもあれば失敗することもある。偉大なるインターネットの力によって、説明の回数を減らし成功率を上げたい。
 というわけで、ずっと前にお問い合わせに応じて書きためていた文章を載せておくことにする。どうかこの記事が、正しい内容を正しいやりかたで説明していますように。そしてgoogleの検索結果の上の方に出てきますように。

 以下の内容はとても単純である。たった一言で要約できる。

 統計ソフトがいう「ウェイト」は、調査でいう「ウェイト」ではない。

 もう一度書きます。SPSSとかSASとか!そういう統計ソフトの! ウェイトとか重み付けとかっていうのは! 調査でいうところの! いわゆるウェイトのことではないです!
 。。。このくらい大声で書いておけば、検索されやすくなるだろうか。

調査における「ウェイティング」

 市場調査の会社にお世話になって驚いたことのひとつは、結果の集計に際してウェイティング(確率ウェイティング)を実に良く使う、ということである。親の仇かというくらいによく使う。

 確率ウェイティングとはこういう話だ。
 ある集団の性質について調べるため、その集団に属する人を対象に調査を行った。調査対象者がこの集団を互いに独立に偏りなく反映していると仮定し、得られたデータから統計的推測を行いたい。
 さて、私たちはこの集団(対象母集団)における男女比が5:5であると知っている。しかしなにかの事情によって、調査対象者の男女比は7:3にならざるを得なかった。この事情にはいろんな種類がありうるが、たとえば調査設計時の事情が考えられる(比例割当でない層別抽出)。

 この調査対象者の回答をそのまま集計し、平均や割合を求めてしまうと、その値は男の回答をより強く反映してしまい、母集団の推測としては歪んだ値になってしまう。
 これを避けるために、男の回答には小さな重み(ウェイト)、女の回答には大きな重みを与えて集計する。これが調査でいうところのウェイティング、すなわち確率ウェイティングである。
 なお英語では、probability weighting, survey weighting, sampling weightingなど、いろいろな呼び方が用いられている。検索するときに困る。

 市場調査の業界団体である日本マーケティング・リサーチ協会がまとめた「マーケティング・リサーチ用語辞典」を見ると、「ウェイトつき集計」という項で確率ウェイティングが説明されている。市場調査に限らず、広義の社会調査全般において、ウェイト、ウェイティング、ウェイト・バックという言葉は確率ウェイティングを指して使われることが多いと思う。

確率ウェイティングと拡大推計

 確率ウェイティングは拡大推計とごっちゃにして語られることもある。

 たとえば、前掲の業界団体による市場調査会社向けガイドラインには、ウェイティングについての注意点を挙げている箇所がある。

一般的にウエイトバックが必要となる理由は
1. 層化サンプリングの段階でウエイト付けする場合
2. 回収率の違いにより標本構成に偏りが出てそれをウエイト付けにより補正する場合
上記2点が多いと考えられるが、[中略] ウエイト付けをした場合、主な表すべてにおいて、ウエイト付けした基数(ベース)としない基数(両者を明確に区別)を明示すること。

前半はどうみても確率ウェイティングに関する説明である。ところが後半では、「ウェイト付けした基数」という文言が登場する。ウェイティングしようがしまいが標本サイズは標本サイズなのであって、これは意味不明である。

 実は、「1000人の対象者のうち3000人がこの製品を買いたいと答えた、市場には潜在顧客が1000万人いるから、発売すればきっと300万人が買ってくれるだろう」というような単純な拡大推計を指して「ウェイティング」ということもある。
 ガイドラインの文言は、確率ウェイティングと拡大推計の両方について同時に述べているようである。

 確率ウェイティングは標本の性質に応じて必要となるもので、あらゆる統計的分析に影響する。いっぽう、拡大推計は標本の性質とは無関係に要請され、要約統計量のなかでも合計にしか影響しない。クロス集計の手順のみについて考えれば似たようなものかもしれないが、このふたつは分けて考えたほうが良いと思う。

ウェイティングした分析はできますか?

 調査データの分析において確率ウェイティングはどんなときに必要か。仮に必要だとして、ウェイトはどうやって求めるか。
 これはなかなか難しい問題で、容易な答えはない。悩みの種である。あまりに悩んだ末、高名な専門家による解説を全訳してしまったことさえある。若気の至り、窓の雪。しかしその話は脇に置いておく。

 本題は、確率ウェイティングのためのウェイトが各ケース(ローデータの行)に対してすでに付与されているとき、そのウェイトを用いながらさまざまな統計的推測を行うことができるか、という点である。つまり、統計ソフトを使い、データに確率ウェイティングを行いつつ、差の検定を行ったり、相関を求めたり、回帰分析や因子分析を行ったりすることができるか(行うべきかどうかは別にして)、という点である。

 私は繰り返し説明してきた。できません。十中八九できません。なぜなら、お使いの統計ソフトがいうところの「ウェイト」は、私たちが思っているような意味でのウェイトではないからです。

手計算による確率ウェイティング

 ごく簡単な例について考えよう。
 母集団における男女比は5:5。事情があって、標本における男女比が7:3となるように設計した。母集団は十分に大きいものとする。非現実的ではあるが、標本サイズは男7, 女3, 計10であるとしよう。
 この10人が、ある変数について次のような値を持っている。

男性1, 2, 2, 3, 3, 4, 4
女性1, 2, 3

 このデータに基づき、母集団におけるこの変数の平均(母平均)について推測したい。

 まずは、統計ソフトを使わずに考えてみよう。ひとまず男性に注目する。標本平均は、(1+2+2+3+3+4+4)/7=2.714。つまり、母集団における男性の平均は2.714と推定される。
 この推定はどのくらいあてになるだろうか。標本平均の分散を求めよう。統計学の初級の教科書に書いてあるように、標本平均の分散は (不偏分散)/(標本サイズ) である。ええと、不偏分散は1.238, 標本サイズは7 だから、標本平均の分散は1.238/7=0.177だ。
 同様に女性についても電卓を叩いてみると、標本平均は2.000, 標本平均の分散は0.333である。

 では、全体の平均について考えよう。いま、男性の平均が2.714、女性の平均が2.000と推定されていて、かつ母集団のうち男性は5割、女性は5割であることがわかっているのだから、男女をあわせた平均は 2.714*0.5 + 2.000*0.5 = 2.357 と推定できる。
 この推定はどのくらいあてになるか。推定量の分散は、{(男の標本平均の分散) * (男の割合)^2 + (女の標本平均の分散) * (女の割合)^2} となる (ご不審の向きは、L. Kish "Survey Sampling"(Wiley) の80-82頁、豊田秀樹「調査法講義」(朝倉書店) の131-132頁をご覧ください)。ええと、0.177*0.25 + 0.333*0.25 = 0.128である。

 以上が「正解」である。表にまとめておく。なお、標本平均の分散の平方根は標準誤差と呼ばれ、広く用いられているので、それも添えておく。

母平均の推定値(標本平均)その分散標準誤差
男性2.7140.1770.421
女性2.000 0.3330.577
全体2.3570.1280.357

手計算による確率ウェイティング(ケース・ウェイトを用いて)

 この計算を、個々のケースにウェイトを与える形で表現しなおしてみよう。

 まず、各ケースにウェイトを与える。このウェイトは、ケースが母集団から抽出される確率の逆数に比例した値であればなんでもよいのだが、一番わかりやすいのはこういう求め方だ。母集団における男の割合は5割、標本では7割だから、男のウェイトは5/7=0.714。同様に、女のウェイトは5/3=1.667。ウェイトの合計は10となる。

 母平均の推定値、すなわち「ウェイティングされた標本平均」を求める際には、個々の値にウェイトを掛けながら合計し、ウェイトの合計で割る。{1*0.714 + 2*0.714 + ... + 2*1.667 + 3*1.667) / 10 = 2.357。
 上記の結果とぴったり一致する。計算式を変形しただけだから、当たり前である。

SPSSで試してみると

 では、統計ソフトを使って試してみよう。みんな大好きな IBM SPSS Statistic での例を紹介する (ver.19で試しました)。下図のデータを用意する。クリックすると大きくなるはずです。

Image.png

SPSSのメニューには[データ]-[ケースの重み付け]という機能がある。こいつはきっと確率ウェイティングの機能にちがいない、というわけで、「ケースの重み付け」を選び、変数weight を指定してみよう。次に、[分析]-[グループの平均]で変数valueの平均を求める。[オブション]で「平均値の標準誤差」を選んでおく。アウトプットは...

Image2.png

 平均値は2.357。「正解」と一致する。ところが、平均の標準誤差は0.332。「正解」よりも少し小さい。

 平均の標準誤差なんてどうでもいいです、などというなかれ。たいていの推測統計手法は標準誤差と関係している。みんなが死ぬほど大好きな検定だって、もちろんそのひとつだ。標準誤差が誤りなら、検定の結果だって誤りだ。これはとても大事な話なのである。

頻度ウェイティング

 なぜこのような結果になるのか。理由は拍子抜けするほどに簡単だ。
 SPSSでいう「ケースの重み付け」は確率ウェイティングのことではない。それは単に「ウェイトの数だけこのケースをコピーしてくれ」ということである。

Image3.png

 たとえば上の3行のデータを与え、変数weightで「ケースの重み付け」を行うと、SPSSはこう考える。変数Xが10である行が100行、20である行が120行、30である行が130行、全部で350行のデータなんですね、わかりました、と。

 こうした機能は「頻度ウェイティング」と呼ばれることがある。頻度ウェイティングはデータ行列のサイズを小さくするための工夫に過ぎない。

 頻度ウェイトによってコンパクトに表現されたデータ行列を一発で扱う方法について考えてみよう。ある変数の標本平均を求める際には、個々の値に頻度ウェイトを掛けながら合計し、頻度ウェイトの合計で割ればよい(なぜなら、頻度ウェイトは本来のローデータにおけるケース数だから)。
 この式は、幸か不幸か、確率ウェイティングのもとで母平均の推定値を求めるための式と、た・ま・た・ま、同じ形になっている。そのせいで、SPSSがウェイトを頻度ウェイトとみなして求めた標本平均は、同じウェイトが確率ウェイトだった場合の母平均の推定値 2.357 と、た・ま・た・ま、一致する。だからといって、SPSSが確率ウェイトを正しく扱えるわけではない。

 なお、SPSSにとってはウェイトは頻度なのだから、それは整数でないとおかしい。さきほどはウェイトとして整数でない値を与えたが、よくもエラーにならないものである。実際、上記で試した[分析]-[グループの平均]は、整数でないウェイトをむりやり頻度として捉えるが、他の機能は、ウェイトを整数に丸めてしまったり、そもそも無視してしまったりすることがあるらしい。

分析ウェイティング

 もうひとつの有名な統計ソフトであるSASでは、また少し話が違う。
 SASでは各プロシジャのweight 文でケース・ウェイトを指定できる。その扱いはプロシジャによっても設定によってもちがうのだけれど、通常は、「このデータの行はサブグループの平均を表していて、それぞれのグループのサイズがウェイトで表されているのね」という意味に解される。

Image3.png

 たとえば、さきほどの3行のデータをSASに与え、weight列をケース・ウェイトに指定すると、通常は次の意味になる。「サブグループ1はサイズ100, Xの平均は10だった。サブグループ2はサイズ120, Xの平均は20だった。サブグループ3はサイズ130、Xの平均は30だった。」
 こうしたウェイティングは「分析ウェイティング」と呼ばれることがある。これもまた、確率ウェイティングではない。

統計ソフトにおける「ウェイティング」とは

 統計ソフトの中には確率ウェイティングを扱うことができるものもある。

 そういうソフトでないかぎり、お使いの統計ソフトでは ... もっと具体的にいえば、SPSS Statistic BaseならびにAdvanced Statistics、SAS/BASEならびにSAS/STATのsurvey系プロシジャを除く全プロシジャ、その他もろもろの統計ソフトでは ... 確率ウェイティングを伴う分析を行うことはできない。メニューに「重み付け」という機能があっても、プロシジャにweight文というのがあっても、それは確率ウェイティングではない。

 考えてみれば、統計ソフトのユーザは多方面にわたるが、確率ウェイティングをしたいと思う人はせいぜい市場調査・社会調査の関係者くらいであろう。一般的な統計ソフトのメニューにある「重み付け」が、調査実務家の考える意味でのウェイティングでないとしても、文句をいえる筋合いではない。

 ややこしいのは、確率ウェイティングに対応していないソフトに誤って確率ウェイトを与えて要約統計量を求めたとき、平均と割合だけは期待通りの値が出力されることがある、という点である。繰り返しになるが、それはた・ま・た・まである。他の結果は誤っている。
 統計ソフトがいう「ウェイト」は、調査でいう「ウェイト」ではない。

あえて教訓を求めるならば

 最初に書いたように、私には世の中に訴えたいことは特にない。とりたてて語るべき知恵もない。ただ単に、日頃繰り返している説明を文章にしておき、仕事を楽にしたいだけである。
 しかし、あえて教訓らしきものを引き出すとすれば、こういうことはいえる。

 私の説明にオリジナリティはない。webを探せば、もっとわかりやすい説明がみつかる。たとえばこちらこちら。統計ソフトのアウトプットをこまかーく調べた面白い研究もある(PDF)。
 さらに、SPSS Statisticsを起動し[ケースの重み付け]ダイアログの[ヘルプ] ボタンを押すと、次の説明が現れる。

[ケースの重み付け] では、統計分析を行うため (複製をシミュレートすることにより) ケースに異なる重みを付けることができます。重み付け変数の値は、データ ファイル内の 1 つのケースが表す観測数を示していなければなりません。
いささか不親切ではあるけれど、簡にして要を得た説明である。この記事で長々と書いてきたことはすべて、ソフトのマニュアルに、はっきりと書いてあることなのだ。

 にもかかわらず、推定や検定やもっと複雑な統計的分析の際に、SPSS Statisticsの[ケースの重み付け]で確率ウェイティングを行おうとしている方に、私はこれまでたくさん出会ってきた。それはもう、びっくりするくらいにたくさん。そのなかには、すごく優秀な分析者もいらっしゃったし、調査実務に携わって何十年というベテランの方もいらっしゃった。

 というわけで、教訓はこうだ。実務家の豊富な経験なるものは、必ずしもあてにならない。統計ソフトの使い方を先輩に教わったり、誰かの書いたブログを読んでいる暇があったら、そのソフトのマニュアルを読んだほうがよいのではないかと思います。

統計ソフトの「ウェイト」は調査の「ウェイト」ではない

2013年4月17日 (水)

 連続変数のベクトル y を指標として持つ潜在クラスモデルについて考える(いわゆる潜在プロファイル・モデル)。話を単純にするため、共変量はなし、指標はすべて局所独立、条件つき分布に多変量正規性を仮定する。
 対象者 i がデータ y_i を持っているとしよう。この人の所属クラス C が k である事後確率は
 P(C=k | y_i) = P(C=k) [y_i | C=k] / [y_i]
ただし [y_i | C=k] は、クラス k の平均ベクトルと共分散行列を持つMVNである。みよ、大自然の単純な美しさを。

 いったいなにを云わんとしているのかというと... 指標がV1とV2の2つしかない、2クラスの潜在クラスモデルを推定したとしよう。mplusが吐いたパラメータ推定値が
Latent Class 1: Means V1 3.248 V2 5.626; Variances V1 1.985 V2 8.243
Latent Class 2: Means V1 8.845 V2 5.298; Variances V1 1.985 V2 8.243
Categorical Latent Variables: Means C#1 0.220
であったとする。さて、V1=6, V2=4のオブザベーションがあったら、そのクラス所属確率は?
 まず事前確率について。mplusはカテゴリカル潜在変数の平均を最後のクラスで 0 とするので、P(C=1) ∝ exp(0.22) = 1.24, P(C=2) ∝ exp(0) = 1。足しあげて割合にすると、0.55, 0.45である。
 次に尤度。分子のほうだけ考える。Excel風に書くと、C=1ならばnormdist(6, 3.248, sqrt(1.985)) * normdist(4, 5.256, sqrt(8.243)) * 0.55 = 0.00351。C=2ならば0.00206。これを足しあげて割合にして、0.63と0.37。
 これが事後確率である。mplusのSAVEDATAで出したものとぴったり、ぴったり一致する。みよ、大宇宙の壮大な神秘を。

 なにが云いたいのかというと... マーケット・セグメンテーションの手法として、消費者調査データに基づく対象者分類を行った場合、あとで別の対象者に同じ調査項目を聴取し、その回答に基づいてその人がどのセグメントに属するかを判別したい、ということがある。ものすごくよくある。細かい商売ですみません。
 で、分類の際にk-means法を使っていた場合は判別関数をつくるのは容易だが、潜在クラスモデルのような確率モデルをつくってしまうと、あとの判別が大変だ。などと思い込んでいる人がいる。愚かなり。実に愚かなり。
 というか、私自身がついついそう思いこんでいたのである。このたび用事があってこの件について考える羽目になり、なんとか自分で考えずに済ませられないものかと考えたが逃げ道はなく、どんよりした気分で机の上に資料を揃え、気分転換用のお菓子も用意し、深呼吸してから考え始めた。15分後、菓子に手をつける間もなく、あまりの簡単さにあっけにとられている私がいた(お菓子は後で食べましたが)。私が愚かでした、反省してます。反省のあまりブログに記録する次第である。いっけんややこしそうだからといって、思考停止してはいけないのだ。

潜在クラスモデルの所属確率について (反省の弁)

2012年11月14日 (水)

 流れ流れて市場調査の会社に拾って頂いた当初,いろいろ面食らうことが多かったのだが,実験計画法? そんなのどうでもいいですよ,科学者じゃないんだから,と真顔でいわれたことがあって,このときも言葉を失った。たぶん実験ということばは,実験室やら白衣やら,おぞましくもアカデミックななにかを連想させ,人を思考停止に導いてしまうのであろう。
 たとえば,調査対象者をランダムに2群に分け,一方の群にはある製品に現行のパッケージラベルを貼って提示し,他方の群には同じ製品に新しいパッケージラベルを貼って提示する。で,好意度なり購入意向を比較する。これだって立派な実験である。
 ここで「パッケージラベル」を要因 factor という。登場するパッケージラベルは「現行ラベル」と「新ラベル」の二種類あるが,このそれぞれを水準 level という。つまり,このごくありふれたパッケージ・テストは,1要因2水準の実験である。

完全実施要因計画と一部実施要因計画

 パッケージラベルの違いは製品への好意度に影響を与えるか。この点について推測する際には,要するに,それぞれのパッケージラベルに対して好意度を聴取し回答を集め,2グループの回答の集まりを比較すればよい。話の都合上,一方の水準について測定したこと(つまり,一方のパッケージラベルについて好意度を聴取したこと) を +1 と表現し,他方の水準について測定したことを -1 と表現することにする。

 要因の数がもっと増えた場合について考えよう。たとえば,パッケージラベルの影響だけでなく,蓋の色による影響も,パッケージの形状による影響も同時に調べたい,というような場合である。3要因(A, B, Cと呼ぶ)がそれぞれ2水準を持つとすれば,水準の組み合わせは2*2*2=8個。この8個の製品をつくり,それぞれについて好意度を聴取すれば良いわけだ。以下,以下の8行の表のイメージで考えることにする。

製品番号ABC
1-1-1-1
2-1-1+1
3-1+1-1
4-1+1+1
5+1-1-1
6+1-1+1
7+1+1-1
8+1+1+1

 この表の行,つまり水準の組み合わせのことを,とりあえず製品と呼ぶことにする。英語ではラン run と呼ぶことがあるが,なんと訳せばいいのかわからない。

 上の表のように,全ての属性の全ての水準の全ての組み合わせからなる全ての製品を使って実験するのが,ひとつの理想である。こういう実験計画を完全実施要因計画 full factorial design という。しかし,すべての製品を用意するのは大変だ。それに,要因の数が4つ, 5つ, ... と増えたり,2水準の要因だけではなく3水準の要因や4水準の要因が登場すると,可能な製品の数は膨大になり,そのすべてを使った実験は現実的でなくなってしまう。いくつかの少数の製品だけをうまく選び出して実験したい。これを一部実施要因計画 fractional factorial design という。

推定精度のシミュレーション

 一部実施要因計画の場合,用いる製品をどうやって選べばよいだろうか。
 四の五の言わずに,かんたんなシミュレーションをしてみよう。まず,「正解」を決めてしまう。好意度を7件法評定(1~7点)で聴取する場合について考えよう。各製品はある決まった好意度の高さ(「正解」)を持っている,と想定する。値はなんでもよいのだが,仮に,8製品を通じた好意度の平均は4.0ポイント,Aが-1のときよりも+1のときのほうが好意度が2.4ポイント高い,同様に B, Cではそれぞれ0.6, 0.3の差がある... ということにしよう。こうしてできた「正解」は以下の通り。

製品番号ABC正解
1-1-1-11.9
2-1-1+12.5
3-1+1-13.1
4-1+1+13.7
5+1-1-14.3
6+1-1+14.9
7+1+1-15.5
8+1+1+16.1

 では,完全実施要因計画に基づいた架空の実験を行おう。対象者120人を8つの製品にランダムかつ均等に割り当てる。対象者には,自分が割り当てられた製品について好意度を評定してもらい,それでお役御免とさせていただく (市場調査ではこういう実験計画をモナディック・デザインと呼ぶことが多い)。各製品に対する人々の好意度はその製品が持っている「正解」の周りをばらつく,と仮定しよう。データ解析の世界でのちょっとした伝統に従い,ばらつきは左右対称なある決まった形状をとる(N(0,2)に従う)ということにしておく。
 以上の仮定に従い,コンピュータ上で疑似的な回答データをつくる。現実世界と違って一瞬にしてデータが集まるので気分が良い。よくよく見ると,7件法で聴取したはずなのに回答は整数になっていないし,負の値や7以上の値も含まれているけれど,この際それは気にしないことにしよう。
 得られた回答の集計結果をご報告しよう。このたびはどうも運が悪かったようで,各組み合わせに対する回答の平均は「正解」からかなりズレている。

製品番号ABC正解回答平均
1-1-1-11.92.1
2-1-1+12.53.4
3-1+1-13.13.3
4-1+1+13.73.3
5+1-1-14.33.9
6+1-1+14.94.8
7+1+1-15.54.9
8+1+1+16.16.7

 こうした実験を行った場合,データの分析には分散分析と呼ばれる手法を用いることが多いので,ここでもそれに従おう。主効果モデルと呼ばれる分散分析モデルを適用する(そのことの是非については,のちほど)。このモデルでは,全製品の平均, Aが+1であるときの平均に対する増加分, Bが+1であるときの増加分, Cが+1であるときの増加分,の4つをデータから推定する。以下ではこれらを「切片」「Aの係数」「Bの係数」「Cの係数」と呼ぶ。その「正解」は,4.0, 1.2, 0.6, 0.3である。
 推定の結果は,それぞれ3.92, 1.16, 0.63, 0.38となった。Aの影響が過小評価され,B, Cの影響が過大評価されている。

 この実験はたまたま運が悪かったのかもしれない。そこで,同じ設定の実験を10000回繰り返す。Excelなり統計ソフトウェアなりに慣れている人なら,さほど難しくはない作業だ。
 各実験において,得られた4つの推定値と「正解」とのズレを調べ,10000回の実験についてズレを集計する。データ解析のちょっとした伝統に従って,ここでは推定値と「正解」の差の二乗の平均を集計値とする。推定の精度が低いほど大きな値になる。結果は... 切片0.033, Aの係数0.034, Bの係数0.033, Cの係数0.034となった。

 完全実施要因計画ならば,ズレの大きさはこれくらいだ,という見当がついた。今度は,一部実施要因計画を用いて実験してみよう。あれこれ考えるのが面倒なので,片っ端から実験してみることにする。
 ええと,8製品のうち7つを用いる実験計画は8パターン,8製品のうち6つ用いるなら28パターン,8製品のうち5つ用いるなら56パターン, 8製品のうち4つ用いるなら70パターン,ここまで計162パターン。さらに,8製品のうち3つなら... いや,それはさすがに虫が良すぎるだろう。実際,ちょっと考えてみればわかるのだが,各製品が持つ好意度の「正解」について予備知識がない状況で,3つの要因がもたらす影響を評価するためには,少なくとも4つの製品について調べる必要がある。

 というわけで,対象者数を120人に固定したまま,8つの製品から4つ以上を用いる計162パターンのそれぞれについて,完全実施要因計画と同様に10000回の架空実験を繰り返してみる。さすがにExcelではちょっと面倒だが,なんらかの統計ソフトウェアの使用方法に慣れていれば,それほど難しい作業ではない。これをSASで試した時は,20分くらいで計算開始にこぎつけた。このたびRで書き下ろしたら,計算までに小一時間かかってしまったが,これは私が下手だからで,上手い人ならあっという間だと思う。

plot_zoom_png.png

 結果を上図に示す(クリックで拡大表示)。162パターンのなかには,とてもじゃないが無理,というようなパターンも登場してしまうので(例, 製品{1, 2, 3, 4}。どれもAが-1である),そのようなパターンは除外してある。
 チャートの底部を水平に走る黒点線が,完全実施要因計画の場合のズレの大きさである。製品数 7 の場合のズレの大きさは,その少し上の緑色。つまり,製品を1つ減らしたせいで,「正解」からのズレが少し大きくなっている。以下,6製品, 5製品, 4製品と,製品数を減らすにつれてズレがますます大きくなることがわかる。
 。。。いや,ちょっと待て! チャートの底部をよく見てほしい。目をこらすと,ああ,なんということだろうか! 黒点線の周りに,かすかな赤い輝きがみえるではないか!!

直交計画の威力

 わざとらしい驚き方で恐縮だが,実は製品数4の一部実施実験計画のなかに,たったふたつだけ,完全実施要因計画に匹敵する推定精度を誇る,特別なパターンが存在するのである。それは,製品{1,4,6,7}(下表のアミカケ部分)を用いる計画,そしてその残りの製品{2,3,5,8}を用いる計画である。

製品番号ABC
1-1-1-1
2-1-1+1
3-1+1-1
4-1+1+1
5+1-1-1
6+1-1+1
7+1+1-1
8+1+1+1

 このふたつの実験計画はなぜ特別なのか。上の表のアミカケ部分をよくみると,このふたつの実験計画は次の2つの特徴を備えていることがわかる。

  1. A,B,Cのどれをみても,+1と-1が2行ずつはいっている。つまり,ある要因のなかで,各水準を持っている製品数が水準間で同じである。この性質を,釣り合いが取れている balanced という。
  2. たとえば要因AとBを取り出してみよう。Aが-1のとき,Bは-1がひとつ,+1がひとつである。Aが+1のときも,Bは-1がひとつ,+1がひとつである。つまり,どの2要因を取り出しても,「一方が+1だったら,そうでない場合に比べて,他方は+1になることが多い(少ない)」という関連性がない。この性質を直交している orthogonal という。

 少々ヤヤコシイのだが,この2つの性質を同時に兼ね備えた実験計画のことを特に直交計画 orthogonal design という。ここでいう「直交」とは,さきほどの「直交」とはちょっと別の意味で用いられている(計画行列における直交性を指している)。

 考えてみれば,これは不思議な話だ。私たちはいま,いくつかの製品についての消費者の好意度の回答を集め,それを分析することで,好意度の背後にあるメカニズムを探ろうとしている。だったら,調べる製品の数は多ければ多いほうがよい,と考えるのが自然であろう。しかし実際には,製品の選び方にはコツがあり,多ければ良いとは限らないのである。
 ね? 実験計画法って必要でしょう?

主効果と交互作用

 うまい話には大抵裏がある。直交計画の推定精度にも,それ相応の代償が伴うのだけれど,その話の前に,まず直交計画の性質について紹介しておこう。
 話をいったん完全実施要因計画に戻す。たいていの場合,私たちは個々の要因が好意度にどの程度影響しているかを知りたい。たとえば,要因Aの2つの水準のあいだで,好意度にどのくらいの差が生じるかどうかを知りたい。これを要因Aの主効果 main effect という。その推測に際しては,Aが+1である4製品に対する好意度の集まりと,-1である4製品に対する好意度の集まりを比較すればよい。
 しかし,世の中は時としてわれわれが期待するよりも複雑に出来ている。私たちがそれを知りたいかどうかは別にして,現実には,「要因AとBが両方+1だったときにのみ好意度が高くなる」... というような影響が生じているかもしれない。これを要因Aと要因Bの交互作用効果 interaction effect という。その推測に際しては,A列とB列の符号が合致している4行と,符号が異なっている4行を比較すればよい。

直交計画の分解能

 交互作用効果の早見表をつくるつもりで,A列とB列を掛けた列(AB列)を表に追加しておこう。AとBの交互作用をチェックしたいときは,AB列が+1である4行と-1である4行とを比較すればよい。ついでに,「要因1と3の交互作用」(AC列),「要因2と3の交互作用」(BC列),「要因1,2,3の交互作用」(ABC列) も作っておくことにしよう。こうしてできた表は,実験計画法の世界でL8直交表と呼ばれる有名な道具となる。

ABCABACBCABC
-1-1-1+1+1+1-1
-1-1+1+1-1-1+1
-1+1-1-1+1-1+1
-1+1+1-1-1+1-1
+1-1-1-1-1+1+1
+1-1+1-1+1-1-1
+1+1-1+1-1-1-1
+1+1+1+1+1+1+1

 さて,話を一部実施要因計画に戻そう。さきほどみつけた直交計画とL8直交表を見比べると,L8直交表でABC列が+1である4行,ないし-1である4行を抜き出したのが,さきほどの直交計画であることがわかる。いいかえれば,さきほどの直交計画をつくるためには,L8直交表でABC列が+1ないし-1である行を抜き出せばよい。つまりこの例で,ABC列は「ありうる製品のうち,一部実施計画に取り入れるべき製品はどれか」を表していることになる。
 この「ABC」をこの一部実施要因計画のジェネレータ generator という。また,ジェネレータの桁数を分解能 resolution と呼ぶ。なぜかわからないがローマ数字で書く習慣がある。さきほどの直交計画は分解能IIIと呼ばれる。英語でいうとかっこいいですね。自分がなにかすごくつまらない仕事をしているのではないかと気がめいるときは、小声で「よし、レゾリューション・スリーで実験だ」などと呟けば、すこしは虚栄心が満たせるかもしれない。

分解能IIIの直交計画,その欠点と役割

 分解能IIIの直交計画にはある深刻な欠点がある。もう一度,よくみてみよう。

製品番号ABC
1-1-1-1
2-1-1+1
3-1+1-1
4-1+1+1
5+1-1-1
6+1-1+1
7+1+1-1
8+1+1+1

 アミカケされた4行では,Aが-1である2行ではBとCの符号が異なり,Aが+1である2行ではBとCの符号が一致している。
 つまりこういうことだ。仮に,実験の結果Aが+1である2製品に対する好意度の集まりが,Aが-1である2製品に対する好意度の集まりよりも高かったとしよう。それは「Aによって好意度が変わる」ことを意味しているのかもしれないし,「BとCの組み合わせによって好意度が変わる」ことを意味しているのかもしれない。つまり,Aの主効果を意味しているのかもしれないし,BとCの交互作用効果を意味しているのかもしれない。
 このように,分解能IIIの直交計画では,Aの主効果とBCの交互作用効果を区別して捉えることができない。同じことが,Bの主効果とACの交互作用効果の間にも,Cの主効果とABの交互作用効果の間にもいえる。
 さきほどのシミュレーションでは,データに主効果モデルと呼ばれる分散分析モデルを適用した。つまり,主効果のみ推定するモデルである。本来,交互作用効果があるかどうかわからない状況下では,交互作用効果に関心があろうがなかろうが,交互作用効果の推定も行うべきである。交互作用効果が存在するのに,それを無視したモデルを使ってしまうと,誤った知見が得られてしまうかもしれないからだ。
 しかし,分解能IIIの直交計画で得たデータでは,そもそも主効果と交互作用効果を切り分けられない。分散分析では主効果モデルしか使えない。得られた主効果の推定値がほんとに主効果を表しているのかどうかは,分析者の解釈に任されている。3要因の主効果をたった4製品で精度よく推定できることの,これが代償である。

 いや,そういう言い方は後向きですね。もっとポジティブに言い直せば,こうだ。いま,測定値に影響しそうな要因がたくさんあるんだけど,どれが影響しているのかわからない,としよう。そんなときは,思いつく要因を片っ端から組み込んだ実験をやって,影響がありそうなものを選び出すことができたら便利である。そんな場合は,とりあえず交互作用効果は無視して,主効果だけに注目するという作戦が有用であろう。分解能 III の計画は,そんな場面で役に立つ。

最適計画への道

 さて,上のアイデアは,完全実施要因計画で必要とされる膨大な数の製品のなかから,釣り合いと直交性を保ったまま現実的な数の製品をうまいこと選び出すというアイデアであった。
 ところが実際には,要因の数や水準の数が多くなると,完全実施要因計画の製品数は膨大になる。そこから少数の製品を抜き出すと,釣り合いと直交性がどうしても崩れてしまうことが多い。どうしたらよいか?

 さきほど試したような簡単なシミュレーション実験を行い,推定精度が相対的にましな計画を選べばよい,というのがひとつの答えである。しかし,ありがたいことに,分散分析のような線形モデルがさまざまな実験計画のもとで持つ推定精度は,数理的に単純な一定の性質を持つことが分かっている(パラメータの推定量の共分散行列は情報行列の逆行列に比例する)。それによれば,推定精度がもっともよいのは完全実施要因計画ないし直交計画であり,釣り合いが崩れたり直交でなくなったりするにつれて推定精度が下がっていく。その低下の程度も簡単に計算できる。
 そこで,直交計画が作れない場合にも,実現可能な範囲内でもっとも推定精度のよい実験計画を探そう,という考え方が登場する。こういうアプローチを最適計画 optimum design という。完全実施要因計画のときに100%となるような推定精度の指標を使い,なるべく100%に近い一部実施実験計画を探そうというわけである。大変な作業だが,コンピュータが一瞬にして良い答えを出してくれる。推定精度を測る指標にはいくつかあるが,D-最適性という指標が用いられることが多い(この指標は,情報行列の逆行列の固有値の幾何平均の小ささを表す)。

 面白いことに,釣り合いの崩れ方のちがいのせいで,要因が直交している計画よりも直交していない計画のほうがD-最適性が高くなったりすることがある。そのせいで,実験計画法の教科書の説明どおりに丁寧に計画するより,専用のソフトを使って力づくで作ったほうが,良い計画ができることがある。
 市場調査の教科書をみると,実験計画の初歩的概念は紹介してあるし,直交計画について説明している親切な教科書もある。でも,ソフトウェアでの実験計画作成や,D-最適性について説明している本は見かけない。他の分野では常識となっている話なのに,なんでですかね? もっと普及してよい考え方だと思う。

最適計画の使い道

 もちろん,実験計画の良し悪しは,あるひとつの基準だけで決められるものでもない。たとえば,交互作用効果を無視した実験ではあるものの,あとでこっそり交互作用の有無もチェックしたい,といったこみいった裏の事情がある場合も珍しくないだろう。D-最適性が多少低くなってもいいから,釣り合いが崩れているのだけは困る,とか。どうしても実現が難しい製品があるとか。どうしてもこの製品を入れたいとか。現実は常に多様である。
 それに,ここまでの話はすべて,データのサイズが一定で,用いる製品数だけが変えられる,という場合の話である。実際には,製品数が増えればデータサイズも増える場合が多いだろう(各対象者に全製品を提示する場合とか)。その場合は,もちろんデータサイズが大きくなった方が推定精度は向上する。しかし,用いる製品数の増加は,実験にかかる時間やコストの増大を意味したり,対象者の疲労や回答の質の低下につながるかもしれない。いろいろな事情が複雑に絡み合い,もはやD-最適性どころの騒ぎではない。実験の設計はどうしたって,複数の条件を睨みながら考える複雑な意思決定プロセスなのだ。

 というわけで,一部実施要因計画をつくる際のお勧めはこうだ。リサーチャーが残業して,手作業で必殺の実験計画をひとつだけ作り込むのは,もうやめよう。さっさと実験計画用のソフトを買い,さまざまな設定の下の実験について、実験計画のいくつかの候補をつくってしまえばよい。何個でも一瞬にして作れるし、追加料金はかからない。で,それぞれの候補について,D-最適性や,その場で求められているさまざまな特性についてチェックし,一番ましな奴を選ぶ。という手順が,よろしいのではないかと思います。

Postscript

 ずいぶん前に,当時の同僚からもらった質問に応じて書きかけていた返事が,途中で放ってあったのをみつけた。捨ててしまうのもなんなので,具体的な業務と関わる部分を削って,ブログに書き留めて供養とする次第である。これを読んだ誰かが,私の愚かな誤りを見つけて馬鹿にするかもしれないが,それはまあ,そのときである。
 要因計画を用いた製品テストなんて市場調査の世界では稀だ,とおっしゃる方がいらっしゃるかもしれない。そのような,良く言えば現実力旺盛な方(悪く言えば想像力を欠いた方)のためには,たとえばBIBDによる製品テスト,応答曲面モデルによる最適化,そして離散選択型コンジョイント分析といった,市場調査においてもっとおなじみの場面においても,最適計画というアイデアがきわめてパワフルであることを示せばよいのだが,ちょっと気力が尽きました。
 本稿の種本はKuhfeld, Tobias & Garrett (1994, J. Mktg. Res.),岩崎(2006)「実験計画法」,などだが,どれもきちんと読んで理解しているとは言い難いので,すべての誤りの責任はもちろん私にあります。

実験計画の最適性について (お問い合わせに答えて)

2012年3月29日 (木)

 被験者に $t$ 個の刺激のうち 2個を呈示し、どちらが好きかを比較させる場合について考える。提示順序を考慮すると刺激対の数は$t(t-1)$個ある。対象者をいずれかの刺激対にランダムに割り当て、たとえば、先に見せたほうの刺激が非常に良いときに3, 同じときに0, 後にみせたほうが非常に良いときには-3... というふうに答えてもらうことにする。各刺激に対する選好の程度に,刺激間でどのような差があるか。

 刺激対$(i, j)$について対象者$k(=1,...,n)$が評価した際の反応を $x_{ijk}$ とし、
 $x_{ijk} = (\alpha_i - \alpha_j) + \gamma_{ij} + \delta_{ij} + \epsilon_{ijk}$
と考える。$\gamma$は組み合わせの効果($\gamma_{ij} = -\gamma_{ji}$)、$\delta$は順序効果($\delta_{ij} = \delta_{ji}$)、$\epsilon$は誤差項。$\alpha$の総和、$\gamma$の$i$ないし$j$を通した総和は0とする。
 この構造モデルに基づく分散分析は難しくない。全平方和は
 $S_T = \sum_i \sum_j \sum_k x_{ijk}^2$
自由度は$nt(t-1)$, 値の個数そのものである。要因は4つあるが、横着して主効果と誤差のことだけ考える。
 先に $\alpha_i$ の推定量について考えておく。これを $a_i$ と書くことにする。また、構造モデルの右辺から誤差項を取ったやつを$\mu_{ij}$、さらに$\delta$も取ったやつを$\pi_{ij}$、それぞれの推定量を $u_{ij}$, $p_{ij}$と書くことにする。順にゆっくり考えていくと、
 $u_{ij}= 1/n \sum_k x_{ijk}$
 $p_{ij} = 1/2( u_{ij} - u_{ji})$
 $a_i = 1/t \sum_j p_{ij}$
 以上を整理すると、$a_i = 1/(2tn) (\sum_j \sum_k x_{ijk} - \sum_j \sum_k x_{jik}) $となる。
 $\alpha_i$ の推定量が手に入ったところで、主効果の平方和について考えると、
 $S_\alpha = n \sum_i \sum_k (a_i - a_j)^2 = 2nt \sum_i a_i + 2n \sum_i \sum_j (a_i a_j)$
よくよくみると第二項は0である。スバラシイ。自由度は$t-1$。
 誤差の平方和は
 $S_\mu = n \sum_i \sum_j u_{ij}^2$
で、$S_E = S_T - S_\mu$ と考えれば楽勝である。自由度は $t(t-1)(n-1)$。このように、実に綺麗に算出できる。刺激間比較のためにはスチューデント化された範囲を使えば良い。

 では、まったく同じデザインで、提示順序を無視して分析したらどうなるか? 構造モデルは
 $x_{ijk} = (\alpha_i - \alpha_j) + \gamma_{ij} + \epsilon_{ijk}$
セル(i,j)とセル(j,i)をコミにして分析することになる。
 $u_{ij}= 1/(2n) (\sum_k x_{ijk} - \sum_k x_{jik})$
全平方和と主効果はさっきと同じで、誤差が変わる。自由度は $t(t-1)(2n-1)/2$となる。なぜそうなのか真正面から考えていたらなんだかわけがわかんなくなってきたんだけど、搦め手から考えると$\gamma$の自由度は$(t-1)(t-2)/2$だから、これでつじつまが合うのは確かだ。

 実のところ、前者のモデルはいわゆる「Scheffeの一対比較法」、後者のモデルは62年に芳賀敏郎が発表した「芳賀の変法」なのである。以上、佐藤信「統計的官能検査法」(日科技連)より。

 学部生の頃だったか、修士1年の頃だったか、当時心理統計を担当していた先生が、本来教えるべき内容をそっちのけにして、サーストンのケースIIIだのVだの、シェッフェ法だの誰々の変法だの、官能検査関連の超・超・超古典的な手法について蛇のように執念深く語り倒し、心底辟易したことがあった。それは呆れるほどにかびくさく,死ぬほど面倒で、何の役に立つのかさっぱりわからない議論であった。そんなのどうでもいいじゃん、もっと普通の講義をしてくださいよ、と思ったものである。あれは冬だったのだろうか、窓の外には枯れた芝生がみえた。
 このたび勤務先で、たまたま一対比較データの分析の話になり、いやその種のデータの分析方法は古くに確立してて、サーストン法とかシェッフェ法とか... と電話口で説明し始めたところで、突然あの時の記憶が、教室の空気の匂いさえも鮮やかに甦り、クラクラとめまいがするような感覚に襲われた。それはほんの数秒のことで、すぐにささやかで平凡な中年男ライフに戻ることができたけれども。
 電話を切ってから、頬杖をついて、手元にあった本をぼんやりめくっていたら、シェッフェ法とその変法についてわかりやすく説明していて、なあんだ、あのややこしい話って、こんなシンプルな線形モデルだったのか... と、気の抜けるような思いであった。先生、あの話、もう少しわかりやすく説明できたかもしれませんよ。そして私も、もっとまじめに勉強するか、あるいはもう少し別のことをしておいたほうがよかったかもしれません。
 とはいえ、いまとなってはなにもかも詮無い話である。あの先生が亡くなってからずいぶん経つ。

突然に、一対比較法について

2009年9月25日 (金)

 片側検定が許されるのはどんな状況だろうか?
 こういうことを延々と考えても,なんら益がないことははっきりしているのだが,それでもいったん考え始めると,途中で投げ出すのは難しい。いや,むしろどうでもいいことだからこそ,ここまで真剣にあれこれ考えるのかもしれない。展望なき人生において,真に考えるべき事柄はおしなべて深刻な事柄であり,そして深刻な事柄を考えるだけで頭のなかにブザー音が響く。ストップ! 考えるな!
 こういうのをなんていうんだっけ? 現実に対して真剣に向き合うだけの力がないことを。「駄目な人」? もっとひどい言い回しもたくさんありそうだ。その辺について熟考するのは,次回の人生にとっておいて,当座の問題は。。。片側検定が許されるのはどんな状況か。

 仮説検定は統計学やデータ解析法の初級コースに登場する基礎的概念だ。片側検定は,仮説検定の説明のなかで登場する考え方で,これまた初歩的内容にはちがいない。
 ところが,こういう話題の困ったところは,基本的な疑問であるにもかかわらず,統計学の教科書に答えが書いてあるとは限らない,という点である。たとえば,手元にある本のなかから中村隆英ほか「統計入門」(1984)をめくってみると,この本でひとり勉強したころの手垢や落書きが目に付いて,もう涙が出そうである。20年近い年月が経ってしまった。いや,そういう感傷は置いておいて。。。

「棄却域が x ≦ a または x ≧ b ( a < b ) の範囲となる検定方式を両側検定といい,[...] 棄却域が x ≧ c といった範囲になる検定方式を片側検定という。」(p.211)

 これだけである。この直後から,説明は所与の対立仮説の下での棄却域設定についての議論へと移っていく(帰無仮説μ=c, 対立仮説μ≠c に対する一様最強力検定は存在しない,とか)。しかし,棄却域が片側になるような対立仮説(たとえば μ>0)を設定してよいのはいったいどんなときなのか,という疑問には,この本は答えてくれない。それこそが,俺にとっての疑問であり,多くの人にとっての疑問であるはずなのだ。
 この「統計入門」は初学者向けの統計学の教科書で,決して数理統計学の専門書ではない。大変わかりやすい,良い本だと思う。それでも,ユーザの肝心の疑問にはなかなか答えてくれない。教科書とは往々にしてそういうものである。一冊の教科書に頼ってはいけない。

 「片側検定が許されるのはどんなとき?」という疑問に対し,可能な答えが3つあるように思う。例として,母平均μと定数cとを比較する検定について考えよう。H1:μ>cという対立仮説を設定し,片側検定を行ってよいのはどんなときか?

 もっと具体的な事例に当てはめて書き直したいのだが,これが案外難しい。平均と定数を比較する検定など,実際にはなかなか用いられないからだ。そこで,俺がかつて考えた素晴らしい事例を紹介したい。以前統計学の講義を担当していたときに考案した名作である。これが人々に知られないまま消えていくのは,あまりにもったいない。

町工場経営の山田さんはパチンコ玉を作っています。愛用のパチンコ玉製造機は,もう何十年も休みなく動き続け,新品のパチンコ玉を吐き出し続けています。パチンコ玉一個の重さの平均はぴったり2g, 寸分の狂いもありません。年月とともにスイッチやツマミの文字は薄れてしまいましたが,あまりに安定的な機械なので,オーバーホールの必要もツマミをまわす必要もなく,山田さんはすっかり安心していました。
 ところがある日,山田さんの愛猫が,パチンコ玉の重さを変えるツマミの上に飛び乗ってしまいました。そのツマミを右に回すと,パチンコ玉は少しだけ重くなり,左に回すと少しだけ軽くなってしまうのです。大変だ!山田さんはあわててツマミを調べましたが,目盛がすっかり消えてしまっているので,ネコがツマミをまわしたのかどうか,まったく見当がつきません。
 そこで山田さんは,ネコが飛び乗ったあとで生産されたパチンコ玉からN個を抜き出し,その重さを測定器で調べることにしました。
 山田さんは次のように考えました。これから生産されるパチンコ玉の重さの集合を母集団と考えよう。私はこれから,無作為抽出したサイズNの標本を手に入れるわけだ。母平均が2gと異なっているといえるか,検定によって調べてみよう。

 。。。おかしいな,名作だと思ったのに。こうして書いてみると,俺の資質と能力,人としての常識,といったあたりに深刻な疑念を感じざるを得ない。
 まあいいや,この例で話を進めると,山田さんが片側検定を行って良い場合について,以下の3つの答えかたがある。

。。。書いていてだんだん頭が痛くなってきた。ネコの神ってなんだよ。
 この3つの説明の違いは,いっけん言い回しの差のようにみえるかもしれないが,よく考えてみるとかなり異なる示唆を持っている。そのことは,山田さんが得た標本平均が,たとえば1.95gだった場合を考えればよくわかる。

 この問題について述べている解説書を探してみると,意外なほどに意見が割れている。手元にある日本語の解説書に限定して書き抜いてみよう。
 まず,(A) μ<cではないと事前に知っているときしか片側検定を認めない解説書。

吉田「本当にわかりやすいすごく大切なことが書いてあるごく初歩の統計の本」(1998), p.164

研究者がある方向の差を研究仮説として設定したというだけでは片側検定を行なうべきではなく,研究仮説とは逆の方向の差がその領域や文脈においてまったく意味をもたないと考えられるような特殊な場合でない限り,基本的には両側検定を行うべきだと思います

市原「バイオサイエンスの統計学」(1990), p.194

あらかじめデータがどちらの方向に偏るかを特定できる場合[...]は,片側の効果だけを調べればよいという理屈になる。しかしそのような場合ですら,例外的に逆の効果をもつこともありえるので,やはり,両方の可能性を考えて両側検定しておくのが妥当と考えられる

実のところ,(A)における片側検定に反対する人はいないのではないかと思う。もしμ<cではないことをあらかじめ知っているのであれば,仮に標本平均がcをはるかに下回ったとしても,帰無仮説μ=cを棄却する理由がない。事前知識を生かして棄却域を設定すること自体には,文句のつけようがない。意見が割れるのは,(A)以外の状況での片側検定を認めるか,という点である。上の2冊は,片側検定が許される状況を(A)に限定している解説書であるといえよう。
 (A)の状況は現実にはほぼあり得ないから(ふつうネコは喋らない),この立場に立てば,片側検定はめったに適用できないことになる。

 次に,(B) μ>cと予測しているときに片側検定を許容する解説書。

山内「心理・教育のための統計法」(1987) p.104

[片側検定の]このような方向性を決定するのは,いままでに知られている情報や,理論的な仮定である。ただし,実験が行われる前に,方向性に関してH1や棄却域は決定しておかなければならない。[...] 例題:心理学者は,ある社会性のテストにおいて,テストの誤答数が通常の小学4年生において平均値25個,標準偏差4.3個であることがわかっているとしよう。いま,リーダーシップの高い子どもは,そのテストでより少ない誤答数を示すであろうと想定した(したがって,H1は方向性のある片側検定を要請する)

いま手元にないのだけれど,名高い古典である岩原信九郎「教育と心理のための推計学」も,たしかこのような説明をしていたと思う。
 うーむ。ど素人の俺がこんなことを書くのは身の程知らずもよいところだが,こういう説明は,そのー,ちょっとまずいんじゃないでしょうか,と思う次第である。
 たいていの分析者は,なんらかの予測なり想定なりを事前に胸に抱いているものだ。リーダーシップの高い子どもは誤答が少ないんじゃないかとか,ネコはツマミを右に回したんじゃないかとか。上記の解説書に従えば,とにかく想定がはっきりしていれば片側検定が使用できる,ということになろう(両側検定の出番などほとんど無くなってしまうのではなかろうか)。さて,健全な実証研究においては,蓋を開けてみると,事前の想定は往々にして裏切られるものである。実のところ,リーダーシップの高い子どものほうが誤答が多い,ということがあきらかになるかもしれない。山田さんの予測とは逆に,パチンコ玉の重さの標本平均は1.8gかもしれない。データに裏切られた分析者は,そのデータを投げ出すべきだろうか? それはデータに失礼というものだ。むしろ,自らの理論を修正するために,そのデータを再分析するのがあるべき姿であろう。というわけで,今度は新たな対立仮説についての検定が有用になるだろう。その新たな対立仮説が,たとえばμ<cならば,片側検定を2回繰り返すことになるわけだ。2回の検定をあわせたType I Errorは有意水準の2倍になってしまう。だったら最初から両側検定でType I Errorを制御しておいたほうが良いのではないでしょうか?
 と,どきどきしながら偉そうなことを書いているわけだが,だんだん気弱になってきたので援軍を呼ぼう。世の中には物好きな人がいるもので,Lombardi&Hurlbert(2009)は統計学の教科書を52冊集め,片側検定が許される状況についての説明を抜き出し,Reasonable(白), Vague(灰色), Bad(黒)の3つに分類している。俺の分類と照らし合わせると,どうやら(A)と(C)は白,(B)は黒とされているようだ。ただし,この著者らはなかなか厳しくて,俺の目から見て(C)タイプの説明であっても,分析者が対立仮説を片側に設定する理由がassumptionとかpredictionなどと表現されているだけで,その教科書はもう黒呼ばわりである。ひええ。

 さて,いよいよ本命の登場である。(C) μ=c か μ<c かというちがいに関心がないときに片側検定を許容する解説書。

南風原「心理統計学の基礎―統合的理解のために」(2002)

いま,帰無仮説をH0:ρ≦0とし,対立仮説をH1:ρ>0と設定したとします。つまり,『母集団相関はゼロまたは負である』という帰無仮説を『母集団相関は正である』という対立仮説に対して検定するということです。このような検定は,たとえばある市において,学級の児童数と,学級のなかで授業がわからないという不満を持つ児童の割合との相関係数について,帰無仮説が棄却されて『母集団相関は正である』という対立仮説が採択されたら,学級定員を減らすための検討を開始する,というような場合などに考えられる可能性のあるものです。このとき,r=-.5のように,負で絶対値が比較的大きな r が得られたとしても,これは明らかに対立仮説を支持するものではありませんから,こうした負の領域には棄却域を設けず,正で値の大きい r の範囲のみを棄却域とする片側検定を用いるのが合理的です。このように,片側検定が有効であるケースは考えられますが,通常の研究において検定が用いられる文脈では,このようなケースはあまりありません

永田「統計的方法のしくみ―正しく理解するための30の急所」(1996), p.107

たとえば,何らかの新しい機械の購入を考えているとしよう。その機械のある特性値の母平均μが"基準値"よりも上回るときのみ購入の価値があると考えるのならば[...](右)片側検定をおこなうことが妥当である

吉田「直感的統計学」(2006), p.245

消費者グループが,A自動車会社の『高速道路では一リットルあたり30キロ走る』という新車に関する主張をテストしている。[...] 消費者は一リットル30キロ走るのかそれ以下なのかに興味がある。もしもA自動車会社の主張に対して,サンプル平均が著しく低ければ,μ=30は到底受け入れることはできない。H0を棄却するということは会社を訴えるといったアクションを取ることになるかもしれない。それに対して,もしもサンプル平均が 30キロ以上だったら,つまり燃費が自動車会社が主張しているより良いときに消費者はアクションを取るだろうか?多分彼らはただH0を受け入れ,黙っているだろう。そういう意味でμ=30キロでもμ=35キロでも同じ受容域に入るのである

俺の手元にある本の中では,このタイプの解説が一番多い。英語なので抜き書きは省略するが,古典的名著として知られるWiner "Statistical Principles In Experimental Design" (2nd ed.)も,一番頼りにしている Kirk "Experimental Design" (3rd ed.) も,このタイプの説明であった。
 告白すると,(C)の状況下での片側検定は間違っている,と俺はこれまで固く固く信じてきた。片側検定について他人に説明する際にも,それは山田さんの関心の持ちようで決まる問題じゃない,ネコがツマミを左に回していないと山田さんが事前に知っているかどうかという,事前知識の問題なんですよ,と強調してきた。時々見かける(B)(C)タイプの説明は,筆が滑ったか,ないし誤りだと考えていたのである。
 このたび,仕事でこの関連の話を考える用事があって,南風原(2002)を読み直し,この名著がなんと(C)の状況での片側検定を認めていることに気づいて,うわあ,また俺は嘘をついていたか,と青ざめた。それがこの文章を書くきっかけである。

 (C)の状況下での片側検定は誤りだ,と俺が考えていたのはなぜか。たぶん,吉田(1998)や市原(1990)のようなタイプの説明をどこかの本で読んだせいだろう。それがなんという本だったのか,思い出せないのがつらいところだ。
 でも,思い返してみれば,自分なりに納得した理由が,大きく二つあったように思う。テクニカルな理由プラクティカルな理由である。

 まずテクニカルな理由。(C)の状況での片側検定を許容すると,物事はもう泣きたいくらいにややこしくなってしまう。
 両側検定の場合について考えよう。帰無仮説をμ=c, 対立仮説をμ≠cとする。帰無仮説の下での平均の標本分布を考え,その両側5%ぶんを塗りつぶす。この5%を有意水準と呼ぶ。有意水準はType I Errorの確率,すなわち「帰無仮説が真のときそれを棄却してしまう確率」に等しい。
 では,片側検定の場合はどうなるか。対立仮説をμ>cとしよう。素直に考えると,(C)の場合,対立仮説が真でないときはμ=cかもしれないしμ<cかもしれないのだから,帰無仮説はμ≦cだ。南風原(2002)もそう述べている。しかし,南風原先生は説明を端折っておられるが,帰無仮説がμ≦cであっても,棄却域を決める際にはμ=cの下での平均の標本分布を用いざるをえない。μ≦cの下での標本分布は一意に決まらないからだ。すると,この分布の右側5%ぶんを塗りつぶしたとして,この5%の意味,すなわち有意水準の意味が変わってくる。それはもはやType I Errorの確率ではない。「Type I Errorの最大確率」とでも呼ぶべきものになってしまうのである。
 これに対し,永田(1996)がそうしているように,片側検定においても帰無仮説はμ=cだ,と言い張る手もある。(A)における片側検定と同じく,μ<cを考慮の外に置いてしまうわけである。クリアな解決策だが,μ<cはどこにいったんだ,という疑問が浮かぶ。
 南風原流の説明と永田流の説明ではどちらが一般的だろうか。また物好きな人がいるもので,Liu&Stone(1999)は44冊の教科書を調べ,うち24冊はH0:μ≦c(南風原路線),20冊はH0:μ=c(永田路線)であると報告している。彼らに言わせると,どっちの説明でもかまわないけど,前者の路線を取る場合には有意水準の意味が変わることを説明しなければならない。しかし,永田路線の20冊のうち11冊がその点をちゃんと説明しているのに対し,南風原路線を取る24冊のうち実に20冊はその説明を端折っている由。これは統計教育における由々しき問題点だ,と彼らは憂えているのだが,うーん,教科書がややこしい話を端折ったからと言って,一概に責めることはできないと思うんだけど。
 それはともかく,このややこしさこそが,(C)の状況下での片側検定が許されない理由なのだ,と俺は考えていた。H0:μ≦0ならば(南風原路線),仮説検定においてもっとも基礎的な概念である有意水準の意味が,両側か片側かというちょっとした選択によって変わってしまうことになる。H0:μ=0とすると(永田路線),直感に反する帰無仮説を設定することになる。
 いまこうして書いてみると,俺が感じていたのは,(C)の状況下での片側検定はいろいろとややこしい問題を抱えている,ということに過ぎないようだ。片側検定が誤りだという根拠とは言い難い。反省。

 というわけで,(C)における片側検定は間違いだ,という俺の確信は哀れにも揺らぎ始めているのだが,完全な納得にまでは至っていない。そのように未練がましく考える背後にあるのは,かつて(C)が誤りだと思った第二の理由,プラクティカルなほうの理由である。
 素朴な言い方で恥ずかしいが,「μ=c か μ<c かというちがいに関心がない」だなんて,そんなの嘘じゃないか,と思うのである。
 これは具体的な例で考えた方がいい話だと思う。二つ例を挙げてみたい。まず非劣性試験の話。聞きじったところによると,ジェネリック医薬品の承認申請のためには,新薬の場合とは異なり,先行薬よりも(ある限界を超えて)劣っていない,ということを示すための臨床試験を行うのだそうだ。で,テクニカルにはいろいろな工夫があるにせよ,基本的発想としては片側検定を用いるのだそうである。なるほど,既存の薬との間に差がなかろうが,ジェネリック薬のほうが優れていようが,それはどうでもいいわけだから,まさに(C)の状況であるといえよう。
 今度は市場調査の話。日本の市場調査業界の基礎を築いた有名な実務家が書いたハンドブックから引用する。製品テストを例に,(C)的な状況を端的に説明している。

後藤「市場調査マニュアル」(1997)

テストの目的が製品A, Bのどちらを選ぶか,または選択を保留するかという場合には両側検定となる。[...] 製品Aの評価が高ければAを選び,そうでなければBを選ぶという場合には[...]片側検定となる。[...たとえば] 改良品Aがよければ現行品Bを改良品に切り替えるが,そうでなければ現行品Bを変えない,あるいは競合銘柄Bよりよければ当社製品Aを発売するが,そうでなければ発売しない場合。

新製品が対照製品と同等であろうが,新製品のほうが劣っていようが,それはどうでもいい,だから片側検定だ,というわけである。理屈としてはさきほどの臨床試験の話と同じなのだが,こんどは違和感を感じる。これは本当に(C)の状況なのか,と疑問に思われてならない。
 なぜなら,もし新製品Aの評価が対照製品Bよりも低かったら,単に「発売しない」と決定するだけでは済まないだろう,と思うからだ。Aは,Bよりも劣った製品になってしまったのだろうか,それとも差がないだけなのだろうか。観察された差が統計的にみて意味を持つものなのかどうか,調べてみたくなりませんか? 再び検定してみたくなりませんか? 新製品が対照製品と同等であろうが劣っていようがどうでもいいなんて,嘘じゃありませんか?
 このふたつの話はどこがちがうのか。まず,統計的推測と意思決定の間の距離がちがうと思う。臨床試験のほうは,その結果がそのまま製造販売の承認という社会的決定に直結する。いっぽう,消費者調査の結果をそのまま新製品上市の決定に直結させる発想は,あまり現実的とは思えない。むしろ,ビジネス上の意思決定は(科学的推論がそうであるように)幅広い情報に基づいてなされるべきものであり,消費者調査はそのデータソースの一つに過ぎない,と考えたほうがよいだろう。
 別の観点からいうと,このふたつの話は,分析主体の「関心」の多様性が異なる。前掲のLombardi&Hurlbert(2009)は面白いことをいっている。現代の仮説検定論を築いたNeyman(1937)は,片側検定が許される状況について"[where] we are interested only [...] in one limit"と表現しているそうなのだが(Cタイプですね),このinterestはただのinterestではなく,個々の分析者から独立な"collective interest"として解するべきものである由。なるほど,上手いことをいうものだ。統計的仮説の設定において問われるのは,個々の分析者の関心ではなく,データにアクセス可能な当事者たちの間のコンセンサスなのだ。臨床試験の場合は,「ジェネリック薬が先行薬と同等であろうが,ジェネリック薬のほうが優れていようが,そのちがいはどうでもいい」という点に関係者すべてが合意するだろう。いっぽう製品テストのほうでは,この合意が得られるかどうか怪しい。調査のステイク・ホルダーは多種多様な関心を持っているはずだ。経営陣は「新製品が対照製品と同等であろうが,新製品のほうが劣っていようが,そのちがいはどうでもいい」と思っているのに対して,R&Dは大いに関心がある,という風に。
 こうして整理してみると,俺が(C)タイプの説明に抵抗を感じるのは,μ<cとμ=cのちがいに「全く関心がない」場合の片側検定が論理的に誤っていると思うからではない。単に,そんな状況は極めて限定的だ,と感じているに過ぎない。いいかえると,片側検定が許される状況を,その場その場の「関心」をキーワードにして定義するのは甘すぎる,と危惧しているだけである。これもまた,(C)タイプの説明を誤りとみなすだけの根拠とはいえない。反省。

 片側検定が許される状況をどのように定義するかという問題について,勝手に悩み,勝手に論破されて勝手に反省しているわけだが,そんな定義にかまけていないで,(D)そもそも片側検定をやめちゃえばいいじゃん,という意見もあるだろう。日本語でみかけた例としては,佐伯・松原編「実践としての統計学」(2000)のなかで,佐藤俊樹がそういう意見を書いている。この先生は「不平等社会日本」で知られる気鋭の社会学者だが,こういうテクニカルな話題にも造詣が深いんですね。
 佐藤によれば,(C)の状況における片側検定には論理的飛躍がある。

帰無仮説θ=0の否定はθ≠0であって,θ>0 (やθ<0)ではない。にもかかわらず,何の説明もなしに『片側検定では帰無仮説θ=0を棄却することでθ>0(あるいはθ<0)という対立仮説がとられる』といえば,それはまさに論理的飛躍である

しかし佐藤は(A)の立場も拒否する。

論理を一貫させようとすれば,いくつかの本でやっているように,片側検定は何らかの根拠でθ<0(あるいはθ>0)でないと確定できる(無視できる)場合に使われる,とするしかない。[...]こういう形にすれば明確だが,[...] 限られたケースでしか片側検定は使えなくなる

 佐藤の議論の眼目はこうだ。この論理的飛躍は片側検定に限ったことではない。両側検定だって同じことである。なぜなら,分析者は両側検定によって帰無仮説μ=cを棄却するが,それを根拠として彼が主張するのは,たいていμ≠cではなく,μ>cなりμ<cなりであるからである。検定ユーザはどのみち論理的に飛躍しているのである。

この飛躍は統計学的論理の一貫性と有用性という二つのメタ手順を同時に考えることでしか対処できない。θ<0 (あるいはθ>0)でないと確定している場合や,『差がある』とだけいえればいい場合には単純である。そうでない場合が問題になる。これについてはいくつかの解があろう。一つの解は,どちらにしても飛躍が発生するのだから,両側検定と片側検定は無差別である,とする。要するに,どちらを使ってもよい。例えば『片側で有意水準5% (両側なら10%)』と明記して,有意水準の実質的値について誤解が起きないようにすればよい。もう一つの解は,両側検定を使った上で,統計学的には根拠のない推論を一部していると認めるほうがよい,とする。結論だけいえば統計学的論理の一貫性を追求するのと同じだが,理由がちがう。飛躍を統計学に押し付けて消去するより,研究の内部にリスクとして明示的におくべきだ,という判断による。私としては第2の解がよいと思う

二つの解が実質的にどう異なるのか,俺はちょっと理解できていないんだけど,魅力的な見解だと思う。こうしてみると,片側検定が許される状況の定義なんて,どうでもいいような気がしてくる。
 しかし,個々の分析者が佐藤に従い,熟慮した上で片側検定を使うのを止めたとしても,他人が片側検定を乱用するのを止めることはできない。誰もが仮説検定論について熟慮するわけではないのだ(仮にそんなことが起きたら,それこそ労力の無駄遣いだと思う)。さらに,片側検定が乱用される動機は十分にある。分析者はなんとか有意差を得たいと思っていることが多い。想定する方向での棄却域を広げてくれる片側検定は,分析者にとっての甘い誘惑なのである。実は私はμ>cだけに「関心」があるのです!と勝手に宣言し,片側検定を乱用する人々の姿が目に浮かぶ。
 片側検定を使うのを止めよう,という呼びかけだけでは不十分だ。片側検定が許容されるのはどんな状況か,ユーザのための操作的ガイドラインが必要だと思う。

 というわけで,あれこれ読みかじってみた結果,論理的には(C)タイプの説明が正しいと,半ば納得するに至った。その反面,(C)が示唆するcollective interestという基準はなかなか理解されにくいだろうなあ,という危惧も捨てきれない。南風原先生がそうしているように,片側検定の乱用を戒める文言を付け加えるのも大事だが,我々統計ユーザに対するガイドラインとしては,いっそ(A)タイプの説明にまで後退しちゃったほうがいいんじゃないか,と未練がましく考える次第である。

片側検定の迷路

2008年5月 3日 (土)

 2要因以上のunbalancedなANOVAでは,ある水準についてのunweightedな平均(SASのproc glmでいうmeans)と重みつき平均(lsmeans)の値が異なる。事後の多重比較で用いるべきはどちらだろう?
 森・吉田の本は前者を比較する方法を紹介しているが,それは前者のほうが手計算しやすいからだと思う。この本の基になっているKirkの本は,不揃いになった理由が処理水準と無関係な欠損である場合は前者でよいとか,その不揃いさが母集団の比率を反映しているのなら後者がよいとか述べていたと思うが(小牧「データ分析法要説」によればWinerもそういっているらしい),この2つのどちらでもないケースも少なくないと思う。
 そもそも,unweightedな平均を比較してよい理由がよくわからん。本来は常にlsmeansを使うべきなのではないか?

 ほんとはこんな古い話題に頭を悩ませていないで,もうちょっと新しい方法を考えるとか(混合モデルを使うとか),あるいはもっとオカネになることを勉強するとか,あるいはもっと楽しいことを考えるとか...いろいろあるだろうに。知りもしない分野で専門家づらしている悲劇。この問題を解決したところで誰が喜ぶわけでもないという悲喜劇。

 考えてみたら,こんな疑問はありふれたFAQのはずだ。どこかに解説論文があるにちがいない。悩んでいないでさっさと検索すべきだな。

覚え書き:今週の悩み事

2007年7月11日 (水)

集中講義で話の種に使おうと思うので,ここでメモしておく。

日経新聞

「日本は美しい」5割強、内閣府世論調査
 内閣府は5日、安倍晋三首相が目指す国家像「美しい国」づくりに関する特別世論調査の結果をまとめた。今の日本を「美しい」とした人は10.6%で「どちらかというと美しい」の42.7%と合わせ5割強に上った。「日本の美しさとは何か」を複数回答で尋ねたところ、1位は山や森などの「自然」で 80.0%。伝統工芸などの「匠(たくみ)の技」が58.5%、田園・里山といった「景観」が52.8%で続いた。[...]

毎日新聞

内閣府調査:過半数が「美しい国」 20、30代では逆転
 内閣府は5日、安倍晋三首相が提唱する「美しい国」づくりに関する世論調査の結果を発表した。現在の日本を「美しい」と回答した人は、「どちらかというと美しい」(42.7%)も含めると53.3%と過半数を占めた。「美しくない」と回答した人は、「どちらかというと」(31.8%)を含めて43%だった。
 ただ年齢別にみると、40代以上は「美しい」が過半数を占めたのに対し、20代の57.6%、30代の50.5%は「美しくない」と回答。若い世代は日本を美しく思っていない側面が浮き上がった。
 一方、日本の美しさを複数回答で尋ねたところ、山や森、海などの「自然」が80%、伝統工芸や宮大工の技術など「匠(たくみ)の技」が58.5%、田園・里山の風景など「景観」が52.8%--などが上位だった。内閣府は「上位の自然や匠の技などに若者はなかなか関心を持ちにくいので、『美しくない』が多くなったのではないか」と分析している。[...]

産経新聞

日本は「美しい」53.3% 内閣府世論調査.
 日本を美しいと考える人が53.3%、美しくないとの回答は43.0%であることが5日、内閣府が発表した「美しい国づくりに関する特別世論調査」で分かった。40~70歳代以上の過半数が「美しい」と答え、20~30歳代では過半数が「美しくない」とした。「重要と思う美しい国の姿」では「文化、伝統、自然、歴史を大切にする国」(47.5%)が最も多く、次いで「規律を知り、凛(りん)とした国」(22.6%)、「世界に信頼され、尊敬されるリーダーシップのある国」(16.7%)の順。[...]

時事通信

「日本は美しい」53%=「美しくない」43%-内閣府調査
 政府は5日、安倍晋三首相が政権構想の基盤に据える「美しい国づくり」に関する特別世論調査の結果を発表した。それによると、現在の日本について「美しい」と回答した人は53.3%だったのに対し、「美しくない」は43.0%に上った。若い世代で「美しくない」と答えた人が過半数を占めており、内閣府は「いじめや虐待などの報道に接し、現状に問題意識を持っている人が多い表れ」としている。
 調査は、日本の「美しさ」に対する国民のイメージを探り、国内外に向けた日本の魅力発信などの施策づくりに役立てるのが狙い。5月24日から6月3日にかけ、3000人の成人男女を対象に実施し、有効回答率は60.9%だった。

共同通信

「日本美しい」は1割強 内閣府が世論調査
 内閣府は5日、安倍晋三首相が唱える「美しい国づくり」に関する特別世論調査の結果を発表した。日本の現状を「美しい」とした回答は10・6%にとどまり、「どちらかというと美しい」を合わせてようやく53・3%と半数を超えた。
 「美しくない」は11・2%、「どちらかというと美しくない」は31・8%。40代以上は「美しい」「どちらかというと美しい」が過半数だったが、20-30代は逆に「美しくない」「どちらかというと美しくない」が過半数。内閣府は「特に若い世代が美しくないととらえている点をよく考慮し、今後の政策を考えたい」としている。
 「日本の美しさとは何か」という質問(複数回答可)では、山や海など「自然」が80・0%でトップ。このほか伝統工芸など「匠の技」が58・5%、田園や街並みなど景観が52・8%、能や歌舞伎など「伝統文化」が50・8%だった。敬語や方言など「日本語」は32・8%にとどまった。

朝日,読売ではみつからなかった。内閣府のリリースはここ

とにかく,誰かを馬鹿にするのだけはやめよう,と思う。それはなにも生まない。ただ自分の頭が良いような気がするだけ,ただ気分がよくなるだけだ。
なぜこんな調査が生まれるのか。それの結果はどのように流通し,どのように利用されていくか。そのことだけを考えたいと思う。

美しい国づくりに関する特別世論調査

2007年6月 4日 (月)

 先週だったか,新聞に子どもの生活についての調査研究が紹介されていた。いわく,小学5年生を対象にした調査で,起床時間と「学校が楽しいかどうか」との間に相関が見られた由。早起きの子どものほうが,学校を楽しいと感じている割合が高いのだそうだ。ふうん,と適当に読み流して,それきり忘れていた。
 ところが,数日前にぼんやりネットを眺めていたら,この調査をぼろくそにやりこめている人がいた。あれこれ見てみると,そういうブログがたくさんあって(一例),ちょっと驚いた。なんというか,世の中には正義感に満ちた人が多い。
 批判や罵倒の数々を眺めているうち,なにやら形容しがたい感情がこみ上げてくるのを感じて,自分でも意外であった。以来数日,このことをあれこれ考えていたのだが,これは誰が悪いという話ではなくて,構造的な問題なのではないか,と思い至った。

 生活サイクルの地域差を指摘するものから,学者はみんなアホばかりだ式の悪口まで,ブログにみられる議論のレベルは様々だが,批判の主旨を煎じ詰めると,(1)標本に無作為性がない,(2)因果関係はわからない(早起きのせいで学校が楽しくなるわけではない),という二点に尽きるように思う。
 (1)はまた別の話になるので横に置いておくとして,(2)については,なるほどその通りだ。起床時間と学校の楽しさに相関があるとしても,どちらが原因でどちらが結果かはわからないし,そもそもこの2つの変数は共になんらかの変数の影響下にあるのかもしれない(たとえば,東京圏よりも地方都市の子どものほうが早起きで,かつ学校を楽しく感じるのかもしれない)。ではこのツッコミは,この調査に対して(ないし報道に対して),果たして正当な批判となっているだろうか。
 たとえば朝日新聞の記事(05/29)はこうだ。

早起きの子どもは学校が好きで、楽しいとも感じている――。早起きと学校好きの間にそんな関係があることが、教育学や食物学の専門家でつくる「子どもの生活リズム向上のための調査研究会」の調査でわかった。

調査の概要が述べられ,最後に

「研究会代表をつとめる明石要一・千葉大教授は「起床が早く、バランスのいい朝食を食べている家庭は、生活にリズムがある。学校も家庭も、もっとこの問題に関心を持ってほしい」と話す。

早起きの「せいで」学校が好きになるとは,どこにも書かれていない。他社の記事では別の知見が紹介されているが(朝食が和食の子どものほうが早起きだ,云々),書き方は同様である。記事で述べられているのはあくまで相関的知見であって,因果的知見ではない。「××な人は○○の傾向がある」という記述を「××であることは○○をもたらす」と勝手に読み替えるのは,あまりに軽率というものだ。そんな読み方が許されてしまうようでは,相関の記述などおよそ不可能になってしまう。
 この調査の報告書を読んだわけではないけれども,生活スタイルと学校への態度の間に因果関係が主張されているとは思えない。横断調査から因果関係はわからない,というのは調査のイロハであり,研究者はもちろん新聞記者だって先刻ご承知である。鬼の首を取ったかのように悪口を書く前に,もう少し丁寧に文章を読んだほうがよいのではなかろうか?

 と,いうのが最初に感じたことなのだが,あれこれ考え合わせると,この人たちの言うことにもそれなりの説得力はあるのかもしれない,と思う。
 たとえば上記のように,早起きの「せいで」学校が好きになるわけじゃないだろう,と批判する人がいる。そんなことは記事のどこにも書かれていないわけで,あきらかに勇み足である。しかしもしこれらの報道を,多くの人が因果的説明として理解しているとしたらどうだろうか。現に,いまwebでこの調査について言及しているブログを検索してみたら,その半分は「くだらない調査だ」というツッコミだが,もう半分は「なるほど,規則正しい生活って大事だ」「やっぱり朝は和食にしようかしら」といった,驚くほど素直な感想なのである。こうしてみると,勇み足にも理由がないとはいえない。
 この学者たちは学校の問題を家庭に押しつけようとしてるんじゃないかしら,と嫌悪感を示す人もいる。もちろんこれは,批判というよりは勘ぐりである。しかし,研究会代表である教育学の先生は,中教審分科会委員を勤める,文教政策に関わりのある学者である。この調査もおそらくは,文科省が只今鋭意推進中の「早寝早起き朝ごはん」国民運動と関係があるのだろう。決して無理な勘ぐりとはいえない。
 
 というわけで,あちこちで目にする批判は,正当とはいえないけれども,社会的文脈を踏まえれば一定の意義を持っているのかもしれない,と思う。
 にもかかわらず --- ここからが本題なのだが --- ブログをあれこれ眺めるうちに,なんともいえない感情がこみ上げてくる。ありのままに言ってしまえば,強い不快感といらだちを感じる。いったいなぜだろう?

 上記調査についての各社の報道を,いまwebで眺めていて感じるのは,ああみんな上手く書くものだなあ,ということである。 
 日経の記事はこうだ:

 主食と主菜、副菜、一汁の4品がそろった朝食を食べている小学5年生の61.8%が「学校がとても楽しい」と感じていることが28日、千葉大の明石要一教授(教育社会学)ら研究者グループが2006年に実施した調査でわかった。

あくまで相関の記述に留まっている。実に正しい文章である。しかし同時にこの記事は,「朝食が学校への態度に影響する」という因果関係を読み手にうまく暗示している。この因果関係は,家庭生活が大事だ,という社会的通念と合致しており,それゆえにわかりやすい。こういうところが,ああ上手いなあ,と思う。
 書き手の立場で考えると,相関の記述だけではつまらないのである。人は常に説明を求める。学校が楽しいと感じる子どもの特徴について知りたい,などという奇特な人はいない。その子たちはなぜそう感じるのか,学校を楽しい場所にするためにはどうしたらよいのか。人は理解可能な原因,実行可能な方策を求める。だから書き手は,たとえそれが横断調査の結果に過ぎないとしても,つまり手元にあるのは相関的知見だけであるとしても,嘘にならない範囲内で,なんとかして因果的なニュアンスを醸し出そうとする。こうしてセクシーな報告ができあがる。
 だから,たとえばこういう文章を書いてはいけない。

 「学校がとても楽しい」と感じている小学5年生には,主食と主菜、副菜、一汁の4品がそろった朝食を食べている傾向があることがわかった。

同じく相関を記述しているのだが,このタイプの文章を読むと,多くの人が「わかりにくい」という(おかげで,俺は延々悩んだことがある)。正確には,わかりにくいのではない。ただセクシーさに欠けるのである。なぜなら,(社会的通念によって想定されるところの)結果側が先,原因側が後に書かれているせいで,因果的なニュアンスを読み取るのが難しくなっているからである。
 さらに高度な例を朝日の記事にみることができる。先生のコメントの素晴らしさといったらない。「起床が早く、バランスのいい朝食を食べている家庭は、生活にリズムがある」 よく読むと,この文には特になにも述べていない。ある種の家庭のありかたを,「生活にリズムがある」という言葉で形容しているだけである。にもかかわらずこの文章は,ああ生活のリズムが大事だなあ,と読み手に感じさせてくれる。内容がないのにわかりやすさだけが感じられる。これが技術である。
 
 では,責めるべきは新聞記者なのだろうか。彼らは文章上のテクニックによって,ただの相関関係を因果関係に見せかけ,「早寝早起き朝ごはん」国民運動のお先棒をかついでいるのだろうか。マクロにみれば,そういう批判も可能かもしれない。しかしミクロにみたとき,そこにあるのは単なる仕事熱心さだと思う。読み手の関心を少しでも惹き付けようとするのは,書き手にとって当然の工夫である。
 教育学者たちを責めるべきだとも思えない。彼らがなにを推進しようとしているにせよ,子どもの家庭生活と学校生活の関係を調べる調査には一定の価値があると思う。そのメディア向けリリースの際に,「朝ごはんが子どもにとって大事だ」という信念なり,社会的通念なりが顔を出しても,それは仕方のないことだ。
 この記事を読んで,ああ朝ご飯は大事だなあ,と素直に信じる人がいたとしても,その人を馬鹿にする気にはとてもなれない。調査結果からなにを読み取るのかはその人の勝手である。
 では,ブログでこの調査を罵倒し,マスコミと学者と無知な大衆を馬鹿にする人たちはどうだろうか。上で考えたように,その(いささか軽率な)批判にも,ある種の意義があるのかもしれない。でも,彼らは結局のところなにも破壊していないし,なにも建設していない。
 相関は因果を意味しない,という指摘は正しい。しかし,相関から因果を読み取ろうとする人々を押しとどめることはできない。新聞記事を読んで早起きは大事だと得心したお母さんは,なるほどあなたの仰る通りかもしれません,この調査結果は早起きが大事だということの直接の証拠にはなっていないんでしょうね,でもわたしやっぱり早起きって大事だと思うわ,というだろう。研究者たちにしたところで,仰るように学校の楽しさを支える因果的メカニズムについては今後の研究課題です,わざわざご指摘頂いてどうも,というのが落ちである。この手の批判者が何人いても,たとえば「早寝早起き朝ごはん」国民運動はとまらないのである。皮肉な言い方をすれば,そうやって批判することで自尊心を満たすことができたんだから,あなたもこの調査の受益者のひとりなんじゃないか,という気がしてしまう。
 どこにも悪い人はいない。それぞれの行いは正しい。でもマクロにみると,つまらない調査が量産され,根拠のない因果的説明が漠然と宙を漂い,社会を不可思議な決定へと導いていく。あえていうならば,調査者も伝達者も受け手も,全員が共犯なのではないだろうか。

 延々と書いてきたけれども,実はこの話はどうでもよくて,書きながら念頭にあるのは,いまの仕事であるところの市場調査のことである。調査会社,発注側の担当者,調査報告の受け手の全員が,誠実に調査にあたり誠実に報告し誠実にそれを批判し,でも結局はあいまいな共犯関係から抜け出ることができない,そんな事態がありうるのではないか,そのループから抜け出すためにはどうしたらいいのか,などと考えこむわけである。

「早起きの子どもは学校が好きだ」

2005年11月30日 (水)

 勤め先の会社で,調査データをどう分析するかという相談をしていたら,相手がふとこんなことを云った。データを調べる前に,あらかじめ「見たいもの」を決めてしまうようなやり方は嫌なんです。なんだかそれって予定調和だという感じがする。
 そうつぶやいた若い女性は,データ解析に関してはほぼnoviceだといってよい。しかし,(会社の先輩に対して大変失礼な言い方だが) とても勘の良い人で,明確で筋道だった問題意識を持っている人だった。
 どんな話の流れの中でそんな言葉がでてきたのか,よく思い出せない。その時は特に強い印象を受けず,なにか適当な受け答えで済ませてしまったような気がする。ずっとあとになって,あの言葉は俺に対する間接的な批判ではなかったか,と考えるようになった。正確にいえば,それは意図的な批判ではない。そのときの文脈やそのあとのやりとりから考えても,その人自身は俺に対して反論する意図を全く持っていなかったし,その言葉も深い意味合いをこめたものではなかったと思う。しかし俺にしてみると,その人に対して俺が繰り返し強調していたことについて,その問題点を一言で言い当てられたような,首筋にそっとナイフを当てられたような気がしたのだった。

 俺がその人に主張していたのはこういうことだ。ここにデータ行列があるから,それを虚心に眺め,なにか面白い事柄をみつけよう --- こういう発見的な方略は,魅力的ではあるけれども,勝算がほとんどない。まず「なにが見たいのか」を決めましょう。どんなナイーブな思いつきでもいい,分析に際しての視点を定めましょう。
 俺の考え方は決して間違っていないと思う。何の予断も持っていないように思えても,実はデータを手にした段階で,我々は現象を捉えるなんらかの枠組みを持っているものだ。我々はその枠組みに沿って問いを立て,その枠組みに沿って結果を解釈する。だから,その枠組みを顕在化・明確化することがまず必要だし,それによって「見たいもの」は自ずから定まるはずだ。強い言い方をすれば,分析に際して「見たいもの」がはっきりしないというのは,当該の現象そのものについてよほど関心がないか,でなければ物事をきちんと考えていないということだと思う。
 さらにいえば,実際の分析手順としてみても,ある程度までは仮説検証的なアプローチを採るのが現実的だ。データ行列が含んでいる情報はあまりに豊かであり,データをインタラクティブに視覚化して仮説探索していく作業は,仮説検証的な分析よりもはるかに高いスキルを必要とする。要するに,仮説を探すよりも検討する作業のほうがラクなのだ。考えてみればそれは当たり前の話で,仮説検証と仮説探索のちがいは,現象を抽象化するというしんどい作業をあらかじめ済ましておくかどうかというちがいでもある。
 でもその一方で,あらかじめ視点を限定すればするほど,そこから得られる知見は痩せていく。現象駆動的な問題意識からスタートしている場合,確たる理論的背景がないこと自体は悪いことではないし,そのときデータとの対話の中で自分の分析枠組みを明確化していこうとするのはひとつの見識だし,そうすることによってしか得られないタイプの視点もあるはずだ。俺はそういう要素をできるだけ減らして,コストと納期を見積もりやすい分析へと向かおうとしていたわけだが,それが一種の出来レースだという批判は,なるほど一面の真実である。

 というのは,ずいぶん奇麗な言い方だ。
 要するに怖いのです。データと際限なく格闘する泥沼に突き進むのが。答えがないかもしれない問いに巻き込まれるのが。容赦なく時間だけが経ち,はっきりするのは自分の無能さだけ,そういう恐怖をあなたは知っていますか(誰に向かって訴えているんだかわからんが)。真実なんてどうだっていい,平和が欲しい,幸せになりたい。そういうお年頃なのです。

予定調和

2005年10月10日 (月)

順序変数yが個人特性ηを測定しているとしよう。 ηのもとでyがc番目のカテゴリより上に落ちる確率P(y≧c|η)は, 比例オッズモデルなら
  P(y≧c|η) = F[(カテゴリcの下側の閾値)+(傾き)×η]
と表現するところだし,IRTなら
  P(y≧c|η) = F[(識別力)×(個人特性η-(カテゴリcの困難度))]
と表現するところだが,いったん連続的な潜在反応変数
  y* = (切片ν) + (因子負荷λ)×(個人特性η) + (誤差ε)
を考えてやって,
  P(y≧c|η) = F[(y* - カテゴリcの下側の閾値τ_c) / (誤差εのSD)]
と表現しても同じことである。ここでy*は,yの背後にある実質的に意味のある変数だと考えてもいいし,ただの計算上の道具だと割り切ってもかまわない。
 y*は潜在変数だから,平均0,分散1と考えるのが普通だ。その場合は
  (切片 ν) = 0
  (因子負荷 λ の二乗) × (ηの分散 ψ) + (誤差の分散 θ) = 1
となる。でも,別にy*の分散は1でなくてもいいわけで,たとえばθが1だと考えたっていい。そこで,
  Δ = 1 / (y*のSD)
をスケール・ファクタと呼ぶことにする。スケール・ファクタは,λとψとθをコミにしてあらわしたものである。単群の分析であれば,どんな値に決めても構わない。ここまでは,まあいいや。
 これが多群の分析となると,話がややこしくなってくる。群によってλだかψだかθだかがちがってくる,と考えないといけないからだ。そこで,ふたつの路線が登場する。ひとつは,ある群のy*の分散を(つまりスケール・ファクタΔを)1にして,ほかの群のΔを推定する,という考え方で,これをdelta approachという。このアプローチの背後には,どうせ誤差分散θは群間でちがうんだからどうでもいいや,という考え方がある。いっぽう,誤差分散θがちがうかどうかに関心がある場合には,ある群のθを1にして他の群のθを推定する,という考え方もできる。これをtheta approachという。
 モデルの推定上は前者のほうが都合が良いんだそうな(簡単だからだろうか)。

以上,Mplus Web Notes #4 より。カテゴリカル変数をつかったSEMについてわかりやすく説明したものが見つけられなくて困っていたのだが,これを読んでやっと少しだけ理解できた。Muthenさんの説明はなんでこんなにわかりやすいのだろうか。お菓子でも送りたい気分だ。

カテゴリ変数をつかったSEMでの多群分析

2005年9月26日 (月)

そんな分析では独立変数が他の変数とごっちゃになってしまっているですよ,ということの説明に小一時間をかけ,しかもまだ完全な理解を得ていない。話が通じない相手ではないのに,なぜだろうか。
(理由1) データの分析に際して問題意識を強く持っているから。それ自体はもちろん良いことなんだけど,往々にして自分の分析手法の欠陥が見えにくくなってしまうことが多いように思う。
(理由2) 学生ならば適当なところで納得してくれるが,オトナはそうではない。俺の説明能力は,俺が思っていたよりもはるかに低いのかもしれない。

話は違うが,屁が臭すぎて頭が痛い。自分の屁だから逃れようもない。どうすればよいのか。甘栗を食うのを止めれば済む問題なのか。

交絡と屁

rebuilt: 2017年5月18日 17:22
validate this page