Thinking about Latent Variables

Example

Let’s see an example with more digestible results. The following data contain 25 personality items regarding the Big 5 scales. There are 2800 total subjects, and five categories of questions¹.

• agreeableness: a tendency to be compassionate and cooperative rather than suspicious and antagonistic towards others
• conscientiousness: a tendency to be organized and dependable, show self-discipline, planning
• openness to experience: appreciation for art, emotion, adventure, unusual ideas, curiosity, and variety of experience.
• extraversion: energy, assertiveness, sociability and the tendency to seek stimulation in the company of others, and talkativeness.
• neuroticism: prone to physiological stress, tendency to experience unpleasant emotions easily, such as anger, anxiety, depression, and vulnerability

We’ll just look at part of this to keep things simple, namely the agreeableness and neuroticism items.

library(psych)
bfi_trim = bfi %>% 
  select(matches('A[1-5]|N[1-5]'))

So let’s do a PCA. We’ll look at two components, though remember that PCA technically returns as many components as there are variables. You’re just seeing part of the results. Note that PCA almost universally requires you to standardize the data, so that each variable has the same variance. Otherwise you’ll just get components reflecting the relative variances. In what follows the PCA is done on the correlation matrix, which amounts to the same thing.

pc = principal(bfi_trim, 2)
pc

Principal Components Analysis
Call: principal(r = bfi_trim, nfactors = 2)
Standardized loadings (pattern matrix) based upon correlation matrix
     RC1   RC2   h2   u2 com
A1  0.09 -0.50 0.25 0.75 1.1
A2  0.03  0.77 0.60 0.40 1.0
A3 -0.01  0.80 0.63 0.37 1.0
A4 -0.07  0.61 0.38 0.62 1.0
A5 -0.16  0.70 0.52 0.48 1.1
N1  0.81 -0.12 0.67 0.33 1.0
N2  0.80 -0.11 0.65 0.35 1.0
N3  0.81 -0.03 0.66 0.34 1.0
N4  0.68 -0.15 0.49 0.51 1.1
N5  0.66  0.05 0.43 0.57 1.0

                       RC1  RC2
SS loadings           2.90 2.40
Proportion Var        0.29 0.24
Cumulative Var        0.29 0.53
Proportion Explained  0.55 0.45
Cumulative Proportion 0.55 1.00

Mean item complexity =  1
Test of the hypothesis that 2 components are sufficient.

The root mean square of the residuals (RMSR) is  0.09 
 with the empirical chi square  1923.02  with prob <  0 

Fit based upon off diagonal values = 0.91

First, focus on the portion of the output where it says SS loadings . The first line is the sum of the squared loadings² for each component (in this case where we are using a correlation matrix as the basis, summing across all 10 possible components would equal the value of 10). The Proportion Var tells us how much of the overall variance the component accounts for out of all the variables (e.g. 2.9 / 10 = 0.29). The Cumulative Var tells us that 2 components make up over 53% the variance. The others are the same thing just based on the 2 retained components rather than all 10 variables (i.e. the cumulative explained variance would = 1). We can see that the second component accounts for less variance, and this would continue with additional components, where each accounts for a decreasing amount of variance.

Some explanation of the other parts of the output:

• h2: the amount of variance in the item/variable explained by the (retained) components. It is the sum of the squared loadings, a.k.a. communality. For example, the first reasoning item shares only 38% of its variance with these components.
• u2: 1 - h2
• com: A measure of complexity. A value of 1 might be seen for something that loaded on only one component, and zero otherwise (a.k.a. perfect simple structure).
• RC: these factors are called RC to denote these are rotated components. They’ve been rotated to aid interpretability of the loadings.

Loadings, also referred to as the pattern matrix, in this scenario represent the estimated correlation of an item/variable with the respective components, and provide the key way in which we interpret the factors. As an example, we can reproduce the loadings by correlating the observed variables with the estimated component scores³.

cor(bfi_trim, pc$scores, use = 'pair') %>% round(2)

     RC1   RC2
A1  0.09 -0.50
A2  0.03  0.77
A3 -0.02  0.80
A4 -0.08  0.62
A5 -0.17  0.70
N1  0.81 -0.12
N2  0.80 -0.12
N3  0.81 -0.04
N4  0.69 -0.15
N5  0.66  0.05

