R言語で構造方程式モデリング その1

R

行動科学などの領域では、重さや長さのように直接観測可能な特性だけでなく、「先高感」「知能」「ブランドイメージ」のように直接測定出来ない特性である構成概念（潜在変数）を扱う必要がある。
構造方程式モデリングは、構成概念や観測変数の性質を調べるために集めた多くの観測変数を同時に分析するための統計的方法。

構造方程式モデリング

  • 多変量解析の1手法
  • 第2世代の多変量解析
  • 伝統的な統計手法を下位モデルとして実行可能
    • 因子分析法
    • 多変量解析
    • パス解析法
    • 時系列解析
    • etc...

変数の分類

方程式モデルには、様々な種類の変数が登場するため、ここでは、3つの観点から分類しその性質を確認する

  • 観測変数と潜在変数
    • 観測変数 : 直接観測出来る変数。
    • 潜在変数 : モデルに導入した直接観測出来ない変数。構成概念
  • 構造変数と誤差変数
    • 構造変数 : 考察の対象となっている構成概念を含んでいる変数
    • 誤差変数 : 考察の対象となっている構成概念では説明出来ない結果側の変数の変動を生み出す要因としてモデルに導入された変数
  • 外生変数と内生変数
    • 外生変数 : モデルの中で1度も他の変数の結果にならない変数
    • 内生変数 : モデルの中で生成される変数

方程式モデルの可視化

パス図を使って、方程式モデルを可視化する。

  • 観測変数は四角で囲む
  • 潜在変数（構成概念）は楕円で囲む
  • 誤差変数は囲まない
  • 影響を与える変数から、与えられる変数に単方向の矢印を書き、矢印に因果の影響力を示す数値を付与
  • 共変動を示す2つの外生変数に因果関係を仮定しない時には双方向の矢印を書き、矢印に共分散を示す数値を付与

モデルの記述

構造方程式モデリングでは、手元のデータの性質や分析社の仮説に基づいてオリジナルの分析モデルをあつらえたり、従来型の多変量解析の多くのモデルを下位モデルとして実行することが出来る。
ここでは、3つのモデルをやってみる

2次因子分析モデル

下図のように、観測変数の背後に仮定した構成概念の背後に、さらに高次の構成概念を仮定したモデル。


MIMICモデル（Multiple Indicator MultIple Cause Model）

下図のように、複数の観測指標によって、構成概念が規定され、その構成概念が観測変数の原因となっているモデル。

縦断的モデル1

下図のように、時間経過に伴って構成概念が変化し、各時点で複数の測定がなされているモデル。
縦断的モデルでは、構成概念の事を「波（wave）」と呼び、例えば下図を3波2指標パネルモデルと呼ぶ。


実行例

以下、3つのモデルについて、R言語での実行例を書く。

準備

パッケージのインストール及び読み込み

> install.packages("sem")
> library(sem)
2次因子分析モデル
  • モデルが識別されるように、いくつかの母数を固定する
    • f1 : 外生変数なので、\sigma^2_{f1}, a1, a2, a3のうちどれかを1つを固定。
    • f2 : 内生変数なので、b1,b2のどちらか1つを固定。
    • f3 : 内生変数なので、b3,b4のどちらか1つを固定。
    • f4 : 内生変数なので、b5,b6のどちらか1つを固定。
  • データ
    • 以下のような6つのテストの相関行列
TEST1 1.00
TEST2 0.703 1.00
TEST3 0.432 0.469 1.00
TEST4 0.424 0.461 0.738 1.00
TEST5 0.330 0.303 0.471 0.427 1.00
TEST6 0.348 0.388 0.433 0.490 0.717 1.00
  • 構成概念の命名
    • f1を知能、f2を動作知能、f3を記憶力、f4を数理知能と命名
  • 実行

相関行列を入力

> TEST.cor <- read.moments(names = c("TEST1", "TEST2", "TEST3", "TEST4","TEST5", "TEST6"))
1:  1.00
2:  0.703 1.00
4:  0.432 0.469 1.00
7:  0.424 0.461 0.738 1.00
11: 0.330 0.303 0.471 0.427 1.00
16: 0.348 0.388 0.433 0.490 0.717 1.00
22: 
Read 21 items

> TEST.cor
      TEST1 TEST2 TEST3 TEST4 TEST5 TEST6
