井出草平の研究ノート

RStanのインストールに関するメモ(2022年7月)[R]

2022年7月現在の、RにStanをインストールする手順を記載する。
現在の注意点はRStanはR4.2に対応していない、ということである。最新バージョンのRとRtoolsではRStanは動かないことに注意が必要である。 なお、Windowsでの方法であり、Macでは試していないのでわからないものの、R4.2がRStanに対応していないのは同じだと思われる。

実行可能な環境はRバージョンは4.1.3、Rtoolsバージョンは4.0

R version: 4.1.3
Rtools version : 4.0 https://discourse.mc-stan.org/t/issues-while-installing-rstan-after-r-update-to-4-2-1/27956/2

Windowsをお使いのStanユーザーは、可能であればR4.2へのアップグレードを控えていただくことをお勧めする。
https://blog.mc-stan.org/2022/04/26/stan-r-4-2-on-windows/

R4.1.3(Windows)

ダウンロードページ。

https://cran.r-project.org/bin/windows/base/old/4.1.3/

Rtoolsの設定

https://cran.r-project.org/bin/windows/Rtools/rtools40.html

Windowsの場合は下記をダウンロードする。

On Windows 64-bit: rtools40-x86_64.exe (includes both i386 and x64 compilers).

RtoolsのPATHの追加

インストールが完了したら、Rパッケージをコンパイルできるようにするために、もう1つのステップを実行する必要がある。それは、Rtoolsのmakeユーティリティ(bash、makeなど)の場所をPATHに置くことです。最も簡単な方法は、Documentsフォルダーに以下の行を含むテキストファイル.Renvironを作成することである。

PATH="${RTOOLS40_HOME}\usr\bin;${PATH}"

これはテキストエディタで行うこともできる、Rから以下のように行うこともできる(Rのコードではバックスラッシュをエスケープする必要があることに注意)。

write('PATH="${RTOOLS40_HOME}}\usrr\bin;${PATH}"', file = "~/.Renviron", append = TRUE)

Rを再起動し、Rtoolsのインストール先が表示され、makeが見つかるか確認する。

Sys.which("make")
## "C:\\rtools40\\usr\\bin\\make.exe"

次に、Rパッケージをソースからインストールしてみる。

install.packages("jsonlite", type = "source")

これで成功すればOK。

RStanのインストール

こちらは日本語のガイドがある。
下記のガイドどおり実行すれば、インストールは可能である。

https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started-(Japanese)

尤度に基づいた検定、IC指標の比較

Nylund、Asparouhov、Muthénの論文からnaive chi-square(NCS)、Lo–Mendell–Rubin、bootstrap likelihood ratio testの比較をした論文の結論部分。BLRTの優位性が示された論文として見かけることがあるが、あくまでもデータ次第であることが指摘されてことが抜け落ちていることが多いのではないかと思われる。

実践的に検定を用いていると、必ずしもBLRTが最適ではないケースに遭遇する。LMR、BICAICを示すことが必要になるだろう。

https://www.statmodel.com/download/LCA_tech11_nylund_v83.pdf

  • Nylund, K. L., Asparouhov, T., & Muthén, B. O. (2007). Deciding on the Number of Classes in Latent Class Analysis and Growth Mixture Modeling: A Monte Carlo Simulation Study. Structural Equation Modeling: A Multidisciplinary Journal, 14(4), 535–569. https://doi.org/10.1080/10705510701575396

尤度に基づいた検定

表8は、3つの尤度ベースの検定(NCS、LMR、BLRT)の性能に関する情報を含んでいる。表7に示した結果と同様に、これらの結果は、一連のモデルについて、異なるクラス数での比較に基づいている。 これらの値を計算するために、各複製について、2クラスから6クラスのモデルを指定した結果のp値を調べた。なお、これらのLMRとBLRTでは、2クラスモデルを指定した場合、1クラスモデルと2クラスモデルを比較している。また、2クラスから6クラスの範囲の代替モデルを指定したことに注意。したがって、表8では、LCAモデルでは1クラスモデルから5クラスモデルまで、FMAモデルとGMMモデルでは1クラスから3クラスまでがLMRとBLRTの可能解となる。NCSについては、伝統的なカイ二乗分布を用いてp値を求める。このp値は,指定されたモデルとクラスが1つ少ないモデルとの間に有意な改善があるかどうかを評価するために使用される。これらのp値を見て、最初の有意でないp値(p > .05)の発生に基づいて選択されたモデルを特定した。.95に近い太字の数値は、テストが完璧に近いパフォーマンスを示したことになる。カテゴリ別LCA8クラスモデルのLMR列の太字の数値は、n=500 は,再現実験の78%が4クラスの解が正しいモデルであると結論づけたことを示している。

