井出草平の研究ノート

暴力的なゲームと攻撃性を結びつける研究に不正が見つかり撤回される

ネット・ゲーム依存

山根信二さんのTwitter経由で知った情報(https://twitter.com/shinjiyamane/status/1362068833641504776)。

www.sciencemag.org

今回、問題になっているのは重慶南西大学の心理学者である张谦Qian Zhang氏。张氏は暴力的ゲームをすると暴力性が増加するといった研究や、暴力的な映画を観ると暴力性が増すといった論文を次々に発表しているが、データが捏造されたものだと指摘され、张氏の論文のうち2本は撤回されている。

どうやら、調査をせずに架空のデータで、分析も捏造し、論文を書いたようである。

ただ、撤回されたのは2本だけで、今年に入ってからも张氏は精力的に論文を発表している。

ざっくりと翻訳してみた。

张氏は不正行為を否定しているが、2つの論文は撤回された。他の論文は雑誌にはまだ残っており、そのデータはメタアナリシスに含まれた形になっていることが「大きな問題」になるとメディアと行動を研究するケンブリッジ大学の認知科学者、エイミー・オーベンAmy Orben氏は述べている。この不正の影響は、学術の中に留まらないという。张氏が不正をした研究は、メディアの警告ラベルであったり、親や医療専門家がゲームをどのように子ども与えるかという根拠になってしまっているのだ。

調査のきっかけとなったのは、イリノイ州立大学の心理学者ジョー・ヒルガードJoe Hilgard氏で、先月、张氏の研究についての懸念点についてブログ記事にして公開した(http://crystalprisonzone.blogspot.com/2021/01/i-tried-to-report-scientific-misconduct.html)。ヒルガードは、彼が3000人の被験者を対象とした別の研究である"Youth & Society"の2018年の論文(https://journals.sagepub.com/doi/abs/10.1177/0044118X18770309)を初めてみた時感銘を受けたという。「何てこった！」と思いましたよ、と彼は言う。この研究では、一部のティーンエイジャーが暴力的なビデオゲームをプレイした後、より攻撃的になっていると書かれていたのだ。膨大なサンプルサイズを考えると、この研究は「強力なエビデンス」になる可能性があったとヒルガード氏は言う。

しかし、ヒルガード氏はこの論文に統計的な誤り、数学的に不可能な部分があることを発見した。张氏と彼の共著者たちは、高い統計的有意差を報告しているが、実際には統計的有意差のない小さな差であった。ヒルガード氏は张氏とジャーナルに警告を発し、张氏は修正を提出し、より信憑性のある論文になったかのように見えたが、それでも数値をごまかしていた。

ヒルガード氏は、3つの異なる論文でほぼ同じ結果が報告されているなど、张氏の他の論文にも問題があることを発見した。彼は张氏に電子メールを送り、彼のデータを見るように頼んだが、Zhangは拒否したという。ヒルガード氏はその後、ノースカロライナ大学チャペルヒル校の心理学者であり、複数の論文で张氏との共著者であるドロシー・エスペラージュDorothy Espelage氏に連絡を取った。ドロシーも张氏にデータを送るように言ったが、拒否されたとのことだった。ヒルガード氏が重慶南西大学に調査を依頼した後に、张氏が映画の暴力に関する"Youth & Society"の論文(https://journals.sagepub.com/doi/10.1177/0044118X18775846)のデータをヒルガード氏に送ってきただけだった。

しかし、そのデータも奇妙なもので、その論文で行われた実験と同じような実験でみられる特徴が欠けていたとヒルガード氏は言う。ヒルガード氏は調査結果を张氏の大学に送ったが、张氏のデータには問題はあったが不正ではなかったとし、张氏には "統計的な知識と研究方法に欠陥がある "と結論付けた。この返答に不満を感じたヒルガード氏は、関係するすべての雑誌に自分の観察結果を送った。

ミシガン大学アナーバー校の心理学者であり、"Youth & Society"の編集者でもあるマーク・ジマーマンMarc Zimmerman氏によると、ヒルガード氏のメールは、別の情報源が张氏のデータに疑問を呈した1週間後に来たという。2カ月後の2019年12月、ジャーナルは张氏の2つの論文を撤回した。2つの告発がなければ、迅速で、かつ、撤回という形にならなかったかもしれないと言う。「間違ったデータを公表したいと思うジャーナル編集者おらず、深刻な告発だった」とジマーマン氏は言う。

张氏はScience誌へ宛ての電子メールで、不正行為を否定し、「データを見るための『正しい』方法は1つだけではない」と書いている。张氏は、ヒルガード氏が暴力的なゲームと攻撃性との関連性に懐疑的であり、彼が他の人の研究に難癖をつけることで、有名になろうとしていると書いている。しかし、ダートマス大学の心理学者であるジェイ・ハルJay Hull氏(暴力的なビデオゲームと攻撃性の間に関連があるという立場である)は、ヒルガード氏を支持すると述べている。ハル氏は科学は不正やいんちきに基づいてはならないと述べている。

张氏と共著者であったドロシー・エスペラージュは"Aggressive Behavior"に掲載された論文の撤回を支持しているが、"Aggressive Behavior"を含む他のジャーナルは、张氏の研究をいまだ撤回していない。アイオワ州立大学の心理学者であり、"Aggressive Behavior"の編集長であるクレイグ・アンダーソンCraig Anderson氏はコメントを辞退した。

张氏の最近の論文は、ヒルガード氏の指摘を避けて書かれているが、ヒルガード氏は満足していない。こういった論文の不正への対策が遅れると「不正が巧妙になり、発見できなくなる恐れがある」とヒルガード氏は述べている。

メンタルドクターSidowさんのツイートの問題点

中谷さんから話を振られたので、少し詳しく書いてみようと思う。

メンタルドクターSidowさんのツイート

f:id:iDES:20210218000042p:plain https://twitter.com/dr_sidow/status/1360527222423998464

彼の認識が特殊だというわけではなく、一般的な精神科医の認識もこのような感じではないかと思う。
ただ、一般的な認識が正しいとは限らないこともあり、このツイートには大きな誤りがある。

ゲーム障害＝精神疾患

世界保健機関WHOのICD-11が発効するのは2022年からなので、厳密にいうと、ゲーム障害(ゲーム症)は現時点ではまだ精神疾患ではない。ただ、ほぼ間違いなく、来年からはICD-11が採択されると思われるので、現時点で精神疾患だと言っていたとしても、特に問題はないだろう。

精神疾患は病気を意味するのか

なんじょうさん(@nnanjoh)がICDの話もされているので用語回りを整理しておこう。
結論から言うと、精神疾患は病気や疾病といった意味を持っていない。

ICDを含め、公式に、英語・日本語とも述語が混乱しているので、このあたりの用語かなりわかりにくい。
専門家も混乱した述語で困っているため、WHOに直接ICD-11での述語の使い方を問いただしたが、回答は、矛盾なく統一しようと試みたがICDのカバーする範囲が広いため、断念した(大意)というものだった。
また、アメリカ精神医学会のDSMの翻訳はDSM-IIIの時から精緻だが、ICDのF項目(精神障害)の翻訳は間違いや矛盾が散見される。

このようにICD原版も混乱があり、日本語訳にも混乱があるので、二重苦なのだが、一応、下記のように捉えると、間違いはない。

  • illness: 病気
  • disease: 疾病
  • disorder: 障害

疾患は対応する英語がない。日本では疾病も含んだ精神障害なども含意する用語である。

精神医学ではMental Disorderを「精神障害」と翻訳する時もあれば「精神疾患」と翻訳する時もある。DSMの本文では精神障害されているが、タイトルでは精神疾患が使用されている。このあたりも混乱を招いている原因の一つである。DSMの翻訳者によれば、精神障害という用語が日本語としてわかりにくいため、タイトルは精神疾患という表現にしたようなのだが、余計に混乱を来した感じはある。

精神障害＝精神疾患＝疾病・病気とは異なる概念

と捉えておけば間違いはない。

ちなみに、ICD-10では比較的登場回数の多かったillnessという用語はICD-11ではおおむね削除することに成功したようである。

過度にゲームを続けると脳や行動に悪影響が及ぶことは研究レベルで証明されています。

ここが問題の部分である。
結論から言うと証明されていない。日本語の文献や久里浜の樋口進さんの本などを読んでいると、そう誤解しても仕方ないが、脳画像研究の論文を読んでいると、この認識が誤りであることがわかるはずである。
現在の脳画像研究の考え方と、その研究がどういう経緯をたどってきたかを簡単ではあるが、紹介しておこう。

脳の仕組み

脳の神経組織では白質と灰白質という用語が使われる。
脳神経細胞そのものは灰白質、その伝達になっているのが白質だと理解すると早い。
白質の画像研究では拡散テンソルイメージング（DTI）法を用いる。
灰白質の画像研究ではfMRIを用い、特にボクセルベースの形態計測法（VBM）が用いられている。

黎明期

白質の研究は1990年代末に始まり、本格的に論文が出るのでは2000年代に入ってからである。灰白質の研究は少し遅れ2000年半ばから後半くらいから始まる。白質の研究の方が少し研究され、それを追う形で灰白質の研究は進んでいったというのがざっくりとした研究の流れである。

今回は神経細胞の話なので、灰白質の研究に限定しよう。
本格的にゲーム障害と灰白質の研究の中で最初にまとまったものが出版されたのは2011年の下記の研究Yuan et al.(2011)だと思う。

journals.plos.org

このブログでも(解説記事)をエントリした論文だ。確かに、この論文には、ゲーム障害＝ゲームのやりすぎは脳にダメージ(特定の部位の脳の容積が低下する)と示唆されるいう趣旨のことが書かれている。久里浜の樋口さんなどの論調はこの2010年代前半の論文によく登場する考え方である。

研究の蓄積による論調の変化

灰白質の画像研究は新しい領域の研究であったため、2000年代、2010年代を通して、多く研究がされた。もちろん、さまざまな分野・疾患で実施された。研究が積み重なっていくうちに、次第に論調が変化してきた。

ゲーム障害はADHDとの併存(ADHDが先んじる)ことが多いが、ADHDとゲーム障害で問題となる箇所が酷似していることが判明している。(解説記事)

pubmed.ncbi.nlm.nih.gov

精神病性障害を除いた精神障害の灰白質の減少箇所は共通していることも判明している。(解説記事)

www.ncbi.nlm.nih.gov

共通してみられるのは前帯状皮質(ACC)と背外側前頭前野(DLPFC)と島(insula)であるが、精神障害によって特定の分野の脳の容積が少ないといった対応関係が同定はされていない。

では、もともとADHDであり、その後ゲーム障害になった児童に脳画像研究をした場合、そこで写しているのは、ADHDの画像なのか、ゲーム障害の画像なのかという疑問が当然浮かぶ。
この疑問を解決するには、多変量解析が必要だが、fMRIで脳画像を撮影する値段が馬鹿高いためサンプル数が集められず、今のところ実現していない。できれば、コホートでゲーム障害が生じる前から実施すると一番理想的だが、こちらも資金面から困難である。

現在は以下の3つの可能性が示唆されている。

  1. ゲームが神経細胞を破壊した可能性
  2. 神経細胞が少ない人がゲーム障害になりやすい可能性
  3. ゲーム障害と併存するADHDやうつ病によって神経細胞が破壊された可能性

現在では2を指示するエビデンスが次第に増えている。例えば、2019年のアルコール使用障害の研究では、灰白質の体積の少なさが原因かもしれないと指摘されるようになっている。(解説記事)

pubmed.ncbi.nlm.nih.gov

2010年代前半の研究は、精神障害が原因で特定部位の脳の容積が低下しているのではなと推論される傾向があったが、現在では逆転した理解がされるのが標準的である。
つまり、生まれつき、もしくは成長の過程で、脳の特定部位の容積が少ないことが、ゲーム障害を含め、他の精神障害のリスクになるのではないか、というものだ。

最近のゲーム障害の脳画像研究

久里浜の樋口さんもよく引用する、ゲーム障害のボクセルベースの形態計測法のレビュー、メタアナリシス論文がある(Yao et al. 2017)。執筆年は2017年である。(解説記事)

pubmed.ncbi.nlm.nih.gov

この論文を読めばわかるが、ゲームのやりすぎによって脳の特定部位の容積低下が起こったなどとは一言も書かれていない。
この論文の結論は1)ゲーム障害では、行動面のみならず脳における神経変異があるということ、2)両側前帯状皮質がバイオマーカーになるのではないかというものだ。

ただ、この論文は日本では、久里浜・樋口式のアレンジが加わり、2010年代前半に議論されていた、ゲームのやりすぎによって脳に影響が生じると喧伝されるされているのが現状である。

時代の変化

2011年の論文(Yuan et al. 2011)と2017年の論文(Yao et al. 2017)を実際に読み比べれば、時代の流れはすぐに把握できる。Yuan et al. (2011)は久里浜の樋口さんやゲーム害悪派の主張に割と近い論調であるが、Yao et al.(2017)は全く異なっている*1。

現在では、Yao et al.(2017)のゲームの論文だけではなく、他の精神障害の論文でも、精神障害が原因で脳の特定部位の容積が低下したと書く論文はない。

脳の特定部位の容積の低下は精神障害の結果であるという見解は2010年代初頭にあった考え方で、最近は、1)原因、もしくは2)精神障害を見つけるバイオマーカーという2つの論調で論文が書かれているのが普通である。

