library(OpenMx) dat <- read.table(file = "LGC1.dat") names(dat) <- c("id", paste0("Y", 1:10), paste0("T", 1:10)) manifests <- paste0("Y", 1:10) latents <- c("eta0", "eta1") #Two-class Growth Mixture Model #---------------------------------------------------------------------------------- dataRaw <- mxData( observed = dat[, -1], type = "raw" ) # residual variances resVars <- mxPath( from = manifests, to = manifests, arrows = 2, free = TRUE, values = rep(1, 10), labels = paste0("r", 1:10)) # latent variances and covariance latVars <- mxPath( from = latents, to = latents, arrows = 2, connect = "unique.pairs", free = TRUE, values = c(25, 1, 2), labels = c("vari1", "cov1", "vars1") ) # intercept loadings intLoads <- mxPath( from = latents[1], to = manifests, arrows = 1, free = FALSE, values = rep(1, 10)) # slope loadings sloLoads <- mxPath( from = latents[2], to = manifests, arrows = 1, free = FALSE, values = 0, labels = paste0("data.T", 1:10)) # latent means latMeans <- mxPath( from = "one", to = latents, arrows = 1, free = TRUE, values = c(100, -8), labels = paste0("mu", latents)) # enable the likelihood vector funML <- mxFitFunctionML(vector = TRUE) class1 <- mxModel("Class1", type = "RAM", dataRaw, manifestVars = manifests, latentVars = latents, resVars, latVars, intLoads, sloLoads, latMeans, funML) # latent variances and covariance latVars2 <- mxPath( from = latents, to = latents, arrows = 2, connect = "unique.pairs", free = TRUE, values=c(25, 1, 2), labels = c("vari2","cov2","vars2") ) # latent means latMeans2 <- mxPath( from = "one", to = latents, arrows = 1, free = TRUE, values = c(5,1), labels = paste0("mu", latents, "2")) class2 <- mxModel(class1, name = "Class2", latVars2, latMeans2) classP <- mxMatrix( type = "Full", nrow = 2, ncol = 1, free=c(TRUE, FALSE), values=1, labels = c("p1","p2"), name = "weights" ) fit <- mxFitFunctionAlgebra("mixtureObj") gmm2 <- mxModel("GMM2Class", class1, class2, classP, mxExpectationMixture(paste0('Class',1:2), scale="softmax"), mxFitFunctionML()) gmm2Fit <- mxRun(gmm2, suppressWarnings=TRUE) summary(gmm2Fit) #----------------------------------------------------------------------------- #Three Class Growth Mixture Model #----------------------------------------------------------------------------- #---------------------------------------------------------------------------------- # residual variances resVars <- mxPath( from = manifests, to = manifests, arrows = 2, free = TRUE, values = rep(1, 10), labels = paste0("r", 1:10)) # latent variances and covariance latVars <- mxPath( from = latents, to = latents, arrows = 2, connect = "unique.pairs", free = TRUE, values = c(25, 1, 2), labels = c("vari1", "cov1", "vars1") ) # intercept loadings intLoads <- mxPath( from = latents[1], to = manifests, arrows = 1, free = FALSE, values = rep(1, 10)) # slope loadings sloLoads <- mxPath( from = latents[2], to = manifests, arrows = 1, free = FALSE, values = 0, labels = paste0("data.T", 1:10)) # latent means latMeans <- mxPath( from = "one", to = latents, arrows = 1, free = TRUE, values = c(100, -8), labels = paste0("mu", latents)) # enable the likelihood vector funML <- mxFitFunctionML(vector = TRUE) class1 <- mxModel("Class1", type = "RAM", dataRaw, manifestVars = manifests, latentVars = latents, resVars, latVars, intLoads, sloLoads, latMeans, funML) # latent variances and covariance latVars2 <- mxPath( from = latents, to = latents, arrows = 2, connect = "unique.pairs", free = TRUE, values=c(25, 1, 2), labels = c("vari2","cov2","vars2") ) # latent means latMeans2 <- mxPath( from = "one", to = latents, arrows = 1, free = TRUE, values = c(5,1), labels = paste0("mu", latents, "2")) class2 <- mxModel(class1, name = "Class2", latVars2, latMeans2) # latent variances and covariance latVars3 <- mxPath( from = latents, to = latents, arrows = 2, connect="unique.pairs", free = TRUE, values = c(25, 1, 2), labels = c("vari3","cov3","vars3") ) # latent means latMeans3 <- mxPath( from = "one", to = latents, arrows = 1, free = TRUE, values = c(5, 1), labels = paste0("mu", latents, "3")) class3 <- mxModel(class1, name="Class3", latVars3, latMeans3) classP <- mxMatrix( type="Full", nrow=3, ncol=1, free=c(TRUE, TRUE, FALSE), values=1, labels = c("p1","p2", "p3"), name="weights" ) gmm3 <- mxModel("GMM3Class", class1, class2, class3, classP, mxExpectationMixture(paste0('Class',1:3), scale="softmax"), mxFitFunctionML()) gmm3Fit <- mxRun(gmm3, suppressWarnings=TRUE) summary(gmm3Fit) #------------------------------------------------------------------------ mxCompare(gmm3Fit, gmm2Fit) mxCompare(gmm3Fit, gmm2Fit, boot=TRUE)