表8の結果は、NCSやLMRに対するBLRTの明確な優位性を示している。ほぼ全てのLCAモデルとサンプルサイズを考慮した場合、BLRTはほぼ95%の確率で真のkクラスモデルを正しく同定することができる。NCSを見て、LCAモデルを考慮すると、サンプルサイズが大きくなるにつれて性能が低下し、一貫して低いままであることに気づくま。FMAモデルでは、NCSの性能は一貫して低く、具体的には、すべてのサンプルサイズにおいて、正しいモデルを約60%の時間でのみ識別している。しかし、GMMモデルについては、BLRTほどではないが、NCSは一貫して良好な結果を示した。カテゴリー別LCA15項目モデルでは、LMRはBLRTとほぼ同等の性能を示し、90%から94%の確率でLMRは正しいクラス数を選択した。同じ設定において、NCSの性能は非常に低く、正しいモデルを識別できたのはせいぜい3%であった。連続的な結果を持つLCAモデルでは、LMRとBLRTはn = 500と1,000で非常に良い性能を発揮する。カテゴリカルと連続のLCAの両方において、15項目のモデルでは、LMRは、すべてのサンプルサイズにおいて、90%以上の正しい割合でモデルを同定し、むしろ良いパフォーマンスを示す。
このシミュレーションでは、3つの尤度検定(NCS、LMR、BLRT)と、カテゴリと連続の両方の結果を持つ混合モデルにおいてクラス数を決定するために使用される一般的な適合指標の性能について検討した。 我々は、これらのモデルにNCS検定を使用した場合、新たに提案されたLMRやBLRTと比較して、どの程度誤解を招くかを理解することに努めた。また、従来のIC指標についても検討した。

尤度に基づく検定の比較

本研究の結果は、k - 1 対 k クラスモデルの検定に NCS を使用した場合、真のモデルを棄却する頻度が高すぎるという事実を支持するものであった。NCSはサンプルサイズに特に敏感であり、LCAモデルについてはサンプルサイズが大きくなるにつれてその性能は実際に悪化することが指摘されている。前述したように、この検定は入れ子型混合モデルの検定には適していないことが分かっている。この仮定違反は、サンプルサイズが大きくなるにつれて増幅される可能性がある。表8で見られるように、NCSがうまくいかないとき、クラスの数を過大評価する傾向があることに注意することが重要である。NCSの代替として、他のLRTであるLMRとBLRTは、より正確に正しいモデルを識別する能力を有することが示された。BLRTはほぼ全てのモデリング設定において、LMRよりも明らかに優れた性能を示した。
NCSテストのType 1エラー率が高いので、LMRとBLRTの検出力の結果のみを考察する。本研究で検討したほぼすべてのモデルで、LMRとBLRTの両方が良好な検出力を有している。BLRTは、LMRよりも、すべてのサンプルサイズにわたって、より一貫した検出力を持つ。BLRTとLMRの検出力のわずかな差は,検出力(すなわち,kクラスとk - 1クラスのモデルを区別すること)の点で,2つのテストの性能に大きな差があることを示すものではない。
表 8 の LMR の結果を考慮すると、一般的に LCA モデルにおいて LMR がモデルを正しく認識できない場合、クラス数を過大評価する傾向があることがわかる。したがって、LMRをクラス列挙のツールとして使用する分析者は、LMRのp値が有意でない場合、最大でその数のクラスが存在すると確信することができるが、実際にはもっと少ないかもしれない。クラス数を過大評価することは、クラスを過小評価することよりも良いと考えることができる。なぜなら、k + 1クラスの解から真のkクラスの解を抽出することができるからである。例えば、真の解がk = 3で、モデルLMRが4クラス解を同定した場合、4クラス解のクラスの1つが実質的に意味をなさないか、同定が困難な非常に小さなクラスが存在する場合がある。その結果、LMRの結果にもかかわらず4クラス解を採用せず、真の3クラス解に落ち着くということもあり得る。しかし、LMRがあった場合、BICの優位性がより明らかになる。BICは修正BICよりも一貫して正しいモデルを識別し、修正BICはn = 200の8項目モデルで49%まで低下する。FMAとGMMの両方において、BICはうまく機能しており、最悪の場合、84%の確率で正しいモデルを識別することができる。これらの結果から、BICはクラス数を正しく特定するために検討されたICの中で最も整合的なICであると結論づけられる。表7より、BICはモデルの種類によらず、サンプルサイズが小さい場合に感度を持つことがわかる。 これらのICの性能に及ぼす構造の影響をより深く理解するためには、さらなる研究が必要であることは間違いない。

BICとBLRTの比較

