R

因子分析＆確証的因子分析

lavaan,psych,gdataパッケージ　gawk graphviz

OSはlinux(zorinOS)

windowsではlavaan:::sem使用時に変数名が日本語だとエラーになりました。

（参考）
lavaanチュートリアル(関西大学 荒木先生)
心理データ解析演習(心理デザインデータ解析演習) 因子分析，確証的因子分析 

1
2
3
4
5
6
7
8
9
library(gdata)
library(psych)
library(lavaan)
data<-read.xls("http://www.educ.kyoto-u.ac.jp/cogpsy/personal/Kusumi/datasem11/osanai_factor.xls",1)
data<-data[,-1]
head(data)
#png("parallel01.png", width =600, height =600)  # 描画デバイスを開く
faev<-fa.parallel(data)
#dev.off()

因子数は2つとする

1
2
3
4
5
factML<- fa(r=data, nfactors=2 ,rotate="promax", fm="ml", scores=T)
factML
#群馬大学青木先生の関数
#source("http://aoki2.si.gunma-u.ac.jp/R/src/factanal2.R", encoding="euc-jp")
#factanal2(data,factors=2,scores="Bartlett")

ML1 ML2 h2 u2 com
積極性 0.83 0.06 0.74 0.26 1.0
陽気 -0.15 0.85 0.62 0.38 1.1
先導 0.75 0.02 0.59 0.41 1.0
無愛想 -0.03 -0.82 0.70 0.30 1.0
話好き 0.03 0.89 0.82 0.18 1.0
やる気 0.88 0.00 0.78 0.22 1.0
躊躇 -0.83 0.10 0.61 0.39 1.0
人気 0.23 0.59 0.54 0.46 1.3

第1因子は「活動性」
第2因子は「社交性」
と名づける。

1
2
3
4
5
6
7
8
9
10
11
12
13
model<-'
活動性 =~ 積極性 + 先導 + やる気 + 躊躇
社交性 =~ 陽気 + 無愛想 + 話好き + 人気
  '
fit <- cfa(model,data=data)
summary(fit, fit.measures=TRUE)
# 適合度指標の算出
inspect(fit,"fit")
# 修正インデックス
inspect(fit,"modindices")
#semPlot:::semPaths関数を使うと簡単にグラフ化できる
#library(semPlot)
#semPaths(fit,"std",style="lisrel",fade=F)

出力は省略

今回もgawkを使ってdot ファイルを出力

関数を定義(cat関数しか使ってない。）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
makeawk <- function(latent){
cat("BEGIN {",sep="\n",file="lavaan_dot.awk")
cat("i = 0",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"digraph  fit  \\{\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"rankdir=TB;\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"size=\\\"8,8\\\";\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"edge [fontname=\\\"sans\\\"  ,fontsize=8,arrowsize = 0.8,penwidth=0.8];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"graph [ordering = out,splines = true,overlap = false];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"center=1;\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"node [shape =ellipse];",sep=" ",append=T,file="lavaan_dot.awk")
cat(latent,sep=" ",append=T,file="lavaan_dot.awk")
cat(";\"",sep="\n",append=T,file="lavaan_dot.awk")
#cat("print \"切片1  [shape =triangle,fixedsize = true, width = 0.3, height = 0.3,fontsize=5,label=1];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"node [fontname=\\\"serif\\\" ,fontsize=14, shape=box];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"{rank=min };\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"=~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 \",dir=back\" $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 \",dir=both\" $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~1\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("i = i + 1",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"-> 切片\" i $3  $4 \",dir=back\" $5",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("END{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"\\}\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
}

潜在変数名をlatentにいれて、makeawk(latent)を実行する。
作業フォルダにlavaan_dot.awkができる。
あとは「sem_lavaanパッケージ」と同じ

1
2
3
4
5
6
latent<-c("活動性","社交性")
makeawk(latent)
pars <- parameterEstimates(fit, standardized = TRUE)
g<-data.frame(pars$lhs,pars$op,pars$rhs," [label=",signif(pars$std.all,digits=2),"];")
write.table(g,"out",sep=" ",col.names=FALSE,row.names=FALSE,quote=FALSE, na="")
system( 'gawk -f "lavaan_dot.awk" "out">out.dot')

out.dotをビューワーで確認。
活動性->活動性[label=1dir=both];
社交性->社交性[label=1dir=both];
をコメントアウトする。

1
system('dot -Tpng out.dot -o out.png')

makeawk関数を下のように変えるとnodeに着色する。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
makeawk <- function(latent){
cat("BEGIN {",sep="\n",file="lavaan_dot.awk")
cat("i = 0",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"digraph  fit  \\{\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"rankdir=TB;\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"size=\\\"8,8\\\";\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"edge [fontname=\\\"sans\\\" ,fontsize=8,arrowsize = 0.8,penwidth=0.8];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"graph [ordering = out,splines = true,overlap = false];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"center=1;\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"node [shape =ellipse, style = filled,fillcolor = \\\"#f4fd78\\\"];",sep=" ",append=T,file="lavaan_dot.awk")
cat(latent,sep=" ",append=T,file="lavaan_dot.awk")
cat(";\"",sep="\n",append=T,file="lavaan_dot.awk")
#cat("print \"切片1  [shape =triangle,fixedsize = true, width = 0.3, height = 0.3,fontsize=5,label=1, style = filled,fillcolor =pink];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"node [fontname=\\\"serif\\\" ,fontsize=14, shape=box, style = filled,fillcolor = \\\"#a9fab1\\\"];\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"{rank=min };\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"=~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 \",dir=back\" $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~~\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"->\" $3 $4 $5 \",dir=both\" $6",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("if($2 == \"~1\")",sep="\n",append=T,file="lavaan_dot.awk")
cat("{",sep="\n",append=T,file="lavaan_dot.awk")
cat("i = i + 1",sep="\n",append=T,file="lavaan_dot.awk")
cat("print  $1 \"-> 切片\" i $3  $4 \",dir=back\" $5",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
cat("END{",sep="\n",append=T,file="lavaan_dot.awk")
cat("print \"\\}\"",sep="\n",append=T,file="lavaan_dot.awk")
cat("}",sep="\n",append=T,file="lavaan_dot.awk")
}