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)

出力される図ひとつめ

出力される図ふたつめ