チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

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は便利なので,もう少し使いこなしていきたいところ。

よく必要になるけど,よく忘れるR(とExcel)の小技メモ

文字(数字)の桁数を揃えたい。
参考URL:
http://d.hatena.ne.jp/teramonagi/20100902/1283429021

formatC("hoge", widh = 5, flag="0")

flagは '0 + - #' だけが指定可能とのこと。


不要な空白を消したい。

# いらない空白を全て削除する
# 全角半角の指定に注意。
gsub(" ", "", x, fixed = TRUE)
# 文字列の前後の空白を削除する
stringr::str_trim(x, side="both")

全角半角を変換したい。
Nipponパッケージを利用するといいらしいが,上手く動作しない場合がある。
そんなときは,Excelを使う。
参考URL:全角を半角に、半角を全角に変換する エクセルの関数技

文字を半角にするにはASC関数を使い,半角文字を全角にするにはJIS関数を使う。

不要な全角半角スペースを削除するには,=SUBSTITUTE(SUBSTITUTE(A1," ","")," ","")という関数を利用するとよい。
(trimは文字列前後の空白を削除するための関数。)


Rで図を出力する場合のあれこれ。
参考URL
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
R のグラフで日本語を使う - 知に至る病
https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html
Rのグラフで軸の目盛りの刻み幅を変更する方法 - Rプログラミングの小ネタ

# 文字のフォントを指定する。
# これを指定しておけば,とりあえず文字化けしない。
pdf(file="hoge.pdf",  family="Japan1GothicBBB")

# これはMac向き
par(family="HiraKakuProN-W3")  # ヒラギノ角ゴシックProN W3

# 空のプロットを作成する。
 plot.new()

# tclで軸を内側に向ける。
# 正の値が内側で,枠の大きさに対する割合で指定で指定する。0.3ぐらいがけっこうちょうどいい。
plot(0,1,tcl = 0.3)

# メモリを独自に振り直したいときに使う。
#メモリを消す
curve(dnorm, xlim=c(-3, 3), xaxt="n") 
# atに書き込みたい目盛りのベクトルを指定する
axis(side=1, at=-3:3) 
# sideは下から時計回りに1,2,3,4になっている。

# ラベルを各軸に並行して描く
las = 0
# ラベルをすべて水平に描く
las = 1
# ラベルを軸に対して垂直に描く
las = 2
# ラベルをすべて垂直に描く
las = 3

Mplusでターゲット行列を指定した回転を行う方法

Analysisコマンドで,rotation = targetを指定し,
modelコマンド内で目標とする値を~で指定する。
また,各因子の最後の(*t)などとつけて,各因子負荷が1つの因子負荷行列で表すことができることを示さないといけない。
Mplusのマニュアルは,versin 8のp681などに記載されている。正直,いろいろと回転がありすぎて何が何だかわからない。
本文をtargetで検索するのがはやい。
これを使えば,任意の値を指定したターゲット行列に向けた回転を行ってくれる。
(つまり,プロクラステス(procrustes)回転ができる。)
確認的因子分析よりもゆるく,探索的因子分析よりも仮説を反映した分析ができるので,もう少し利用されてもいい方法のように感じる。


使用したデータ
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat

TITLE:
DATA:
FILE IS ex4.1a.dat;
VARIABLE:
NAMES ARE y1-y12;
ANALYSIS: 
ESTIMATOR IS ML;
 ROTATION=TARGET;
 MODEL:
f1 BY y1-y12 y1~0 y2~0 y3~0 y4~0 y5~1 y6~1(*t);
f2 BY y1-y12 y7~0 y8~0 y9~0 y10~0 y11~1 y12~1(*t);
  OUTPUT:
MODINDICES CINTERVAL;
  

参考PDF
http://www.restore.ac.uk/appliedpsychometrics/Tutorial.materials/Introduction.to.Mplus/Introduction%20to%20Mplus%20session%202.pdf

Wordの小技(大文字小文字の変換,テキストボックス内の文字列の検索)

