# # Copyright 2007-2009 The OpenMx Project # # 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 data <- read.table("lazarsfeld.ord",na.string=".", col.names=c("Armyrun", "Favatt", "squaredeal", "welfare", "freq")) freq <- data[,5] vars <- data[,1:4] nclass <- 2 nvar <- 4 nthresh <- 1 Mx1Threshold <- rbind(c(-.9, .3, -.3, .4)) nameList <- names(vars) # Define the model lcamodel <- mxModel("lca", Class1 <- mxModel("Class1", mxMatrix("Iden", name = "R", nrow = nvar, ncol = nvar, free=F), mxMatrix("Iden", name = "M", nrow = nvar, ncol = nvar, free=F), mxMatrix("Full", name = "ThresholdsClass1", nrow = nclass, ncol = 1, free=TRUE), mxFIMLObjective(covariance="R", means="M", dimnames=nameList, thresholds="ThresholdsClass1",vector=TRUE), mxData(vars, type="raw")), Class2 <- mxModel("Class2", mxMatrix("Iden", name = "R", nrow = nvar, ncol = nvar, free=TRUE), mxMatrix("Full", name = "ThresholdsClass2", nrow = nclass, ncol = 1, free=TRUE), mxFIMLObjective(covariance="R", means="M", dimnames=nameList, thresholds="ThresholdsClass2",vector=TRUE), 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, names = c(paste("pclass", 1:nclass, sep=""))), mxBound(c(paste("pclass", 1:nclass, sep="")),0,1), mxConstraint("1","=","sum(ClassMembershipProbabilities)"), # Define the objective function mxAlgebra(-2*sum(freq * log(pclass1%x%Class1.objective + pclass2%x%Class2.objective)), name="lca"), mxAlgebraObjective("lca")) # Run the job model <- mxRun(lca) summary(model) # Results from old Mx: #omxCheckCloseEnough(mxEval(thresh, model)[,1], Mx1Threshold[,1], 0.01) #omxCheckCloseEnough(mxEval(thresh, model)[1,2], Mx1Threshold[1,2], 0.01) #omxCheckCloseEnough(mxEval(R, model), Mx1R, 0.01) #omxCheckCloseEnough(model@output$Minus2LogLikelihood, 4081.48, 0.02)