It can be difficult to sort it out just by looking at the values, so we’ll look at it visually. In the following plot, stronger and more positive loadings are indicated by blue (negative are indicated by red), and we can see the different variables associated with different components. Note that the smaller negative loading comes from a reverse scored item (‘I am indifferent to the feelings of others’)⁴.

Many would use what’s called a biplot, which shows the both observations and variables on the same plot, in with components/loadings on the axes, the component scores of the data, and the loadings as vectors whose length and direction specify the strength of the loading on a particular component.

The sheering of the data here is just due to the fact we have variables that only take on 6 values, and we are looking at them in terms of continuous PC scores. Biplots often take a bit to figure out just for two components, especially when one doesn’t have a clean result, and one can’t really get a sense of more than two at time. However, using a graphical model makes clear what’s going on for all retained components. Low loadings have been suppressed for additional clarity.

Interpretation is the fun, but commonly difficult, part. In this case, the variables group like we’d expect⁵. So, we have one component that represents agreeableness, and one that represents neuroticism.

It’s worth mentioning the naming fallacy at this point. Just because we associate a factor with some concept, doesn’t make it so. In addition, for PCA the goal is not interpretation, but to reduce the data while retaining as much of the variance as possible.

Example

Now let’s do a factor analysis using the same data as we did for the PCA.

fa_model = fa(bfi_trim, 2)
fa_model

Factor Analysis using method =  minres
Call: fa(r = bfi_trim, nfactors = 2)
Standardized loadings (pattern matrix) based upon correlation matrix
     MR1   MR2   h2   u2 com
A1  0.07 -0.36 0.14 0.86 1.1
A2  0.05  0.69 0.47 0.53 1.0
A3  0.03  0.76 0.56 0.44 1.0
A4 -0.05  0.47 0.24 0.76 1.0
A5 -0.12  0.60 0.39 0.61 1.1
N1  0.78 -0.03 0.61 0.39 1.0
N2  0.76 -0.02 0.58 0.42 1.0
N3  0.77  0.05 0.58 0.42 1.0
N4  0.58 -0.08 0.36 0.64 1.0
N5  0.54  0.08 0.29 0.71 1.0

                       MR1  MR2
SS loadings           2.44 1.78
Proportion Var        0.24 0.18
Cumulative Var        0.24 0.42
Proportion Explained  0.58 0.42
Cumulative Proportion 0.58 1.00

 With factor correlations of 
      MR1   MR2
MR1  1.00 -0.19
MR2 -0.19  1.00

Mean item complexity =  1
Test of the hypothesis that 2 factors are sufficient.

The degrees of freedom for the null model are  45  and the objective function was  2.82 with Chi Square of  7880.99
The degrees of freedom for the model are 26  and the objective function was  0.23 

The root mean square of the residuals (RMSR) is  0.04 
The df corrected root mean square of the residuals is  0.05 

The harmonic number of observations is  2759 with the empirical chi square  396.78  with prob <  5.7e-68 
The total number of observations was  2800  with Likelihood Chi Square =  636.27  with prob <  1.6e-117 

Tucker Lewis Index of factoring reliability =  0.865
RMSEA index =  0.092  and the 90 % confidence intervals are  0.085 0.098
BIC =  429.9
Fit based upon off diagonal values = 0.98
Measures of factor score adequacy             
                                                   MR1  MR2
Correlation of (regression) scores with factors   0.92 0.88
Multiple R square of scores with factors          0.84 0.77
Minimum correlation of possible factor scores     0.68 0.54

The results in terms of loadings are similar to the PCA. But we can see some differences. The PCA components accounted for 53% the variance, while here we only account for 42%. In addition, we have an estimate of the correlation of the latent factors, -0.19, where the PCA components are uncorrelated by definition. And finally, there are measures of ‘fit’ that come from the structural equation modeling context. See this for more detail on those measures. I will note however that BIC is provided, which allows for a way to compare different factor solutions.

“Confirmatory” Factor Analysis

Some distinguish between exploratory factor analysis and confirmatory factor analysis, but all the latter is doing is constraining some loadings/paths to be zero. Otherwise they are identical techniques, or can be made to be. For example, if we only had a single factor, the results from ‘confirmatory’ and ‘exploratory’ FA are identical if they use the same estimation techniques (e.g. maximum likelihood).

Why would you set loadings or factor correlations to zero? Theory would suggest it. For example, the agreeableness items shouldn’t be a measure of neuroticism after accounting for the correlation of the factors, so we can set things up to coincide with our substantive assumptions. In other cases, theory might be available, but it might just also be a simplifying assumption.

