Mplusメモ

Mplusはミクスチャーモデルの分析もらくらくできるし，計算も結構速い。しかし，人間業が深いもので，もっと早く分析して欲しい，なんていう欲望が出てくる。
そんな時はamalysisコマンドの中にprocessors = 4;などと記述すると，複数コアで計算を回してくれる。こんな簡単に並列計算が出来るなんていい世の中になったものだ。
このPDFにいろいろと書いてある。
https://www.statmodel.com/download/handouts/MuthenV7Part1.pdf

追記：PROCESSORS = a b; ! a = # processors, (b = # threads)らしい。
コア数とスレッド数の意味については以下のページが参考になる。
CPU選びの基本：マルチコアとマルチスレッド
Windows OS入門：第3回 プロセスとスレッド (1/2) - ＠IT

Windowsではsetコマンドを使うとか。
Windows機で搭載CPU数がわかるコマンド - その他([技術者向] コンピューター) 解決済 | 教えて！goo
Macでコア数などを調べる場合にはターミナルでsystem_profiler SPHardwareDataTypeとすればよい。

ちなみに，mixtureモデルを分析する場合には，outputでtech14を指定するとBootstrap Likelihood Ratio Test(BLRT)をやってくれる。詳しくは文献を参照してもらいたいが，潜在クラス数を決めるのには有効な方法である。

jagsUIなるパッケージ

MCMCサンプラーのうちの1つであるJAGSをRから利用する新しいパッケージができたようだ。基本はrjagsのラッパーのようだが，並列計算の指定などがやりやすそうだ。今後，これを中心に使っていくかもしれないのでヘルプにあった例をコピペしておく。初期値の指定が関数になっているが，いいのだろうか。詳細はのちほど検証のこと。

マニュアル
https://cran.r-project.org/web/packages/jagsUI/jagsUI.pdf
作者のgithub
GitHub - kenkellner/jagsUI: Run JAGS (Just Another Gibbs Sampler) analyses from within R.

#install.packages("jagsUI")
library("jagsUI")

#Analyze Longley economic data in JAGS
#Number employed as a function of GNP
######################################
##   1. Collect and Package Data    ##
######################################
#Load data (built into R)
data(longley)
head(longley)
#Separate data objects
gnp <- longley$GNP
employed <- longley$Employed
n <- length(employed)
#Input data objects must be numeric, and must be
#scalars, vectors, matrices, or arrays.
#Package together: several possible ways

#1. A named list of the objects
data <- list(gnp=gnp,employed=employed,n=n)
#2. A character vector of the names of the objects
data <- c( gnp , employed , n )
#3. A list of names of the objects
data <- list( gnp , employed , n )

######################################
##      2. Write model file         ##
######################################
#Write a model in the BUGS language
#Generate model file directly in R
#(could also read in existing model file)
writeLines("
           model{
           #Likelihood
           for (i in 1:n){
           employed[i] ~ dnorm(mu[i], tau)
           mu[i] <- alpha + beta*gnp[i]
           }
           #Priors
           alpha ~ dnorm(0, 0.00001)
           beta ~ dnorm(0, 0.00001)
           sigma ~ dunif(0,1000)
           tau <- pow(sigma,-2)
           }
           ", con="model.txt")
#Identify filepath of model file;
#in this case in the working directory
modfile <-  "model.txt"

######################################
##    3. Initialize Parameters      ##
######################################
#Best to generate initial values using function
inits <- function(){
  list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1),sigma=runif(1,0,3))
}
inits
#In many cases, JAGS can pick initial values automatically;
#you can leave argument inits=NULL to allow this.

######################################
##  4. Set parameters to monitor    ##
######################################
#Choose parameters you want to save output for
#Only parameters in this list will appear in output object
#(deviance is added automatically if DIC=TRUE)
#List must be specified as a character vector
params <- c('alpha' , 'beta' , 'sigma' )

######################################
##        5. Run Analysis           ##
######################################
#Call jags function; specify number of chains, number of adaptive iterations,
#the length of the burn-in period, total iterations, and the thin rate.
out <- jags(data = data,
            inits = inits,
            parameters.to.save = params,
            model.file = modfile,
            n.chains = 3,
            n.adapt = 100,
            n.iter = 1000,
            n.burnin = 500,
            n.thin = 2,
            parallel = T,
            n.cores=3)
