You are here

2-class growth mixture model: adding time-invariant covariates & viewing class probabilities

5 posts / 0 new
Last post
PCFLaw's picture
Offline
Joined: 08/26/2020 - 03:15
2-class growth mixture model: adding time-invariant covariates & viewing class probabilities
AttachmentSize
Binary Data script.R5.22 KB
Plain text icon dataset.txt3.67 KB

Hi,

I am wondering if I can get some help on a 2-class growth mixture model I'm having trouble with (script and example dataset attached). The model describes a single linear trajectory across 8 time-points (t1-t8). I'm new to OpenMx and I've looked through this forum/other sites but I'm at a loss. I'm trying to (i) add 3 time-invariant covariates ('Sex', 'Age', 'pooled_VA') to the 2-class growth mixture model and (ii) view the class probabilities for each case.

Regarding problem (i), after I run the code below to fit the model with time-invariant covariates, the output I get when I run 'gmmFit$classProbs' is NULL. This surely can't be right but I don't know where I messed up.

classPM <- mxMatrix(type = "Full", nrow = 2, ncol = 4,
                    free=c(TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), values=1,
                    labels = c("p11","p21","p12", "p22",
                               "p13", "p23","p14", "p24"),
                    name = "weightsM")
classPV <- mxMatrix(nrow=4, ncol=1, labels=c(NA, "Sex", "pooled_VA", "Age"), values=1, name="weightsV")
# class proportions 
classP <- mxAlgebra(weightsM %*% weightsV, name="weights")
expec <- mxExpectationMixture(c('Class1', 'Class2'), weights='weights', scale='softmax')
objective <- mxFitFunctionML()
gmm <- mxModel("Growth Mixture Model", dataRaw, class1, class2, classPM, classPV, classP, expec, objective) # specify model to estimate
gmmFit <- mxRun(gmm) # run model
summary(gmmFit) # model parameter results

Regarding problem (ii), when I run the code below to print class probabilities, the output after running 'print(icp)' is "1" for one class and "0" for the other class for every case (see output below). This seems wrong but I don't know where I went wrong.

gmmFit$classProbs
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
  if (missing(model) || !(isS4(model) && is(model, "MxModel"))) {
    stop("'model' argument must be an MxModel object")
  }
  if (missing(classProbs) || !(is.character(classProbs) && length(classProbs) == 1)) {
    stop("'classProbs' argument must be a character string")
  }
  cp <- eval(substitute(mxEval(x, model), list(x = as.symbol(classProbs))))
  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)
}
 
icp <- indClassProbs(gmmFit, "classProbs")
print(icp)

[,1] [,2]
[1,] 1 0
[2,] 1 0
[3,] 1 0
[4,] 1 0
[5,] 1 0

I also get the above output when I fit the model without covariates using the code below.

# class proportions 
classP <- mxMatrix(type="Full", nrow=2, ncol=1, free=c(TRUE,FALSE), values=1, lbound=0.001,labels=c("p1","p2"),name="Props")
# class probabilities
classS <- mxAlgebra(Props%x%(1/sum(Props)), name="classProbs") 
algFit <- mxAlgebra(-2*sum(log(classProbs[1,1]%x%Class1.fitfunction+classProbs[2,1]%x%Class2.fitfunction)), name="mixtureObj")
fit <- mxFitFunctionAlgebra("mixtureObj") #This specifies the log likelhood function used to estimate class probabilities
gmm <- mxModel("Growth Mixture Model", dataRaw, class1, class2, classP, classS, algFit, fit) # specify model to estimate
gmmFit <- mxRun(gmm) # run model
summary(gmmFit) # model parameter results

Any assistance would be great. Thanks in advance!

Ryne's picture
Offline
Joined: 07/31/2009 - 15:12
So I do think the problem is

So I do think the problem is that you're hitting a boundary condition (probabilities of 1 or 0 are bad), but I'm not entirely sure where you went wrong. But I've never seen a matrix (rather than a vector) of class probabilities, and I don't understand the role that classPV is trying to play. But let's start with the model without covariates, as that should run successfully prior to adding more complexity.

I notice that your first model uses the old manual calculation, while the second using the mixture objective. Can you run the first model using the mixture objective and confirm that the class probabilities aren't 0 and 1?

PCFLaw's picture
Offline
Joined: 08/26/2020 - 03:15
Hi,

Hi,

I ran the first model using the mixture objective (code below; it ran successfully so hopefully I did it right), and the output is still 0 and 1.

classPM <- mxMatrix(type = "Full", nrow = 2, ncol = 1,
                    free=c(TRUE,FALSE), values=1,
                    labels = c("p1","p2"),
                    name = "weightsM")
classPV <- mxMatrix(nrow=1, ncol=1, labels=c(NA), values=1, name="weightsV")
# class proportions 
classP <- mxAlgebra(weightsM %*% weightsV, name="weights")
expec <- mxExpectationMixture(c('Class1', 'Class2'), weights='weights', scale='softmax')
objective <- mxFitFunctionML()
gmm <- mxModel("Growth Mixture Model", dataRaw, class1, class2, classPM, classPV, classP, expec, objective) # specify model to estimate
gmmFit <- mxRun(gmm) # run model
summary(gmmFit) # model parameter results

I'm not sure if it's relevant, but when gmmFit is run it gives me the following message:
Running Growth Mixture Model with 19 parameters
Warning message:
In model 'Growth Mixture Model' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Hi Could you please post the

Hi Could you please post the output from this command on your system?
mxVersion()

I tried running your script under R 4.0 with OpenMx and found an issue with the function indClassProbs() as follows:

Browse[2]> 
debug at #13: cp <- eval(substitute(mxEval(x, model), list(x = as.symbol(classProbs))))
Browse[2]> 
debug at #14: cp2 <- as.vector(cp)
Browse[2]> 
debug at #15: cps <- diag(length(cp2))
Browse[2]> cp
     [,1]        
[1,] "classProbs"
Browse[2]> ?substitute
Browse[2]> cp2
[1] "classProbs"
Browse[2]> cps
Error: object 'cps' not found
Browse[2]> diag(length(cp2))
     [,1]
[1,]    1

Which isn't what I think the code would be expecting, but might be wrong. Can you confirm that this is what was expected?

PCFLaw's picture
Offline
Joined: 08/26/2020 - 03:15
Hi,

Hi,

The following output when I run mxVersion() is:

OpenMx version: 2.17.4 [GIT v2.17.4]
R version: R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32
Default optimizer: CSOLNP
NPSOL-enabled?: No
OpenMP-enabled?: No

I'm sorry, I'm not familiar with the function indClassProbs(), so I don't know what I'm supposed to look for (and confirm whether it is what was expected). I just know that the output I'm supposed to see after running function indClassProbs() is not 1 and 0. I'm still very new to OpenMx.