構成概念間のグラフィカルモデリング


2007年の心理学研究に掲載されている論文「[構成概念間のグラフィカルモデリング（豊田・福中・尾崎・川端, 2007）](http://t.umblr.com/redirect?z=https%3A%2F%2Fwww.jstage.jst.go.jp%2Farticle%2Fjjpsy1926%2F78%2F4%2F78_4_416%2F_pdf&t=MGYyNzFjMTdiMjljZmQ1ODgyZmM5NDNkNzk4NThkMTJhODc3M2VmMyxnTE4yTjVlMg%3D%3D)」に使用されていた分析例のデータとスクリプト（Mx用）を見つけたので（[この論文](http://t.umblr.com/redirect?z=https%3A%2F%2Fdspace.wul.waseda.ac.jp%2Fdspace%2Fbitstream%2F2065%2F37670%2F3%2FHonbun-5970.pdf&t=MjNhMGJmZGFiNzQ3Y2RmNTJhMjMzNmQ5ZTdjNjgxMWQ0MDBiMjI5MCxnTE4yTjVlMg%3D%3D)の最後の方），とりあえずこれをRのOpenMx 2.0で動かしてみました。

データはharman.csvという名前で別ファイルに保存してあるものとします。

### 下準備
まず，構造方程式モデリングのための準備作業を行います。

```r:harman_prep.R
require(OpenMx)

#### データの読み込み ####  
cmat<-read.csv("Harman.csv")
varname=paste("x",c(1:17),sep="")
rownames(cmat)<-colnames(cmat)<-varname

cmat<-as.matrix(cmat)

#### 行列A：潜在変数から観測変数への係数 ####
A<-  mxMatrix("Full",17,4,name="A",
              free=  c(F,F,F,F,  T,F,F,F,  T,F,F,F,  T,F,F,F,  
                       F,F,F,F,  F,T,F,F,  F,T,F,F,  F,T,F,F,  F,T,F,F,  
                       F,F,F,F,  F,F,T,F,  F,F,T,F,  F,F,T,F,  
                       F,F,F,F,  F,F,F,T,  F,F,F,T,  F,F,F,T), byrow=T)

#### 各潜在変数の最初の観測変数の係数を1に固定 ####
A@values[1,1]  <- 1.0; A@values[2,1]  <- 0.5
A@values[3,1]  <- 0.5; A@values[4,1]  <- 0.5
A@values[5,2]  <- 1.0; A@values[6,2]  <- 0.5
A@values[7,2]  <- 0.5; A@values[8,2]  <- 0.5
A@values[9,2]  <- 0.5
A@values[10,3] <- 1.0; A@values[11,3] <- 0.5
A@values[12,3] <- 0.5; A@values[13,3] <- 0.5
A@values[14,4] <- 1.0; A@values[15,4] <- 0.5
A@values[16,4] <- 0.5; A@values[17,4] <- 0.5

#### 行列E：誤差分散 ####
E<-  mxMatrix("Diag",17,17,free=T,name="E", values=0.5)

#### 行列R：共分散行列の逆行列 ####
R<-  mxMatrix("Symm",4,4,name="R",
              values=c(  3.051, 
                        -0.950,  2.511, 
                        -0.892, -0.463,  3.810, 
                        -0.908, -0.751, -1.239, 4.470),
              byrow=T, free=T)

#### 行列I：（偏）相関係数処理用の単位行列 ####
I<- mxMatrix("Iden",4,4,name="I")

#### 行列COV：分析モデル ####
COV<-mxAlgebra(A %&% solve(R)+ E * E , name="COV")

#### 行列S：相関行列 ####
S<-mxAlgebra(solve(sqrt(I*solve(R))) %*% solve(R) %*% 
       solve(sqrt(I*solve(R))), name="S")

#### 行列P：偏相関行列 ####
P<-mxAlgebra(solve(sqrt(I*solve(S))) %*% solve(S) %*%
       solve(sqrt(I*solve(S))), name="P")

#### 分析モデルを指定 ####
Func<-mxExpectationNormal("COV",dimnames=varname)
fit<-mxFitFunctionML()

#### 観測データの指定 ####
data <- mxData(cmat, 'cor', numObs = 145)
```

### 分析手順
まず，偏相関を固定しないモデルを作成して実行します。

```r
model.0<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.0<-mxRun(model.0)
```

相関行列をデータとする場合には適合度がうまく計算出来ないとの警告メッセージが出るが，とりあえず結果は求められます。

#### 警告メッセージ

```r
In FUN(X[[1L]], ...) :
  OpenMx does not yet correctly handle mxData(type='cor') standard errors and fit statistics. See Steiger (1980), "Tests for comparing elements of a correlation matrix".
```

結果の要約を出力します。

```r
summary(Res.0)
```

以下，モデル適合度に関する部分のみ掲載すると次のようになります。

```r
observed statistics:  136 
estimated parameters:  40 
degrees of freedom:  96 
-2 log likelihood:  1590.964 
saturated -2 log likelihood:  1388.592 
number of observations:  145 
chi-square:  X2 ( df=113 ) = 202.3716,  p = 4.986103e-07
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       10.37158               282.3716                       NA
BIC:     -275.39486               401.4409                 274.8667
CFI: 0.9032155 
TLI: 0.883516 
RMSEA:  0.07385442  [95% CI (0.05366435, 0.09308408)]
```

偏相関係数を表示させます。

```r
Res.0$P
```

ただし，この行列の偏相関係数は符号を逆転させる前のものである点に注意が必要です。

偏相関行列の中で，絶対値が最も小さい箇所を探します。

因子2と3の偏相関（0.15）が最小なので，共分散行列の逆行列`R`の因子2と3の箇所を0と置いて固定します。

```r
R@values[3,2]<-R@values[2,3]<-0
R@free[3,2]<-R@free[2,3]<-F
```

修正モデルを作成して実行します。

```r
model.32<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.32<-mxRun(model.32)
summary(Res.32)
```

結果の適合度以降の部分は次のようになります。

```r
observed statistics:  136 
estimated parameters:  39 
degrees of freedom:  97 
-2 log likelihood:  1592.55 
saturated -2 log likelihood:  1388.592 
number of observations:  145 
chi-square:  X2 ( df=114 ) = 203.9575,  p = 4.69223e-07
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       9.957483               281.9575                       NA
BIC:    -278.785690               398.0501                 274.6403
CFI: 0.902581 
TLI: 0.8837808 
RMSEA:  0.07377041  [95% CI (0.05366051, 0.09292332)]
```

`AIC.0=282.3716 > AIC.32=281.9575`と，モデルの適合度が上がっているので，さらに偏相関を固定してモデルを修正します。

このモデルの偏相関行列は以下のとおりです。

```r
Res.32$P

$result:
           [,1]          [,2]          [,3]       [,4]
[1,]  1.0000000 -3.893661e-01 -3.274052e-01 -0.1938470
[2,] -0.3893661  1.000000e+00  3.400425e-17 -0.2820814
[3,] -0.3274052 -1.020127e-16  1.000000e+00 -0.3401498
[4,] -0.1938470 -2.820814e-01 -3.401498e-01  1.0000000
dimnames: NULL
```

絶対値が最小の偏相関係数は因子１と４の間なので，行列`R`で因子１と４の部分を0に固定します。

```r
R@values[4,1]<-R@values[1,4]<-0
R@free[4,1]<-R@free[1,4]<-F

model.41<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.41<-mxRun(model.41)
summary(Res.41)
```

結果は次のとおりです（関連部分のみ記載）。

```r
observed statistics:  136 
estimated parameters:  38 
degrees of freedom:  98 
-2 log likelihood:  1594.207 
saturated -2 log likelihood:  1388.592 
number of observations:  145 
chi-square:  X2 ( df=115 ) = 205.6151,  p = 4.343026e-07
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       9.615086               281.6151                       NA
BIC:    -282.104820               394.7310                 274.4855
CFI: 0.9018688 
TLI: 0.8839492 
RMSEA:  0.07371695  [95% CI (0.0536906, 0.09279173)]
```

`AIC.32=281.9575 > AIC.41=281.6151`と，モデルの適合度が上がっているので，再度偏相関を固定してモデルを修正します。

```r
Res.41$P

$result:
              [,1]       [,2]          [,3]          [,4]
[1,]  1.000000e+00 -0.4231979 -3.799717e-01 -1.339203e-16
[2,] -4.231979e-01  1.0000000 -8.002529e-17 -3.403886e-01
[3,] -3.799717e-01  0.0000000  1.000000e+00 -3.970419e-01
[4,] -1.281201e-16 -0.3403886 -3.970419e-01  1.000000e+00
dimnames: NULL
```

偏相関係数の絶対値が最小であるのは因子２と４の間なので，そこを０に固定します。

```r
R@values[4,2]<-R@values[2,4]<-0
R@free[4,2]<-R@free[2,4]<-F

model.42<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.42<-mxRun(model.41)
summary(Res.42)
```

得られた結果（適合度関連部分のみ）は以下のとおりです。

```r
observed statistics:  136 
estimated parameters:  37 
degrees of freedom:  99 
-2 log likelihood:  1606.951 
saturated -2 log likelihood:  1388.592 
number of observations:  145 
chi-square:  X2 ( df=116 ) = 218.3588,  p = 2.900845e-08
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       20.35884               292.3588                       NA
BIC:     -274.33780               402.4980                 285.4169
CFI: 0.889151 
TLI: 0.8700391 
RMSEA:  0.07800989  [95% CI (0.05864812, 0.09668359)]
```

`AIC.41=281.6151 < AIC.42= 292.3588`で，修正後のモデル適合度が修正前より下がってしまったので，このモデルは破棄し，修正前のモデルを最終モデルとして採択します。

この分析で得られた最終的なモデルの因子間相関は，

```r
Res.41$S

$result:
          [,1]     [,2]      [,3]      [,4]
[1,] 1.0000000 0.572729 0.5493207 0.4206697
[2,] 0.5727290 1.000000 0.4218840 0.5167430
[3,] 0.5493207 0.421884 1.0000000 0.5511866
[4,] 0.4206697 0.516743 0.5511866 1.0000000
```

この間のRMSEAの変化は

```r
model.0  = 0.07385442
model.32 = 0.07377041
model.41 = 0.07371695
model.42 = 0.07800989
```

となり，原論文とほぼ同じ結果を得ることができました。