#Arguments will be passed to JAGS; you will see progress bars
#and other information
#Examine output summary
out
#Look at output object elements
names(out)
#Plot traces and posterior densities
plot(out)
#Plot traces
traceplot(out)
#Update model another 1000 iterations
out <- update(out,n.iter = 1000)

追記
JAGSの並列計算を楽にしてくれる，runjagsパッケージもあるようだ。
https://cran.r-project.org/web/packages/runjags/runjags.pdf
このへんの使い方も理解しておきたい。
MCMCをするときの計算時間が尋常じゃないので，こういうものを使って少しでも楽になればいい。並列計算のやりかたを調べたり何だりしているときの虚しさたるや。。。それはもう悲惨なので。。。

学会に行く意味？

 学会に行く意味とはなんだろうか。今回はこのことに関して自分にとって何が良いことなのか考えてみる。

　まず，第1に締め切りがあることによる効果が挙げられる。締め切りがあることによって，適当なところで研究をまとめなければならない。コレによって，1つの区切りをつけることができる。1つの区切りがつけば，次にやることは1つの纏まった研究にするために何が必要なのかを考えることにシフトする。この作業なしに論文を書ければいいのだが，なかなかそうはいかないのが辛いところである。

　次に，学会発表によって人からのフィードバックを得ることが出来る点が挙げられる。内輪の議論ではなかなかできない外からの視点を得るためにも，学会発表は行うのが好ましいように思う。

　そして，単純な業績になるということである。学会発表自体はそこまで大きなものではないが，ないよりはマシである。もちろん学会発表だけではなく，キチンと論文化まで行わなければならないのだが。

　最後に，これが一番大事だと思うのだが，研究のモチベーションを保つというのがあるのではないかと思う。発表してフィードバックをもらうのが憧れの先生であるならば，これは非常に嬉しいことだ。さらに，そのような先生と懇親会で話ができるというのも楽しみなことである。先生以外にも，同じあるいは隣接した領域の知り合いが増えて，交流関係が広がるのも魅力だ。こういった，人との直接的な関わりというのは，論文を読むだけでは得られない。学会発表に足を運んでこういった交流を行うことが大事なように思う。

　おまけの拾い物としては，面白い発表を聞いて，最新(?)の研究動向を知ることが出来るというのもある。

　デメリットとしては，時間とお金がかかることである。しかし，ある意味では大学院生にとっての就活の一環であるので，サボるほうが将来的な損失は大きい可能性がある。最も，意味のない学会に出る必要はないのだが。

趣味を始める

　いきなり新しい趣味を始めたくなった。昔から見るのが好きだったジャグリングをやろうと思う。一番基本的なボールというやつから手をつける。3つのボールを小気味良くポコポコと投げ続けられるようになりたいものだ。
　ネットでやり方を調べてみたり，youtubeで実際の動画を見て見ながら，見よう見まねでやるものの，連続して投げ続けるのは難しい。だが，新しいことを始めるのは楽しいものだ。体を使うのも心地よい。普段使わない脳みその部分を使っている気分なのも大変によい。
　道具の場所もとらないし，持ち運びも楽なので，時間を見つけてちょくちょく練習することにする。ただし，三日坊主になってはいけないので，気合の入れ過ぎには注意する。

文章を書く上で気をつけること

自分で文章を書くときにやりがちなことをまとめておく。文章を書くのは苦手なので，注意点をキチンと書いてみる。この文章の目的は，自分にとっての現在把握している気をつけるべき点を明確にし，書き残すことで自分の作文能力を向上させるための資料を作成することである。

以下の点をあつかう。

1. 目的を明確にする。
2. 論理は不可欠
3. ”具体的に”という言葉の具体的さを考える。
4. 長い文章は書かない。
5. ギャップに注意する。
6. はじめと終わりを合わせる。
7. 接続詞に注意する。
8. 人にチェックをお願いする。
9. 自分で読み上げる。

目的を明確にする。

