################################################################################ # Growth Mixture Model Example # Author: Ryne Estabrook # Date: October 18, 2010 # Description: The following code outlines a growth mixture model with two # classes fit to 'myGrowthMixtureData' from the OpenMx library ################################################################################ library(OpenMx) data(myGrowthMixtureData) # build lgc templace lgc <- mxModel("LGC", mxMatrix("Full", 5, 2, values=c(1, 1, 1, 1, 1, 0:4), name="lambda"), mxMatrix("Symm", 2, 2, TRUE, labels=c("vi", "cis", "vs"), name="phi"), mxMatrix("Full", 1, 2, TRUE, labels=c("mi", "ms"), name="alpha"), mxMatrix("Diag", 5, 5, TRUE, values=1, labels="e", name="epsilon"), mxAlgebra(lambda %*% phi %*% t(lambda) + epsilon, name="cov"), mxAlgebra(alpha %*% t(lambda), name="mean"), mxFIMLObjective("cov", "mean", dimnames=names(myGrowthMixtureData), vector=TRUE) ) # use the template to make the first class model class1 <- mxModel(lgc, name="Class1") class1 <- omxSetParameters(class1, labels=c("vi", "cis", "vs", "mi", "ms"), values=c(1, .5, 1, 5, 1), newlabels=c("vi1", "cis1", "vs1", "mi1", "ms1") ) # use the template to make the second class model class2 <- mxModel(lgc, name="Class2") class2 <- omxSetParameters(class2, labels=c("vi", "cis", "vs", "mi", "ms"), values=c(1, .5, 2, 0, -1), newlabels=c("vi2", "cis2", "vs2", "mi2", "ms2") ) # make the class proportion matrix, # fixing one parameter at a non-zero constant (one) classP <- mxMatrix("Full", 2, 1, free=c(TRUE, FALSE), values=1, lbound=0.001, labels = c("p1", "p2"), name="Props") # rescale the class proportion matrix into # a class probability matrix by dividing by their sum # (done with a kronecker product of the class props and 1/sum) classS <- mxAlgebra(Props%x%(1/sum(Props)), name="classProbs") # calculate the sum of the class-specific objectives, # weighted by the class probabilities algObj <- mxAlgebra(-2*sum( log(classProbs[1,1]%x%Class1.objective + classProbs[2,1]%x%Class2.objective)), name="mixtureObj") # make an mxAlgebraObjective obj <- mxAlgebraObjective("mixtureObj") # put it all in a model gmm <- mxModel("Growth Mixture Model", mxData( observed=myGrowthMixtureData, type="raw" ), class1, class2, classP, classS, algObj, obj ) # run it gmmFit <- mxRun(gmm) # view the results summary(gmmFit) # view the class probabilities gmmFit$classProbs ################################################################################ # Model Stats and Interpretation ################################################################################ # view the class probabilities gmmFit$classProbs indClassProbs <- function(model, classProbs, round=NA){ # this function takes a mixture model in OpenMx # and returns the posterior class probabilities # using Bayes rule, individual person-class likelihoods # and the model class probability matrix, as described in # Ramaswamy, Desarbo, Reibstein, and Robinson, 1993 if (missing(model) || !(isS4(model) && is(model, "MxModel"))) { stop("'model' argument must be an MxModel object") } if (missing(classProbs) || !(is.character(classProbs) && length(classProbs) == 1)) { stop("'classProbs' argument must be a character string") } cp <- eval(substitute(mxEval(x, model), list(x = as.symbol(classProbs)))) cp2 <- as.vector(cp) cps <- diag(length(cp2)) diag(cps) <- cp2 subs <- model@submodels if(min(dim(cp))!=1)stop("Class probabilities matrix must be a row or column vector.") if(max(dim(cp))==1)stop("Class probabilities matrix must contain two or more classes.") of <- function(num) { return(mxEval(objective, subs[[num]])) } rl <- sapply(1:length(names(subs)), of) raw <- (rl %*% cps) tot <- 1 / apply(raw, 1, sum) div <- matrix(rep(tot, length(cp2)), ncol=length(cp2)) icp <- raw * div if (is.numeric(round)){icp <- round(icp, round)} return(icp) } entropy <- function(classProbs) { # this function takes a matrix of class probabilities from a # mixture model with people in rows and classes in columns # and returns the entropy statistic, as given by # Ramaswamy, Desarbo, Reibstein, and Robinson 1993 n <- dim(classProbs)[1] k <- dim(classProbs)[2] e <- 1 - sum(-classProbs * log(classProbs))/(n * log(k)) return(e) } icp <- indClassProbs(gmmFit, "classProbs") print(icp) entropy(icp) hist(icp[,1], main="Class Probabilities", xlab="P(Class==1)", ylab="Count", breaks=20) ################################################################################ # Serial Respecification ################################################################################ # how many trials? trials <- 20 # place all of the parameter names in a vector parNames <- names(omxGetParameters(gmm)) # make a matrix to hold all of the input <- matrix(NA, trials, length(parNames)) dimnames(input) <- list(c(1: trials), c(parNames)) output <- matrix(NA, trials, length(parNames)) dimnames(output) <- list(c(1: trials), c(parNames)) fit <- matrix(NA, trials, 4) dimnames(fit) <- list(c(1: trials), c("Minus2LL", "Status", "Iterations", "pclass1")) # populate the class probabilities input[,"p1"] <- runif(trials, 0.1, 0.9) input[,"p1"] <- input[,"p1"]/(1-input[,"p1"]) # populate the variances v <- c("vi1", "vs1", "vi2", "vs2", "e") input[,v] <- runif(trials*5, 0, 10) # populate the means m <- c("mi1", "ms1", "mi2", "ms2") input[,m] <- runif(trials*4, -5, 5) # populate the covariances r <- runif(trials*2, -0.9, 0.9) scale <- c( sqrt(input[,"vi1"]*input[,"vs1"]), sqrt(input[,"vi2"]*input[,"vs2"])) input[,c("cis1", "cis2")] <- r * scale # keep time serial <- proc.time()[3] for (i in 1: trials){ temp1 <- omxSetParameters(gmm, labels=parNames, values=input[i,] ) temp1@name <- paste("Starting Values Set", i) temp2 <- mxRun(temp1, unsafe=TRUE, suppressWarnings=TRUE, checkpoint=TRUE) output[i,] <- omxGetParameters(temp2) fit[i,] <- c( temp2@output$Minus2LogLikelihood, temp2@output$status[[1]], temp2@output$iterations, round(temp2$classProbs@result[1,1], 4) ) } serial <- proc.time()[3]-serial fit table(round(fit[,1], 3), fit[,2]) ################################################################################ # Parallel Respecification ################################################################################ library(snowfall) sfInit(parallel=TRUE, cpus=4) sfLibrary(OpenMx) topModel <- mxModel("Top") makeModel <- function(modelNumber){ temp <- mxModel(gmm, independent=TRUE, name=paste("Iteration", modelNumber, sep="")) temp <- omxSetParameters(temp, labels=parNames, values=input[modelNumber,]) return(temp) } mySubs <- lapply(1:20, makeModel) topModel@submodels <- mySubs results <- mxRun(topModel) fitStats <- function(model){ retval <- c( model@output$Minus2LogLikelihood, model@output$status[[1]], model@output$iterations, round(model$classProbs@result[1,1], 4) ) return(retval) } resultsFit <- t(omxSapply(results@submodels, fitStats)) sfStop() #### compare times parallel <- topModel@output$wallTime serial parallel serial-parallel