Wordで便利な小技を見つけた。
どちらも,「高度な検索と置換」を使う。
特に,ワイルドカードと,検索対象を指定するのがポイント。
「[0-9]」で検索を行うと,大文字の数字が選択される。その上で,「A->a」で「半角」を選択すればOK。

テキストボックス内の文字列の選択では,ワイルドカードの「*」を使って,検索対象を「メイン文書」ではなく「メイン文書のテキストボックス」にして検索をする。その後,検索ウィンドを閉じて「F9」をつかって,フィールドの更新をすればよい。

Wordで数字を大文字から小文字に変換する。
office-qa.com


Wordで複数のテキストボックスの文字列を選択する。
detail.chiebukuro.yahoo.co.jp


そして,フィールドについて,何処かでまとめて勉強したほうがいい気がする。
あとで,まとめて読んでみること。

pasonyu.com

pasonyu.com

フィールドについて

フィールドについて:Word(ワード)2010基本講座

平均値の信頼区間のシミュレーションプロット

Rで,信頼区間を書いていいプロットするコードのメモ。
layout関数で,上手くプット画面を分割して,さも1枚のグラフかのように見せかけるのがポイント。
けっこう可愛い感じのグラフを描画できる。

# 
# 縦向きの,信頼区間の図。
# ポイントは,レイアウトとparの使い方。
# 
layout(matrix(c(1,2),ncol=1), heights=c(1/5,4/5))
par(mar=c(0,3,1,1))
mu <- 0
sigma <- 1
n <- 5
size <- 100