　まず，何か文章を書くということには目的が無くてはならない。特に自分のやっていることを人に伝えるための文章を書く場合には，何が伝わったらいいのか（十分なのか）を正確に定義しなければならない。目的のない”つれづれなるままに”書かれた文章が許されるのは，日記やかなり狭いコミュニティでの交流に限られると考えられる。
　多くの人に自分自身のやっていること，アイデア，研究内容を伝える際には目的が必要であるし，目的があるから次に何をしなければならないかが自然に出てくる。数学のように定義から積み上げていく書き方も大事ではあるが，限られた紙面・時間で概要を伝えるためには目的を先に明示するのが望ましい。
　目的には，何か社会的，分野的な問題意識も関わっている。”そもそも意味ないんじゃない？”と思われるのを防ぐためには，解決すべき問題がある，ということを明示しなければならない。基本的には，理想状態と現実のギャップがある，そのギャップを埋める，というのがわかりやすいテンプレのように思う。例えば，手法Aを使えばこんなに世界が良くなるが，使うためにはいくつかの問題がある，そこで，その問題を解決する，といった形式だ。目的は，理想の状態，問題はギャップとするとわかりやすい。
　このような大きな目的を設定できるかが勝負である。隠れたニーズや欲求を察知したり，適切に嗅ぎ分けられる嗅覚が必要である。これは多分にセンスに依存している気もするが，訓練次第で勘が働くようになりそうである。そもそも，解決する必要がある問題や目的を正確に定めるのが何より大事になる。これがずれると後の解決策も意味が無い。正確な問題の設定のために，先輩，上司，教授，友人，同僚，の手を借りることを躊躇してはいけない。問題のオリジナリティというは大事だが，それはそれとして，すごい人というのがどのように世の中を見ているのかに触れるためにも，自分の狭い世界に閉じこもっていてはいけない。

論理は不可欠

　論理的に書くというのはこれまた難しい。数式の証明ならば，わかりやすいが，論理的な文章となると，それだけで多くの本が書かれているくらい大変なテーマである。あんまり難しいことは考えずに，私が学んだ小さい一歩を書き留めておく。
　やるべきことは，同じ概念には同じ用語をキチンと割り当てる，独自の用語を使う場合にはその定義を明確に記述する，ということである。言い換えをしてもいいのは日常用語のみ。テクニカルタームは具体的に説明することは行ってもいいが，言い換えを行ってはいけない。これは基本として忘れないようにする。
　この上で，Aが問題だ，なぜならば●●，解決には●●が必要，なぜ解決できるか，というのを書いていけば良い。型を学ばなければならない。変にオリジナリティを出して行けはいけない。型にあてはめて書く。型にはめられないことは書かない。それは書いても伝わらないことだから。目的に合わせて，1つずつ必要な要素を解決していくことが必要。目的のために問題を分割する。分割した問題をさらに分割する。この作業を十分に繰り返していく。そうするとやるべきことが見えてくるハズ。あとは，それを解決すれば，ドミノ式にはじめの問題の解決に向かうはず。だが，念のため，それが解決した場合に，実際にどうなるかを想像することが必要である。

”具体的に”という言葉の具体的さを考える。

　具体的にはどれくらい具体的なのか。個人的には，具体的の加減が分からず非常に困っていた。いや，未だに困っているのだが。数式であれば，一般的な式に具体的な値を入れたり，具体例というものは作れるのだが。なんとなくわかってきたのが，いわゆる，5W1Hというものがハッキリしていることが具体的と言われた時に満たすべき条件として重要なのかもしれない。相手が，自分の文章を読んである程度正確にやることがわかるくらいに５W1Hの内容を書くことを目指す。
　研究であれば，仮説が明確（独立変数と従属変数が何で，分析が何か），対象は誰か，いつやるか，どのような実験をやるか，どのような理論を使うか，（どこでやるか）などがあれば結構具体的になるような気がする。というか，論文で求められている項目をそのまま言えばいいのか。なるほど。

長い文章は書かない。

　文章というか，一文を長くしない。できるだけ短文にする。イメージとしては英語を書く時のように，主語，述語，目的語など書いて，修飾は少なめにする。論理的に繋がる部分は，2文に分けても良い。というか，私の場合はわけないと分けの分からない文章に成ってしまう。基本は，短めの文章を書くことを心がける。

ギャップに注意する。

　論理のギャップは常に気をつけないといけない。意外と実証されていないこと，未定義のことを暗に仮定していることもあるので注意が必要だ。今言及していることは，直前の内容を正確に受けているか，内容はこれまでに説明していないことが入っていないかなど注意する。

