You are here

Latent Profile Analysis / Finite Mixture Model

20 posts / 0 new
Last post
KB's picture
KB
Offline
Joined: 12/19/2018 - 13:53
Latent Profile Analysis / Finite Mixture Model

Hi,

I have been following the other post about latent class analysis using EasyMx https://openmx.ssri.psu.edu/comment/reply/4282/7227. I still have some questions and concerns. My goal is to run a two component finite mixture model (aka Latent Profile analysis with two classes).

Is the optimization being done gradient descent or some other algorithm like EM algorithm. I would like to be using the EM algorithm.

If Easy MX is not using the EM Algorithm is there a way to use the EM algorithm with EasyMx? Or is there a way to set up a latent profile analysis using the EM Algorithm. Thank you for your help!

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
EM

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

Veronica_echo's picture
Offline
Joined: 02/23/2018 - 01:57
A follow-up question

If I am using mxExpectationMixture() and mxFitFunctionML() to construct a finite mixture model, which kind algorithm I am applying? I assume it still should be sorts of iterative methods, but may I know more details, like the name or steps? Thanks in advance.

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
optimizer

The default optimizer is used. That's usually CSOLNP which does gradient descent.

Veronica_echo's picture
Offline
Joined: 02/23/2018 - 01:57
Got it. Thank you!

Got it. Thank you!

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
Note of clarification

Note of clarification: in some disciplines, "gradient descent" means "steepest descent", i.e. referring to algorithms that use only the first derivatives of the objective function. None of CSOLNP, SLSQP, or NPSOL is a steepest-descent optimizer. Rather, all 3 are quasi-Newton optimizers that use the BFGS approximation of the objective-function Hessian.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
mxComputeEM

Yeah. All of the test scripts that use mxComputeEM() are IRT models that use the IFA (item factor analysis) module. But it's possible to use mxComputeEM() with an arbitrary MxAlgebra, isn't it?

KB's picture
KB
Offline
Joined: 12/19/2018 - 13:53
OpenMx and EM Algorithm

Thanks for your feedback!

Is the general process the same? I build my model using pathways (as opposed to matrix algebra) then I just need to use mxComputeEM() to execute the optimization as opposed to something like mxExpectationMixture() and mxFitFunctionML(). I guess I am struggling to see where mxComputeEM fits in. Or is the process different and I actually wouldn't need to specify the model using pathways ? Thanks for your help!

KB's picture
KB
Offline
Joined: 12/19/2018 - 13:53
MxAlgebra

In order to use mxComputeEM() my model needs to be specified in terms of matrix algebras as opposed to the pathways? Is that correct?

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
compute plans

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.

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
what is the E step?

The first thing I want to know is what constitutes the E step in this model? Typically you'd create an mxAlgebra to perform the E step with fixed=TRUE. Then CompteEM can recompute this algebra with mxComputeOnce.

KB's picture
KB
Offline
Joined: 12/19/2018 - 13:53
E-step: an MxAlgebra Objects

The E step for a latent profile model is the calculating of the subject-specific posterior class probabilities where k is the number of classes. These probabilities are a function of the prior class probabilities (mixing proportions) and the subject normal densities/likelihoods which are a function of the class means and standard deviations, as shown below:
P(class=k|i, p_k, L_k) = p_k * L_ik / sum_k=1^m (p_k * L_ik)

So if I am understanding correctly, I need to create an MxAlgebra object which specifies the algebraic equation for these subject specific posterior probabilities. If that is correct, I have a few practical questions about creating this MxAlgebra (please stop me if I am way off base!) ... below.

I have created the prior probabilities using an MxMatrix object and I know I can use MxMatrix objects in the MxAlgebra ... so that is not an issue. What I am unsure of is how to refer to or call the data points in the likelihood (L_ik) as well as the means and standard deviations. My means and standard deviations are specified using MxPathways. So, should I refer to them in the equation by their labels? Can I refer to the data points in the likelihood using data.x (this is a guess)?

Thank you!

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
can you show the code?

I think the most efficient way to proceed is to exchange working code. Can you reply with your preliminary code attached and comments where you have specific questions?

Veronica_echo's picture
Offline
Joined: 02/23/2018 - 01:57
Thank you

Dr. Kirk,

Thank you very much! It's really very helpful. I need to add this page as my reference.

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
some progress

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)
AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
I just want to point out a

I just want to point out a crucial detail concerning Joshua's script: it gets different results from an MxModel that uses mxExpectationMixture() and mxFitFunctionML(). See the syntax below (which also avoids the label-switching problem by ensuring that class 2 is always the "responder" class):

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",
    data4mx,
    mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=mean(x),labels="mu1",name="Mu1"),
    mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=var(x),labels="sigmasq1",name="SigmaSq1"),
    mxExpectationNormal(covariance="SigmaSq1",means="Mu1",dimnames="x"),
    mxFitFunctionML(vector=T)
)
 
class2 <- mxModel(
    "Class2",
    data4mx,
    mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=mean(x),labels="mu1",name="Mu1"),
    mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.1,lbound=0,labels="deltamu",name="DeltaMu"),
    mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=var(x),labels="sigmasq2",name="SigmaSq2"),
    mxAlgebra(Mu1 + DeltaMu, name="Mu2"),
    mxExpectationNormal(covariance="SigmaSq2",means="Mu2",dimnames="x"),
    mxFitFunctionML(vector=T)
)
 
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)
mxEval(mean(Posteriors),mmfit,F)
 
 
mm2 <- mxModel(
    "Mixture2",
    data4mx,
    class1, class2,
    mxMatrix(type="Full",nrow=2,ncol=1,free=c(T,F),values=1,labels=c("p1","p2"),name="Props"),
    mxExpectationMixture(components=c("Class1","Class2"),weights="Props",scale="sum"),
    mxFitFunctionML()
)
mm2fit <- mxRun(mm2)
summary(mm2fit)
mm2fit$expectation$output
ChrisMcManus's picture
Offline
Joined: 06/18/2020 - 00:59
"fitfunction must be in probability units"

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

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Good question

Hello Chris - my undergrad prof from Bedford College no less! I have kicked your model around a bit and I don't see anything obvious yet. The error message makes me think that the likelihood of one of the classes (probably the first) is negative at some point, but both class1 and class2 run fine on their own, e.g., mxRun(ModelC2$class1).

Not sure if I ever stumped you on a question in class, but if I did you may be enjoying mentor's revenge :). I'll bring it up with the rest of the team, hopefully in tomorrow's meeting.

Best wishes
Mike

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Use vector=TRUE

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.

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Thanks Mike

Good stuff. I think we should try to improve the error message... "Very sorry, but it looks like you are trying to use a mixture distribution function, but at least some of the components of the model don't have the required argument to mxFitFunctionML(), i.e., they don't say mxFitFunctionML(vector=TRUE)." This is probably some variety of check prior to attempting to evaluate the likelihood.