Example 8.23: Expanding latent class model results

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.

Example 8.21: latent class analysis

Latent class analysis is a technique used to classify observations based on patterns of categorical responses. Collins and Lanza's book,"Latent Class and Latent Transition Analysis," provides a readable introduction, while the UCLA ATS center has an online statistical computing seminar on the topic.

We consider an example analysis from the HELP dataset, where we wish to classify subjects based on their observed (manifest) status on the following variables: 1) on street or in shelter in past 180 days [homeless], 2) CESD score above 20, 3) received substance abuse treatment [satreat], or 4) linked to primary care [linkstatus]. We arbitrarily specify a three class solution.

SAS
Support for this method in SAS is available through the proc lca and proc lta add-on routines created and distributed by the Methodology Center at Penn State University. While it's customary in R to use researcher-written routines, it's less so for SAS; the machinery which allows independently written procs thus has the potential to mislead users. It bears explicitly stating that third-party procs probably don't have the same level of robustness or support as those distributed by SAS Institute.

The proc lca code assumes that the data exist in the dataset ds. The current coding of 0's and 1's needs to be changed to 1's and 2's.

data ds_0; set "c:\book\help.sas7bdat"; run;

data ds; set ds_0;
   homeless = homeless+1;
   cesdcut = (cesd > 20) + 1;
   satreat = satreat+1;
   linkstatus = linkstatus+1;
run;

The call to the LCA procedure specifies the number of classes, the variables to include, the number of categories per variable, and information about the starting values and random starts. It's highly recommended to run a "large" number of random starts to ensure that the true maximum likelihood estimate is reached (the 20 we used is likely too few for more complex models).

proc lca data=ds;
   title '3 class model';
   nclass 3;
   items homeless cesdcut satreat linkstatus;
   categories 2 2 2 2;
   seed 42;
   nstarts 20;
run;

The output begins with diagnostic information, and indicates that 40% of the seeds were associated with the best fitting model.

Data Summary, Model Information, and Fit Statistics (EM 
Algorithm)

Number of subjects in dataset:              431
Number of subjects in analysis:             431

Number of measurement items:             4
Response categories per item:            2 2 2 2
Number of groups in the data:            1
Number of latent classes:                3
Rho starting values were randomly generated (seed = 42).

No parameter restrictions were specified (freely estimated).

Seed selected for best fitted model:    1486228051
Percentage of seeds associated with best fitted model:   40.00%

The model converged in 3241 iterations.

Maximum number of iterations: 5000
Convergence method: maximum absolute deviation (MAD)
Convergence criterion:  0.000001000

A number of fit statistics are provided to help with model comparison (e.g. number of classes, constraints in more complex models).

=============================================
Fit statistics:
=============================================
Log-likelihood:     -1032.48
G-squared:              1.22
AIC:                   29.22
BIC:                   86.15
CAIC:                 100.15
Adjusted BIC:          41.72
Entropy R-sqd.:         0.94
Degrees of freedom:        1

The results indicate that 22% of subjects are in class 1, just 8% in class 2, and 70% in class 3.

                            Parameter Estimates
Gamma estimates (class membership probabilities): 
Class:                     1          2          3
                      0.2163     0.0785     0.7052

The next set of output describes the classes. The prevalence for each level of each variable is described for each class. The last response category is redundant (equal to 1 minus the sum of the other probabilities).

Rho estimates (item response probabilities): 
  Response category  1:
Class:                     1          2          3
  homeless    :       0.2703     1.0000     0.5625
  cesdcut     :       0.1154     0.4214     0.1678
  satreat     :       0.0004     0.0000     1.0000
  linkstatus  :       0.6029     1.0000     0.5855

  Response category  2:
Class:                     1          2          3
  homeless    :       0.7297     0.0000     0.4375
  cesdcut     :       0.8846     0.5786     0.8322
  satreat     :       0.9996     1.0000     0.0000
  linkstatus  :       0.3971     0.0000     0.4145

Members of class 1 were primarily homeless subjects with a larger proportion of high scores on the CESD, with substance abuse treatment history, and 40% of whom linked to primary care. Class 2 (the smallest group) was comprised of non-homeless subjects with lower CESD scores, substance abuse treatment, but no linkage. Class 3 was 44% homeless, had high levels of CESD, did not report substance abuse treatment, and 41% linked to primary care.

R

We begin by reading in the data, Then we use the within() function (section 1.3.1) to generate a dataframe with the variables of interest.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
ds = within(ds, (cesdcut = ifelse(cesd>20, 1, 0)))

The poLCA package supports estimation of latent class models in R. The poLCA() function, like proc lca, can incorporate polytomous categorical variables, but also like proc lca requires the variables to be coded starting with positive integers. We specify 10 repetitions (with random starting values).

library(poLCA)
res2 = poLCA(cbind(homeless=homeless+1, 
   cesdcut=cesdcut+1, satreat=satreat+1, 
   linkstatus=linkstatus+1) ~ 1, 
   maxiter=50000, nclass=3, 
   nrep=10, data=ds)

This generates the following output:

Model 1: llik = -1032.889 ... best llik = -1032.889
Model 2: llik = -1032.889 ... best llik = -1032.889
Model 3: llik = -1032.484 ... best llik = -1032.484
Model 4: llik = -1032.889 ... best llik = -1032.484
Model 5: llik = -1032.889 ... best llik = -1032.484
Model 6: llik = -1032.484 ... best llik = -1032.484
Model 7: llik = -1032.484 ... best llik = -1032.484
Model 8: llik = -1032.889 ... best llik = -1032.484
Model 9: llik = -1032.889 ... best llik = -1032.484
Model 10: llik = -1032.889 ... best llik = -1032.484
Conditional item response (column) probabilities,
 by outcome variable, for each class (row) 
 
$homeless
           Pr(1)  Pr(2)
class 1:  0.2703 0.7297
class 2:  1.0000 0.0000
class 3:  0.5625 0.4375

$cesdcut
           Pr(1)  Pr(2)
class 1:  0.1154 0.8846
class 2:  0.4213 0.5787
class 3:  0.1678 0.8322

$satreat
          Pr(1) Pr(2)
class 1:      0     1
class 2:      0     1
class 3:      1     0

$linkstatus
           Pr(1)  Pr(2)
class 1:  0.6029 0.3971
class 2:  1.0000 0.0000
class 3:  0.5855 0.4145

Estimated class population shares 
 0.2162 0.0785 0.7053 
 
Predicted class memberships (by modal posterior prob.) 
 0.181 0.1137 0.7053 
 
========================================================= 
Fit for 3 latent classes: 
========================================================= 
number of observations: 431 
number of estimated parameters: 14 
residual degrees of freedom: 1 
maximum log-likelihood: -1032.484 
 
AIC(3): 2092.967
BIC(3): 2149.893
G^2(3): 1.221830 (Likelihood ratio/deviance statistic) 
X^2(3): 1.233247 (Chi-square goodness of fit) 

The results are consistent with those found in proc lca. We note that, also similar to proc lca the global maximum likelihood estimates were reached 3 times out of 10-- this can be discerned by examination of the 10 model results. It's always a good idea to fit a large number of iterations to ensure that the global maximum likelihood estimates have been reached.