《结构方程模型及其应用》(侯, 温, & 成,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));

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>