メンタルドクターSidowさんの問題

シンプルに勉強不足である。また、精神科医という肩書は専門家とみなされるので、発言をするときには、しっかりと英文で最新情報を押さえたうえで、発言するべきである、ということだろうか。

Twitterの名前には脳のアイコンもついており、脳にこだわりも持っているようなので、是非、脳画像研究の論文もしっかりと読んでいってもらいたいものである。

*1:2011年の論文(Yuan et al. 2011)でも「ゲームのやりすぎで脳に悪影響がある」と言い切っているわけではない。

鹿児島県警サイバー犯罪対策課のツイートはどこが問題か

ネット・ゲーム依存

鹿児島県警サイバー犯罪対策課がゲーム障害についてツイートして、批判が寄せられたため、取り下げるということが起こったらしい

news.livedoor.com

Twitterのみなさんのスピードには追い付かないが、僕の立場からできる指摘を行っておこう。 当のツイートは以下のようなものだった。

f:id:iDES:20210215234114p:plain

スマホでゲームをする時間を自分でコントロールできますか？

おそらく、ゲームをやり始めて、ついつい当初の時間より長くやってしまった、という意味だろう。使用時間が増えることが、ゲーム障害やネット依存症の証拠だと考えている人もいるようだ。しかし、それは間違いである。

単に「ゲームをする人」にとどまらない「ゲーマー」と目され、ゲームマニアと呼ばれることもある熱狂的なプレイヤーたちの一群も存在する。これらの人々は、週に20-30時間をゲームに費やしていながら、IGD(インターネット・ゲーム障害)の基準を満たしていない。他のゲーマーとの競技に参加し、観衆を楽しませて給与を得るプロゲーマーたちのリーグもあるが、こうした人々もIGD基準に適合しないことが多い。(ダニエル・キング『ゲーム障害 ゲーム依存の理解と治療・予防』:60)

依存と使用時間には関連があるはずだ、と考える人は多いが、現実はそうではない。例えば以下のような研究がある。

ネット依存症とネット使用時間は無関連。関連があるのは不安症状。

ides.hatenablog.com

私たちは行動をコントロールできているのか？

ゲームをついついやりすぎてしまうことは異常なのことなのだろうか。おそらく「よくあること」である。ついついごはんを食べ過ぎてしまったり、ついつい観るつもりのなかったテレビ番組をみてしまったり、するものだ。

