elsur.jpn.org >

2019年1月28日 (月)

年末年始で本の記録を整理していて、この2年ほどは小説から遠ざかっているという危機感を持ち、1月は移動時間中に極力フィクションを読むよう心掛けていた。といっても、コートのポケットにはいる海外ミステリに偏ってしまうのだが...

Bookcover 償いの雪が降る (創元推理文庫) [a]
アレン・エスケンス / 東京創元社 / 2018-12-20
えーと、青春ミステリっていうのかしらん? こういうのだと知ってたら読まなかったタイプの小説なんだけど、ま、面白かった。

Bookcover 凍てつく街角 (ハヤカワ・ポケット・ミステリ) [a]
ミケール カッツ クレフェルト,Michael Katz Krefeld / 早川書房 / 2017-02-09
デンマークの人気作家だそうだ。心に傷を抱え酒浸りの刑事が、失踪事件を追うなかで自分を取り戻していく、というような話。

Bookcover あなたの人生の物語 (ハヤカワ文庫SF) [a]
テッド・チャン / 早川書房 / 2003-09-30
えーと、これはSF短編集で、表題作はおととし公開の映画「メッセージ」の原作。映画を観た直後にだいたい読んでいたんだけど、整理のためこのたび読了。「バビロンの塔」がいちばんおもしろかった。

Bookcover あなたを愛してから (ハヤカワ・ミステリ1933) [a]
デニス・ルヘイン / 早川書房 / 2018-07-05
人気作家ルヘインの新作。昨年読んだ「ザ・ドロップ」が全然印象に残ってないので、あんまり期待してなかったんだけど、いやー、これはまさにノンストップ・スリラーだわ...

フィクション - 読了:「あなたを愛してから」「あなたの人生の物語」「凍てつく街角」「償いの雪が降る」

Bookcover 観音さま (講談社学術文庫) [a]
鎌田 茂雄 / 講談社 / 2018-11-11

哲学・思想(2011-) - 読了:「観音さま」

Bookcover 豆腐と生きる 2 今、伝えたいリーダーの言葉 豆腐メーカー経営者インタビュー集 [a]
フードジャーナル編集部 / 株式会社フードジャーナル社 / 2016-07-29
豆腐好きが高じてついつい買ってしまった。食品業界誌による、豆腐メーカー経営者のインタビュー集。

Bookcover オスマン帝国-繁栄と衰亡の600年史 (中公新書) [a]
小笠原 弘幸 / 中央公論新社 / 2018-12-19

Bookcover 日本の「中国人」社会 (日経プレミアシリーズ) [a]
中島 恵 / 日本経済新聞出版社 / 2018-12-11

Bookcover 俳句の世界 (講談社学術文庫) [a]
小西 甚一 / 講談社 / 1994-12-27
偉い先生に対して失礼だけど、柔らかい文章ってのは古くなるんだなあ、と痛感...

Bookcover ドキュメント 候補者たちの闘争――選挙とカネと政党 [a]
井戸 まさえ / 岩波書店 / 2018-12-14
著者は岩波新書「日本の無国籍者」を書いた人で、経済誌出身のライターの方として名前を憶えていたのだけれど、実は民主党の衆院議員だったのだそうだ(現在は落選中で立憲に所属)。政党の候補者選定に焦点を当てたノンフィクション。個人的体験としてもひどい目にあっている由で、いやあ、政治家って大変だなあ。

ノンフィクション(2018-) - 読了:「豆腐と生きる2」「オスマン帝国」「日本の『中国人』社会」「ドキュメント 候補者たちの闘争」「俳句の世界」

Bookcover 壁の向こうの住人たち――アメリカの右派を覆う怒りと嘆き [a]
A.R.ホックシールド / 岩波書店 / 2018-10-26
「感情労働」概念で知られる社会学者ホックシールドが、トランプ支持層を取材して書いた本。非常に面白かった。これの大坂維新支持層版とか、どなたか書いて下さらないかしらん。それにしても、こういう研究者がインタビューを通じて定性的解釈(著者のいう「ディープ・ストーリー」)に至る思考のプロセスって、どうなってんだろう...

Bookcover 戦略の世界史(下) 戦争・政治・ビジネス [a]
ローレンス・フリードマン / 日本経済新聞出版社 / 2018-09-26
年末年始に気合入れて読んだ分厚い上下巻本。読んでるときは面白かったけど、どんな内容だったかよく覚えてない...

Bookcover 聖遺物崇敬の心性史 西洋中世の聖性と造形 (講談社学術文庫) [a]
秋山 聰 / 講談社 / 2018-10-12

Bookcover フランス現代史 (岩波新書) [a]
小田中 直樹 / 岩波書店 / 2018-12-21

Bookcover 花の命はノー・フューチャー: DELUXE EDITION (ちくま文庫) [a]
ブレイディみかこ / 筑摩書房 / 2017-06-06

ノンフィクション(2018-) - 読了:「花の命はノーフューチャー」「戦略の世界史」「壁の向こうの住民たち」「フランス現代史」「聖遺物崇敬の心性史」

Bookcover パンダ探偵社 1 (torch comics) [a]
澤江ポンプ / リイド社 / 2018-12-27
徐々に人間性を失っていく奇病をめぐる人間ドラマ。これは面白いなあ..

Bookcover アルテ 10 (ゼノンコミックス) [a]
大久保圭 / 徳間書店 / 2019-01-19

Bookcover 昴とスーさん 1巻 (ハルタコミックス) [a]
高橋 那津子 / KADOKAWA / 2017-08-10

Bookcover 北北西に曇と往け 3 (ハルタコミックス) [a]
入江 亜季 / KADOKAWA / 2019-01-15

Bookcover 僕らの色彩(1) (アクションコミックス) [a]
田亀 源五郎 / 双葉社 / 2019-01-12

