自閉スペクトラム症の原因の探索

自閉スペクトラム症関連の記事の紹介を2つほど。

blog.livedoor.jp

抗生物質が原因ではないかという論文が紹介されている。
場末P科病院の精神科医のblogがまた再開されるのかわからないが、再開するなら喜ばしいニュースだ。

日本語ではPPA(プロピオン酸)に関する論文の紹介がNewsweekでされていた。

www.newsweekjapan.jp

どちらの仮説の段階である。
ワクチン仮説のように否定される可能性も十分にあるので、何かの行動に移すのではなく、少し静観しておいた方が良いように思う。

ちなみに、現在時点で、自閉スペクトラム症の増加原因として、間違いなく確定しているのは、母親・父親のが高齢で出産するようになったことである。

食品添加物がない時代にも自閉スペクトラム症があっただろうという声が上がっているようだが、増加というのは促進要因が加わったという意味である。
遺伝的要因+環境的要因(受精から出産まで)でほとんどの精神障害は8割程度が説明できるのだが、環境要因に新たなものが加わったため、発病しやすくなったという理解が正しい。

統合失調症が過去は人口の1%程度あったが、1) 栄養状態の改善(妊娠後期の母体に十分な栄養が供給されるようになった)、2) 住環境の改善(特に冬に暖房で寒い部屋で母親が過ごさなくてもよくなった)ことによって、有病率は著しく減少した。遺伝的に素因があったとしても、発病するか否かは、(物理的な)環境によって意外と左右されるのである。

ひきこもりはロスジェネ世代に多いのか――ロスジェネ仮説を検証する

synodos.jp

シノドスに論考が掲載されました。
ロスジェネ仮説の検証を前面に押し出した方がいいというのは芹沢さんのアドバイスで、趣旨が明確になったように思います。

統計技法的には区間推定にagresti-coull法を初めて使いました。

今学期のおすすめ学習アプリ2019年秋

大阪大学でICT教育を取り上げている授業でのまとめ。2018秋、2019春と以前に2回エントリを入れているのでそちらも参照のこと。

ides.hatenablog.com

ides.hatenablog.com

スタディプラス

play.google.com

今までなぜか登場していなかったアプリ。超有名。
受験勉強の際に自身の勉強時間の管理ができる。臨床心理学でいう「認知療法」である。一般的には「見える化」などの一部だろうか。

また、同じ志望校の人たちがどのくらい勉強しているかなどもわかるので、合格するための勉強時間のベンチマークとして使えるアプリである。受験生にオススメ。
大学になったら要らないかも。

HiNative

hinative.com

onna と onna no hito はどう違いますか?
説明が難しい場合は、例文を教えて下さい。

という質問に対して下記の回答が返ってくる。

onna = woman

no = of
hito = human| people | person
onna no hito = a polite word for "woman"

【onna】 と 【onna no hito】 はどう違いますか? | HiNative

知恵袋などよりもよほど信頼性が高い、とのこと。

Spotify

www.spotify.com

言わずと知れた超有名アプリだが、英語学習用のポッドキャストがあるらしい。

www.fluentu.com

これはかなり興味深い。

Reverso

www.reverso.net

Google翻訳的なもの。スマホのアプリも存在している。
Google翻訳は言語によってはダメダメなので、Reversoが役立つときもあるとのこと。
機械翻訳は毎日使っているので、試してみたい。

DROPS

languagedrops.com

多言語に対応した言語習得アプリ。おもしろそうなので少し使ってみたい。

アラビア語検索エンジン アラジン

アラビア語専攻の学生が2名いたため、アラビア語関係のアプリ・サイトがいくつかあがった。

www.linca.info

アラジンというのは、オンラインのアラビア語-日本語の辞典のようである。

アラジンのすばらしさは、アラビア語との関係がない人にはあまりピンとこないので、こちらのブログエントリを参照してもらうとわかるかもしれない。

blog.goo.ne.jp

日本語-英語は学習環境が整っているが、アラビア語となると、あまり環境が整っていないので、いろいろと苦労するのだろう。

Arabic Almanac

Arabic Almanac - Hans Wehr

Arabic Almanac - Hans Wehr

  • Omar Jahangir
  • 教育
  • 無料
apps.apple.com