表7と表8の結果を比較して、ICの中で最も性能が良いBICと、LRTの中で最も性能が良いBLRTの性能を理解することができる。カテゴリカルな結果を持つLCAでは、BICよりもBLRTの方が正しいクラス数を特定することに一貫性がある。なぜなら、最悪の場合、BLRTは49%の確率で正しいクラス数を特定するからである。これは、不等クラス、n = 200を持つ8項目のカテゴリ結果LCAモデルにおいて、最悪の場合、正しいクラス数を全く識別できないBICよりも優れている。連続的な結果を持つLCAを考慮すると、BICは良いパフォーマンスを示すが、n = 200に対して74%しか正しいモデルを同定することができない。この設定において、BLRTは一貫しており、94%の確率で正しいモデルを同定している。表7と表8の結果をFMAとGMMモデルについて比較すると、BICとBLRTはともに良好な結果を示している。この設定において、BICは最悪の場合、n = 200のGMMに対して84%の割合で正しいモデルを同定し、BLRTは最悪の場合、n = 200のFMAに対して87%の割合で正しいモデルを同定している。このように、本研究のすべてのモデルを考慮し、表7と表8の結果を比較すると、BLRTは正しいクラス数を示す最も一貫した指標として際立っていることがわかる。
ここで示した結果は、すべての混合モデルのごく一部であるが、FMA と GMM の結果は、それぞれ 1 つのモデルのみに基づいていることに注意することが重要である。したがって,表 7 と表 8 に示した結果を比較する場合,FMA と GMM の結果は LCA モデルほどには重視されない.Tofighi and Enders (2006) による最近のシミュレーション研究では,より広範な GMM モデルについて,クラス列挙の問題がより詳細に検討されている.
研究の結果をまとめると、BLRT テストは NCS や LMR と比較して明らかに有利であり、LCA モデルのクラス数を決定するための信頼できるツールとして使用できることが示された。BICは、LCAモデル、FMAおよびGMMモデルにおいて、検討した他のICよりもクラス数を正しく同定できることがわかった。LMRはNCSよりも良い性能を示したが、BLRTほどではなかった。もし、IC指標を一つ選ばなければならないとしたら、BICはクラス数の最も良い指標と思われるツールであろう。BLRTは、正しいクラスモデルを選択する一貫性から、LMRよりも選択されるであろう。全体として、すべてのモデルとサンプルサイズにわたって表7と8の結果を比較すると、BLRTは、この論文で検討したすべてのインデックスとテストの中で最もよく機能する統計ツールである。
しかし、BLRTには欠点もある。BLRTを使用した場合、我々の例では計算時間が5倍から35倍に増加したことに留意されたい。BLRTアプローチのもう一つの欠点は、分布とモデルの仮定に依存することである。複製されたデータセットは、推定されたモデルから生成され、モデルで使用されたものと正確な分布を持つ。そのため、モデルや変数の分布に誤った仕様があると、複製されたデータセットが元のデータセットと類似した性質を持たなくなり、誤ったp値の推定につながる。例えば、あるクラス内のデータが歪んでいるにもかかわらず、正規分布としてモデル化されている場合、BLRTのp値は不正確となる可能性がある。外れ値もまた、誤ったp値の推定につながる可能性がある。さらに、BLRTは現在、複雑な調査データに対応することができない。同様に、様々なICは、モデル、分布、サンプリングの仮定に依存する。一方、LMRはパラメータ推定値の分散に基づいており、様々なモデルや分布の仮定の下で頑健かつ有効であり、複雑な調査データに対応することができる。したがって、このような状況では、LMRが望ましいと考えられる。しかし、我々のシミュレーションは検定の頑健性を評価しておらず、このトピックに関するさらなる研究が必要である

実践のための推奨事項

BLRTの計算時間が長くなるため、モデル探索の最初のステップではBLRTを要求しない方がよいかもしれない。その代わりに、BICとLMRのp値をガイドとして使って可能な解に近づき、いくつかのもっともらしいモデルが特定されたら、これらのモデルをBLRTを要求して再分析することができる。さらに、尤度の値は標準的な差の検定には使えないことが知られているが、尤度の実際の値は、次のようにクラス数を決定するための探索的診断ツールとして使用することが可能である。

図5は、クラス数の異なるモデルについて、クラス数を決定するための説明変数として利用できる対数尤度値のプロットを4つ表示したものである。左上図は、8項目のLCAで、カテゴリ別の結果(n = 1,000)を見てみると、2クラスから3クラスに移るときに尤度が大きく増加し、3クラスから4クラスに移るときにも尤度が大きく増加するパターンがあることがわかる。そして、4クラスから5クラスに移るとき、また5クラスから6クラスに移るとき、同様に平坦になる。4クラスから5クラスの間で線が平らになっているのは、4クラスから5クラスになったときに、可能性が実質的に増加しないことを示唆している。さらに、8項目のモデルについては、それが真のk = 4クラスモデルであることが分かっているので、正しいポイントで起こる平坦化が観察された。ほぼすべてのモデリング設定において、同様のパターンが見られた。10項目、n = 200のモデルは、他のモデルほど劇的なフラットニングはありませんが、先に述べたように、これは最も難しいモデリング設定である。これは、GMMモデルのプロットでも観察される。4つのプロットしかないが、一般的な知見として、対数尤度プロットはn = 500とn = 1,000のLCAモデルに対してかなり一貫して正しいモデルを同定することが示唆される。従来の方法では、これらの対数尤度の値を差の検定のために検定できないことは分かっているが、このプロットは、クラス数を探索する際に尤度を記述的なツールとして使用する方法である。
本試験のモデルの選択により、特定のモデルおよびサンプルの属性が結果に及ぼす影響について結論が得られなかったことに留意することが重要である。例えば、単純構造モデルのような特定のタイプのモデル構造に対するこれらのテストとインデックスの性能について一般化することはない。むしろ、モデルやサンプルの選択は、LRTやICの性能を理解するために、様々な混合モデルを探求したいことが動機となっている。今後の研究では、モデルの構造、結果の性質(カテゴリーか連続か)、項目数との相互関係、クラスの列挙のための指標の性能をよりよく理解することを目指すことができる。

