Posterior Class Probabilities and Entropy: Example

The first function calculates posterior class probabilities for a mixture model based on the following equation:

P(class=k|i, p_k, L_k) = p_k * L_ik / sum_k=1^m (p_k * L_ik)

where classes 1 through M are indexed by k, persons are indexed by i, p_k is the model estimated probability of being in class k and L_k is a vector of person specific likelihood for class k. Posterior class probabilities are given by Bayes theorum, with the equation above coming from Ramaswamy, Desarbo, Reibstein, and Robinson, 1993. In words, the individual class probabilities are the ratio of that person's raw likelihood for any one class divided by the sum of their likelihoods across all classes, each weighted by the model's estimated class probabilities for the sample.

The function below takes two arguments: an MxModel that was used to fit a mixture model, and the name of the matrix or algebra in that model that defines the class probabilities. It is assumed that the submodels in the provided model are the models for each of the classes that use FIML or RowObjective and return a vector of likelihoods. For instance, if you fit a mixture model with two classes, your structure should be that some object (''gmm'') that contains a model ("GrowthMixtureModel") with two submodels ("class1" and "class2") representing the models for each class. There should be a matrix in the parent model ("GrowthMixtureModel") that contains class probabilities.


indClassProbs <- function(model, classProbs, round=NA){
# this function takes a mixture model in OpenMx
# and returns the posterior class probabilities
# using Bayes rule, individual person-class likelihoods
# and the model class probability matrix, as described in
# Ramaswamy, Desarbo, Reibstein, and Robinson, 1993
cp <- mxEval(classProbs, model)
cp2 <- as.vector(cp)
cps <- diag(length(cp2))
diag(cps) <- cp2
subs <- model@submodels
if(min(dim(cp))!=1)stop("Class probabilities matrix must be a row or column vector.")
if(max(dim(cp))==1)stop("Class probabilities matrix must contain two or more classes.")
of <- function(num){
return(mxEval(objective, subs[[num]]))
}
rl <- sapply(1:length(names(subs)), of)
raw <- (rl%*%cps)
tot <- 1/apply(raw, 1, sum)
div <- matrix(rep(tot, length(cp2)), ncol=length(cp2))
icp <- raw * div
if (is.numeric(round)){icp <- round(icp, round)}
return(icp)
}