チラシの裏の落書き日記

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

研究について思うこと。

「研究とは何か」ということについて。

研究とはどんなゲームか。
分野にとって解く事ができる,価値のある重要な問題に効率的に解を与えるゲーム。

  1. 人を納得させればよい
  2. ゴールの見極めが大事。
  3. 報告はフォーマットに則って行う。

ルールをキチンと理解することがゲームを上手く進めるコツ。

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)

出力される図ひとつめ
f:id:kazzstat:20160620111433p:plain

出力される図ふたつめ
f:id:kazzstat:20160620112040p:plain

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ステップで潜在クラスの数を決める。

  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について。

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