《结构方程模型及其应用》(侯, 温, & 成,2004)部分章节R代码

John FOX教授的sem包试写了这本教材的几个例子,结果都与LISREL8报告的Minimum Fit Function Chi-Square吻合。不过LISREL其它拟合指标用的是Normal Theory Weighted Least Squares Chi-Square?,所以看上去比Minimum Fit Function Chi-Square报告的结果要好那么一些些。LISREL历史上推的GFI/AGFI曾因为经常将差报好被批评,圈内朋友私下嘲笑这样LISREL就更好卖了,买软件的用户也高兴将差报好,只有读Paper的人上当。

目前只写了chap3_1..到Chap3_2_...,一共5(或6)个例子。尺有所短,寸有所长--sem包不会自动报告所有修正指数,不能做样本量不一样多的多组模型,对复杂的模型要写的代码太多。不过,在已经尝试的几个例子里,有一个是LISREL跑不出来但sem包能跑出结果的。目前sem包还没有达到结构方程众多商业软件的成熟水准,但R庞大的义工武器库已经使sem包至少已经在Missing Data Multiple Imputation、Bootstrapping等等应用上胜出一筹。所以例子中还附带了缺失数据一讲Multiple Impuation的R示范代码。

欢迎各位有兴趣作类似尝试的同学将结果email我, 可以陆续更新到下面的GPL版权代码集合中。

代码下载:lixiaoxu.googlepages.com (中大镜像)

[update] AMos, Mplus 与 SAS/STAT CALIS 缺省报告与使用的都是 Minimum Fit Function Chi-Square,可通过各种软件(Albright & Park, 2009)的结果对比查验。Normal Theory Weighted Least Squares Chi-Square并非总是比Minimum Fit Function Chi-Square报告更“好”的拟合结果(虽然常常如此)。Olsson, Foss, 和 Breivik (2004) 用模拟数据对比了二者,确证Minimum Fit Function Chi-Square计算得到的拟合指标在小样本之外的情形都比Normal Theory Weighted Least Squares Chi-Square的指标更适用。

R的sem包目前在迭代初值的计算上还做得不够好。我遇到的几个缺省初值下迭代不收敛案例,将参数设定合理范围的任意初值(比如TX,TY都设成0-1之间的正实数)之后都收敛了。

--

Albright, J. J., & Park., H. M. (2009). Confirmatory factor analysis using Amos, LISREL, Mplus, and SAS/STAT CALIS. The University Information Technology Services Center for Statistical and Mathematical Computing, Indiana University. Retrieved July 7, 2009, from http://www.indiana.edu/~statmath/stat/all/cfa/cfa.pdf.

Olsson, U. H., Foss, T., & Breivik, E. (2004). Two equivalent discrepancy functions for maximum likelihood estimation: Do their test statistics follow a non-central chi-square distribution under model misspecification? Sociological Methods Research, 32(4), 453-500.

[update]R中的sem包妙处之一是可在线实现结构方程应用界面。下面这个最粗陋的例子对应原书Chap3_1_4_CFA_MB.LS8的结果。
--荣耀属于sem package的作者Rweb的作者、以及服务器运算资源的提供者。