こちらも辞書。詳しくないので不確かだが、アプリの説明によるとアラビア語の有名な辞書Hans Wehr辞書のオンライン版だそうだ。アラビア語を勉強する人にとっては、紙の辞書を持ち歩かなくてもいいので、かなり便利だろう。


僕のおすすめ

前回、前々回として出しているので枯れてきているが、思い当たるものをあげてみよう。

Dropbox Paper

www.dropbox.com

共同作業のプラットフォーム。文書作成が主たる機能。GoogleのG Suiteが最も使われている分野だが、個人的にはDropbox Paperの方が優れていると思う。

Markdown Preview Enhanced

shd101wyy.github.io

AtomエディタとVSCodeのエクステンション。作者はイリノイ大学のYiyi Wangさん。 AtomエディタやVSCodeを使っていないと、使えないのだが、このエクステンションのためにAtomの導入をしてもよいくらいだと思う。このエクステンションを使うようになってから、Wordで文章を書くのが、ひどく煩雑で、バカらしくなる感じている。

Wangさんによれば、RStudioを参考に、このプラグインを作ったのだとか。Rのプログラムが走らないことを除けば、このアプリの方が優れたものに仕上がっている。
RStudioは文章が打ちづらく、エディタに作られている他のエクステンションも使えないので、編集はAtomの方が圧倒的に使いやすい。
AtomからRが操作できれば最強なのではないかと思う。

Thesaurus.com

https://www.thesaurus.com/

オンラインの英語の類語辞典。割とよく利用している。英文を書く上でシソーラスは必須。

Pronounce Names

www.pronouncenames.com

人物の読み方を調べるサイト。英語圏ではない人の名前が論文の著者にあるときが多いので、割と使っている。同種のサイトもあるので、併用すると良いと思う。Google検索で名前 Pronounce名前 How to PronounceとかでもOK。

Tables Generator

www.tablesgenerator.com

Markdown利用者は限定されるかもしれないが、Markdownのエディタは表づくりがあまり得意ではない。Excelのセルをコピーして貼って使えるところがよい。非常に便利。LaTex、HTML、Wiki形式にも対応している。

R-bloggers

www.r-bloggers.com

統計学による分析にRを使っている人限定だが、かなりレベルが高いサイト。よく知らないが、プラグラミングのQiitaみたいな存在なのかもしれない。

guru

www.guru.com

アウトソーシングのマッチングサイト。人力なのでもちろん有料。

勉強というより、仕事に使うサイトだろう。一般的な英語の校正であったり、アカデミック・ライティングの校正を請け負っている人がいるので、英語のアウトプットを出すときには利用するといいかもしれない。日本の業者の校正があまりよいと思ったことがないので、それよりは、guruの中でスキルの高い人を探してみる方がよいと思う。

