#
# Copyright 2007-2019 by the individuals mentioned in the source code history
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
require(OpenMx)
#Latent class analysis example
#------------------------------------------------------+
# Marginal ml |
# data from http://www.people.vcu.edu/~nhenry/LSA50.htm|
# see |
# Latent Class Probs Item 1 Item 2 Item 3 Item 4 |
# .4243 .9240 .6276 .5704 .5125 |
# .5757 .4324 .1871 .1008 .0635 |
# Mx recovers: |
# .4440 0.90736 0.61444 0.55828 0.50176 |
# .5560 0.42838 0.18218 .09394 .05627 |
#------------------------------------------------------+
data("lazarsfeld", package ="OpenMx", verbose= TRUE)
data <- lazarsfeld
freq <- data[,5]
data[,1] <- as.ordered(data[,1])
data[,2] <- as.ordered(data[,2])
data[,3] <- as.ordered(data[,3])
data[,4] <- as.ordered(data[,4])
vars <- data[,1:4]
nclass <- 2
nvar <- 4
nthresh <- 1
nameList <- names(vars)
class1 <- mxModel("Class1",
mxMatrix("Iden", name = "R", nrow = nvar, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "M", nrow = 1, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "ThresholdsClass1", nrow = 1, ncol = nvar,
dimnames = list("Threshold",nameList), free=TRUE),
mxExpectationNormal(covariance="R", means="M", dimnames=nameList,
thresholds="ThresholdsClass1"),
mxFitFunctionML(vector=TRUE))
class2 <- mxModel("Class2",
mxMatrix("Iden", name = "R", nrow = nvar, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "M", nrow = 1, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "ThresholdsClass2", nrow = 1, ncol = nvar,
dimnames = list("Threshold",nameList), free=TRUE),
mxExpectationNormal(covariance="R", means="M", dimnames=nameList,
thresholds="ThresholdsClass2"),
mxFitFunctionML(vector=TRUE))
# Define the model
lcamodel <- mxModel("lcamodel", class1, class2, mxData(vars, type="raw"),
# Create class membership probabilities, constrain them to be in the range 0-1, and make their sum equal 1.0
mxMatrix("Full", name = "ClassMembershipProbabilities", nrow = nclass, ncol = 1, free=TRUE,
labels = c(paste("pclass", 1:nclass, sep=""))),
mxBounds(c(paste("pclass", 1:nclass, sep="")),0,1),
mxMatrix("Iden", nrow = 1, name = "constraintLHS"),
mxAlgebra(sum(ClassMembershipProbabilities), name = "constraintRHS"),
mxConstraint(constraintLHS == constraintRHS),
# Define the objective function
mxAlgebra(-2*sum(freq * log(pclass1%x%Class1.objective + pclass2%x%Class2.objective)),
name="lca"),
mxFitFunctionAlgebra("lca")
)
# Run the job
model <- mxRun(lcamodel, suppressWarnings=TRUE)
summary(model)
model$matrices
# Check results against those hard-coded from old Mx:
negativeThresholds <- as.vector(mxEval(Class1.ThresholdsClass1, model))
positiveThresholds <- as.vector(mxEval(Class2.ThresholdsClass2, model))
if (negativeThresholds[[1]] > 0) {
temp <- negativeThresholds
negativeThresholds <- positiveThresholds
positiveThresholds <- temp
}
omxCheckCloseEnough(negativeThresholds, c(-1.3247, -0.2909, -0.1466, -0.0044),.01)
omxCheckCloseEnough(positiveThresholds, c(0.1805, 0.9071, 1.3169, 1.5869),.01)
omxCheckCloseEnough(sort(as.vector(mxEval(ClassMembershipProbabilities, model))), c(0.4440, 0.5560),.01)
omxCheckCloseEnough(model$output$Minus2LogLikelihood, 4696.444, 0.01)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.