install.packages("OpenMx")
library(OpenMx)
library(haven)

# Script 7.2. OpenMx Script for the Two-Class Means Model

plan <- omxDefaultComputePlan()
plan$steps <- list(
  SA=mxComputeSimAnnealing(),
  GD=plan$steps$GD,
  ND=plan$steps$ND,
  SE=plan$steps$SE,
  HQ=plan$steps$HQ,
  RD=plan$steps$RD,
  RE=plan$steps$RE
)

class1 <- mxModel('Class1', type='RAM',
                  manifestVars=c('health_2','health_3','health_4','health_5','health_6','health_7','health_8','health_9'),
                  latentVars=c('eta_1','eta_2'),
                  # residual variances
                  mxPath(from=c('health_2','health_3','health_4','health_5','health_6','health_7','health_8','health_9'), arrows=2,
                         free=TRUE, values=1, labels='theta'),
                  # latent variances and covariance
                  mxPath(from=c('eta_1','eta_2'), arrows=2, connect='unique.pairs',
                         free=TRUE, values=c(1, 0.5, 1), labels=c('psi_11', 'psi_12', 'psi_22')),
                  # intercept loadings
                  mxPath(from='eta_1', to=c('health_2','health_3','health_4','health_5','health_6','health_7','health_8','health_9'), arrows=1,
                         free=FALSE, values=1),
                  # slope loadings
                  mxPath(from='eta_2', to=c('health_2','health_3','health_4','health_5','health_6','health_7','health_8','health_9'), arrows=1,
                         free=FALSE, values=c(0, 1,2,3,4,5,6,7)),
                  # latent variable means
                  mxPath(from='one', to=c('eta_1', 'eta_2'), arrows=1,
                         free=TRUE, values=c(4,0.5), labels=c('alpha_1_c1', 'alpha_2_c1')),
                  # ML Fit Function to run for each individual
                  mxFitFunctionML(vector=TRUE)
) # close model for class 1

class2 <- mxModel(class1,
                  # latent means
                  mxPath(from='one', to=c('eta_1', 'eta_2'), arrows=1,
                         free=TRUE, values=c(3,0.3), labels=c('alpha_1_c2', 'alpha_2_c2')),
                  name='Class2') # close model


classRP <- mxMatrix('Full', 2, 1, free=c(FALSE, TRUE),
                    values=1, lbound=0.001, labels=c('p1', 'p2'), name='RProps')
classP <- mxAlgebra(RProps %x% (1 / sum(RProps)), name='Props')
algObj <- mxAlgebra(-2* sum(log(Props[1,1] %x% Class1.objective + Props[2,1] %x% Class2.objective)),
                    name='mixtureObj')
obj <- mxFitFunctionAlgebra('mixtureObj')
gmm.2.means <- mxModel('2Class Means Growth Mixture Model',
                       mxData(observed=OpenMx_dfwide_poster1, type='raw'),
                       class1, class2, classRP, classP, algObj, obj,  mxComputeSimAnnealing(plan = plan, method='tsallis1996', control=list(stepsPerTemp=3)))

gmm.2.means.fit <- mxRun(gmm.2.means)
summary(gmm.2.means.fit)