require(OpenMx) setwd("N:/Alyssa") myGrowthMixtureData <- read.csv("LCGA.csv") class1 <- mxModel("Class1", mxMatrix( type="Full", nrow=15, ncol=15, free= c(F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F), values=c(0,0,0,0,0,0,0, 1, 0, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 1, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0, 1,0, 0, 0,0, 0,0,0,0,0,0,0,0.5, 0, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0,0.5, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0,0.5, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0,0.5,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,1, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,1,0.8,0.64,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,1,1.1,1.21,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,1,1.4,1.96,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 0, 0, 0,0, 0, 0,0, 0,0,0,0,0,0,0, 0, 1, 0, 0,0, 0, 0,0), labels=c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "pcl", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "pcl", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "pcl", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "pcl", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), lbound=c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, -0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, -0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, -0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, -0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), ubound=c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 0.999, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), byrow=TRUE, name="A" ), mxMatrix( type="Symm", nrow=15, ncol=15, free= c(T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T), values=c(5,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,5,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,5, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0,0.1, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0,0.1, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0,0.1, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0,0.1,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,1,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,1,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,1,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,1,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,5), labels=c("resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var01", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var11", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var21", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var31", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS"), byrow= TRUE, name="S" ), mxMatrix( type="Full", nrow=7, ncol=15, free=F, values=c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0), byrow=T, name="F" ), mxMatrix( type="Full", nrow=1, ncol=15, values=c(1,1,1,3,3,3,3,0,0,0,0,1,1,-1,1), free=c(F,F,F,T,T,T,T,F,F,F,F,T,T,T,F), labels=c(NA,NA,NA,"meanpcl","meanpcl","meanpcl","meanpcl",NA,NA,NA,NA,"meani1","means1","meanq1",NA), name="M" ), mxRAMObjective("A","S","F","M", vector=TRUE, dimnames = c(names(myGrowthMixtureData),"PTSD0","PTSD1","PTSD2","PTSD3","i","s","q","phantom")) ) # close model #Create an MxModel object # ----------------------------------------------------------------------------- class2 <- mxModel(class1, mxMatrix( type="Symm", nrow=15, ncol=15, free= c(T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, F, T), values=c(5,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,5,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,5, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0,0.1, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0,0.1, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0,0.1, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0,0.1,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,1,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,1,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,1,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,1,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0, 0, 0,0,0,0,0,0,0,0,5), labels=c("resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resPCL", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var02", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var12", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var22", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "var32", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "resCAPS"), byrow= TRUE, name="S" ), mxMatrix( type="Full", nrow=1, ncol=15, values=c(1,1,1,3,3,3,3,0,0,0,0,1,1,-1,1), free=c(F,F,F,T,T,T,T,F,F,F,F,T,T,T,F), labels=c(NA,NA,NA,"meanpcl","meanpcl","meanpcl","meanpcl",NA,NA,NA,NA,"meani2","means2","meanq2","phantom"), name="M" ), name="Class2" ) # close model #Create an MxModel object # ----------------------------------------------------------------------------- classP <- mxMatrix("Full", 2, 1, free=c(TRUE, FALSE), values=1, lbound=0.001, labels = c("p1", "p2"), name="Props") classS <- mxAlgebra(Props%x%(1/sum(Props)), name="classProbs") # make a matrix of class probabilities # ----------------------------------------------------------------------------- algObj <- mxAlgebra(-2*sum( log(classProbs[1,1]%x%Class1.objective + classProbs[2,1]%x%Class2.objective)), name="mixtureObj") obj <- mxAlgebraObjective("mixtureObj") gmm <- mxModel("Growth Mixture Model", mxData( observed=myGrowthMixtureData, type="raw" ), class1, class2, classP, classS, algObj, obj ) gmmFit <- mxRun(gmm, suppressWarnings=TRUE) summary(gmmFit) gmmFit$classProbs omxCheckCloseEnough(gmmFit@output$Minus2LogLikelihood, 8739.05, 0.01) omxCheckCloseEnough(max(mxEval(classProbs, gmmFit)), 0.6009, 0.01) omxCheckCloseEnough(min(mxEval(classProbs, gmmFit)), 0.3991, 0.01) # Check results to see if they are within specified bounds