3 Auxiliary Variable Analyses

Latent class and latent profile models are often only the first step in an analysis. After identifying a useful class solution, researchers commonly want to know whether the classes differ on variables that were not used to define the classes. These variables are often called auxiliary variables, distal outcomes, or external variables.

In this vignette you will learn how to:

distinguish proper auxiliary variable analysis from naive approaches assuming perfect class assignment;
inspect classification uncertainty before conducting auxiliary analyses;
use the Bolck-Croon-Hagenaars (BCH) method implemented in tidySEM::BCH();
use pseudo-class draws with tidySEM::pseudo_class() and tidySEM::append_class_draws();
test overall and pairwise class differences in BCH models using tidySEM::lr_test();
extract and report model parameters using tidySEM::table_results(); and
communicate auxiliary-variable results without treating latent class membership as perfectly observed.

To follow along and run the code yourself, open the file auxiliary_variables.qmd in the lca_course project, as explained in Appendix A.

3.1 Why auxiliary-variable analyses require special care

Suppose a latent profile model identifies three classes of students based on several indicators of study behavior. You may then want to ask whether the classes differ in academic achievement, well-being, or another external variable.

At face value, it seems sensible to assign each student to their most likely class, and then run an ANOVA or regression to examine class differences in academic achievement. The difficulty is that latent class membership is not directly observed, so there will be error/uncertainty in these class assignments. Treating the most likely class as if it were known without error discards this uncertainty, leading to biased estimates and standard errors in downstream analyses.

Proper methods for auxiliary variable analysis account for this classification uncertainty. This vignette focuses on two documented approaches:

BCH analysis. BCH() estimates an auxiliary model as a multi-group model. Every case contributes to every class-specific group, with contributions determined by group-specific BCH weights.
Pseudo-class analysis. pseudo_class() repeatedly draws a class value for each person according to that person’s estimated class-membership probabilities, estimates the auxiliary model in each resulting dataset, and pools the results across datasets.

3.2 Example Data and Model

This vignette assumes that you have already selected and interpreted a latent class solution estimated with a tidySEM mixture-model function such as mx_profiles(), mx_lca(), or mx_mixed_lca(), see Chapter 1. You can use your own, or follow along with the examples below using the zegwaard_carecompass data included with tidySEM, see Chapter 2. We can load the final model constructed in that previous chapter from a file. This gives us four (theoretically expected) profiles of caregivers of close others (e.g., the elderly, or those suffering from chronic illness). The profiles are constructed based on self-reported burden, feelings of being trapped, negative affect, and loneliness. The dataset contains several other properties of the caregiver and the person being cared for, which we could analyze as auxiliary variables:

library(tidySEM)
library(OpenMx)
lpa_model <- readRDS("res_zegwaard.rds")
df <- zegwaard_carecompass
names(df)

 [1] "burdened"     "trapped"      "negaffect"    "loneliness"   "sex"         
 [6] "sexpatient"   "cohabiting"   "distance"     "freqvisit"    "relationship"

First examine the latent profile model.

table_fit(lpa_model)

Minus2LogLikelihood	n	Parameters	observedStatistics	df	RMSEASquared	RMSEANull	modelName	AIC	BIC	saBIC	Classes	Entropy	prob_min	prob_max	n_min	n_max	LL
3832.507	513	23	2030	2007	0	0.05	mix4	3878.507	3976.033	3903.028	4	0.75487	0.8142145	0.9151939	0.1559454	0.3430799	-1916.253

table_results(lpa_model)

