You are here

Latent Profile Analysis / Finite Mixture Model

16 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
Offline
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
Offline
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
Offline
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
Offline
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