Example 8.25: more latent class models (plus a graphical display)

In recent entries (here, here, here and here), we've been fitting a series of latent class models using SAS and R. One of the most commonly used and powerful package for latent class model estimation is Mplus. In this entry, we demonstrate how to use the MplusAutomation package to automate the process of fitting and interpreting a series of models using Mplus.

The first chunk of code needs to be run on a Windows computer with Mplus installed. We undertake the same analytic steps as before, then run the prepareMplusData() function to create the dataset, then createModels() to create the Mplus input files.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
library(MplusAutomation)
cesdcut = ifelse(cesd>20, 1, 0)
smallds = na.omit(data.frame(homeless, cesdcut, 
  satreat, linkstatus))
prepareMplusData(smallds, file="mplus.dat")
createModels("mplus.txt")

Once the preliminaries have been completed, we can run each of the models, then scarf up the results.

runModels()
summary=extractModelSummaries()
models=readModels()

To help create graphical summaries of the classes, I crafted some code that utilizes the Mplus output files. These include functions to calculate the antilogit (alogit()), determine class probabilities in terms of model parameters (calcclasses() and findprobs()), and a routine to plot the resulting values (plotres()). Unfortunately, these routines require tweaking for models with different number of predictors (as well as care if some of the predictors have more than 2 levels).

alogit = function(x) {
        return(exp(x)/(1+exp(x)))
}

calcclasses = function(vals) {
   numvals = length(vals)
   classes = numeric(numvals+1)   
   for (i in 1:numvals) {
      classes[i] = exp(vals[i])/
         (1+sum(exp(vals[1:numvals])))
   }
   classes[numvals+1] = 1/(1+sum(exp(vals[1:numvals])))
   return(classes)
}

findprobs = function(df) {
   numclass = length(levels(as.factor(df$LatentClass)))-1
   v1 = numeric(numclass)
   v2 = numeric(numclass)
   v3 = numeric(numclass)
   v4 = numeric(numclass)
   for (i in 1:numclass) {
      v1[i] = alogit(-df$est[4*(i-1)+1])
      v2[i] = alogit(-df$est[4*(i-1)+2])
      v3[i] =  alogit(-df$est[4*(i-1)+3])
      v4[i] =  alogit(-df$est[4*(i-1)+4])
   }
   if (numclass>1) {
      classes=calcclasses(df$est[(4*numclass+1):(4*numclass+numclass-1)])
   } else classes=1
   return(list(prop=cbind(v1=v1, v2=v2, v3=v3, v4=v4), 
     classes=classes))
}

plotres = function(mylist, roundval=1, cexval=.75) {
# can only plot at most 4 classes!
   reorder = order(mylist$classes)
   probs = mylist$classes[reorder]
   results = mylist$prop
   dimen = dim(results)
   cols = c(1,4,2,5) # black, blue, red, and turquoise
   ltys = c(3,1,2,4) # dotted, solid, dash, dash-dot
   ltyslines = ltys[rank(-mylist$classes)]
   colslines = cols[rank(-mylist$classes)]
   ltysrev = rev(ltys[1:dimen[1]])
   colsrev = rev(cols[1:dimen[1]])
   plot(c(0.9, dimen[2]), c(-0.08,1), xlab="", 
           ylab="estimated prevalence", xaxt="n", type="n")
   abline(h=0, col="gray")
   abline(h=1, col="gray")
   for (i in 1:(dimen[1])) {
      lines(1:(dimen[2]), results[i,], lty=ltyslines[i], 
         col=colslines[i], lwd=2)
      points(1:(dimen[2]), results[i,], col=colslines[i])
   }
   text(1,0,"homeless", pos=1, cex=cexval)
   text(2,0,"cesdcut", pos=1, cex=cexval)
   text(3,0,"satreat", pos=1, cex=cexval)
   text(4,0,"linkstat", pos=1, cex=cexval)
   legendval = character(dimen[1])
   for (i in 1:(dimen[1])) {
      legendval[i] = paste("class ",i , " (", round(100*probs[i], 
         roundval), "%)", sep="")
   }
   legend(2, 0.4, legend=legendval, 
      lty=ltysrev, col=colsrev, cex=cexval, lwd=2)
}

files=list.files()
files=files[grep(".out",files)
]
par(mfrow=c(4,1), mar=c(1, 2, 2, 1) + 0.1)
for (i in 1:length(files)) {
   cat("file=",files[i],"\n")
   res = extractModelParameters(target=files[i])
   newres = findprobs(res$unstandardized)
   plotres(newres)
   title(substr(files[i], 1, nchar(files[i])-4))
}

Finally, by spelunking through the output files in the current directory (using list.files()) the results from each of the four models can be collated and displayed as seen in the Figure above.

The single class model simply reproduces the prevalences of each of the predictors. The two class solution primarily separates those not receiving substance abuse treatment from those that do. The three class solution further splits along substance abuse treatment, as well as homeless status. The four class solution is somewhat jumbled, with a group of homeless subjects comprising the largest class. None of the classes distinguish linkstatus (the primary outcome of the RCT).

Example 8.22: latent class modeling using randomLCA

In Example 8.21 we described how to fit a latent class model to data from the HELP dataset using SAS and R. Subjects were classified based on their observed (manifest) status on the following variables (on street or in shelter in past 180 days [homeless], CESD scores above 20, received substance abuse treatment [satreat], or linked to primary care [linkstatus]). We arbitrarily specify a three class solution.

In this example, we fit the same model using the randomLCA() function within the package of the same name.

R

We begin by reading in the data.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
library(randomLCA)

We start by creating a dichotomous variable with high scores on the CESD, and put this together as part of a dataframe to be given as input.

cesdcut = ifelse(cesd>20, 1, 0)
smallds = na.omit(data.frame(homeless, cesdcut, satreat, linkstatus))
results = randomLCA(smallds, nclass=3, notrials=1000)
summary(results)

This generates the following output:

  Classes      AIC      BIC    logLik
        3 2092.968 2149.893 -1032.484
Class probabilities 
Class  1 Class  2 Class  3 
 0.07846  0.70534  0.21620 
Outcome probabilities 
          homeless cesdcut   satreat linkstatus
Class  1 9.465e-06  0.5786 1.000e+00  9.538e-06
Class  2 4.375e-01  0.8322 9.988e-06  4.145e-01
Class  3 7.297e-01  0.8846 1.000e+00  3.971e-01

The results are equivalent to the results from the prior example, though the scientific notation for the observed prevalences in each class are hard to read. Other objects are available from the returned value, though they are also not in an easily digestible form:

> names(results)
 [1] "fit"        "nclass"     "classp"     "outcomep"  
 [5] "se"         "np"         "nobs"       "logLik"    
 [9] "observed"   "fitted"     "deviance"   "classprob" 
[13] "bics"       "random"     "level2"     "byclass"   
[17] "blocksize"  "call"       "probit"     "quadpoints"
[21] "patterns"   "notrials"   "freq"      
> results$patterns
   homeless cesdcut satreat linkstatus
1         0       0       0          0
2         0       0       0          1
3         0       0       1          0
4         0       0       1          1
5         0       1       0          0
6         0       1       0          1
7         0       1       1          0
8         0       1       1          1
9         1       0       0          0
10        1       0       0          1
11        1       0       1          0
12        1       0       1          1
13        1       1       0          0
14        1       1       0          1
15        1       1       1          0
16        1       1       1          1
> results$freq
 [1] 17 10 16  1 82 62 33  9 15  9  4  4 64 45 37 23

We'll address workarounds for these shortcomings in a future entry.