label	est_sig	se	pval	confint	class
Means.burdened.class1	3.27***	0.04	0.00	[3.18, 3.36]	class1
Means.trapped.class1	1.28***	0.05	0.00	[1.18, 1.38]	class1
Means.negaffect.class1	2.31***	0.06	0.00	[2.20, 2.42]	class1
Means.loneliness.class1	2.73***	0.04	0.00	[2.64, 2.82]	class1
Variances.burdened.class1	0.23***	0.02	0.00	[0.19, 0.27]	class1
Variances.trapped.class1	0.17***	0.02	0.00	[0.14, 0.20]	class1
Variances.negaffect.class1	0.31***	0.02	0.00	[0.27, 0.36]	class1
Variances.loneliness.class1	0.24***	0.02	0.00	[0.20, 0.28]	class1
Means.burdened.class2	3.40***	0.06	0.00	[3.28, 3.52]	class2
Means.trapped.class2	2.27***	0.06	0.00	[2.15, 2.38]	class2
Means.negaffect.class2	2.81***	0.06	0.00	[2.70, 2.93]	class2
Means.loneliness.class2	2.79***	0.06	0.00	[2.66, 2.91]	class2
Means.burdened.class3	4.25***	0.07	0.00	[4.12, 4.38]	class3
Means.trapped.class3	2.67***	0.05	0.00	[2.58, 2.77]	class3
Means.negaffect.class3	2.92***	0.06	0.00	[2.80, 3.03]	class3
Means.loneliness.class3	2.01***	0.06	0.00	[1.89, 2.14]	class3
Means.burdened.class4	2.38***	0.06	0.00	[2.26, 2.50]	class4
Means.trapped.class4	0.38***	0.05	0.00	[0.28, 0.49]	class4
Means.negaffect.class4	1.78***	0.07	0.00	[1.65, 1.91]	class4
Means.loneliness.class4	3.18***	0.06	0.00	[3.07, 3.30]	class4
mix4.weights[1,1].NA	1.00	NA	NA	NA	NA
mix4.weights[1,2].NA	0.86***	0.15	0.00	[0.56, 1.15]	NA
mix4.weights[1,3].NA	0.66***	0.11	0.00	[0.44, 0.88]	NA
mix4.weights[1,4].NA	0.47***	0.08	0.00	[0.32, 0.63]	NA

Quiz

How many classes does the model contain?
The class-specific model parameters include means and variances.
Auxiliary variables influence the parameter estimates of the class solution
The fit table suggests that class separation is:

3.3 Inspect Classification Uncertainty

Before testing external differences, inspect how clearly observations are classified. Note that the fit table already gave us two metrics: entropy, and the range of the probabilities that observations in each class will be correctly classified (prob_min and prob_max). Note that these are all reasonably low, indicating poor classification accuracy. Thus, this is a prime example of an analysis where it’s important to account for classification uncertainty! If entropy and prob_min and prob_max were all near 1, you could treat class membership as observed without introducing substantial bias, but that is not the case here.

The function class_prob() returns posterior classification probabilities for each individual, and several summaries and diagnostic statistics based on those:

cp <- class_prob(lpa_model)

Examine the individual posterior probabilities:

format_numeric(head(cp$individual))

     class1 class2 class3 class4 predicted
[1,] "0.87" "0.13" "0.00" "0.00" "1.00"   
[2,] "0.28" "0.00" "0.00" "0.72" "4.00"   
[3,] "0.81" "0.19" "0.00" "0.00" "1.00"   
[4,] "0.96" "0.02" "0.00" "0.02" "1.00"   
[5,] "0.92" "0.07" "0.00" "0.01" "1.00"   
[6,] "0.01" "0.05" "0.93" "0.00" "3.00"

These probabilities make clear why an auxiliary analysis should not normally be based only on a hard class assignment. Whereas you might be relatively certain that person 4 belongs in class 1, there is substantial uncertainty about person 2, for example.

If you perform auxiliary analyses, the basic intuition is that each person contributes to estimating the parameters (e.g., auxiliary variable means) in each class in proportion to their posterior classification probability of belonging to that class. Thus, person 1 will contribute most to class 1, but also a little bit to class 2.

3.4 The BCH method

The BCH method performs a weighted multi-group analysis. All cases contribute to all groups, but their contributions are weighted by the most likely class membership matrix:

format_numeric(cp$mostlikely.class)

          assigned.1 assigned.2 assigned.3 assigned.4
avgprob.1 "0.88"     "0.08"     "0.00"     "0.03"    
avgprob.2 "0.11"     "0.81"     "0.07"     "0.00"    
avgprob.3 "0.01"     "0.13"     "0.86"     "0.00"    
avgprob.4 "0.08"     "0.00"     "0.00"     "0.92"