Bookcover プリンセスメゾン (6) (ビッグコミックス) [a]
池辺 葵 / 小学館 / 2019-01-11

コミックス(2015-) - 読了:「パンダ探偵社」「プリンセスメゾン」「僕らの色彩」「北北西に曇と往け」「アルテ」「昴とスーさん」

溜めてしまうと面倒になるので、こまめに記録することにしよう... というわけで、年明けから読んだ本。まずはコミックスから。

Bookcover 伊豆漫玉ブルース (ビームコミックス) [a]
桜 玉吉 / KADOKAWA / 2019-01-12
: 孤高のエッセイマンガ家・桜玉吉さんのコミック・ビームの掲載作をまとめたもの。伊豆に「ボ泣き石」っていうのがあるんですね...

Bookcover メタモルフォーゼの縁側(1) (単行本コミックス) [a]
鶴谷 香央理 / KADOKAWA / 2018-05-08

Bookcover こはぜ町ポトガラヒー ~ヒト月三百文晦日払~ (1) (ビッグコミックス) [a]
昌原 光一 / 小学館 / 2018-12-27

Bookcover 結ばる焼け跡 (ニチブンコミックス) [a]
雨瀬 シオリ / 日本文芸社 / 2018-12-28

Bookcover かっぱのねね子 こうの史代小品集 [a]
こうの史代 / 朝日新聞出版 / 2018-12-13
こうの史代さんの「夕凪の街 桜の国」以前、まだマイナーなほのぼのギャグマンガ家であった頃の作品を収録。週刊モーニングの懸賞ページのイラスト、なんていうのも載っている。すでに絵柄は完成されていたんだなあ...

Bookcover 天地創造デザイン部(3) (モーニング KC) [a]
たら子 / 講談社 / 2019-01-23

コミックス(2015-) - 読了:「天地創造デザイン部」「かっぱのねね子」「伊豆漫玉ブルース」「メタモルフォーゼの縁側」「こはぜ町ポトガラヒー」「結ばる焼け跡」

2019年1月21日 (月)

 統計的因果推論の分野における神々のひとりJudea Pearlの主著 Causality は、そのとてつもない難解さで有名である... 少なくとも私のなかでは有名である。訳書でさえ限りなくギリシャ語に近い(訳者の先生すいません)。もっとも、最近の著書"The Book of Why"はもう少し親しみやすそうなんだけど、あいにく、オフィスの机の脇に聳える必読書たちの山の、二合目あたりに埋もれている。
 Pearlさんが打ち出した重要な概念のひとつにdo(x)オペレータというのがある。説明は端折るけど、まあとにかく、そういうのがあるんです。その意味するところは、なかなか難しい。少なくとも私にとっては、わりかしギリシャ語っぽい。
 このたびPearlさんご本人が、do(x)オペレータの解釈についての短い文章をご発表になった。で、みずから「書いたから読んでね」とメールでアナウンスなさっていた(こういうところ、Pearl先生は実にマメである)。えーと、そのメールによれば、RubinとかRobinsとかHeckmanとかHermanとか(←Pearlさんに並ぶ神々たち)の「操作できない変数に因果的特性を帰属させることはできない」という主張は、もはや過去の遺物である、とのこと。ご覧ください、相変わらずのこの戦闘性を。ここだけの話、畏れながら私は、Pearlさんを心中ひそかにマッド・ドッグ先生とお呼びしている(これはインドネシアのアクション・スター、ヤヤン・ルヒアンさんの通り名であり、ファンたちの心からの畏敬の念を表しているのです)。

 というわけで、飯のついでに読んでみた。
 そんなに真剣なわけじゃないけど、ちょっと個人的な関心もあって... 以前なにかで、「Aさんが肺がんになった原因のひとつはAさんが喫煙者だったからだ」というのは理解できるけど、「Aさんが前立腺がんになった原因のひとつはAさんが男性であることだ」というのは「原因」という言葉の使い方としておかしい、というような議論を読んで、そうかなあ、とずっと心にひっかかっていた。

Pearl, J. (2019) On the interpretation of do(x). Technical Report R-486, UCLA Cognitive Systems Laboratory.

 いわく。
 この論文では、次の2つの問いに答える。

  • 問1, $Q = E[Y|do(X=x)]$ というときの$Q$の、数学的でない意味。モデルの数学的な特性は、現実の特徴を表す場合もあれば、表しているとは言いにくい場合もある(例, 「このモデルに含まれる変数の数は素数だ」)。$Q$にはどんな現実的な意味があるのか。
  • 問2, その意味を実証的に検証する方法はあるのか?

 問1について。
 モデル$M$に基づき観察研究をやって$Q$を識別できたとしよう。その使い方には次の3つがある。

 その1, $Q$は、将来利用可能であろう操作可能な介入がもたらす因果効果の理論的限界を表現している。
 $I = \{I_1, \ldots, I_n\}$を操作可能な介入の集合とする。アウトカム$Y$にもたらすそれらの効果を比べたい。これらの介入は、それ自体に直接介入することはできない$X$に効果を与え、それを通じて$Y$に影響するのだとしよう。
 さて、観察研究において$Q = E[Y|do(X=x)]$を識別し推定できたとしよう。このとき、$Q$はいずれかの介入が$Y$にもたらしうる究極的な効果である。つまり、個別の介入$I_i$についてその効果を評価したいという人にとっては、$Q$はすぐに使えるものではないかもしれない。しかし、$X$をより強く制御するために新しい介入方法を開発すべきかどうか決めたいという人にとっては、$Q$は価値を持つ。
 なお、$Q$は$X$を通じた$I_i$の影響の天井なのであって、$I_i$の影響そのものの天井ではないということに注意。$X$を経由しない効果もあるかもしれない。

 その2, $Q$は、現在操作可能な変数の因果的効果に対する制約として働く。
 $I$ → $X$ → $Y$というシンプルな線形モデルを考える。構造係数をそれぞれ$a$, $b$とする。交絡変数はなく、$I$から$Y$へのパスもないものとする。$X$は操作不能とする。
 人々はいう、$X$は操作可能じゃないから、$b$を因果効果と呼ぶのは変だと。そんなこたあない。介入$I$の$Y$に対する平均因果効果$ACE(I)$は、適切な正規化の下で
 $ACE(I) = E[Y|do(I+1)] - E[Y|do(I)] = a \times b$