私たちは自身を完全にコントロールできているはずがない。まぁ、もちろん、中には、当初の予定時間通りに「まるで機械のように」動くことのできる人もいるだろうが、多くの人たちはそうではないはずだ。

精神医学では「過度に(excessive)」という言葉をよく使う。例えば、うつ病の抑うつ気分という症状は、気分の落ち込みを意味する。誰だって気分が乗らない時や落ち込むとはあるだろう。しかし、ベッドから立てず、ごはんを食べるのもままならず、仕事にもいけず、といった状態が2週間も続けば「質的に」別の問題となる。質的に別の問題、つまり、過度な落ち込みである抑うつ状態、うつ病と呼ばれる状態だ。

ゲームをついつい長くやってしった人は、そのうちゲーム障害になるのか、というと、その可能性は低い。ショックなことで落ち込んでしまった人がかならずうつ病にならないのと同じだ。落ち込んでしまったり、ゲームをついつい長くやり過ぎることは精神疾患ではなく、人間らしい、よくあることである。

日常生活よりゲームを優先する

日常生活という意味が幅広いので、どういう意図で書かれたのかはわからないが、結論からいえば、人の趣味趣向の問題で、人の生き方である。

MMORPGなどではゲームの中でいろいろな人とコミュニケーションをとっていて、リアルな関係よりオンラインの関係の方が魅力的であることもあるだろう。

日常生活に問題が生じる

精神医学では社会的機能と呼んでいるものである。一般的に浸透していない概念だけに、好きなように使われてしまっているのが現状だ。

日常生活に支障とは、通常、GAFが50以下もしくはSOFASが50以下といったものを指す。例えば、うつ病でベットから動けず、会社に行くことができないといった状態であればGAFもSOFASも50以下となる。食物を食べない、排せつの処理が不適切(垂れ流し)などになってくると、GAF、SOFASとも20以下と評価する。

ゲーム障害で、日常生活に支障というのは、学生であれば、少なくとも不登校か、他の社会的活動ができないレベルまで社会的機能が低下しなくてはならない。社会人であれば、仕事に行けない、仕事をすることができない、他の社会的活動もないことが続くことを意味する。

鹿児島県警を含め一般的にはもっとゆるい条件のように解釈されているのではないだろうか。

ゲームに熱中していたので、宿題をするのを忘れてしまったとか、そういうレベルのことを想像しているように思えて仕方ない。

新しいテクノロジーは悪者にされる

ゲームやネットをやり過ぎると、生活に悪影響が出ると考えられがちだが、間違いといってよいだろう。
韓国で行われたインターネット利用した調査ではオンラインゲームをするネット使用時間の長いグループの方が学業成績が良いという結果が出ている。

ides.hatenablog.com

ゲームのせいで成績が落ちたと嘆く親御さんもいるかもしれないが、それはゲームのせいではないかもしれない。学術論文の指摘は、ゲームがなくとも、自分の子どもの成績は落ちていた可能性があったことを気づかせてくれるだろう。ゲームのせいにしたい気持ちはわからなくもないが。

LINEばかりしていて家族の会話が減ったと、親や保護者がボヤくこともあるが、それはLINEのせいなのだろうか。親子関係が冷えていたりするので、子どもがLINEばかりしている可能性はないだろうか。
それに、反抗期に親や保護者と会話を嫌がるのは極めて正常なことでもある。

新しいテクノロジーは悪者扱いされやすいのでテクノロジーは叩かれがちである。しかし、そうやっているうちは、真の原因はみつからず、ただ解決が先延ばしにされているだけだということを忘れない方がいい。

ゲーム障害という精神疾患

ICD-11の採択は2022年の予定なので、将来的にはゲーム障害が精神疾患になるので、間違いとも言い切れない。
今の段階でゲーム障害が精神疾患だ、というのは不正確だろう。

とはいえ、精神障害の構成要件が間違えているので、やはり不適切である。

ゲームを止めましょう

警察という立場で「ゲームを止めろ」と言えば、ゲームの売り上げが下がるなどの効果が予測できる。これは、ゲーム制作会社やゲーム関連の仕事をついている人たちへの威力業務妨害である。

精神医学の知識が足らないのも問題だが、鹿児島県警という立場で発言するべきか否かの判断がついていないことも大きな問題のように感じるのだが、いかがだろうか。

prevalenceパッケージまわりの補足[R]

R 計量

こちらのエントリの補足。

ides.hatenablog.com

感度・特異度がわかっていれば、スクリーニング調査から真の有病率の推定ができる、という手法である。
もちろん、統計的推計なので、診断の有無をしっかりと調べていくようなものを再現できる訳ではなく、いろいろと欠点は存在するが、その点はまた別の機会に。

神戸の新型コロナウィルス抗体検査からの推計

検査人数(神戸): 1000
抗体陽性者数(神戸): 33
感度(クラボウ): 0.7638
特異度(クラボウ): 1

library("prevalence")
Prev <- truePrev(33, 1000, SE = 0.7638, SP = 1, nchains = 4, 
                            burnin = 10000, update = 10000, verbose = FALSE)
summary(Prev)

以前に小数点3位以下が表示できない方法しか示していなかったがsummary()を使えば、普通に使えることがわかった。

$TP
                 mean     median       mode          sd          var       2.5%      97.5% samples
chain 1    0.04433895 0.04404722 0.04423811 0.007389166 5.459978e-05 0.03095667 0.05996641   10000
chain 2    0.04433545 0.04396456 0.04359424 0.007482319 5.598510e-05 0.03073782 0.06040244   10000
chain 3    0.04469223 0.04432558 0.04262175 0.007613095 5.795921e-05 0.03085897 0.06043127   10000
chain 4    0.04431771 0.04396270 0.04456388 0.007538102 5.682298e-05 0.03075770 0.06056437   10000
all chains 0.04442108 0.04407672 0.04416214 0.007507470 5.636211e-05 0.03081621 0.06034275   40000

$SE
             mean median   mode sd var   2.5%  97.5% samples
chain 1    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 2    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 3    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 4    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
all chains 0.7638 0.7638 0.7638  0   0 0.7638 0.7638   40000

$SP
           mean median mode sd var 2.5% 97.5% samples
chain 1       1      1    1  0   0    1     1   10000
chain 2       1      1    1  0   0    1     1   10000
chain 3       1      1    1  0   0    1     1   10000
chain 4       1      1    1  0   0    1     1   10000
all chains    1      1    1  0   0    1     1   40000

　感度・特異度の信頼区間の計算