3.4.1 Comparing class means on a continuous outcome

If the research question is simply whether latent classes differ in the mean of a continuous external variable, for example, to test whether the classes differ in their distance to the caree, we can call the BCH function and simply supply the auxiliary variable to the data argument, omitting the model argument:

bch_mean <- BCH(x = lpa_model,
               data = df$distance)
table_results(bch_mean)

label	est_sig	se	confint	group
Means.y.class1	155.25***	3.80	[147.80, 162.69]	class1
Variances.y.class1	2464.08***	266.68	[1941.39, 2986.77]	class1
Means.y.class2	159.52***	2.79	[154.05, 164.98]	class2
Variances.y.class2	1144.27***	133.43	[882.76, 1405.78]	class2
Means.y.class3	147.24***	6.09	[135.30, 159.18]	class3
Variances.y.class3	4200.66***	558.51	[3105.99, 5295.32]	class3
Means.y.class4	167.01***	7.06	[153.17, 180.84]	class4
Variances.y.class4	3989.68***	630.60	[2753.71, 5225.64]	class4

Quiz

Which class, on average, lives closest to the caree?
The parameters in this model are:

A BCH analysis is a multi-group OpenMx model but not a test of differences between groups. The function lr_test() can be used to compare parameters across groups. The function conducts overall and pairwise likelihood-ratio tests for a multi-group MxModel whose submodels have the same structure:

lr_test(bch_mean, compare = "M")

BCH test for equality across classes

Overall likelihood ratio test:
 LL_baseline LL_restricted   LL_dif df         p
    5426.959      5432.425 5.465735  3 0.1407025

Pairwise comparisons using likelihood ratio tests:
 Model1 Model2 LL_baseline LL_restricted    LL_dif df          p
 class1 class2    5426.959      5427.779 0.8197093  1 0.36526517
 class1 class3    5426.959      5428.203 1.2440167  1 0.26469835
 class2 class3    5426.959      5430.281 3.3217853  1 0.06836755
 class1 class4    5426.959      5429.092 2.1334872  1 0.14411257
 class2 class4    5426.959      5427.933 0.9735965  1 0.32378486
 class3 class4    5426.959      5431.396 4.4373903  1 0.03515995

3.4.2 Comparing Proportions Across Classes

For a single ordinal variable, like the sex of patients, we could call:

aux_sex <- BCH(lpa_model, data = ordered(zegwaard_carecompass$sexpatient, levels = c("male", "female")))

To obtain an omnibus likelihood ratio test of the significance of these sex differences across classes, as well as pairwise comparisons between classes, use:

lr_test(aux_sex)

BCH test for equality across classes

Overall likelihood ratio test:
 LL_baseline LL_restricted   LL_dif df          p
    695.2568       703.954 8.697276  3 0.03359867

Pairwise comparisons using likelihood ratio tests:
 Model1 Model2 LL_baseline LL_restricted    LL_dif df           p
 class1 class2    695.2568      695.7158 0.4590770  1 0.498055675
 class1 class3    695.2568      700.2993 5.0425058  1 0.024732690
 class2 class3    695.2568      703.0550 7.7982664  1 0.005229638
 class1 class4    695.2568      695.9637 0.7068917  1 0.400477634
 class2 class4    695.2568      697.1501 1.8933247  1 0.168827407
 class3 class4    695.2568      696.4700 1.2132661  1 0.270686249

The results can be reported in probability scale:

table_prob(aux_sex)

Variable	Category	Probability	group
y	1	0.4596910	class1
y	2	0.5403090	class1
y	1	0.4216089	class2
y	2	0.5783911	class2
y	1	0.5959575	class3
y	2	0.4040425	class3
y	1	0.5164841	class4
y	2	0.4835159	class4

Quiz

There are significant sex differences across classes.
Which class differs significantly from class 1?
Which class cares for the largest proportion of female patients?

3.4.3 Comparing Auxiliary Models Across Classes

We can also compare a simple model between classes. For example, we previously compared mean distances across classes - but this model also left variances free across classes. If we want to constrain those, we could use the following code:

df_aux <- zegwaard_carecompass[, "distance", drop = FALSE]
aux_mean <- BCH(lpa_model, model = "distance~1; distance ~~c*distance",
                 data = df_aux)
table_results(aux_mean)

label	est_sig	se	confint	group
Means.distance.class1	155.25***	3.98	[147.44, 163.05]	class1
Variances.distance.class1	2707.70***	169.39	[2375.70, 3039.70]	class1
Means.distance.class2	159.52***	4.29	[151.11, 167.93]	class2
Means.distance.class3	147.24***	4.89	[137.65, 156.83]	class3
Means.distance.class4	167.01***	5.82	[155.61, 178.40]	class4

For a slightly more complex model, let’s examine whether the distance predicts the frequency of visits differently across classes (treated as continuous).

df_aux <- zegwaard_carecompass[, c("freqvisit", "distance")]
df_aux$freqvisit <- as.numeric(df_aux$freqvisit)
aux_model <- BCH(lpa_model, model = "freqvisit ~ distance",
                 data = df_aux)
table_results(aux_model)

label	est_sig	se	pval	confint	group
Regressions.freqvisit.ON.distance.class1	0.00	0.00	0.80	[-0.00, 0.00]	class1
Means.freqvisit.class1	3.99***	0.18	0.00	[3.64, 4.35]	class1
Means.distance.class1	155.25***	3.80	0.00	[147.80, 162.70]	class1
Variances.freqvisit.class1	0.53***	0.06	0.00	[0.42, 0.64]	class1
Variances.distance.class1	2464.07***	266.70	0.00	[1941.34, 2986.79]	class1
Regressions.freqvisit.ON.distance.class2	0.00	0.00	0.77	[-0.00, 0.01]	class2
Means.freqvisit.class2	3.66***	0.43	0.00	[2.81, 4.51]	class2
Means.distance.class2	159.52***	2.79	0.00	[154.05, 164.98]	class2
Variances.freqvisit.class2	1.19***	0.14	0.00	[0.92, 1.46]	class2
Variances.distance.class2	1144.27***	133.42	0.00	[882.77, 1405.76]	class2
Regressions.freqvisit.ON.distance.class3	-0.00	0.00	0.35	[-0.00, 0.00]	class3
Means.freqvisit.class3	3.95***	0.27	0.00	[3.43, 4.47]	class3
Means.distance.class3	147.24***	6.09	0.00	[135.30, 159.18]	class3
Variances.freqvisit.class3	1.29***	0.17	0.00	[0.95, 1.63]	class3
Variances.distance.class3	4200.66***	558.63	0.00	[3105.75, 5295.56]	class3
Regressions.freqvisit.ON.distance.class4	-0.00	0.00	0.91	[-0.00, 0.00]	class4
Means.freqvisit.class4	3.32***	0.38	0.00	[2.57, 4.08]	class4
Means.distance.class4	167.02***	7.06	0.00	[153.18, 180.86]	class4
Variances.freqvisit.class4	1.48***	0.23	0.00	[1.02, 1.93]	class4
Variances.distance.class4	3989.68***	630.52	0.00	[2753.87, 5225.48]	class4

To obtain an omnibus likelihood ratio test of the difference in regression coefficients across classes and pairwise comparisons between classes, use:

lr_test(aux_model, compare = "A")

BCH test for equality across classes

Overall likelihood ratio test:
 LL_baseline LL_restricted    LL_dif df        p
    6852.492      6853.473 0.9812781  3 0.805782

Pairwise comparisons using likelihood ratio tests:
 Model1 Model2 LL_baseline LL_restricted     LL_dif df         p
 class1 class2    6852.492      6852.521 0.02932549  1 0.8640297
 class1 class3    6852.492      6853.328 0.83645365  1 0.3604130
 class2 class3    6852.492      6853.043 0.55116219  1 0.4578432
 class1 class4    6852.492      6852.540 0.04766352  1 0.8271800
 class2 class4    6852.492      6852.582 0.08964307  1 0.7646314
 class3 class4    6852.492      6852.721 0.22892426  1 0.6323226

The results indicate that there are no significant sex differences across classes, \(\Delta LL(3) = 0.98, p = .81\).

