# Example 8.23: Expanding latent class model results

**SAS and R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In Example 8.21 we described how to fit a latent class model to data from the HELP dataset using SAS and R (using `poLCA()`, and then followed up in example 8.22 using `randomLCA()`. In both entries, we classified subjects 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 write a function to augment the default output of `randomLCA()` to make it easier for the analyst to interpret the results.

**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 will write a function wrapper for `randomLCA` that does some additional work in a generic fashion. This will allow easier estimation of other models. We annotate the function to explain what we’re doing. The resulting objects are `outcomep`, which contains the outcome probabilities, and `classp`, with the class probabilities.

runlca = function(df, nclass=2, names=c("item"), verbose=FALSE) { nvars = dim(df)[2] # create a list of names for the items if (length(names)==1) { names = rep(names, nvars) } # include only complete cases bigtable = table(na.omit(df)) allpatterns = as.data.frame(ftable(bigtable)) # keep only the patterns that occur nonzeropatterns = allpatterns[allpatterns$Freq > 0,] # fit the model results = randomLCA(nonzeropatterns[,1:nvars], nonzeropatterns$Freq, nclass=nclass, calcSE=FALSE) # display available sample size cat("nobs=", results$nobs, "\n") oldopt = options(digits=2) if (verbose==TRUE) { # display patterns whichclass = apply(results$classprob, 1, which.max) nonzeropatterns$class = whichclass print(nonzeropatterns[order(whichclass),]) } print(summary(results)) resvals = cbind(results$outcomep, results$classp) # label the margins with our desired variable names (plus class probability) colnames(resvals) = c(names, "classprob") # annotate standard output with rounded values print(round(resvals, 2)) options(oldopt) return(results) }

Now let’s apply the function. 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 to the function. Then we call the `runlca()` function. By specifying the `verbose` option the code displays each of the patterns, sorted by which class it is in (based on the highest predicted probability).

cesdcut = ifelse(cesd>20, 1, 0) smallds = data.frame(homeless, cesdcut, satreat, linkstatus) results = runlca(smallds, nclass=3, names=c("homeless", "cesd", "satreat", "linkstatus"), verbose=TRUE)

This generates the following output:

nobs= 431 homeless cesdcut satreat linkstatus Freq class 5 0 0 1 0 16 1 7 0 1 1 0 33 1 6 1 0 1 0 4 2 8 1 1 1 0 37 2 13 0 0 1 1 1 2 14 1 0 1 1 4 2 15 0 1 1 1 9 2 16 1 1 1 1 23 2 1 0 0 0 0 17 3 2 1 0 0 0 15 3 3 0 1 0 0 82 3 4 1 1 0 0 64 3 9 0 0 0 1 10 3 10 1 0 0 1 9 3 11 0 1 0 1 62 3 12 1 1 0 1 45 3 Classes AIC BIC logLik 3 2093 2150 -1032 Class probabilities Class 1 Class 2 Class 3 0.07846 0.21621 0.70534 Outcome probabilities homeless cesd satreat linkstatus classprob [1,] 0.00 0.58 1 0.00 0.08 [2,] 0.73 0.88 1 0.40 0.22 [3,] 0.44 0.83 0 0.41 0.71

The results are equivalent to the results from the prior example, but the predicted classes are listed, and the class probabilities (and proportion endorsing the item) are more clearly discernible. It might be useful in a later iteration of the function to add some blank lines and the proportion of the seeds that resulted in the maximum likelihood.

**leave a comment**for the author, please follow the link and comment on their blog:

**SAS and R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.