2007年の心理学研究に掲載されている論文「構成概念間のグラフィカルモデリング（豊田・福中・尾崎・川端, 2007）」に使用されていた分析例のデータとスクリプト（Mx用）を見つけたので（この論文の最後の方），とりあえずこれをRのOpenMx 2.0で動かしてみました。
データはharman.csvという名前で別ファイルに保存してあるものとします。
下準備
まず，構造方程式モデリングのための準備作業を行います。
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)
分析手順
まず，偏相関を固定しないモデルを作成して実行します。
model.0<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.0<-mxRun(model.0)
相関行列をデータとする場合には適合度がうまく計算出来ないとの警告メッセージが出るが，とりあえず結果は求められます。
警告メッセージ
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".
結果の要約を出力します。
summary(Res.0)
以下，モデル適合度に関する部分のみ掲載すると次のようになります。
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)]
偏相関係数を表示させます。
Res.0$P
ただし，この行列の偏相関係数は符号を逆転させる前のものである点に注意が必要です。
偏相関行列の中で，絶対値が最も小さい箇所を探します。
因子2と3の偏相関（0.15）が最小なので，共分散行列の逆行列Rの因子2と3の箇所を0と置いて固定します。
R@values[3,2]<-R@values[2,3]<-0
R@free[3,2]<-R@free[2,3]<-F
修正モデルを作成して実行します。
model.32<-mxModel(A,E,R,I,COV,S,P,fit,Func,data)
Res.32<-mxRun(model.32)
summary(Res.32)
結果の適合度以降の部分は次のようになります。
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と，モデルの適合度が上がっているので，さらに偏相関を固定してモデルを修正します。
このモデルの偏相関行列は以下のとおりです。
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@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)
結果は次のとおりです（関連部分のみ記載）。
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と，モデルの適合度が上がっているので，再度偏相関を固定してモデルを修正します。
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@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)
得られた結果（適合度関連部分のみ）は以下のとおりです。
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で，修正後のモデル適合度が修正前より下がってしまったので，このモデルは破棄し，修正前のモデルを最終モデルとして採択します。
この分析で得られた最終的なモデルの因子間相関は，
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の変化は
model.0 = 0.07385442
model.32 = 0.07377041
model.41 = 0.07371695
model.42 = 0.07800989
となり，原論文とほぼ同じ結果を得ることができました。