Rの関数のソースコードなど見たいときの方法

一番有名な，関数名をタイプするものでは見れないものについて。
あまり日本語の記事が無いようなきがしたので，メモを残してておく。

library("Matrix")
Matrix
# Matrix関数の中にある謎の関数。このままだとエラーが返ってくる。
spV2M
# コロン３つなのがポイント
Matrix:::spV2M
# あるいは以下の関数を使うと中身がみれる。
getAnywhere("spV2M")

untar(download.packages(pkgs = "Matrix", 
                  destdir = ".",
                  type = "source")[,2])

とすると，作業ディレクトリにパッケージがダウンロードされてuntarされて中身を見ることができるようになる。
srcフォルダにC++のソースコードなどがあるのでみてみること。

Rの組み込みのプリミティブな関数についてはpryrパッケージのshow_c_souceなんかで見れるらしい。

このページが非常に勉強になる。
stackoverflow.com

MplusでMCMCした結果を使ったSEMでの情報仮説の評価

ちょっと気になったので，近年一部界隈で話題になっているベイズ情報仮説の評価をやってみる。
普通に計算すると大変なベイズファクターも，情報仮説の元では簡単に計算できる。

元ネタは以下の資料
http://www.tandfonline.com/doi/abs/10.1080/17405629.2012.732719
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjZ9-6K-PXOAhVGpJQKHf9-B2kQFgghMAA&url=https%3A%2F%2Frensvandeschoot.files.wordpress.com%2F2012%2F07%2F02-bayes-informative-hypothesis.ppt&usg=AFQjCNHLedgBUV_FFYDOK-6AQVfL-55v3A&sig2=0kAg9LL2V3VDKPBcIF3yfg&bvm=bv.131783435,d.dGo
http://www3.psy.senshu-u.ac.jp/~ken/JJPS2014.pdf
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=4&cad=rja&uact=8&ved=0ahUKEwjP6YGA-_XOAhUHsJQKHVFFASQQFggyMAM&url=http%3A%2F%2Fir.acc.senshu-u.ac.jp%2Findex.php%3Faction%3Dpages_view_main%26active_action%3Drepository_action_common_download%26item_id%3D9914%26item_no%3D1%26attribute_id%3D15%26file_no%3D1%26page_id%3D13%26block_id%3D52&usg=AFQjCNFFbxBCCZtGqypBtdGwYOs0BKsRmw&sig2=BKGhm-TClAlSrhKU4aKDSw&bvm=bv.131783435,d.dGo
などなど。
あまり，この手の話題に詳しくはないが，MplusでMCMCをする練習にもなるので，よいかと思う。

# 使用するパッケージ
#　install.packages("MplusAutomation")
library(MplusAutomation)
# Mplusの出力から，R上でトレースプロット等を描くためのもの。
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")

folder <- "/Users/hoge/Dropbox/mplus/infohypo"

