井出草平の研究ノート

Rで平行分析

Mplusでの平行分析の実行の仕方は以前に書いたが、Rでも実行できる。

今回もbfiデータの1~25列目を使用する(参考: http://ides.hatenablog.com/entry/2019/03/19/093726)。

library("psych")
data(bfi)
d1<-bfi[1:25]

平行分析は次のように指定する。

fa.parallel(d1) 

"Parallel analysis suggests that the number of factors = 6 and the number of components = 6"と出力され6因子、6成分と提案があった。 今まで行ってきた他の方法よりも因子数が多い印象である。

f:id:iDES:20190722165507p:plain

ちなみに図にあるPCとは" principal components"のことで、主成分分析FAはprincipal axis factor analysisで主因子法のことである。

オプションとしてfa="minres"のように表記して、minres(最小残差法), ml(最尤法), uls(重みづけのない最小二乗法), wls(重み付き最小2乗法), gls(一般化最小二乗法), pa(主成分)が指定できる。デフォルトはminresである。

対角SMC平行分析はオプションで指定を行う。

fa.parallel(d1, SMC=TRUE) 

こちらは8因子を提案がされてきた。

堀啓造先生によると下記のような特徴があるらしい。

対角に SMC 入れた分析(SMC-EIGEN)の場合、ほかの推定法に比べ極端に多くの因子数を推測する。
対角SMC の平行分析はその点それ以上の因子をとってはいけないという指標となりうるのである。
欠点としてはマイナー因子も拾うことがあるので,まったく因子を想定しない項目群では多めの因子数を推定する。
https://www.ec.kagawa-u.ac.jp/~hori/yomimono/nfac.ppt

たしかに今回のシミュレーションでも最も多い因子数を提案してきた。

因子負荷と相関係数の違い

因子分析をしていて混同しがちなのは因子負荷(因子パターン)と相関係数ではないかと思う。

誰も教えてくれなかった因子分析: 数式が絶対に出てこない因子分析入門

誰も教えてくれなかった因子分析: 数式が絶対に出てこない因子分析入門

この2つを直観的に区別するのに役に立つのが松尾・中村(2002)にある図2.5.6である。

f:id:iDES:20190719201919p:plain

これは2因子を想定、斜交回転をした図である。 第1因子と第2因子は45°くらいで交わっている。

因子負荷は軸と並行に下ろしてきたところの数値である。 一方で、相関係数は軸と90°直角に下ろしてきたところの数値である。

このため、因子負荷と相関係数は異なった値となわけだが、直交行回転の場合は同じ値となる。

因子負荷とは技法的には相関係数ではあるが、実際に相関係数とは別の値になることがこの図でわかるのではないだろうか。この図も含めてこの本は直観的に因子分析が理解できる良書だと思う。

MAP(最小平均偏相関)

探索的因子分析において因子数を決定する基準として使われるMAPについて。

MAP:Minimum Average Partial correlationであり、日本語だと最小平均偏相関になる。
Velicerによつて開発された方法である。

link.springer.com

因子数を決める基準はどれが適切かという議論は比較的多くされているようだ。

下記の論文ではMAPについて次のように書かれている。

http://www.glmj.org/archives/articles/Pearson_v39n1.pdf

試験1では、VelicerのMAPはある程度良好であった;しかし、PA-R、BICLRTほどではありませんが、BICPA-Rと同様に、変動率が小さい場合には過小評価される可能性があった。また、過小評価はBICPA-Rよりも変動率に影響をうける。その傾向は研究2でより厳しく、この知見はZwickとVelicer(1986)のそれと類似していた。 したがって、MAPを用いている研究者は、特に変動因子比が低い場合には、この方法が示唆するよりも1つまたは2つ多い因子を用いたモデルを試した方が良いかもしれません。一般的な統計パッケージではMAPも並行分析も標準的な選択肢ではないため、これらのいずれか外部ツールを使うことになる。分析者は代わりに共通因子並行分析を試みるべきである。

この論文ではAICの結果が良かったとされている。
少なくとも、AIC>BICであるのは確かなようだ。

ともあれRを用いたMAPの出し方を書いておこう。
psychパッケージのVSS()で計算することができる。

personality-project.org

psychパッケージでのMAPの求め方

library(psych)
data(bfi)
d1 <- bfi[1:25]

IPIP-NEOのデータを使用する。IPIP-NEOの詳細はこちら 5因子5問ずつのデータを使うので1~25列目までをd1に格納した。

res01<-vss(d1)
print(res01,digits =4) 

偏相関の値がわりと同じような数値になりやすいため、小数点以下を4桁に調整してみた。

VSSはどれが提案されたのかがちゃんと文章で出るようになっている。

The Velicer MAP achieves a minimum of 0.0148 with 5 factors

ということで5因子が提案されたようだ。

追記:

MAP基準の推定方法を指定できるようだ。

vss(data, fm="minres")  #最小残差法
vss(data, fm="mle") #最尤法
vss(data, fm="pc") #主成分

www.rdocumentation.org

Mplusで探索的因子分析を行う

Rでの方法はこちらを参照。

ides.hatenablog.com

データはRのpsychパッケージに含まれるものを利用するので、作成方法は下部に記載する。 (追記: 2021/12/06) データの内容についてはこちらを参照。

Mplusのコード

TITLE:
    Exploratory Factor Analysis by Mplus, using IPIP-NEO data.

DATA:
    FILE = "bfi.dat";

VARIABLE:
    NAMES = A1-A5 C1-C5 E1-E5 N1-N5 gender education age;
    USEVARIABLES = A1-A5 C1-C5 E1-E5 N1-N5;
    MISSING = .;

ANALYSIS:
    TYPE = EFA 3 6;
    ESTIMATOR = ML;
    ROTATION = oblimin;
    PARALLEL = 50;

PLOT:
   TYPE = plot2;

OUTPUT:
    residual;

DATA:
データを置いた場所を指定する。日本語のフォルダ内に入れない方が良い。データのエンコードUTF-8にしておこう。

VARIABLE:
今回使うのはA1からN5の25個のデータ。

ANALYSIS:
EFAとはExploratory Factor Analysisの省略。
3 6となっているのは、3因子モデルから6因子モデルを検討するという意味。
推定法はML、回転はオブリミオン回転を指定している。回転を指定しない場合のデフォルト値はオブリミオン回転なので、指定はしなくてもよい。
PARALLELは平行分析のサンプリング数。この数を大きくすると、計算が長くなるのでほどほどにしておこう。Mplusのマニュアルでは50回と例示されてある*1が、 Finch & Bolin(2017)では1000という例示がある(p.177)のでマシンパワーがあれば大きい数字でもよいのだろう。

PLOT:
スクリープロットの出力を指定。

まず因子数を決めることから始める。 とはいっても、このデータの妥当な因子数はよくわからない。

ガットマン基準

RESULTS FOR EXPLORATORY FACTOR ANALYSIS


           EIGENVALUES FOR SAMPLE CORRELATION MATRIX
                  1             2             3             4             5
              ________      ________      ________      ________      ________
                4.729         2.581         2.019         1.468         0.979

固有値1を超えるのは4因子モデルである。5因子モデルになると0.979と1を切るため、固有値を1以上とするガットマン基準からは4因子モデルが提案される。

モデル改善の検定

モデルのフィッティングを見てみよう。

SUMMARY OF MODEL FIT INFORMATION


                                  Degrees of
     Model           Chi-Square    Freedom     P-Value

     3-factor          2794.270       133       0.0000
     4-factor          1376.733       116       0.0000
     5-factor           903.210       100       0.0000
     6-factor           649.008        85       0.0000

                                               Degrees of
     Models Compared              Chi-Square    Freedom     P-Value

     3-factor against 4-factor      1417.537        17       0.0000
     4-factor against 5-factor       473.523        16       0.0000
     5-factor against 6-factor       254.203        15       0.0000

Models Comparedのところをみると、モデルが改善されたかがわかる。改善しなくなったところで打ち止めをするのが通常である。具体的にはカイ二乗検定のP値をみる。例えば、6因子モデルのP値が0.360とかになっていると、モデルが改善していないとして、一つ前の5因子モデルを採用する。

ただ、実際には6因子モデルも0.000なので、5→6因子モデルにすることによってフィッティングは改善していることになっている。

追記(2021/12/05)
Bengt O. Muthenは下記の記述があった。 http://www.statmodel.com/discussion/messages/8/26371.html?1539133417

このような矛盾したメッセージは、実際のデータでもよく起こる。なぜなら、データは、ある数の因子を持つ完全な因子モデルを表していない可能性があるからである。例えば、m個の因子があっても、いくつかの残差相関がある場合がある。この場合、カイ二乗はmよりも多くの因子を持つべきだとでますが、これは本当のモデルではない。これは、いくつかの因子が2つの有意な負荷量しか持たないことでよく現れる。残差相関が存在するかどうかを確認するために、修正指数を見てほしい。一般的には、平行分析の方がこの問題に強いかもしれないが、これが徹底的に研究されているかどうかはわからない。そしてもちろん、どのような解決策がトピックの理論とどのように関連しているかを確認する必要がある。

MAP

Velicer(1976)のMinimum Average Partial(MAP)はMplusのオプションには含まれないようだ*2

MAPの値を出せるのは、服部先生の「忍者はっとりくん」、清水先生のHAD、RのpsychパッケージにあるVSS()関数のようだ。MAPのデータが必要な時には別のプログラムが必要である。

BIC

Model BIC
3-factor 184190.261
4-factor 182907.658
5-factor 182561.134
6-factor 182425.992

BICは小さいほうがフィッティングがよいとされている。微妙に数値が現象して行っているが、6因子が最も小さく、値も劇的に改善しているとも言いづらいため、あまり参考にはならなさそうである。

ちなみに、BICはプロマックス回転では出力されない。

スクリープロット

下の図はスクリープロットと平行分析の図である。

f:id:iDES:20190712174514p:plain

スクリープロットは推移がなだらかになる前までの数値を採用するので、4因子モデルが良さそうである。

平行分析

黒色の斜め横点線より上側の因子数を採用する方法である。清水先生のブログでちゃんと説明されているのでま詳しい解説はそちらで。

このデータの場合は4因子モデルが支持されている。

結論

基準を総合すると4因子モデルがよさそうである。
このデータはNEOなのでもともとは5因子モデルのデータであったはずなのだが、統計量は4因子モデルを指示している。

データ作成 追記: 2021/12/06

Rで作成する。

データ呼び出し

library("psych")
data(bfi)
df1 <- bfi

実数型に変数型を変更

df1$A1<-as.numeric(df1$A1)
df1$A2<-as.numeric(df1$A2)
df1$A3<-as.numeric(df1$A3)
df1$A4<-as.numeric(df1$A4)
df1$A5<-as.numeric(df1$A5)
df1$C1<-as.numeric(df1$C1)
df1$C2<-as.numeric(df1$C2)
df1$C3<-as.numeric(df1$C3)
df1$C4<-as.numeric(df1$C4)
df1$C5<-as.numeric(df1$C5)
df1$E1<-as.numeric(df1$E1)
df1$E2<-as.numeric(df1$E2)
df1$E3<-as.numeric(df1$E3)
df1$E4<-as.numeric(df1$E4)
df1$E5<-as.numeric(df1$E5)
df1$N1<-as.numeric(df1$N1)
df1$N2<-as.numeric(df1$N2)
df1$N3<-as.numeric(df1$N3)
df1$N4<-as.numeric(df1$N4)
df1$N5<-as.numeric(df1$N5)
df1$O1<-as.numeric(df1$O1)
df1$O1<-as.numeric(df1$gender)
df1$O1<-as.numeric(df1$education)
df1$O1<-as.numeric(df1$age)

この作業はMplusAutomation()がnumericでないとうけつけないため。

Mplusデータの書き出し

library(MplusAutomation)
variable.names(df1) # 変数名を書き出し
prepareMplusData(df1, filename="bfi.dat", overwrite=T)

ユースケ・サンタマリアのうつ病

ユースケ・サンタマリアはある時期、うつ病であったようだ。

文庫版でユースケ・サンタマリアのインタビューが追加されている。他のインタビューに答えている人はうつ病かよくわからないが、ユースケ・サンタマリアに関しては比較的重症のうつ病であったようだ。

ユースケ 病院に行くと異常なしで、「じゃあストレスだ」って心療内科とかに行かされるわけですよ。そうすると、欝じゃないのにそういうとこに行くことになっちゃった人って、先生からいろいろ言われたときに芝居しちゃうっていうか、そっちに引きずられちゃうんですよ。先生が僕の目とかじっと見ながらしゃべってるんだけど、ちょっと目線をキョドらせてみたりとか(笑)。人間のわけわかんないところなんだけど、その場に乗っかっちゃう。
吉田 なぜか期待に応えようとしちゃうわけですね(笑)。
ユースケ そう。で、やっぱり「ちょっと欝の気があるね」って抗欝剤とか出されちゃって、知り合いにその薬を見せたら「絶対に飲まないほうがいい。こんなもん飲んだらホントに病気になっちゃうよ」って言われて。だから飲んでないんだけど、そんな薬でこの症状が治まるんなら飲もうかな、みたいな。藁にもすがるような思いで、占い師みたいな人に診てもらったりね。宗教に走る人の気持ちがちょっとわかったというか。「わかります。誰でもそうなりますよ」とか言われると信じちゃうんだから。弱ってると信じちゃうんですよね。だからちょっと高い授業料を払って。で、ホントに悪かったのは五年ぐらいです、ゲロゲロで。

占い師には見てもらったが、抗うつ剤は飲まなかったらしい。

実際のところはわからないが、抗うつ剤を入れていたら、簡単に治っていたように思う。抗うつ剤を入れるとずっと飲まないといけないというのは誤りで、うつ病がなおっていれば、比較的簡単に抜くことができる。

ずっと抗うつ薬を飲んでいる人というのは、うつ病が治っていないから飲んでいるだけにすぎない。惰性で飲んでいる人もいるとは思うが。

抗うつ剤を飲んだら「ホントに病気になっちゃうよ」というのは、原因と結果が逆転している。病気だから薬を飲むのである。

抗うつ剤への抵抗はまだまだ大きいようだ。

うつ状態にあるときのエピソード。

ユースケ でも、その頃って周りみんな元気だから、「どうしたの?ユースケさん!美味い寿司屋あるから行こうよ!」って、とにかく俺を元気づけたいから。「気分悪いから寿司なんて一番食えないんだよ」って話なんだけど、でもそれ言うと向こうが「ああ・・・そう」みたいな感じで二度と連絡ないですから。「せっかく人が元気づけようとしてやってるのに断りやがった」「後輩のくせに」とか思われちゃう。俺は飯が食えないし、人と会ったりするテンションじゃないんだよ、動けないんだもん。
吉田 「そんなときは美味しい御飯食べれば元気が出るよ」と言われでも・・・。
ユースケ そう、「肉行こう肉!」みたいな。それが一番食えないのに。

よくある話。

宮本輝「パニック障害がもたらしたもの」

パニック障害になった有名人の一人として宮本輝がいる。

宮本輝のエッセイに「パニック障害がもたらしたもの」というものがある。

いのちの姿 完全版 (集英社文庫)

いのちの姿 完全版 (集英社文庫)

宮本輝パニック発作を初めて経験した時は次のように書かれてある。

もうじき淀駅に着くというころ、私は自分がいつもの自分とは違うことに気づいた。自分が自分でないような、神経の焦点が合っていないような、理由もなくかすかな恐怖感が心のなかでかすかに波立っているような、そんな感覚がつづいたあと、私は不意に全身が地の底に沈んでいきそうになって、慌てて電車の座席に両手を突いて支えた。
電車は何の支障もなく走りつづけている。何等かの事故が起こったのではなく、異変が生じたのはこの自分なのだとすぐにわかった。
その途端、どうにも抑えようのない、耐えられないほどの不安と恐怖がせりあがってきた。視界が白っぽかった。
自分にいったい何事が生じたのか、さっぱりわからなくて、とにかくいっときも早く電車から降りたいと思っているうちに、眩暈(めまい)と耳鳴りがして、不安感はさらに強まり、動悸(どうき)が頭の芯にまで響き始めた。全身から血の気が引き、掌が汗で濡れていくのも感じた。

電車の中でパニック発作を初めて経験したようだ。

翌日、私は会社の近くにある病院で診てもらった。
(中略)
「何か悩み事はない?」
「借金がありますねェ。秋に結婚するんですが、その費用を使い込んで、結婚相手にばれんうちに穴埋めしようと思って、きのう競馬場に行こうとしてたんです」

なかなかいい感じのクズっぷりである。

ここで終わると宮本輝をディすってるるようにも読めるので、一応全体の骨子を述べると、パニック障害になってサラリーマンが続けられなくなり、小説家になるしかないと考えた。必死に頑張って小説を仕上げて、世の中に評価されるようになったという感じの流れである。パニック障害があるから小説家として大成したという良い?話。

ひきこもりと暴力・犯罪

川崎の事件、元農水省事務次官による長男の刺殺事件を受けて、ひきこもりと暴力の関係性があるのではないかと考える人が多いようで、メディアの取材等も来ているので、少しまとめてみようと思う。

ひきこもりというものを有名にした一つは新潟少女監禁事件(wikipedia)である。事件が発覚した時に犯人である佐藤宣行の年齢は39歳であり、長期間ひきこもり状態にあったようである。ひきこもりと暴力・犯罪の問題はこの問題が社会問題化した当初から関心が持たれていた。

調査結果から見えること

ひきこもりの調査で暴力にいて質問したものがいくつかあるため、まずはそのデータを整理してみたい。

まず内閣府の若年(~39歳)ひきこもり調査(2010)では下記のような結果になっている*1。ランダムサンプリングであるため、日本全国の平均値が知ることができる調査である。

質問番号 質問文 ひきこもり 非ひきこもり
Q28-13 家族を殴ったり蹴ったり
してしまうことがある
5.1% 2.2%
Q28-14 壁や窓を蹴ったりたたいたり
してしまうことがある
10.2% 6.2%
Q28-15 食器などを投げて壊すことがある 8.5% 0.9%
Q28-16 大声を上げて怒鳴り
散らすことがある
8.5% 10.2%

これらの項目すべてが家庭内暴力にあるが、最も気になるのは対人暴力であろう。結果は5.1%である。この調査の一般人口(ひきこもり以外)の一般人口の対人暴力がどのくらいかという正確な値はわからない(対象者がウソをつくので正確な値は把握できない)が、だいたい倍近い結果となっている。

次に40~64歳を調査した内閣府のひきこもり(2019)調査*2見てみよう。

質問番号 質問文 ひきこもり 非ひきこもり
Q36-13 家族を殴ったり蹴ったり
してしまうことがある
4.3% 1.0%
Q36-14 壁や窓を蹴ったりたたいたり
してしまうことがある
8.5% 3.1%
Q36-15 食器などを投げて壊すことがある 4.3% 0.8%
Q36-16 大声を上げて怒鳴り
散らすことがある
19.1% 7.3%

似たような結果だが、ひきこもり群で大声を上げて怒鳴り散らすことがあるの割合が高くなっている。 対人暴力は4.3%と若年調査(~39歳)と大差がない値が出ている。

家族外への暴力の噴出

家族外の暴力を聞いている調査は知っている中では大分県の調査*3だけである。

これまであった問題行動では、家族用アンケートで「社会参加ができないこと以外に深刻な問題がない」が23人(41%)、「家族以外の暴力がある」0人であった反面、「家族への暴力」7人(13%)、「家庭内で物を壊す」12人(21%)と家族内の問題も少なくなかった。家族の調整や家庭内暴力に対する対応等の支援も必要である。

割合
家庭内暴力(対人) 12.5%
家庭内暴力(対物) 21.4%
家庭外への暴力 0%

家族以外の暴力は確認できなかった。

この調査はランダムサンプリングした一般人口を基にしたものではなく、相談機関を対象したものである*4。対象機関以下のものである。

内閣府の調査に比べて家庭内暴力の値が高い。おそらく相談機関で票を集めたためであろう。というのは、比較的深刻なケース(家庭内暴力があるなども含め)が相談機関に集まりやすいからである。

それでも、家庭外へ暴力が噴出することは確認できなかったようである。

少なくとも一般の犯罪率よりは低いというくらいは断言できそうである。

ひきこりはありふれた現象である

以前シノドスで書いたが、ひきこもりを経験したか、現在ひきこもり状態にある人は10人に1人いる。

現在、ひきこもり状態にある若者は1.6%であり、経験がある若者8.4%を合わせると9.7%なる(注2)。およそ10人に1人の若者が、過去にひきこもり経験があるか、現在ひきこもりであるということになる。現在ひきこもり状態にある若者は100人に1~2人程度であるが、経験者も含めると10人に1人という割合になる。少なくとも日本において、ひきこもりは決して稀な現象ではないのである。

f:id:iDES:20190605164952p:plain

自宅にいることが多いため、あまり目立たないが、ひきこもりは日本にあふれるほど多くいるのだ。ひきこもり経験も自ら進んで話す話題でもないため、耳にすることも少ないかもしれないが、特別な現象ではない。日本人の10人に1人に起きていることであれば、ひきこもりを特殊な集団として取り上げる意味はほとんどない。

自閉スペクトラム症と犯罪

全例把握できているわけではないが、メディアで報道される目立つ事件で、ひきこもり状態にあるものの犯罪は、知る限りすべて自閉スペクトラム症の疑いがある。最近起こった事件については別ブログのエントリを参照のこと。

ひきこりと犯罪についてのまとめ

  1. ひきこもり状態にある者が犯罪を犯すことはあるが、珍しい。調査でも確認できていない。
  2. ひきこもりはありふれた現象であり、異質で特別な集団と扱うのは誤り。
  3. ひきこもりによる犯罪は、ひきこもりの中の自閉スペクトラム症の疑いがあるグループに見られる。
  4. ひきこもりは一様なグループではなく、その中に犯罪との親和性が高いサブグループがあり、そのサブグループに対して対応をすべきである。
  5. そのサブグループとして現在判明しているのは自閉スペクトラム症の疑いがあるグループである。

*1:https://www8.cao.go.jp/youth/kenkyu/hikikomori/pdf_gaiyo_index.html

*2:https://www8.cao.go.jp/youth/kenkyu/life/h30/pdf-index.html

*3:大分県精神保健福祉センター、2004、『ひきこもりの実態調査報告書』

*4:調査の対象になったのは、保健所14カ所・市町村57カ所・福祉事務所(県・市)17カ所・児童相談所2カ所医療機関(精神科・心療内科)56カ所・教育事務所6カ所・高等学校57カ所・警察(フレンドリーサポートセンター)・フリースクール2カ所・精神保健福祉センター・教育センター各1カ所(調査対象機関209施設)