データ解析B：第5回授業 因子分析と構造方程式モデリング（Structural Equation Modeling）
因子分析
研究の設計
- 研究の問い：人々の不安と格差に関する態度は、幸福度にどのような影響を与えるか？
- 利用する変数：不安感と格差／平等に関する因子、社会経済的属性に関する変数
- 利用する分析方法：因子分析とそれを含んだ構造方程式モデリング
因子分析とは何か？
主成分分析を思い出してみよう。
主成分分析の特徴：変数の相関係から、分散が最大となるように軸を設定。その軸に対して、無相関の軸を特定すること多変量の情報を縮約する
⇒因子分析も多変量の情報を縮約するという点では、主成分分析と同じ目的を持つ。しかしその理論的背景は、因子分析とは対照的といえるほどに異なる。
- 因子分析の特徴：相関をもとに、共通の因子を特定することで多変量の情報を縮約する
⇒言い換えるなら、複数の観測された多変量データの背後にある潜在的な因子(latent factor)を特定することが因子分析の目的。主成分を合成することで縮約しようとする主成分分析とは異なり、背景にある因子を絞り込むことで、多変量を縮約する発想に立っている。
- とはいえ、以下で見るように、主成分分析と因子分析の結果はかなりの程度類似したものとなることにも注目しよう。
因子分析の2つの方向性
探索的因子分析:多変量間に存在する複数の相関関係が、どの程度の数の、どのような潜在因子によって説明可能であるかを探索的に調べることを目的にした因子分析
確認的因子分析：因子数と観測された多変量との間に仮説を設定し、それが妥当であるかを確認することを目的にした因子分析
多くの場合、因子分析は前者の方向性を取ることが多い
因子分析における推定値
因子分析では、以下の模式図にあるように、観測されている多変量の背景にある共通因子と独自因子によって観測変数を説明しようとする。
観測変数＝共通因子＋独自因子 （共通因子も独自因子もともに潜在因子）
各観測変数は次式によって表されるが、以下の表記を*測定方程式と呼ぶ。
\[y_{防衛強化}=a_{1}f_{保守}+e_{防衛強化}\] \[y_{先制攻撃}=a_{1}f_{保守}+e_{防衛強化}\] \[\vdots\] \[y_{公共事業}=a_{1}f_{保守}+a_{2}f_{リベラル}+e_{公共事業} \] \[y_{公共事業}=a_{1}f_{保守}+a_{2}f_{リベラル}+e_{公共事業} \]
因子分析の手順: 主成分分析でも利用した「政策態度変数」を利用しよう
固有値ベクトルと平行分析
(2)固有値ベクトルを求める。
- 固有値ベクトルとは何か？
\[ Ax=\lambda x\] ベクトル\(x\)に関して、\(n\times n\)の行列を乗しても、長さが変化するだけで方向は同じとなる場合の\(\lambda\)が固有値。
(3)固有値ベクトルをもとに、スクリープロットすることで因子数を特定する。
(4)固有値が「1」以上の数の因子数の設定が妥当といわれるが、さらに正確に因子数を検討するために「平行分析（parallel analysis）」を実行。
## Parallel analysis suggests that the number of factors = 5 and the number of components = 3- 平行分析の解釈の仕方
PC Actual Data：主成分分析のための実際のデータを使った固有値のプロット
PC Simulated Data:主成分分析のためのシミュレーションされたデータを使った固有値のプロット
PC Resampled Data: 主成分分析のためのリサンプリングされたデータを使った固有値のプロット
FA Actual Data: 因子分析のための実際のデータを使った固有値のプロット
\(\vdots\)
平行分析では、実際のデータとリサンプリングされたデータの交点に当たるところの因子数をもって、適切な因子数であるとみなす。この場合、因子分析における適切な因子数は5となる。
因子分析と回転
(5)因子分析をfactanalコマンドによって実行するが、まずは初期解を求める。
factor_policy_init<-factanal(partyr,factors=5,scores = "Bartlett",rotation="none")
factor_policy_init##
## Call:
## factanal(x = partyr, factors = 5, scores = "Bartlett", rotation = "none")
##
## Uniquenesses:
## politician.Q4_1 politician.Q4_2 politician.Q4_3 politician.Q4_4
## 0.174 0.478 0.263 0.373
## politician.Q4_5 politician.Q4_6 politician.Q4_7 politician.Q4_8
## 0.377 0.561 0.484 0.445
## politician.Q4_9 politician.Q4_10 politician.Q4_11 politician.Q4_12
## 0.599 0.608 0.307 0.883
## politician.Q4_13 politician.Q4_14 politician.Q4_15 politician.Q4_16
## 0.242 0.213 0.354 0.694
## politician.Q4_17
## 0.662
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5
## politician.Q4_1 0.854 0.287 0.109
## politician.Q4_2 0.683 -0.129 0.172
## politician.Q4_3 0.789 0.242 0.233
## politician.Q4_4 -0.356 0.643 -0.288
## politician.Q4_5 0.769 0.165
## politician.Q4_6 0.562 0.312 -0.100 0.111
## politician.Q4_7 0.382 -0.488 0.309 0.166
## politician.Q4_8 0.291 -0.492 0.111 0.409 0.222
## politician.Q4_9 -0.444 0.191 0.393
## politician.Q4_10 -0.554 -0.215 0.182
## politician.Q4_11 0.793 0.168 -0.178
## politician.Q4_12 0.181 0.251 -0.126
## politician.Q4_13 0.804 0.198 0.177 -0.195
## politician.Q4_14 -0.765 0.334 0.125 0.263
## politician.Q4_15 -0.703 0.204 0.152 0.294
## politician.Q4_16 -0.318 0.344 0.225 0.181
## politician.Q4_17 0.271 0.377 -0.309 0.156
##
## Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings 6.119 1.183 0.908 0.645 0.427
## Proportion Var 0.360 0.070 0.053 0.038 0.025
## Cumulative Var 0.360 0.429 0.483 0.521 0.546
##
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 107.53 on 61 degrees of freedom.
## The p-value is 0.000221初期解をもとにすると因子負荷量は次のようにプロットされる。
しかし初期に得られたこの解よりも、軸の位置を「変化＝回転」させることで、より因子負荷量が明確になる軸を探すことできる。
(6)次にバリマックス回転(varimax rotation)による解を求める。
##
## Call:
## factanal(x = partyr, factors = 5, rotation = "varimax")
##
## Uniquenesses:
## politician.Q4_1 politician.Q4_2 politician.Q4_3 politician.Q4_4
## 0.174 0.478 0.263 0.373
## politician.Q4_5 politician.Q4_6 politician.Q4_7 politician.Q4_8
## 0.377 0.561 0.484 0.445
## politician.Q4_9 politician.Q4_10 politician.Q4_11 politician.Q4_12
## 0.599 0.608 0.307 0.883
## politician.Q4_13 politician.Q4_14 politician.Q4_15 politician.Q4_16
## 0.242 0.213 0.354 0.694
## politician.Q4_17
## 0.662
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5
## politician.Q4_1 0.863 -0.116 0.247
## politician.Q4_2 0.550 -0.313 0.183 0.274 0.113
## politician.Q4_3 0.848
## politician.Q4_4 0.320 -0.120 -0.698 -0.130
## politician.Q4_5 0.635 -0.210 0.294 0.235 -0.185
## politician.Q4_6 0.602 -0.103 0.242
## politician.Q4_7 0.214 -0.134 0.633 -0.170 -0.150
## politician.Q4_8 0.735
## politician.Q4_9 -0.279 0.562
## politician.Q4_10 -0.552 0.199 0.133 0.155
## politician.Q4_11 0.667 -0.391 0.210 0.170 0.151
## politician.Q4_12 0.328
## politician.Q4_13 0.701 -0.357 0.348 0.126
## politician.Q4_14 -0.495 0.486 -0.268 -0.169 0.454
## politician.Q4_15 -0.485 0.480 -0.115 -0.157 0.378
## politician.Q4_16 0.496 -0.178 0.128
## politician.Q4_17 0.292 -0.188 0.458
##
## Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings 4.473 1.638 1.428 1.105 0.638
## Proportion Var 0.263 0.096 0.084 0.065 0.038
## Cumulative Var 0.263 0.359 0.443 0.508 0.546
##
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 107.53 on 61 degrees of freedom.
## The p-value is 0.000221バリマックス回転のもとでの因子負荷量は次のようにプロットできる。
バリマックス回転では、軸を直交させた（直角）のまま回転させる。
軸の位置は回転により変化するが、各因子の相対的な位置は変化しない。また、初期解では有意な相関が認められなかった変数についても、相関が認められる場合もある。
(7)最後にプロマックス回転（斜交回転）による解を求める。
factor_policy_prom<-factanal(partyr,factors=5,rotation="promax",scores ="Bartlett")
factor_policy_prom##
## Call:
## factanal(x = partyr, factors = 5, scores = "Bartlett", rotation = "promax")
##
## Uniquenesses:
## politician.Q4_1 politician.Q4_2 politician.Q4_3 politician.Q4_4
## 0.174 0.478 0.263 0.373
## politician.Q4_5 politician.Q4_6 politician.Q4_7 politician.Q4_8
## 0.377 0.561 0.484 0.445
## politician.Q4_9 politician.Q4_10 politician.Q4_11 politician.Q4_12
## 0.599 0.608 0.307 0.883
## politician.Q4_13 politician.Q4_14 politician.Q4_15 politician.Q4_16
## 0.242 0.213 0.354 0.694
## politician.Q4_17
## 0.662
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5
## politician.Q4_1 0.909
## politician.Q4_2 0.465 -0.205 0.225
## politician.Q4_3 0.944 -0.164
## politician.Q4_4 0.158 0.248 -0.115 -0.781
## politician.Q4_5 0.544 0.222 0.116 -0.218
## politician.Q4_6 0.634 -0.171
## politician.Q4_7 0.147 0.630 -0.156 -0.103
## politician.Q4_8 0.169 0.829 0.178
## politician.Q4_9 0.701 0.279
## politician.Q4_10 -0.526 0.181 0.273 0.148
## politician.Q4_11 0.588 -0.291 0.105
## politician.Q4_12 0.132 0.330
## politician.Q4_13 0.668 -0.236 0.196 -0.111 0.122
## politician.Q4_14 -0.172 0.496 0.557
## politician.Q4_15 -0.192 0.529 0.486
## politician.Q4_16 0.256 0.595 -0.145 0.241
## politician.Q4_17 0.323 0.140 -0.169 0.386
##
## Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings 3.992 1.715 1.352 1.060 0.836
## Proportion Var 0.235 0.101 0.080 0.062 0.049
## Cumulative Var 0.235 0.336 0.415 0.478 0.527
##
## Factor Correlations:
## Factor1 Factor2 Factor3 Factor4 Factor5
## Factor1 1.000 -0.2527 0.3520 -0.2554 0.5294
## Factor2 -0.253 1.0000 0.1241 0.0527 -0.5218
## Factor3 0.352 0.1241 1.0000 -0.0732 0.1181
## Factor4 -0.255 0.0527 -0.0732 1.0000 0.0654
## Factor5 0.529 -0.5218 0.1181 0.0654 1.0000
##
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 107.53 on 61 degrees of freedom.
## The p-value is 0.000221プロマックス回転のもとでの因子負荷量は次のようにプロットできる。
因子分析の結果の解釈
(8)因子負荷量の解釈のために、負荷量の強弱を可視化しよう。因子負荷量を以下のコマンドで導出し、データフレーム化する。
loadings<-factor_policy_prom$loadings
#因子負荷量を抽出し、データフレームを作成
#[m,]はm列目、[,n]はn行目を表す
#また[m:n,]はm列目からn列目までを指定している
df_loadings<-data.frame(loadings[1:17,1], loadings[1:17,2],loadings[1:17,3],
loadings[1:17,4],loadings[1:17,5])
#行の名前の変数を追加
df_loadings$policy<-c("防衛強化","先制攻撃","北朝鮮","非核化","靖国","小さな政府",
"公共事業","景気対策","教育","課税強化","治安維持","外国人労働",
"原発容認","夫婦別姓","同性婚","選挙権","一院制")
#変数名を変換
colnames(df_loadings) <- c("第1因子","第2因子","第3因子","第4因子","第5因子","policy")tidyverseによって、tidy(きちんと成形された)データを構築しよう。
df_loadings.m <- df_loadings %>%
tidyr::gather(key = factor, value =loading, -policy)
#policy変数を基準にして
#Factorという新しい変数を作成
df_loadings.m #もとのdf_loadingとかたちを見比べてみてほしい## policy factor loading
## 1 防衛強化 第1因子 0.908971793
## 2 先制攻撃 第1因子 0.464757708
## 3 北朝鮮 第1因子 0.943556862
## 4 非核化 第1因子 0.158333294
## 5 靖国 第1因子 0.543748343
## 6 小さな政府 第1因子 0.633693064
## 7 公共事業 第1因子 0.146783359
## 8 景気対策 第1因子 -0.003248962
## 9 教育 第1因子 -0.041818340
## 10 課税強化 第1因子 -0.526072157
## 11 治安維持 第1因子 0.587754368
## 12 外国人労働 第1因子 0.091796823
## 13 原発容認 第1因子 0.668159202
## 14 夫婦別姓 第1因子 -0.172445408
## 15 同性婚 第1因子 -0.191567648
## 16 選挙権 第1因子 0.255858827
## 17 一院制 第1因子 0.323264590
## 18 防衛強化 第2因子 0.061491496
## 19 先制攻撃 第2因子 -0.204810755
## 20 北朝鮮 第2因子 0.072996740
## 21 非核化 第2因子 0.248359252
## 22 靖国 第2因子 -0.056300373
## 23 小さな政府 第2因子 0.048696409
## 24 公共事業 第2因子 0.006319035
## 25 景気対策 第2因子 0.169496633
## 26 教育 第2因子 0.700923087
## 27 課税強化 第2因子 0.180945248
## 28 治安維持 第2因子 -0.291473893
## 29 外国人労働 第2因子 0.063576378
## 30 原発容認 第2因子 -0.236059983
## 31 夫婦別姓 第2因子 0.495624061
## 32 同性婚 第2因子 0.529393901
## 33 選挙権 第2因子 0.594664954
## 34 一院制 第2因子 0.140281477
## 35 防衛強化 第3因子 -0.021241322
## 36 先制攻撃 第3因子 0.080906479
## 37 北朝鮮 第3因子 -0.007293509
## 38 非核化 第3因子 -0.115387535
## 39 靖国 第3因子 0.221624861
## 40 小さな政府 第3因子 -0.170533538
## 41 公共事業 第3因子 0.630349910
## 42 景気対策 第3因子 0.828697342
## 43 教育 第3因子 0.278582757
## 44 課税強化 第3因子 0.068968911
## 45 治安維持 第3因子 0.058565504
## 46 外国人労働 第3因子 -0.019851242
## 47 原発容認 第3因子 0.195919406
## 48 夫婦別姓 第3因子 -0.098146724
## 49 同性婚 第3因子 0.074210504
## 50 選挙権 第3因子 -0.037952021
## 51 一院制 第3因子 -0.169148146
## 52 防衛強化 第4因子 0.056115575
## 53 先制攻撃 第4因子 0.224693755
## 54 北朝鮮 第4因子 -0.164395289
## 55 非核化 第4因子 -0.781277557
## 56 靖国 第4因子 0.116409416
## 57 小さな政府 第4因子 0.079503777
## 58 公共事業 第4因子 -0.155939021
## 59 景気対策 第4因子 0.178099504
## 60 教育 第4因子 -0.036057268
## 61 課税強化 第4因子 0.272741826
## 62 治安維持 第4因子 0.105288197
## 63 外国人労働 第4因子 0.132158240
## 64 原発容認 第4因子 -0.111088315
## 65 夫婦別姓 第4因子 -0.065276855
## 66 同性婚 第4因子 -0.046168266
## 67 選挙権 第4因子 -0.144527694
## 68 一院制 第4因子 0.386124596
## 69 防衛強化 第5因子 -0.046088907
## 70 先制攻撃 第5因子 0.050451318
## 71 北朝鮮 第5因子 0.006299160
## 72 非核化 第5因子 0.030514580
## 73 靖国 第5因子 -0.217697875
## 74 小さな政府 第5因子 -0.095646698
## 75 公共事業 第5因子 -0.103313405
## 76 景気対策 第5因子 0.002171319
## 77 教育 第5因子 0.067530071
## 78 課税強化 第5因子 0.147777821
## 79 治安維持 第5因子 0.096850364
## 80 外国人労働 第5因子 0.329670835
## 81 原発容認 第5因子 0.122470041
## 82 夫婦別姓 第5因子 0.556507583
## 83 同性婚 第5因子 0.485704746
## 84 選挙権 第5因子 0.240881576
## 85 一院制 第5因子 0.020783302#各政策分野についての因子負荷量をまとめる
# 各バーの長さは因子負荷量の絶対値・強弱を表す
# 色は因子負荷量の正負の符号を表す
loading_plot<-ggplot(df_loadings.m, aes(policy, abs(loading), fill=loading)) +
facet_wrap(~ factor, nrow=1) + #「因子ごとに図を折り返す」の意味
geom_bar(stat="identity") +
coord_flip() + #図を縦方向に変換
#色を指定
scale_fill_gradient2(name = "loading",
high = "blue", mid = "white", low = "red",
midpoint=0, guide=F) +
ylab("因子負荷量の強弱") +
xlab("政策分野変数")+
theme_bw(base_size=10)
loading_plot
このレジュメの中で、上記の図は日本語上手く表示されていないけれども、画面では日本語の表記が確かめられるはず。
【第1因子】
正の相関（中でも太字は高い正の相関）
- Q4_１(防衛強化)
- Q4_2(先制攻撃)
- Q4_3(北朝鮮圧力)
- Q4_5(靖国参拝)
- Q4_6(小さな政府)
- Q4_7(公共事業)
- Q4_11(治安維持・プライバシー侵害)
- Q4_13(原発容認)
- Q4_17(一院制)
負の相関
- Q4_9(教育無償化)
- Q4_10(課税強化)
- Q4_14(夫婦別姓)
- Q4_15(同性婚)
⇒これらの相関関係をもとに考えると第一因子は「保守軸」と定義づけることが可能
【第2因子】
正の相関
- Q4_4(非核)
- Q4_9(教育無償化)
- Q4_10(課税強化)
- Q4_14(夫婦別姓)
- Q4_15(同性婚)
- Q4_16(選挙権年齢引き下げ)
負の相関
Q4_1(防衛強化)
Q4_2(先制攻撃)
Q4_5(靖国参拝)
Q4_7(公共事業)
Q4_11(治安維持)
Q4_13(原発容認) ⇒第2因子は「リベラル軸」と定義することが可能
この分析結果を前回の主成分分析の結果と比較検討してみよう。
###因子得点のプロット
(9)主成分分析の際に主成分得点をプロットしたように、因子得点もプロットしてみよう。
#データフレームの設定
politician<-read.csv("politician2017.csv", fileEncoding = "SJIS")
attach(politician)
partyr<-data.frame(politician$Q4_1,politician$Q4_2,politician$Q4_3,politician$Q4_4,politician$Q4_5,politician$Q4_6,politician$Q4_7,politician$Q4_8,politician$Q4_9,politician$Q4_10,politician$Q4_11,politician$Q4_12,politician$Q4_13,politician$Q4_14,politician$Q4_15,politician$Q4_16,politician$Q4_17,politician$NAME)
for(i in 1:17){
reason <- partyr[,i]
reason[partyr[,i] %in% 99] <- NA
reason[partyr[,i] %in% 5] <- 1
reason[partyr[,i] %in% 4] <- 2
reason[partyr[,i] %in% 3] <- 3
reason[partyr[,i] %in% 2] <- 4
reason[partyr[,i] %in% 1] <- 5
partyr[,i] <- reason
}
partyr<-na.omit(partyr)
df_fa.score<-data.frame(factor_policy_prom$scores[1:409,1],factor_policy_prom$scores[1:409,2],
factor_policy_prom$scores[1:409,3],factor_policy_prom$scores[1:409,4],
factor_policy_prom$scores[1:409,5],partyr$politician.NAME)
#mutate()関数を利用し、因子得点を含めたデータセットを構築
tb_fact<-mutate(df_fa.score, score1 = df_fa.score[,1], score2 = df_fa.score[, 2])
tb_fact<-na.omit(tb_fact)
#政治家名を明示しながら、因子得点をプロット
library(ggrepel)
library(ggplot2)
#散布図をプロット
factscore_fig<-ggplot(tb_fact, aes(x=tb_fact$score1, y=tb_fact$score2,label=tb_fact$partyr.politician.NAME))
factscore_fig <- factscore_fig + geom_point() +
geom_text(aes(y=tb_fact$score2 + .1, label=tb_fact$partyr.politician.NAME),size = 2, vjust=0)+
xlab("第1因子：保守軸")+
ylab("第2因子：リベラル軸")
plot(factscore_fig)
- 主成分分析の回のレジュメを参考に、党首名だけのプロットを作ってみよう。
因子分析の実践的学習：人々の「不安感」と「格差感」
データと変数の設定
(1)「A1」の質問項目をもとに、変数の設定をしていこう。
df_fac1<-data.frame(qa1_01, qa1_02, qa1_03, qa1_04, qa1_05, qa1_06, qa1_07, qa1_08, qa1_09, qa1_10, qa1_11, qa1_12, qa1_13,qa1_14, qa1_15, qa1_16, qa1_17, qa1_18, qa1_19, qa1_20, qa1_21, qa1_22, qa1_23, qa1_24)
for(i in 1:24){
fac <- df_fac1[,i]
fac[df_fac1[,i] %in% 99] <- NA
fac[df_fac1[,i] %in% 5] <- 1
fac[df_fac1[,i] %in% 4] <- 2
fac[df_fac1[,i] %in% 3] <- 3
fac[df_fac1[,i] %in% 2] <- 4
fac[df_fac1[,i] %in% 1] <- 5
df_fac1[,i] <- fac}固有値ベクトルと平行分析
(2)平行分析を実行し、適切な因子数を割り出そう。
## Parallel analysis suggests that the number of factors = 7 and the number of components = 4因子分析の実行
(3)因子分析をfactanalコマンドによって実行し、各種の回転を試行してよいと思われる結果を特定していこう。
##
## Call:
## factanal(x = df_fac1, factors = 7, scores = "Bartlett", rotation = "promax")
##
## Uniquenesses:
## qa1_01 qa1_02 qa1_03 qa1_04 qa1_05 qa1_06 qa1_07 qa1_08 qa1_09 qa1_10 qa1_11
## 0.481 0.819 0.771 0.005 0.490 0.729 0.656 0.747 0.679 0.800 0.178
## qa1_12 qa1_13 qa1_14 qa1_15 qa1_16 qa1_17 qa1_18 qa1_19 qa1_20 qa1_21 qa1_22
## 0.345 0.495 0.351 0.621 0.865 0.713 0.686 0.824 0.285 0.421 0.680
## qa1_23 qa1_24
## 0.681 0.612
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 Factor7
## qa1_01 -0.360 0.191 0.557
## qa1_02 0.316 0.188 0.159
## qa1_03 0.235 -0.118 0.136 0.230
## qa1_04 1.013
## qa1_05 0.776 -0.216
## qa1_06 0.496 0.160
## qa1_07 0.164 -0.294 0.192 0.105 0.495 0.131
## qa1_08 0.605
## qa1_09 0.370 0.210 0.113 -0.104
## qa1_10 0.147 0.139 0.323
## qa1_11 0.589 0.667
## qa1_12 0.712 0.319
## qa1_13 0.723
## qa1_14 0.832
## qa1_15 0.192 0.277 0.217 -0.194
## qa1_16 0.140 0.272 -0.114
## qa1_17 0.131 0.334 -0.100 0.170
## qa1_18 0.196 0.371
## qa1_19 0.183 0.392
## qa1_20 0.835 0.133
## qa1_21 -0.114 0.741
## qa1_22 -0.108 0.129 0.464 0.148
## qa1_23 -0.108 0.113 0.550 -0.145
## qa1_24 -0.120 -0.133 0.635 0.202
##
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 Factor7
## SS loadings 2.546 1.962 1.532 1.133 1.070 1.000 0.649
## Proportion Var 0.106 0.082 0.064 0.047 0.045 0.042 0.027
## Cumulative Var 0.106 0.188 0.252 0.299 0.343 0.385 0.412
##
## Factor Correlations:
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 Factor7
## Factor1 1.000 0.3903 0.1654 0.411 0.1683 -0.0970 0.2674
## Factor2 0.390 1.0000 0.1895 0.586 0.0990 0.0311 0.4269
## Factor3 0.165 0.1895 1.0000 0.228 0.0967 0.3488 0.5106
## Factor4 0.411 0.5858 0.2283 1.000 0.2064 0.2781 0.5857
## Factor5 0.168 0.0990 0.0967 0.206 1.0000 -0.2481 0.0257
## Factor6 -0.097 0.0311 0.3488 0.278 -0.2481 1.0000 0.3632
## Factor7 0.267 0.4269 0.5106 0.586 0.0257 0.3632 1.0000
##
## Test of the hypothesis that 7 factors are sufficient.
## The chi square statistic is 401.63 on 129 degrees of freedom.
## The p-value is 9.51e-30- これらの因子はどのような性質のものとして定義できるだろうか？
(4)次に、平等／格差感にかかわる因子(A2)を質問2群の中から特定していこう。まずは上記と同じように変数を設定する。
(5)平行分析の後に、因子分析を実行し結果について解釈しよう。
## Parallel analysis suggests that the number of factors = 6 and the number of components = 4##
## Call:
## factanal(x = df_fac2, factors = 6, scores = "Bartlett", rotation = "promax")
##
## Uniquenesses:
## qa2_01 qa2_02 qa2_03 qa2_04 qa2_05 qa2_06 qa2_07 qa2_08 qa2_09 qa2_10 qa2_11
## 0.710 0.552 0.648 0.683 0.746 0.558 0.482 0.367 0.490 0.237 0.252
## qa2_12 qa2_13 qa2_14 qa2_15 qa2_16
## 0.705 0.597 0.565 0.468 0.643
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
## qa2_01 -0.103 0.600
## qa2_02 -0.119 0.752
## qa2_03 0.371 -0.105 0.360
## qa2_04 0.361 0.255
## qa2_05 0.102 0.452
## qa2_06 0.730 -0.168
## qa2_07 0.524 0.101 0.213 -0.152 -0.162
## qa2_08 0.931 -0.181
## qa2_09 0.419 0.198 0.185
## qa2_10 0.938
## qa2_11 0.837
## qa2_12 0.141 0.384 -0.280
## qa2_13 0.186 0.157 0.593
## qa2_14 -0.177 0.142 0.664 0.119
## qa2_15 0.138 -0.105 0.639
## qa2_16 0.659 0.192
##
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
## SS loadings 2.064 1.685 1.425 0.977 0.783 0.574
## Proportion Var 0.129 0.105 0.089 0.061 0.049 0.036
## Cumulative Var 0.129 0.234 0.323 0.384 0.433 0.469
##
## Factor Correlations:
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
## Factor1 1.000 0.607 0.4415 0.291 -0.228 0.3595
## Factor2 0.607 1.000 0.5494 0.402 -0.276 0.3941
## Factor3 0.441 0.549 1.0000 0.178 0.064 0.0532
## Factor4 0.291 0.402 0.1783 1.000 -0.163 0.3446
## Factor5 -0.228 -0.276 0.0640 -0.163 1.000 -0.3834
## Factor6 0.359 0.394 0.0532 0.345 -0.383 1.0000
##
## Test of the hypothesis that 6 factors are sufficient.
## The chi square statistic is 86.62 on 39 degrees of freedom.
## The p-value is 1.8e-05構造方程式モデリング(Structural Equation Modeling: SEM)
構造方程式モデリングとは何か
複数の変数間、つまり多変量間の因果性・相関性を検証可能なモデル
因子の特定もモデルの中に含めることで、潜在因子（＝構成概念）も含めたモデル化が可能
このモデルを、パス図(paty diagram)によって明示することが可能
→ここで今回の授業における「研究の問い」を思い出してみよう。「不安や格差／平等感は幸福度に影響を与えているのか？」を分析することが目的。これを図で表すと、おおよそ以下のようになるが・・・
- SEMはまさに上記のような模式図を「モデル」として、これをそのまま推定することが特徴
SEMにおける用語とモデルの記法（決まりごと）
用語
観測変数：直接的に観察された変数
潜在変数：直接的に観察されない変数であり、因子を指す。構成概念を意味する
誤差変数：結果側の変動生み出す変数であり、原因となる変数によっては説明されない変動を表す変数
外生変数：何からも影響を受けない変数。モデルの中では、他の変数からの矢印が入らない変数
内生変数：原因であり結果である、双方向の因果性がある変数。モデルの中では、その変数から発する矢印も、その変数に入る矢印も双方が存在する
モデルの記法
観測変数は四角で囲む
潜在変数（因子）は丸または楕円で囲む
誤差変数は丸または楕円で囲むか、何でも囲まずに表記する
矢印の始点は原因、矢印の終点は結果を表す
双方向矢印は相関を表す
片方矢印を受けた結果側の変数には、必ず誤差変数が付属する
構造方程式モデリングの手順：幸福度をめぐる因果関係を検証しよう
構造方程式モデリングの準備、package「lavaan」の利用
(1)上記のように、分析者が何らかのモデルを設定する。例えば私は、こんなモデルを当初は考え、手書きしていた。
(2)SEMにはいくつかのパッケージが必要になるので、あらかじめ読みこむ。
(3)avaanパッケージによってSEMを実行する場合、大きく分けて以下の3段階のプロセスからなる。
モデルの設定
モデルのあてはめ(fit)
分析結果をもとにしたパス図の出力
第1段階目：モデルの設定
(4)モデルを設定する
- 回帰モデルの場合は従属変数と独立変数群を、以下のように\(\sim\)（チルダ）でつなぐ
\[ y\sim x1+x2+x3 \]
\[ x1\sim f1+f2+f3 \]
- 潜在変数についてのモデルは、以下のように\(\sim=\)で表す
\[ f1=\sim x4+x5+x6 \]
- 分散共分散については、以下のように\(\sim\sim\)で表す
\[ x1\sim\sim x1 \] \[ x2\sim\sim x3 \]
- モデルの設定は、モデル名のオブジェクトを作成しその中にモデルを格納するという手順
#幸福度
happy<-qa10
happy[qa10==99]<-NA
#性別
gender<-qb1
gender[qb1==2]<-0
#年齢
age<-2018-qb3y
#失業
unemployment<-qb8
unemployment[qb8==888888|qb8==999999]<-0
unemployment[unemployment>0]<-1
#非正規雇用
albeit<-qb9_2
albeit[qb9_2==1]<-0
albeit[qb9_2==2|qb9_2==3| qb9_2==4|qb9_2==5]<-1
albeit[qb9_2==6|qb9_2==8|qb9_2==9]<-NA
#収入
income<-qb11
income[qb11==99]<-NA
#結婚
marriage<-qb2
marriage[qb2==1]<-1
marriage[qb2==2]<-1
marriage[qb2!=1]<-0
#離婚
divorce<-qb2
divorce[qb2==3]<-"di"
divorce[divorce!="di"]<-0
divorce[divorce=="di"]<-1
divorce<-as.numeric(divorce)
#子ども
qb18a[qb18a==99]<-NA
qb18b[qb18b==99]<-NA
children<-qb18a+qb18b
#1人親
single<-qb15
single[qb15==1]<-11
single[qb15==5]<-1
single[single!=1]<-0sem_model<-"
#潜在変数モデル1
anxiety =~ qa1_01+qa1_07+qa1_09+qa1_10+qa1_11+qa1_12+qa1_12+
qa1_13+qa1_14+qa1_15+qa1_16+qa1_17+qa1_18+qa1_22+qa1_23
#潜在変数モデル2
disparity =~ qa2_01+qa2_02+qa2_04+qa2_06+qa2_07+qa2_08+qa2_09+
qa2_14+qa2_15
#回帰モデル
happy~ anxiety+disparity+gender+age+unemployment+albeit+
income+marriage+divorce+children+single
"
#モデルの最初と最後は「'」(shift+7)でくくる第2段階目：モデルのあてはめ
(5)先に定義したモデルを、構造方程式モデリングにあてはめる
#データフレームの設定
df_sem<-data.frame(happy,gender,age,unemployment,albeit,income,marriage,
divorce,children,single,qa1_01,qa1_07,qa1_09,qa1_10,
qa1_11,qa1_12,qa1_13,qa1_14,qa1_15,qa1_16,qa1_17,
qa1_18,qa1_22,qa1_23,qa2_01,qa2_02,qa2_04,qa2_06,
qa2_07,qa2_08,qa2_09,qa2_14,qa2_15)
df_sem<-na.omit(df_sem)
#あてはめ
fit <- lavaan::sem(sem_model, data = df_sem)
summary(fit)## lavaan 0.6-12 ended normally after 232 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 59
##
## Number of observations 679
##
## Model Test User Model:
##
## Test statistic 2204.773
## Degrees of freedom 457
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## anxiety =~
## qa1_01 1.000
## qa1_07 0.186 0.102 1.826 0.068
## qa1_09 -1.355 0.160 -8.453 0.000
## qa1_10 -0.256 0.116 -2.194 0.028
## qa1_11 -2.375 0.227 -10.470 0.000
## qa1_12 -2.842 0.262 -10.845 0.000
## qa1_13 -1.858 0.194 -9.576 0.000
## qa1_14 -2.020 0.201 -10.061 0.000
## qa1_15 -1.648 0.189 -8.724 0.000
## qa1_16 -0.152 0.170 -0.898 0.369
## qa1_17 -0.626 0.112 -5.604 0.000
## qa1_18 -0.831 0.123 -6.756 0.000
## qa1_22 0.311 0.115 2.711 0.007
## qa1_23 -0.564 0.101 -5.591 0.000
## disparity =~
## qa2_01 1.000
## qa2_02 5.836 8.258 0.707 0.480
## qa2_04 13.421 18.825 0.713 0.476
## qa2_06 15.434 21.621 0.714 0.475
## qa2_07 22.372 31.322 0.714 0.475
## qa2_08 21.735 30.430 0.714 0.475
## qa2_09 21.848 30.594 0.714 0.475
## qa2_14 8.420 11.860 0.710 0.478
## qa2_15 15.509 21.742 0.713 0.476
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## happy ~
## anxiety -2.511 0.290 -8.653 0.000
## disparity 0.800 3.058 0.262 0.794
## gender -0.366 0.137 -2.672 0.008
## age -0.007 0.006 -1.219 0.223
## unemployment 2.472 1.468 1.684 0.092
## albeit 0.212 0.151 1.401 0.161
## income 0.232 0.047 4.909 0.000
## marriage 0.369 0.202 1.826 0.068
## divorce 0.180 0.390 0.463 0.644
## children -0.038 0.066 -0.572 0.567
## single -0.193 0.367 -0.526 0.599
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## anxiety ~~
## disparity -0.003 0.005 -0.711 0.477
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .qa1_01 0.590 0.033 17.920 0.000
## .qa1_07 0.762 0.041 18.412 0.000
## .qa1_09 0.845 0.048 17.776 0.000
## .qa1_10 0.986 0.054 18.406 0.000
## .qa1_11 0.509 0.034 14.932 0.000
## .qa1_12 0.296 0.029 10.137 0.000
## .qa1_13 0.781 0.046 17.089 0.000
## .qa1_14 0.605 0.037 16.359 0.000
## .qa1_15 1.074 0.061 17.668 0.000
## .qa1_16 2.158 0.117 18.422 0.000
## .qa1_17 0.708 0.039 18.262 0.000
## .qa1_18 0.731 0.040 18.145 0.000
## .qa1_22 0.935 0.051 18.395 0.000
## .qa1_23 0.578 0.032 18.263 0.000
## .qa2_01 0.653 0.035 18.420 0.000
## .qa2_02 0.512 0.028 18.188 0.000
## .qa2_04 0.596 0.034 17.342 0.000
## .qa2_06 0.369 0.023 16.093 0.000
## .qa2_07 0.420 0.030 14.074 0.000
## .qa2_08 0.363 0.027 13.671 0.000
## .qa2_09 0.518 0.034 15.074 0.000
## .qa2_14 0.626 0.035 18.021 0.000
## .qa2_15 0.616 0.036 17.023 0.000
## .happy 2.023 0.116 17.471 0.000
## anxiety 0.126 0.023 5.396 0.000
## disparity 0.001 0.002 0.357 0.721第3段階目：分析結果の表への出力
(6)分析結果を表として見やすいように出力する。semTalbleというパッケージの中の```SemTable()’’’というコマンドを使ってみよう。以下のようにコマンドを入力すると、難解なコードが出力されてくる。
library(semTable)
fit1.t1 <- semTable (fit , type="latex",columns = c("estsestars", "p"),
columnLabels = c(estsestars = "Est w/stars", p = "p-value"),longtable=TRUE)## \begin{longtable}{@{}r*{2}{S[
## input-symbols = ( ) +,
## group-digits = false,
## table-number-alignment = center,
## %table-space-text-pre = (,
## table-align-text-pre = false,
## table-align-text-post = false,
## table-space-text-post = {***},
## parse-units = false]}@{}}
##
## \endfirsthead
## \endhead
## & \multicolumn{2}{c}{Model}\tabularnewline \hline
## & \multicolumn{1}{c}{Est w/stars}& \multicolumn{1}{c}{p-value}\tabularnewline\hline
## & \multicolumn{2}{c}{\underline{Factor Loadings}}\tabularnewline \multicolumn{1}{l}{\underline{anxiety}}\tabularnewline
## qa1.01& 1.00$^+$& \tabularnewline
## qa1.07& 0.19{(0.10)}$\phantom{{^{***}}}$& .068\tabularnewline
## qa1.09& -1.36{(0.16)}$^{***}$& .000\tabularnewline
## qa1.10& -0.26{(0.12)}$^{*}\phantom{{^{**}}}$& .028\tabularnewline
## qa1.11& -2.38{(0.23)}$^{***}$& .000\tabularnewline
## qa1.12& -2.84{(0.26)}$^{***}$& .000\tabularnewline
## qa1.13& -1.86{(0.19)}$^{***}$& .000\tabularnewline
## qa1.14& -2.02{(0.20)}$^{***}$& .000\tabularnewline
## qa1.15& -1.65{(0.19)}$^{***}$& .000\tabularnewline
## qa1.16& -0.15{(0.17)}$\phantom{{^{***}}}$& .369\tabularnewline
## qa1.17& -0.63{(0.11)}$^{***}$& .000\tabularnewline
## qa1.18& -0.83{(0.12)}$^{***}$& .000\tabularnewline
## qa1.22& 0.31{(0.11)}$^{**}\phantom{{^{*}}}$& .007\tabularnewline
## qa1.23& -0.56{(0.10)}$^{***}$& .000\tabularnewline
## \multicolumn{1}{l}{\underline{disparity}}\tabularnewline
## qa2.01& 1.00$^+$& \tabularnewline
## qa2.02& 5.84{(8.26)}$\phantom{{^{***}}}$& .480\tabularnewline
## qa2.04& 13.42{(18.82)}$\phantom{{^{***}}}$& .476\tabularnewline
## qa2.06& 15.43{(21.62)}$\phantom{{^{***}}}$& .475\tabularnewline
## qa2.07& 22.37{(31.32)}$\phantom{{^{***}}}$& .475\tabularnewline
## qa2.08& 21.74{(30.43)}$\phantom{{^{***}}}$& .475\tabularnewline
## qa2.09& 21.85{(30.59)}$\phantom{{^{***}}}$& .475\tabularnewline
## qa2.14& 8.42{(11.86)}$\phantom{{^{***}}}$& .478\tabularnewline
## qa2.15& 15.51{(21.74)}$\phantom{{^{***}}}$& .476\tabularnewline
## & \multicolumn{2}{c}{\underline{Regression Slopes}}\tabularnewline \multicolumn{1}{l}{\underline{happy}}\tabularnewline
## anxiety& -2.51{(0.29)}$^{***}$& .000\tabularnewline
## disparity& 0.80{(3.06)}$\phantom{{^{***}}}$& .794\tabularnewline
## gender& -0.37{(0.14)}$^{**}\phantom{{^{*}}}$& .008\tabularnewline
## age& -0.01{(0.01)}$\phantom{{^{***}}}$& .223\tabularnewline
## unemployment& 2.47{(1.47)}$\phantom{{^{***}}}$& .092\tabularnewline
## albeit& 0.21{(0.15)}$\phantom{{^{***}}}$& .161\tabularnewline
## income& 0.23{(0.05)}$^{***}$& .000\tabularnewline
## marriage& 0.37{(0.20)}$\phantom{{^{***}}}$& .068\tabularnewline
## divorce& 0.18{(0.39)}$\phantom{{^{***}}}$& .644\tabularnewline
## children& -0.04{(0.07)}$\phantom{{^{***}}}$& .567\tabularnewline
## single& -0.19{(0.37)}$\phantom{{^{***}}}$& .599\tabularnewline
## & \multicolumn{2}{c}{\underline{Residual Variances}}\tabularnewline
## qa1.01& 0.59{(0.03)}$^{***}$& .000\tabularnewline
## qa1.07& 0.76{(0.04)}$^{***}$& .000\tabularnewline
## qa1.09& 0.84{(0.05)}$^{***}$& .000\tabularnewline
## qa1.10& 0.99{(0.05)}$^{***}$& .000\tabularnewline
## qa1.11& 0.51{(0.03)}$^{***}$& .000\tabularnewline
## qa1.12& 0.30{(0.03)}$^{***}$& .000\tabularnewline
## qa1.13& 0.78{(0.05)}$^{***}$& .000\tabularnewline
## qa1.14& 0.61{(0.04)}$^{***}$& .000\tabularnewline
## qa1.15& 1.07{(0.06)}$^{***}$& .000\tabularnewline
## qa1.16& 2.16{(0.12)}$^{***}$& .000\tabularnewline
## qa1.17& 0.71{(0.04)}$^{***}$& .000\tabularnewline
## qa1.18& 0.73{(0.04)}$^{***}$& .000\tabularnewline
## qa1.22& 0.93{(0.05)}$^{***}$& .000\tabularnewline
## qa1.23& 0.58{(0.03)}$^{***}$& .000\tabularnewline
## qa2.01& 0.65{(0.04)}$^{***}$& .000\tabularnewline
## qa2.02& 0.51{(0.03)}$^{***}$& .000\tabularnewline
## qa2.04& 0.60{(0.03)}$^{***}$& .000\tabularnewline
## qa2.06& 0.37{(0.02)}$^{***}$& .000\tabularnewline
## qa2.07& 0.42{(0.03)}$^{***}$& .000\tabularnewline
## qa2.08& 0.36{(0.03)}$^{***}$& .000\tabularnewline
## qa2.09& 0.52{(0.03)}$^{***}$& .000\tabularnewline
## qa2.14& 0.63{(0.03)}$^{***}$& .000\tabularnewline
## qa2.15& 0.62{(0.04)}$^{***}$& .000\tabularnewline
## happy& 2.02{(0.12)}$^{***}$& .000\tabularnewline
## gender& 0.25$^+$& \tabularnewline
## age& 104.94$^+$& \tabularnewline
## unemployment& 0.00$^+$& \tabularnewline
## albeit& 0.24$^+$& \tabularnewline
## income& 2.67$^+$& \tabularnewline
## marriage& 0.14$^+$& \tabularnewline
## divorce& 0.03$^+$& \tabularnewline
## children& 1.14$^+$& \tabularnewline
## single& 0.03$^+$& \tabularnewline
## & \multicolumn{2}{c}{\underline{Residual Covariances}}\tabularnewline
## gender w/age& 0.35$^+$& \tabularnewline
## gender w/unemployment& 0.00$^+$& \tabularnewline
## gender w/albeit& -0.10$^+$& \tabularnewline
## gender w/income& 0.42$^+$& \tabularnewline
## gender w/marriage& 0.01$^+$& \tabularnewline
## gender w/divorce& -0.01$^+$& \tabularnewline
## gender w/children& -0.08$^+$& \tabularnewline
## gender w/single& -0.01$^+$& \tabularnewline
## age w/unemployment& 0.01$^+$& \tabularnewline
## age w/albeit& 1.04$^+$& \tabularnewline
## age w/income& -1.59$^+$& \tabularnewline
## age w/marriage& 0.83$^+$& \tabularnewline
## age w/divorce& -0.19$^+$& \tabularnewline
## age w/children& 2.58$^+$& \tabularnewline
## age w/single& 0.03$^+$& \tabularnewline
## unemployment w/albeit& 0.00$^+$& \tabularnewline
## unemployment w/income& -0.00$^+$& \tabularnewline
## unemployment w/marriage& 0.00$^+$& \tabularnewline
## unemployment w/divorce& -0.00$^+$& \tabularnewline
## unemployment w/children& 0.00$^+$& \tabularnewline
## unemployment w/single& -0.00$^+$& \tabularnewline
## albeit w/income& -0.49$^+$& \tabularnewline
## albeit w/marriage& 0.00$^+$& \tabularnewline
## albeit w/divorce& 0.00$^+$& \tabularnewline
## albeit w/children& 0.06$^+$& \tabularnewline
## albeit w/single& 0.01$^+$& \tabularnewline
## income w/marriage& 0.04$^+$& \tabularnewline
## income w/divorce& 0.00$^+$& \tabularnewline
## income w/children& -0.06$^+$& \tabularnewline
## income w/single& -0.01$^+$& \tabularnewline
## marriage w/divorce& -0.03$^+$& \tabularnewline
## marriage w/children& 0.20$^+$& \tabularnewline
## marriage w/single& -0.02$^+$& \tabularnewline
## divorce w/children& 0.01$^+$& \tabularnewline
## divorce w/single& 0.02$^+$& \tabularnewline
## children w/single& 0.02$^+$& \tabularnewline
## & \multicolumn{2}{c}{\underline{Latent Variances}}\tabularnewline
## anxiety& 0.13{(0.02)}$^{***}$& .000\tabularnewline
## disparity& 0.00{(0.00)}$\phantom{{^{***}}}$& .721\tabularnewline
## & \multicolumn{2}{c}{\underline{Latent Covariances}}\tabularnewline
## anxiety w/disparity& -0.00{(0.00)}$\phantom{{^{***}}}$& .477\tabularnewline
## & \multicolumn{2}{c}{\underline{Fit Indices}}\tabularnewline
## $\chi^{2}(\mathrm{df})$& 2204.77(457)$^{***}$& .000\tabularnewline
## CFI& 0.59& \tabularnewline
## TLI& 0.56& \tabularnewline
## RMSEA& 0.08& \tabularnewline
## \hline\multicolumn{3}{l}{$^+$Fixed parameter}\tabularnewline\multicolumn{3}{l}{$^{*}\phantom{{^{**}}}$p<0.05, $^{**}\phantom{{^{*}}}$p<0.01, $^{***}$p<0.001}\tabularnewline
## \end{longtable}
##
## - このコードは、「LaTex」(テフ／ラテックス)と呼ばれる文書作成ソフト上で表を編集（compile）するためのコードである。そして、このレジュメは、Rmarkdownでの冒頭部分yamlで、次のように文書の出力をpdf_documentに設定している。
LaTexでは数式を上手に表記できたり、見やすい分析結果の表を作成することができる。Latexをtinylatexというパッケージを使って、Rmarkdown上に連動させながら、pdf_documentを作成していく。ここでは、見やすい各種分析結果の表を
semTable()によって出力し、pdfに出力している。以下では、実際に編集(compile)した分析結果の表を挙げているので、上のコードと見比べてみてほしい。表にする際に、有意水準をアスタリスクで表すように指示している。また表が長いので、引数の中で、
longtable=TRUEを指定し、長い表に対応したものになっている。
第4段階目：パス図の作成
(7)パス図によって因果性などを可視化する。
semPaths(fit,"par", style="lisrel",
mar=c(6,1,3,1), # 余白の指定、下、左、上、右の順
edge.label.cex=.8, # 矢印の係数の文字の大きさ
fade=F,
theme = 'gray')
いまいちわかりにくいので・・・
(8)少し難しいけれども、自作の関数をもとにDiagrammeRパッケージを用いて、見やすい図を作成する。DiagrammeRで出力されてくるのは、.dotファイル形式のものなので、これを.png形式に変換してレポートやペーパーに挿入できるように変換する作業まで頑張ってみよう。
#"lavaan_plot.R"というスクリプト・ファイルにはDiagrammeRを
#使ってより分かりやすい図を出力するための関数my_lavaan_plot()を
#作成してある。その関数をsource()関数を使って呼び出す
source("lavaan_plot.R")
#自作の関数をもとに、fitに格納した分析結果をパス図として出力する
fit_plot<-my_lavaan_plot(fit)
#上記の操作の時点でviewerで図を確認できるが、その図を.pngファイルで出力したい。
#以下のパッケージは図をDiagrammeRで出力した図を
#.png形式で出力するために必要な関数を実装したパッケージ
library(DiagrammeRsvg)
library(magrittr)
library(svglite)
library(rsvg)
library(png)
#fit_plotをpngファイルに変換する
export_svg(fit_plot) %>%
charToRaw %>% rsvg %>% png::writePNG('sem_diagramme.png')モデルの再設定・修正と再解釈
(9)モデルの当てはまりが悪いので、それを改善するためにdiparity変数を除外し、潜在変数モデル内で統計的に有意ではない変数を除外する。
sem_model_anx<-"
#潜在変数モデルを不安感に関するもののみに絞る
anxiety =~ qa1_11+qa1_12+qa1_12+qa1_13+qa1_14
#回帰モデル
happy~ anxiety+gender+age+unemployment+
income+marriage+children
"
#あてはめ
fit_anxiety <- lavaan::sem(sem_model_anx, data = df_sem)
summary(fit_anxiety)## lavaan 0.6-12 ended normally after 42 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 16
##
## Number of observations 679
##
## Model Test User Model:
##
## Test statistic 201.189
## Degrees of freedom 29
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## anxiety =~
## qa1_11 1.000
## qa1_12 1.248 0.057 21.831 0.000
## qa1_13 0.765 0.052 14.841 0.000
## qa1_14 0.838 0.049 17.129 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## happy ~
## anxiety 1.013 0.079 12.762 0.000
## gender -0.425 0.137 -3.100 0.002
## age -0.006 0.006 -1.065 0.287
## unemployment 2.637 1.492 1.767 0.077
## income 0.220 0.041 5.318 0.000
## marriage 0.395 0.175 2.251 0.024
## children -0.029 0.063 -0.463 0.644
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .qa1_11 0.513 0.036 14.371 0.000
## .qa1_12 0.213 0.034 6.221 0.000
## .qa1_13 0.802 0.047 17.146 0.000
## .qa1_14 0.623 0.038 16.338 0.000
## .happy 2.106 0.120 17.600 0.000
## anxiety 0.704 0.064 11.018 0.000library(semTable)
fit.anxiety <- semTable (fit_anxiety , type="latex",columns = c("estsestars", "p"),
columnLabels = c(estsestars = "Est w/stars", p="p-value"),longtable=TRUE)lavaanのsemplot()によるパス図。変数が減ったことで、少し見やすくなった。
semPaths(fit_anxiety,"par", style="lisrel",
mar=c(6,1,3,1), # 余白の指定、下、左、上、右の順
edge.label.cex=.8, # 矢印の係数の文字の大きさ
fade=F,
theme = 'gray')
自作関数may_lavaan_plot()によるパス図は以下の通り。
fit_plot_anx<-my_lavaan_plot(fit_anxiety)
#fit_plotをpngファイルに変換する
export_svg(fit_plot_anx) %>%
charToRaw %>% rsvg %>% png::writePNG('sem_diagramme_anx.png')- 再分析の結果以下の変数について、幸福度に対してどのような解釈できるだろうか。
- 不安感
- 性別
- 無職
- 所得
- 既婚
- モデルの当てはまり具合は、どの程度改善されただろうか？
構造方程式モデリング実践的学習
(1)本授業の受講生の皆さんは有権者です。「有権者調査コードブック」を一読しながら、有権者の投票選択（安倍政権下の自民党）を決める要因について、各自でモデルを組み立てください。
(2)モデルを「モデルの設定」に従って設定してください。
(3)semを実行してみましょう。
(4)表を出力してください。
(5)パス図を書いてください。