A Bayesian Approach to Multilevel Latent Class Analysis

A BAYESIAN APPROACH TO MULTILEVEL LATENT CLASS ANALYSIS Ole Schacht Student ID: 01911379 Promotor: Prof. Dr. Koen De Turck Co-Promotor: Prof. Dr. Bert Weijters Tutors: Prof. Dr. Frank goedertier, Mr. Joeri Van Den Bergh A dissertation submitted to Ghent University in partial fulfillment of the requirements for the degree of Master of Science in Statistical Data Analysis. Academic year: 2024-2025 The author and promoters give permission to consult this master dissertation and to copy it or parts of it for personal use. Every other use falls under the restrictions of the copyright, in particular concerning the obligation to mention explicitly the source when using results of this master dissertation. Ghent, August 25, 2025 The promotors, The author, Prof. Dr. Koen De Turck Prof. Dr. Bert Weijters Ole Schacht Acknowledgements This thesis marks the final step to obtaining my degree in Statistical Data Analysis. I would have not been able to complete this work without the help of several important people. First, I would like to thank my supervisor, Prof. Dr. Koen De Turck, for introducing me to Bayesian statistics and for his help throughout this project. I am equally thankful to my co-supervisor, Prof. Dr. Bert Weijters, who first introduced me to modeling and consumer research. I further wish to thank my tutors, Prof. Dr. Frank Goedertier and Mr. Joeri Van Den Bergh, with whom I have had the pleasure of collaborating on several research projects over the past few years alongside Prof. dr. Bert Weijters. I am also grateful to the market research company, Human8, for providing me with access to some of their data. It has helped me to improve this thesis from being purely methodological to also demonstrate the practical application on a real data set. I would also like to express my appreciation towards the department of Data-Analysis at the faculty of Psychology. They have established a work environment where continuous learning is not only supported but encouraged. In particular, I am grateful to (former) PhD colleagues Sara Dhaene, Julie De Jonckere, Stijn Debrouwere, Lara Vankelecom, and Jasper Bogaert. When things got tough, I often thought of them, who walked the same road a few years earlier. I am thankful to my PhD supervisor, Prof. Dr. Tom Loeys, and my co-supervisor, Prof. Dr. Beatrijs Moerkerke, for their flexibility and understanding. They always gave me the time I needed to complete this degree. Of course, I am also thankful to my classmates, and now (soon to be) colleague(s), Felipe Fontana Vieira and Nathan Laroy. It was by no means an easy journey but we made sure to keep it fun and helped each other when needed. Finally, I want to express my thanks to several people outside the academic setting. To my parents, whose belief in me have always been a source of motivation. To my closest friends for being there: Viktor, Basile, Jari, Andy, Rutger, Dauwe, Jilles, Jona, Fouke, Thibault, Dorien, Laurien, and Sare. And most of all, to my girlfriend Charlotte. Without your presence, I would not have been able to finish this degree during my PhD. iii Table of Contents Acknowledgments iii Abstract vii 1 Introduction 1 1.1 LatentClassAnalysis .................................. 1 1.1.1 Parameters in a Latent Class Model . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Estimation .................................... 2 1.2 Bayesian Statistics and MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Bayesian Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Monte Carlo Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Gibbs Sampling and Metropolis-Hastings . . . . . . . . . . . . . . . . . 6 1.2.3 ConjugatePriors................................. 7 1.3 NestedData........................................ 7 1.4 Objectives and Empirical Application . . . . . . . . . . . . . . . . . . . . . . . . 8 2 A Bayesian Approach to LCA 9 2.1 The Likelihood Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Bayesian Estimation by Bayes’ Theorem . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 ConjugatePriors................................. 12 Beta Density for Item Response Probabilities . . . . . . . . . . . . . . . . 12 Dirichlet Density for Class Membership Probabilities . . . . . . . . . . 12 2.2.2 Fully Specified Posterior Distribution . . . . . . . . . . . . . . . . . . . . 13 2.3 Gibbs Sampling Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Sampling Class Probabilities . . . . . . . . . . . . . . . . . . . . . . . . . 14 Sampling Item Response Probabilities . . . . . . . . . . . . . . . . . . . . 15 Sampling Class Membership . . . . . . . . . . . . . . . . . . . . . . . . . 15 A Note on Numerical Stability . . . . . . . . . . . . . . . . . . . . . . . . 16 3 A Bayesian Approach to Multilevel LCA 17 3.1 Multilevel Structure: Likelihood Function . . . . . . . . . . . . . . . . . . . . . . 17 3.2 ConjugatePriors..................................... 21 3.3 Posterior Distribution by Bayes’ Theorem . . . . . . . . . . . . . . . . . . . . . . 22 3.4 Gibbs Sampling Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 v vi 3.4.1 A Note on Numerical Stability . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4.2 Some RCode................................... 27 3.5 ShinyApp......................................... 28 4 Empirical Application 29 4.1 Data Description and Modeling Strategy . . . . . . . . . . . . . . . . . . . . . . 29 4.1.1 Data........................................ 29 4.2 ModelSelection ..................................... 30 4.3 Formal Comparison and Convergence . . . . . . . . . . . . . . . . . . . . . . . . 32 4.3.1 Posterior Means and Uncertainty . . . . . . . . . . . . . . . . . . . . . . . 33 4.3.2 Diagnostics and Label Switching . . . . . . . . . . . . . . . . . . . . . . . 36 Our Proposed Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Substantive Label Switching . . . . . . . . . . . . . . . . . . . . . . . . . 37 TracePlots .................................... 39 4.4 Substantive Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.4.1 Cluster Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.4.2 ClassDescription ................................ 40 4.5 CovariateAdjustment.................................. 41 4.5.1 Predicted Probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5 Discussion 43 References 45 Abstract Latent Class Analysis (LCA) is a popular classification tool to categorize individuals into distinct subgroups based on a set of binary features, such as yes/no survey items or the presence/absence of disease symptoms. This classification typically happens at an unobservable level and is learned from the data itself. When these latent classes are substantively meaningful, researchers are often interested in their relative size, as well as in predicting class membership from covariates, such as age or gender. Recent methodological advancements have extended LCA to account for clustering of the data, such as individuals nested within different countries. This multilevel modeling approach may lead to more interpretable class structures while also relaxing the independence assumption on individual data. The parameters of a latent class model are typically estimated using the Expectation-Maximization (EM) algorithm. In contrast, Bayesian methods offer full posterior distributions for all parameters, which quantify the uncertainty on estimated parameters while allowing to take into account prior knowledge. While Bayesian approaches for single-level LCA are well established in the literature, often using Gibbs sampling, no detailed or accessible analog exists for multilevel LCA. This gap is largely due to the analytical and computational complexity of deriving conjugate priors and implementing an efficient sampler for multiple hierarchical levels that accounts for label switching of classes and clusters. In this dissertation, we address this gap by deriving a Gibbs sampling algorithm for multilevel LCA. The sampling approach uses Dirichlet priors for cluster and class probabilities, combined with beta priors for the individual item probabilities. We describe the likelihood function, prior densities, and full conditional distributions in detail, and compare the Bayesian posterior means to their frequentist Maximum Likelihood counterparts. We provide Rcode and a user-friendly Shiny applet that allows researchers to use the sampler for their own work (https://oleschacht.shinyapps.io/shinyapp/). To illustrate the approach, we apply it to an empirical data set consisting of 13,028 consumers from 17 countries reporting on sustainable consumption behaviors. We conclude with a discussion and provide directions for future research. vii 1Introduction Abstract. This thesis develops a Bayesian approach for estimating multilevel latent class models. To do so, we first discuss some important concepts in this introductory chapter. We begin by presenting Latent Class Analysis (LCA) and describe a typical application. Estimation of LCA model parameters via the Expectation-Maximization algorithm is briefly explained. Next, we turn to the basics of Bayesian statistics and discuss how LCA models can be estimated using Markov Chain Monte Carlo algorithms. We then extend this framework to multilevel LCA, which is useful when individuals are nested within higher-level units. This nested structure comes with additional complexity, especially in the Bayesian setting. Finally, we outline the objectives of this thesis and present a running example based on data from 13,028 consumers nested within 17 countries. 1.1 Latent Class Analysis Latent Class Analysis (LCA) is a popular multivariate data analysis technique that classifies individuals into latent segments or classes based on a set of binary features [Nylund-Gibson and Choi, 2018, Sinha et al., 2020]. These features might represent typical yes/no survey items, but they can also reflect the presence or absence of other attributes, such as disease symptoms [Formann and Kohlmann, 1996]. In any case, the defining characteristic of the data analysis is that a set of binary responses is available for each individual and that a latent decomposition is performed [Goodman, 1974]. By decomposing these response patterns into a smaller number of latent classes, researchers can study similarities between individuals at a covert level [Weller et al., 2020]. For example, consider a market research firm aiming to understand consumers’ environmentally sustainable behavior. The firm collects a representative sample of consumers and presents them with a list of eco-friendly behaviors, such as “I engage in recycling” and “I reduce my meat consumption.” Based on the responses given to these items, the firm seeks to segment the market into distinctive and meaningful groups. Additionally, the firm may be interested in predicting class membership from covariates included within the same model [Lyrvall et al., 2024], or in relating these latent classes to distal (outcome) variables [Bakk and Vermunt, 2015, Clark and Muthén, 2009, Nylund-Gibson et al., 2019]. 1 8Introduction nested within higher-level units, such as students within schools [Mayworm et al., 2023], patients within treatment centers [Harrison et al., 2013], or consumers within countries [Bassi, 2023, Yang et al., 2025]. Ignoring this hierarchical structure leads to violations of the independence assumption between observations [Park and Yu, 2016], and also discards valuable information at the group-level [Lyrvall et al., 2025]. In the context of LCA, this nested structure can be explicitly modeled by grouping higher-level units into latent clusters, which gives rise to the multilevel latent class model [Vermunt, 2003]. This approach is often called non-parametric, because it models heterogeneity at the group-level without imposing parametric assumptions on random effects [Finch and French, 2013]. 1.4 Objectives and Empirical Application While multilevel LCA has been developed within a frequentist framework, it has not been discussed in detail from a Bayesian perspective. However, it deserves notice that several papers have come fairly close to the current setting (such as Di and Bandeen-Roche [2010], Lee et al. [2025], and Vidotto et al. [2018]), but they differ with respect to several important specifics. To the best of our knowledge, no study has written an accessible Bayesian estimation approach for multilevel LCA with binary responses. The main complexity is situated in the hierarchical structure of the data: the Gibbs sampling algorithm will need to be extended accordingly. Also, this approach is currently not available in free-to-use statistical software. In summary, this thesis considers the following objectives: 1. Derive a Gibbs algorithm for non-parametric multilevel latent class analysis. 2. Program this sampling algorithm in Rand provide a Shiny application. 3. Assess the convergence between the Bayesian and frequentist approach. 4. Perform a real data analysis to test and illustrate the proposed method. The data and analysis code for this thesis are available online through the following link: https://osf.io/wxcyd/. For the last objective, I have been granted access by Human8 to analyze data from a large multi-country survey on environmental sustainability which includes a total of 13,028 individuals across 17 countries. The Bayesian approach to multilevel LCA will be validated on this data set in Chapter 4. For completeness (not but part of the main objectives), it will also be shown how both individual-level and grouplevel covariates can be included in a frequentist multilevel latent class model. The multinomial regression coefficients of this analysis will be briefly discussed, as well as the predicted probabilities for a range of sensible covariate values. 2A Bayesian Approach to LCA Abstract. In this chapter we focus on the Bayesian approach to fitting latent class models with binary responses and non-nested data. This exposition necessitates a more mathematical tone, as the concepts from the introductory chapter are now formalized. Despite this slightly more technical style, we go over the respective elements in detail and start from the very beginning. We first explain how the likelihood function arises in LCA, and how this translates into the posterior distribution by using conjugate priors. Since the resulting posterior distribution is high-dimensional, we resort to Gibbs sampling as the preferred MCMC algorithm. For numerical stability, we also introduce the algorithm on the logarithmic scale. The reader should be informed that this chapter draws reference to the paper by Li et al. [2018], which provides the necessary stepping stones for deriving the multilevel Bayesian framework in the subsequent chapter. 2.1 The Likelihood Function In latent class analysis, the primary goal is to uncover unobserved (latent) classes that explain the patterns in the observed binary responses of the individuals. The result of this analysis yields two key sets of parameter estimates: (i) the class membership probabilities, and (ii) the item response probabilities conditional on class membership. In what follows, we begin by deriving the likelihood function for the observed data, which will later come into play when constructing the posterior distribution in our Bayesian framework. Suppose we have full data available from Nindividuals. These individuals are asked to respond to Kbinary items, further referenced as k=1,2, . . . , K. The set of binary responses (taking only values 0or 1) for any given individual can be represented by a vector of length K, denoted as yifor individual i=1,2, . . . , N. Specific responses are referenced as yik. In this chapter, the latent class probabilities will be denoted by vector notation π, and individual elements are indexed by πjfor classes j=1,2, ..., J. The set of item response probabilities per class will be denoted by pj. Because there are Ksuch probabilities per latent class, individual elements are indexed by pjk, with element jreferencing to the specific latent class and kto the specific item. The complete matrix can be written as p. The likelihood function of observing the vector yifor individual iis a function of the 9 10 A Bayesian Approach to LCA latent class probabilities and item response probabilities. In particular, we can write: f(yi∣π,p)=J ∑ j=1 πj K ∏ k=1 pyik jk (1−pjk)1−yik .(2.1) To make this likelihood for individual ivery specific, suppose there are in total K=3 items and we observe for individual ithe following pattern: yi=(1,1,0). Assume further that there are J=2latent classes, each having their own set of item response probabilities. Table 2.1 shows the true underlying parameter values for this hypothetical example. pj1pj2pj3 π1=0.7 0.7 0.6 0.8 π2=0.3 0.4 0.2 0.4 Table 2.1 Hypothetical latent class probabilities (π) and item-response probabilities (p). Assuming that item response probabilities are conditionally independent given class membership, these can be multiplied across all items to yield the conditional probability of observing yi. If individual ibelongs to the first class, then this probability equals 0.7×0.6× (1−0.8)=0.084. As such, the conditional chance of observing yi=(1,1,0)can be written compactly as follows: ∏3 k=1pyik 1k(1−p1k)1−yik =0.084. However, if individual ihappens to belong to class 2, then this vector of probabilities, and hence the conditional probability of observing yi, will be different. Considering that π1=0.7and π2=0.3, the overall probability of observing yiis a weighted average. From introductory statistics, we know that this average is simply the sum of respective probabilities to observe yigiven class membership. As such, we arrive at the marginal likelihood to observe yi[Li et al., 2018]. If we work this out using the entries from Table 2.1, we obtain: f(yi∣π,p)=0.7×[(0.710.30)×(0.610.40)×(0.800.21)] +0.3×[(0.410.60)×(0.210.80)×(0.400.61)] =0.059 +0.014 =0.073 (2.2) Hence, assuming that the true parameters of the data-generating process are known, the marginal probability of observing yiis 0.073. This is also called the overall likelihood because it marginalizes (sums) the likelihood of observing the vector of responses for individual iover the latent classes. If we denote the class membership one-hot vector1of individual iby ci, then we can write this conditional probability down more explicitly, as follows: P(yi=(1,1,0),ci=(1,0)∣π1, p1k)=0.059 1A one-hot vector is a vector object in which a single entry equals 1 and all other entries equal 0. A Bayesian Approach to LCA 11 In words, the probability of observing yiin class 1is equal to 0.059. The same notation can be written down for class 2. For later, it will be convenient to incorporate this latent class membership notation into the likelihood function. In particular, we can write: f(yi,ci∣πj, pjk)=J ∏ j=1 ⎡ ⎢ ⎢ ⎢ ⎣πj K ∏ k=1 pyik jk (1−pjk)1−yik )⎤ ⎥ ⎥ ⎥ ⎦ cij (2.3) For example, c1(class membership vector for person i=1) can take on a specific set of values (such as c1=(1,0)), such that the exponents in the above expression yield a likelihood based solely on π1and p1. Indeed, if we work out the equation for the element c12, and knowing that this individual belongs to latent class 1, then this second part of the product will be raised to the power of 0, which always returns 1. Finally, the joint likelihood to observe the full matrix y, assuming mutually independent observations, is given by: f(y,c∣π,p)=N ∏ i=1 J ∏ j=1 ⎡ ⎢ ⎢ ⎢ ⎣πj K ∏ k=1 pyik jk (1−pjk)1−yik )⎤ ⎥ ⎥ ⎥ ⎦ cij (2.4) This concludes the likelihood function of the data [Li et al., 2018]. The final equation of this section has to do with conditional probabilities. Assuming that class membership is known, then response yi=(1,1,0)is most likely given by a person in class 1 because π= (0.7,0.3). More generally, the ratio between the proportions and the summed likelihood defines the conditional probability of an individual belonging to that specific class: P(cij ∣yi, πj, pjk)=J ∏ j=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ πj∏K k=1pyik jk (1−pjk)(1−yik) ∑J j=1πj∏K k=1pyik jk (1−pjk)(1−yik) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ cij (2.5) 2.2 Bayesian Estimation by Bayes’ Theorem The likelihood function says something about the likelihood of the data, given the parameters. For discrete responses, the likelihood function is numerically equal to the probability of the observed data for given parameters. In contrast, Bayesians seek the probability of the parameters, given the data. From introductory statistics, we know that: P(A∣B)=P(B∣A)P(A) P(B),(2.6) popularized as Bayes’ Theorem [Cornfield, 1967]. Applied to the current problem, we have: P(πj, pjk∣yi,ci)=P(yi,ci∣πj,jk )P(πj)P(pjk) P(yi,ci).(2.7) In the above expression, the numerator is the product of the summed likelihood of individual iand the prior distributions of the parameters. To complete the Bayesian estimation approach, we need to specify the prior distribution for each parameter of the model. In the 12 A Bayesian Approach to LCA context of latent class analysis, these priors are distributions of the latent class probabilities and the item response probabilities [Qiu et al., 2023]. In principle, the researcher is free to choose any distribution, as long as the parameter space is respected. For instance, it would not make sense to use a standard normal distribution as a prior for pjk, as these parameters cannot be negative and exceed 1. We now move to a discussion of conjugate priors. 2.2.1 Conjugate Priors The Bayesian approach yields a posterior distribution over the parameters of interest by combining the likelihood function with the prior distributions on the unknown parameters [Gelman et al., 1995]. These priors represent the researcher’s beliefs or knowledge about the parameters before observing any data. A prior is said to be conjugate to the likelihood if the resulting posterior distribution belongs to the same family as the prior [Fink, 1997]. This property is computationally desirable since it allows for closed-form posterior updates. In latent class analysis, there are two sets of parameters for which we need to specify priors: (i) item response probabilities, and (ii) class membership probabilities [Li et al., 2018]. Beta Density for Item Response Probabilities Each pjk represents the probability of a binary outcome and can be interpreted as the probability of a success in a Bernoulli trial, conditional on latent class j. A natural conjugate prior for a Bernoulli likelihood is therefore the Beta distribution. It is defined on the interval [0,1]and is parameterized by two positive hyperparameters, denoted αand β: pjk ∼Beta(αjk, βjk) f(pjk)=Γ(αjk +βjk) Γ(αjk)Γ(βjk)pαjk−1 jk (1−pjk)βjk −1∝pαjk −1 jk (1−pjk)βjk −1. (2.8) Here, Γ(⋅)denotes the gamma function. The mean of this prior is α/(α+β), and the parameters αand βcan be interpreted as pseudo-counts for prior successes and failures, respectively. For instance, choosing α=2and β=8implies a prior belief that pjk is likely small, with a mean of 0.2. In the absence of prior information, a non-informative or uniform prior can be used by setting α=β=1. As such, it assigns equal weight to all values in [0,1]. In this thesis, we will assume that the item response probabilities pjk are mutually independent given class membership, and we use the same hyperparameters for each item. Dirichlet Density for Class Membership Probabilities The Beta priors deals with parameters that have only two possible values, such as individual item probabilities. The vector of class probabilities π=(π1, . . . , πJ)lies on the (J−1)- A Bayesian Approach to LCA 13 dimensional probability simplex. The conjugate prior for a multinomial or categorical distribution is the Dirichlet distribution [Di and Bandeen-Roche, 2010]. It is parameterized by a vector of positive hyperparameters u=(u1, . . . , uJ), with density: π∼Dirichlet(u) f(π)=Γ(u1+⋯+uJ) Γ(u1)⋯Γ(uJ) J ∏ j=1 πuj−1 j∝J ∏ j=1 πuj−1 j. (2.9) As with the Beta distribution, the hyperparameters ucan be interpreted as prior counts for the number of individuals expected in each class. A non-informative prior is obtained by setting all uj=1. We now have the necessary information to derive the posterior distribution for the latent class analysis model where observations are assumed independent. 2.2.2 Fully Specified Posterior Distribution The product of likelihood times prior, divided by the probability of observing the data, results in the posterior distribution, as formalized in Equation (2.7). This posterior represents the complete state of knowledge on the vector of unknown parameters θ=(π,p). However, it deserves noting that we do not actually need the denominator in Equation (2.7). That is, we can always recover it from the available information. To make this clear, recall that the total probability in Equation (2.2) was 0.073 =(0.059+0.014), identical in all cases of ci. But the conditional probability of belonging to a given class was different, which relates to conditional probability as discussed in Equation (2.5). To relate this to the posterior distribution, the denominator merely serves as a scaling constant to convert the numerator vector [0.059,0.014]into a normalized vector [0.808,0.192], so that it sums to one and represents an actual probability distribution [Li et al., 2018]. It has no effect on how the numerators are compared with one another. Consequently, the posterior distribution is said to be proportional to the numerator. If we work this out in detail, we find: P(π,p∣y,c)∝P(y,c∣π,p)×P(π)×P(p) =N ∏ i=1 J ∏ j=1 ⎡ ⎢ ⎢ ⎢ ⎣πj K ∏ k=1 pyik jk (1−pjk)(1−yik)⎤ ⎥ ⎥ ⎥ ⎦ cij ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Summed likelihood J ∏ j=1 πuj−1 j ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Dirichlet prior K ∏ k=1 pαjk −1 jk (1−pjk)βjk −1 ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Beta priors =N ∏ i=1 J ∏ j=1 ⎡ ⎢ ⎢ ⎢ ⎣πcij j K ∏ k=1 pyikcij jk (1−pjk)(1−yik)cij ⎤ ⎥ ⎥ ⎥ ⎦ J ∏ j=1 πuj−1 j K ∏ k=1 pαjk −1 jk (1−pjk)βjk −1 =N ∏ i=1 J ∏ j=1 πcij +uj−1 j ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Dirichlet posterior K ∏ k=1[pyij cij +αjk −1 jk (1−pjk)(1−yij )cij +βjk −1] ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Beta posteriors (2.10) 14 A Bayesian Approach to LCA The above equation shows, using Bayes’ Theorem, that the posterior is proportional to likelihood times prior. Moreover, because we select conjugate priors, the resulting posterior remains in the same distributional form as the likelihood derived in Equation (2.4). A reader may be confused at this point, and ask him-/herself the question why a rather complex sampling algorithm will be needed to perform posterior inference. Even though we have a known distributional form, it is not possible to directly compute posterior summaries because this density is high-dimensional (for this toy example, we already have 7 unknown parameters). Therefore, simulation methods have to be used to approximate this posterior density. The next section discusses a Gibbs sampling algorithm to do exactly this. We again note that this section is quite similar to Li et al. [2018]. We use it do gently work our way up to the multilevel setting, which has not yet been discussed in full detail. 2.3 Gibbs Sampling Algorithm In this section, we briefly discuss a Gibbs sampling algorithm that can be used for sampling from the posterior distribution of the parameters in ordinary latent class analysis. A Gibbs sampler iteratively draws samples from the full conditional distributions of each parameter by conditioning on the current values of the others. For our posterior density, we need to sample successively: (i) class membership probabilities, (ii) item response probabilities, and (iii) latent class membership indicators. We will now go over the respective steps in detail. Sampling Class Probabilities To sample π, we need to condition on all other parameter values and draw from the Dirichlet posterior distribution. To make this more intuitive, consider a single individual iwith observed responses yi=(1,1,0)and latent class indicator ci=(1,0). Suppose the starting values of the item probabilities pjk are as given in Table 2.1, and assuming prior parameters u=(1,1). Since ciis assumed known in this step (or through random starting values), the conditional distribution of π(assuming only a single individual in our sample) becomes: π∣rest ∼Dirichlet(u1+ci1−1, u2+ci2−1)=Dirichlet(1,0)(2.11) Another individual with ci=(0,1)would yield Dirichlet(0,1). But evidently, because πreflects the population-level class proportions, we need to pool across all individuals to obtain the full conditional posterior (this step is explained in more detail by Li et al. [2018]): π∣rest ∼Dirichlet ⎛ ⎝ N ∑ i=1 ci1+u1−1, . . . , N ∑ i=1 ciJ +uJ−1⎞ ⎠(2.12) A Bayesian Approach to LCA 15 The result of this step will be a vector of length Jthat contains class membership probabilities that sum to one. Conditioning on all parameter values that are treated as constant in iteration t, then the newly sampled latent class probabilities are denoted by π(t+1). In the next step, these newly sampled values of πare used in the sampling of the other parameters. As such, in each iteration of the algorithm, the current values are treated as known and new samples are drawn from the full conditional distributions. This process is repeated many times to generate enough samples to approximate the full posterior distribution. Because some readers may be unfamiliar with Gibbs sampling, we will proceed by explaining the next steps of the algorithm. Sampling Item Response Probabilities Assume we have drawn π(t+1)=(0.430,0.570)based on, say, 100 observations and we now move on to update the item response probabilities. Given the current values of cand y, we sample each pjk from their full conditional posterior. Considering an individual with yi=(1,1,0)and ci=(1,0), its relevant contribution to the likelihood for class 1 becomes: K ∏ k=1 pyik 1k(1−p1k)1−yik =p1 11 ×p1 12 ×(1−p13)1(2.13) For class 2, this entire term is raised to the power 0, and thus contributes 1. Since we ignore scaling constants in this step (e.g., πj), we focus solely on the binomial structure of the likelihood. The newly sampled pjk, conditional on c(t), can be written as follows: p(t+1) jk ∣rest ∼Beta ⎛ ⎝ N ∑ i=1 yikc(t) ij +αjk −1, N ∑ i=1(1−yik)c(t) ij +βjk −1⎞ ⎠(2.14) Sampling Class Membership Assume that, in a next step, the newly sampled item response probabilities are as follows: p(t+1)=⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0.743 0.431 0.634 0.213 0.536 0.862 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Together with π(t+1), these can now be used to update our latent class membership vector c(t+1). To do so, we calculate the posterior probability that individual ibelongs to each latent class j=1, . . . , J (by using Equation 2.5), given their response vector yi, the updated class probabilities π(t+1), and the updated item probabilities p(t+1): p(c(t+1) ij =1∣rest)=π(t+1) j∏K k=1p(t+1)yik jk (1−p(t+1) jk )1−yik ∑J j=1π(t+1) j∏K k=1p(t+1)yik jk (1−p(t+1) jk )1−yik (2.15) 16 A Bayesian Approach to LCA These normalized class probabilities now define a multinomial distribution from which we can sample a one-hot vector per individual. For example, if the conditional class probabilities for individual iare (0.908,0.092), we can simply draw the class indicator vector in Rusing the command: t(rmultinom(1, size = 1, prob = c(0.908, 0.092))). Now, the updated class membership vectors serve as the new values in a next iteration, and the process repeats until convergence. Once the iteration loop begins, we only need to keep track of the values in the current iteration (t), the previous iteration, and the latest iteration (t+1). This concludes all the necessary mathematical steps in Bayesian LCA. A Note on Numerical Stability To increase the efficiency of the sampler, we can apply the logarithmic transformation to the joint likelihood function [Demmel, 1984]. This will be especially important for the multilevel Gibbs sampler in Chapter 3. The corresponding posterior distribution then becomes: log P(π,p∣y,c)=C+log P(y,c∣π,p)+log P(π)+log P(p) =C+N ∑ i=1 J ∑ j=1 cij ⎡ ⎢ ⎢ ⎢ ⎣log πj+K ∑ k=1(yik log pjk +(1−yik)log(1−pjk))⎤ ⎥ ⎥ ⎥ ⎦ ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-likelihood +J ∑ j=1(uj−1)log πj ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-Dirichlet prior +J ∑ j=1 K ∑ k=1[(αjk −1)log pjk +(βjk −1)log(1−pjk)] ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-Beta priors =C+J ∑ j=1 ⎛ ⎝ N ∑ i=1 cij +uj−1⎞ ⎠log πj +J ∑ j=1 K ∑ k=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎛ ⎝ N ∑ i=1 yikcij +αjk −1⎞ ⎠log pjk +⎛ ⎝ N ∑ i=1(1−yik)cij +βjk −1⎞ ⎠log(1−pjk)⎤ ⎥ ⎥ ⎥ ⎥ ⎦, (2.16) where Cis now an additive constant in stead of a proportionality constant. The Gibbs sampling algorithm on the logarithmic scale follows accordingly. To save space in the main text, we do not demonstrate how this Gibbs sampler can be implemented in R. The corresponding OSF repository (https://osf.io/wxcyd/) for this thesis contains a function with full annotation that implements this sampler. In general, it is quite similar to the Rcode provided by Li et al. [2018], but we made their sampler somewhat more efficient, reduced unnecessary storage of parameter chains, and implemented the batch means method. There are also several user-friendly packages in Rthat support Bayesian latent class analysis, such as those written by White and Murphy [2014] and Plummer [2003]. 3A Bayesian Approach to Multilevel LCA Abstract. This third chapter develops a Gibbs sampler for multilevel latent class models that accounts for nested binary data. Although we are definitely not the first to describe Bayesian estimation methods for multilevel mixture models with categorical data (such as Vidotto et al. [2018] and Lee et al. [2025]), we believe that currently no paper describes our setting. We begin by revising the likelihood function to accommodate for group-level dependencies and demonstrate how group-level class membership can be inferred jointly with individual-level class membership. This multilevel approach leads to multiple latent class probability distributions. The algorithm is first presented formulaic but we also illustrate this step-by-step and provide computer code. To facilitate application, we additionally introduce a Shiny app that implements the full procedure. 3.1 Multilevel Structure: Likelihood Function In this chapter, we relax the assumption of mutual independence between individuals. Instead, we assume that each person belongs exclusively to one of Gdifferent groups. These groups may represent, for example, countries, organizations, schools, or neighborhoods. This multilevel structure implies that individuals (level 1) are nested within groups (level 2), and a latent class decomposition is performed at both levels of analysis [Vermunt, 2003]. In the current setup, each group is assumed to belong to one of Llatent clusters. Although clusters are occasionally referred to as group-level classes in the literature, we will avoid this terminology whenever possible. That is, we distinctively use classes (level 1) and clusters (level 2) as the standard terminology. Each cluster will share a common set of latent classes, but their frequency distribution may vary across these clusters [Lyrvall et al., 2025]. Figure 3.1 provides a visual depiction of this hierarchical structure. At the top of the figure are the observed responses given to Kitems, indexed by k=1, . . . , K. These items are linked to a set of Jlatent classes, represented as circles, and constitute the lower level of analysis. At the bottom of the graphic, we introduce Llatent clusters to which groups may belong. These clusters represent latent classes at the group-level, with each group having a probability distribution of belonging these clusters. For completeness, Figure 3.1 also includes covariates: X1and X2at the individual-level, as well as Z1and Z2at the group-level. These predictor variables can be included to model class and cluster member17 24 A Bayesian Approach to Multilevel LCA 3.4.1 A Note on Numerical Stability As discussed, working on the logarithmic scale when evaluating likelihoods and posteriors is recommended to avoid arithmetic underflow when probabilities become extremely small [Lee et al., 2025]. We now follow the same approach as in Chapter 2 to transform the posterior to the logarithmic scale, with corresponding terms grouped by parameter. log P(ψ,π,p∣y,c,w)=C+log P(y,c,w∣ψ,π,p)+log P(ψ)+log P(π)+log P(p) =C+G ∑ g=1 L ∑ ℓ=1 wgℓ log ψℓ+G ∑ g=1∑ i∈Ng J ∑ j=1 L ∑ ℓ=1 cij wgℓ log πj∣ℓ+G ∑ g=1∑ i∈Ng J ∑ j=1 cij K ∑ k=1[yik log pjk +(1−yik)log(1−pjk)] ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-likelihood +L ∑ ℓ=1(uℓ−1)log ψℓ ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-prior for ψ +L ∑ ℓ=1 J ∑ j=1(vjℓ −1)log πj∣ℓ ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-priors for πℓ +J ∑ j=1 K ∑ k=1[(αjk −1)log pjk +(βjk −1)log(1−pjk)] ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Log-priors for p =C+L ∑ ℓ=1 ⎛ ⎝ G ∑ g=1 wgℓ +uℓ−1⎞ ⎠log ψℓ ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Dirichlet log-posterior for ψ +L ∑ ℓ=1 J ∑ j=1 ⎛ ⎝ G ∑ g=1∑ i∈Ng cij wgℓ +vjℓ −1⎞ ⎠log πj∣ℓ ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Dirichlet log-posteriors for πℓ +J ∑ j=1 K ∑ k=1 ⎛ ⎝ G ∑ g=1∑ i∈Ng cij yik +αjk −1⎞ ⎠log pjk +⎛ ⎝ G ∑ g=1∑ i∈Ng cij (1−yik)+βjk −1⎞ ⎠log(1−pjk) ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ Beta log-posteriors for p In the above expression, terms that involve the same parameter are grouped in an intermediary step but this is omitted for brevity. Also, Cagain serves as an additive constant. As can be observed, we transform products of many small probabilities into sums of logprobabilities which are numerically more stable to work with. We now discuss a tailored Gibbs sampler that generates posterior samples for all unknown parameters of the model. Initialization (t=0)In the initialization step, the researcher specifies starting values from the respective prior distributions. Specifically, one draws the following values: • Cluster, class, and item response probabilities: ψ(0)∼Dirichlet(u1, . . . , uL),π(0) .∣ℓ∼Dirichlet(v1ℓ, . . . , vJℓ), p(0) jk ∼Beta(αjk, βjk) • Latent membership indicators for each group and individual: w(0) g=(w(0) g1, . . . , w(0) gL ), L ∑ ℓ=1 w(0) gℓ =1, w(0) gℓ ∈{0,1} c(0) i=(c(0) i1, . . . , c(0) iJ ), J ∑ j=1 c(0) ij =1, c(0) ij ∈{0,1} where g=1, . . . , G, and i∈Ng. A Bayesian Approach to Multilevel LCA 25 Iterative Sampling (t=1,2,...,T)At each successive iteration, the Gibbs sampler updates parameters sequentially by sampling from their full conditional distributions. In particular, the algorithmic steps are as follows: 1. Sample cluster membership w(t+1) g: For each group g, the posterior probability of membership in cluster ℓ(on the logarithmic scale) is updated based on the current cluster probability ψ(t) ℓ, the current conditional class probabilities π(t) j∣ℓ, and the current class membership vectors c(t) ij : log P(w(t+1) gℓ =1∣⋅)=C+log ψ(t) ℓ+∑ i∈Ng log ⎛ ⎝ J ∑ j=1 π(t) j∣ℓ K ∏ k=1[p(t) jk yik (1−p(t) jk )1−yik ]⎞ ⎠. (3.16) These log-transformed probabilities (up to a additive constant) are then normalized using the log-sum-exp trick1to improve numerical stability: P(w(t+1) gℓ =1∣⋅)=exp (log ψ(t) ℓ+∑i∈Nglog (∑J j=1π(t) j∣ℓ∏K k=1[p(t) jk yik (1−p(t) jk )1−yik ])) ∑L ℓ=1exp (log ψ(t) ℓ+∑i∈Nglog (∑J j=1π(t) j∣ℓ∏K k=1[p(t) jk yik (1−p(t) jk )1−yik ])). We then sample w(t+1) gfrom the resulting multinomial distribution over clusters: w(t+1) g∼Multinom (1; P(w(t+1) g1=1∣⋅), . . . , P(w(t+1) gL =1∣⋅)).(3.17) 2. Sample latent class membership c(t+1) i: For each individual in group g, conditional on cluster membership w(t+1) gℓ =1, we compute the log-probability (up to a proportionality constant) of belonging to latent class j, using the current latent class probabilities π(t) j∣ℓand item response probabilities p(t) jk : log P(c(t+1) ij =1∣⋅)=C+log π(t) j∣ℓ+K ∑ k=1[yik log p(t) jk +(1−yik)log(1−p(t) jk )].(3.18) These log-probabilities are efficiently normalized using the same trick as before: P(c(t+1) ij =1∣⋅)=exp (log π(t) j∣ℓ+∑K k=1[yik log p(t) jk +(1−yik)log(1−p(t) jk )]) ∑J j=1exp (log π(t) h∣ℓ+∑K k=1[yik log p(t) jk +(1−yik)log(1−p(t) jk )]). (3.19) The updated class membership vector c(t+1) iis then sampled accordingly: c(t+1) i∼Multinom (1; P(c(t+1) i1=1∣⋅), . . . , P(c(t+1) iJ =1∣⋅)) (3.20) 1In Bayesian statistics, the log-sum-exp trick is a clever way to normalize a vector of log probabilities by shifting the values to avoid extremely large or small numbers during exponentiation. See Blanchard et al. [2020], or https://gregorygundersen.com/blog/2020/02/09/log-sum-exp/ for details. 26 A Bayesian Approach to Multilevel LCA 3. Sample cluster probabilities ψ(t+1): Given the cluster assignments w(t+1), we update the cluster probability vector ψby sampling from its posterior Dirichlet distribution: ψ(t+1)∼Dirichlet ⎛ ⎝ G ∑ g=1 w(t+1) g1+u1−1, . . . , G ∑ g=1 w(t+1) gL +uL−1⎞ ⎠.(3.21) 4. Sample class probabilities π(t+1): For each cluster ℓ, the algorithm updates the class membership probabilities using the counts of individuals assigned to each class within that cluster: π(t+1) .∣ℓ∼Dirichlet ⎛ ⎝ G ∑ g=1∑ i∈Ng c(t+1) i1w(t+1) gℓ +v1ℓ−1, . . . , G ∑ g=1∑ i∈Ng c(t+1) iJ w(t+1) gℓ +vJℓ −1⎞ ⎠. (3.22) 5. Sample item response probabilities p(t+1): Finally, for each latent class jand item k, the algorithm samples the item response probabilities from the corresponding Beta distributions that are updated by the observed data and weighted by class membership indicators: p(t+1) jk ∼Beta ⎛ ⎝ G ∑ g=1∑ i∈Ng c(t+1) ij yik +αjk −1, G ∑ g=1∑ i∈Ng c(t+1) ij (1−yik)+βjk −1⎞ ⎠. (3.23) 6. Prepare for the next iteration: In the final step, the algorithm reassigns the old values with the current values: ψ(t)←ψ(t+1),π(t)←π(t+1),p(t)←p(t+1),c(t)←c(t+1),w(t)←w(t+1) Samples from the burn-in period are typically discarded because the algorithm is then still influenced by the starting values and the posterior has not yet reached its stationary distribution. The duration of burn-in can be based on trace plots of the posterior samples as they arrive [Lee et al., 2025]. A common heuristic is to discard the first 1000 iterations, but complex (mixture) models may require a longer burn-in period [Lu et al., 2011]. A second point of discussion concerns the estimation of Monte Carlo variance. Because the Gibbs sampler generates samples sequentially, each sample depends on the previous one, and consecutive draws are auto-correlated. As a result, the variance will be underestimated when all samples are treated as independent. In this thesis, we use the batch means method [Flegal and Jones, 2010], which divides the Tsamples into T∗=T/Mbatches of size Meach. These averages behave as almost-independent samples when Mis chosen large enough. A final point of discussion is related to the label switching problem [Stephens, 2000, Jasra et al., 2005]. We discuss this issue in Chapter 4 and propose a workable solution. A Bayesian Approach to Multilevel LCA 27 3.4.2 Some RCode The algorithm is now ready to be put into a computer program. For brevity, we only provide the actual sampling step in the following Rlisting. The complete function with preand post-processing steps (such as computing the posterior means and addressing label switching) is omitted but can be accessed on the OSF page: https://osf.io/wxcyd/. 1while (iter <= n_iter) { 2log_p_jk <- log(p_jk); log_1mp_jk <- log(1 - p_jk) 3 4# Sample cluster membership indicators w_{gl} 5for (g in 1:G) { 6indices <- group_map[[as.character(groups[g])]] 7y_g<- y[indices, , drop = FALSE] 8 9log_lik_j<- y_g %*% t(log_p_jk) + (1 - y_g) %*% t(log_1mp_jk) 10 log_probs_l<- numeric(L) 11 12 for (l in 1:L) { 13 log_lik_jl <- sweep(log_lik_j, 2, log(pi_jl[l, ]), "+") 14 log_probs_l[l] <- log(psi[l]) + sum(rowLogSumExps(log_lik_jl)) 15 } 16 17 # Normalize probabilities 18 probs <- exp(log_probs_l-max(log_probs_l)) 19 probs <- probs / sum(probs) 20 w_gl[g, ] <- rmultinom(1, 1, probs) # Sample w_{gl} 21 } 22 23 # Sample class membership indicators c_{ij} 24 log_lik_c <- y %*% t(log_p_jk) + (1 - y) %*% t(log_1mp_jk) 25 26 for (i in 1:N) { 27 g<- which(groups == group[i]) # Group of individual i 28 l<- which(w_gl[g, ] == 1) # Cluster assignment 29 log_lik_c[i, ] <- log_lik_c[i,]+log(pi_jl[l, ]) 30 } 31 32 c_probs <- exp(log_lik_c - rowMaxs(log_lik_c)) 33 c_probs <- c_probs /rowSums(c_probs) 34 c_ij <- t(apply(c_probs, 1, function(prob) rmultinom(1, 1, prob))) 35 36 # Sample cluster probabilities 28 A Bayesian Approach to Multilevel LCA 37 psi <- rdirichlet(1, colSums(w_gl) + prior_psi) 38 39 # Sample class probabilities 40 for (l in 1:L) { 41 g_in_l<- which(w_gl[, l] == 1) # Groups in cluster l 42 i_in_l<- unlist(group_map[as.character(groups[g_in_l])]) 43 pi_jl[l, ] <- rdirichlet(1, colSums(c_ij[i_in_l,,drop = F]) + prior_pi[l, ]) 44 } 45 46 # Sample item response probabilities p_{jk} 47 p_jk <- t(apply(c_ij, 2, function(z) 48 rbeta(K, alpha + colSums(z *y), beta + colSums(z *(1 - y))))) Listing 3.1 Multilevel Gibbs sampler 3.5 Shiny App To facilitate application, we also developed a user-friendly wrapper in the form of a Shiny App such that applied researchers can easily use the sampler. The applet is available online at https://oleschacht.shinyapps.io/shinyapp/ or offline through OSF. The Shiny app features a left sidebar to configure the algorithm and a main panel with two tabs to view the results. After uploading a data set and selecting the suitable variables for analysis, the sidebar includes input fields for specifying the number of latent classes and clusters, the number of Gibbs sampling iterations, burn-in period, and the size of the batch means. A “Run Bayesian mLCA” button then triggers the analysis. While running, the user also gets an indication of how long the sampler is estimated to continue based on the current iteration. The main screen features an ’about section’, and underneath that, the user finds two main panels: tables and plots. Upon completion, these panels are loaded with relevant posterior means which can be downloaded as separate .RData files. One specific plot that deserved further explanation is the dynamic panel plot. It is an interactive graph that plots the posterior means of item response probabilities for each latent class (similar to Figure 1.1 from Chapter 1). Users can dynamically group similar items into faceted panels, which is useful when the number of items is large or when items can be grouped meaningfully (e.g., by content domain, response scale, or difficulty). 4Empirical Application Abstract. In this chapter we present an empirical application to showcase the Bayesian approach to multilevel latent class analysis. It helps to make the technical concepts from the previous chapters more tangible. The analysis also serves to show that data is often far from ideal and that important choices have to be made when selecting the number of classes and clusters. Another goal of this chapter is to evaluate the convergence of Bayesian posterior means and frequentist MLEs. As such, it serves as a formal validation for our sampler. Finally, while a Bayesian treatment of covariate adjustment is beyond the scope of this thesis, we include a brief section that demonstrates how covariates at both levels can be included into the model within a purely frequentist framework. 4.1 Data Description and Modeling Strategy 4.1.1 Data Data were collected in 2024 by global market research firm Human81. The survey included responses from 13,028 individuals across 17 countries (about 800 respondents per region), and quotas were used for gender and age. A demographic breakdown is omitted for brevity. The survey included questions relating to a variety of topics. In one of these, consumers were asked about pro-environmental behaviors, and specifically about the set of behaviors they already engage in as of today: “Please indicate which of the following you are already doing” (1 = selected, 0 = not selected). The following 15 items were presented based on academic relevancy and stakeholder interest: (a) consciously buying local products; (b) consciously buying seasonal products; (c) trying to consume less in general; (d) refusing plastic bags when shopping; (e) avoiding single-use plastic items; (f) avoiding flying; (g) using a more sustainable form of transport ; (h) removing meat from my diet; (i) avoiding buying fast fashion; (j) avoiding buying water in plastic bottles; (k) looking for ethical or sustainable products; (l) sharing and renting rather than buying goods; (m) buying secondhand goods; (n) following a plant-based diet; and (o) consciously buying from inclusive brands. A final option, (p) “none of the above”, was also provided. The goal of this analysis is to perform a multilevel latent class decomposition on the sustainability responses. 1https://www.wearehuman8.com/ 29 30 Empirical Application Table 4.1 presents the estimated proportion of respondents selecting each sustainability item, disaggregated by country. Across all countries, the most frequently endorsed behavior is refusing plastic bags when shopping (Item (d); 56.4%), followed by avoiding singleuse plastics (Item (e); 50.8%). In contrast, the least commonly reported behavior is sharing or renting goods (Item (l); 12.7%), followed by removing meat from one’s diet (Item (h); 14.2%). Country-specific differences further demonstrate variation in eco-friendly behavior. For instance, German consumers report the highest prevalence on five sustainability actions among all countries surveyed. 4.2 Model Selection Arguably the most important step in multilevel latent class analysis is selecting the optimal model in terms of the number of classes and clusters. Following recommendations from the literature [Di Mari et al., 2023, Lukoˇcien˙ e et al., 2010, Lyrvall et al., 2024, Yu and Park, 2014, Lyrvall et al., 2025], we use a three-step strategy to identify the most suitable model. Step 1: Selecting the Number of Classes In the first step, we temporarily ignore the hierarchical structure of the data and estimate a series of latent class models for all individuals pooled across countries. The goal of this preliminary step is to find a baseline class solution that does not yet account for betweencountry heterogeneity. We estimate models with an increasing number of classes and evaluate them using multiple criteria: the Bayesian Information Criterion (BIC), Consistent Akaike Information Criterion (CAIC), entropy, and average proportion of classification error [Lyrvall et al., 2025]. We also consider substantive interpretation of these models. As shown in Figure 4.1, the entropy and average classification error suggest that an 8-class solution is preferred in this first step. However, the BIC and CAIC either favor a parsimonious 5-class solution or a slightly more complex 9-class specification. Given our objective of retaining a baseline classification that is as clear as possible, we decide to proceed with an 8-class solution based on the entropy and classification error. Step 2: Selecting the Number of Clusters In the second step, we fix the number of latent classes at 8 (as determined in Step 1) and introduce the nested structure by estimating multilevel models with an increasing number of country clusters. As seen in Figure 4.1, the improvement in model fit in terms of CAIC and BIC seems to level off after three clusters. A 5-cluster solution offers somewhat better fit (also in terms of class entropy) but introduces a singleton cluster (USA only). We decide to proceed with a 3-cluster solution because it is more parsimonious and interpretable. Empirical Application 31 Table 4.1 Estimated item proportions per country. Item Country Nga b c d e f g h i j k l m n o ARG 706 0.289 0.244 0.341 0.660 0.564 0.166 0.493 0.153 0.356 0.418 0.312 0.110 0.320 0.139 0.140 AUS 801 0.358 0.306 0.393 0.573 0.517 0.196 0.267 0.161 0.320 0.449 0.290 0.095 0.338 0.119 0.129 BRA 810 0.436 0.298 0.447 0.460 0.489 0.126 0.405 0.117 0.210 0.379 0.422 0.137 0.259 0.205 0.296 CHI 520 0.279 0.227 0.385 0.613 0.579 0.156 0.438 0.162 0.415 0.377 0.285 0.144 0.479 0.150 0.133 CHN 989 0.350 0.406 0.369 0.444 0.567 0.149 0.576 0.158 0.256 0.299 0.340 0.174 0.158 0.260 0.294 COL 772 0.325 0.192 0.346 0.685 0.615 0.137 0.551 0.141 0.323 0.482 0.396 0.153 0.255 0.119 0.133 FR 800 0.278 0.460 0.455 0.540 0.499 0.379 0.409 0.148 0.348 0.399 0.242 0.135 0.441 0.119 0.064 GER 805 0.482 0.491 0.417 0.694 0.458 0.426 0.422 0.219 0.389 0.396 0.191 0.088 0.325 0.205 0.144 HK 655 0.209 0.209 0.469 0.623 0.594 0.096 0.537 0.092 0.284 0.434 0.237 0.093 0.255 0.137 0.191 MEX 796 0.325 0.219 0.329 0.642 0.559 0.156 0.422 0.126 0.339 0.474 0.304 0.131 0.339 0.107 0.113 PHI 805 0.530 0.224 0.439 0.612 0.627 0.140 0.514 0.109 0.283 0.391 0.528 0.173 0.384 0.224 0.142 SAF 805 0.538 0.296 0.360 0.451 0.451 0.204 0.398 0.119 0.306 0.266 0.441 0.111 0.298 0.195 0.202 SIN 804 0.282 0.178 0.394 0.552 0.434 0.137 0.501 0.147 0.254 0.373 0.256 0.136 0.219 0.138 0.144 THA 750 0.457 0.472 0.260 0.540 0.340 0.159 0.323 0.061 0.185 0.145 0.379 0.111 0.359 0.263 0.245 UAE 609 0.368 0.278 0.443 0.555 0.534 0.125 0.397 0.189 0.261 0.282 0.394 0.194 0.258 0.256 0.215 UK 800 0.289 0.264 0.398 0.650 0.496 0.252 0.365 0.175 0.374 0.408 0.269 0.058 0.380 0.128 0.135 USA 801 0.367 0.218 0.400 0.358 0.348 0.213 0.213 0.135 0.273 0.336 0.288 0.125 0.383 0.125 0.144 TOTAL 13 028 0.366 0.297 0.390 0.564 0.508 0.192 0.426 0.142 0.302 0.371 0.329 0.127 0.316 0.171 0.170 Note. Proportion of consumers ticking off a given sustainability element (e.g., 0.90 implies that an estimated 90% of surveyed consumers in that country endorse that item); Sustainability items: (a) consciously buying local products, (b) consciously buying seasonal products, (c) trying to consume less in general, (d) refusing plastic bags when shopping, (e) avoiding single-use plastic items, (f) avoiding flying, (g) using a more sustainable form of transport (e.g., bicycle, public transport), (h) removing meat from their diet, (i) avoiding buying fast fashion (i.e., mass production, low-cost fashion), (j) avoiding buying water in plastic bottles, (k) looking for ethical or sustainable products when shopping, (l) sharing and renting rather than buying goods, (m) buying second-hand goods, (n) following a plant-based diet, (o) consciously buying from inclusive brands; ARG=Argentina, AUS=Australia, BRA=Brazil, CHI=Chile, CHN=China, COL=Colombia, FR=France, GER=Germany, HK=Hong Kong, MEX=Mexico, PHI=Philippines, SAF=South Africa, SIN=Singapore, THA=Thailand, UAE=United Arab Emirates, UK=United Kingdom, USA= United States of America 32 Empirical Application Step 3: Re-Evaluating the Number of Classes In the third step, we fix the number of clusters at 3 and re-fit a series of models with an increasing numbers of latent classes. This final step evaluates whether the optimal number of classes changes once the hierarchical structure is accounted for. Based on the elbow plots, the (C)AIC and BIC point to a 7-class solution. Adding more classes only yields marginal gains in model fit but at the expense of interpretation. In conclusion, we retain a model with 7 latent classes and 3 country clusters. It should, however, be noted that other choices would have led to different final models. Of these choices, two alternative modeling options included the following: (i) choosing five latent classes in Step 1 (based on the elbow plot), and/or choosing five country clusters in Step 2. We effectively tried all these possibilities but concluded that these alternatives did not provide as substantively meaningful results. 216000 217000 218000 219000 2 3 4 5 6 7 8 9 10 11 12 Number of Classes AIC / BIC / CAIC Step 1 : Model Fit Criteria 214000 215000 216000 217000 1 2 3 4 5 6 7 8 9 10 Number of Clusters AIC / BIC / CAIC Step 2 : Model Fit Criteria 214000 216000 218000 2 3 4 5 6 7 8 9 10 11 12 Number of Classes AIC / BIC / CAIC Step 3 : Model Fit Criteria 0.2 0.4 0.6 2 3 4 5 6 7 8 9 10 11 12 Number of Classes Entropy / Class Error Probability Step 1 : Classification Metrics 0.4 0.5 0.6 12345678910 Number of Clusters Entropy / Class Error Probability Step 2 : Classification Metrics 0.45 0.50 0.55 0.60 0.65 0.70 2 3 4 5 6 7 8 9 10 11 12 Number of Classes Entropy / Class Error Probability Step 3 : Classification Metrics AIC BIC CAIC ClassErrProb Entropy Figure 4.1 Elbow plots of relevant fit statistics per model selection step. 4.3 Formal Comparison and Convergence Having fixed the number of classes and clusters, we now turn to a Bayesian estimation of the model parameters. Specifically, we use the Gibbs sampler derived in Chapter 3 to approximate the joint posterior distribution of all unknown quantities. We run 23,000 iterations and discard the first 3,000 samples as burn-in to allow the chain to converge to its stationary distribution. To address auto-correlation in the Markov chain, we apply the batch means method with a batch size of M=100. This yields a sample size of 200 almost-independent draws for computing the posterior means. It should be emphasized that all post-burn-in samples are effectively used in the computation (i.e., the average of many local averages). Empirical Application 33 For comparison, we estimate the same model using the frequentist EM algorithm implemented in the MultilevLCA package in R[Lyrvall et al., 2025, Di Mari and Lyrvall, 2024]. This library provides maximum likelihood estimates for multilevel latent class models, with or without covariates. As discussed in Chapter 1, the EM algorithm is a deterministic method that iteratively updates parameter estimates to maximize the log-likelihood function [Dempster et al., 1977]. When starting values are initialized, its trajectory is fully determined. In contrast, the Gibbs sampler is stochastic by nature: at each iteration, new samples are drawn randomly from the conditional distributions [Casella and George, 1992]. 4.3.1 Posterior Means and Uncertainty We now compare Bayesian posterior means with their frequentist counterparts to evaluate convergence. Figure 4.2 summarizes this comparison. Panel A displays the posterior means of the item response probabilities, color-coded per latent class, and grouped across five thematic categories. These panels are purely for visualization and were not part of the estimation procedure. Line-type indicates the estimation method: solid lines are used for plotting frequentist MLEs and dashed lines for Bayesian posterior means. The comparison clearly points to strong agreement between both methods after 23,000 iterations of the sampler. Minor discrepancies are observed only for the blue-colored segment (vegetarians; discussed below), specifically for items related to plant-based dieting. In this case, the Bayesian posterior means tend to be slightly lower than their frequentist analogs. For completeness, Table 4.2 provides a numerical comparison of all parameter estimates. Avoid plastic Consume less Ethical consumption Plant−based diet Transport item_d item_e item_j item_c item_i item_l item_m item_a item_b item_k item_o item_h item_n item_f item_g 0.00 0.25 0.50 0.75 1.00 Probability Dark greens Conscious shoppers Vegetarians Minimalists Plastic avoiders Minimally engaged consumers Disengaged consumers Bayesian Frequentist 2% 8% 2% 45% 7% 21% 15% 2%2%4% 9% 23% 13% 45% 2% 14% 2% 12% 2% 61% 7% 2% 8% 2% 45% 7% 20% 16% 3%2% 5% 10% 23% 13% 46% 2% 14% 2% 12% 2% 60% 8% Frequentist Bayesian 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Cluster 3 Cluster 2 Cluster 1 A B Figure 4.2 Comparison of Bayesian and frequentist estimates. 40 Empirical Application Cluster 2: "Hybrid engagement economies". These are countries with high urbanization rates and increasing exposure to global sustainability issues (Argentina, Chile, Colombia, Hong Kong, Mexico, the Philippines, and Singapore). Theoretically, these countries may reflect a “hybrid” consumer culture: sustainability awareness is growing, but macroeconomic barriers (such as affordability and access) may still limit widespread adoption. The cultural diversity of this cluster also suggests that sustainability engagement may be shaped by local norms in addition to economic conditions. Cluster 3: "Post-industrial economies". This final cluster includes high-income, postindustrial, post-materialistic countries (Australia, France, Germany, the United Kingdom, and the United States). These countries often have strong environmental awareness and public discourse around sustainability. Consumers may be more likely to encounter incentives, social norms, and macro-economic support for sustainability behaviors such as recycling, sustainable transport, ethical shopping, and dietary shifts. 4.4.2 Class Description Class 1: “Dark Greens”. The first group shows very high engagement across nearly all behaviors. Their posterior probabilities of engaging in actions like avoiding fast fashion, reducing consumption in general, refusing plastic bags, and buying ethical products are all close to 1. Their behavioral pattern makes them the prototypical eco-friendly consumers in the sample, therefore we label these the “Dark Greens”. Class 2: “Conscious Shoppers”. These individuals are selectively engaged in the shopping domain and by using a more sustainable form of transport. They also report relatively high engagement with avoiding plastic, buying local and seasonal products, and looking for ethical brands. This class is focused on making sustainable choices within more traditional behaviors such as shopping, therefore labeled the “Conscious Shoppers”. Class 3: “Vegetarians”. This class is characterized by a strong commitment to vegetarianism and plant-based dieting, which clearly differentiates them from all other latent classes. Their engagement on other behaviors (e.g., ethical consumption) is moderate. Sustainability identity seems to center on dietary ethics, therefore labeled the “Vegetarians”. Class 4: “Minimalists”. Minimalists demonstrate moderate engagement across a range of behaviors; particularly in consumption-reduction and plastic avoidance. They are not as intensely committed as the Dark Greens, but they show a tendency toward doing more with less: avoiding fast fashion, refusing plastic, and reducing consumption overall. This suggests a lifestyle grounded in simplicity, therefore labeled the “Minimalists”. Class 5: “Plastic Avoiders”. This group is characterized by relatively high engagement with plastic-related behaviors: refusing bags, avoiding single-use items, and not buying Empirical Application 41 bottled water. These consumers show relatively low probabilities on most other items, therefore labeled the “Plastic Avoiders”. Class 6: “Minimally Engaged Consumers”. These individuals show low-to-moderate engagement across all items. Their highest endorsements are in refusing plastic bags and using more sustainable transport, but these are relatively weak overall. This group of consumers may be passively aware of sustainability but not strongly motivated to act, or lacking the resources to do so; therefore labeled the “Minimally Engaged consumers”. Class 7: “Disengaged Consumers”. Finally, we discuss the least engaged group overall. Posterior means for nearly all sustainability behaviors hover around or below 0.25, with particularly low scores on buying from sustainable brand, avoiding single-use plastics, or following a plant-based diet. Their profile points to a near-complete disengagement from sustainability, whether due to disinterest, lack of awareness, or competing priorities. For this reason, we label this class the “Disengaged Consumers”. 4.5 Covariate Adjustment Finally, we discuss how covariates can be included in multilevel latent class models, both at the country-level and individual-level [Lyrvall et al., 2025, Vermunt, 2003, Di Mari et al., 2023]. But importantly, we remind the reader that we now switch to a frequentist framework. This structural part is implemented as a two-level multinomial logistic regression model. At the lower level, latent class membership serves as the dependent variable, with individual-level covariates predicting class assignment within each cluster. At the higher level, the probability of cluster membership is modeled as a function of group-level predictors. More details on covariate-adjustment are discussed in [Lyrvall et al., 2024]. At the between-group level, we use secondary data for post-materialism and Gross Domestic product (PPP). The most recent post-materialism score (ranging from very materialistic (1) to very post-materialistic (4)) for each country in our sample are taken from the World Values Survey. Post-materialism is approximately normally distributed (mean = 2.17, SD = 0.34, range = 1.32–3.04). GDP (PPP) for each country is taken from the 2023 Worldometer We log-transform these values to account for extreme right-skewness in the original distribution (log-transformed mean = 10.57, SD = 0.70, range = 9.28–11.86). On the individual level, we consider gender (male as reference, 50.2% female), meancentered age (original mean = 41.61, SD = 13.98, range = 19-72), and their interaction as predictors for class membership within each of the three country clusters. Class 1 serves as the reference category across all three clusters. For brevity, we will not discuss the results on the logit scale but directly move on to the fitted probabilities for a range of sensible covariate values. The detailed results can again be found in the OSF-repository. 42 Empirical Application 4.5.1 Predicted Probabilities We now visualize the predicted class and cluster probabilities for this covariate-adjusted model for a range of covariate values. In particular, Figure 4.5, panel A, shows the fitted probabilities from the multinomial regression (separately per cluster) with class membership as the dependent variable. Panel B, then, shows the fitted probabilities from the multinomial regression with cluster membership as the dependent variable. Within the latter panel, we visualize these fitted probabilities for a range of log GDP values, separately for low (< 2.08) and high (> 2.26) post-materialism. This binary split corresponds to quartile 1 and quartile 3 of the sample, respectively. For conciseness, we only discuss the cluster-level predicted probabilities; the class-level results follow a similar interpretation. The left-hand plot in panel B shows the results for countries with low post-materialism scores. It is observed that, as log GDP increases, the predicted probability of cluster 2 (cluster 3) membership increases (decreases). This suggests that, among low post-materialism countries, log GDP primarily separates clusters 2 and 3. In the right plot of panel B (high post-materialism), the probability of cluster 1 (cluster 3) membership increases (decreases) for higher GDP levels. For cluster 2, the association seems less pronounced. Dark Greens −20 −10 0 10 20 0.01 0.02 0.03 0.04 Age (mean−centered) Predicted Probability Conscious Shoppers −20 −10 0 10 20 0.00 0.05 0.10 0.15 Age (mean−centered) Predicted Probability Vegetarians −20 −10 0 10 20 0.02 0.04 0.06 Age (mean−centered) Predicted Probability Minimalists −20 −10 0 10 20 0.0 0.1 0.2 0.3 Age (mean−centered) Predicted Probability Plastic Avoiders −20 −10 0 10 20 0.1 0.2 0.3 0.4 0.5 Age (mean−centered) Predicted Probability Minimally Engaged Consumers −20 −10 0 10 20 0.2 0.4 0.6 Age (mean−centered) Predicted Probability Disengaged Consumers −20 −10 0 10 20 0.0 0.1 0.2 0.3 0.4 0.5 Age (mean−centered) Predicted Probability Low Post−Materialism 9 10 11 12 0.0 0.2 0.4 0.6 0.8 1.0 Log GDP (PPP) Predicted Probability High Post−Materialism 9 10 11 12 0.0 0.2 0.4 0.6 0.8 1.0 Log GDP (PPP) Predicted Probability A B Cluster 3 Cluster 2 Cluster 1 Male Female Figure 4.5 Fitted probabilities for class/cluster membership as a function of covariates. 5Discussion In this thesis we worked out a Gibbs sampling algorithm that can be used for Bayesian multilevel latent class analysis. In addition, we provided a Shiny app that allows applied researchers to use this method without requiring extensive knowledge about Bayesian methods. Finally, we provided a real data analysis to showcase and validate the method. In Chapter 1, we started with an overview of LCA and its typical applications. It was explained how this modeling tool can be approached both from a Bayesian and frequentist standpoint. In doing so, we discussed the main ideas behind Bayesian updating and MCMC algorithms. In Chapter 2, we started by explaining the Bayesian approach to standard (non-nested) latent class analysis. In particular, we discussed the likelihood function in great detail and illustrated how the posterior distribution is proportional to the likelihood of the data times the prior distribution. The chapter ended by providing a Gibbs sampler which can be found on the OSF page of this project: https://osf.io/wxcyd/. As discussed many times, we followed closely the work by Li et al. [2018]. However, we decided to work on the logarithmic scale and also made the Gibbs sampler more efficient by avoiding redundant computations and storage use. It should also be noted that there are several good Rpackages that can be used for (Bayesian) latent class analysis, such as BayesLCA [White and Murphy, 2014] and JAGS [Qiu, 2022]. We decided to work out the sampler from scratch to make the transition to multilevel LCA somewhat more gradual. In Chapter 3 we extended the likelihood function to accommodate for the hierarchical structure that is often present in behavioral data. That is, individuals that provide responses to a survey are often embedded within larger entities, and individuals within entities respond more similarly than individuals between entities [Vermunt, 2003, Lukoˇcien˙ e et al., 2010]. As a result of bringing in this extra layer of information, the multilevel LCA models a set of class probabilities per cluster to which one or more groups (entities) belong [Lyrvall et al., 2025]. Following the same build-up from Chapter 2, we constructed the likelihood of the data and provided an intuitive example to make the derivations easier to follow. We then introduced the set of conjugate priors for the likelihood based on similar arguments from the non-nested case [Li et al., 2018]. Consequently, we showed how the posterior distribution (up to a proportionality constant) is again obtained by collecting exponents. After having derived the log-transformed posterior distribution (up to an additive constant), we explained how a Gibbs sampling algorithm can be constructed from the full conditional distributions [Casella and George, 1992]. As such, we showed the algorithmic 43 44 Discussion steps to estimate the parameters from a multilevel LCA within a full Bayesian framework. In the main text, we also briefly mentioned some papers that considered similar Bayesian approaches to similar settings. One of these attempts is a very recent paper on Bayesian multilevel latent class profile analysis [Lee et al., 2025], which was published after the current sampler was derived. The difference with their approach is that they consider a setting where individuals respond to the same set of items across time, while simultaneously being nested within groups. As such, they actually consider three levels of latent categorization, including a decomposition of latent change profiles. It may well be that our sampler is a special case of their work, where only a single time occasion is considered. Future research should assess whether our approach can be mapped onto the work by [Lee et al., 2025]. After having discussed the Gibbs sampler, we provided corresponding algorithmic code in the Rlanguage, as well as a user-friendly Shiny app. To save space in the main text, we only explained the heart of the algorithm; the interested reader is referred to the OSF repository for the full implementation. The Shiny App enables researchers with limited Rexperience to use the Gibbs sampler with ease. Upon completion, posterior means for all model parameters can be downloaded and plotted. Finally, in Chapter 4, we provided a data analysis from over 13,000 consumers from 17 countries to test and showcase our sampler. As a validation, all posterior means were very close to their corresponding maximum likelihood estimates after label switching was accounted for [Jasra et al., 2005, Stephens, 2000]. We visualized the uncertainty of parameter estimates by comparing 95% confidence/credible intervals for the latent class probabilities. The current thesis leaves important extensions for future work. In particular, we did not discuss covariate adjustment in a Bayesian setting. In the frequentist realm, covariate adjustment is a well-established and active area of research [Lyrvall et al., 2025]. Extending this to the Bayesian context is conceptually straightforward, but it introduces considerable computational challenges. In particular, this would require a Metropolis-Hastings-withinGibbs algorithm, which is not easy to program. We also did not explore alternative Bayesian implementation tools such as JAGS [Plummer, 2003] or Stan [Gelman et al., 2015]. In principle, these tools could the simplify the sampling process and improving efficiency. In conclusion, this thesis presented a full Bayesian approach to multilevel latent class analysis. We discussed the derivation of the likelihood, specification of conjugate priors, construction of a Gibbs sampler, and a solution to label switching. The proposed model allows researchers to handle both withinand between-group heterogeneity in binary responses. Beyond our theoretical work, we also translated this sampler into an Rfunction and a Shiny app to make the approach accessible to researchers. The empirical application provided a check for the convergence of Bayesian and frequentist parameter estimates. Bibliography J. Albert. Bayesian computation with R. Springer, New York, 2007. doi: 10.1007/ 978-0-387-92298-0. Z. Bakk and J. K. Vermunt. Robustness of stepwise latent class modeling with continuous distal outcomes. Structural Equation Modeling: A Multidisciplinary Journal, 23(1):20–31, May 2015. ISSN 1532-8007. doi: 10.1080/10705511.2014.955104. URL http://dx.doi. org/10.1080/10705511.2014.955104. F. Bassi. European consumers’ attitudes towards the environment and sustainable behavior in the market. Sustainability, 15(2):1666, Jan. 2023. ISSN 2071-1050. doi: 10.3390/su15021666. URL http://dx.doi.org/10.3390/su15021666. T. Bijmolt, L. Paas, and J. Vermunt. Country and consumer segmentation: Multi-level latent class analysis of financial product ownership. International Journal of Research in Marketing, 21(4):323–340, 2004. ISSN 0167-8116. Pagination: 17. P. Blanchard, D. J. Higham, and N. J. Higham. Accurately computing the log-sum-exp and softmax functions. IMA Journal of Numerical Analysis, 41(4):2311–2330, Aug. 2020. ISSN 1464-3642. doi: 10.1093/imanum/draa038. URL http://dx.doi.org/10.1093/ imanum/draa038. G. Casella and R. Berger. Statistical Inference. Chapman and Hall/CRC, Apr. 2024. ISBN 9781003456285. doi: 10.1201/9781003456285. URL http://dx.doi.org/10.1201/ 9781003456285. G. Casella and E. I. George. Explaining the gibbs sampler. The American Statistician, 46 (3):167–174, Aug. 1992. ISSN 1537-2731. doi: 10.1080/00031305.1992.10475878. URL http://dx.doi.org/10.1080/00031305.1992.10475878. H. Chen, L. Han, and A. Lim. Beyond the em algorithm: constrained optimization methods for latent class model. Communications in Statistics - Simulation and Computation, 51(9): 5222–5244, May 2020. ISSN 1532-4141. doi: 10.1080/03610918.2020.1764034. URL http: //dx.doi.org/10.1080/03610918.2020.1764034. S. L. Clark and B. Muthén. Relating latent class analysis results to variables not included in the analysis. Technical report, Muthén & Muthén, 2009. URL https: //www.statmodel.com/download/relatinglca.pdf. J. Cornfield. Bayes theorem. Revue de l’Institut International de Statistique / Review of the International Statistical Institute, 35(1):34, 1967. ISSN 0373-1138. doi: 10.2307/1401634. URL http://dx.doi.org/10.2307/1401634. 45 46 Discussion J. Demmel. Underflow and the reliability of numerical software. SIAM Journal on Scientific and Statistical Computing, 5(4):887–919, Dec. 1984. ISSN 2168-3417. doi: 10.1137/0905062. URL http://dx.doi.org/10.1137/0905062. A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society Series B: Statistical Methodology, 39(1):1–22, Sept. 1977. ISSN 1467-9868. doi: 10.1111/j.2517-6161.1977.tb01600.x. URL http://dx.doi.org/10.1111/j.2517-6161.1977.tb01600.x. C.-Z. Di and K. Bandeen-Roche. Multilevel latent class models with dirichlet mixing distribution. Biometrics, 67(1):86–96, Jun 2010. doi: 10.1111/j.1541-0420.2010.01448.x. R. Di Mari and J. Lyrvall. multilevLCA: Estimates and plots single-level and multilevel latent class models, 2024. [Computer software manual]. R. Di Mari, Z. Bakk, J. Oser, and J. Kuha. A two-step estimator for multilevel latent class analysis with covariates. Psychometrika, 88(4):1144–1170, Aug 2023. doi: 10.1007/s11336-023-09929-2. W. H. Finch and B. F. French. Multilevel latent class analysis: Parametric and nonparametric models. The Journal of Experimental Education, 82(3):307–333, Sept. 2013. ISSN 1940-0683. doi: 10.1080/00220973.2013.813361. URL http://dx.doi.org/10.1080/ 00220973.2013.813361. D. Fink. A compendium of conjugate priors, 1997. URL https://citeseerx.ist.psu. edu/viewdoc/summary?doi=10.1.1.157.5540. J. M. Flegal and G. L. Jones. Batch means and spectral variance estimators in markov chain monte carlo. The Annals of Statistics, 38(2), Apr. 2010. ISSN 0090-5364. doi: 10.1214/ 09-aos735. URL http://dx.doi.org/10.1214/09-AOS735. A. K. Formann and T. Kohlmann. Latent class analysis in medical research. Statistical Methods in Medical Research, 5(2):179–211, June 1996. ISSN 1477-0334. doi: 10.1177/ 096228029600500205. URL http://dx.doi.org/10.1177/096228029600500205. A. Gelman and D. B. Rubin. Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457–511, 1992. A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin. Bayesian Data Analysis. Chapman and Hall/CRC, 1995. A. Gelman, D. Lee, and J. Guo. Stan: A probabilistic programming language for bayesian inference and optimization. Journal of Educational and Behavioral Statistics, 40(5):530–543, Oct. 2015. ISSN 1935-1054. doi: 10.3102/1076998615606113. URL http://dx.doi. org/10.3102/1076998615606113. Discussion 47 S. Ghosal and A. W. van der Vaart. Entropies and rates of convergence for maximum likelihood and bayes estimation for mixtures of normal densities. The Annals of Statistics, 29(5), Oct. 2001. ISSN 0090-5364. doi: 10.1214/aos/1013203452. URL http://dx.doi. org/10.1214/aos/1013203452. L. A. Goodman. Exploratory latent structure analysis using both identifiable and unidentifiable models. Biometrika, 61(2):215–231, 1974. ISSN 1464-3510. doi: 10.1093/biomet/61. 2.215. URL http://dx.doi.org/10.1093/biomet/61.2.215. W. J. Harrison, M. S. Gilthorpe, A. Downing, and P. D. Baxter. Multilevel latent class modelling of colorectal cancer survival status at three years and socioeconomic background whilst incorporating stage of disease. International Journal of Statistics and Probability, 2 (3), July 2013. ISSN 1927-7032. doi: 10.5539/ijsp.v2n3p85. URL http://dx.doi.org/ 10.5539/ijsp.v2n3p85. A. Jasra, C. C. Holmes, and D. A. Stephens. Markov chain monte carlo methods and the label switching problem in bayesian mixture modeling. Statistical Science, 20(1), Feb. 2005. ISSN 0883-4237. doi: 10.1214/088342305000000016. URL http://dx.doi.org/ 10.1214/088342305000000016. A. Jones. Posterior consistency, 2021. URL https://andrewcharlesjones.github. io/journal/posterior-consistency.html. Technical blog post. J. Lee, D. B. McCoach, O. Harel, and H. Chung. Bayesian multilevel latent class profile analysis: Inference and estimation for exploring the diverse pathways to academic proficiency. Multivariate Behavioral Research, page 1–19, May 2025. ISSN 1532-7906. doi: 10.1080/00273171.2025.2501341. URL http://dx.doi.org/10.1080/00273171. 2025.2501341. Y. Li, J. Lord-Bessen, M. Shiyko, and R. Loeb. Bayesian latent class analysis tutorial. Multivariate Behavioral Research, 53(3):430–451, Feb 2018. doi: 10.1080/00273171.2018.1428892. W. A. Link and M. J. Eaton. On thinning of chains in mcmc. Methods in Ecology and Evolution, 3(1):112–115, June 2011. ISSN 2041-210X. doi: 10.1111/j.2041-210x.2011.00131.x. URL http://dx.doi.org/10.1111/j.2041-210X.2011.00131.x. Z. L. Lu, Z. Zhang, and G. Lubke. Bayesian inference for growth mixture models with latent class dependent missing data. Multivariate Behavioral Research, 46(4):567–597, July 2011. ISSN 1532-7906. doi: 10.1080/00273171.2011.589261. URL http://dx.doi.org/ 10.1080/00273171.2011.589261. O. Lukoˇcien˙ e, R. Varriale, and J. K. Vermunt. The simultaneous decision(s) about the number of lowerand higher-level classes in multilevel latent class analysis. Sociological Methodology, 40(1):247–283, Aug 2010. doi: 10.1111/j.1467-9531.2010.01231.x. 48 Discussion J. Lyrvall, Z. Bakk, J. Oser, and R. Di Mari. Bias-adjusted three-step multilevel latent class modeling with covariates. Structural Equation Modeling: A Multidisciplinary Journal, 31(4): 592–603, Feb 2024. doi: 10.1080/10705511.2023.2300087. J. Lyrvall, R. Di Mari, Z. Bakk, J. Oser, and J. Kuha. Multilevel latent class analysis: State-ofthe-art methodologies and their implementation in the r package multilevlca. Multivariate Behavioral Research, 60(4):731–747, Mar. 2025. ISSN 1532-7906. doi: 10.1080/00273171. 2025.2473935. URL http://dx.doi.org/10.1080/00273171.2025.2473935. A. M. Mayworm, J. D. Sharkey, and K. Nylund Gibson. An exploration of the authoritative school climate construct using multilevel latent class analysis. Contemporary School Psychology, 27:283–302, 2023. D. Mindrila. Bayesian latent class analysis: Sample size, model size, and classification precision. Mathematics, 11(12):2753, Jun 2023. doi: 10.3390/math11122753. I. J. Myung. Tutorial on maximum likelihood estimation. Journal of Mathematical Psychology, 47(1):90–100, Feb. 2003. ISSN 0022-2496. doi: 10.1016/s0022-2496(02)00028-7. URL http://dx.doi.org/10.1016/S0022-2496(02)00028-7. K. L. Nylund, T. Asparouhov, and B. O. Muthén. Deciding on the number of classes in latent class analysis and growth mixture modeling: A monte carlo simulation study. Structural Equation Modeling: A Multidisciplinary Journal, 14(4):535–569, Oct. 2007. ISSN 1532-8007. doi: 10.1080/10705510701575396. URL http://dx.doi.org/10.1080/ 10705510701575396. K. Nylund-Gibson and A. Y. Choi. Ten frequently asked questions about latent class analysis. Translational Issues in Psychological Science, 4(4):440–461, Dec 2018. doi: 10.1037/ tps0000176. K. Nylund-Gibson, R. P. Grimm, and K. E. Masyn. Prediction from latent classes: A demonstration of different approaches to include distal outcomes in mixture models. Structural Equation Modeling: A Multidisciplinary Journal, 26(6):967–985, Apr. 2019. ISSN 15328007. doi: 10.1080/10705511.2019.1590146. URL http://dx.doi.org/10.1080/ 10705511.2019.1590146. J. Park and H.-T. Yu. The impact of ignoring the level of nesting structure in nonparametric multilevel latent class models. Educational and Psychological Measurement, 76(5):824–847, Jul 2016. doi: 10.1177/0013164415618240. J. Park and H.-T. Yu. Recommendations on the sample sizes for multilevel latent class models. Educational and Psychological Measurement, 78(5):737–761, July 2017. ISSN 1552-3888. doi: 10.1177/0013164417719111. URL http://dx.doi.org/10.1177/ 0013164417719111. Discussion 49 M. Plummer. Jags: A program for analysis of bayesian graphical models using gibbs sampling. In K. Hornik, F. Leisch, and A. Zeileis, editors, Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003), page –, Mar. 2003. M. Qiu. A tutorial on bayesian latent class analysis using jags. Journal of Behavioral Data Science, 2(2), Dec. 2022. ISSN 2574-1284. doi: 10.35566/jbds/v2n2/qiu. URL http: //dx.doi.org/10.35566/jbds/v2n2/qiu. M. Qiu, S. Paganin, I. Ohn, and L. Lin. Bayesian nonparametric latent class analysis for different item types. Multivariate Behavioral Research, 58(1):156–157, Jan 2023. doi: 10. 1080/00273171.2022.2160958. C. Robert and G. Casella. Monte Carlo Statistical Methods. Springer New York, 2004. C. P. Robert. The metropolis-hastings algorithm, 2015. URL https://doi.org/10. 48550/arXiv.1504.01896. Accessed: 2025-08-14. P. Sinha, C. S. Calfee, and K. L. Delucchi. Practitioner’s guide to latent class analysis: Methodological considerations and common pitfalls. Critical Care Medicine, 49(1): e63–e79, Nov. 2020. ISSN 0090-3493. doi: 10.1097/ccm.0000000000004710. URL http: //dx.doi.org/10.1097/CCM.0000000000004710. M. Stephens. Dealing with label switching in mixture models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 62(4):795–809, Nov. 2000. ISSN 1467-9868. doi: 10.1111/1467-9868.00265. URL http://dx.doi.org/10.1111/1467-9868.00265. C. J. Van Lissa, M. Garnier-Villarreal, and D. Anadria. Recommended practices in latent class analysis using the open-source r-package tidysem. Structural Equation Modeling: A Multidisciplinary Journal, 31(3):526–534, Oct 2023. doi: 10.1080/10705511.2023.2250920. D. van Ravenzwaaij, P. Cassey, and S. D. Brown. A simple introduction to markov chain monte–carlo sampling. Psychonomic Bulletin amp; Review, 25(1):143–154, Mar. 2016. ISSN 1531-5320. doi: 10.3758/s13423-016-1015-8. URL http://dx.doi.org/10.3758/ s13423-016-1015-8. J. K. Vermunt. Multilevel latent class models. Sociological Methodology, 33(1):213–239, Aug 2003. doi: 10.1111/j.0081-1750.2003.t01-1-00131.x. D. Vidotto, J. K. Vermunt, and K. van Deun. Bayesian multilevel latent class models for the multiple imputation of nested categorical data. Journal of Educational and Behavioral Statistics, 43(5):511–539, Apr 2018. doi: 10.3102/1076998618769871. Y. Wang, E. Kim, S.-H. Joo, S. Chun, A. Alamri, P. Lee, and S. Stark. Reconsidering multilevel latent class models: Can level-2 latent classes affect item response probabilities?