である。$b$は実現可能で操作可能な介入$I$が持っている理論的特性である。
 非線形システムでもそうだ。$X=f(I, \varepsilon_x), Y = g(X, \varepsilon_y)$とすると
 $E(Y | do(I)) = \sum_x P(x|do(I)) E[Y|do(x)]$
である。$Y$への$X$の理論的効果$E(Y|do(x))$がゼロならば、介入$I$の因果的効果もゼロだ。
 線形システムに戻って、$X$→$Y$に媒介変数とか交絡変数とかがあるとしよう。そのときでも
 $c = E[Y|do(x+1)] - E[Y|do(x)]$
とすれば$ACE(I)=a \times c$である。要するに、$Q = E[Y|do(X=x)]$を識別し推定できれば、$I$の因果的効果も同定できる。
 いや、もちろん$ACE(I)$そのものが推定できるんなら話は別で、その際には$b$や$c$には意味なかろう。でも$I$の開発中には、まだ$I'$しか開発できてないなんてこともあるわけだ。そのとき、$b$や$c$がわかったら大いに役に立つ。もしそれが小さかったら、もう開発なんてやめちゃおうかという話になるじゃないですか。

 その3, $Q$は、操作可能な変数群の因果効果を導出するための補足的な数学的構成物として働く。
 かつてDawid(2000 JASA)はこう述べた。我々は「反事実的な量」なるものを観察したことがない、「反事実的な量」についての我々のモデリング上の想定について、その妥当性を実証的評価できたこともない。反事実という概念に頼ることは哲学的にいって間違っておる、と。
 実証的妥当性についてはそのとおりだ。しかし、実用主義的な実証主義と教条主義的な実証主義を区別する必要がある。前者は、問いが実証的に検証できることを求めるが、その検証の際にどんなツールを使うかについては、便利さと想像力の問題と捉える。後者は、分析に登場するすべての補足的シンボル、すべての中間段階が、実証的な厳密さを持つ用語だけを含むことを求める。極端にいえば、負の値での割り算も認めないってことになるだろう。
 因果推論の文脈で言うと、実用主義者は、個々の測定単位の観察されない反事実(たとえば$Y_x(u)$)を喜んでみとめるだろう。それが集団効果についての実証的に検証可能な推定に役立つならば。Rosenbaum & Rubin(1983)の「潜在的アウトカム」というフレームワークでも、こうやって反事実が使われている。
 話を$Q$に戻そう。人々はいう、$X$が操作不能なとき、$Q$は実証性がないと。百歩譲って、$Q$は因果的効果でもなければ因果的効果の限界でもないとしよう。でも、少なくとも純粋に数学的な構成物ではある。それはそれ自体に経験的内容はないけれど、実証的に意味のある結果を導き出すことを可能にしてくれる構成物である。そういうのって、科学においては決して珍しくない。「複素数」とかね。

 問2について。
 $X$を直接操作することはできなくても、$Q(x) = q(x)$という主張を反証できるような研究計画を組み、データが主張と矛盾しないことを確かめることはできる。つまり、$Q(x) = q(x)$という主張をもたらしたモデリング上の想定から得られる、検証可能な含意を確認すればよい。
 たとえば、$Q(x)$がバックドア基準を通じて識別可能であり、バックドア基準を満たす共変量セットが手に入るから、それらの共変量を調整しよう、とか。
 $I$ → $X$ → $Y$であり、かつ$I$と$Y$に刺さる未知の共変量$U$があるので、あいにく観察研究で検証可能な含意はないが、RCTで$I$をランダム化して$P(y|do(I))$を推定しよう、とか。

 このように、操作不能な変数も因果的効果を持つことができる。それを検証することもできる。その変数が操作可能かどうかなんて、(実験計画はともかく)分析のフェーズでは、気にせんでよろしい。

論文:データ解析(2018-) - 読了: Pearl (2019) わがdo(x)オペレータの正しい解釈

2019年1月20日 (日)

 構造方程式モデリング界の帝王、泣く子も黙るMplus様には分厚いマニュアルがついていて、膨大な数のモデル例が紹介されている。そのなかには私には到底理解できないものもあれば、いっけん易しそうだがよく考えてみるとわけがわからないものもある。
 後者のひとつがExample 9.26「Cross-classifiedデータのランダム二値項目IRTモデル」で、よんどころない都合により、本日はずっとこの事例についてあれこれ思い悩んでいた。そりゃあ出世しないわ...!
 
 というわけで、この事例について考えたことをメモしておく。文字通り、私の・私による・私のためのメモであって、いちいちブログに書くんじゃないよ、情けないな、という気もする。

 Example 9.26とはこういうモデルである。Mplus 7.2から追加された TYPE=CROSSCLASSIFIEDを使っている。

