Latent Profile Analysis / Finite Mixture Model

OpenMx offers mxComputeEM, but we don't have any demonstration scripts to show how to use it for latent profile analysis.

Knowing about compute plans would probably be useful for both Keighly and Jin. For the most commonly-used analyses done in OpenMx, the MxModel's compute plan is created automatically, without user involvement. Jin, you can see what the default compute plan would be, in the case of ML estimation and requesting confidence intervals, with

plan <- omxDefaultComputePlan(intervals=T)

and inspect the compute steps of the plan. They are, in order:
1. MxComputeGradientDescent, which is the primary quasi-Newton optimization to find point estimates.
2. MxComputeConfidenceInterval, which is a series of additional optimizations to find the confidence limits of profile-likelihood confidence intervals.
3. MxComputeNumericDeriv, during which the fitfunction gradient and Hessian are numerically calculated using central differences and several iterations of Richardson extrapolation.
4. MxComputeStandardError, during which standard errors of the point estimates are estimated from the diagonals of the inverted Hessian.
5. MxComputeHessianQuality, during which the condition number of the Hessian is computed.
6. MxComputeReportDeriv, during which the results of the numeric derivatives are reported to the frontend.
7. MxComputeReportExpectation, during which the MxModel's expectation object is populated with its results from the backend.

So, to use `mxComputeEM()`, the syntax might be something like

plan <- omxDefaultComputePlan()
plan$steps$GD <- mxComputeEM() #<--Need appropriate arguments to mxComputeEM()

in order to replace the gradient-descent step with an EM step, and then, put `plan` into your MxModel. There is probably more to it, though, since I do not think all of the default steps are compatible with EM. I've never actually used `mxComputeEM()` myself, so I'll have to ask Joshua, who implemented it, for the details.

There are some bugs in the v2.12.x that prevent this kind of model from working. However, things are working better in the master branch. Here's an example,


library(OpenMx)
set.seed(190127)

mu1 <- -1.5
mu2 <- 1.5

N <- 200
x <- matrix(c(rnorm(N/2,mu1,1),
rnorm(N/2,mu2,1)),ncol=1,dimnames=list(NULL,"x"))
data4mx <- mxData(observed=x,type="raw")

class1 <- mxModel("Class1",
mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=-.5,name="Mu"),
mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=4,name="Sigma"),
mxExpectationNormal(covariance="Sigma",means="Mu",dimnames="x"),
mxFitFunctionML(vector=T))

class2 <- mxRename(class1, "Class2")
class2$Mu$values[1,1] <- .5

mm <- mxModel(
"Mixture", data4mx, class1, class2,
mxAlgebra((1-Posteriors) * Class1.fitfunction, name="PL1"),
mxAlgebra(Posteriors * Class2.fitfunction, name="PL2"),
mxAlgebra(PL1 + PL2, name="PL"),
mxAlgebra(PL2 / PL, fixed=TRUE,
initial=matrix(runif(N,.4,.6), nrow=N, ncol = 1), name="Posteriors"),
mxAlgebra(-2*sum(log(PL)), name="FF"),
mxFitFunctionAlgebra(algebra="FF"),
mxComputeEM(
estep=mxComputeOnce("Mixture.Posteriors"),
mstep=mxComputeGradientDescent(fitfunction="Mixture.fitfunction")))

mmfit <- mxRun(mm)
summary(mmfit)

confusion <- table(mmfit$Posteriors$result < .1, c(rep(TRUE,N/2),rep(FALSE,N/2)))
print(confusion)


Hi -- Forgive a total newbie who has not used this site before if I have done things the wrong way, and offended any local etiquette... In particular I may be entirely in the wrong place, I realise.

I am using OpenMx to fit mixtures of normal distributions to data, and with univariate data it seemed to work well. Essentially, I started with the example of a mixture of growth models on p.89 of the User Guide (v2.17.2) and slowly altered things until I got the results I wanted with univariate data.

All was fine until I started working on some bivariate data and I kept running into a new error which said, "Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings, : ModelC2.fitfunction: component class1.fitfunction must be in probability units"

This I didn't understand at all as the example talks about fitting proportions and then making them into probabilities later. More frustratingly I have scoured every detail of my model and cannot find where on earth it is different from the GMM model, and neither can I find anything on how to make the fitfunction be in probability units.

At this point I should confess that I find the internal workings of OpenMx a little opaque, despite having used LISREL and laavan for many years (but they don't do mixture modelling!). So forgive once more.

I think I have attached a file which generates dummy data to fit two bivariate normal distributions to some data (and the plot shows they are clearly a mixture of two components).

Many thanks for any help you can offer,

Chris

PS OpenMx version: 2.17.4 [GIT v2.17.4]
R version: R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32
Default optimizer: CSOLNP
NPSOL-enabled?: No
OpenMP-enabled?: No

When using the mixture model you need to have mxFitFunctionML(vector=TRUE) in the component models. So change


class1 <- mxModel(model="class1", type="RAM", dataRaw, manifestVars=c("x","y"),
meanx1,meany1,varx1,vary1,cov1)

to

class1 <- mxModel(model="class1", type="RAM", dataRaw, manifestVars=c("x","y"),
meanx1,meany1,varx1,vary1,cov1, mxFitFunctionML(vector=TRUE))

and similarly for the class2 model. Everything should work fine then.

This means that you can't run them individually, but they will work in the mixture.