library(OpenMx)

# put the data in OpenMx format
dataset <- read.delim(file.choose(), header = TRUE, sep = "\t")
dataset.gmm <- dataset[,c("Age", "Sex", "pooled_VA", "t1", "t2", "t3", "t4", "t5", "t6", "t7", "t8")] 
names(dataset.gmm) <- c("Age", "Sex", "pooled_VA", paste0("t", 1:8))
dataRaw <- mxData(observed=dataset.gmm, type="raw")

# residual variances
resVars <- mxPath(from=paste0("t", 1:8), to=paste0("t", 1:8), arrows=2,
                  free=TRUE, values=rep(1, 8),labels = paste0("residual", 1:8))
# latent variances and covariance (class 1)
latVars <- mxPath(from=c("intercept", "slope"), to=c("intercept", "slope"), arrows=2, connect="unique.pairs",
                  free=TRUE, values = c(1,.4,1), labels=c("vari1", "cov1", "vars1"))
# intercept loadings
intLoads <- mxPath(from="intercept", to=paste0("t", 1:8), arrows=1,
                   free=FALSE, values=rep(1, 8))
# slope loadings
sloLoads <- mxPath(from="slope", to=paste0("t", 1:8), arrows=1,
                   free=FALSE, values=seq(from=0,to=7))
# manifest means
manMeans <- mxPath(from="one", to=paste0("t", 1:8), arrows=1,
                   free=FALSE, values=rep(0, 8))
# latent variable means
latMeans <- mxPath(from="one", to=c("intercept","slope"), arrows=1,
                   free=TRUE, values=c(0,-1), labels=c("meanintercept", "meanslope"))
# enable the likelihood vector
funML <- mxFitFunctionML(vector=TRUE)
class1 <- mxModel("Class1", type="RAM",
                  manifestVars=paste0("t", 1:8),
                  latentVars=c("intercept", "slope"),
                  resVars, latVars, intLoads, sloLoads, manMeans, latMeans, funML)

# latent variances and covariance (class 2)
latVars2 <- mxPath(from=c("intercept", "slope"), to=c("intercept", "slope"), arrows=2, connect="unique.pairs",
                   free=TRUE, values=c(1, .5, 1), labels=c("vari2", "cov2", "vars2"))
latMeans2 <- mxPath(from="one", to=c("intercept", "slope"), arrows=1,
                    free=TRUE, values=c(5,1), labels=c("meani2", "means2"))
class2 <- mxModel(class1, name="Class2", latVars2, latMeans2)

---------------------------------------------------------------------
  # Model fitting without covariates
  
  # class proportions 
  classP <- mxMatrix(type="Full", nrow=2, ncol=1, free=c(TRUE,FALSE), values=1, lbound=0.001,labels=c("p1","p2"),name="Props")
# class probabilities
classS <- mxAlgebra(Props%x%(1/sum(Props)), name="classProbs") 
algFit <- mxAlgebra(-2*sum(log(classProbs[1,1]%x%Class1.fitfunction+classProbs[2,1]%x%Class2.fitfunction)), name="mixtureObj")
fit <- mxFitFunctionAlgebra("mixtureObj") #This specifies the log likelhood function used to estimate class probabilities
gmm <- mxModel("Growth Mixture Model", dataRaw, class1, class2, classP, classS, algFit, fit) # specify model to estimate
gmmFit <- mxRun(gmm) # run model
summary(gmmFit) # model parameter results

---------------------------------------------------------------------
  # Model fitting with covariates
  
  classPM <- mxMatrix(type = "Full", nrow = 2, ncol = 4,
                      free=c(TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), values=1,
                      labels = c("p11","p21","p12", "p22",
                                 "p13", "p23","p14", "p24"),
                      name = "weightsM")
classPV <- mxMatrix(nrow=4, ncol=1, labels=c(NA, "Sex", "pooled_VA", "Age"), values=1, name="weightsV")
# class proportions 
classP <- mxAlgebra(weightsM %*% weightsV, name="weights")
expec <- mxExpectationMixture(c('Class1', 'Class2'), weights='weights', scale='softmax')
objective <- mxFitFunctionML()
gmm <- mxModel("Growth Mixture Model", dataRaw, class1, class2, classPM, classPV, classP, expec, objective) # specify model to estimate
gmmFit <- mxRun(gmm) # run model
summary(gmmFit) # model parameter results

---------------------------------------------------------------------
  # 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)
}

icp <- indClassProbs(gmmFit, "classProbs")
print(icp)