VARIABLE: NAMES = u subject item; 
CATEGORICAL = u;
CLUSTER = item subject;
ANALYSIS: TYPE = CROSSCLASSIFIED RANDOM;
ESTIMATOR = BAYES;
PROCESSORS = 2;
MODEL: %WITHIN%
%BETWEEN subject%
s | f BY u;
f@1;
u@0;
%BETWEEN item%
u; [u$1];
s; [s];

 Mplusのサイトで公開されているデータを確認したところ、クラスタ変数 item は50水準, subjectは100水準。で、subject * item ごとに1行を持つ5000行のデータである。

 このモデルを数式で表現してみよう。
 
 コードに登場する変数は、u, item, subjectの3つ。そのうちitem と subject はCLUSTER=に指定されている。以下、subject を$i$, item を$j$とする。データの構造をコードから知るのは難しいが、上で確認したように、u は subject * item ごとにひとつ得られている二値の測定値である。これを$u_{ij}$と書く。

 VARIABLEセクションの注目点は2つ。

  • CLUSTER= item subject; となっている。あれこれ試していて気が付いたのだが、後述するように、CLUSTER=で指定する順番は無意味ではない。Mplusのマニュアルでは、cross-classifiedモデルにおける2つのbetweenレベルは2A, 2Bと表現されていて、CLUSTER=で最初に指定したのが2B, 次が2Aになる。ここではsubjectが2A, itemが2Bである。
  • WITHIN=もBETWEEN=もない。ここで、u の指定は次の4通りがあり得た。(1)WITHIN=に指定する、(2)BETWEEN=で subject レベルの変数として指定する、(3)BETWEEN=で item レベルの変数として指定する、(4)WITHIN=でもBETWEEN=でも指定しない。ここでは(4)が選択されているわけで、$u_{ij}$が、subjectのレベルでもitemのレベルでもモデル化される、ということがわかる。
     なお、$u_{ij}$がWITHINレベルでモデル化されるかどうかは、これだけではわからない(このことに気が付くまでに時間がかかった)。通常の階層モデル(TYPE = TWOLEVEL)では、ある変数がWITHIN=にもBETWEEN=にも指定されていないとき、その変数はどちらのレベルでもモデル化されるということがわかる。しかし、TYPE=CROSSCLASSIFIEDの場合は一概にそうともいえないのだ。このExample 9.26がよい例で、実は$u_{ij}$はWITHINレベルではモデル化されていない($u_{ij}$のモデルを式で表現したとき、その右辺に、$ij$を添え字に持つ項がない)。

 MODELセクションに移ろう。

 Mplusでは、カテゴリカル従属変数の背後には連続潜在変数が想定される。これを$u^*_{ij}$と書く。デフォルトではプロビット・リンクになる。つまり、標準正規分布関数を$\Phi$として
 $Prob(u_{ij}=1 | u^*_{ij}) = \Phi(u^*_{ij} - \tau)$
である。%BETWEEN item%に [u$1] と指定されているから、ついつい$\tau$が item ごとに動くのかと考えてしまったのだが、そうではない。(ここも私が混乱した点のひとつであった。考えてみれば、通常のTWOLEVELモデルで%BETWEEN%に [u$1]と指定したら、それは動かないよね。落ち着け俺)
 
 なお、Mplusのアウトプット上では、[u$1]は item レベルのパラメータとして扱われる。%BETWEEN item%に[u$1]と書こうが書くまいが結果は変わらない。%BETWEEN subject%に[u$1]と書くとエラーになり、ふたつあるBETWEENのどっちかで固定しろといわれる。そこで

MODEL:    %WITHIN% 
%BETWEEN subject%
s | f BY u;
f@1;
u@0; [u$1];
%BETWEEN item%
u; [u$1@0];
s; [s];
と書くと、推定は走るけど、パラメータu$1は0に固定され、結果として自由パラメータがひとつ減る。想像するに、MODEL=CROSSCLASSIFIEDにおいてカテゴリカル従属変数の閾値を推定する際は、レベル2B(つまり、CLUSTER=で1つめに指定したクラスタ変数のレベル)でしか指定できない、というルールがあるのかもしれない。

 ここからは$u^*_{ij}$のモデル。

 まずは、説明変数は無視して残差項にだけ注目しよう。uは%WITHIN%には登場せず、%BETWEEN subject%では u@0, %BETWEEN item%では u となっている。つまり、残差はsubject間では動かず、item間で動く。itemを「テスト項目」と捉えるなら、残差という言い方はなんだか変だな。むしろ切片というべきであろう。これを$a_j$と書き、その分散を$\sigma_a^2$と書こう。

 いよいよ説明変数について考える。はい深呼吸!
 %BETWEEN subject%に s | f BY u と指定されている。つまり、われらが $u^*_{ij}$は謎の潜在変数 f の指標なのだ。この潜在変数を $F$と書こう。[f]という記載がどこにもないので平均は0, f@1とあるので分散は1である。

 sは因子負荷である。おさらいしよう。Mplusでは、縦棒を使って傾きや因子負荷に名前を付けると、それはただの傾き・因子負荷ではなく、BETWEENレベルで変動するランダム傾き・ランダム因子負荷になる。通常の階層モデル(TYPE=TWOLEVEL RANDOM)では、%WITHIN%において s | y ON x; とかs | f BY y; と指定するのが常道である。このとき、WITHINレベルでは(つまりクラスタの中では) s は変動しない。ただし、ランダム傾きに限り、アスタリスクをつけて s* | y ON x; と書けば、sはクラスタ内でも変動する。
 しかし TYPE = CROSSCLASSIFIED RANDOMでは勝手が違う。まず、ランダム傾きだろうがランダム因子負荷だろうがアスタリスクは使えない。つまり、そのレベルでは変動しない。さらに、%WITHIN%だけでなく%BETWEEN%でも指定できてしまう。%BETWEEN%で指定した場合、そのレベルでは変動しないけど、隣のレベルでは変動しうるのだ。びっくり!
 ここでは、sは%BETWEEN subject%で定義され、%BETWEEN item%で s; [s]; とある。ってことは、このモデルでは、因子負荷 s は subject によっては動かないが、itemによって動くのだ。sを $s_j$と書き、その平均を$\mu_s$, 分散を$\sigma_s^2$と書こう。

 というわけで、この事例を数式で書き下ろすと、こういうモデルだと思う。
 $Prob(u_{ij}=1 | u^*_{ij}) = \Phi(u^*_{ij} - \tau)$
 $u^*_{ij} = a_j + s_j F_i $
 $E(F) = 0, \ \ Var(F) = 1$
 $a_j \sim N(0, \sigma_a^2), s_j \sim N(\mu_s, \sigma_s^2)$
