MplusのMCMCサンプルをGetする方法(+Mplusの結果をRでプロットする方法)その1

気が向いたらStanやJAGSなどとの対応関係をきちんと取りならが検算を行っていきたい。

MCMCサンプルは「.gh5」ファイルに入っているので，これを読み込む。

以下の資料のp.65あたりの「mplus.get.bayesian.」を利用すると，ギブスサンプラーのMCMCサンプルが得られるはず。
http://www.statmodel.com/mplus-R/Mplus%20R%20tutorial.pdf

以下のinpファイルを作成して，Mplusで推定を行う。
データは
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat
にあるものを，ダウンロードしてくる。

TITLE:	this is an example of an exploratory 
		factor analysis with continuous 
		factor indicators
DATA:		FILE IS ex4.1a.dat;
VARIABLE:	NAMES ARE y1-y12;
ANALYSIS:

TYPE = general;
ESTIMATOR = BAYES;
PROCESSORS = 2;
BITERATIONS = (2000);

model:

y1-y12 with y1-y12;
y1-y12;
[y1-y12];


OUTPUT: CINTERVAL;
plot: type = plot1 plot2 plot3;

library("MplusAutomation")
# MacではRを使ってプロットをする必要がある。
# 以下のsourceを読み込むと必要な関数がいろいろと読み込まれる。
source("https://www.statmodel.com/mplus-R/mplus.R")

# アウトプットの読み込み
hoge1 <- readModels("ex4.1bayes.out")
# 結果の表示
hoge1$parameters

# Rstudioのコンソールでmplusと入力して，タブを押すと色々な関数が候補として出てくる。

# プロットできるもの一覧を表示する。
mplus.view.plots("ex4.1bayes.gh5")
# 変数名を見る
(v_list <- mplus.list.variables("ex4.1bayes.gh5"))
# ヒストグラムの作成
mplus.plot.histogram("ex4.1bayes.gh5","Y1")
# 密度曲線を追加
mplus.plot.densityplot("ex4.1bayes.gh5","Y1")
# qqplotの追加
mplus.plot.qqnorm("ex4.1bayes.gh5","Y1")
# 散布図の作成
mplus.plot.scatterplot("ex4.1bayes.gh5","Y1","Y2")
# データの取得：1変数だけ。
mplus.get.data("ex4.1bayes.gh5","Y1")

ベイズ推定の結果をプロットするためのもの。
まだ利用しきれていない関数もあるので，今後追記する。

# 
# ベイズ推定の結果をいろいろと得るためのもの
# mplus.get.bayesian.●●
#
# パラメタ名を取得する
# ついでに，パラメタ名が必要な関数も表示してくれる。
(par_list <- mplus.list.bayesian.parameters("ex4.1bayes.gh5"))
#
# 特定のパラメタのMCMCサンプルを取得する。
# 列がチェインに対応している。
#
par_data <- mplus.get.bayesian.parameter.data("ex4.1bayes.gh5", par_list[1])
# トレースプロット
matplot(par_data, type="l")
mplus.plot.bayesian.traceplot("ex4.1bayes.gh5",paramstr = par_list[1])
# 自己相関
mplus.plot.bayesian.autocorrelation("ex4.1bayes.gh5",paramstr = par_list[1])
# 事後分布
mplus.plot.bayesian.distribution("ex4.1bayes.gh5",paramstr = par_list[1])

Mplusは便利なので，もう少し使いこなしていきたいところ。