TEST1 1.000 0.000 0.000 0.000 0.000     0
TEST2 0.703 1.000 0.000 0.000 0.000     0
TEST3 0.432 0.469 1.000 0.000 0.000     0
TEST4 0.424 0.461 0.738 1.000 0.000     0
TEST5 0.330 0.303 0.471 0.427 1.000     0
TEST6 0.348 0.388 0.433 0.490 0.717     1

モデルを指定

> TEST.model <- specify.model()
知能 -> 動作知能, a1, NA
知能 -> 記憶力,   a2, NA
知能 -> 数理知能, a3, NA
動作知能 -> TEST1, NA, 1
動作知能 -> TEST2, b2, NA
記憶力   -> TEST3, NA, 1
記憶力   -> TEST4, b4, NA
数理知能 -> TEST5, NA, 1
数理知能 -> TEST6, b6, NA
知能    <-> 知能, NA, 1
動作知能 <-> 動作知能, d2, NA
記憶力   <-> 記憶力,  d3, NA
数理知能 <-> 数理知能, d4, NA
TEST1 <-> TEST1, e1, NA
TEST2 <-> TEST2, e2, NA
TEST3 <-> TEST3, e3, NA
TEST4 <-> TEST4, e4, NA
TEST5 <-> TEST5, e5, NA
TEST6 <-> TEST6, e6, NA

Read 19 records

モデル構築

> TEST.sem <- sem(ram = TEST.model, S = TEST.cor, N = 1000)
> summary(TEST.sem)

 Model Chisquare =  55.935   Df =  6 Pr(>Chisq) = 3.0007e-10
 Chisquare (null model) =  2918.6   Df =  15
 Goodness-of-fit index =  0.98236
 Adjusted goodness-of-fit index =  0.93827
 RMSEA index =  0.091273   90% CI: (0.0703, 0.11386)
 Bentler-Bonnett NFI =  0.98083
 Tucker-Lewis NNFI =  0.957
 Bentler CFI =  0.9828
 SRMR =  0.015852
 BIC =  14.488 

 Normalized Residuals
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1.25000 -0.11500  0.00001 -0.00386  0.24000  0.97900 

 Parameter Estimates
   Estimate Std Error z value Pr(>|z|)                         
a1 0.56105  0.034836  16.1055 0.0000e+00 動作知能 <--- 知能    
a2 0.76424  0.036971  20.6712 0.0000e+00 記憶力 <--- 知能      
a3 0.57208  0.034329  16.6645 0.0000e+00 数理知能 <--- 知能    
b2 1.07429  0.052826  20.3363 0.0000e+00 TEST2 <--- 動作知能   
b4 1.00951  0.040304  25.0478 0.0000e+00 TEST4 <--- 記憶力     
b6 1.06620  0.051607  20.6599 0.0000e+00 TEST6 <--- 数理知能   
d2 0.33961  0.033132  10.2501 0.0000e+00 動作知能 <--> 動作知能
d3 0.14698  0.038053   3.8625 1.1224e-04 記憶力 <--> 記憶力    
d4 0.34521  0.034049  10.1385 0.0000e+00 数理知能 <--> 数理知能
e1 0.34561  0.031352  11.0237 0.0000e+00 TEST1 <--> TEST1      
e2 0.24478  0.033320   7.3463 2.0384e-13 TEST2 <--> TEST2      
e3 0.26896  0.026163  10.2802 0.0000e+00 TEST3 <--> TEST3      
e4 0.25498  0.026277   9.7034 0.0000e+00 TEST4 <--> TEST4      
e5 0.32752  0.031346  10.4486 0.0000e+00 TEST5 <--> TEST5      
e6 0.23553  0.033210   7.0921 1.3205e-12 TEST6 <--> TEST6      

 Iterations =  50

因子得点の算出

> fscores(TEST.sem, data=df, centering=T, scale=F) #今回はデータがないため未実行

パス図の描画（graphvizで出力する）

> con <- file("TEST1.dot", encoding="UTF-8")
> capture.output(file = con,
 path.diagram(TEST.sem, ignore.double=F,
 edge.labels = "both", standardize = F,
 node.font = c("Osaka", 12),
 edge.font = c("Osaka", 10),
 rank.direction = "TB"))
> system("dot -Tpng TEST.dot -o TEST.png")


MIMICモデル
  • モデルが識別されるように、母数を固定する
    • f1: 内生変数なので、αb41,αb51,αb61のうちどれか一つを固定
  • データ
    • 以下のような相関行列
