You are here

The OpenMx website will be down for maintenance from 9 AM EDT on Tuesday, September 17th, and is expected to return by the end of the day on Wednesday, September 18th. During this period, the backend will be updated and the website will get a refreshed look.

five-factor model

3 posts / 0 new
Last post
loveyu3317's picture
Offline
Joined: 07/25/2024 - 14:38
five-factor model

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))
AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Data?

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 algebra Lambda %*% 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.

AdminHunter's picture
Offline
Joined: 03/01/2013 - 11:03
Make your life easier

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.

require(OpenMx)
nv <- 14
mv <- paste0('y', 1:nv)
 
# Make fake data for demonstration
phenfile0 <- as.data.frame(matrix(rnorm((nv+2)*100), nrow=100, ncol=nv+2))
names(phenfile0) <- c(letters[1:2], mv)
 
phenofacmodel <- mxModel(
  name="PhenFactor",
  manifestVars=list(exo=mv),
  latentVars=list(exo=paste0('F', 1:5)),
  type='LISREL',
 
  # Data
  mxData(observed=phenfile0[,3:16],type="raw"),
 
  # Means
  mxPath('one', mv, labels=paste0('m', 1:nv), values=colMeans(phenfile0[,3:16],na.rm=T)),
 
  # Loadings
  mxPath('F1', mv[12:14], labels=paste0('l1_', 12:14), values=1),
  mxPath('F2', mv[9:11], labels=paste0('l2_', 9:11), values=1),
  mxPath('F3', mv[6:8], labels=paste0('l3_', 6:8), values=1),
  mxPath('F4', mv[3:5], labels=paste0('l4_', 3:5), values=1),
  mxPath('F5', mv[1:2], labels=paste0('l5_', 1:2), values=1),
 
  # Factor Variances
  mxPath(paste0('F', 1:5), arrows=2, free=FALSE, values=1),
 
  # Unique variances
  mxPath(mv, arrows=2, labels=paste0('u', 1:14), lbound=0.0001,
         values=diag(cov(phenfile0[,3:16],use="pair")))
)
 
# Loadings
mxEval(LX, phenofacmodel)
 
# Resid
mxEval(TD, phenofacmodel)
 
# Intercepts/means
mxEval(TX, phenofacmodel)
 
# Factor variances
mxEval(PH, phenofacmodel)
 
phenofacrun <- mxRun(phenofacmodel)
summary(phenofacrun, refModels=mxRefModels(phenofacrun, run = TRUE))