##Input Correlation Matrix
R.DHP<-matrix(0,ncol=17,nrow=17);
R.DHP[col(R.DHP) >= row(R.DHP)] <- c(
1,
.34,1,
.38,.35,1,
.02,.03,.04,1,
.15,.19,.14,.02,1,
.17,.15,.20,.01,.42,1,
.20,.13,.12,.00,.40,.21,1,
.32,.32,.21,.03,.10,.10,.07,1,
.10,.17,.12,.02,.15,.18,.23,.13,1,
.14,.16,.15,.03,.14,.19,.18,.18,.37,1,
.14,.15,.19,.01,.18,.30,.13,.08,.38,.38,1,
.18,.16,.24,.02,.14,.21,.21,.22,.06,.23,.18,1,
.19,.20,.15,.01,.14,.24,.09,.24,.15,.21,.21,.45,1,
.18,.21,.18,.03,.25,.18,.18,.18,.22,.12,.24,.28,.35,1,
.08,.18,.16,.01,.22,.20,.22,.12,.12,.16,.21,.25,.20,.26,1,
.12,.16,.25,.02,.15,.12,.20,.14,.17,.20,.14,.20,.15,.20,.50,1,
.20,.16,.18,.04,.25,.14,.21,.17,.21,.21,.23,.15,.21,.22,.29,.41,1
);
R.DHP<-t(R.DHP);
colnames(R.DHP)<-rownames(R.DHP)<-paste('X',1:17,sep='');
print('Inputted Correlation Matrix');
print(R.DHP);
##
##
##Input Model Specification of Chap3_1_4_CFA_MB.LS8
require(sem);
model.B <- matrix(ncol=3,byrow=TRUE,data=c(
'X1 <-> X1' , 'TD1_1'  , NA  ,
'X2 <-> X2' , 'TD2_2'  , NA  ,
'X3 <-> X3' , 'TD3_3'  , NA  ,
'X5 <-> X5' , 'TD5_5'  , NA  ,
'X6 <-> X6' , 'TD6_6'  , NA  ,
'X7 <-> X7' , 'TD7_7'  , NA  ,
'X8 <-> X8' , 'TD8_8'  , NA  ,
'X9 <-> X9' , 'TD9_9'  , NA  ,
'X10<-> X10', 'TD10_10', NA  ,
'X11<-> X11', 'TD11_11', NA  ,
'X12<-> X12', 'TD12_12', NA  ,
'X13<-> X13', 'TD13_13', NA  ,
'X14<-> X14', 'TD14_14', NA  ,
'X15<-> X15', 'TD15_15', NA  ,
'X16<-> X16', 'TD16_16', NA  ,
'X17<-> X17', 'TD17_17', NA  ,
'xi1<-> xi1', NA       , '1' ,
'xi2<-> xi2', NA       , '1' ,
'xi3<-> xi3', NA       , '1' ,
'xi4<-> xi4', NA       , '1' ,
'xi5<-> xi5', NA       , '1' ,
'xi1<-> xi2', 'PH12'   , NA  ,
'xi1<-> xi3', 'PH13'   , NA  ,
'xi1<-> xi4', 'PH14'   , NA  ,
'xi1<-> xi5', 'PH15'   , NA  ,
'xi2<-> xi3', 'PH23'   , NA  ,
'xi2<-> xi4', 'PH24'   , NA  ,
'xi2<-> xi5', 'PH25'   , NA  ,
'xi3<-> xi4', 'PH34'   , NA  ,
'xi3<-> xi5', 'PH35'   , NA  ,
'xi4<-> xi5', 'PH45'   , NA  ,
'X1 <- xi1' , 'LX1_1'  , NA  ,
'X2 <- xi1' , 'LX2_1'  , NA  ,
'X3 <- xi1' , 'LX3_1'  , NA  ,
'X5 <- xi2' , 'LX5_2'  , NA  ,
'X6 <- xi2' , 'LX6_2'  , NA  ,
'X7 <- xi2' , 'LX7_2'  , NA  ,
'X8 <- xi1' , 'LX8_1'  , NA  ,
'X9 <- xi3' , 'LX9_3'  , NA  ,
'X10<- xi3' , 'LX10_3' , NA  ,
'X11<- xi3' , 'LX11_3' , NA  ,
'X12<- xi4' , 'LX12_4' , NA  ,
'X13<- xi4' , 'LX13_4' , NA  ,
'X14<- xi4' , 'LX14_4' , NA  ,
'X15<- xi5' , 'LX15_5' , NA  ,
'X16<- xi5' , 'LX16_5' , NA  ,
'X17<- xi5' , 'LX17_5' , NA  )
);
class(model.B)<-'mod';
##
##
N=350;##sample size;
##R.DHP[-4,-4] excludes X4
## Result
(summary(sem.B<-sem(model.B, R.DHP[-4,-4], N)));
## Residuals
(round(residuals(sem.B),3));
#####################
boxplot.matrix = function(M,ylim=c(-1,1)) {
  M = as.matrix(M);
  boxplot(c(M[row(M)>col(M)]),at=1,xlab='',ylab='',ylim=ylim);
  points(rep(1,length(c(M[row(M)>col(M)]))),c(M[row(M)>col(M)]),pch='-',col='red');
  stem(c(M[row(M)>col(M)]));
  boxplot.stats(c(M[row(M)>col(M)]));
}
####################
boxplot(residuals(sem.B));

Automatize LISREL jobs

LISREL routine can run in DOS or in command line mode of windows (windows-key + R -> CMD) . The command line is just like --

D:\My Documents>"C:\Program Files\lisrel87\lisrel87.exe" "C:\Program Files\lisrel87\LS8EX\EX61.LS8" D:\myOutput.out

1. You only need edit and input the bold part.
2. Quotation marks are used wherever the paths or filenames include blanks.
3. The 2nd argument is the output file. You can still specify more output options in your .ls8 file.
4. A .bat file can automatize batches of such lisrel jobs.