oneFactorThresholdModel <- mxModel("oneFactorThresholdModel", ### MALES #### mxModel("Males", ## Factor Loadings ## mxMatrix( type="Full", nrow=n, ncol=k, free=TRUE, name="LambdaM" ), ## Factor Variance ### mxMatrix( type="Full", nrow=k, ncol=k, free=FALSE, value=1, name="PsiM" ), ## Identity matrix ### mxMatrix( type="Iden", nrow=n, ncol=n, name="I" ), ## Residuals ## mxAlgebra( expression=vec2diag(diag2vec(I) - diag2vec((LambdaM %*% PsiM %*% t(LambdaM)))) , name="ThetaM" ), ## Expected Covariances ### ## Factor Variance ### mxAlgebra( expression=LambdaM %*% PsiM %*% t(LambdaM) + ThetaM, name="SigmaM" ), ## Observed means ## mxMatrix( type="Full", nrow=1, ncol=n, free=FALSE, values=c(meM), name="MuM"), ## Observed Covariances ## mxMatrix("Full", ncol=n, nrow=n, free=FALSE, values=c(covM), name="Sm"), ## Diff obs - Expected mxAlgebra( Sm - SigmaM, name="D"), ## Diff obs - Expected as vector mxAlgebra(vech(D), "vech_D"), # Weight matrix multiplied by n mxMatrix("Full", ncol=(n*n-n)/2+n, nrow=(n*n-n)/2+n, values=Nm*c(solve(acovM)), free=FALSE, name="acovM_S"), # Objective function mxAlgebra( t(vech_D) %*% acovM_S %*% vech_D , name="obj"), mxAlgebraObjective('obj') ), ### FEMALES ### mxModel("Females", ## Factor Loadings ## mxMatrix( type="Full", nrow=n, ncol=k, free=TRUE, name="LambdaF" ), ## Factor Variance ### mxMatrix( type="Full", nrow=k, ncol=k, free=FALSE, value=1, name="PsiF" ), ## Factor Variance ### mxMatrix( type="Diag", nrow=n, ncol=n, free=TRUE, lbound=0.01, name="ThetaF" ), ## Expected Covariances ### mxAlgebra( expression=LambdaF %*% PsiF %*% t(LambdaF) + ThetaF, name="SigmaF" ), ## Observed means ## mxMatrix( type="Full", nrow=1, ncol=n, free=FALSE, values=c(meF), name="MuF"), ## Observed Covariances ## mxMatrix("Full", ncol=n, nrow=n, free=FALSE, values=c(covF), name="Sf"), ## Diff obs - Expected mxAlgebra( Sf - SigmaF, name="D"), ## Diff obs - Expected as vector mxAlgebra(vech(D), "vech_D"), # Weight matrix multiplied by n mxMatrix("Full", ncol=(n*n-n)/2+n, nrow=(n*n-n)/2+n, values=Nf*c(solve(acovF)), free=FALSE, name="acovF_S"), # Objective function mxAlgebra( t(vech_D) %*% acovF_S %*% vech_D , name="obj"), mxAlgebraObjective('obj') ), mxAlgebra(Males.objective + Females.objective, name="FullModel"), mxAlgebraObjective("FullModel") ) fit <- mxRun(oneFactorThresholdModel) summary(fit)