推定するパラメータは、$\tau, \mu_s, \sigma_s^2, \sigma_a^2$の4つである。

 つまりこのモデルは、50指標1因子のカテゴリカル因子分析モデル、ないしIRTモデルといってよいと思う。ちがいは、因子負荷$s_1, \ldots, s_{50}$, 項目切片$a_1, \ldots, a_{50}$をモデルのパラメータとみるのではなく、確率変数の実現値として捉えているという点である。

雑記:データ解析 - 覚え書き: Cross-classifiedデータのランダム二値項目IRTモデルをMplusで表現する

2019年1月18日 (金)

 cross-classifiedマルチレベルモデルについて大急ぎで勉強する必要が生じ、よく引用されている参考書である
Raudenbush, S.W., Bryk, A. (2002) Hierarchical Linear Models: Applications and Data Analysis Methods. Second Edition.
を急遽注文し、12章だけ読んだ。さぞや難しかろうと覚悟していたのだが、意外にもわかりやすい内容であった。
 若くてやる気のある学生のみなさんに私がひとつだけ勝てるとしたら、このように、高めの本をためらいなくバカスカ買っちゃうということだけなのである... 悲しい...

 以下はそのメモ。

 本章ではデータが階層的ではあるんだけどきれいな階層になっていなくて、複数のユニットにcross-classifyされているようなデータについて扱う。子供が地域と学校に属しているとか、子供が中学と高校に属しているとか、従業員を職業と産業で分類するとか。
 なお、子供がある学校に属していて、学校の中で実験群か統制群のどちらかに割り当てられている、というような場合と、構造は似ているけど全然違う。「子どもが地域と学校に属している」場合、分析者は結果をいろんな地域・学校に一般化したいと思っているので、地域と学校はランダム効果である。「子どもが学校と実験条件に属している」場合、いろんな地域に一般化したいとは思うが「実験条件の母集団」というのはないので、学校はランダム効果で実験条件は固定効果となり、実験条件をlevel-1変数にとったtwo-levelとなる。

 以下、子ども$i$, 地域$j$, 学校$k$に属する子ども$i$の到達度を$Y_{ijk}$とする。

 まずunconditonalなモデル(予測子のないモデル)について。
 withinモデルは
 $Y_{ijk} = \pi_{0jk} + e_{ijk}$
 $e_{ijk} \sim N(0, \sigma^2)$
 betweenモデルは
 $\pi_{0jk} = \theta_0 + b_{00j} + c_{00k} + d_{0jk}$
 $b_{00j} \sim N(0, \tau_{b00})$
 $c_{00k} \sim N(0, \tau_{c00})$
 $d_{00j} \sim N(0, \tau_{d00})$
 $b_{00j}$は地域$j$の効果, $c_{00k}$は学校$k$の効果, $d_{0jk}$は学校と地域の交互作用効果。セル(学校x地域)のサイズが小さい時は交互作用効果の推定は無理。
 このモデルはこう書き換えられる。
 $Y_{ijk} = \theta_0 + b_{00j} + c_{00k} + d_{0jk} + e_{ijk}$
要するに二元配置分散分析である。

 分散が$\sigma^2, \tau_{b00}, \tau_{c00}, \tau_{d00}$に分解されたので、ユニット内相関係数も3種類定義できることになる[...メモ省略]

 ではconditionalなモデルについて。
 地域特性を$W_j$, 学校特性を$X_k$とする。話を単純にするため、どちらも1変数、固定効果、2水準(実験群と 統制群)とする。個人特性を$a_{ijk}$とする。話を単純にするため、1変数2水準(男と女)とする。
 withinモデルは
 $Y_{ijk} = \pi_{0jk} + \pi_{1jk} a_{ijk} + e_{ijk}$
 $e_{ijk} \sim N(0, \sigma^2)$
 さて、betweenモデルは...

 まずは簡単に、地域の予測子の効果も学校の予測子の効果も固定と考えよう。以下、$\pi_{0jk}$と$\pi_{1jk}$を一発で書くために、$p=0,1$とする。
 $\pi_{pij} = \theta_p + b_{p0j} + \gamma_p W_j + c_{p0k} + \beta_p X_k + d_{pjk}$

 このモデルだと、地域の効果は全学校で共通、学校の効果は全地域で共通である。
 これを緩和すると
 $\pi_{pij} = \theta_p + b_{p0j} + (\gamma_p + c_{p1k}) W_j + c_{p0k} + (\beta_p + b_{p1j}) X_k + d_{pjk}$
地域特性と学校特性の交互作用をいれると
 $\pi_{pij} = \theta_p + b_{p0j} + (\gamma_p + c_{p1k}) W_j + c_{p0k} + (\beta_p + b_{p1j}) X_k + \delta_p X_k W_j + d_{pjk}$

 実際には、理論と照らし合わせて、もうちょっと簡略にしないといけないわけで...
 ここからは実例の紹介。

 実例1, 子どもの達成に地域と学校がどう効くか研究。
 withinモデルは
 $Y_{ijk} = \pi_{0jk} + \pi_{1jk} a_{ijk} + e_{ijk}$
 $e_{ijk} \sim N(0, \sigma^2)$
 betweenでは主効果だけ考えた。
 $\pi_{0jk} = \theta_0 + b_{00j} + c_{00k}$
 $b_{00j} \sim N(0, \tau_{b00})$
 $c_{00k} \sim N(0, \tau_{c00})$
 このモデルを推定して級内相関を求めると...[略]

 予測子をいれてみます。withinレベルで、入学前の言語能力VRQと読み到達度Reading, 父親の職業DADOCCと教育DADED, 母親の教育MOMEDと働いているかDADUNEMP, 性別SEX。地域特性は、社会的剥奪レベルDEPRIVATION。ただし、学校によって違ってくるかもしれないので
 $\pi_{0jk} = \theta_0 + (\gamma_{01} + c_{01k})DEPRIVATION_j + b_{00j} + c_{00k} $