クラボウのパンフレット(https://www.kurabo.co.jp/bio/English/support/download.php?M=PL&CID=9#catalog111)には感度・特異度が信頼区間の計算をしていないのでやってみよう。

f:id:iDES:20210212150417p:plain

まず分割表を作成する。

dat.kurabo <- as.table(matrix(c(97,0,30,394), nrow = 2, byrow = TRUE))
colnames(dat.kurabo) <- c("PCR+","PCR-")
rownames(dat.kurabo) <- c("Test+","Test-")
dat.kurabo

クラボウの抗体検査キットのテストの分割表である。

      PCR+ PCR-
Test+   97    0
Test-   30  394

感度読緯度はepiRパッケージのepi.testsを使用する。

www.rdocumentation.org

library(epiR)
dat.kurabo <- epi.tests(dat.kurabo, conf.level = 0.95)
print(dat.kurabo); summary(dat.kurabo)

結果は以下。

Point estimates and 95 % CIs:
---------------------------------------------------------
Apparent prevalence                    0.19 (0.15, 0.22)
True prevalence                        0.24 (0.21, 0.28)
Sensitivity                            0.76 (0.68, 0.83)
Specificity                            1.00 (0.99, 1.00)
Positive predictive value              1.00 (0.96, 1.00)
Negative predictive value              0.93 (0.90, 0.95)
Positive likelihood ratio              Inf (NaN, Inf)
Negative likelihood ratio              0.24 (0.17, 0.32)
---------------------------------------------------------

95%信頼区間を考慮したtruePrev()の推定

Prev2 <- truePrev(33, 1000, SE = ~dunif(0.68, 0.83), SP = ~dunif(0.99, 1.00), nchains = 4, 
                            burnin = 10000, update = 10000, verbose = FALSE)
summary(Prev2)

少しまともな推計になっている気がする。

$TP
                 mean     median       mode          sd          var       2.5%      97.5% samples
chain 1    0.03855326 0.03814766 0.03755704 0.008713162 7.591920e-05 0.02252821 0.05702676   10000
chain 2    0.03878964 0.03831414 0.03807106 0.008852792 7.837193e-05 0.02273063 0.05734680   10000
chain 3    0.03872776 0.03827282 0.03691512 0.008820912 7.780849e-05 0.02310364 0.05757549   10000
chain 4    0.03874573 0.03819929 0.03651368 0.008811758 7.764709e-05 0.02292265 0.05766032   10000
all chains 0.03870410 0.03822552 0.03722087 0.008799941 7.743896e-05 0.02279400 0.05738077   40000

$SE
                mean    median      mode         sd         var      2.5%     97.5% samples
chain 1    0.7535388 0.7531481 0.6970268 0.04316914 0.001863575 0.6839307 0.8256182   10000
chain 2    0.7526508 0.7514917 0.7112665 0.04338263 0.001882053 0.6834678 0.8259837   10000
chain 3    0.7521326 0.7505534 0.6963293 0.04311100 0.001858558 0.6833736 0.8257220   10000
chain 4    0.7529574 0.7512809 0.6955498 0.04319945 0.001866192 0.6838645 0.8262247   10000
all chains 0.7528199 0.7516303 0.6949127 0.04321705 0.001867714 0.6836405 0.8258802   40000

$SP
                mean    median      mode          sd          var      2.5%     97.5% samples
chain 1    0.9949075 0.9948642 0.9913735 0.002900040 8.410230e-06 0.9902096 0.9997145   10000
chain 2    0.9950091 0.9949967 0.9986009 0.002892036 8.363873e-06 0.9902785 0.9997547   10000
chain 3    0.9949031 0.9948460 0.9913536 0.002881859 8.305113e-06 0.9902280 0.9997297   10000
chain 4    0.9949990 0.9949875 0.9911999 0.002906795 8.449454e-06 0.9902602 0.9997894   10000
all chains 0.9949547 0.9949262 0.9914553 0.002895512 8.383993e-06 0.9902446 0.9997498   40000

別の書式で

SEとSPを別々に指定し、分布をdistで指定することもできるようだ。

SE <- list(dist = "uniform", min = 0.68, max = 0.83)
SP <- list(dist = "uniform", min = 0.99, max = 1.00)
TP <- truePrev(x = 33, n = 1000, SE = SE, SP = SP, nchains = 4, 
               burnin = 10000, update = 10000, verbose = FALSE)
summary(TP)

プロットも下記のように記述して出力できる。

par(mfrow = c(2, 2))
plot(TP)

f:id:iDES:20210212150441p:plain

密度プロットとトレースプロットだけを表示させるには以下のコードでOKである。

densplot(Prev, exclude_fixed = TRUE)
traceplot(Prev, exclude_fixed = TRUE)

LaTeX形式の数式コードを書くパッケージequatiomatic[R]

R 回帰分析 計量

中澤港さんの日記に書かれていたパッケージ(https://minato.sip21c.org/im3r/20210131.html)。

これを使うと，lmer()で使ったモデルが，そのままTeXのコードになる。LaTeXで論文を書いている人はそのまま取り込めるし，WordやLibreOfficeで論文を書いている人でも，(1)ヘッダをつけてpandocでWordのdocxやWriterのodtに変換するとか，(2)MS Office＋IguanaTexを使ってWordに取り込むとか，(3)LibreOffice＋TexMaths拡張を使ってLibreOffice WriterやImpressに取り込むなどすれば，LaTeXがわからなくても使えてしまう。forecastパッケージのArima()関数にも使えるようだ。

www.rdocumentation.org

主に恩恵を受けるのは、lme4::lmer Modelsとforecast::Arima Modelsを使う人のようだ。

適用できる分析

Model Packages/Functions
linear regression stats::lm
logistic regression stats::glm(family = binomial(link = 'logit'))
probit regression stats::glm(family = binomial(link = 'probit'))
ordinal logistic regression MASS::polr(method = 'logistic'); ordinal::clm(link = 'logit')
ordinal probit regression MASS::polr(method = 'probit'); ordinal::clm(link = 'probit')
auto-regressive integrated moving average forecast::Arima
regression with auto-regressive integrated moving average errors forecast::Arima

回帰分析

mod1 <- lm(mpg ~ cyl + disp, mtcars)
summary(mod1)

結果。

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.66099    2.54700  13.609 4.02e-14 ***
cyl         -1.58728    0.71184  -2.230   0.0337 *  
disp        -0.02058    0.01026  -2.007   0.0542 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.055 on 29 degrees of freedom
Multiple R-squared:  0.7596,    Adjusted R-squared:  0.743 
F-statistic: 45.81 on 2 and 29 DF,  p-value: 1.058e-09

重回帰分析

数式の表示にはextract_eqコマンドを使用する。

library(equatiomatic)
extract_eq(mod1)

$$
\operatorname{mpg} = \alpha + \beta{1}(\operatorname{cyl}) + \beta{2}(\operatorname{disp}) + \epsilon
$$

f:id:iDES:20210211005124p:plain

推定値を入れた式を描く

extract_eq(mod1, use_coefs = TRUE)

$$
\operatorname{\widehat{mpg}} = 34.66 - 1.59(\operatorname{cyl}) - 0.02(\operatorname{disp})
$$

f:id:iDES:20210211005136p:plain

Texコードで一部をアレンジ

extract_eq(mod1, 
           wrap = TRUE, 
           intercept = "\\hat{\\phi}", 
           greek = "\\hat{\\gamma}",
           raw_tex = TRUE)

切片をphi hatにして、βではなくγ hatにした場合。

$$
\begin{aligned}
\operatorname{mpg} &= \hat{\phi} + \hat{\gamma}{1}(\operatorname{cyl}) + \hat{\gamma}{2}(\operatorname{disp}) + \epsilon
\end{aligned}
$$

f:id:iDES:20210211005151p:plain

ロジスティック回帰

model_logit <- glm(sex ~ bill_length_mm + species, 
                   data = penguins, family = binomial(link = "logit"))
extract_eq(model_logit, wrap = TRUE, terms_per_line = 3)

wrap = TRUEで改行した形で表記、terms_per_lineは一行入る数式の数。

$$
\begin{aligned}
\log\left[ \frac { P( \operatorname{sex} = \operatorname{male} ) }{ 1 - P( \operatorname{sex} = \operatorname{male} ) } \right] &= \alpha + \beta{1}(\operatorname{bill_length_mm}) + \beta{2}(\operatorname{species}{\operatorname{Chinstrap}})\ + \ &\quad \beta{3}(\operatorname{species}_{\operatorname{Gentoo}})
\end{aligned}
$$

f:id:iDES:20210211005247p:plain

プロビット回帰

model_probit <- glm(sex ~ bill_length_mm + species, 
                    data = penguins, family = binomial(link = "probit"))
extract_eq(model_probit, wrap = TRUE, terms_per_line = 3)

$$
\begin{aligned}
P( \operatorname{sex} = \operatorname{male} ) &= \Phi[\alpha + \beta{1}(\operatorname{bill_length_mm}) + \beta{2}(\operatorname{species}{\operatorname{Chinstrap}})\ + \
&\qquad\ \beta{3}(\operatorname{species}_{\operatorname{Gentoo}})]
\end{aligned}
$$

f:id:iDES:20210211005257p:plain

順序回帰分析

library(MASS)
set.seed(1234)
df <- data.frame(outcome = factor(rep(LETTERS[1:3], 100),
                                  levels = LETTERS[1:3],
                                  ordered = TRUE),
                 continuous_1 = rnorm(300, 100, 1),
                 continuous_2 = rnorm(300, 50, 5))

model_ologit <- MASS::polr(outcome ~ continuous_1 + continuous_2, 
                           data = df, Hess = TRUE, method = "logistic")
model_oprobit <- MASS::polr(outcome ~ continuous_1 + continuous_2, 
                            data = df, Hess = TRUE, method = "probit")

extract_eq(model_ologit, wrap = TRUE)

$$
\begin{aligned}
\log\left[ \frac { P( \operatorname{A} \geq \operatorname{B} ) }{ 1 - P( \operatorname{A} \geq \operatorname{B} ) } \right] &= \alpha{1} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2}) \
\log\left[ \frac { P( \operatorname{B} \geq \operatorname{C} ) }{ 1 - P( \operatorname{B} \geq \operatorname{C} ) } \right] &= \alpha{2} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2})
\end{aligned}
$$