はじめと終わりを合わせる。

　論理をつらつらつなげていくと，間違ってはいないけど，目的からはずれているということがしばしばある。そのため，目的を書いたら，その解答・結論を先に書いてしまったほうがいいかもしれない。目的と結論から議論がずれないように気をつける。

接続詞に注意する。

　接続詞のリストは何度でも確認したほうがいい。
http://pothos.main.jp/setuzokusi.htm
こういうものを参考にする。机の前に貼って自分のものになるまで何度も参照する。

人にチェックをお願いする。

　恥を忍んで，人に自分の文章を晒してく。基本的に誤字脱字などはなくして見てもらうのがマナーだと思うが，私自身は誤字脱字が多いので見せるのに躊躇してしまう。誤字脱字のチェック方法はどうにか確立しなければならない。
　人に診てもらうことによって，自分では考えていなかった観点から意見が出る。自分では自明だと思っていたことが相手にとっては全く自明でないことも多々ある。意見を聞いたら，取り入れるかどうかは自分で判断すること。複数人に見せていると色々な人の意見が入り混じってむしろ文章が混乱することがある。しかし，見せないよりは見せたほうがよい（少なくとも今の私には）。また，どのように見てもらいたいのかを明確にして依頼するのも生産的だろう。論旨を見て欲しい，てにをはを直して欲しい，正確かどうかを見て欲しい，かっこいい言い回しはないだろうか見て欲しい，などなど。今の自分の文章に必要なことを考えて頼む。そこまで思い至らない場合でもとにかく頼む。直してもらう。直される中で，自分に欠けている部分がわかったり，よい言い回しを学んだり，いろいろとメリットが有る。また，依頼した人には懇ろにお礼を忘れないこと。相手が困っていたら相手を助けること。お互いに良い関係を育むこと。

自分で読み上げる。

　誤字脱字チェック法。しかし，私は自分の文章に自身がないため読み上げるのが恥ずかしい。恥ずかしいが，非常に効果のある方法のようだ。読めない部分，読みにくい部分は一目瞭然。自分の声にならないのだから。ということで，私自身の課題は読むのを恥ずかしがらずにできるかどうかである。これは今後克服しなければならない重要な課題である。

　ほとんど書いた後にふと思ったが，自動の読み上げ機能を使うのもいいのではないかと思った。私は，Macを使用しているのでデフォルトでついているテキスト読み上げ機能を付けて読み上げを行ってみる。自分で読むのはなんだか恥ずかしい人種なので，今後はこれを活用してみたい。
選択したテキストの音声化 - Word for Mac

このような機能が実用に堪えるのかを検証するために，文章を書いたら読み上げ機能を使ってみることとする。
実際に使ってみたが，かなり間違いを発見できることがわかった。今後，これを活用する。一文字でも入っていないと読み上げる際にかなり変に聞こえるので，一発で分る。

この辺の記事も参考にする。Wordの校正機能をMAXにする。
penya.jp

このサイトも便利そうである。
www.value-press.com

今後，このような自分のための注意点を充実させていく。

Mplusのメモ

Mplusは潜在変数モデリングに適した解析ソフトで，心理学，社会学，疫学などの分野で近年人気が出ているソフトウェアである。ホームページは次のものがある。
Muthén & Muthén, Mplus Home Page

このソフトは非常に便利で，文法も簡潔なのであるが，文法が特に覚えられないため，自分のための備忘録として徐々にまとめを書き始める。基本的にはユーザーガイドを真似して，その都度似た分析のものをコピペするという使い方が基本。これまでにやった分析の書き方を基本に。