結論と今後の方向性

本研究では、混合モデリングにおける正しいクラス数を特定するためのICおよび尤度ベースのテストの性能について検討した。本研究で検討したツールの中で、結果は、BLRTが他のツールより優れていることを示した。2番目に良かったのはBICで、次いで修正BICであった。先行研究では、クラス数を決定するための様々な適合指標やテストが検討されているが、この論文は、この種の混合モデルに対するBLRT法の性能を詳細に検討した最初のものの1つである。限られた数のFMAとGMMモデルだけでなく、カテゴリと連続の両方の結果を持つLCAモデルを考慮することで、ICと尤度ベースのテストの性能の理解を、これまで検討されてきたものよりも拡大することができる。しかしながら、これらの結果は単なるプレビューに過ぎない。BLRTの有用性について広く述べる前に、より広い範囲のLCAモデル(例えば、共変数を含む、異なる構造のモデル、連続とカテゴリカルの結果の組み合わせなど)に対するBLRTとその性能についてより多くの研究を検討する必要がある。とはいえ、これらの結果は、これらの指標の性能をより深く理解することで、混合モデルにおけるクラス数の決定方法に関するさらなる理解に貢献するものである。

相関係数の高い説明変数があれば多重共線性を考慮して片方を除くべき?

重回帰分析はとても分かりやすく有効な分析法だが、説明変数間の相関が高すぎる場合は、 パラメータの推定が不安定になるという問題点がある。 これは、説明変数間にすでに別の線型回帰関係が含まれているということであり、その意味でこのような現象を多重共線性という。多重共線性を回避する最も単純な方法は、 まず説明変数間の相関行列を見て、 相関がとても強いものがあれば、片方を説明変数から除くことである。 (p.183)

久々に読んでみたが、この記述はよろしくない。 前半は正しいが、後半の相関係数の記述は前半と整合的になっていない。

McElreath(2020)は次のように述べている。

科学文献の中には、多重共線性に対処するための様々な方法が紹介されている。その中で因果関係を考慮したものはほとんどない。実際に、モデルを当てはめる前に一対の相関を調べ、相関の高い予測因子を特定して削除するよう学生に教えている分野もある。これは間違いである。一対相関が問題なのではない。相関関係ではなく、条件付きの関連性が問題なのである。 そして、その場合でも、何をするのが正しいかは、何が共線性を引き起こしているかによる。データ内の関連性だけでは、何をすべきかを決めることはできない。(Statistical Rethinking-E-Book: A Bayesian Course with Examples in R and STAN (Chapman & Hall/CRC Texts in Statistical Science) 2nd Edition, Kindle Edition (English Edition))

相関係数を求めるのではなく、ペアプロットを描画するというのが正しい対策である。

より詳しい説明は以下の2つのエントリを参照のこと。

ides.hatenablog.com

ides.hatenablog.com

多重共線性を授業で習った記憶はないのでどのように教えられているかはわからないが、誰かに伝えるのは難しいな、とは思う。

多重共線性のシミュレーション

下記エントリーの続き。

ides.hatenablog.com

こちらの教科書から多重共線性について

  • Richard McElreath - Statistical Rethinking_ A Bayesian Course with Examples in R and STAN

6.1. 多重共線性

一般に、回帰モデルに追加する潜在的な予測変数が多くあることは事実である。たとえば、霊長類のミルク・データの場合、我々が結果として選ぶどの列も予測するために利用可能な7つの変数がある。なぜ、7つすべてを含むモデルを構築しないのか? いくつかの危険性がある。

まず、最も心配のない「多重共線性」から始めよう。多重共線性とは、2つ以上の予測変数の間に非常に強い関連があることを意味する。生の相関は重要ではない。むしろ重要なのは、モデル中の他の変数に対する条件付きで、関連性があることである。多重共線性の結果は,たとえすべての変数が実際には結果と強く関連していても,事後分布は,どの変数も結果と確実に関連していないことを示唆しているように見えるということである。 このイライラする現象は、重回帰がどのように機能するかの詳細から発生する。実際、多重共線性には何の問題もない。モデルは予測には問題なく機能するす。ただ、それを理解しようとするとフラストレーションがたまるだけである。多重共線性を理解すれば、一般的な回帰モデルをよりよく理解できるようになることが期待できる。

まずは簡単なシミュレーションから始めよう。その後、再び霊長類のミルクのデータに目を向け、実際のデータセットで多重共線性を見てみよう。

6.1.1. 多重共線性 脚

個人の足の長さを予測変数として、その人の身長を予測しようとすることを想像しよう。確かに身長は脚の長さと正の相関がある、少なくとも我々のシミュレーションではそうだと仮定する。しかし、両足(右足と左足)をモデルに入れると、やっかいなことが起こる。

以下のコードは、100人の身長と脚の長さをシミュレートしている。それぞれについて、まず身長がガウス分布からシミュレートされる。次に、各個体の脚の高さの比率を0.4から0.5までの範囲でシミュレートする。最後に、それぞれの脚は、実際の集団で典型的に見られるように、左右の脚の長さが全く同じにならないように、若干の測定誤差や発育誤差を加味してある。

