2016-05-03

今日のできごとなど

雑記

今日やったこと

某ポスター作成
発表資料の作成
MCMCの観察
なす購入

1.2.については直前に再度確認を行うこと。

3. について，バーンインは結構早い。サンプリングがえらく時間がかかる。いらないパラメータはサンプリングしないでいたほうがいいか？
欲しいところだけとる感じ。
thiningを大きめにして，自己相関を小さくすること。
4.なすは適当なときに植える。

明日やること

シミュレーションを実行し直すこと
申請書No.1の内容を修正すること
申請書No.2の草稿を作成すること。
某発表のための資料の暫定版を作成すること
対応の続き。
某お仕事について。

1.すっかり忘れていた。
2. わりかし大事。イントロの部分を中心に直す。印刷してみて考えること。
3.全く手をつけていないが，No1の流れを参考にしてみる？
4. 暫定版を作成。必要なシミュを追加する？
5.引き続き。ただ，相談は必要。
6.誰がやるのだろうか。。。？

2016-05-03

Rの使いかたメモ（purrrの使いかたについて）

R 統計

自分が使うとしたらどう使うか？

一番やりそうだなぁと思ったのは最後の@sinhrksさんが紹介してくれているような，複数のモデルを作成して，まとめて分析するというやり方だと思う。
例えば，共分散構造分析を行うときに複数のモデルを立てて分析する，というような使いかたや，探索的因子分析を行うときに，因子数や回転方法を変えながら分析するという感じ。

確認的因子分析の例はこんな感じ。
このように書くメリットは，モデルの構築の部分の分析を分けられるから，記述がすっきりすることかと。

library("purr")
library("lavaan")

# lavaanの例を流用
model1<- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 '

model2<- ' f1  =~ x1 + x2 + x3 + x4 + x5
          f2 =~ x6 + x7 + x8 + x9 '

model3<- ' f1  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 '

model4<- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 
f1  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 
f1 ~~ 0* visual
f1 ~~ 0* textual
f1 ~~ 0* speed'

model5 <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 
f1 =~ visual +  textual +speed'

models <- list(model1,model2, model3, model4, model5)

# 推定できないモデルもあるが，とりあえずは気にしない。
fits <- purrr::map(models,~lavaan::cfa(., data=HolzingerSwineford1939))
fits %>% map(.,~summary(., fit.measures=TRUE))
fits %>% map(.,~lavaan::fitMeasures(.))

summary(fits[[1]])

sapply(fits,lavaan::fitMeasures)

psychパッケージのfa関数をつかって分析をおこなう。本来は因子数を予め決める必要があるが，ここでは，因子数を変える。さらに，回転方法も回転無し，直交回転，斜交回転の３種類に対して分析を行なってみる。それぞれの回転方法にたいして，１から7因子の場合を当てはめている。一応screeplotや平行分析からは５，６因子あたりが良さそう。推定方法は最尤推定を用いる。

library("purrr")
library("psych")
library("GPArotation")

# Big fiveの25項目
dat <- bfi[1:25]

# 因子数の確認。
scree(dat)
fa.parallel(dat)

# 関数をつくる関数を作成する。
myfa <-function(rotate){
    function(nf){
      fa(dat,nf,rotate=rotate,fm="ml")
    }
  }

rotate <- list("none", "varimax",  "promax")
nf <- c(1:7)
# 
results <- purrr::map(rotate,myfa) %>% map(.,~lapply(nf,function(x).(x))) %>% flatten()
# どんな出力があるのかを確認
names(results[[1]])

results %>% map(.,~.[c("RMSEA")])
results %>% map_df(.,~.[c("BIC","EBIC","TLI")])