関西学院大学の清水先生がわかりやすいまとめを書いてくださっているので，この辺も参考に。
[http://www.slideshare.net/simizu706/mplus-lecture-1:embed:cite]

ユーザーガイドから幾つか例を書いておく。まず，MIMICモデル。これだけでもMplusの文法の重要な部分がかなり盛り込まれている。

DATA: FILE IS mimic.dat;
VARIABLE: NAMES ARE y1-y6 x1-x3;
MODEL:
f1 BY y1-y3;
f2 BY y4-y6;
f1 f2 ON x1-x3;

利用するデータファイルの名前をDATAコマンドに記述する。一文の終わりはセミコロン";"で締める。変数の名前はvariableコマンドで指定する。"-"を使って，一気に定義することができるので楽ちん。尺度ごとにまとめて変数を定義するといいかもしれない。最後にmodelコマンドでモデルを記述する。byは潜在変数の定義，onは回帰を意味する。これだけわかればかなりいろいろなモデルが記述出来ることができる。普通の人がAmosでやることは概ねできる。ちなみに，コメントは"!"で記述する。

次に成長曲線モデル。iは切片，sは傾きパラメータを意味する。この書き方は省略形である。

DATA: FILE IS growth.dat;
VARIABLE: NAMES ARE y1-y4 x1 x2;
MODEL:
i s | y1@0 y2@1 y3@2 y4@3;
i s ON x1 x2;

個人的にはbyを使ったi,sの定義の方がわかりやすい。

DATA: FILE IS growth.dat;
VARIABLE: NAMES ARE y1-y4 x1 x2;
MODEL:
s by y1@0 y2@1 y3@2 y4@3;
i by y1@1 y2@1 y3@1 y4@1;
i s ON x1 x2;

独立変数ありの探索的な潜在クラスモデリング。

DATA: FILE IS lcax.dat;
VARIABLE: NAMES ARE u1-u4 x;
CLASSES = c (2);
CATEGORICAL = u1-u4;
ANALYSIS: TYPE = MIXTURE;
MODEL:
%OVERALL%
c ON x;
u4 ON x;

ポイントは潜在クラスの数をvariableコマンド内で定義し，顕在変数はカテゴリカルな変数であると定義しているてんである。また，analysisコマンド内で，type=mixtureを指定するのがポイント。
モデルの%overall%でクラス全体に関わる部分を記述する。クラスごとにモデリングを行いたい場合には，%c#1%の中に具体的なモデルを記述する。詳細はまた後日。

最後は，マルチレベルモデリング。また難しい記述の仕方である。

DATA: FILE IS reg.dat;
VARIABLE:NAMES ARE clus y x w;
CLUSTER = clus;
WITHIN = x;
BETWEEN = w;
MISSING = .;
DEFINE:
CENTER x (GRANDMEAN);
ANALYSIS: TYPE = TWOLEVEL RANDOM;
MODEL:
%WITHIN%
s | y ON x;
%BETWEEN%
y s ON w;

新しいところは，variableコマンドにcluster, within, betweenが現れ，defineというコマンドが出てきたところである。さらに，analysisのtypeも何やら新しい。順番に見てみよう。

clusterはデータが何によってネストされているかを定義する。例えば，学校の中に個人がネストされていたり，個人の中に測定時点がネストされている構造がある。このとき，withinでレベル1(学校の中の個人，あるいは個人の中の測定時点)の変数を定義している。さらに，betweenではレベル2(学校，個人)の変数を定義している。
defineコマンドによって，レベル1の変数を中心化しており，全体平均によって中心化するか，レベル2の単位で中心化する場合が多い。"ANALYSIS: TYPE = TWOLEVEL RANDOM;"はとりあえず，マルチレベルモデリングを行うときのおまじないと思っておけばよい。

modelコマンドの中では，%within%で個人レベルの回帰のモデルを作成している。”s |”によって，回帰係数を"s"としており，この"s"と従属変数"y"がレベル2の変数によって説明されるモデルを構築している。

ここまでで，相当色々な分析が出来ることがわかる。今後も気が向いた時に，モデルの記述の方法や，出力の読み方(Rでも読み込ませ方)などをメモ書きに残していく。monte carloコマンドやbayes推定の方法など，実用的なところまでたどり着ければよいかと思う。

lmerでの3レベルモデルの実行

3レベルの階層線形を行う時のlmerのメモ

学校のなかに，個人がネストされていて，個人の中に測定時点がネストされている場合。学校IDをschool_id, 個人IDがind_id，あとは独立変数をx, 従属変数をyとする。簡単のため，中心化は行われているとして，レベル2以上の変数は仮定しない。基本は level3/level2などのように上のレベルの変数からスラッシュで区切れば良いようだ。

result <- le4::lmer(y ~ 1 + x + (1 + x|school_id/ind_id), data = dat)

参考にした記事
rpsychologist.com

残差の分散共分散行列の指定の方法が色々と書かれているので，かなり使える記事。
マルチレベルで必要な分析はここでほとんどまかなえる。