最後に、このコードは身長と2本の足の長さを共通のデータフレームに入れる。

N <- 100   # 個体数 
set.seed(909)
height <- rnorm(N,10,2)  # 各脚の高さの合計をシミュレートする
leg_prop <- runif(N,0.4,0.5) # 脚長比 
# 左脚を比例+誤差としてシミュレート
leg_left <- leg_prop*height +    rnorm( N , 0 , 0.02 )
# 右足を比率+誤差としてシミュレート
leg_right <- leg_prop*height + rnorm( N , 0 , 0.02 ) 
# 結合してデータフレームにする
d <- data.frame(height,leg_left,leg_right)

さて、これらのデータを分析して、leg_leftとleg_rightの両方の予測変数で結果の高さを予測しようしかし、事後推定を行う前に、我々が何を期待するかを考えてみよう平均して、個人の脚は身長の45%である(このシミュレーション・データでは)。したがって、脚と身長の関連性を測るベータ係数は、平均身長(10)を平均身長の45%(4.5)で割ったあたりになると予想される。これは、10/4.5≒2.2です。では、代わりにどうなるか見てみよう。ここでは、これから起こることの責任が事前分布にないないことを確認するために、非常に曖昧で悪い事前分布を使う。

library(rethinking)
m6.1 <- quap(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right , a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) , br ~ dnorm( 2 , 10 ) , sigma ~ dexp( 1 )
  ) , data=d )
precis(m6.1)
      mean   sd  5.5% 94.5%
a     0.98 0.28  0.53  1.44
bl    0.21 2.53 -3.83  4.25
br    1.78 2.53 -2.26  5.83
sigma 0.62 0.04  0.55  0.69

この事後平均と標準偏差は狂っているように見える。これはprecis出力のグラフ表示の方が有用なケースである。なぜなら、事後平均と89%区間を一目見て、ここで何かが間違っていることがわかるように表示されるからである。

plot(precis(m6.1))

set.seedの行を省いて、もう何回かシミュレーションをしてほしい。両足の長さがほぼ同じで、身長が足の長さと強く関連しているのなら、この事後分布はなぜこんなに変なのだろうか?事後近似は正しく機能したのだろううか?

正しく動作している。そして、この事後分布は私たちが質問した内容に対する正しい答えである。問題は、その質問である。重回帰は質問に答えるものであることを思い出してほしい他の予測変数がすべてわかっている状態で、各予測変数を知ることの価値は何だろうか? 従って、この場合、質問は次のようになる。各脚の長さを知っていることの価値は何か?

post <- extract.samples(m6.1)
plot( bl ~ br , post , col=col.alpha(rangi2,0.1) , pch=16 )

図6.2. 左:各脚と身長の関連性の事後分布(モデルm6.1より)。両変数はほぼ同じ情報を含むので、事後は負の相関のある値の狭い隆起である。右図。2つのパラメータの和の事後分布は、どちらかの脚と身長の適切な関連性で中央化される。

その結果、図6.2の左側に示すようなプロットが得られた。\ この2つのパラメーターの事後分布は非常に高い相関を持ち、blとbrのもっともらしい値はすべて狭い尾根に沿って横たわっている。blが大きければ、brも小さくなるはずである。ここで起こったことは、両方の脚の変数がほとんど同じ情報を含んでいるので、もしあなたがモデルに両方を含めることにこだわるなら、同じ予測を生み出すblとbrの組み合わせが実質的に無限に存在することになる、ということなのである。

この現象を、このモデルに近似させたと考えるのも一つの方法である

 y_{i} \sim \operatorname{Normal}(\mu i, \sigma)

 \mu_{i} = \alpha + \beta_1 x_i  +  \beta_2 x_i

変数yは例の身長のような結果で、xは例の脚の長さのような1つの予測変数であるここで,x は2回使われており,これは両方の脚の長さを使うことによって起こる問題の完璧な例である.ゴーレムの視点からは,μi のモデルは次のようになる.

 \mu_{i} = \alpha+ \left(\beta_1 + \beta_2 \right) x_i

私が行ったのは、各項からxiを取り除くことだけである。パラメータβ1 とβ2は別々に平均値µに影響を与えることはないので、分離することはできない。

このシミュレーションの例では,事後分布がまさにそうなっている以下は、それらの和の事後分布を計算し、それをプロットする方法である。

sum_blbr <- post$bl + post$br
dens( sum_blbr , col=rangi2 , lwd=2 , xlab="sum of bl and br" )

結果の密度プロットが図6.2の右側に示されている。事後平均は2強と右肩上がりで、標準偏差は和の成分であるblとbrのどちらよりもずっと小さくなっている。脚長変数の片方だけで回帰を当てはめても、ほぼ同じような事後平均が得られる。

m6.2 <- quap(
alist(
height ~ dnorm( mu , sigma ) , mu <- a + bl*leg_left,
a ~ dnorm( 10 , 100 ) , bl ~ dnorm( 2 , 10 ) , sigma ~ dexp( 1 )
) , data=d )
precis(m6.2)
      mean   sd 5.5% 94.5%
a     1.00 0.28 0.54  1.45
bl    1.99 0.06 1.89  2.09
sigma 0.62 0.04 0.55  0.69