The following shows how to do this with the structural equation modeling package lavaan. One just has to specify each latent variable and the corresponding observed indicators that go with it.

modelCode = "
  agree =~ A1 + A2 + A3 + A4 + A5
  neurot =~ N1 + N2 + N3 + N4 + N5

  # not necessary
  # agree ~~ neurot

  # if you don't want them correlated
  # agree ~~ 0*neurot
"

cfa_model = cfa(modelCode, data=bfi_trim)
summary(cfa_model, standardized=T, rsq=T, nd=2)

Here is a visualization of the results. Note how the arrows flow from the factor to the observed variables.

In general, our results are very similar to the previous factor analysis, but we didn’t have to estimate all the cross loadings, so this approach is a simpler model. As before, the factor explains what the items have in common, and what they do not explain is partly due to measurement error, or in other words, a lack of reliability of measurement. For more on confirmatory factor analysis, take a look here.

FA

\[ X \sim \mathcal{N}(ZW' + \mu, \Psi) \] \(\mu\) are the intercepts, \(\Psi\) is a \(DxD\) covariance matrix.

Probabilistic PCA

\[\Psi = \sigma^2I\]

The following plot shows the projection of original data onto the principal component, where we have two variables and one component. By default, PCA results in an orthogonal projection, where PPCA does not require this, and shrinks estimates to the overall mean.

PCA

With PCA we are assuming a noiseless process, and constraining \(W\) to be orthogonal.

\[\sigma^2 \rightarrow 0\]

Count-based Matrix Factorization

Commonly our data regards counts or otherwise compositional data. For example, we might have counts of species seen at different sites, data regarding genetic expression, or the number of times words are seen in various documents. In this case we can use techniques that are better suited to such data, and provide a more interpretable solution.

load('data/threes_mnist.RData')

# many of the columns are useless
zero_idx = colSums(threes) == 0

library(NMF)
# nmf_50 = nmf(threes[,!zero_idx], rank = 50) # if you don't mind waiting
nmf_10 = nmf(threes[,!zero_idx], rank = 10)

# a current bug may cause these from displaying in RStudio
coefmap(nmf_10)   
basismap(nmf_10) 

t(rmultinom(1, size = 20, prob=c(.5, .2, .15, .1, .05))) # first term is most likely, 5th term is least likely

     [,1] [,2] [,3] [,4] [,5]
[1,]    9    6    2    1    2

library(lda)
data(poliblog.documents)
data(poliblog.vocab)

set.seed(8675309)

K <- 10 ## Num topics

# this will take at least a few seconds to run
result <- lda.collapsed.gibbs.sampler(poliblog.documents,
                                       K,  ## Num topics
                                       poliblog.vocab,
                                       2000,  ## Num iterations
                                       0.1,
                                       0.1) 

## Get the top words in the cluster
top.words <- top.topic.words(result$topics, 5)

SEM

Structural Equation Modeling is a very general technique that allows one to incorporate latent variables and observed variables within a single modeling framework, allowing for indirect effects and multiple outcomes. Practically any standard model you’ve conducted can be seen as a special case. We actually already did an SEM when we did the confirmatory factor analysis. Consider the following depiction of a mediation model with three latent variables, each of which would have their respective observed variables (not shown).

We can use factor analysis in the SEM context as a way to assess measurement error. Whatever variance a factor doesn’t explain is a combination of measurement error and random variation. Thus, we can use factor analysis as a means to assess a group of variables’ ability to measure the underlying construct they are supposed to measure.

This document is actually a subset of the SEM document. Check it out for lots more detail and examples.

IRT

Item response theory can be seen as latent variable modeling for categorical data, and with specific parameters with theoretical underpinnings. It is especially applied to testing environments where the items are binary correct vs. not outcomes. An underlying logistic model relates the latent variable, which can be seen as ability in the testing context, to the observed items which have unique difficulty. The conceptual latent variable model is identical to confirmatory factor analysis.

As an example, the one-parameter, a.k.a. Rasch, model (1PM) can be expressed as follows:

\[P(y=1|\theta, \delta) = \mathrm{logis}(\theta_i-\delta_j)\]

In this setting, the probability of endorsement (or getting an item correct), \(\pi_{ij}\), is a function of the difficulty of item \(j\), \(\delta_j\) above, and the latent trait (ability) of person \(i\), \(\theta_i\). In other words, it’s a specific type of logistic regression model. In the testing context, a person with more ‘ability’ relative to the item difficulty will answer correctly. In terms of the logit:

\[\mathrm{logit_{ij}} = \mathrm{log}(\frac{\pi_{ij}}{1-\pi_{ij}}) = \theta_i-\delta_j\]

The following shows an example. With the psych package a factor analysis is done on the tetrachoric correlation matrix of the binary verbal reasoning items. Compare with the standard IRT ‘two parameter’ model from ltm. While the parameterizations are different, they come to the same conclusions about ‘ability’, as can be seen in the plot. See my IRT chapter for more details.

data("ability", package = 'psych') 
verbal = ability %>% 
  data.frame() %>% 
  dplyr::select(contains('reason')) %>% 
  drop_na()

library(psych); library(ltm)
irt_fa = fa2irt(fa(verbal, fm='ml'), rho=tetrachoric(verbal), plot = F)

irt_model = ltm(verbal ~ z1)
# irt_fa$fa  # to compare to the standard fa output
irt_fa$irt

$difficulty
$difficulty[[1]]
  reason.4  reason.16  reason.17  reason.19 
-0.5649361 -0.7011017 -0.8416509 -0.4474843 


$discrimination

Loadings:
          ML1  
reason.4  0.746
reason.16 0.559
reason.17 0.894
reason.19 0.590

                 ML1
SS loadings    2.018
Proportion Var 0.505

summary(irt_model)

Call:
ltm(formula = verbal ~ z1)

Model Summary:
   log.Lik      AIC      BIC
 -3058.632 6133.263 6175.102

Coefficients:
                   value std.err   z.vals
Dffclt.reason.4  -0.6156  0.0549 -11.2151
Dffclt.reason.16 -0.9684  0.0798 -12.1367
Dffclt.reason.17 -0.7647  0.0533 -14.3422
Dffclt.reason.19 -0.6010  0.0631  -9.5229
Dscrmn.reason.4   1.8781  0.1897   9.9001
Dscrmn.reason.16  1.3847  0.1356  10.2099
Dscrmn.reason.17  2.5396  0.3141   8.0852
Dscrmn.reason.19  1.4163  0.1357  10.4370

Integration:
method: Gauss-Hermite
quadrature points: 21 

Optimization:
Convergence: 0 
max(|grad|): 0.002 
quasi-Newton: BFGS 

ICA

The latent linear model versions of PCA and factor analysis assume the observed variables are normally distributed (even standard PCA won’t work nearly as well if the data aren’t. This is not required, and ICA, independent components analysis, does not. This visualization duplicates that seen in Murphy (2012) where we have two (uniform) independent sources. We can see that the ICA correctly recovers those components.

When would you use ICA? If you believe that truly independent sources of signal underlie your data, and that those sources are not normally distributed. It is commonly applied to deal with images or sound⁸.

Multidimensional Scaling

Multidimensional scaling⁹ is a technique for converting a matrix of observation dissimilarities into a low-dimensional map that preserves those distances as well as possible.

The first step in any MDS is choosing a similarity measure. Often this will be a metric or distance, i.e.:

• Euclidean Distance for continuous variables
• Manhattan distance or Jacaard dissimilarity for binary variables
• Similarity measures can often become dissimilarity measures by inverting (x → 1/x) or subtracting from 1 (x → 1−x).

The basic idea is to come up with what we might call components or factors that minimize the difference in factor scores for two observations from the actual distance.

data(UScitiesD)  # base R data set
UScitiesMDS <- cmdscale(UScitiesD)

We can use MDS to obtain a 2-dimensional map preserving distances as well as possible.

t-SNE

t-Distributed Stochastic Neighbor Embedding (t-SNE) is a newer technique in dimension reduction, one that can be seen as a modern version of multidimensional scaling, that tries to put similar rows of data together. Similar in this case can be thought of as follows. If you had a single variable, similarity to data point \(x_i\) would be defined as points with high probability under a normal distribution with mean of the value \(x_i\). However, with several variables we’d be thinking of a multivariate normal distribution with a mean vector. t-SNE will then try to create a new coordinate system of points (call them components or factors if you like) that tries to maintain those similarities. It is particularly well suited to high-dimensional data, or at least some implementations of it are¹⁰. People have used it well for visualization purposes, but it’s no different than what people have been doing to visualize the results of MDS and PCA for decades, just that the underlying technique is different.

The t-SNE approach is more difficult to implement, as there are particular parameters to set beyond just the number of components, and they will definitely affect the visual result. For example, I could only duplicate the MDS result for the US cities data with some fiddling and luck. In addition, there is some randomness in the process, so you can get slightly different results each time. None of this really matters, as the only real stated goal with this technique is good visualization. You might also could use it as a prelude to determining the number of clusters before cluster analysis.

library(tsne)
tsne_result_us = tsne(UScitiesD, k = 2, initial_dims=6, max_iter = 2000, perplexity = 4, epoch = 250)

Here’s what it came up with for the MNIST data set.

Recommender Systems

Practically everyone has been exposed to recommender systems such as collaborative filtering and related models. That’s how Netflix, Amazon and others make their recommendations to you given the information you’ve provided about likes and dislikes, what other similar people have provided, and how similar the object of interest is to others.

The following image, taken from Wikipedia (click the image to go there), conceptually shows how a user-based collaborative filtering method would work, where a recommendation is given based on what other similar users have given.

<img src=“../img/collaborative_filtering.gif” style=“display:block; margin: 0 auto;” width:35%>

Let’s go with movies as an example. You might only rate a handful, and indeed most people will not rate most movies. But at some point most movies will have been rated. So how can one provide a recommendation for some movie you haven’t seen? If we group similar movies into genres, and similar people into demographic categories and based on taste, one can recommend something from a similar genre of movies that you like, that people in the same demographic category seem to like as well.

If you think of genres as latent categories of movies, now you can employ the techniques we’ve talked about¹¹. Similarly, we can find clusters of people using cluster analytic techniques. In short, collaborative filtering/recommender systems are using latent variable techniques to a specific type of data, e.g. ratings. More modern approaches will incorporate user and item characteristics, recommendations from other systems, and additional information. The following provides some code for you to play with, using a straightforward singular value decomposition on movie ratings, which is the same technique used in Base R’s default prcomp function for PCA. You might compare it with method = 'POPULAR'.

library(recommenderlab)
data("MovieLense")
MovieLense
barplot(table(getRatings(MovieLense)), col = '#ff5500', border='#00aaff') # rating frequencies
recommenderRegistry$get_entries(dataType = "realRatingMatrix")  # the methods available
recom_svd <- Recommender(MovieLense, method = "SVD")
recom_svd

After running it try some predictions.

# predicted ratings for two users
recom <- predict(recom_svd, MovieLense[2:3], type="ratings")
recom
as(recom, "matrix")[,1:10]


# comparison model
recom_popular <- Recommender(MovieLense, method = "POPULAR")
getModel(recom_popular)$topN
recom <- predict(recom_popular, MovieLense[2:3], type="topNList")
recom
as(recom, "list")

All in all, thinking about your data in terms of a recommendation system might not be too far-fetched, especially if you’re already considering factor analysis of some sort.

Time-based

Hidden Markov models can be used to model latent discrete states that result in a sequence of observations over time. In terms of a graphical model, we can depict it as follows:

In this sense, we have latent variable \(z\) that represents the hidden state of a system, while the outcome is what we actually observe. There are three latent states above, and the relative width of the edge reflects the transition probability of moving from one state to the other. The icons are the categories of observations we can potentially see, and there is some probability, given a latent state of seeing a particular category. Such a situation might lead to the following observed sequence of observations:

We start in state 1 where the heart and yellow smiley are most probable. Let’s say we observe ♥. The second state is most likely so we get there, and because it has some probability of staying in that state, we observe § and §. We finally get to state 3 and see ☻, where we go back to state 1, where we see ☺, jump to latent state 3 etc. And such a process continues for the length of sequence we see.

This can be seen as a mixture model/latent class situation that we’ll talk about more later. The outcomes could also be continuous, such that the latent state determines the likelihood of the observation in a manner more directly akin to the latent linear model for standard factor analysis.

Random Effects Models

For those familiar with mixed effects models that incorporate random effects, the random effects can be seen as a latent variable in the same manner as we have been discussing here. In fact, with latent growth curve models widely use in the SEM context, they are explicitly modeled as such, but it turns out these are just a different way to run a mixed model. For typical mixed models, they are essentially a latent construct that represents unknown sources of variance that aren’t captured by the model. In addition, the residual covariance structures in mixed models, and covariance functions in Gaussian process models more generally, can actually be modeled with a factor analytic approach where loadings are estimated to produce the covariances. And finally, certain IRT models can be reformulated as mixed models. Thinking about random effects as latent variables possibly opens the door for other interpretive opportunities when dealing with them in other contexts.

Bayesian

There are Bayesian versions for any of the model-based techniques we’ve discussed, essentially just adding priors for the parameters of interest. As mentioned, LDA can be seen as latent semantic analysis with priors for the topics. With standard factor analysis, one can actually estimate the factor scores directly as additional parameters in the model. In other approaches, it can require a bit of finagling to estimate the factor analytic result, as it is inherently unidentified and requires constraints.

Beyond that, Bayesian non-parametric techniques, such as Dirichlet process approaches, allow for the possibility of not having to specify the number of clusters or latent variables in advance. Instead number is allowed to grow as the data grows according to the prior, and so can be applied to online learning situations as well. See my chapter on this for some examples.

Mixture Models

A basic approach for categorical latent variable analysis from a model-based perspective could be seen as follows:

1. Posit the number of clusters you believe there to be
2. For each observation, estimate those probability of coming from either cluster assuming some underlying distribution
3. Assign observations to the most likely class (i.e. the one with the highest probability)

Typical methods are assuming an underlying normal distribution, such that each observation is the weighted combination of its likelihood under the parameters for each proposed cluster with their respective estimated parameters. For example, with two clusters, each with an estimated mean and variance, we get the likelihood for a particular observation for each cluster. We weight them by its estimated probability of membership in that cluster, which are also estimated. The parameters of interest are those of the assumed distribution (e.g. means and (co-)variances with a normal distribution).

As an example, we’ll use mclust on the old faithful data, positing 2 clusters. As such we’ll get the proportion of observations classified to each cluster, and the means and covariance matrix (not shown) for each cluster.

library(mclust)
mod = Mclust(faithful, G = 2)
## summary(mod, parameters=TRUE)  # more detail

We can see that things come out more or less as we might expect, though there is a bit more uncertainty with observations in the middle.

The nice thing about model based approaches is that if we want to compare different solutions we can use the same tools we do in regression settings. Here, mclust does a search over many models (not just differences in cluster sizes, but also things like equal variances, shape, etc.) and selects the ‘best’ model according to BIC¹². Interestingly, the result here suggests three clusters for the Old Faithful data.

mod = Mclust(faithful, verbose=F)

library(flexmix)
flexmix(eruptions ~ 1,  data=faithful, k = 2) %>% parameters()
Mclust(faithful$eruptions, G=2, verbose=F) %>% summary(parameters=T)

mod = flexmix(eruptions ~ waiting,  data=faithful, k = 2)
summary(mod)

Call:
flexmix(formula = eruptions ~ waiting, data = faithful, k = 2)

       prior size post>0 ratio
Comp.1 0.498  148    272 0.544
Comp.2 0.502  124    272 0.456

'log Lik.' -191.5743 (df=7)
AIC: 397.1485   BIC: 422.3891 

parameters(mod)    # means (Intercept) and std dev (sigma) for each group

                      Comp.1      Comp.2
coef.(Intercept) -1.73204354 -2.06336751
coef.waiting      0.07822496  0.07374856
sigma             0.34866389  0.39692190

“Cluster analysis”

Aside from mixture models, when people use the term ‘cluster analysis’ they are typically referring to distance-based methods. Given a distance matrix that informs how dissimilar observations are from one another, the methods try to create clusters of observations that are similar to one another.

Latent Class Analysis

The term latent class analysis (LCA) refers to the case where we assume a latent categorical variable, but is commonly applied to the observed data being categorical as well (usually binary). In fact, some disciplines have a name for every categorical continuous combination you can have.

We could simply call them by the names we have already, mixture models in the case of categorical latent variables, and factor analysis otherwise. The IRT demonstrated above can be seen as an example of latent trait analysis.

The following example comes from survey data from the 2000 American National Election Study. Two sets of six questions with four responses each, asking respondents’ opinions of how well various traits (moral, caring, knowledgeable, good leader, dishonest, intelligent) describe presidential candidates Al Gore and George W. Bush. The responses are measured on a four point scale from ‘Extremely well’ to ‘Not well at all’. In what follows, we are assuming three latent classes can explain these responses. We can also model the latent class as an outcome. Which we do here by likely vote cast (the respondent’s 2000 vote choice when asked).

library(poLCA)
data("election")
# ?election  # for details
f <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
           MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB) ~ VOTE3
lca_mod <- poLCA(f, election, nclass=3, verbose = F)
# lca_mod

The problem with LCA, and this is regardless of tool used, is that the output is ridiculously messy, because you get estimated proportions of response for every latent class (3) for every observed class (4), for every variable (8). So we’ll just plot the result instead.

Darker colors represent higher probability of response. In general, we have classes representing naysayers, more positive responders, and those with a more mixed response (these would also be known as low, high and moderate scores if we’d just left the latent variable continuous). The following shows the proportions for each latent class by vote type.

So the more negative responding types were also those likely to vote for Bush, while the other two are more likely to vote for Gore.

How many latent factors/classes?

Determining the most interpretable number of factors from a factor analysis is perhaps one of the greatest challenges in factor analysis. There are many solutions to this problem, none of which is uniformly the best.

W. Revelle

“Solving the number of factors problem is easy, I do it everyday before breakfast. But knowing the right solution is harder.”

H. F. Kaiser

The only rule you have to follow is to never use a scree plot to make your decision. Elbows indeed.

Anonymous

Many tools have been offered to assist with this, and if one were able to definitively say how many factors were ‘best’, you’d have heard about it by now. For FA and PCA, see the help file for the nfactors in psych (also unidim), which will provide eight or so metrics, and probably cause confusion more than clarity.

In general, if you’re using a latent variable or factor analysis approach for interpretative purposes, theory should be a strong, but not overwhelming guide. If your data suggest three are clearly better than two (e.g. via AIC comparison), then go with that, because that may just be how it is for your sample. For measurement models, some get hung up on ‘validated’ scales, and are reluctant to even explore other possibilities. But if the scale was validated on different types of people than your sample, you should probably be more exploratory. I can also mention that scales have been used for decades that were never adequately shown to be very reliable in the first place.

Interpretation of latent constructs

Interpretation comes largely from the weights/loadings or related coefficients. With PCA and common FA, rotation goes a long way to improving interpretability. I will say that if you interpret latent clusters in an ordinal fashion, you should strongly consider simply leaving the latent variable continuous. Even a mixture of two continuous latent variables (and/or their interaction) may uncover the latent classes you would have found, but leave the data in a more expressive context.

How to treat the data

In many situations, you may be able to treat the observed data as discrete or continuous. Unless there is a strong need otherwise, I would suggest leaving continuous, even for Likert type scales. The issue is that categorical outcome models are generally harder to estimate, and often without providing much interpretability benefit. PCA wouldn’t be your best option, but it wouldn’t necessarily be wrong to use it on binary data. Remember also, that even if your data is binary that doesn’t mean the latent variables have to be discrete (think IRT).

Predictive modeling

Many often use the clusters or factor scores obtained from an initial analysis for prediction. As such, they perform a two-step process, one where they obtain the cluster membership or scores, and then perform some other model using it, say standard regression. Unfortunately, this ignores the uncertainty in those estimates of the class membership or factor score. Essentially your standard errors will be a bit optimistic.

While personally I’m less worried about this, using SEM or other approaches will allow for simultaneous estimation, and should be employed if not too much trouble. Even then it might not be straightforward. For example, several suggestions have been made for doing so with latent classes in the SEM context to keep the latent class determination free of contamination by the predictive model. They generally just find a different issue arising with each, and using the ‘appropriate’ technique may mean other standard post-processing is not available. When the latent variables are continuous, there are concerns about smaller sample sizes with SEM, but you can use fsr (factor score regression) in lavaan to correct the standard errors. I would also recommend looking at the help file for the factor.scores function in psych.

There are many disciplines that will just do the two-step and not blink an eye, so you shouldn’t feel bad just using factor analysis and related as an initial pre-processing dimension reduction step.

If prediction is the primary goal of the analysis however, you should actually be wary of dimension reduction techniques. Having more columns than rows or mountains of data is commonly dealt with by using other techniques that can handle it, not by reducing the data so you can do a regression. Dimension reduction would still perhaps be useful in small data settings though, allowing you to retain information from multiple sources without making the ultimate model overly complex.

Issues with Latent Classes