list操作めも
Rのリスト操作について改めてまとめる。
特に,purrrなどを用いた操作を理解すること。
purrrの使い方についてのページ
purrr: ループ処理やapply系関数の決定版 - Heavy Watal
Hadley神のListについてのページ
Lists · R for Data Science
forループ的なlapplyやらmapの使い方が便利。
forループを使う場合には,
df <- df <- data.frame( a = rnorm(10), b = rnorm(10), c = rnorm(10), d = rnorm(10) ) temp <- numeric(length(df)) for(i in 1:length(df)){ temp[i] <- mean(df[[i]]) } temp
などと事前に結果を入れるオブジェクトを容易しないといけない。
一方,lapplyやsapplyを使うと,わりと簡潔に書ける。
lapply(df, mean) sapply(df, mean)
purrrパッケージのmapなどはlapplyを今風に処理するためのもの。
map_dfなどは型を指定して,結果をベクトルで返すもの。結果はベクトルなので注意。
library("purrr") map(df,mean) map_dbl(df, mean) map_df(df, mean)
さらに,mapやlapplyをforループ的に使うこともできる。これは言われれば納得だが,知らないと確かに書けない。今後使って行きたい。
# for loop的な使い方 map(1:3, ~df[,.]) map(1:5, runif)
また,特定の条件を満たす要素に関数を適用する,map_ifとmap_atも地味に便利。
# 0より大きい要素に対して,なんか処理をする。 -2:2 %>% map_if(.,.>0, as.factor) # 部分的に適用して,残りのリスト全体を返す。 ? map_if iris %>% map_if(is.factor, as.character) %>% str() # 全部数値 mtcars %>% str() # 1,4,5列画文字に mtcars %>% map_at(c(1, 4, 5), as.character) %>% str()
purrrパッケージにはまだ良くわかってない使い方があるので,今後もまた調べていきたいところ。
Rの関数のソースコードなど見たいときの方法
一番有名な,関数名をタイプするものでは見れないものについて。
あまり日本語の記事が無いようなきがしたので,メモを残してておく。
library("Matrix") Matrix # Matrix関数の中にある謎の関数。このままだとエラーが返ってくる。 spV2M # コロン3つなのがポイント 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
3ステップで潜在クラスの数を決める。
- 最適な対数尤度探し
- optseedとtech 11の挿入
- 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)をやってくれる。詳しくは文献を参照してもらいたいが,潜在クラス数を決めるのには有効な方法である。