収入 1.00
学歴 0.434 1.00
職業威信 0.508 0.459 1.00
人脈 0.359 0.142 0.248 1.00
知名度 0.220 0.234 0.215 0.816 1.00
社会参加 0.259 0.299 0.220 0.868 0.779 1.00
  • 実行

相関行列の入力

> social.cor <- read.moments(names = c("収入", "学歴", "職業威信", "人脈","知名度", "社会参加"))
1.00
0.434 1.00
0.508 0.459 1.00
0.359  0.142 0.248 1.00
0.220 0.234 0.215 0.816 1.00 
0.259 0.299 0.220 0.868 0.779 1.00
22: 
Read 21 items

モデルの指定

> social.model <- specify.model()
収入      -> 社会的地位, αd1
学歴      -> 社会的地位,αd2
職業威信   -> 社会的地位, αd3
社会的地位 -> 人脈, NA, 1
社会的地位 -> 知名度,αb51, NA
社会的地位 -> 社会参加,αb61, NA
社会的地位 <-> 社会的地位,d1, NA
収入 <-> 学歴,   σ21, NA
収入 <-> 職業威信,σ31, NA
学歴 <-> 職業威信,σ32, NA
人脈    <-> 人脈,    e4,  NA
知名度  <-> 知名度,   e5, NA
社会参加 <-> 社会参加, e6, NA
収入 <-> 収入, NA, 1
学歴 <-> 学歴, NA, 1
職業威信 <-> 職業威信, NA, 1
17: 
Read 16 records

モデル構築

> social.sem <- sem(ram = social.model, S = social.cor, N = 1000)
> summary(social.sem)

 Model Chisquare =  324.19   Df =  9 Pr(>Chisq) = 0
 Chisquare (null model) =  3610.8   Df =  15
 Goodness-of-fit index =  0.9131
 Adjusted goodness-of-fit index =  0.79722
 RMSEA index =  0.18723   90% CI: (0.17006, 0.20498)
 Bentler-Bonnett NFI =  0.91022
 Tucker-Lewis NNFI =  0.8539
 Bentler CFI =  0.91234
 SRMR =  0.035230
 BIC =  262.02 

 Normalized Residuals
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-2.0100 -0.0436  0.0000  0.0413  0.1460  3.3300 

 Parameter Estimates
     Estimate Std Error z value Pr(>|z|)                             
αd1  0.253322 0.036540   6.9327 4.1296e-12 社会的地位 <--- 収入      
αd2  0.050914 0.036033   1.4130 1.5767e-01 社会的地位 <--- 学歴      
αd3  0.091187 0.035735   2.5518 1.0718e-02 社会的地位 <--- 職業威信  
αb51 0.889684 0.021464  41.4504 0.0000e+00 知名度 <--- 社会的地位    
αb61 0.946008 0.019856  47.6432 0.0000e+00 社会参加 <--- 社会的地位  
d1   0.804465 0.040311  19.9563 0.0000e+00 社会的地位 <--> 社会的地位
σ21  0.434000 0.023368  18.5726 0.0000e+00 学歴 <--> 収入            
σ31  0.508000 0.020846  24.3696 0.0000e+00 職業威信 <--> 収入        
σ32  0.459000 0.022561  20.3446 0.0000e+00 職業威信 <--> 学歴        
e4   0.081529 0.011974   6.8089 9.8324e-12 人脈 <--> 人脈            
e5   0.272995 0.015040  18.1518 0.0000e+00 知名度 <--> 知名度        
e6   0.178031 0.012924  13.7755 0.0000e+00 社会参加 <--> 社会参加    

 Iterations =  35

因子得点の算出

> fscores(social.sem, data=df, centering=T, scale=F) #今回はデータがないため未実行

パス図の描画（graphvizで出力する）

> con <- file("social.dot", encoding="UTF-8")
> capture.output(file = con,
 path.diagram(social.sem, ignore.double=F,
 edge.labels = "both", standardize = F,
 node.font = c("Osaka", 12),
 edge.font = c("Osaka", 10),
 rank.direction = "LR"))
> system("dot -Tpng social.dot -o social.png")


※ 見やすいように、dotファイルをちょっと修正してします。

縦断的モデル
  • モデルが識別されるように、いくつかの母数を固定する
    • f1 : 外生変数なので、\sigma^2_{f1}, a1, a2, a3のうちどれかを1つを固定。
    • f2 : 内生変数なので、αb32,αb42のどちらか1つを固定。
    • f3 : 内生変数なので、αb53,αb63のどちらか1つを固定。
  • データ
    • 以下のような相関行列
