Hi,
I'm using OpenMx to run a five-factor CFA model for 14 variables. "phenfile0" is the dataset and I used its 3rd to 16th columns as a input with the colnames of "y1" to "y14".
"y12" "y13" "y14" belong to "F1"; "y9" "y10" "y11" belong to "F2"; "y6" "y7" "y8" belong to F3; "y3" "y4" "y5" belong to F4; and "y1" "y2" belong to F5.
Issue: the output gave a larger SE for some loadings. I'm not sure if the code is correct. Thanks for your help!
phenofacmodel <- mxModel( name="PhenFactor", # data mxData(observed=phenfile0[,3:16],type="raw"), mxMatrix(type="Full",nrow=1,ncol=14,free=T,values=colMeans(phenfile0[,3:16],na.rm=T), labels=c("m1","m2","m3","m4","m5","m6","m7","m8","m9","m10","m11","m12","m13","m14"),name="Mu", dimnames=list(NULL,c('y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8', 'y9', 'y10', 'y11', 'y12', 'y13', 'y14'))), mxMatrix(type="Full",nrow=14,ncol=1,name="Lambda1", free= c(F,F,F,F,F,F,F,F,F,F,F,T,T,T), values=c(0,0,0,0,0,0,0,0,0,0,0,1,1,1), labels= c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,"l1_12","l1_13","l1_14"), dimnames=list( c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,'y12','y13','y14'),"F1")), mxMatrix(type="Full",nrow=14,ncol=1,name="Lambda2", free= c(F,F,F,F,F,F,F,F,T,T,T,F,F,F), values=c(0,0,0,0,0,0,0,0,1,1,1,0,0,0), labels= c(NA,NA,NA,NA,NA,NA,NA,NA,"l2_9","l2_10","l2_11",NA,NA,NA), dimnames=list( c(NA,NA,NA,NA,NA,NA,NA,NA,'y9','y10','y11',NA,NA,NA),"F2")), mxMatrix(type="Full",nrow=14,ncol=1,name="Lambda3", free= c(F,F,F,F,F,T,T,T,F,F,F,F,F,F), values=c(0,0,0,0,0,1,1,1,0,0,0,0,0,0), labels= c(NA,NA,NA,NA,NA,"l3_6","l3_7","l3_8",NA,NA,NA,NA,NA,NA), dimnames=list( c(NA,NA,NA,NA,NA,'y6','y7','y8',NA,NA,NA,NA,NA,NA),"F3")), mxMatrix(type="Full",nrow=14,ncol=1,name="Lambda4", free= c(F,F,T,T,T,F,F,F,F,F,F,F,F,F), values=c(0,0,1,1,1,0,0,0,0,0,0,0,0,0), labels= c(NA,NA,"l4_3","l4_4","l4_5",NA,NA,NA,NA,NA,NA,NA,NA,NA), dimnames=list( c(NA,NA,'y3','y4','y5',NA,NA,NA,NA,NA,NA,NA,NA,NA),"F4")), mxMatrix(type="Full",nrow=14,ncol=1,name="Lambda5", free= c(T,T,F,F,F,F,F,F,F,F,F,F,F,F), values=c(1,1,0,0,0,0,0,0,0,0,0,0,0,0), labels= c("l5_1","l5_2",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA), dimnames=list( c('y1','y2',NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),"F5")), #The common factor's variance, fixed to 1.0 to identify the scale: mxMatrix(type="Symm",nrow=1,free=F,values=1,labels="fvar1",lbound=0,name="Fvar1"), mxMatrix(type="Symm",nrow=1,free=F,values=1,labels="fvar2",lbound=0,name="Fvar2"), mxMatrix(type="Symm",nrow=1,free=F,values=1,labels="fvar3",lbound=0,name="Fvar3"), mxMatrix(type="Symm",nrow=1,free=F,values=1,labels="fvar4",lbound=0,name="Fvar4"), mxMatrix(type="Symm",nrow=1,free=F,values=1,labels="fvar5",lbound=0,name="Fvar5"), #The unique variances: mxMatrix( type="Diag",nrow=14,free=T, labels=c("u1","u2","u3","u4","u5","u6","u7","u8","u9","u10","u11","u12","u13","u14"),lbound=0.0001,name="U", values=diag(cov(phenfile0[,3:16],use="pair"))), #Compute mxAlgebra(Lambda1%*%Fvar1%*%t(Lambda1),name="Comm1"), mxAlgebra(Lambda2%*%Fvar2%*%t(Lambda2),name="Comm2"), mxAlgebra(Lambda3%*%Fvar3%*%t(Lambda3),name="Comm3"), mxAlgebra(Lambda4%*%Fvar4%*%t(Lambda4),name="Comm4"), mxAlgebra(Lambda5%*%Fvar5%*%t(Lambda5),name="Comm5"), mxAlgebra(Comm1+Comm2+Comm3+Comm4+Comm5+U,name="Sigmay",dimnames=list( c('y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8', 'y9', 'y10', 'y11', 'y12', 'y13', 'y14'), c('y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8', 'y9', 'y10', 'y11', 'y12', 'y13', 'y14'))), mxExpectationNormal(covariance="Sigmay",means="Mu"), mxFitFunctionML() ) phenofacrun <- mxRun(phenofacmodel) summary(phenofacrun, refModels=mxRefModels(phenofacrun, run = TRUE))
Hi
It's quite possible for factor loadings to have different standard errors. This will occur if the variables have different variances, or different response probabilities if the data are binary or ordinal. BTW are they continuous, normally distributed variables or something else?
To simplify model specification for those who get path diagrams, the Onyx-sem software is really handy - it will spit out OpenMx code from a diagram.
I note that there are also some possible simplifications of your code, where the lambda vectors are
cbind()
(column-wise glued together) into a matrix Lambda (or just set up that way initially), and the variances specified as a diagonal matrix FCov- actually an identity matrix since all the variances are 1 and no covariances between factors are specified. Then the covariance algebraLambda %*% FCov %*% t(Lambda)
stays the same however many factors you have vary the number of columns, and you also have the option of letting them correlate.I agree with Mike Neale about the standard errors.
Below is a version of the same model that is just easier (in my opinion) to read and write.