以前紹介したgrammarly(https://app.grammarly.com/)にも、アウトソーシングで英文校正をしてもらう機能がある(有料)。grammarlyでの校正は非常に手軽という利点はあるが、人が選べないのと、クオリティ的に多少、難ありと感じることもあるので、guruも利用した方が良いと思う。日本語論文のアブストレベルであれば、grammarlyでも十分。

僕は英語しか使わないので英語に関してしか知らないが、guruであれば、その他言語への翻訳や校正をしてくれる人もいるように思う。

Rのmclustパッケージで潜在プロファイル分析

今回はRのmclustパッケージで潜在プロファイル分析を行う方法について解説したい。 潜在クラス分析はカテゴリカル変数、潜在プロファイル分析は(基本的には)連続変数という違いはあるものの、ケースごとのパターンからクラス分けをするのは同じである。

使い分けは次のようにするとよい。
潜在クラス分析は2値データの解釈が行いやすいため、基本的には2値データでの使用に留めるのが賢明である。
順序の無い3カテゴリの場合も潜在クラス分析をすることになる。しかし、解釈がややこしくなる。二項ロジスティック回帰分析の解釈が簡単なのに比べて多項ロジスティック回帰分析の解釈がややこしいのと同じである。とはいえ、解釈ができないわけではないので頑張ってクリアカットな結果を出すように頑張るれば使えないわけではない。

もし、順序のついているもの、つまり順序変数であれば、潜在プロファイル分析をした方がよい。
リコードをして潜在クラス分析をしている論文もあるが、潜在プロファイルを使うべきである。 例えば、リッカート尺度で計測したようなものである。「あてはまらない」-「どちらでもない」-「あてはまる」といった、5件法、7件法で計測した質問項目である。

正確に言えば、順序変数と連続変数の水準は異なるのだが、ある程度の順序があれば、連続変数とみなして計算しても、それほど問題は起こらないことがわかっている。シミュレーションによって諸説あるが、だいたい8カテゴリある順序変数であれば、連続変数と同等の結果が得られる。5件法だとさすがに連続変数と同等とはいえないのだが、潜在クラス分析をするためにリコードを行うの方が情報が失われるため、潜在プロファイル分析を行う方が、よほどまし、ということである。

今回取り上げるのはmclustパッケージの開発者たちが書いた論文の例題である。

www.ncbi.nlm.nih.gov

この例題がこちらのブログで解題されていたので、このブログ記事を参照している。

www.r-bloggers.com

データ

データは、Kaggle.comで無料で使用できるYoung People Surveyである。

library(tidyverse)
survey <- read_csv("https://raw.githubusercontent.com/whipson/tidytuesday/master/young_people.csv") %>% select(History:Pets)

Young People Surveyはだいたい以下のような感じである。

  History Psychology Politics Mathematics Physics Internet    PC `Economy Manage… Biology Chemistry Reading Geography `Foreign langua…
    <dbl>      <dbl>    <dbl>       <dbl>   <dbl>    <dbl> <dbl>            <dbl>   <dbl>     <dbl>   <dbl>     <dbl>            <dbl>
1       1          5        1           3       3        5     3                5       3         3       3         3                5
2       1          3        4           5       2        4     4                5       1         1       4         4                5
3       1          2        1           5       2        4     2                4       1         1       5         2                5
4       4          4        5           4       1        3     1                2       3         3       5         4                4
5       3          2        3           2       2        2     2                2       3         3       5         2                3
6       5          3        4           2       3        4     4                1       4         4       3         3                4
# … with 19 more variables

32の興味/趣味に関する項目が含まれていて、1(関心がない)から5(非常に関心がある)の5段階の順序変数で計測されている。

外れ値を除外する

KaggleによるとSatisficerがいるようである(Satisficerに関してはこちらを参照)。Satisficerとは同じ値を何度も選択した参加者など不適切な回答をしている者である。 carelessパッケージを使用して「文字列応答」(string responding)を識別する。マハラノビス距離で多変量の外れ値を探す。

library(careless)
library(psych)
interests <- survey %>%
  mutate(string = longstring(.)) %>%
  mutate(md = outlier(., plot = FALSE))

最大10に応答する文字列に上限を付け、α=.001のマハラノビス距離をカットオフとして外れ値を除外する。

cutoff <- (qchisq(p = 1 - .001, df = ncol(interests)))

interests_clean <- interests %>%
  filter(string <= 10,
         md < cutoff) %>%
  select(-string, -md)

外れ値を除外した分析データはinterests_cleanの中に格納された状態になっている。

潜在プロファイル分析

潜在プロファイル分析はmclustパッケージを使用して実行する。

library(mclust)
interests_clustering <- interests_clean %>%
  na.omit() %>%
  mutate_all(list(scale))
BIC <- mclustBIC(interests_clustering)

BICをプロットする。

plot(BIC)

f:id:iDES:20191210162844p:plain

この分析の場合、プロットだけでは、ラインが近接しているので見にくいので、上位3つのBICsummary()を使って表示させる。

summary(BIC)

下記のように表示される。

Best BIC values:
            VVE,3       VEE,3      EVE,3
BIC      -75042.7 -75165.1484 -75179.165
BIC diff      0.0   -122.4442   -136.461

最も良好な値を取っているのはVVE,3である。とはいえ、この3つの統計量にそれほどの違いはない。
各モデルを視覚的に表現すると下記のようになる。こちらの論文に掲載されている図である。

f:id:iDES:20191210162932p:plain

VEE,3は分布の形状が等しくなるように制約をかけているので、実際に分析するのはこのモデルがよいと考えられる。
VEE,3は以下のようなモデルである。

mod1 <- Mclust(interests_clustering, modelNames = "VEE", G = 3, x = BIC)
summary(mod1)

結果は以下。

---------------------------------------------------- 
Gaussian finite mixture model fitted by EM algorithm 
---------------------------------------------------- 

Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: 

 log-likelihood   n  df       BIC       ICL
      -35455.83 874 628 -75165.15 -75216.14

Clustering table:
  1   2   3 
137 527 210 

これらの統計量単体で適切性は判断できないが、クラスタリングの結果、極端に少ないグループがないことが確認できる。

ICL: Integrated Completed Likelihood (ICL)

BICのほかにICLでもモデルのフィッティングを見る。

ICL <- mclustICL(interests_clustering)
plot(ICL)
summary(ICL)

結果をプロットすると以下のようになる。

f:id:iDES:20191210163018p:plain

Best ICL values:
             VVE,3        VEE,3      EVE,3
ICL      -75134.69 -75216.13551 -75272.891
ICL diff      0.00    -81.44795   -138.203

グラフからも、ICLのサマリからも、BICと同じ傾向が確認できる。 ICLで最もフィッティングがよいのはVEE,3である。

以上の分析からVEE,3が最も適切であると判断できる。根拠となるのは以下の3点である。

  1. 分布の形状が同じになる制約
  2. BICでVVE,3とほとんど変わらず2番目にフィッティングが良好なモデル
  3. ICLで最もフィッティングが良好なモデル

BLRT

潜在プロファイル分析でもBLRTを行う。 モデルの選択ではなく、クラス数の選択であり、BICやICLで確認していることとは異なるので注意が必要である。

mclustBootstrapLRT(interests_clustering, modelName = "VEE")

結果は以下。

------------------------------------------------------------- 
Bootstrap sequential LRT for the number of mixture components 
------------------------------------------------------------- 
Model        = VEE 
Replications = 999 
              LRTS bootstrap p-value
1 vs 2    197.0384             0.001
2 vs 3    684.8743             0.001
3 vs 4   -124.1935             1.000

3 vs 4でP-Valueが有意ではなくなっており、4クラスモデルで改善が見られなかったということなので、3クラスモデルが良いという提案である。
BIC、ICLはどちらも3クラスモデルを提案していたので、3クラスモデルで間違いなさそうである。

ちなみにBLRTが良いと言っているのは、Mplus開発グループの論文Nylund,Asparouhov, Muthén(2007)である。
https://www.statmodel.com/download/LCA_tech11_nylund_v83.pdf

プロット

各プロファイルの平均値を抽出し、+1 SDを超える値をトリミングする。この作業をしてないとプロットで枠外へはみ出るので必要になる作業である。

library(reshape2)

means <- data.frame(mod1$parameters$mean, stringsAsFactors = FALSE) %>%
  rownames_to_column() %>%
  rename(Interest = rowname) %>%
  melt(id.vars = "Interest", variable.name = "Profile", value.name = "Mean") %>%
  mutate(Mean = round(Mean, 2),
         Mean = ifelse(Mean > 1, 1, Mean))

ggplot2でプロットを行う。

library(ggplot2)
p1<- means %>%
  ggplot(aes(Interest, Mean, group = Profile, color = Profile)) +
  geom_point(size = 2.25) +
  geom_line(size = 1.25) +
  scale_x_discrete(limits = c("Active sport", "Adrenaline sports", "Passive sport",
                              "Countryside, outdoors", "Gardening", "Cars",
                              "Art exhibitions", "Dancing", "Musical instruments", "Theatre", "Writing", "Reading",
                              "Geography", "History", "Law", "Politics", "Psychology", "Religion", "Foreign languages",
                              "Biology", "Chemistry", "Mathematics", "Medicine", "Physics", "Science and technology",
                              "Internet", "PC",
                              "Celebrities", "Economy Management", "Fun with friends", "Shopping", "Pets")) +
  labs(x = NULL, y = "Standardized mean interest") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")
p1

f:id:iDES:20191210163120p:plain

赤のグループは科学に興味があり、青のグループは芸術と人文科学に大きな関心があり、緑のグループは、科学と芸術の両方に関心がないように見えるが、他のことにある程度関心があるように捉えられる。

そこで、1クラスを科学(Science)、2クラスを無関心(Disinterest)、3クラスを人文(Arts & Humanities)とする。
それぞれのクラス名を代入したプロットに書き換える。上のコードとの違いはmutateでクラス名を入れているだけである。

p2 <- means %>%
  mutate(Profile = recode(Profile, 
                          X1 = "Science: 16%",
                          X2 = "Disinterest: 60%",
                          X3 = "Arts & Humanities: 24%")) %>%
  ggplot(aes(Interest, Mean, group = Profile, color = Profile)) +
  geom_point(size = 2.25) +
  geom_line(size = 1.25) +
  scale_x_discrete(limits = c("Active sport", "Adrenaline sports", "Passive sport",
                              "Countryside, outdoors", "Gardening", "Cars",
                              "Art exhibitions", "Dancing", "Musical instruments", "Theatre", "Writing", "Reading",
                              "Geography", "History", "Law", "Politics", "Psychology", "Religion", "Foreign languages",
                              "Biology", "Chemistry", "Mathematics", "Medicine", "Physics", "Science and technology",
                              "Internet", "PC",
                              "Celebrities", "Economy Management", "Fun with friends", "Shopping", "Pets")) +
  labs(x = NULL, y = "Standardized mean interest") +
  theme_bw(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "top")
p2

f:id:iDES:20191210163135p:plain

図で確認すると、各クラスの特徴がわかりやすくなる。
いろいろと傾向が読み取れるが、面白いのは、科学(Science)のクラスに属している人は映画館/演劇とセレブに興味がないというのは世情に疎いということだろうか。思い当たる節は確かにある。バリバリの理系の人と話をするときには、芸能話題を選ぼうとは思わない。日本であれば何かのオタクである確率が高い気がする。人文(Arts & Humanities)クラスが経済管理(Economy Management)に興味がない、と出ているが、確かに人文系と経済は相性が悪い。人文系には、経済学的な発想がないというのもあるし、個人的な金銭管理にもあまり興味がないように思う。

余談だが、mclustのブートストラップの反復回数は999回のようである。マシンパワーにもよるが分析は多少時間がかかる。
以前、99について書いたが最近は999ルールなのかもしれない。

RからMplusを使う Rですべて完結させる

RからMplusが使えるようになるMplusAutomationパッケージについてのエントリ2つ目。
前回はデータをMplus形式に変換する方法について解説した。

ides.hatenablog.com

今回も下記の条件で書いている。

  1. RStudioを使っている
  2. ファルダパスに日本語が混じっていない

今回はR上からMplusの分析を行う方法について解説する。

分析データと分析コード

パッケージのガイドに掲載されている簡単な例を実行する。 使用するデータはmtcarsである。Rにデフォルトで入っているので比較的有名なデータである。

mtcars function | R Documentation

head(mtcars)

以下のようなデータである。

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

今回使用するのは4つの変数。

  • mpg...マイル/ガロン(燃費)
  • hp...Gross horsepower...総馬力
  • wt...重さ
  • disp...排気量

Mplusの分析コードは下記のように書く。

TITLE: Path Mode;
DATA: FILE = "Path_model.dat";
MODEL: 
     mpg ON hp;
     wt ON disp;
OUTPUT:"CINTERVAL;

これをR上から実行する。

Mplusの分析コマンドをR上から実行

library(MplusAutomation)
pathmodel <- mplusObject(
   TITLE = "Path Mode;",
   MODEL = "
     mpg ON hp;
     wt ON disp;",
   OUTPUT = "CINTERVAL;",
   rdata = mtcars)
fit <- mplusModeler(pathmodel, modelout = "Path_model.inp", run = 1L, hashfilename = FALSE)

このスクリプトを走らせると、Rstudioの作業フォルダ内にデータの.datファイル、分析コードファイル.inp、分析結果のファイル.outが出力される。

変更するのは下記の点である。

  • コロンはイコールにする
  • ダブルクォーテーションで囲う
  • 文末にコンマ

hashfilename= FALSEにすると、inpファイル名と同じデータファイル(.dat)が作成される。

結果をRに読み込む

Mplusの結果はpath_model.outに出力されている。
結果を見たい場合は、直接見てもよいだろう。

Rへ結果を返したいときには、下記の関数を使う。

res<-readModels()

結果をresに格納する。 走らせた分析が1つでout形式ファイルが1つの時にはこの方法が取れる。

結果は下記のように引数を指定すると表示できる。
まずは統計量。

res$summaries

下記のようにやや混沌とした形で出力される。

  Mplus.version       Title AnalysisType   DataType Estimator Observations NGroups NDependentVars NIndependentVars NContinuousLatentVars
1             8  Path Mode;      GENERAL INDIVIDUAL        ML           32       1              2                2                     0
  Parameters ChiSqM_Value ChiSqM_DF ChiSqM_PValue ChiSqBaseline_Value ChiSqBaseline_DF ChiSqBaseline_PValue       LL UnrestrictedLL   CFI
1          7       15.564         2         4e-04             106.604                5                    0 -101.059        -93.277 0.867
    TLI     AIC     BIC    aBIC RMSEA_Estimate RMSEA_90CI_LB RMSEA_90CI_UB RMSEA_pLT05  SRMR     AICC       Filename
1 0.666 216.117 226.377 204.555           0.46         0.266         0.685       0.001 0.136 220.7837 path_model.out

BICだけ取り出したいときは更にBICをつける。

res$summaries$BIC

こちらはシンプルでよい。

226.377

推定結果は次のように書くと返される。

res$parameters$unstandardized

出力。

         paramHeader param    est    se est_se  pval
1             MPG.ON    HP -0.065 0.009 -6.889 0.000
2              WT.ON  DISP  0.006 0.001  8.739 0.000
3            WT.WITH   MPG -1.020 0.384 -2.657 0.008
4         Intercepts   MPG 29.586 1.529 19.346 0.000
5         Intercepts    WT  1.816 0.180 10.112 0.000
6 Residual.Variances   MPG 14.045 3.524  3.986 0.000
7 Residual.Variances    WT  0.209 0.056  3.751 0.000

95%信頼区間ci.を挟むんで$parameters$ci.unstandardizedと書くと返される。

texregパッケージできれいに表示

きれいに表示するにはtexregパッケージが使用できるようだ。

library(texreg)
screenreg(fit, summaries = c("Observations", "CFI", "SRMR"), single.row=TRUE)

下記のように整えて出力される。

==================================
                  Model 1         
----------------------------------
 MPG<-HP          -0.06 (0.01) ***
 WT<-DISP          0.01 (0.00) ***
 WT<->MPG         -1.02 (0.38) ** 
 MPG<-Intercepts  29.59 (1.53) ***
 WT<-Intercepts    1.82 (0.18) ***
 MPG<->MPG        14.04 (3.52) ***
 WT<->WT           0.21 (0.06) ***
----------------------------------
Observations      32              
CFI                0.87           
SRMR               0.14           
==================================
*** p < 0.001, ** p < 0.01, * p < 0.05

モデルを比較する

Rのlmやglmで使用するupdate関数が利用できる。 下記のように新しいパスを追加する。

pathmodel2 <- update(pathmodel, MODEL = ~ . + "
    mpg ON disp;
    wt ON hp;")
fit2 <- mplusModeler(pathmodel2, modelout = "model2.inp", run = 1L, hashfilename = FALSE)

この結果を先ほどのモデルと比較する。

screenreg(list(
  extract(fit, summaries = c("Observations", "CFI", "SRMR")),
  extract(fit2, summaries = c("Observations", "CFI", "SRMR"))),
  single.row=TRUE)

出力。

====================================================
                  Model 1           Model 2         
----------------------------------------------------
 MPG<-HP          -0.06 (0.01) ***  -0.02 (0.01)    
 WT<-DISP          0.01 (0.00) ***   0.01 (0.00) ***
 WT<->MPG         -1.02 (0.38) **   -0.73 (0.26) ** 
 MPG<-Intercepts  29.59 (1.53) ***  30.74 (1.27) ***
 WT<-Intercepts    1.82 (0.18) ***   1.68 (0.19) ***
 MPG<->MPG        14.04 (3.52) ***   8.86 (2.21) ***
 WT<->WT           0.21 (0.06) ***   0.19 (0.05) ***
 MPG<-DISP                          -0.03 (0.01) ***
 WT<-HP                             -0.00 (0.00)    
----------------------------------------------------
Observations      32                32              
CFI                0.87              1.00           
SRMR               0.14              0.00           
====================================================
*** p < 0.001, ** p < 0.01, * p < 0.05

95%信頼区間を併記する場合にはcis=TRUEをオプションに追加しておく。

screenreg(list(
  extract(fit, cis=TRUE, summaries = c("Observations", "CFI", "SRMR")),
  extract(fit2, cis=TRUE, summaries = c("Observations", "CFI", "SRMR"))),
  single.row=TRUE)

結果省略。

パッケージの限界

上記のように書くと非常に便利な機能のように感じられるが、実際の所、すべての分析に対応しているわけではなさそうである。
因子分析を走らせてみたところ、結果のMplusが分析するとこまではできるが、Rへの結果のリターンができなかった。 Mplusに分析させるまでは、おそらくシンプルではないモデルであっても、可能だがRへ結果を返せないのだ。

Rは図表を整えるパッケージが揃っているので、RにMplusがデータを渡すことができれば、かなり重宝するだろう。
僕が理解できる範囲では一部の分析結果にとどまっている。

MplusAutomationパッケージでMplusコマンドを書くとややこしくなる。
それであれば前回の変換編で書いたように、データだけ書いてMplusで実行ファイルを書いた方がよいと思う。

RでMplusの分析コードを書いてMplusに計算させるというのは、実用性はそれほど高くなく、やはり、結果をRで読み込めることの方が重要性は高い。
結果の読み込みができるとパーフェクトに近いパッケージになるように思われた。

参考までに こちらで行った因子分析をMplusAutomationパッケージの形式で書き換えたものを掲載しておこう。

library(psych)
data(bfi)
FAmodel <- mplusObject(
  TITLE = "Exploratory Factor Analysis by Mplus via R, using NEO data.;",
  VARIABLE = "USEVARIABLES = A1-A5 C1-C5 E1-E5 N1-N5;",
  ANALYSIS = "TYPE = EFA 4 6;
              ESTIMATOR = ML;
              ROTATION = oblimin;
              PARALLEL = 50;",
  PLOT = "TYPE = plot2;",
  OUTPUT = "residual;",
  usevariables = c("A1", "A2", "A3", "A4", "A5", "C1", "C2", "C3", "C4", "C5",
                   "E1", "E2", "E3", "E4", "E5", "N1", "N2", "N3", "N4", "N5", 
                   "O1", "O2" ,"O3", "O4", "O5"),
  rdata = bfi)
FA_fit <- mplusModeler(FAmodel, modelout = "bfi.inp", run = 1L, hashfilename = FALSE)

RからMplusを使う ファイル変換編

RからMplusが使えるようになるMplusAutomationパッケージというものがある。

cran.r-project.org

有用な関数が複数含まれているので使い方がいくつかあるが、今回はRを使ってMplusのデータとコードの一部を作る方法を紹介しよう。

RStudioを利用する

RとMplusのユーザーであれば、RStudioを使っていると思うので、RStudioが導入されている前提で話を進める。

RStudioの使い方はこちらを参照のこと。

www.idesohei.net

Mplusの分析ファイルはどこに置けばいいの?

日本語で得られる情報では、Mplusの作業ディレクトリは絶対パスで設定しているものばかりだがRStudioと同じく、パスの指定は不要である。 ただ、2バイト文字をパスに含めないという違いがあるだけである。 詳しくはリンク先のMplusの作業ディレクトリ解説を参照のこと。

ides.hatenablog.com

Mplusの形式でのデータ書き出し

psychパッケージにあるbfiというデータをMplus形式で書きだす。

library(MplusAutomation)
library(psych)
data(bfi)

bfiはこのブログではよく使っている5因子パーソナリティのデータIPIP-NEOである。 データの詳細はこちらを参照のこと。 http://ides.hatenablog.com/entry/2019/03/19/093726

変数名を抽出する

変数すべてを書き出すことはおそらく稀で、分析に必要なデータだけを書き出すことが多い。 そこで、書き出す変数の指定をするために変数名を抽出しておいた方が便利だ。

variable.names(bfi)

以下のように出力される。

 [1] "A1"        "A2"        "A3"        "A4"        "A5"        "C1"        "C2"        "C3"        "C4"        "C5"       
[11] "E1"        "E2"        "E3"        "E4"        "E5"        "N1"        "N2"        "N3"        "N4"        "N5"       
[21] "O1"        "O2"        "O3"        "O4"        "O5"        "gender"    "education" "age"     

Mplusデータ形式として出力

因子分析だけする場合には、"gender","education","age"の変数は使わないのでオミット(ドロップ)させる。保持する変数をしている場合はkeepColsと書く。

prepareMplusData(bfi, filename="bfi.dat", 
                dropCols=c("gender","education", "age"),
                overwrite=T)
  • bfi...出力するデータ
  • filename...出力するデータ名
  • keepCols...保持する変数
  • dropCols...保持しない変数...keepとdropはどちらかのみ指定

Mplusで使用するdatが同じディレクトリ(フォルダ)上に出力される。 Rの出力は下記のように出る。

DATA: FILE = "bfi.dat";
VARIABLE: 
NAMES = A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4 O5; 
MISSING=.;

これは、Mplusの分析コードであり、inp形式データの前半部分である。
つまり、 Mplusのコードの前半部分が自動生成されるのだ。

意外にMplusのコードでミスが出やすいのがこの箇所である。ANALYSISやMODELの部分はしっかり考えて書くので、割と間違えない。

データ出力だけ使用する場合には、inpファイルは通常の方法で作成する。

ちなみに、Mplusで因子分析をするには以下のような分析コードになる。

TITLE:""
DATA: FILE = "bfi.dat";
VARIABLE: 
    NAMES = A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4 O5; 
    MISSING=.;
ANALYSIS:
    TYPE = EFA 3 6;
    ESTIMATOR = ML;
    ROTATION = oblimin;
    PARALLEL = 50;
PLOT: TYPE = plot2;
OUTPUT: residual;

MplusAutomationの使い方としてはいくつかあるが、現実的にはここまでをRで行うのがいいのではないかと思う。
僕は今までExcelでMplus用のデータを作成していた。手間がかかるし、手間がかかればミスも出る可能性がある。なによりも、何をどのように処理したか、記録が取れないのが最も問題だった。個人的にはかなりの福音である。

1台のパソコンで異なったバージョンのMplusを使うには

今日は小技の紹介。
1台のパソコンで異なったバージョンのMplusを使うにはどうしたらよいか、という話である。

このようなことを試みた人はあまりいないかもしれないが、Mplusの異なったバージョンを同じPCにインストールしてもうまく使えない。 今回のはその解決法である。

Mplusの仕組み

Mplusのエディタは下記のものだ。

f:id:iDES:20191203005506p:plain:w400

Mplusの正式な名称はわからないが、いわゆるエディタの類なので、本エントリではエディタと呼ぶことにする。

実はMplusの本体はこれではない。エディタはフロントエンドにすぎない。
本体はこちらである。

f:id:iDES:20191203005841p:plain:w400

分析を実行すると、黒いウィンドウが出て計算をしているやつである。

Mplusはエディタで書いた分析の命令を本体に送って計算させて、終了するとエディタに返す、という仕組みである。

異なったバージョンのMplusを使うために必要な2つのこと

以下の2つを実施すると、異なったバージョンのMplusが利用可能である。

  1. MplusをそれぞれMplus8やMplus7といったフォルダにインストールする
  2. Pathを使うバージョンを変更するごとにMplus8やMplus7というフォルダに変更する

インストール場所をデフォルトから変更する(インストール時)

64bit版の場合、Mplusのデフォルトのインストール場所は下記の場所である。

C:\Program Files\Mplus

これを"Mplus8"というようにバージョンの番号をいれる。

C:\Program Files\Mplus8

Pathを変更する

つぎにPathを変更作業を行う。 Mplusをインストールするとpathが自動的に設定される。

f:id:iDES:20191203010520p:plain

この図ではMplusのPathはC:\Program Files\Mplusに通っている。
この場所がMplusの本体の格納場所を示している。
エディタから本体へ分析命令をどこに送ればいいか、ということを指定する。

したがって、このPathはMplus8に変えたり、Mplus7に変えれば、異なったMplusで分析ができる。

Pathの変え方は以下のページで解説がされている。

proengineer.internous.co.jp

どのバージョンで分析したか

分析結果のファイルである".out"の先頭行に書かれてある。
例えばこんな感じである。

Mplus VERSION 7
MUTHEN & MUTHEN
12/03/2019  08:40 AM

Mplus7で分析しているつもりでも、Mplus6に分析が送られていた、ということも起こりうるので、outファイルの先頭行には注意しておこう。