f:id:iDES:20210211005527p:plain

順序ロジスティック回帰

$$
\begin{aligned}
\log\left[ \frac { P( \operatorname{A} \geq \operatorname{B} ) }{ 1 - P( \operatorname{A} \geq \operatorname{B} ) } \right] &= \alpha{1} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2}) \
\log\left[ \frac { P( \operatorname{B} \geq \operatorname{C} ) }{ 1 - P( \operatorname{B} \geq \operatorname{C} ) } \right] &= \alpha{2} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2})
\end{aligned}
$$

f:id:iDES:20210211005308p:plain

順序プロピット回帰

extract_eq(model_oprobit, wrap = TRUE)

$$ \begin{aligned}
P( \operatorname{A} \geq \operatorname{B} ) &= \Phi[\alpha{1} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2})] \
P( \operatorname{B} \geq \operatorname{C} ) &= \Phi[\alpha{2} + \beta{1}(\operatorname{continuous_1}) + \beta{2}(\operatorname{continuous_2})]
\end{aligned}
$$

f:id:iDES:20210211005413p:plain

lavaanを用いた測定不変性の確認[R]

R lavaan 構造方程式モデリング

lavaanで測定普遍性を確認する方法。

rstudio-pubs-static.s3.amazonaws.com

ある構成要素についてグループ間で比較を行うとき、暗黙のうちに測定の不変性を仮定している。回帰分析、t検定、混合効果モデルなどを行っている時も、その構成要素が同じように機能するであろうという仮定の下での分析になっている。この仮定を確かなものにするための一つの方法が、測定不変性を確認することである。ただし、間接的な方法なので、根本的なアプローチかというとそうでもないところには注意が必要であるが、測定普遍性を確認し、グループ間の違いについて論じる際には、言及が必要となる。

データ

古典的なHolzinger and Swineford (1939)のデータセット。2つの異なる学校（PasteurとGrant-White）の7年生と8年生の子供たちの知能テストのスコアである。オリジナルのデータセット（MBESSパッケージで利用可能）では、26のテストのスコアがある。lavaanで利用できるのは、9つの変数を持つより小さなサブセットであり、このデータセットがより広く使われている（たとえば、Joreskogの1969年の論文では、Grant-Whiteの学校からの145人の被験者のみを使用している）。

library("lavaan")
data("HolzingerSwineford1939")
head(HolzingerSwineford1939)
  id sex ageyr agemo  school grade       x1   x2    x3       x4   x5        x6       x7   x8       x9
1  1   1    13     1 Pasteur     7 3.333333 7.75 0.375 2.333333 5.75 1.2857143 3.391304 5.75 6.361111
2  2   2    13     7 Pasteur     7 5.333333 5.25 2.125 1.666667 3.00 1.2857143 3.782609 6.25 7.916667
3  3   2    13     1 Pasteur     7 4.500000 5.25 1.875 1.000000 1.75 0.4285714 3.260870 3.90 4.416667
4  4   1    13     2 Pasteur     7 5.333333 7.75 3.000 2.666667 4.50 2.4285714 3.000000 5.30 4.861111
5  5   2    12     2 Pasteur     7 4.833333 4.75 0.875 2.666667 4.00 2.5714286 3.695652 6.30 5.916667
6  6   2    14     1 Pasteur     7 5.333333 5.00 2.250 1.000000 3.00 0.8571429 4.347826 6.65 7.500000

x1, x2, x3 は視覚的要因、x4, x5, x6 はテキスト要因、x7, x8, x9はスピード要因である。

確証的因子分析(CFA)

最もシンプルなCFAは以下のものである。

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '
fit <- cfa(HS.model, data=HolzingerSwineford1939)
summary(fit, fit.measures=TRUE)
fitMeasures(fit, c("chisq","df","pvalue","gfi","agfi","cfi","tli","rmsea", "srmr"))

フィッシング指標のみ示しておこう。

 chisq     df pvalue    gfi   agfi    cfi    tli  rmsea   srmr 
85.306 24.000  0.000  0.943  0.894  0.931  0.896  0.092  0.065

悪くはないがいまいちなところもある。

構成モデル Configural model

最初のステップ。モデルが、グループのそれぞれに適合することを示すことである。ここでのグループは学校ごとである。group = "school"を追加する。

configural <- cfa(HS.model, data=HolzingerSwineford1939, group = "school")
summary(configural, fit.measures=TRUE)
fitMeasures(configural, c("chisq","df","pvalue","gfi","agfi","cfi","tli","rmsea", "srmr"))

こちらもフィッシング指標のみ示しておこう。

  chisq      df  pvalue     gfi    agfi     cfi     tli   rmsea    srmr 
115.851  48.000   0.000   0.995   0.989   0.923   0.885   0.097   0.068

やはりいまいちである。

弱不変性 Weak Invariance

測定基準不変性 metric invarianceのことである。

因子負荷がグループ間で等しくなるように制約するモデルである。これは因子がグループ間で同じ有意性を持っていることを仮定している。

weak.invariance <- cfa(HS.model, data=HolzingerSwineford1939, group = "school", group.equal = "loadings")
summary(weak.invariance, fit.measures = TRUE)
fitMeasures(weak.invariance, c("chisq","df","pvalue","gfi","agfi","cfi","tli","rmsea", "srmr"))

こちらもフィッシング指標のみ示しておこう。

  chisq      df  pvalue     gfi    agfi     cfi     tli   rmsea    srmr 
124.044  54.000   0.000   0.995   0.989   0.921   0.895   0.093   0.072

尤度比検定　構成モデルvs.弱不変性

lavTestLRT(weak.invariance, configural)

少し前まではanova(a,b)で動作していたようだが、現在はlavTestLRTに変更されている。

Chi-Squared Difference Test

                Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)
configural      48 7484.4 7706.8 115.85                              
weak.invariance 54 7480.6 7680.8 124.04     8.1922       6     0.2244

帰無仮説は、構成モデルは弱い不変性モデルよりも適合性が向上していないというものだ。P=0.2244であったため有意だとは認められず、帰無仮説は棄却されない。これは、弱い不変性があることの証拠である。測定普遍性を示す場合、P<.05になるとマズいということだ。

フィッシング指標をまとめるコードもあったので掲載しておこう。