3.5 Pseudo-class analysis

The pseudo-class approach represents classification uncertainty through repeated stochastic class assignments. For each generated dataset, a person’s class value is drawn according to that person’s posterior probability of belonging to each latent class. The auxiliary model is then fitted in each dataset, and the estimates are pooled across datasets.

For example, we can regress the external outcome on pseudo-class membership (treated as a factor via I(factor()), and dropping the intercept via -1:

pseudo_mean <- pseudo_class(
  x = lpa_model,
  model = lm(distance ~ I(factor(class)) -1, data = data),
  data = df[, "distance", drop = FALSE],
  m = 100)
pseudo_mean

term	estimate	std.error	statistic	df
I(factor(class))1	155.7188	4.309138	36.13689	406.1078
I(factor(class))2	158.3141	4.709850	33.61340	384.6969
I(factor(class))3	149.2713	5.356927	27.86510	387.0363
I(factor(class))4	165.2482	6.258045	26.40573	407.3928

Note that the results are very similar, but not identical, to those of the BCH model with fixed variances:

tab_auxmean <- table_results(aux_mean, columns = NULL)
tab_auxmean[tab_auxmean$Category == "Means", c("est", "se")]

	est	se
1	155.25	3.98
3	159.52	4.29
5	147.24	4.89
7	167.01	5.82

The pooled coefficients summarize the estimated relationship between pseudo-class membership and the external outcome while incorporating the uncertainty produced by repeatedly drawing class membership.

The interpretation of regression coefficients depends on how the generated class variable is parameterized in the analysis. See ?pseudo_class for more examples.

3.5.1 Using OpenMx Models

The pseudo_class() function is flexible about how the auxiliary model is supplied. The manual documents that model may be a model expression, a function that performs the analysis on each generated dataset, or a character string that can be interpreted as a structural equation model through as_ram().

If this structural equation model does not explicitly refer to the class variable, it will be used to construct a multi-group model:

pct_mx <- pseudo_class(
  x = lpa_model,
  model = "distance ~ 1",
  data = df_aux,
  m = 20
)
pct_mx[grepl("^Mean", pct_mx$term), ]

	term	estimate	se	statistic	df	df_complete
1	Means.distance.1	155.6841	4.198924	37.07714	235.9822	505
3	Means.distance.2	158.4626	3.905128	40.57810	250.1561	505
5	Means.distance.3	149.2216	5.802345	25.71746	384.3806	505
7	Means.distance.4	165.3696	7.272674	22.73849	343.0726	505

Note that the results are, again, very close to the BCH analysis.

3.6 Generating pseudo-class draws explicitly

Sometimes you may want access to the generated pseudo-class datasets themselves. append_class_draws() generates multiple datasets containing a variable named class. Each person’s class value is sampled according to their posterior class-membership probabilities.

draws <- append_class_draws(
  lpa_model,
  data = df_aux,
  m = 20
)

head(draws)

id_dataset	freqvisit	distance	class
1	3	0.000000	1
1	4	225.103494	4
1	5	115.505252	1
1	6	104.147144	2
1	3	65.699952	1
1	5	5.915863	3

The returned object has class class_draws and can itself be supplied to workflows that support pseudo-class analyses.

A crucial practical requirement is that the rows of the auxiliary dataset must be in exactly the same order as the observations used to estimate the latent class model. The class probabilities belong to specific individuals; reordering the auxiliary data would associate probabilities with the wrong cases.

This is a good reason to prepare the indicator and auxiliary datasets together before model fitting, rather than independently sorting or filtering them later.

3.7 BCH or Pseudo-class?

Both methods are designed to avoid treating modal class membership as perfectly observed, but they do so differently.

The BCH method expresses the auxiliary analysis as a weighted multi-group model. This is particularly convenient when the external question can be represented as an OpenMx/RAM model and when you want direct class-specific estimates plus likelihood-ratio comparisons across groups. Each case contributes proportionally to all groups’ estimates.

The pseudo-class method instead creates repeated plausible class assignments and pools the resulting analyses. Its main practical advantage is flexibility: the manual allows the model argument to be an expression, function, or structural-model syntax. This can make it useful for auxiliary analyses that are naturally handled by familiar R modeling functions. Each case contributes to one group’s estimates only, which introduces bias; pooling the results across multiple sets of class draws averages out this bias, but introduces Monte Carlo error.

The choice between them should therefore be guided by the inferential question and the model you want to fit, rather than by which method gives the most favorable p-value.

3.8 Predicting class membership

So far, we have asked whether latent classes differ on auxiliary variables. We can also reverse this question: Which variables predict class membership?

In principle, a test of class differences on an auxiliary variable also provides evidence about the association between that variable and class membership. If an auxiliary variable is significantly associated with class membership, this association will be evident when tested in either direction, whether class -> auxiliary or auxiliary -> class. The statistical parameters and their interpretation differ, but the significance test is informative. The interpretation of this significant effect is determined by the assumed direction of causality. For example, if age is assumed to be a predictor of class membership, then a significant BCH analysis of class differences in age tells you that the reverse effect is also significant.

In some cases, researchers might prefer to model class membership as the outcome and age as the predictor. With tidySEM, this type of analysis can conveniently be conducted using pseudo_class(). For example, we can use the nnet package to perform a multinomial regression, predicting class membership while accounting for uncertainty in class assignment.

library(nnet)
pct_pred <- pseudo_class(x = lpa_model,
                         model = function(data){
                           multinom(class ~ freqvisit + distance, data = data)
                           },
                         data = df_aux,
                         m = 10)

# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 682.410257
final  value 682.290626 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 680.095372
final  value 679.179913 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 680.304712
final  value 680.230835 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 676.930976
final  value 675.812143 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 680.859207
final  value 680.010316 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 675.710307
final  value 675.326779 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 675.670269
final  value 675.227329 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 675.313995
final  value 673.562526 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 668.921827
final  value 668.648656 
converged
# weights:  16 (9 variable)
initial  value 704.237535 
iter  10 value 669.424039
final  value 669.097802 
converged

pct_pred

term	estimate	std.error	statistic	df	p.value
(Intercept)	0.2161202	0.6447646	0.3351924	122.34783	0.7380540
freqvisit	-0.1289636	0.1324206	-0.9738941	77.79097	0.3331278
distance	0.0009323	0.0023445	0.3976424	215.94726	0.6912866
(Intercept)	0.7375894	0.6311327	1.1686756	301.01459	0.2434591
freqvisit	-0.2081532	0.1230617	-1.6914542	382.63818	0.0915643
distance	-0.0023268	0.0025741	-0.9039239	173.36214	0.3672901
(Intercept)	0.4671871	0.7289171	0.6409331	167.12011	0.5224441
freqvisit	-0.4656306	0.1563363	-2.9783915	85.82719	0.0037665
distance	0.0031832	0.0029506	1.0788542	153.23856	0.2823480

Note that the results are a bit messy, but essentially, we get four sets of regression coefficients, predicting membership of classes 1:4. In this case, only freqvisit is a significant predictor of membership of class 4, nothing else.

An important advantage of this approach is that several predictors can be included simultaneously. Their coefficients then represent partial associations with class membership: the association of each predictor with class membership while holding the other predictors constant.

This is a real difference from testing class differences separately for each auxiliary variable via the BCH. For example, testing whether classes differ in age can establish a univariate association between age and class membership. However, a multinomial regression can address a more specific question: Does age predict class membership after accounting for sex and education?

3.9 Forecasting Class for New Cases

This LCA model was developed to help classify care providers in a clinical context, so that mental healthcare professionals can provide tailored support to those who take care of their clients. In tidySEM, it is possible to predict class membership for new data. Imagine that we administer the care compass questionnaire to a new individual. We can assign their scale scores to a data.frame, and supply it to the predict_class() function (in previous versions, we overloaded the predict() function) via the newdata argument. The result includes the individual’s most likely class, as well as posterior probabilities for all classes.

df_new <- data.frame(
  burdened = 2,
  trapped = 0.5,
  negaffect = 1.5,
  loneliness = 4
)
predict_class(lpa_model, newdata = df_new)