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

`proc`s thus has the potential to mislead users. It bears explicitly stating that third-party

`proc`s 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.

## 7 comments:

Is there a way to save class membership for each subject in poLCA ?

The "poLCA()" function returns an object of class "poLCA": one of the entries is called "posterior". This consists of a matrix of posterior class membership probabilities, and should do the trick as a way of merging back class membership for each subject.

How I should perform the trick? for merging back?

is there a way to save class membership using Proc LCA for each subject?

According to the FAQ for Proc LCA (http://methodology.psu.edu/ra/lcalta/faq):

"Can I obtain the predicted probability of membership in each latent class/status for each individual?

Yes. Posterior probabilities of class or status membership are available by using the OUTPOST option in PROC LCA and PROC LTA. See the user's guide for details on the syntax. If the goal is to assign individuals to a latent class based on their predicted probabilities and link class membership to outcomes, note that this approach does not incorporate the uncertainty of class membership into the analysis, thus biasing inference."

This is in response to April's "Anonymous" comment, and references the poLCA manual: http://userwww.service.emory.edu/~dlinzer/poLCA/poLCA-manual-1-3-1.pdf.

As Nick mentioned, the "poLCA()" function returns an object of class "poLCA", a list of 28 elements. One of these elements is "posterior", "an N by R matrix containing each observation’s posterior class membership probabilities". N is defined as the “number of cases used in the model,” and R is the number of classes specified in the call to poLCA(). It is worth noting that poLCA has a logical argument, na.rm, that specifies how to deal with missing values – if TRUE (the default), then these cases are listwise deleted before estimating the model.

Another interesting element of a “poLCA” object is “predclass,” a vector with integer values, of length N, “of predicted class memberships.”

I fully recommend reading the poLCA manual for more details; there is some latent class model theory but everything after chapter 5 should be of practical interest.

I have a few questions about adding complication to a baseline model:

1) Adding covariates:

a. I do not see G2 anymore in the output. The only fit stat I get it Log-likelihoo. Any way to get all the fit stats?

b. When running a model with covariates (even with the OUTPAR option) one can obtain a different solution compared to the baseline model without covariates. It seems that in this case arguing about the effects of covariates cannot be done with reference to classes identified in a baseline model, but the meaning of classes should be discussed from the model with covariate, right?

c. Other programs (e.g Latent Gold) allow to use covariates as ‘inactive’ so that they do not affect the solution, but they are only used to describe class membership probabilities. Is this possible in PROC LCA? And how can one judge absolute model fir when df is high (btw, what is the threshold after which we can say that G2 doesn’t follow a chi-square distribution?)

2) Multi-group LCA

a. I thought that adding a variable to define multiple-groups ,imposing measurement invariance, would be mathematically equivalent to adding that variable as covariate. However, from my empirical tries, I obtain different results.

b. when running a model with multiple groups and imposing measurement invariance, class membership probabilities are given for each group separately. Is it possible to get in the output the class membership probabilities for the whole sample?

Post a Comment