とし、$(c_{00k}, c_{01k})'$の分散共分散を3つすべて推定した。で、カイ二乗検定で「$c_{01k}$の分散も$c_{00k}$との共分散も0」という帰無仮説が棄却されなかったので、$c_{01k}$を取っ払った。[←なるほど...]

 実例2はめんどくさいので読み飛ばした。

 勉強になりましたです...
 実例1で、地域特性の効果が学校によって違うってのをモデル化しているけど($(\gamma_{01} + c_{01k})W_j$という項をいれる。なるほどねえ)、こういうのMplusだとどうやって指定するんだろう?

 そうか、私がいまやりたいのは
 $Y_{ijk} = \theta_0 + b_{00j} + b_{01j} c_{00k} + e_{ijk}$
というモデルだということに気が付いた。うーむ、こういうのはMplusだとどうなるんだ???

雑記:データ解析 - 覚え書き:cross-classifiedなマルチレベルモデル

2019年1月14日 (月)

Chow, S., Ho, M.R., Hamaker, E., Dolan, C.V. (2010) Equivalence and Differences Between Structural Equation Modeling and State-Space Modeling Techniques. Structural Equation Modeling, 17, 303-332.

 仕事の都合で泣く泣く読んだ奴。いま気がついたのだが、著者に動的因子分析のHamakerさん, 状態空間モデルのDolanさんがはいっている。先に著者名をちゃんとみておこうよ、俺...

 SEM(構造方程式モデリング)と状態空間モデリングのどこが同じでどこがちがうか説明しましょう、という論文。
 壮大な問題設定なのでびびってしまうが、あとでみるように、この論文がSEMと呼んでいるのはいわゆるSEMよりちょっと狭くて、測定変数間にパスがなく、因子に外生的な説明変数がないようなやつのことである(つまり、CFAの因子間にパス引いた奴)。SEM関連の本は全部オフィスの本棚に並べておるのでいまわからんが、普通そういうモデルはなんて呼ぶんですかね。多重指標モデル?

 イントロは省略して...

 まずSEMについて説明しよう。[原文の添字を簡略化してメモする]
 $y_i = \tau + \Lambda \eta_i + \epsilon_i$
 $\eta_i = \alpha + B \eta_i + \zeta_i$