fit.stats <- rbind(fitmeasures(configural, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
fitmeasures(weak.invariance, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")))
rownames(fit.stats) <- c("configural", "weak invariance")
fit.stats

内容は必要なものを入れて使うとよいだろう。

                   chisq df      rmsea       tli       cfi      aic
configural      115.8513 48 0.09691486 0.8850976 0.9233984 7484.395
weak invariance 124.0435 54 0.09283654 0.8945646 0.9209235 7480.587

強不変性 strong invariance

尺度妥当性scalar invarianceのことである。

強不変性は切片がグループ間で等しいという制約を追加する。これは、構成要素（因子負荷）の有意味性と、基礎となる顕在変数（切片）のレベルが、両方のグループで等しいというモデルである。具体的にはgroup.equal = c( "loadings", "intercepts")を追加する。loadingsの代わりにinterceptsを追加するのでは、両方とも入れる。

strong.invariance <- cfa(HS.model, data=HolzingerSwineford1939, group = "school", 
                         group.equal = c( "loadings", "intercepts"))
summary(strong.invariance, fit.measures = TRUE)
fitMeasures(strong.invariance, c("chisq","df","pvalue","gfi","agfi","cfi","tli","rmsea", "srmr"))

こちらもフィッシング指標のみ示しておこう。

  chisq      df  pvalue     gfi    agfi     cfi     tli   rmsea    srmr 
164.103  60.000   0.000   0.993   0.987   0.882   0.859   0.107   0.082

視覚、テキスト、スピードの潜在平均が推定されている。強不変性があるかどうかをテストするまでは、グループ間で比較すべきではないとされている。

尤度比検定　弱不変性vs. 強不変性

lavaanはベータ版と宣言しているようにコマンドが時々変わり、尤度比検定のコマンドが変更されている(参照)

lavTestLRT(strong.invariance, weak.invariance)
Chi-Squared Difference Test

                  Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
weak.invariance   54 7480.6 7680.8 124.04                                  
strong.invariance 60 7508.6 7686.6 164.10     40.059       6  4.435e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

P<.05であるため、帰無仮説が棄却される。
弱不変モデルが強不変モデルよりもよく適合しているようだ。弱不変モデルの方がAICとBICが低いことからも明らかである。

次にすることは部分不変性partial invarianceがあるかどうかを調べることである。lavaanのlavTestScore()を使用することで確認できる。この関数を使うことで、グループ間の等値制約を解放するべき場所を見つけることができる。modindices()は、新しいパスに関連して新たに追加されたパラメータの修正インデックスのみを表示する。

lavTestScore(strong.invariance)
$test

total score test:

   test     X2 df p.value
1 score 46.956 15       0

$uni

univariate score tests:

     lhs op   rhs     X2 df p.value
1   .p2. == .p38.  0.306  1   0.580
2   .p3. == .p39.  1.636  1   0.201
3   .p5. == .p41.  2.744  1   0.098
4   .p6. == .p42.  2.627  1   0.105
5   .p8. == .p44.  0.027  1   0.871
6   .p9. == .p45.  0.004  1   0.952
7  .p25. == .p61.  5.847  1   0.016
8  .p26. == .p62.  6.863  1   0.009
9  .p27. == .p63. 19.193  1   0.000
10 .p28. == .p64.  2.139  1   0.144
11 .p29. == .p65.  1.563  1   0.211
12 .p30. == .p66.  0.032  1   0.857
13 .p31. == .p67. 15.021  1   0.000
14 .p32. == .p68.  4.710  1   0.030
15 .p33. == .p69.  1.498  1   0.221

pはparTable()を使うと具体的に確認できる。

options(max.print=100000) 
parTable(strong.invariance)
   id     lhs op     rhs user block group free ustart exo label plabel start    est    se
1   1  visual =~      x1    1     1     1    0      1   0         .p1. 1.000  1.000 0.000
2   2  visual =~      x2    1     1     1    1     NA   0  .p2.   .p2. 0.769  0.576 0.101
3   3  visual =~      x3    1     1     1    2     NA   0  .p3.   .p3. 1.186  0.798 0.112
4   4 textual =~      x4    1     1     1    0      1   0         .p4. 1.000  1.000 0.000
5   5 textual =~      x5    1     1     1    3     NA   0  .p5.   .p5. 1.237  1.120 0.066
6   6 textual =~      x6    1     1     1    4     NA   0  .p6.   .p6. 0.865  0.932 0.056
7   7   speed =~      x7    1     1     1    0      1   0         .p7. 1.000  1.000 0.000
8   8   speed =~      x8    1     1     1    5     NA   0  .p8.   .p8. 1.227  1.130 0.145
9   9   speed =~      x9    1     1     1    6     NA   0  .p9.   .p9. 0.827  1.009 0.132
10 10      x1 ~~      x1    0     1     1    7     NA   0        .p10. 0.698  0.555 0.139
11 11      x2 ~~      x2    0     1     1    8     NA   0        .p11. 0.752  1.296 0.158
12 12      x3 ~~      x3    0     1     1    9     NA   0        .p12. 0.673  0.944 0.136
13 13      x4 ~~      x4    0     1     1   10     NA   0        .p13. 0.660  0.445 0.069
14 14      x5 ~~      x5    0     1     1   11     NA   0        .p14. 0.854  0.502 0.082
15 15      x6 ~~      x6    0     1     1   12     NA   0        .p15. 0.487  0.263 0.050
16 16      x7 ~~      x7    0     1     1   13     NA   0        .p16. 0.585  0.888 0.120
17 17      x8 ~~      x8    0     1     1   14     NA   0        .p17. 0.476  0.541 0.095
18 18      x9 ~~      x9    0     1     1   15     NA   0        .p18. 0.489  0.654 0.096
19 19  visual ~~  visual    0     1     1   16     NA   0        .p19. 0.050  0.796 0.172
20 20 textual ~~ textual    0     1     1   17     NA   0        .p20. 0.050  0.879 0.131
21 21   speed ~~   speed    0     1     1   18     NA   0        .p21. 0.050  0.322 0.082
22 22  visual ~~ textual    0     1     1   19     NA   0        .p22. 0.000  0.410 0.095
23 23  visual ~~   speed    0     1     1   20     NA   0        .p23. 0.000  0.178 0.066
24 24 textual ~~   speed    0     1     1   21     NA   0        .p24. 0.000  0.180 0.062
25 25      x1 ~1            0     1     1   22     NA   0 .p25.  .p25. 4.941  5.001 0.090
26 26      x2 ~1            0     1     1   23     NA   0 .p26.  .p26. 5.984  6.151 0.077
27 27      x3 ~1            0     1     1   24     NA   0 .p27.  .p27. 2.487  2.271 0.083
28 28      x4 ~1            0     1     1   25     NA   0 .p28.  .p28. 2.823  2.778 0.087
29 29      x5 ~1            0     1     1   26     NA   0 .p29.  .p29. 3.995  4.035 0.096
30 30      x6 ~1            0     1     1   27     NA   0 .p30.  .p30. 1.922  1.926 0.079
31 31      x7 ~1            0     1     1   28     NA   0 .p31.  .p31. 4.432  4.242 0.073
32 32      x8 ~1            0     1     1   29     NA   0 .p32.  .p32. 5.563  5.630 0.072
33 33      x9 ~1            0     1     1   30     NA   0 .p33.  .p33. 5.418  5.465 0.069
34 34  visual ~1            0     1     1    0      0   0        .p34. 0.000  0.000 0.000
35 35 textual ~1            0     1     1    0      0   0        .p35. 0.000  0.000 0.000
36 36   speed ~1            0     1     1    0      0   0        .p36. 0.000  0.000 0.000
37 37  visual =~      x1    1     2     2    0      1   0        .p37. 1.000  1.000 0.000
38 38  visual =~      x2    1     2     2   31     NA   0  .p2.  .p38. 0.896  0.576 0.101
39 39  visual =~      x3    1     2     2   32     NA   0  .p3.  .p39. 1.155  0.798 0.112
40 40 textual =~      x4    1     2     2    0      1   0        .p40. 1.000  1.000 0.000
41 41 textual =~      x5    1     2     2   33     NA   0  .p5.  .p41. 0.991  1.120 0.066
42 42 textual =~      x6    1     2     2   34     NA   0  .p6.  .p42. 0.962  0.932 0.056
43 43   speed =~      x7    1     2     2    0      1   0        .p43. 1.000  1.000 0.000
44 44   speed =~      x8    1     2     2   35     NA   0  .p8.  .p44. 1.282  1.130 0.145
45 45   speed =~      x9    1     2     2   36     NA   0  .p9.  .p45. 0.895  1.009 0.132
46 46      x1 ~~      x1    0     2     2   37     NA   0        .p46. 0.659  0.654 0.128
47 47      x2 ~~      x2    0     2     2   38     NA   0        .p47. 0.613  0.964 0.123
48 48      x3 ~~      x3    0     2     2   39     NA   0        .p48. 0.537  0.641 0.101
49 49      x4 ~~      x4    0     2     2   40     NA   0        .p49. 0.629  0.343 0.062
50 50      x5 ~~      x5    0     2     2   41     NA   0        .p50. 0.671  0.376 0.073
51 51      x6 ~~      x6    0     2     2   42     NA   0        .p51. 0.640  0.437 0.067
52 52      x7 ~~      x7    0     2     2   43     NA   0        .p52. 0.531  0.625 0.095
53 53      x8 ~~      x8    0     2     2   44     NA   0        .p53. 0.547  0.434 0.088
54 54      x9 ~~      x9    0     2     2   45     NA   0        .p54. 0.526  0.522 0.086
55 55  visual ~~  visual    0     2     2   46     NA   0        .p55. 0.050  0.708 0.160
56 56 textual ~~ textual    0     2     2   47     NA   0        .p56. 0.050  0.870 0.131
57 57   speed ~~   speed    0     2     2   48     NA   0        .p57. 0.050  0.505 0.115
58 58  visual ~~ textual    0     2     2   49     NA   0        .p58. 0.000  0.427 0.097
59 59  visual ~~   speed    0     2     2   50     NA   0        .p59. 0.000  0.329 0.082
60 60 textual ~~   speed    0     2     2   51     NA   0        .p60. 0.000  0.236 0.073
61 61      x1 ~1            0     2     2   52     NA   0 .p25.  .p61. 4.930  5.001 0.090
62 62      x2 ~1            0     2     2   53     NA   0 .p26.  .p62. 6.200  6.151 0.077
63 63      x3 ~1            0     2     2   54     NA   0 .p27.  .p63. 1.996  2.271 0.083
64 64      x4 ~1            0     2     2   55     NA   0 .p28.  .p64. 3.317  2.778 0.087
65 65      x5 ~1            0     2     2   56     NA   0 .p29.  .p65. 4.712  4.035 0.096
66 66      x6 ~1            0     2     2   57     NA   0 .p30.  .p66. 2.469  1.926 0.079
67 67      x7 ~1            0     2     2   58     NA   0 .p31.  .p67. 3.921  4.242 0.073
68 68      x8 ~1            0     2     2   59     NA   0 .p32.  .p68. 5.488  5.630 0.072
69 69      x9 ~1            0     2     2   60     NA   0 .p33.  .p69. 5.327  5.465 0.069
70 70  visual ~1            0     2     2   61     NA   0        .p70. 0.000 -0.148 0.122
71 71 textual ~1            0     2     2   62     NA   0        .p71. 0.000  0.576 0.117
72 72   speed ~1            0     2     2   63     NA   0        .p72. 0.000 -0.177 0.090
73 73    .p2. ==   .p38.    2     0     0    0     NA   0              0.000  0.000 0.000
74 74    .p3. ==   .p39.    2     0     0    0     NA   0              0.000  0.000 0.000
75 75    .p5. ==   .p41.    2     0     0    0     NA   0              0.000  0.000 0.000
76 76    .p6. ==   .p42.    2     0     0    0     NA   0              0.000  0.000 0.000
77 77    .p8. ==   .p44.    2     0     0    0     NA   0              0.000  0.000 0.000
78 78    .p9. ==   .p45.    2     0     0    0     NA   0              0.000  0.000 0.000
79 79   .p25. ==   .p61.    2     0     0    0     NA   0              0.000  0.000 0.000
80 80   .p26. ==   .p62.    2     0     0    0     NA   0              0.000  0.000 0.000
81 81   .p27. ==   .p63.    2     0     0    0     NA   0              0.000  0.000 0.000
82 82   .p28. ==   .p64.    2     0     0    0     NA   0              0.000  0.000 0.000
83 83   .p29. ==   .p65.    2     0     0    0     NA   0              0.000  0.000 0.000
84 84   .p30. ==   .p66.    2     0     0    0     NA   0              0.000  0.000 0.000
85 85   .p31. ==   .p67.    2     0     0    0     NA   0              0.000  0.000 0.000
86 86   .p32. ==   .p68.    2     0     0    0     NA   0              0.000  0.000 0.000
87 87   .p33. ==   .p69.    2     0     0    0     NA   0              0.000  0.000 0.000

total score test:として出力されているのはラグランジュ乗数検定である。これは、等値制約を解放することによって、ベース・モデルに対する適合性の改善を表すかどうかの検定である。カイ2乗差検定で帰無仮説が棄却されるものを探す。カイ2乗差（X2）で最大の変化を持つのはx3（.p27 == .p63）である（カイ二乗値19.193であり、P=0.000である)。ここの切片を解放すべきであることがわかる。最初にこの点をモデルへ反映させよう。

部分的強不変性を導くステップ1

先の切片の等値制約を外すのはx3~x1であるため、group.partial = c("x3 ~ 1")を指定する。

strong.invariance.x3 <- cfa(HS.model, data=HolzingerSwineford1939, group = "school", 
                            group.equal = c( "loadings", "intercepts"), 
                            group.partial = c("x3 ~ 1"))
lavTestScore(strong.invariance.x3)

結果。

$test

total score test:

   test     X2 df p.value
1 score 27.528 14   0.016

$uni

univariate score tests:

     lhs op   rhs     X2 df p.value
1   .p2. == .p38.  0.734  1   0.392
2   .p3. == .p39.  0.485  1   0.486
3   .p5. == .p41.  2.760  1   0.097
4   .p6. == .p42.  2.630  1   0.105
5   .p8. == .p44.  0.026  1   0.872
6   .p9. == .p45.  0.002  1   0.960
7  .p25. == .p61.  2.833  1   0.092
8  .p26. == .p62.  2.833  1   0.092
9  .p28. == .p64.  2.136  1   0.144
10 .p29. == .p65.  1.560  1   0.212
11 .p30. == .p66.  0.032  1   0.857
12 .p31. == .p67. 15.023  1   0.000
13 .p32. == .p68.  4.727  1   0.030
14 .p33. == .p69.  1.492  1   0.222

ラグランジュ検定が棄却されているので、さらに等値制約の解放が必要であることがわかる。
次に最大のカイ2乗差（X2）を持つのはx7 (.p31. == .p67.)であるため、モデルをさらに改変する。

ステップ2

strong.invariance.x3x7 <- cfa(HS.model, data=HolzingerSwineford1939, group = "school", 
                              group.equal = c( "loadings", "intercepts"),
                              group.partial = c("x3 ~ 1", "x7 ~ 1"))
lavTestScore(strong.invariance.x3x7)
$test

total score test:

   test     X2 df p.value
1 score 12.583 13   0.481

$uni

univariate score tests:

     lhs op   rhs    X2 df p.value
1   .p2. == .p38. 0.734  1   0.391
2   .p3. == .p39. 0.492  1   0.483
3   .p5. == .p41. 2.769  1   0.096
4   .p6. == .p42. 2.631  1   0.105
5   .p8. == .p44. 0.013  1   0.910
6   .p9. == .p45. 0.062  1   0.803
7  .p25. == .p61. 2.832  1   0.092
8  .p26. == .p62. 2.832  1   0.092
9  .p28. == .p64. 2.135  1   0.144
10 .p29. == .p65. 1.563  1   0.211
11 .p30. == .p66. 0.032  1   0.858
12 .p32. == .p68. 0.053  1   0.818
13 .p33. == .p69. 0.053  1   0.818

P=0.481となったため、ラグランジュ検定が棄却されていないことを確認。

次に行うのは、弱不変性モデルとの比較である。切片に等値制約を掛けた場合、P<.05になっていたので、部分的強不変モデルを作成したわけだが、この作成したモデルが弱不変性モデルと比べて改善できているかを確認する。

lavTestLRT(strong.invariance.x3x7, weak.invariance)

結果。

Chi-Squared Difference Test

                       Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)
weak.invariance        54 7480.6 7680.8 124.04                              
strong.invariance.x3x7 58 7478.0 7663.3 129.42     5.3789       4     0.2506

P<.05ではないため、測定不変性が確認できている。ちなみに、このコマンドはanova(a,b)がまだ通るようである(不思議)。

fit.stats2 <- rbind(fitmeasures(strong.invariance, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")),
fitmeasures(strong.invariance.x3x7, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")))
rownames(fit.stats2) <- c("strong", "strong with x3 x7")
fit.stats <- rbind(fit.stats, fit.stats2)
round(fit.stats, 4)

フィッティング指標をモデルごとに比較してみよう。

                     chisq df  rmsea    tli    cfi      aic
configural        115.8513 48 0.0969 0.8851 0.9234 7484.395
weak invariance   124.0435 54 0.0928 0.8946 0.9209 7480.587
strong            164.1028 60 0.1074 0.8590 0.8825 7508.646
strong with x3 x7 129.4225 58 0.0905 0.8999 0.9194 7477.966
summary(strong.invariance.x3x7)

切片の一部だけ表示。

Group 1 [Pasteur]:
Intercepts:
                   Estimate  Std.Err  Z-value  P(>|z|)
    visual            0.000                           
    textual           0.000                           
    speed             0.000                        

Group 2 [Grant-White]:
Intercepts:
                   Estimate  Std.Err  Z-value  P(>|z|)
    visual            0.051    0.129    0.393    0.695
    textual           0.576    0.117    4.918    0.000
    speed            -0.071    0.089   -0.800    0.424

テキストについては強不変性があり、視覚と速度の両方の要因については部分的に強不変性があることがわかった。学校2(Grant-White)の生徒が、学校1(Pasteur)よりもテキストの潜在的構成要素で優れていることが判明した。彼らはPasteurの生徒よりも0.576高くなると予想されます。部分的な強不変性を許容すると、Grant-Whiteの生徒とPasteurの生徒の間には、視覚と速度に関して差がない。

厳密不変性 Strict Invariance

厳密不変性はグループ間の顕在変数の残差分散が等しいという制約を追加する。今までの分析で部分的不変性しか確認できなかったため、あるべきパラメータ（たとえば、x3 と x7 の切片）を自由にして、残差分散と同様に、他のすべてのパラメータ（すなわち、因子負荷量と切片）を制限することになる。raw score に基づく比較を行いたい場合（例えば，変数のスコアの合計に基づく）我々は厳密不変性を必要になる。 これは，観察された分散が真のスコア分散と残差/誤差分散の合計であるからである。残差分散が同じであれば、それらは同量の真のスコア分散を持ちます。前と同様に潜在平均を推定する。

等値制約に"residuals"を追加する。

strict.invariance.x3x7 <- cfa(HS.model, data=HolzingerSwineford1939, group = "school", 
                              group.equal = c( "loadings", "intercepts", "residuals"), 
                              group.partial = c("x3 ~ 1", "x7 ~ 1"))
anova(strong.invariance.x3x7, strict.invariance.x3x7)

こちらはlavTestLRT(a,b)は通らず、anova(a,b)しか通らないので注意。そのうち通るようになるのだろうが。

Chi-Squared Difference Test

                       Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)  
strong.invariance.x3x7 58 7478.0 7663.3 129.42                                
strict.invariance.x3x7 67 7477.8 7629.8 147.26     17.838       9     0.0371 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

カイ二乗差検定は有意だが、P値が.037というのは、それほど違和感があるわけではない。他の適合尺度を見よう。

fit.stats <- rbind(fit.stats, fitmeasures(strict.invariance.x3x7, fit.measures = c("chisq", "df", "rmsea", "tli", "cfi", "aic")))
rownames(fit.stats)[5] <- "strict invariance"
round(fit.stats, 4)
                     chisq df  rmsea    tli    cfi      aic
configural        115.8513 48 0.0969 0.8851 0.9234 7484.395
weak invariance   124.0435 54 0.0928 0.8946 0.9209 7480.587
strong            164.1028 60 0.1074 0.8590 0.8825 7508.646
strong with x3 x7 129.4225 58 0.0905 0.8999 0.9194 7477.966
strict invariance 147.2605 67 0.0892 0.9026 0.9094 7477.804

厳密不変性は、x3 と x7 の切片が自由である強い不変性よりも、限定的な改善であることがわかる。完全な厳密不変性があれば、顕在変数がグループ間で等しく信頼できると言えるだろう。これは、平均の比較やグループ間の回帰を実行するために、合計スコアを使用できることを意味する。完全な厳密不変性を持っていないので、この性質で比較できる唯一の尺度は、テキスト要因である。構造モデルについて述べるときは、この問題に戻り、構造的不変性、具体的には因子分散、因子共分散、SEMにおける回帰係数について議論する。

実用例

こちらの文献では次のような表記がされている。

Adaptation of the Bangla Version of the COVID-19 Anxiety Scale

f:id:iDES:20210207160130p:plain

関谷徳泰さんのブログ記事(測定不変性(measurement invariance)、測定等価性(measurement equivalence) | Dr.Clover's Computer Clinic)では次のように例が挙げられていた。

f:id:iDES:20210207160403p:plain

ERIC - EJ1038101 - Multiple-Group Confirmatory Factor Analysis in R--A Tutorial in Measurement Invariance with Continuous and Ordinal Indicators, Practical Assessment, Research & Evaluation, 2014-Jul

統計量をどこまで表記するかは場合に拠りけりのようだが、CFIを優先すべきであるようだ。

lavaanのモデルの尤度比検定[R]

R lavaan 構造方程式モデリング

lavTestLRT関数を用いる。昔のバージョンではanova(fit1, fit0)という書き方をしたようだが、現在はlavTestLRT関数に変更されている。

www.rdocumentation.org

LRTはLikelihood-ratio testつまり尤度比検定である。ブートストラップをかけて行うとBLRTである。lavaanで実施できるかはよく知らない。

例題

library(lavaan)
HS.model <- '
    visual  =~ x1 + b1*x2 + x3
    textual =~ x4 + b2*x5 + x6
    speed   =~ x7 + b3*x8 + x9
'
fit1 <- cfa(HS.model, data = HolzingerSwineford1939)
fit0 <- cfa(HS.model, data = HolzingerSwineford1939, 
            orthogonal = TRUE)
lavTestLRT(fit1, fit0)

結果。

Chi-Squared Difference Test

     Df    AIC    BIC   Chisq Chisq diff Df diff Pr(>Chisq)    
fit1 24 7517.5 7595.3  85.305                                  
fit0 27 7579.7 7646.4 153.527     68.222       3  1.026e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

orthogonal = TRUEというオプションはCFAモデル内の潜在変数のすべての共分散が直交するように制約するときのオプションである。
fit1とfit0を比較したときに、制約をかけたfit0の法が有意になっているため、モデルが改善されたことがわかる。

type = "Chisq"で、検定統計量がスケーリングされている場合、特別なスケーリングされた差の検定統計量が計算される。。手法が "satorra.bentler.2001" の場合，Satorra & Bentler (2001) で記述されている単純な近似が使用される。いくつかの設定では，これは負の検定統計量につながることがある。正の検定統計量を確保するために、Satorra & Bentler (2010)によって提案された方法を使用することができる。あるいは、方法が"satorra.2000 "の場合は、Satorra (2000)のオリジナルの公式を使用する。

オプションはlavaanのドキュメントを参照のこと。