curve(dnorm(x,mu,sigma/sqrt(n)),from=-3, to = 3,axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(mu, 0, mu, dnorm(mu,mu,sigma/sqrt(n)), lwd=2, angle = 0,lty =1, col =  rgb(1,0,1,alpha = 0.7))

x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,0,1))
plot(x_mean,1:size, xlim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(x_mean, 1:size, upp, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(x_mean, 1:size, low, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(v = mu , lwd = 2,lty=2,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。
f:id:kazzstat:20171020204451p:plain

#
# 横向きの,信頼区間の図。
#

layout(matrix(c(2,1),nrow=1), widths =c(4/5,1/5))
par(mar=c(3,0,1,1)) # 下,左,上,右の順で,余白を設定する。上が0なので,いい感じでくっつけられる。
mu <- 0
sigma <- 1
n <- 10
size <- 100

# 右に90度,図を転置させる。
t_dnorm <- function(x, mu, sigma){
  cbind(x, dnorm(x,mu,sigma)) %*% matrix(c(0,-1,1,0),ncol=2,byrow=T)
}
plot(t_dnorm(seq(from=-3, to =3, by=0.1), mu, sigma/sqrt(n)),type="l",axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(0, mu, dnorm(mu,mu,sigma/sqrt(n)), mu, lwd=2, angle = 0,lty =3, col =  rgb(1,0,1,alpha = 0.7))


x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,1,0)) # この余白の設定で,絶妙に図をくっつけるのがポイント。
plot(x_mean, ylim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(1:size,x_mean, 1:size, upp, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(1:size,x_mean, 1:size, low,lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(h = mu , lwd = 2,lty=3,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。
f:id:kazzstat:20171020204524p:plain

日々淡々と進むこと

最近思うこととしては、とにかく毎日進み続けることが何よりも大切。

ルーチンワークでやらなくてはならない雑務も多く、体力や気力も削られるかもしれないが、その中でも研究をし続け、論文をなんとしても書き続ける必要がある。

書き続けていれば、書くことの労力も減っていく。自分の中に書くということを続けることが大事。

 

スポーツをやるとこの辺りの感覚もわかる。

自分が練習しないと記録も伸びない。結局、自分との戦いだ。

Mplusによるモンテカルロ・シミュレーション覚書2

今回は前回紹介しなかったコマンドなどについて,紹介する。しかしあまり気合が入っていないのでグダグダな感じになっている。

SCALE OF DEPENDENT VARIABLES FOR ANALYSIS

CENSORED,CATEGORICAL,NOMINAL,COUNTオプションが分析の従属変数に利用できる。指定方法は前述したgenerateと類似していいて,「CENSORED = y1(a);」などと指定する。
この指定によって,モデルと推定での変数の扱いが決まる。
独立変数としては,2値あるいは連続値が指定できる。

OPTIONS FOR DATA ANALYSIS

CLASSES, AUXILIARY, SURVIVAL のみがanalysisコマンドで利用できる。

CLASSES

CLASSESはカテゴリカル潜在変数の名前とクラス数を指定する。
これを使うためにはTYPE =MIXTUREが必要。
Betweenレベルのカテゴリカル潜在変数を使う場合には,betweenオプションを使って,宣言をする必要がある。
使い方としては,「CLASSES = c1 (2) c2 (2) c3 (3);」などと指定する。

AUXILIARY

補助変数は分析モデル内では利用されない変数である。
TYPE=MIXUREとともに,AUXILIARYオプションは自動的に3ステップアプローチ(Asparouhov and Muthén, 2014a, b 参照)を実行するために使用される。
詳細はマニュアルp.880-881を参照。

SURVIVAL

このオプションはTIMECENSOREDオプションと同時に使われないといけない。
気が向いたら追記します。
p.881-883

VARIABLES WITH SPECIAL FUNCTIONS FOR DATA GENERATION AND ANALYSIS

WITHIN

「TYPE=TWOLEVEL:」の場合は,個人レベルの変数を「WITHIN = y1 y2 x1;」と指定する。ここで指定しない変数は,within, betweenの両方でモデル化される。
「TYPE=THREELEVEL:」の場合は,「WITHIN = y1-y3 (level3) y7-y9 (level2) y4-y6;」とする。レベル1,3, 2の順に指定することに注意。y1-y3は個人レベルで,レベル1でのみモデル化される。レベル1変数は他の変数に先立ってモデル化される必要がある。y4-y6はレベル1とレベル2でモデル化される。y7-y9はレベル1とレベル3でモデル化される。
あるいは
「WITHIN = y1-y3;
WITHIN = (level2) y4-y6;
WITHIN = (level3) y7-y9;」
などと指定する。この場合には順番は気にしない。
「TYPE=CROSSCLASSIFIED:」の場合には「WITHIN = y1-y3 (level2a) y4-y6 (level2b) y7-y9;」と指定する。

BETWEEN

クラスターレベルで測定された変数の指定。
withinオプションと同様の指定方法になる。

POPULATION, COVERAGE, AND STARTING VALUES

POPULATIONと COVERAGEは前回の記事参照。
また,「STARTING = estimates.dat;」によって初期値を指定することも出来る。

SAVING DATA AND RESULTS

「REPSAVE = 1 10-15 100;」で保存するデータを指定できる。
単一のデータを保存する場合には,「SAVE = rep1.dat;」とする。
「SAVE = rep*.dat;」とすると,レプリケーションの個別のデータを指定できる。

「RESULTS = results.sav;」で各推定結果を利用できる。
しかし,すこぶる使いにくい。outファイルに以下のような出力が得られ,パラメタの順がわかる。(PARAMETER SPECIFICATIONの部分でパラメタの順も得られるので,それを参照する必要がある。)

RESULTS SAVING INFORMATION

  Order of data

    Replication number
    Parameter estimates
     (saved in order shown in Technical 1 output)
    Standard errors
     (saved in order shown in Technical 1 output)
    Chi-square : Value
    Chi-square : Degrees of Freedom
    Chi-square : P-Value
    Wald Test : Value
    Wald Test : Degrees of Freedom
    Wald Test : P-Value
    H0 Loglikelihood
    H1 Loglikelihood
    Number of Free Parameters
    Akaike (AIC)
    Bayesian (BIC)
    Sample-Size Adjusted BIC
    RMSEA : Estimate
    SRMR

推定法にbayesを指定した場合に,全てのイテレーションのパラメタの事後分布のパラメタ値が得られる。
「BPARAMETERS = bayes.dat;」

時系列解析では,ラグを指定できる。
例:「LAGGED = y (1);」
最大のラグは2である。