# 
# Mplusのinpファイルの作成
# 
cat("DATA: FILE = ex3.1.dat;
! User guideの例を利用
VARIABLE: 
NAMES ARE y1-y6 x1-x4;
USEVARIABLES ARE y1 x1-x3;
MISSING=.;

! ここの指定が大事
analysis:
! 推定方法をベイズに
estimator = bayes;
! チェイン数と連動する
processors =7;
chains = 7;
! イテレーション数。
! 何も指定しないと，半分がバーンインで捨てられる。
fbiter = 20000;
thin = 2;
! 初期値をMLで決める
stvalues = ML;
! アルゴリズムメトロポリスヘイスティングズも指定可能
algorithm = gibbs(px3);
!algorithm = MH; 

model:
y1 on x1-x3*(beta1-beta3);

! 事前分布の指定。
! なくても勝手に決めてくれる。
model priors:
beta1-beta3 ~ N(0, 1);

! 
output:
tech1 tech8 stdy;

! トレースプロットなど
plot:
type = plot2;

savedata:
! MCMCサンプルを保存
bparameters are bprm.dat;
!results are params.dat;
!ESTIMATES ARE estimate.dat;
  ",file=paste0(folder,"lm.inp"))

# モデルの実行
runModels(paste0(folder))


setwd(folder)
getwd()

# MplusAutomationの利用
btest <- getSavedata_Bparams(outfile = "lm.out", discardBurnin = TRUE)

# 変数名が独特なので，colnamesで確認
colnames(btest[[1]] )

# 仮説の表記の仕方。
# beta1 >beta2 かつ beta2 > beta3
testBParamCompoundConstraint(btest,"(Parameter.2_Y1.ON.X1..equality.label.>Parameter.3_Y1.ON.X2..equality.label.)&
  (Parameter.3_Y1.ON.X2..equality.label.>Parameter.4_Y1.ON.X3..equality.label.)")

# 制約がmetなものの割合を見ればよい？
0.2419/(1/6)

?testBParamCompoundConstraint
# beta1 < beta2 かつ　beta2 < beta3
testBParamCompoundConstraint(btest,"(Parameter.2_Y1.ON.X1..equality.label.<Parameter.3_Y1.ON.X2..equality.label.)&
  (Parameter.3_Y1.ON.X2..equality.label.<Parameter.4_Y1.ON.X3..equality.label.)")
# f_i に当たるものが算出出来る
# c_i に当たるものは，3!のなかの組合せのうち，当該の仮説の個数を計算すればよい。
# ここでは1通り。仮説が，beta1 > beta2, beta3 なら，2通りある。
0.1015571 /(1/6)


# 仮説が正しい場合と正しくない場合の比
(0.1015571 /(1/6)) /((1-0.1015571) /(1-1/6))


0.2419/(1/6) /(0.1015571 /(1/6))

# その他関数。
extractModelParameters()
extractModelSummaries()
getSavedata_Fileinfo("lm.out")

# トレースプロットをかいたりできる。
model <- readModels("lm.out")
mplus.traceplot(model, rows = 2, cols = 2, parameters_only = TRUE)
library(coda)
plot(model$bparameters$valid_draw)

Lavaanによるパラメータ推定結果を見やすくグラフにする。（Gelman's secret weaponのパクリ）

尊敬するHoxo_mさんらが翻訳された”みんなのR”にあるGelman's secret weaponを真似して自分の武器を作成してみようかと思う。
d.hatena.ne.jp
https://www.amazon.co.jp/みんなのR-データ分析と統計解析の新しい教科書-Jared-P-Lander-ebook/dp/B012Z86Q6Q

やっていることは簡単。lavaanパッケージでのパラメータ推定結果を流用して，ggplotのgeom_pointrangeを使って点推定値と95%信頼区間を描画しているだけ。
おまけで，5%水準で有意かどうか，固定パラメータかで，色の塗り分けなんかをやってみている。
これで，推定値とそのブレ，有意か否かがさくっと把握できる。
分散・共分散を除きたい場合なんかは投入するデータセットを変更すればよろしい。
色の好みなどは適宜変更すること。

# 使用するパッケージ
library(ggplot2)
library("lavaan")
# help(package="lavaan")
# sem関数の例題を利用
## The industrialization and Political Democracy Example 
## Bollen (1989), page 332
model <- ' 
# latent variable definitions
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + a*y2 + b*y3 + c*y4
dem65 =~ y5 + a*y6 + b*y7 + c*y8

# regressions
dem60 ~ ind60
dem65 ~ ind60 + dem60

# residual correlations
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'
# モデルのあてはめと要約
fit <- sem(model, data=PoliticalDemocracy)
summary(fit, fit.measures=TRUE)
sem_params1 <-  parameterEstimates(fit,standardized = T)


my_weapon <- function(sem_params){
  # パラメータのラベルを作成
  sem_params$coef_name <- apply(sem_params[,1:3],1,function(x)paste0(x,collapse = ""))
  # 係数が有意か否かで色を帰るためのもの。
  sem_params$is_sig <- factor(sem_params$pvalue >= 0.05, levels=c("FALSE","TRUE","Fixed"))
  
  sem_params$is_sig[is.na(sem_params$is_sig)] <- factor("Fixed")
  # SEと信頼区間の描画のためのもの
limits <- with(sem_params, aes(ymax = est+se, ymin= est-se))
  limits2 <- with(sem_params, aes(ymax = ci.upper, ymin= ci.lower))
    
  # ggplotの描画
  p <- ggplot2::ggplot(sem_params, aes(y=est,x=coef_name,colour=is_sig))
  p <- p + ggplot2::geom_point() + ggplot2::geom_pointrange(limits2, alpha=0.5)
  p <- p +  ggplot2::geom_pointrange(limits, size=0.7, alpha=0.3)
  # 転置させるもの。変数名が長くなりがちなので，あったほうがいいかもしれない。
  ggplot2::last_plot() + ggplot2::coord_flip()
  
}

# 実際にplotを作成してみる。
my_weapon(sem_params= sem_params1)


#
# 他のモデルでもやってみる
# 
HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 '

fit <- lavaan(HS.model, data=HolzingerSwineford1939,
              auto.var=TRUE, auto.fix.first=TRUE,
              auto.cov.lv.x=TRUE)
summary(fit, fit.measures=TRUE)
sem_params2 <-  parameterEstimates(fit,standardized = T)

my_weapon(sem_params = sem_params2)

出力される図ひとつめ

出力される図ふたつめ

RからJuliaを使う

Rから今流行のスクリプト言語Juliaを使いたい。
普通にJuliaにパスを通してRのsystem関数を使えばいいのだけれど，パスを通さずにRで下のような関数を書いて作業をする。jlfileにはJuliコードを書いた拡張子.jlであるファイル名を入れる。jlfilepathは作業フォルダ的なもの。jlpathはバージョンが変わったら変わる。

runjl <- function(jlfile, jlfilepath ="hogehoge"){
  jlpath <- "/Applications/Julia-0.4.5.app/Contents/Resources/julia/bin/julia" 
  system(paste0(jlpath," ",jlfilepath,jlfile))
}
# こんな感じで使う。
runjl("example.jl")

わざわざ記事を書くほどのものでも無いけど，せっかくなので。

追記：普通に.bash_profileにpathを書いて，juliaへのパスを通した方が楽。

Mplusで潜在クラス分析を行うときのめも

このノートの荒っぽいまとめ。
https://www.statmodel.com/examples/webnotes/webnote14.pdf

３ステップで潜在クラスの数を決める。

1. 最適な対数尤度探し
2. optseedとtech 11の挿入
3. optseedとtech 14の挿入

ブートストラップ尤度比検定を行うには結構時間がかかるので，いきなりやるより，手順を踏んだほうが良いようだ。
まず，クラス数kとk-1での最適な対数尤度を探す。
この時，start = 600 120;などとすると，初期値を600個ランダムに生成して，初期の対数尤度がいけてるもののうち120個を推定を行い，その中でベストなものを最終的なアウトプットにするというもの。
普通はそこまででいいのだけれど，潜在クラス数を決める場合には，誤ったseedだと検定がうまくできない。そこで，このステップではイケてる乱数のSEEDを探すということになる。
推定が終わったら，outputの中身を"seed"などで検索して，対数尤度のランキング表を見つけること。そこにイケてるseedがある。

つぎに，このseedを使って次のステップに入る。具体的には，
starts =0;!これは乱数の初期値を使わないという意味
opted = 484687;!seedを決める
として， OUTPUT: TECH1 TECH8 TECH11;とtech 11を加える。tech 11ではLo-Mendell Rubin(LMR)検定が出来る。詳細は, Lo et al.(2001)にあるらしい。ブートストラップ尤度比検定(BLRT)より計算量が少ないので，先にやってあたりをつける感じ。

最後はBLRTの登場このseedを使って次のステップに入る。具体的には，先ほどとどうように。
starts =0;
opted = 484687;
として， OUTPUT: TECH1 TECH8 TECH14;とtech 14を加える。tech 11はもはや入れなくて良い。
これで推定しておしまい。
ところが警告が出ることがある。
そんな時は，LRTSTARTS = 0 0 100 20;などとして初期値の個数を増やす。前半2つはクラス数k-1について，後半2つはクラス数kについて。

これで，いけるはず。駄目だったらモデル自体を再考するところから始める。

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をするときの計算時間が尋常じゃないので，こういうものを使って少しでも楽になればいい。並列計算のやりかたを調べたり何だりしているときの虚しさたるや。。。それはもう悲惨なので。。。