検査1 1.00
面接1 0.610 1.00
検査2 0.422 0.297 1.00
面接2 0.268 0.419 0.660 1.00
検査3 0.466 0.223 0.400 0.228 1.00
面接3 0.221 0.458 0.241 0.412 0.585 1.00
  • 構成概念の命名
    • 初期自己像、中期自己像、終期自己像と命名
  • 実行

相関行列を入力

> kensa.cor <- read.moments(names = c("検査1", "面接1", "検査2", "面接2","検査3", "面接3"))
1.00
0.610 1.00
0.422 0.297 1.00
0.268  0.419 0.660 1.00
0.466 0.223 0.400 0.228 1.00
0.221 0.458 0.241 0.412 0.585 1.00
22: 
Read 21 items

モデルを指定

> kensa.model <- specify.model()
初期自己像 -> 検査1, NA, 1
初期自己像 -> 面接1, αb21, NA
中期自己像 -> 検査2, NA, 1
中期自己像 -> 面接2, αb41, NA
終期自己像 -> 検査3, NA, 1
終期自己像 -> 面接3, αb61, NA
初期自己像 -> 中期自己像, αa21, NA
中期自己像 -> 終期自己像, αa32, NA
初期自己像 <-> 初期自己像, NA, 1
中期自己像 <-> 中期自己像, d2, NA
終期自己像 <-> 終期自己像, d3, NA
検査1 <-> 検査1, e1, NA
面接1 <-> 面接1, e2, NA
検査2 <-> 検査2, e3, NA
面接2 <-> 面接2, e4, NA
検査3 <-> 検査3, e5, NA
面接3 <-> 面接3, e6, NA
18: 
Read 17 records

モデル構築

> kensa.sem <- sem(ram = kensa.model, S = kensa.cor, N = 1000)
> summary(kensa.sem)

 Model Chisquare =  788.47   Df =  8 Pr(>Chisq) = 0
 Chisquare (null model) =  2627.1   Df =  15
 Goodness-of-fit index =  0.7859
 Adjusted goodness-of-fit index =  0.43798
 RMSEA index =  0.3125   90% CI: (0.29426, 0.33114)
 Bentler-Bonnett NFI =  0.69987
 Tucker-Lewis NNFI =  0.43976
 Bentler CFI =  0.7012
 SRMR =  0.11342
 BIC =  733.21 

 Normalized Residuals
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.330  -1.620  -0.209   0.766   1.820  10.200 

 Parameter Estimates
     Estimate Std Error z value Pr(>|z|)                             
αb21 0.68102  0.044440  15.3243 0.0000e+00 面接1 <--- 初期自己像     
αb41 0.80549  0.050278  16.0209 0.0000e+00 面接2 <--- 中期自己像     
αb61 0.76603  0.078965   9.7009 0.0000e+00 面接3 <--- 終期自己像     
αa21 0.48021  0.038657  12.4223 0.0000e+00 中期自己像 <--- 初期自己像
αa32 0.49902  0.041492  12.0269 0.0000e+00 終期自己像 <--- 中期自己像
d2   0.59471  0.061660   9.6449 0.0000e+00 中期自己像 <--> 中期自己像
d3   0.56325  0.081439   6.9162 4.6381e-12 終期自己像 <--> 終期自己像
e1   0.10789  0.052846   2.0416 4.1189e-02 検査1 <--> 検査1          
e2   0.57738  0.038787  14.8860 0.0000e+00 面接1 <--> 面接1          
e3   0.19515  0.047282   4.1274 3.6685e-05 検査2 <--> 検査2          
e4   0.47780  0.034629  13.7979 0.0000e+00 面接2 <--> 面接2          
e5   0.23632  0.075062   3.1484 1.6418e-03 検査3 <--> 検査3          
e6   0.55187  0.050110  11.0133 0.0000e+00 面接3 <--> 面接3          

 Iterations =  33

因子得点の算出

> fscores(kensa.sem, data=df, centering=T, scale=F) #今回はデータがないため未実行

パス図の描画（graphvizで出力する）

> con <- file("kensa.dot", encoding="UTF-8")
> capture.output(file = con,
 path.diagram(kensa.sem, ignore.double=F,
 edge.labels = "both", standardize = F,
 node.font = c("Osaka", 12),
 edge.font = c("Osaka", 10),
 rank.direction = "TB"))
> system("dot -Tpng kensa.dot -o kensa.png")



※ 見やすいように、dotファイルをちょっと修正してします。

参考