ただし$y_i$, $\tau$, $\epsilon_i$は長さ$p$のベクトル, $\eta_i$, $\alpha$, $\zeta_i$は長さ$w$のベクトル。$\epsilon_i$の共分散行列を$\Theta$, $\zeta_i$の共分散行列を$\Psi$とする。
 母共分散行列と母平均は
 $\Sigma = \Lambda(I-B)^{-1} \Psi (I-B)^{-1'} \Lambda_{'} + \Theta$
 $\mu = \tau + \Lambda (I-B)^{-1} \alpha$
 SEMの場合、複数時点で測定しているデータについては、ふつうは$y_i$を二次元にするのではなくて長さを伸ばす。$T$時点の縦断データだったら、$T$時点がすべて$y_i$, $\eta_i$に叩き込まれるわけである。
 パラメータ推定にあたっては、ふつはMVNを仮定して尤度を求める。この場合は標本共分散行列$S$と標本平均ベクトル$m$があればよい。でもローデータから尤度を求めることもできる(以下RMLと呼ぶ。FIMLともいう)。時点数が個人数を超えている反復測定データなんかの場合はRMLを使わざるを得ない。
 RMLでは、個々の個人のデータベクトルの対数尤度を足し上げた
 $LL(\theta) = (1/2)\sum_i [-p \log(2\pi) - \log|\Sigma| - (y_i - \mu)^{'} \Sigma^{-1} (y_i - \mu)]$
を最大化する。MVNを仮定すると、この式はこう書き換えられる。
 $F(\theta) = (1/2)[\log|\Sigma|+tr(S\Sigma^{-1}) - \log|S| - p] + (1/2)[(m-\mu)^{'} \Sigma^{-1} (m-\mu)]$
 第1項はモデルの対数尤度、第1項は完全に飽和したモデル(つまり$\Sigma = S, \mu=m$であるモデル)の対数尤度である。これに$-2 (N-1)$を掛けるとカイ二乗統計量になる。[途中から写経モードに... まあ信じますよ先生]

 今度は状態空間モデルについて。
 $y_{it} = \tau + \Lambda \eta_{it} + \epsilon_{it}$
 $\eta_{it} = \alpha + B \eta_{i,t-1} + \zeta_{it}$
標準的な状態空間表記ではすべてのパラメータを時変とみるのだが、ここでは不変とみる。また外生変数もないことにする。
 パラメータ推定はふつうカルマン・フィルタ(KF)でやる。
 まず$t=0$における状態ベクトル$\eta_{i,0|0}$とその共分散行列$P_{0|0}$を決める。たとえば前者は0のベクトル、後者は$\kappa I$($\kappa$は十分に大きな数)とか。さらに、$\epsilon_{it} \sim N(0, \Theta)$, $\zeta_{it} \sim N(0, \Psi)$と仮定する。
 次に、各時点における$i$の状態ベクトルとその共分散行列を、それまでの観察によって予測する。
 $\eta_{i,t|t-1} = \alpha + B \eta_{i,t-1|t-1}$
 $P_{t|t-1} = B P_{t-1|t-1}B^{'} + \Psi$
以下では
 $e_{i,t|t} = y_{it}-(\tau + \Lambda \eta_{i,t|t-1})$
をイノベーションと呼ぶ(one-step-ahead予測誤差ともいう)。期待値は0である。
 さて、この2つの推定値を、その時点の観察を使って更新する。
 $\eta_{i,t|t} = \eta_{i,t|t-1} + K_{t|t} e_{i,t|t}$
 $P_{t|t} = [I-K_{t|t}\Lambda] P_{t|t-1}$
ただし
 $K_{t|t} = P_{t|t-1} \Lambda^{'} [\Lambda P_{t|t-1} \Lambda^{'} + \Theta]^{-1}$
これをカルマン・ゲインという。
 さて、全員について$T$時点の推定を行った暁には、イノベーション$e_{i,t|t-1}$とその共分散行列
 $G_t = E[e_{i,t|t-1}e_{i,t|t-1}^{'}]$
が手に入る。後者はすべての$i$について$\Lambda P_{t|t-1}\Lambda^{'}+\Theta$である。
 ここから対数尤度関数が手に入る。
 $LL(\theta) = (1/2)\sum_i \sum_t [-p \log(2\pi) - \log|G_t| -(e_{i,t|t-1}^{'})G_t^{-1} (e_{i,t|t-1}) ]$
これを予測誤差分解(PED)関数という。これを最大化すればいい。
 要約しよう。まずパラメータに初期値を与える。つぎにカルマンフィルタで個々人の各時点での状態を推定する(このときパラメータは定数とみなす)。次に状態推定を既知と見なしてパラメータをPED関数で推定する。このパラメータを定数とみなして新しい状態を推定する。これを繰り返す。

 いよいよ本題。SEMと状態空間モデリングは、どこが同じでどこがちがうか。
 ひとことでいうと次の通りである。どちらも「プールされた個人内ダイナミクス」を表現している。SEMは潜在変数間の同時的な構造的関係と、その関係の個人間の差を捉えるのに適している。状態空間モデルは、もっと複雑な個人内ダイナミクスを表現するのに適している。

 もっと細かくいうと次のとおり。以下、SEMのパラメータに添字$S$, 状態空間モデルのパラメータに添字$K$をつける。
 状態推定について。$B_K=0$とすれば、カルマンフィルタは回帰法による因子得点推定と同じである。また$B_K=0$で$\Psi_K$が無限大に近いとき(つまり事前情報がないとき)、カルマンフィルタはバートレット法による因子得点推定に近づく。[...中略...]
 パラメータ推定について。$T=1$の場合、$B_S = B_K = 0$となり、2つのモデルは等しい。$T > 1$の場合、$p_S = T \times p_K$ならば、対数尤度関数が同じになる。さらに$B_S = B_K = 0$なら、パラメータ推定値は同じになる。
 問題は$B_S, B_K$が0でない場合。特別な制約を付けない限りパラメータは変わってくる。その制約がどんなものかは、実例をみたほうが早かろう。

 実例1, 横断共通因子モデル。ここでみなさんに、「状態空間モデルはある個人の反復測定データを扱う手法である」というよくある思い込みから脱して頂きます。

 SEMの場合。$B$が0になるので、
 $y_i = \tau + \Lambda \eta_i + \epsilon_i$
 $\eta_i = \alpha + \zeta_i$
共分散行列を$\Theta = E(\epsilon_i \epsilon_i^{'}), \ \Psi = E(\zeta_i^2)$として、
 $\Sigma = \Lambda \Psi \Lambda^{'} + \Theta$
 $\mu = \tau + \Lambda \alpha$
 $LL = (1/2) \sum_i [-p \log(2\pi) - \log |\Sigma| - (y_i - \mu)^{'} \Sigma^{-1} (y_i - \mu)]$
対数尤度関数は状態空間モデルと同一になる。

 [ほんとなの? 試してみよう。
 $y_{it} = \tau + \Lambda \eta_{it} + \epsilon_{it}$
 $\eta_{it} = \alpha + \zeta_{it}$
イノベーションは本来は$e_{i,t|t-1} = y_{it}-(\tau + \Lambda \eta_{i,t|t-1})$だけど、$\eta_{i,t|t-1} = \alpha$だから
 $e_{i,t|t-1} = y_{it}-(\tau + \Lambda \alpha)$
 $G_t = \Lambda \Psi \Lambda^{'}+\Theta$
 $LL = (1/2) \sum_i \sum_t [-p \log(2\pi) - \log|G_t| -(e_{i,t|t-1}^{'})G_t^{-1} (e_{i,t|t-1}) ]$
なるほど、同じだわ]

 シミュレーション。$N=200, T=1, p=4$のデータについて$w=1$の共通因子モデルを組みます。パラメータの真値を適当に決めて架空データをつくり、LISRELとmkfm6(という状態空間モデルのプログラムがある由)で分析してみる。結果はほぼおなじ。

 実例その2, 単変量LDSモデル(Nが大きくてTが小さい場合)。
 McArdle & Hamagamiのlatent difference scoreモデルについて考える。(ああ、読んでおいて良かった... なにがいつ役に立つかわからんねぇ...)
 LDSモデルでは、観察の後ろに潜在変数$\eta_{it}$を考え、さらに「潜在変数の前期からの差」という潜在変数$\Delta \eta_{it}$を考える。
 $y_{it} = \eta_{it} + e_{it}$
 $\eta_{i, t+1} = \eta_{it} + \Delta \eta_{i,t+1}$
でもって、
 $\Delta \eta_{it} = \eta_{si} + \beta \eta_{i,t-1}$
と考える。$\eta_{si}$は$i$さんの定数項で変化の強さを表す。$\beta$は全員共通の定数項で、自己フィードバックを表す。以上をまとめると
 $\eta_{i,t+1} = \eta_{si} + (1+\beta) \eta_{it}$
$i$さんの時点0における切片項と変化の強さを、平均と偏差に分解して
 $\eta_{0i} = \mu_{\eta 0} + \zeta_{\eta 0 i}$
 $\eta_{si} = \mu_{\eta s} + \zeta_{\eta s i}$
$E(\zeta_{\eta 0 i}^2) = \sigma^2_{\eta 0}$, $E(\zeta_{\eta s i}^2) = \sigma^2_{\eta s}$とする。
 LDSモデルは要するに成長曲線モデルだと考えてよい。ちがいは、潜在傾きをさらにモデル化しているという点。

 SEMで表現してみましょう。
 長さ$T$のベクトル$y_i$を考える。測定モデル
 $y_i = \Lambda \eta_i + \epsilon_i$
で、潜在変数$\eta_i$を、長さ$T+1$のベクトル$[\eta_{i1}, \ldots, \eta_{iT}, \eta_{si}]'$とする。$\Lambda$は$T \times T$の単位行列の右に、すべて0の列をくっつけた奴にする。要するに、$\eta_{si}$は測定モデルには出てこないんだけど、話の都合上、潜在変数ベクトルにいれておくわけである。トリックですね。
 構造モデル
 $eta_i = \alpha + B \eta_i + \zeta_i$
では、$\alpha$は先頭が$\mu_{\eta 0}$、末尾が$\mu_{\eta s}$、ほかはすべて0のベクトルにする。$\zeta_i$も同様で、先頭が$\zeta_{0i}$、末尾が$\zeta_{si}$、ほかはすべて0である。
 $B$はどうなるか。時点1については(つまり$B$の1行目)、すべて0でよい。時点$2,\ldots,T$については(つまり$B$の2行目から$T$行目)、対角に$1-\beta$をいれた$(T-1) \times (T-1)$の対角行列の右に、すべて0の列($\eta_{iT}$にかかる係数)、すべて1の列($\eta_{si}$にかかる列)をくっつければ良い。最後の行は$\eta_{si}$をつくるための係数だから、すべて0。ね? 実は$\Delta \eta_{it}$なんていらないのだ。[←なるほど...]

 状態空間で表現してみましょう。[原文とはちがうが、話がややこしくなるので、状態空間モデルでの潜在変数はアスタリスクをつけて$\eta^{*}$とする]
 観察方程式
 $y_{it} = \Lambda \eta_{it}^{*} + \epsilon_{it}$
で、潜在変数$\eta_{it}^{*}$をベクトル$[\eta_{it}, \eta_{si}]'$にする。$\eta_{si}$は時間不変だし、そもそも観察変数には出てこない奴なので、変な感じだが、これもトリックですね。$\Lambda = [1,0]'$である。
 遷移方程式
 $\eta_{it}^{*} = B \eta_{i,t-1}^{*} + \zeta_{it}$
はどうなるか[ここ、原文に誤植があると思うので勝手に直した]。$\zeta_{it}$はすべて0。$B$は$2 \times 2$行列で、1行目($\eta_{it}$をつくる係数)は$[1-\beta, 1]'$。2行目は、時間不変だから$[0,1]'$。
 最後に、初期値$\eta^{*}_{i1|1}$を$\mu_{\eta 0}, \mu_{\eta 1}$, その共分散行列$P_{i1|1}$の対角要素を$\sigma^2_{\eta 0}, \sigma^2{\eta s}$とする。

 シミュレーション。$N=400, T=10$のデータについてLDSモデルを組む。パラメータの真値を適当に決めて架空データをつくり、LISRELとOx(という状態空間モデルのプログラムがある由)で分析してみる。結果はほぼおなじ。
 というわけで、LDSモデルはSEMとしても状態空間モデルとしても推定できるんだけど、後者のほうが便利。なぜなら、SEMだと(1)$T>N$のときに困るし、(2)$T$が大きいだけでえらいことになる。

 実例その3,自己回帰動的因子モデル($T$が大きくて$N=1$)。
 長さ$p$の観察ベクトル$y_{it}$について
 $y_{it} = \Lambda \eta_{it} + u_{it}$
 $\eta_{it} = A_1 \eta_{i,t-1} + A_2 \eta_{i,t-2} + \ldots + \zeta_{it}$
というのがこのモデル[原文に誤植があると思うので勝手に直した]。以下、2因子、6指標、ラグ1までのモデルを考えます。

 状態空間表現の場合。観察方程式
 $y_{it} = \Lambda \eta_{it} + \epsilon_{it}$
の$\eta$は長さ2のベクトル。$\Lambda$は$6 \times 2$の行列で、1列目を$[1, \lambda_{21}, \lambda_{31}, 0, 0, 0]'$, 2列目を$[0,0,0,1, \lambda_{52}, \lambda_{62}]'$とする(識別のためにひとつ1をいれている)。
 遷移方程式
 $\eta_{it}^{*} = B \eta_{i,t-1} + \zeta_{it}$
の$B$は$2 \times 2$の行列で、ラグ1の自己回帰係数。$\zeta_{it}$はこのVAR(1)プロセスで説明できない因子得点である。

 SEMの場合。データ行列は$N \times T p$となるわけで、共分散行列をつかった最尤推定はできない。ひとつの路線はRMLを使うことだが、計算が大変である。もうひとつはブロック・トープリッツ行列を使うという路線。こっちで考えます。
 [あれ、ほかにも方法があるのではなかろうか。まあいいけどさ]

 $N=10, T=200$でモンテカルロシミュレーションしてみよう。[...メモ中略...] とこのように、ブロック・トープリッツ行列路線はパラメータ推定値がけっこう歪む。それに個人レベルの因子得点も推定できない。お勧めできない。

 考察。
 心理学のみなさんは状態空間モデルをあまり使わない。もったいない。SEMにも状態空間モデルにもそれぞれ美点があるので、使い分けましょう。
 云々、云々。

論文:データ解析(2018-) - 読了:Chow, Ho, Hamaker, Dolan (2010) SEMと状態空間モデルのどこが同じでどこが違うか

<< 読了: Im, et al.(2016) cross-classifiedマルチレベルデータをただの階層データとみなしたときの弊害 (多群CFAで測定不変性を検証する編)
 
validate this page / CSS