その1.99はsum_blbrの平均値とほとんど同じである。
基本的な教訓は、これだけです。2つの予測変数が非常に強く相関しているとき(モデル中の他の変数の条件付きで)、モデルに両方を含めると混乱が生じるかもしれないこのような場合、事後分布は間違っていない。これらのデータでは事後分布は間違っていないのである。そして、「答えられない」と言うことは、モデルにとって素晴らしいことなのである。そして、もしあなたが予測に興味があるだけなら、この足のモデルは素晴らしい予測をすることがわかる。ただ、どちらの脚がより重要であるかは主張しないのである。

この脚の例は明確でキュートである。しかし、これは純粋に統計的なものでもある。私たちはここで深刻な因果関係を問うているのではない。次はもっと因果関係のある面白い例で試してみましょう。

多重共線性 ミルク

脚の長さの例では、両脚をモデルに含めることが少し愚かであることは容易に理解できる。しかし、実際のデータセットで生じる問題は、相関の高い予測変数の間の衝突を予期していないかもしれないということである。そして、それゆえ、我々は、どちらの予測変数も重要でないと言うように事後分布を誤って読んでしまうかもしれない。このセクションでは、実際のデータでこの問題の例を見る。

先ほどの霊長類の乳のデータに戻ろう。

library(rethinking)
data(milk)
d <- milk
d$K <- standardize( d$kcal.per.g )
d$F <- standardize( d$perc.fat )
d$L <- standardize( d$perc.lactose )

この例では、perc.fat (脂肪率) と perc.lactose (乳糖率) という変数に注目するこれらを使って、総エネルギー量である kcal.per.g をモデル化する。上のコードでは、これら3つの変数をすでに標準化している。あなたは、多重共線性の自然なケースを探索するために、これらの3つの変数を使用することになる。これらの列には、NAという欠損値がないので、ここでは完全なケースを抽出する必要はないことに注意してほしい。しかし、quapはlmのような無謀な関数とは異なり、黙ってケースを落とすようなことは決してありませんので安心してほしい。
まず、kcal.per.g を perc.fat と perc.lactose の関数として、2つの二変量回帰でモデル化することから始めまよう。これらの事前分布については第5章(147ページ)を参照。

