読書日記

« 2018年11月 | メイン | 2019年4月 »

2019年1月28日 (月)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

天地創造デザイン部(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$を識別できたとしよう。$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$が大きいだけでえらいことになる。

　実例その３，自己回帰動的因子モデル($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と状態空間モデルのどこが同じでどこが違うか

2019年1月12日 (土)

　目の前の仕事に日々追い立てられる毎日であるが、副業というか研究開発的な案件のほうもそれはそれで追い詰められていて、年末、階層モデルについて考えていて頭が混乱してしまい、ちょっとしたパニックに陥った。とりあえず、抱えている問題と関係ありそうなPDFを拾って印刷し、年末年始でだだだーっと読んでみた。そのなかの一本。要は測定不変性に関する話で、とっつきやすい。

Im, M.H., Kim, E.S., Kwok, O.M., Yoon, M., Willson, V.L. (2016) Impact of Not Addressing Partially Cross-Classified Multilevel Structure in Testing Measurement Invariance: A Monte Carlo Study. Frontiers in Psychology, 7, 328.

　いわく。
　教育や社会科学においてはマルチレベルデータを扱うことが多いけど、たいていの場合階層線形モデルが用いられている。そこでは、下のレベルの事例は上のレベルのあるひとつのクラスタに属する。たとえば、ある生徒はある教室に属し、ある教室はある学校に属する。
　しかし現実はそうでない。生徒はふたつの学校に通っていたりするし、学校と地域のように、一方が他方の下にあるわけではないこともある。こういう非階層的なマルチレベルデータをcross-classifiedマルチレベルデータという。
　近年ではcross-classifiedマルチレベルデータを正しく分析することが大事だという理解が広まっている。Goldstein(1986 Biometrika, 1995書籍), Rasbash & Goldstein (1994 J.Educ.Behav.Statist.), Raudenbush & Bryk(2002書籍) をみよ。教科書への記載も増えておるぞ。実証研究も増えておるぞ。Fielding (2002 Sch.Eff.Sch.Improv.), Jayasinghe et al.(2003 JRSS), Marsh et al.(2008 J.Edu.Psy.)をみよ。
　しかるに、一部の者どもはcross-classifiedデータをいまだ階層線形モデルで分析しておる。実に嘆かわしいことである。
　cross-classifiedデータを階層データとして扱うとどうなるか。研究は少ないんだけど、分散成分の推定にバイアスがかかると指摘されている。

　さて。近年の社会科学においては、異なる集団間なり時点間なりでの測定不変性(MI)の検証が一般的になっている。階層的マルチレベルデータでにおけるMIの検証にはマルチレベル確認的因子分析(MCFA)が広く用いられている。それはいいんだけど、実際のマルチレベルデータは必ずしも階層になっていないのに、それでもふつうのMCFAが使われているのである。一般的なSEMのソフトではcross-classifiedマルチレベルデータのMCFAができないからである。[←これを読んで慌てて調べたが、これはちょっと吹かし気味のご発言で、Mplusには2012年秋リリースのVer.7の時点でがTYPE=CROSSCLASSIFIEDが搭載されている。現在もっとも一般的なSEMソフトといえばMplusだと思うんですが。ま、要するにMplus買いましょうってことっすね]

　というわけで、本研究は、cross-classifiedマルチレベルデータが測定不変性を持っているかどうかを検討するとき、cross-classifiedであることを正しく扱わなかったらどうなるかを調べる。

　まずは、ふつうのMCFAとcross-classified MCFAのちがいについて説明しよう[←ありがとうありがとう。それが読みたかったのよー]。

　2レベル, 1因子, 4指標の場合を考える。
　ふつうのMCFAの場合、ある生徒 $i$ はある学校 $j$ にのみ属する。観察ベクトルを$X[i,j]$と呼ぼう。いっぽうcross-classified MCFAでは、$i$は学校$j_1$と地域$j_2$に属する。観察ベクトルを$X[i, (j_1, j_2)]$と書こう。[原文では$X_{ij}$, $X_{i(j_1, j_2)}$なのだが、添字が深くなりすぎるので表記を変える。以下同様]

　ふつうのMCFAでは、Withinレベルで$x_1, x_2, x_3, x_4$に因子FWからパスが刺さり、さらにBetweenレベルで、$x_1, x_2, x_3, x_4$に因子FBからパスが刺さる。いっぽうcross-classified MCFAでは、Betweenレベルが学校のパートと地域のパートに分かれる。学校のパートでは$x_1, x_2, x_3, x_4$に因子FB1からパスが刺さり、地域のバートでは$x_1, x_2, x_3, x_4$に因子FB2からパスが刺さる。

　式で書こう。まず単純に、潜在変数を$\eta$として、普通のMCFAでは
　$X[i,j] = \tau + \Lambda \eta[i,j] + \varepsilon[i,j]$
cross-classified MCFAでは
　$X[i,(j_1,j_2)] = \tau + \Lambda \eta[i,(j_1,j_2)] + \varepsilon[i,(j_1,j_2)]$
観察値はMVNに従うものとする。(IID正規である必要はない。階層データだから)

　クラスタ間で変動するランダム効果をいれよう。ふつうのMCFAの場合、$\eta[i,j]$が2つに分かれる。
　$\eta[i,j] = \alpha + \eta_w[i,j] + \eta_b[j]$
$\alpha$は$\eta[i,j]$の期待値というか全体平均である。
　いっぽうcross-classified MCFAの場合、$\eta[i,(j_1,j_2)]$は3つに分かれる。
　$\eta[i,(j_1, j_2)] = \alpha[j_1] + \alpha[j_2] + \eta_w[i,j] + \eta_b[j_1] + \eta_b[j_2]$
$\alpha[j_1], \alpha[j_2]$はFB1, FB2の期待値である。[FWの期待値は0ってことね]

　元の式のほうも分解してみよう。ふつうのMCFAでは
　$X[i,j] = \tau_b + \Lambda_w \eta_w[i,j] + \Lambda_b \eta_b[j] + \varepsilon_w[i,j] + \varepsilon_b[j]$
　MCFAでは切片はBetweenレベルにしかないという点に注意。
　これがcross-classified MCFAだと
　$X[i,(j_1, j_2)]$
　$= \tau_b[j_1] + \tau_b[j_2]$
　$+ \Lambda_w \eta_w[i,j] + \Lambda_b[j_1] \eta_b[j_1] + \Lambda_b[j_2] \eta_b[j_2]$
　$+ \varepsilon_w[i,j] + \varepsilon_b[j_1] + \varepsilon_b[j_2]$
となる。切片は2つになる。

　因子分散も分解できる。ふつうのMCFAだと
　$V(\eta[i,j]) = \Psi_w + \Psi_b$
で、全分散に占める$\Psi_b$の割合を級内相関(ICC)という。いっぽうcross-classified MCFAだと
　$V(\eta[i,(j_1, j_2)]) = \Psi_w + \Psi_b[j_1] + \Psi_b[j_2]$
となる。ICCが2つできることになる。

　独自分散も分解できる。ふつうのMCFAだと
　$V(\varepsilon[i,j]) = \Theta_w + \Theta_b$
　$V(\varepsilon[i,(j_1,j_2)]) = \Theta_w + \Theta_b[j_1] + \Theta_b[j_2]$

　最後に、共分散行列$\Sigma_T$の分解。ふつうのMCFAの場合、共通因子$\eta$と独自因子$\varepsilon$が独立だという仮定の下で、
　$\Sigma_b = \Lambda_b \Psi_b \Lambda^{'}_b + \Theta _b$
　$\Sigma_w = \Lambda_w \Psi_w \Lambda^{'}_w + \Theta _w$
として$\Sigma_T = \Sigma_b + \Sigma_w$である。同様にcross-classified MCFAでは...[面倒なのでメモは省略するが、$\Sigma_T$が3つに分かれる]

　さて、測定不変性の話。
　測定不変性というのは因子-指標間の非線形的関係を含んだ広い言葉だが、線形因子モデルにおける測定不変性(因子不変性, FI)とは、パラメータが群間で等しいことをいう。その検証のためには、MCFAモデルにパラメータ等値制約を掛けていく。因子パターンが同じ(配置不変)、負荷が同じ(メトリック不変)、潜在変数の切片が同じ(スカラー不変)、独自分散が同じ(厳密不変)、というふうに。[ここで因子パターンといっているのは、いわゆる因子パターン行列のことじゃなくて、因子負荷行列のどこが0でないか、ということであろう]
　cross-classified MCFAの場合だとこうなる。以下、群を$g$で表す。モデルは
　$X_[i,(j_1,j_2),g]$
　$= \tau_b[j_1,g] + \tau_b[j_2,g]$
　$+ \Lambda_w[g] \eta_w[i,j,g] + \Lambda_b[j_1,g] \eta_b[j_1,g] + \Lambda_b[j_2,g] \eta_b[j_2,g]$
　$+ \varepsilon_w[i,j,g] + \varepsilon_b[j_1,g] + \varepsilon_b[j_2,g]$
　共分散行列を分解して
　$\Sigma_b[j_1,g] = \Lambda_b[j_1,g] \Psi_b[j_1,g] \Lambda^{'}_b[j_1,g] + \Theta_b[j_1,g]$
　$\Sigma_b[j_2,g] = \Lambda_b[j_2,g] \Psi_b[j_2,g] \Lambda^{'}_b[j_2,g] + \Theta_b[j_2,g]$
　$\Sigma_w[g] = \Lambda_w[g] \Psi_w[g] \Lambda^{'}_w[g] + \Theta_w[g]$
　まず配置不変性を検討する(因子数と各因子の指標数が群間で等しいといえるかを調べる)。
　次にメトリック不変性を調べる。帰無仮説は
　$\Lambda_b[j_1,1]=\Lambda_b[j_1,2]= \cdots = \Lambda_b[j_1,G]$,
　$\Lambda_b[j_2,1]=\Lambda_b[j_2,2]= \cdots = \Lambda_b[j_2,G]$,
　$\Lambda_w[1]=\Lambda_w[2]= \cdots = \Lambda_w[G]$。
　次にスカラー不変性を調べる。帰無仮説は
　$\tau_b[j_1,1]=\tau_b[j_1,2]= \cdots =\tau_b[j_1,G]$,
　$\tau_b[j_2,1]=\tau_b[j_2,2]= \cdots =\tau_b[j_2,G]$。
　最後に厳密不変性を調べる。帰無仮説は...めんどくさいから省略するけど、$\Theta$が等価だという3つの仮説ね。

　ここからは実験をふたつ。

　実験1。学校が20個、地域が50個あり、一方が決まると他方がある程度まで決まるというデータ(部分的cross-classifiedデータ)をつくる。cross-classified MCFA、1因子4指標で、3つのどのパートでも負荷は0.7から1のあいだであり、とりあえずはレベル間で共通とする(因子不変)。独自分散はwithinレベルで0.25, betweenレベルで0.05。
　動かす要因は２つ。(1)因子不変か。学校を2群にわけ、片方の群についてのみ、学校レベルの因子負荷を変える。大きく変える、小さく変える、変えない、の3水準。(2)ICC。因子分散を動かす。3水準。
　で、地域を無視したふつうのMCFAで測定不変性を調べる。
　結果。カイ二乗検定では、本当は因子不変である場合のType Iエラーには問題がないんだけど、本当は因子不変でない場合でも全然検出できない。学校レベルの因子負荷$\lambda_b[j_1]$の群間差(DIF)は過小評価される。

　実験2ではcross-classified MIMICモデルで測定不変性を調べる(群間で切片が異なるかどうかを、群を因子の共変量にいれることで調べられるかという話であろう)。これもMplus 7.4で実験している。面倒になったのでスキップ。

　考察。
　実験1では、FB2を無視したせいで、$\Psi_b[j_2]$が$\Psi_w$に再配分されてしまい、ICCが低く評価され、DIFが検出できなくなってしまった。[←ああ、そういうことなのか、なるほどね。2要因のANOVAをやるべきところを1要因でやったら、そのぶんセル内の分散が大きくなって、要因の主効果が有意じゃなくなっちゃいました、と云うような話なのであろう]
　このように、cross-classifiedデータで、一方のbetweenレベルの因子負荷が群間で異なるかどうか調べる場合、もう一方の因子を無視したふつうのマルチレベルCFAをやるのはお勧めできない。
　云々。

　ちゃんと読んでないけど、実験1は「注目しているクラスタ変数とクロスしているクラスタ変数を無視して階層多群分析するとクラスタ間での測定変動性を見落とすよ」という話、実験2は「群間で切片が異なることを見つけるだけでよければ多群分析じゃなくて群を共変量にいれたMIMICモデルにするといいよ」という話なのではないかと思う。でもこの実験、Mplus 7.4でやったんでしょ？ なんでTYPE=CROSSCLASSIFIEDを使わないの？と不思議に思いながら読んでいたのだが、あー、そうか、CROSSCLASSIFIEDだとGROUPING=オプションが使えないのかも。

　ともあれ、cross-classified SEMというのがどういうものなのかが理解できたので良しとしよう。
　ほんとはさあ... Muthen一家の誰かが解説を書いてくれるといいんだけどさあ... とぶつぶつ不満を漏らすわけだが、別に新しい話題ってわけじゃないのだろうな。いまのホットトピックはきっと強縦断データとかだろうし。すいません、自分で勉強します。とりあえずはRaudenbush & Brykを手に入れるか...

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

Tellis, G. J. (2006) Modeling Marketing Mix. in Grover, R. & Vriens (eds), The Handbook of Marketing Research: Uses, Misuses, and Future Advances, Sage.

　年末に仕事の都合で大急ぎで目を通した奴。Marketing Mix Modeling (通称MMM。なんでも略すればいいってもんではないだろうに) について急遽整理する必要が生じ、参ったなと思ったら、このハンドブックに1章が割かれていた。困った時のGrover & Vriens である。
　著者はマーケティングの先生で、「意志とビジョン: マーケット・リーダーの条件」という邦訳書がある。しらんけど。

　いくつかメモ:

　著者曰く、広告に対する反応にみられる重要なパターンが7つある。(1) current effect, (2) carryover effect, (3) shape effect (広告-売上関係が曲線的だという話), (4) competitive effects, (5) dynamic effect (wearin, wearout, hysteresis), (6) content effects, (7) media effects. レベルが揃ってなくてすごく気持ち悪い...

　広告反応のモデルとして著者は次の5つを挙げている。

• 基本的な線形モデル: 
  $ Y_t = \alpha + \beta_1 A_t + \beta_2 P_t + \beta_3 R_t + \beta_4 Q_t + \varepsilon_t$
  $Y_t$は売上かなにか、$A_t$は広告かなにか(以下「かなにか」を省略), $B_t$は価格, $R_t$はプロモーション, $Q_t$は品質。$\varepsilon_t$はIID normalだと思っているわけですね。
• 乗法モデル: 
  $ Y_t = \exp(\alpha) \times A_t^{\beta_1} \times P_t^{\beta_2} \times R_t^{\beta_3} \times Q_t^{\beta_4} \times \varepsilon_t $
  推定の際には全部対数にする。利点: (1) シナジーをモデル化できている。(2)いろんな曲線を表現できる。(3)$\beta$が弾力性になっている。欠点: (1)carryover, competitive effects, dynamic effects, content effects, media effectsを表現できていない。(2)S字型の曲線を表現できない。(3)広告の弾力性は一定だと仮定していることになる。[いちいち書かないけど、頭のなかに疑問符が乱舞する... もっと厳密に書かれた説明を読みたいのに]
• 指数魅力モデルと多項ロジットモデル。めんどくさいので途中端折るけど、市場シェアを各ブランドのマーケティング努力の指数におけるシェアとし、マーケティング努力について線形モデルを組んで
  $ M_i = \exp(\sum_k \beta_k X_{ik} + e_i) / \sum_j \exp(\sum_k \beta_k X_{jk} + e_j)$
  という感じ(めんどくさいので独立変数は$X_{ik}$と書く...)。しかしこれだと係数を解釈しにくいし推定も大変なので対数中心化する。結局
  $ \log(M_i M^{-}) = \alpha^{*}_i + \sum_k \beta_k X_{ik}^{*} + e_i^{*}$
  となる。$M^{-}$ってのは幾何平均ね。これは結局多項ロジットモデルになっている。利点: (1)競合も考慮している、(2)関数形がS字型になっている、(3)市場シェアの独立変数に対する弾力性が、マーケティング努力にたいしてベル型になる。
• Koyckモデルと分散ラグモデル。まず、
  $Y_t = \alpha + \lambda Y_{t-1} + \sum_k \beta_k X_{ik} + \varepsilon_t$
  これがKoyckモデル。限界: (1)carryoverは減衰だと決めてかかっている。(2)個別の変数についてのcarryoverを推定するのは大変。(3)データの時間間隔によって結果が大きく変わる。そこでもっと一般化して
  $Y_t = \alpha + \lambda_1 Y_{t-1} + \lambda_2 Y_{t-2} + \ldots + \beta_{10} A_{t} + \beta_{11} A_{t-1} + \ldots $
  というように広げるのが分散ラグモデル。でも多重共線性がひどくなり推定が安定しない。データが足りない場合はKoyckモデルのほうが良い。
• 階層モデル。そのまえにダミー変数回帰について考える。たとえば広告の内容が2種類ある場合、一方であることを表すダミー変数$A_{2t}$をつくり、
  $Y_t = \alpha + \beta_1 A_t + \delta A_t A_{2t} + \ldots $ 
  という風に交互作用項をいれる。当然ながら、広告の種類が増えるとめっちゃ複雑になる。そこで階層モデルの登場...[めんどくさいのでメモ省略。というか、数式で書いてよ！かえってわけわかんなくなるから！]。データがリッチなら、指数魅力モデルとかKoyckモデルとかと併用してもよい。

　後半は価格反応の話とモデル例。死ぬほど眠くなってきたので、ばーっとめくって読了ということに。

　やれやれ。著者の先生には悪いけど、これ、素人向けにわかりやすく書こうとしてかえって曖昧模糊となる、という典型例なのではないかと思う。

論文：マーケティング - 読了: Tellis (2006) マーケティング・ミックス・モデリングについて初心者諸君のために解説してあげよう

2019年1月 6日 (日)

　昨年暮れ、学会発表のための分析をやっててわけがわからなくなり、頭を整理するために、Mplusユーザーズ・マニュアル9章(complex survey dataの階層モデリング)の事例40個を全部読んだ。そのときにとったメモ。
　事例9.30-9.40はAR(1)構造がはいって難しさもさらにグレードアップするのだが、当面使う予定がないので(←言い訳)、メモも省略している。

• TYPE=TWOLEVELで、WITHINレベルでもBETWEENレベルでも x が y に刺さっているとき、データとしてxのクラスタ別平均 xm をつくっておいて、 MODELで

  %WITHIN%
  
  y ON X; 
  
  %BETWEEN% 
  
  y ON xm; 

  とする手と、単に

  %WITHIN%
  
  y ON x;
  
  %BETWEEN% 
  
  y on x; 

  とする手がある。前者の場合はDATAでWITHIN=x; BETWEEN=xm;と宣言すること。いっぽう後者の場合、xはWITHINレベルの潜在変数とBETWEENレベルの潜在変数に分解されていることになる。階層回帰モデルの一般的アプローチではないけれど、xmの信頼性が低いときに後者がお勧め。[←これを読んで、むしろ前者がアリだということに驚いた。個票があるのに標本平均をデータにするなんて邪道という気がするが、クラスタのサイズが十分に大きければ構わんってことか...] なお、後者のモデルにおけるBETWEENレベルの係数とWITHINレベルの係数の差を、Raudenbush & Bryk (2002)は「文脈効果」と呼んでいる[←どういう意味だろう？](事例9.1)
• TYPE=RANDOMのときの縦棒("|")の左はランダム回帰係数の名前を表すが、TYPE=TWOLEVEL RANDOMで%WITHIN%での指定だったらそのランダム回帰係数はBETWEENレベルの変数になりますわね。ところがアスタリスクをつけると(「s* | y on x」という風に)、これはWITHINでも分散を持つようになる。[←これはほんとに、いままでぜんっぜん知らなかった。こういうことがあるからMplusは怖い...] (事例9.14)
• TWOLEVEL RANDOMでは、

  %WITHIN% 
  
  logv | y; 
  
  %BETWEEN% 
  
  logv on w;  

  というようなモデルを組める。logvはyの残差分散の対数を表す。[そうそう、これも謎の機能で... ESTIMATOR=BAYESが登場してから追加された機能だと思う。理屈はわかるけど、どういうときに使えばいいのかがいまいちわからない。係数がクラスタによって違うモデルをつくりたいというのはわかるけど、残差分散がクラスタによって違うモデルをつくりたいってことはあるのか...あるか... そうか...] (事例9.28)

　ところで、Mplusは結構長く使っているはずなのだが、それでも「えええ、そんなことができるの」と呆れたモデルがあったので、呆れついでにメモしておく。事例9.26。変数はu(カテゴリカル)、subject, itemの3つ、subjectとitemがクラスタのCROSSCLASSIFIED RANDOMモデル(THREELEVEL RANDOMではないという点に注意)。

%WITHIN%

%BETWEEN subject%

s | f BY u; f@1; u@0; 

%BETWEEN item%

u; [u$1]; s; [s]; 

つまり、subjectレベルではuは分散1の潜在変数fに残差無しで規定されていて、その負荷というかプロビット回帰係数がsである(uをテスト項目の正誤とすると、fは対象者の能力、sは項目の識別性のようなものであろう)。uはitemレベルで切片が動き、残差分散も動く。さらにsも、subjectのなかでは動かないけどitemレベルで動く。ええええ... そういうのを組みたくなる動機はわかるけど、まさかほんとに組めるとは。 sはsubjectレベルで定義されているからitemレベルではモデル化できないかと思った。CROSSCLASSIFIEDモデルおそるべし。

　いやあ、やっぱり使っているソフトのマニュアルは読むべきだ... 読むべきものが多すぎて辛い...
　Mplusについては、Muthen一家の手によるものを含め、多種多様なレベルの解説が手に入るのだが、CROSSCLASSIFIEDモデルについてのまとまった解説は見当たらない。どこかにないかしらん。

2019/01/20追記: 事例9.26のコードについて誤解していたことにあとで気が付いた。CROSSCLASSIFIEDおそるべし。

雑記：データ解析 - 覚え書き: Mplusがご提供する階層モデリングのための謎機能たち

　ここ数日ずっと机に向かってて、もううんざりしてきて、人生どこで間違えた？というあてのない問いを先程から繰り返している状態なので、気分を変えて、Evernoteに溜まっている昨年中の雑多なメモを整理することに。これが終わったら神社に初詣にでも行こう、もう1月も6日だけど...
　これは昨年6月頃に取ったメモ。読み返すと、これ書いてたとき、たしかなにかの締切のようなものから逃避していたなあ、という情緒的イメージだけが蘇ってきて、心臓がバクバクしてきた。

　New York Timesの電子版をぼけーっと眺めていたら、「この25年間の米国演劇を25本選ぶ」というような記事があった。ふと思ったんだけど、この25本のなかで何本くらいが日本で上演されているんだろうか？ 全然知らないけど、アメリカの芝居の翻訳上演のチラシって、結構よく見かけるような気がするんですよね。
　というわけで、google様で検索してはメモをとる、というのを25回繰り返してみた。書籍と違って舞台というのはあっさり消えていくものなので、たぶん見落としがあるだろうと思うけど。

• Topdog/Underdog (Suzan-Lori Parks, 2001) : 2012年に世田谷シアタートラムで上演されている。小川絵梨子演出。
• An Octoroon (Branden Jacobs-Jenkins, 2014) : 見当たらない。
• The Flick (Annie Baker, 2013) : 2016年, 新国立劇場(「フリック」, マキノノゾミ演出)。
• Mr.Burns, a Post Electric Play(Anne Washburn, 2013): 見当たらない。
• Clybourne Park (Bruce Norris, 2010): 見当たらない。
• Ruined (Lynn Nottage, 2008): 見当たらない。
• How I Learned to Drive (Paula Vogel, 1997): 2009年に青年座で、2014年に下北沢(小川絵梨子演出)で上演されている。題名はいずれも「運転免許 私の場合」。これ、面白そう...
• Seven Guitars (August Wilson, 1996): 見当たらない
• Twilight: Los Angeles, 1992 (Anna Deavere Smith, 1994): 見当たらない
• The Designated Mourner (Wallace Shawn, 2000): 見当たらない。ウォレス・ショーンって、なにかの映画にでていた俳優だと思うのだが、劇作家でもあったのか。
• The Humans (Stephen Karam, 2015): 見当たらない。
• This Is Our Youth (Kenneth Lonergan, 1996): 見当たらない。作者は映画「マンチェスター・バイ・ザ・シー」の監督だそうだ。へー。
• Three Tall Women (Edward Albee, 1994): 作者は「バージニア・ウルフなんてこわくない」のオールビー、晩年の作品だそうだ。黒柳徹子さんがパルコのプロデュースで喜劇の舞台を長く続けていたが(1989-2018年)、1996年と2003年の演目がこれだったらしい。題名は「幸せの背くらべ」。
• Jesus Hopped the 'A' Train (Stephen Adly Guirgis, 2000): 見当たらない。
• Eurydice (Sarah Ruhl, 2006): 見当たらない。
• House/Lights (Wooster Group, 1999): 見当たらない。なんかよくわからんが、マルチメディアを駆使したパフォーマンスなのだそうです。
• The Lamamie Project (Moises Kaufman, 2000): 2001年と2003年に燐光群が上演したらしい。
• Yellow Face (David Henry Hwang, 2007): 見当たらない。面白そうなんだけどなあ...
• August: Osage County(Tracy Letts, 2007): 2016年にシアター・コクーンで上演(「8月の家族たち」, ケラリーノ・サンドロヴィッチ演出)。なお、2013年にメリル・ストリープ主演で映画になっている。
• The Vagina Monologues (Eve Ensler, 1996): これは私でも知っている。訳書も出ているし、日本でも何度も上演されている模様。
• Underground Railroad Game (Jennifer Kidwell and Scott R. Sheppard, 2016): 見当たらない。
• The Wolves (Sarah DeLappe, 2016): 見当たらない。
• The Realistic Joneses (Will Eno, 2012): 見当たらない。
• The Apple Family Plays (Richard Nelson, 2010): 見当たらない。
• The Elaborate Entrance of Chat Deity (Kristoffer Diaz, 2010): 見当たらない。

というわけで、日本での翻訳上演が確認できたのは25本中7本であった。多いような、少ないような。

雑記 - NYTが選ぶところの過去25年間のアメリカ演劇25

2019年1月 4日 (金)

ザ・空気 [a]
永井 愛 / 而立書房 / 2018-07-25
2017年上演、大変評判になった戯曲。

フィクション - 読了：「ザ・空気」

市場の倫理 統治の倫理 (ちくま学芸文庫) [a]
ジェイン ジェイコブズ / 筑摩書房 / 2016-02-09
この秋に読んでいた本。2018年に出会った大当たりの一冊であった。原題は"Systems of Survival: A Dialogue of the Moral Foundations of Commerce and Politics"。「倫理的基盤」といわれるとちょっと引いてしまうけれど、社会の矛盾をふたつの概念的フレームワークの対立という観点から描き出す知的エンターテイメント、という感じの内容。よく考えてみると厳密な議論じゃないような気もするんだけど、面白い本であった。

徹底検証 神社本庁 (ちくま新書) [a]
明, 藤生 / 筑摩書房 / 2018-10-05

海賊の日本史 (講談社現代新書) [a]
山内 譲 / 講談社 / 2018-06-21

異端の時代――正統のかたちを求めて (岩波新書) [a]
森本 あんり / 岩波書店 / 2018-08-22

マルクスと批判者群像 (平凡社ライブラリー) [a]
良知 力 / 平凡社 / 2009-02-01

現代社会はどこに向かうか－高原の見晴らしを切り開くこと (岩波新書) [a]
見田 宗介 / 岩波書店 / 2018-06-20

ルポ 不法移民とトランプの闘い 1100万人が潜む見えないアメリカ (光文社新書) [a]
田原徳容 / 光文社 / 2018-10-16

ノンフィクション(2018-) - 読了：「市場の倫理 統治の倫理」「神社本庁」「異端の時代」「マルクスと批判者群像」「現代社会はどこに向かうか」「不法移民とトランプの闘い」「海賊の日本史」

往生要集―日本浄土教の夜明け (2) (東洋文庫 (21)) [a]
源信 / 平凡社 / 1964-06-01

新版 歎異抄 現代語訳付き (角川ソフィア文庫) [a]
/ 角川書店 / 2001-07-25

親鸞・普遍への道―中世の真実 (ちくま学芸文庫) [a]
阿満 利麿 / 筑摩書房 / 2007-04-01

幸福とは何か - ソクラテスからアラン、ラッセルまで (中公新書) [a]
長谷川 宏 / 中央公論新社 / 2018-06-20

チベット仏教入門 (ちくま新書) [a]
均, 吉村 / 筑摩書房 / 2018-11-06

哲学・思想(2011-) - 読了：「チベット仏教入門」「往生要集」「歎異抄」「親鸞・普遍への道」「幸福とは何か」

FEAR 恐怖の男 トランプ政権の真実 [a]
ボブ・ウッドワード / 日本経済新聞出版社 / 2018-12-01

いいビルの世界 東京ハンサムイースト [a]
東京ビルさんぽ / 大福書林 / 2017-10-09

戦略の世界史(上) 戦争・政治・ビジネス [a]
ローレンス・フリードマン / 日本経済新聞出版 / 2018-09-26

運慶のまなざし――宗教彫刻のかたちと霊性 [a]
金子 啓明 / 岩波書店 / 2017-11-08

漂流児童 [a]
石井 光太 / 潮出版社 / 2018-10-05

ノンフィクション(2018-) - 読了：「いいビルの世界」「運慶のまなざし」「漂流児童」「戦略の世界史」「恐怖の男」

さよならミニスカート 1 (りぼんマスコットコミックス) [a]
牧野 あおい / 集英社 / 2018-11-22
講談社の少女誌「りぼん」が大見得を切って開始した連載。どんなものかと手に取ってみて、なるほどなあ、と納得...

はじめてのひと 3 (マーガレットコミックス) [a]
谷川 史子 / 集英社 / 2018-09-25
名匠・谷川史子の漫画であるからして、素敵な恋愛話になってますけど、でもこれ不倫ですよね！この16歳年上のチェリスト、こいつ要はクズですよね！とひとりで憤りながら読んでいた。おっさんはね！こういう男には厳しいよ！！

ゆりあ先生の赤い糸（１） (BE・LOVEコミックス) [a]
入江喜和 / 講談社 / 2018-07-13
ゆりあ先生の赤い糸(2) (BE LOVE KC) [a]
入江 喜和 / 講談社 / 2018-10-12
私が尊敬してやまないマンガ家・入江喜和さんの新作。前作「たそがれたかこ」もそうだったけど、ちょっと類例が思い浮かばない作品なので、果たして読者がついてくるだろうか...打ち切りになりませんように...と勝手に気を揉んでいる。

コミックス(2015-) - 読了：「さよならミニスカート」「はじめてのひと」「ゆりあ先生の赤い糸」

ダルちゃん: 1 (1) (コミックス単行本) [a]
はるな 檸檬 / 小学館 / 2018-12-06
ダルちゃん (2) (コミックス単行本) [a]
はるな 檸檬 / 小学館 / 2018-12-06
資生堂「花椿」web版での連載を固唾を飲んで追いかけていたマンガが無事終了、このたび書籍になった。現代社会で若い女性が抱える疎外感を鮮烈に描いた作品。正直いって最後の展開は拍子抜けというか、ちょっと納得がいかないところがあったんだけど(私たちが抱えている葛藤と軋轢はそれほど簡単には解決しないのではないかと思う)、それはそれとして、とても面白い物語であった。

ランド(7) (モーニング KC) [a]
山下 和美 / 講談社 / 2018-09-21

神は細部に宿るのよ(5) (ワイドKC) [a]
久世 番子 / 講談社 / 2018-10-12

邦画プレゼン女子高生 邦キチ! 映子さん (ホーム社書籍扱コミックス) [a]
服部 昇大 / ホーム社 / 2018-08-24

きょうのカプセル (ワイドKC) [a]
黒田 硫黄 / 講談社 / 2018-11-22

コミックス(2015-) - 読了：「ダルちゃん」「ランド」「神は細部に宿るのよ」「邦キチ！映子さん」「きょうのカプセル」

銃座のウルナ 6 (ビームコミックス) [a]
伊図透 / KADOKAWA / 2018-09-12

ふしぎの国のバード 5巻 (ハルタコミックス) [a]
佐々 大河 / KADOKAWA / 2018-10-15

二月の勝者 ー絶対合格の教室ー (3) (ビッグコミックス) [a]
高瀬 志帆 / 小学館 / 2018-10-12

にわにはににん (ビームコミックス) [a]
中野 シズカ / KADOKAWA / 2018-09-12

独身OLのすべて(8) (KCデラックス) [a]
まずりん / 講談社 / 2018-09-21

赤狩り THE RED RAT IN HOLLYWOOD (3) (ビッグコミックス) [a]
山本 おさむ / 小学館 / 2018-09-28

コミックス(2015-) - 読了：「赤狩り」「独身OLのすべて」「二月の勝者」「にわにはににん」「銃座のウルナ」「不思議の国のバード」

あげくの果てのカノン 1 (ビッグコミックス) [a]
米代 恭 / 小学館 / 2016-06-10

小路花唄(3) (アフタヌーンKC) [a]
麻生 みこと / 講談社 / 2018-10-05

レイリ(5)(少年チャンピオン・コミックス・エクストラ) [a]
岩明 均,室井 大資 / 秋田書店 / 2018-11-08

孤食ロボット 5 (ヤングジャンプコミックス) [a]
岩岡 ヒサエ / 集英社 / 2018-10-25

めしにしましょう(6) (イブニングKC) [a]
小林 銅蟲 / 講談社 / 2018-10-23

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

コミックス(2015-) - 読了：「BLUE GIANT SUPREME」「あげくの果てのカノン」「小路花唄」「レイリ」「孤食ロボット」「めしにしましょう」

私の少年（５） (ヤングマガジンコミックス) [a]
高野ひと深 / 講談社 / 2018-11-06

闇金ウシジマくん (44) (ビッグコミックス) [a]
真鍋 昌平 / 小学館 / 2018-10-30

大奥 16 (ヤングアニマルコミックス) [a]
よしながふみ / 白泉社 / 2018-10-29

海街diary 9 行ってくる (フラワーコミックス) [a]
吉田 秋生 / 小学館 / 2018-12-10

猫のお寺の知恩さん (8) (ビッグコミックス) [a]
オジロ マコト / 小学館 / 2018-09-28
猫のお寺の知恩さん (9) (ビッグ コミックス) [a]
オジロ マコト / 小学館 / 2018-12-27

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

コミックス(2015-) - 読了：「めしばな刑事タチバナ」「猫のお寺の知恩さん」「海街diary」「大奥」「闇金ウシジマくん」「私の少年」

大砲とスタンプ(8) (モーニング KC) [a]
速水 螺旋人 / 講談社 / 2018-12-21
ウクライナあたりをモデルにした架空の国の「兵站軍」を描く。分野としてはミリタリーマンガということになると思うのだけど、そういうのに全く関心がない私にとっても、ドラマとして面白い。演出はコミカルだけど、よく考えてみるとシビアな物語なのだ。いよいよ大詰めとなる模様...

血の轍 (4) (ビッグコミックス) [a]
押見 修造 / 小学館 / 2018-09-28

７人のシェイクスピア　ＮＯＮ　ＳＡＮＺ　ＤＲＯＩＣＴ（７） (ヤングマガジンコミックス) [a]
ハロルド作石 / 講談社 / 2018-12-06

重版出来! (12) (ビッグコミックス) [a]
松田 奈緒子 / 小学館 / 2018-12-12

ゴールデンゴールド(5) (モーニング KC) [a]
堀尾 省太 / 講談社 / 2018-12-21

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

コミックス(2015-) - 読了：「血の轍」「大砲とスタンプ」「７人のシェイクスピア」「重版出来！」「ゴールデンゴールド」「レベレーション」

2018年中に読んだ本を記録しておく。まずはマンガから。

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

オリンピア・キュクロス 2 (ヤングジャンプコミックス) [a]
ヤマザキ マリ / 集英社 / 2018-12-19

3月のライオン 14 (ヤングアニマルコミックス) [a]
羽海野チカ / 白泉社 / 2018-12-21

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

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

ブラック・ラグーン (11) (サンデーGXコミックス) [a]
広江 礼威 / 小学館 / 2018-11-19

聖☆おにいさん(16) (モーニング KC) [a]
中村 光 / 講談社 / 2018-11-22

コミックス(2015-) - 読了：「聖☆おにいさん」「BLACK LAGOON」「乙嫁語り」「オリンピア・キュクロス」「３月のライオン」「イノサン」「あさひなぐ」

　年末に仕事の都合で、大急ぎで読んだ奴。他にやんなきゃいけないことがいっぱいあるこのときに、調査データの標本ウェイティングなどという話を...

DeVell, M., Krosnick, J. (2010) Computing Weights for American National Election Study Survey Data. ANES Technical Report series, no. nes012427.

　複雑な標本設計の調査において標本にウェイトを付与するためにどうするか、という話を調べていると、時々みかけるのがANES(American National Election Studies)である。これはそのテクニカルペーパーで、Rのanesrakeパッケージの元文献でもある。
　今気がついたけど、第二著者は調査法研究のKrosnickだ... スタンフォードの政治学者だったのか。てっきりミシガンの社会調査の人だと思ってた。

　ちょっと意外だったんだけど、このペーパーは「ANESの公開データに付与されている標本ウェイトを我々がどうつくったか説明する」という話というより、「ANESデータを分析するみなさんのために、標本ウェイトの作り方のガイドラインを示す」という感じの内容である。著者らいわく、正しいウェイトというのはひとつではない、ウェイトは分析目的によっても異なる、とのこと。
　(突き詰めていえば仰るとおりなんだけど、でもふつう標本ウェイティングを使うのは、データ分析の目的を問うことなく、標本抽出確率の不均等性をデータ収集デザインの問題として捉えようとしているときだと思う。もし個別具体的なデータ分析目的の下で考えるんなら、標本抽出確率に影響する変数を説明変数にいれたモデルを組むんじゃなかろうかと思うんだけど、そんなことないですかね？)

　業界の人を読み手として想定した愛想のない文章で、たとえばレイキングという概念についての事前説明はない。
　えーと、レイキングとは次のような話である。たとえば、20-60代男女の調査データにおいて、男女の分布が6:4, 年代の分布が2:2:1:1:1となっている。いっぽう母集団における男女の分布は1:1、年代の分布は1:1:1:1:1である、としよう。このまま集計しちゃうと歪んじゃうので、男や若者には小さめのウェイトを与えたい。さて、ウェイトを計算するにはどうするか。通常は、性x年代の10セルについて人数を数え、母集団における分布と比べてウェイトを決める。これがセル・ウェイティング。ところがこのやり方では、母集団における性x年代の分布が必要になるし、変数がこの2つじゃなくて3つ, 4つ...と増えていったときに大変な目に遭う。そこで、まず性に適切なウェイト値を振って、集計後の男女の分布を1:1にする。次に年代に適切なウェイト値を振って、集計後の年代の分布を1:1:1:1:1にする。すると、集計後の男女の分布が1:1ではなくなるので、ウェイト値をちょっと修正して、再び1:1にする。すると集計後の年代の分布が1:1:1:1:1でなくなるので...という作業を繰り返す。これがiterative proportional fitting, 愛称rakingである。
　以前、お取引様のご担当者(実は前職の同僚)とお電話でご相談していて、上記のレイキングについて説明したところ、それって掻き集めるってことですかね？という。「英語でrakeっていうんですよ、落葉とかをこう掻き集めるの」(彼は米国暮らしが長かったのである)。「あー！そういえばホームセンターでクマデ探してたら、店員さんにレーキですねっていわれた！」「あー！」というわけで、電話口でおっさん二名、ジャニーズ並みの美しいユニゾンとなった。後で調べたところ、rakingは熊手で縦に横に掻き集めるというところから由来している呼び方らしい。

　本文では、単一の横断調査、2ウェーブ・パネル、多ウェーブ・パネルの3つについてウェイティングの手順を説明しているんだけど、今関心があるのは横断調査なので、そこだけメモする。面接調査について説明するけど、RDDでも同様、とのこと。
　箇条書きで21項目もある...さあ深呼吸...

1. まず、層別抽出とかオーバーサンプリングとかに由来する世帯選択確率の不等性に対処するためのウェイトをつくる。各単位がそれぞれの選択段階[←多段抽出のことね]で選択される確率をかけ算し、その逆数をウェイトにする。
2. 世帯内の対象者選択確率の不等性に対処するウェイトをつくる。対象者を世帯内でランダムに選んでいる場合なら、世帯の人数をweighting factorにする。
3. 無回答率の不等性に対処するウェイトをつくる。無回答者についてわかるのはふつうクラスタとか地域とかだけなので、そういう変数の水準ごとに回答率をみて、違っていたら、回答率の逆数をweighting factorにする。
4. ここからは、既知の人口パラメータで事後層別するという話。まずベンチマーク比較をする。測定誤差と項目無回答率が小さいと思われる変数群を決め、その分布(もちろん、1-3で求めたウェイトを使う)を、Current Population Survey(CPS)みたいなベンチマークと比べる。年齢、性、人種、学歴、地域、等々がよろしかろう。選挙後調査なら、得票の分布、投票率も使うべし。投票率のベンチマークはvoting-age人口における投票率じゃなくてvoting-eligible人口における投票率にすべし[ふうん...なんか事情があるんでしょうね]。
5. 上の結果を解釈する。もしずれが小さかったら気にしない、大きかったらその変数で事後層別する。ずれが5%ポイントを超えたら後者、2%ポイントを下回ったら前者とすることが多い。その間なら...たとえば得票率や投票率のような重要なアウトカムと関連している変数とか、それ自体に関心がある変数とかだったら、ウェイティングした方がいいんじゃないですかね。
6. とこのようにして、事後層別に用いる変数を選びます。ただし、
   • 測定が正確で、無回答率が低い変数でないとだめ。特に年収は要注意、項目回答率が95%を切ってたらやめたほうがいい。
   • 得票率と投票率をいれる場合は、いれないバージョンといれるバージョンをつくるべし。
   • 事後層別ウェイティングは一元周辺分布、ないし二元周辺分布で行う[つまり、年代の分布や性x年代の分布はベンチマークに揃えようとするけど、性x年代x人種の分布を揃えようとはしない]。二元周辺分布で行うのは、一元周辺分布のレイキング(後述)ではうまく揃わない場合、ないし、下位グループでの分析に最適化したウェイトをつくりたい場合。
   • 欠損値がある場合、もし単純な補完方法がすぐに使えて、それが良い性質を持っていることがわかっているならばそれで補完すればよい。ふつうはそんな方法はないだろうから、レイキングの際には欠損があるケースもそのまま含めておき、最終的には無回答調整のためのベースウェイトを与える。[うーん、確信が持てないけど、たぶんこういうことではないかと思う。たとえば、性x年代の事後層別ウェイトをレイキングで求めたい、しかし年代には欠損がある、としよう。レイキングでは、年代が欠損である行を、年代の分布をベンチマークと揃えるステップでは無視し、性の分布をベンチマークと揃えるステップでは無視しない。このやりかたでもすべての行になんらかのウェイトが付与されることになる。で、年代が欠損だった行からはそのウェイトを取っ払い、ステップ1-3で求めたウェイトに戻す。あれれ？ ということは、性についても年代についても、最終的な分布はベンチマークとちょっとずれることになるわけ？]
   • 連続変数、6カテゴリ以上のカテゴリカル変数、5%を切るカテゴリを持つカテゴリカル変数の場合はリコード。5カテゴリまで、いずれも5%以上とする。[あああ...ここはちょっと耳が痛い...]
7. 事後層別ウェイティングの際には、セルウェイティングじゃなくてレイキングを使うこと。
8. レイキングのためのコードはスクラッチから書いてもいいが、StataのWESVARとかSURVWGTモジュール、RのSurveyパッケージ, SASのRAKINGマクロを使ってもよい。
9. レイキングの手順。
   • 無回答調整済ベース・ウェイトに値を掛けて、最初の事後層別ファクターの母集団分布に合致させる。たとえば、標本では女性が60パーセント、母集団では52パーセントだったら、ウェイトは0.52/0.60=0.87。
   • さらに値を掛けて、二番目の事後層別ファクターの母集団分布に合致させる。
   • これを残りのファクターについて繰り返す。
   • それぞれのステップで、極端なウェイトはcapする(truncateする)。ウェイトの平均を1として、5を超える奴を5にする。capした行にはその旨のフラグを振っておく。なお、上はこうしてcapするけど、下はcapしなくてよろしい。
   • 最初のファクターに戻って繰り返す。ベンチマークに一致するか、変化がなくなるまで繰り返す。
10. 修正してもう一度レイキングする。
    • capのフラグが特定の特性を持つ行に偏っている場合は、capの閾値を上げてやり直す。
    • レイキングが終わったら、つくったウェイトの下でベンチマークと比較する
    • ベンチマークとずれていたら、capの閾値を上げてやりなおす。
    • レイキングのファクターになっていない変数について、レイキングの影響を調べる。ベンチマークとのずれが大きくなっているなら、後述する修正アプローチを試す。
    • ウェイトの変動係数と、新しいウェイトの下でのデザイン効果を調べる。レイキング後のデザイン効果がレイキング前のデザイン効果を(たとえば)0.5以上超えていたら、後述する修正アプローチを試す。[デザイン効果についての説明がないので困ってしまうが、おそらくはKishの定義、「標本の分散の、同じ要素数の単純無作為標本の分散に対する比」でよいのだろう]
    • 何か特定のマイノリティ集団について調べたいのなら、そこでの変動係数も調べておく。
    • 修正アプローチ: レイキングに使うカテゴリを併合する、事後層別するファクターを削るか入れ替える、capを調整する。
11. 最後に、ウェイトの平均を1にする。
12. 標本全体をベンチマークと一致させるウェイトが、標本を特定の集団に絞ったときにそれをベンチマークと一致させるウェイトになるとは限らない。もし下位集団の分析が研究目的なら、その下位集団のためのウェイトをつくるべきである。手っ取り早いのは、最初から二元周辺分布に対してレイキングしておくことである(たとえば、あとで男性だけの分析をするときにも年代の分布をベンチマークに合わせたい場合は、あらかじめ性x年代の二元周辺分布に対してレイキングしておく)。
13. うまくウェイトがつくれたら、ウェイトはふつう正確性を増大させるが、分散は大きくなる。特定のデモグラフィック・ファクターについての正確性を最大化したい場合や、頑張ってウェイトをつくってみたはいいものの結局ベンチマークとのずれが残っているという場合には、デザイン効果が大きくなっても気にしないウェイトをつくり、これをsupplimentary weightとする。つまり、capは高くし、カテゴリは細かくし、なんなら事後層別ファクターを増やしてレイキングする。[この場合、分析報告では2種類のウェイトについての集計を報告することになるのだろう。はっはっは、めんどくせえなあ]
14. ここからはドキュメンテーションの話。データには、層別・クラスタリングの全水準がわかるように、層とかPUSとか地域とかの変数を含めること。地域は恣意的な番号でよい(その地域が実際にどの地域かがわからない形でよい)。
15. レイキング前のベースウェイトもデータに入れておくこと。
16. ウェイトをつくる各段階で何を使ったかを述べること。
17. 複数種類のウェイトをデータに含める場合には、どのウェイトを使うべきか素人にもわかるように説明すること。
18. 目標母集団の特性とサイズ、そしてレイキングで使ったすべてのベンチマークについて記述すること。ウェイトの合計を母集団サイズに揃えたい場合にはウェイト値を何倍すれば良いかを述べること。
19. レイキングに使った変数を列挙し、それらの変数においては標本抽出誤差が意味を持たないということを説明すること。[あ、なるほど...これは当然だけど大事だ]
20. 関心ある統計量に対するデザイン効果を報告すること。
21. ウェイトをtruncateした場合はその方法を正確に述べること。

　...なにか目の覚めるような斬新なアドバイスがあったりしないかと思ったけど、そういうものではなかった。
　まあとにかく、ウェイト値のトリミング(ここでいうcap)に関して明確な基準がないということがわかった。なにかの目的関数を最適化するような閾値を決めるというアプローチがありそうなものだけれど、レイキングとは相性が悪いのかもしれない。

論文：データ解析(2018-) - 読了：Devell & Krosnick (2010) American National Election Studyにおける標本ウェイティング

« 2018年11月 | メイン | 2019年4月 »