# 6.9 kcal.per.g regressed on perc.fat
m6.3 <-quap(
  alist(
    K ~dnorm(mu,sigma),
    mu <-a+bF*F,
    a ~dnorm(0,0.2),
    bF ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,data=d)

# kcal.per.g regressed on perc.lactose
m6.4 <-quap(
  alist(
    K ~dnorm(mu,sigma),
    mu <-a+bL*L,
    a ~dnorm(0,0.2),
    bL ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,data=d)

precis( m6.3)
precis( m6.4)
      mean   sd  5.5% 94.5%
a     0.00 0.08 -0.12  0.12
bF    0.86 0.08  0.73  1.00
sigma 0.45 0.06  0.36  0.54


       mean   sd  5.5% 94.5%
a      0.00 0.07 -0.11  0.11
bL    -0.90 0.07 -1.02 -0.79
sigma  0.38 0.05  0.30  0.46

bF と bL の事後分布は、基本的に互いに鏡像である。bFの事後平均は、bLの平均が負であるのと同様に正である。両方とも、ほとんどゼロの片側にある狭い事後分布である。それぞれの予測因子と結果の強い関連性を考えると、両方の変数が、種を問わず、乳中の総エネルギーの信頼できる予測因子であると結論づけられるかもしれない。脂肪が多いほど、牛乳のキロカロリーは多くなる。 乳糖が多いほど、牛乳のキロカロリーは少なくなる。しかし、両方の予測変数を同じ回帰モデルに入れたらどうなるかを見てみよう。

m6.5 <- quap(
  alist(
    K ~ dnorm( mu , sigma ) , mu <- a + bF*F + bL*L ,
    a ~ dnorm( 0 , 0.2 ) ,
    bF ~dnorm(0,0.5),
    bL ~dnorm(0,0.5),
    sigma ~dexp(1)
  ) ,
data=d )
precis( m6.5)
       mean   sd  5.5% 94.5%
a      0.00 0.07 -0.11  0.11
bF     0.24 0.18 -0.05  0.54
bL    -0.68 0.18 -0.97 -0.38
sigma  0.38 0.05  0.30  0.46

このとき、bFとbLの事後平均は0に近くなっている。また、両パラメータの標準偏差は二変量モデル(m6.3、m6.4)のときの2倍になっている。 これは脚の例と同じ統計的現象である。何が起こったかというと、変数perc.fatとperc.lactoseが同じ情報を多く含んでいることである。これらは、ほとんどお互いに代替しているのである。その結果、両方を回帰分析に含めた場合、事後分布はbFとbLの組み合わせの長い尾根を記述することになり、等しく妥当であることがわかる。 脂肪と乳糖の場合、この2つの変数は基本的に1つの変動軸を形成している。これを見る最も簡単な方法は、ペアプロットを使うことである。

pairs( ~kcal.per.g+perc.fat+perc.lactose,data=d,col=rangi2)

このプロットを図 6.3 に示す。対角線に沿って、変数がラベル付けされている。対角線から外れた各散布図では、縦軸の変数は同じ行にラベル付けされた変数で、横軸の変数は同じ列にラベル付けされた変数である。 例えば、図6.3の最初の行の2つの散布図は、kcal.per.g(垂直)対perc.fat(水平)、そしてkcal.per.g(垂直)対perc.lactose(水平)である。 脂肪の割合は結果と正の相関があり、乳糖の割合は結果と負の相関があることに注意。次に、中央の列の一番右の散布図を見てほしい。このプロットは、乳糖(水平)に対する脂肪分(垂直)の散布図です。点がほぼ直線上に並んでいることに注目してください。 この2つの変数には負の相関があり、非常に強い相関があるため、ほとんど冗長になっています。どちらかがkcal.per.gを予測するのに役立ちますが、もう一方が分かっていれば、どちらもそれほど役に立たない。

科学文献の中には、多重共線性に対処するための様々な方法が紹介されている。その中で因果関係を考慮したものはほとんどない。実際に、モデルを当てはめる前に一対の相関を調べ、相関の高い予測因子を特定して削除するよう学生に教えている分野もある。これは間違いである。一対相関が問題なのではない。相関関係ではなく、条件付きの関連性が問題なのである。 そして、その場合でも、何をするのが正しいかは、何が共線性を引き起こしているかによる。データ内の関連性だけでは、何をすべきかを決めることはできない。
ミルクの例で起こっていることは、哺乳類の母親が従わなければならない、ミルクに関するトレードオフが存在することです。ある種が頻繁に授乳する場合、そのミルクは水っぽく、エネルギーが低い傾向がある。このようなミルクには糖分(ラクトース)が多く含まれている。そのかわり、授乳回数が少ない種は、母乳は高エネルギーである必要がある。そのようなミルクは脂肪分が非常に高い。これは、次のような因果関係を意味する。

中心的なトレードオフは、ミルクの濃さ(D)を決めることである。この変数は観測されていないので、丸で囲んである次に脂肪Fと乳糖Lが決定される。最後に、F と L の組み合わせでキロカロリー K が決まる。もし D を測定できたり、種の他の側面に基づいて Dを予測する進化的・経済的モデルがあれば、回帰分析でつまずくよりましだろう。
多重共線性の問題は、モデルの適合に関する問題の一つで、非同定性(non-identifiability)と呼ばれるものである。あるパラメータが識別不可能な場合、データとモデルの構造上、そのパラメータの値を推定することができないことを意味する。 この問題はモデルのコーディングの間違いに起因することもあるが、多くの重要なタイプのモデルは、たとえ完全に正しくコーディングされていても、識別不可能な、あるいは弱く識別可能なパラメータを有している。自然は、たとえモデルが正しくても、我々に簡単に推論する義務を負ってはいないのである。
一般的に、利用可能なデータが目的のパラメータについて多くの情報を含んでいるという保証はない。 そのような場合、ベイズ計算機は事前分布に非常に近い事後分布を返す。 したがって、事後分布を事前分布と比較することは、モデルがデータからどの程度の情報を抽出したかを確認する良い方法となる。事後分布と事前分布が似ているということは、計算が間違っているということではない。しかし、もっと良い質問をするきっかけになるかもしれない。

Rethinking: 識別は保証する。理解はあなた次第。技術的に言えば、ベイズモデルにおいて同定性は問題ではない。なぜなら、事後分布が適切である限り、つまり、1へ積分する限り、すべてのパラメータが同定されるからである。しかし、この技術的な事実は、事後分布の意味を理解できることを意味するものではない。だから、ベイズ的な文脈では、パラメータの同定が弱いという言い方をした方がいいかもしれない。しかし、この違いは技術的なものに過ぎないかもしれない。実は、DAGで因果関係が特定できるはずだと言われても、統計的に特定できない場合がある。私たちは、デザインを考えるのと同じように、統計にも力を入れなければならないのである。

Overthinking: 共線性のシミュレーション。2つの予測変数の間の関連性によって事後分布の不正確さがどのように増加するかを見るために、シミュレーションを使おう。以下のコードは、相関のある予測変数を生成し、モデルを当てはめ、perc.fatとkcal.per.gを関連付ける傾きの事後分布の標準偏差を返す関数を作成する。

library(rethinking)
data(milk)
d <- milk
sim.coll <- function( r=0.9 ) {
  d$x <- rnorm( nrow(d) , mean=r*d$perc.fat , 
    sd=sqrt( (1-r^2)*var(d$perc.fat) ) )
  m <- lm( kcal.per.g ~ perc.fat + x , data=d )
  sqrt( diag( vcov(m) ) )[2] # stddev of parameter
}
rep.sim.coll <- function( r=0.9 , n=100 ) {
  stddev <- replicate( n , sim.coll(r) ) 
  mean(stddev)
}

r.seq <- seq(from=0,to=0.99,by=0.01)
stddev <- sapply( r.seq , function(z) rep.sim.coll(r=z,n=100) )
plot( stddev ~ r.seq , type="l" , col=rangi2, lwd=2 , xlab="correlation" )

r.seqの各相関値に対して、このコードは100回の回帰を生成し、それらからの平均標準偏差を返す。このコードでは、暗黙の平行事前分布を使用しているが、これは悪い事前分布である。そのため、共線性変数の効果を誇張してしまう。情報量豊富な事前分布を用いると、標準偏差の増大はより緩やかになる。

多重共線性の根本的な問題はモデルが答えようとする問題の方にありモデルそのものにはない

多重共線性があるかどうかを調べるのにVIFの値を出して4以下/10以下ならOKといった運用をしがちだが、この運用には問題がある。
第一にVIFの値で機械的にカットオフを設けることに問題があること、第二に多重共線性をモデルに求めているという根本的な誤りが存在している。

ResearchGateの質疑より。

https://www.researchgate.net/post/Multicollinearity_issues_is_a_value_less_than_10_acceptable_for_VIF

Alejandro Ros-Gálvez
多重共線性の問題:VIFの値は10未満でよいのか?

Roman Mathias
しかし、私の考えでは、「VIFはXYZ以下であるべきだ」というような儀式的な経験則はすべて、Gigerenzer(2004)に倣って「頭の悪い統計」を行うよう人々を誘うので、潜在的に危険なものである。 回帰モデル中の2つ以上の予測変数がより強く相関している場合、対応する回帰係数の標準誤差は大きくなる。これはモデルの問題である必要はなく、予測変数は本質的に同じ情報を含んでいるので、モデルが応答との共有分散をどちらかに帰することができないことを思い出させるだけなのである。大きな標準誤差は、片方の係数を増やしてもう片方を減らしても、共有された情報のためにほとんど同じ結果が得られることを反映している。
高いVIFのために相関のある予測変数のどちらかをモデルから取り除くと、モデルは実際には良くならず、残りの回帰係数の標準誤差が小さくなるだけである。
相関のある予測変数のいずれかを操作したときに応答が変わるかどうかを実際に見つけるには、操作実験や変数間の因果関係についての専門知識が必要である。一方、回帰係数の高い標準誤差は、ある程度まで予測変数のどちらかがレスポンスの分散を等しく説明できることを示すだけであることを念頭に置いていれば、VIFが6であると報告しても、私には全く問題ないように思われる。根本的な問題は、通常、モデルが答えようとする問題であって、モデルそのものにあるのではない。
McElreath (2020)の6.1.1章では、両足の長さで身長を予測しようとするとどうなるか、明らかに共線性の問題があることが示されている。

  • Gigerenzer G. 2004. Mindless statistics. The Journal of Socio-Economics 33: 587–606.
  • McElreath R. 2020. Statistical rethinking: A Bayesian course with examples in R and Stan. CRC Press.

スマートフォン・アディクション・インベントリー Smartphone Addiction Inventory (SPAI)

journals.plos.org

  • Lin, Y.-H., Chang, L.-R., Lee, Y.-H., Tseng, H.-W., Kuo, T. B. J., & Chen, S.-H. (2014). Development and Validation of the Smartphone Addiction Inventory (SPAI). PLoS ONE, 9(6), e98312. https://doi.org/10.1371/journal.pone.0098312

SPAIの開発

物質関連障害とインターネット中毒の経験を持つ2人の精神科医、LinとChangは、26項目のChen Internet Addiction Scale(CIAS)を「スマートフォン嗜癖」評価用に修正した。CIASの修正版の心理測定研究は、LinがChangの許可を得て行い、探索的因子分析により5つの下位尺度が同定された[5]。「インターネット」という用語は「スマートフォン」に変更された。北京語版の尺度は、専門家委員会によって最終的に修正された。最終的な改訂は以下の通りである。(1) 項目4と6は、12項目のProblematic Cellular Phone Use Questionnaire [10] の意味的に類似した項目2と3に置き換えられた。これは、元の項目が単に「スマートフォン使用」を「インターネット使用」の代わりに使うだけでは意味を成さないためである(2)。また、スマートフォン利用の特殊性から、項目21の「道路横断時にスマートフォンを見る、運転中や待ち時間にスマートフォンをいじり、危険にさらされる」を尺度の末尾に追加した(3)。項目 23 では、元の「睡眠時間を削ってネットする時間を増やすことが習慣になっている」を「スマートフォンを使うことが習慣になり、睡眠の質と総睡眠時間が減った」として文章を修正した。(4)項目25については、"インターネットをするため、いつもの時間に食事ができない "を "スマートフォンの携帯性 "という従来のパソコンによるインターネット利用とは異なる特性に合わせて修正したものである。SPAIは、26点から104点までの範囲で、1=強く反対、2=やや反対、3=やや賛成、4=強く賛成の4段階で評価するよう被験者に求めた。

日本語

日本語翻訳は見つけることはできなかった。

CIAS

Chen Internet Addiction Scale(CIAS)がベースになっている。CIASはこちら。

https://psycnet.apa.org/record/2004-10292-005

ゲーム依存症 低年齢化進む…小4男子の割合高いという調査結果も(ヨミドクター)

yomidr.yomiuri.co.jp

塩島祐子さんの執筆生地。山田さんも指摘のとおり、誤りが多くあり、あまりにも適当すぎる。