Estimating Heterogeneous Treatment Effects Within Latent Class Multilevel Models: A Bayesian Approach

Abstract

This article presents a latent class model for multilevel data to identify latent subgroups and estimate heterogeneous treatment effects. Unlike sequential approaches that partition data first and then estimate average treatment effects (ATEs) within classes, we employ a Bayesian procedure to jointly estimate mixing probability, selection, and outcome models so that misclassification does not obstruct estimation of treatment effects. Simulation demonstrates that the proposed method finds the correct number of latent classes, estimates class-specific treatment effects well, and provides proper posterior standard deviations and credible intervals of ATEs. We apply this method to Trends in International Mathematics and Science Study data to investigate the effects of private science lessons on achievement scores and then find two latent classes, one with zero ATE and the other with positive ATE.

Keywords

causal inference clustered data finite mixture models latent subgroups Bayesian joint estimation double robustness TIMSS science data

Introduction

The Discrete Conceptualization of Treatment Effect Heterogeneity

The variation of treatment effects is an important consideration for planning and evaluating a treatment. There exists an increasing interest in studying and examining treatment effect heterogeneity (Imai & Strauss, 2011; Künzel et al., 2019; Nie & Wager, 2021; Pearl, 2017; Wendling et al., 2018; Xie et al., 2012). Treatment can be effective or ineffective depending on factors, such as psychological or health status. For example, a certain intervention may be helpful for individuals at moderate risk but unsuccessful or even harmful for others at severe risk. Although researchers wish to control for such variables, it is often impossible or impractical to collect all possible baseline or pretreatment covariates that are to be jointly associated with the latent profiles. Machine learning-based treatment effect estimators are increasingly prevalent for causal inference (Athey & Imbens, 2016; Hill, 2011; Imai & Ratkovic, 2013; Wager & Athey, 2018), but these methods do not account for the differences due to unmeasured or unobservable variables. Moreover, estimating fine-grained treatment effects moderated by all possible combinations of the observed covariates is generally unfeasible due to the curse of dimensionality. It can also suffer from the lack of statistical power to detect distinct subgroup-specific effects and may overfit the data and have low predictive accuracy for new samples (Loh & Kim, 2020).

Identifying latent subgroups and estimating different subgroup treatment effects help us investigate the important “what works for whom” question. For conceptual and practical reasons, researchers in the social, behavioral, and health sciences have investigated how the causal effect of a treatment on an outcome differs among individuals in qualitatively different latent subgroups (Kim et al., 2019; Shahn & Madigan, 2017; Spilt et al., 2013; Willke et al., 2012). Latent subgroups are often referred to as “latent classes” but are also called latent components or segmentations. Individuals’ probabilities of belonging to each class have been estimated using discrete latent variable approaches, such as finite mixture models, latent class (regression) models, latent profile analysis, and model-based clustering approaches (Clogg, 1995; McLachlan et al., 2019; Vermunt & Magidson, 2002).

Latent Class Models in Causal Inference

In the context of causal inference, heterogeneous treatment effects can be consequences of different joint covariate distributions across latent classes, different selection processes to choose or be assigned to a treatment, different outcome mechanisms, or any combinations of these. As an example of latent classes in causal inference, suppose a state offered a mathematics summer camp for high school students, and the effectiveness of the program was evaluated by the gain scores of the camp participants’ pre- and postassessment scores. Assume that the state-sponsored camp was offered to a limited number of students from each high school in the state, and some schools selected their high achieving students to enhance their growth, while other schools provided the opportunity to their low performing students to remediate their deficiencies in mathematics. In this scenario, the selection process is positively related to student ability for some students (Class 1) but negatively related to others (Class 2). The selection process can be known (manifest classes) or unknown (latent classes) to the organizer. If the summer camp covered advanced topics, the effect of the treatment (attending the summer camp) would be positive for Class 1 but close to zero for Class 2. On the other hand, if the camp reviewed basic content, the treatment effect will be close to zero for Class 1 but positive for Class 2.

In this example, it is likely that the joint covariate distributions are different between the two classes, but some important variables may not be collected (e.g., family socioeconomic status, parent education) or are unobservable (students’ motivation and self-efficacy). Furthermore, the propensity of being treated is positively correlated with the pretest for students in Class 1 but negatively correlated with the same variable (pretest) for students in Class 2. For the valid evaluation of the treatment effects, we need a statistical method that handles unobserved covariate differences and distinct selection processes, such as latent class models.

Latent class selection models can detect different associations between predictors and outcome by estimating a separate set of regression coefficients for each class, as shown in Kim and Steiner (2015), Kim et al. (2015), and Suk et al. (2021). In the real data analysis section of this article, we also found that substantially different selection and outcome models between two classes, including some significant coefficients in one class, have opposite signs in the other class.

The continuous versus discrete conceptualization of heterogeneity is important in choosing proper methodology for estimating treatment effects. To demonstrate the strength of latent class models in identifying different latent subgroups, Suk et al. (2021) constructed a hypothetical two-class latent structure, where students belong to one of the two latent classes where the average treatment effect (ATE) is two in the first class and zero in the second class. Suk and colleagues fit latent class models first and then used causal forests (Wager & Athey, 2018) within latent classes to estimate individual treatment effects. They also used causal forests without latent class models. The results (Suk et al., 2021, p. 325, figure 1) show that the individual treatment effects estimated by causal forests alone formed a smooth asymmetrical distribution with one mode, and causal forests masked the treatment effects of the two latent classes. In contrast, fitting latent class models first revealed a mixture of two distributions and their respective treatment effects.

As a discrete latent variable method, the distribution for latent classes is not predetermined but approximated by the estimated number of latent classes and class sizes based on a particular data set. Therefore, the nodes and weights of a latent class distribution do not depend on a parametric probability density distribution, in contrast to latent trait models and other continuous latent variable models that assume a parametric (often normal) distribution. Due to its flexibility in terms of its nodes and weights, latent class models can accommodate both continuous and discrete distributions, asymmetric distributions, highly skewed distributions, and multimodal distributions. Latent class modeling is often considered as a semiparametric approach (Heinen, 1996; Lindsay et al., 1991). We extend the standard latent class model for single-level data to clustered data in this article and refer to our method as latent class multilevel models.

Treatment Effects in Clustered Data

To formalize the treatment effects of interest, we use the Neyman–Rubin causal model (Neyman & Iwaszkiewicz, 1935; Rubin, 1974) with its potential outcome notation (Holland, 1986; Rosenbaum & Rubin, 1983; Rubin, 1974, 1978) and extension to multilevel settings (Hong & Raudenbush, 2006). Assume that there are N individuals nested within M clusters. Let $Y_{i j} (1)$ denote the potential outcome if individual i within cluster j was actually treated ( $T_{i j} = 1$ ) and $Y_{i j} (0)$ denote the potential outcome if individual i within cluster j was untreated ( $T_{i j} = 0$ ), where $i = 1, \dots, n_{j}$ in cluster $j = 1, \dots, M$ and $\sum_{j = 1}^{M} n_{j} = N$ . The observed outcome can be presented as

Y_{i j} = T_{i j} Y_{i j} (1) + (1 - T_{i j}) Y_{i j} (0),

under the stable unit treatment value assumption (SUTVA; Rubin, 1986), where the potential outcomes of each individual are not affected by others’ treatment assignments, and there is only a single version of treatment. For more details on SUTVA for multilevel settings, see Hong and Raudenbush (2006, 2013).

As the two potential outcomes $Y_{i j} (0)$ and $Y_{i j} (1)$ are never observed simultaneously, the treatment effects cannot be directly estimated without further assumptions. In general, we can estimate ATEs, the average linear contrast between two potential outcomes $τ = E [Y_{i j} (1) - Y_{i j} (0)]$ , if the pair of potential outcomes ( $Y_{i j} (0), Y_{i j} (1)$ ) is independent of treatment assignment $T_{i j}$ . Block randomized experiments or multisite randomized trials achieve this independence through the randomization of treatment assignment within blocks or sites. For observational multilevel data, we require unconfoundedness or conditional independence, also referred to as strong ignorability in the causal inference literature (Rosenbaum & Rubin, 1983; Rubin, 1978), which implies that the potential outcomes are independent of treatment assignment, given the observed vector of covariates $X_{i j}$

Unconfoundedness : Y_{i j} (1), Y_{i j} (0) ⊥ T_{i j} | X_{i j},

and the probability of each individual being assigned to a particular treatment given a set of observed covariates is strictly between 0 and 1

Positivity : 0 < e (X_{i j}) = Pr (T_{i j} = 1 | X_{i j}) < 1,

where $⊥$ denotes independence between two random variables and $e (X_{i j})$ is the propensity score (Rosenbaum & Rubin, 1983). Propensity scores are commonly estimated by logistic regression with single-level data and by random or fixed effects logistic regression with multilevel data (Fuentes et al., 2021; Leite, 2016). Propensity scores are often used in matching methods to match treated and control units or to weigh cases via inverse probability weighting.

Subgroup-Specific Treatment Effects

This article considers ATEs of unobserved subgroups or latent classes, where each subgroup may have different ATEs from other subgroups. Let $C_{i j} \in {1, \dots, K}$ denote the latent class membership of individual i in cluster j. A class-specific ATE can be viewed as CATE for individuals from latent class $C_{i j} = k$ and can be formalized as follows:

τ_{k} = E [Y_{i j} (1) - Y_{i j} (0) | C_{i j} = k] .

When CATE is defined as conditional ATE on observed covariates $X_{i j}$ (Imbens & Rubin, 2015), $τ (x) = E [Y_{i j} (1) - Y_{i j} (0) | X_{i j} = x$ ], various machine learning methods including causal forests can provide consistent and asymptotically normal estimators of $τ (x)$ under strong ignorability and SUTVA (Wager & Athey, 2018). For latent subgroup k, however, the estimand $τ (x)$ cannot be estimated by machine learning methods directly as $C_{i j}$ is unobserved. This study focuses on the investigation of heterogeneous treatment effects by latent classes. Our proposed method can classify individuals into K latent classes and estimate $τ_{k}$ simultaneously.

There is a major difference between the conditional ATEs estimated by machine learning methods and class-specific treatment effects estimated by latent class modeling. Whereas popular machine learning methods, such as causal forests (Wager & Athey, 2018), Bayesian additive regression trees (Zeldow et al., 2019) and targeted maximum likelihood estimation (Luque-Fernandez et al., 2018) can be used to describe the variability of treatment effects based on observed characteristics of the target population, latent class models can be implemented to better understand the underlying processes of selection or outcome mechanisms due to unobserved heterogeneity that may not be explained by observed variables (Kim & Steiner, 2015). Another important difference between machine learning and latent class methods for heterogeneous treatment effects is that heterogeneity is conceptualized as gradual variation by the former, while latent classes represent heterogeneity as a finite number of distinctive and often qualitatively different subpopulations. The quantitative versus qualitative nature of machine learning versus latent class approaches was illustrated in Suk et al. (2021) with hypothetical two-class data, where the true mixture distribution was approximated as a skewed unimodal distribution for the individual conditional ATE estimates using causal forests, indicating that the machine learning method did not identify the two distinct classes.

Comparisons With Prior Work

A mixture modeling method by Kim and Steiner (2015) provides an efficacious tool for investigating unobserved heterogeneity in causal effects with observational clustered data and has shown promising results in both simulations studies and applications with real data (Kim et al., 2017; Suk et al., 2021). However, there is one important caveat regarding these findings: When mixture selection models are fit to identify latent classes, classification may not be perfect and potential misclassification can be carried over to the subsequent outcome models.

To mitigate the effects of potential misclassification, this study presents a simultaneous procedure for mixing probability models, selection models, and outcome models, such that all the information in the data and all the three models can be used jointly for optimal classification of individuals and model estimation. Perhaps most importantly, misclassification is not carried over from one stage of analysis to another. We argue that this is a significant improvement in the use of latent classes, both theoretically and practically. As another benefit of the simultaneous procedure, our method does not regard latent class membership estimated from the selection model as fixed in the outcome model and provides proper estimates of variation for class-specific treatment effects by accounting for classification uncertainty. We achieved these advances by implementing a Bayesian estimation procedure where the likelihood function is defined as a product of three parts representing a mixing probability model, a selection model, and an outcome model.

Another important development distinguishing our method from other studies using latent classes for selection or outcome models is that, whereas previous work conducted classification at the cluster level (i.e., assigning clusters into one of the latent classes), the current study classifies units at the lowest level in the data. We emphasize that classification at the lowest level in multilevel data allows more flexible forms of heterogeneity and provides richer information about unobserved heterogeneity. As an example, consider a common multilevel data structure where students are nested within schools, and it is of interest to evaluate the effects of an optional online program on student learning outcomes in a distance learning environment. One can anticipate that the effects of a program vary substantially across students depending on students’ characteristics and their environments, such as technical equipment, parent and teacher support, quality of the Internet, and distractions at home, as well as students’ attributes and preference toward the particular program and distance learning itself. It would be beneficial for students, parents, teachers, educational directors, and policymakers to understand to whom and under what circumstances a program is successful. Classification at the individual level can identify a group of students who receive meaningful benefits from a treatment as well as help understand their selection processes to the treatment.

From a methodological point of view, classification at the lowest level allows the investigation of heterogeneity with weaker restrictions than classification at the cluster level. Continuing to use the distance learning program example, if classification is closely related to home environment and students’ characteristics, students within the same schools will often be divided into different latent classes if the schools consist of diverse students. If classification is largely determined by regions, district support, or school resources, class separation will occur mostly at the school level, and results will be similar to the cluster-level classification. Importantly, classification at the lowest level allows investigation of these two scenarios and scenarios in-between (e.g., some schools belong to one class, while other schools consist of students in multiple classes) without assuming particular patterns of classification at any level, in contrast to those methods that do not allow multiple class memberships within clusters. Therefore, we believe that the current research substantially contributes to the literature in the field by accounting for classification uncertainty directly in the estimation of treatment effects and also classifying units at the lowest level in the data.

Bayesian Joint Procedure for Classification and Estimation

Model Specification

Suppose there are M clusters and n_j individuals within cluster j, and each individual comes from one of K latent classes labeled as $1, 2, \dots, K$ . For individual i from cluster j, let $X_{i j}, Z_{i j}, Y_{i j}, T_{i j}$ , and $C_{i j}$ denote its covariates with fixed effects, covariates with random effects, outcome, treatment indicator, and class membership, respectively. Here, $C_{i j}$ is unobserved while other variables are observed. $X_{i j}$ includes all the covariates: Level 1, Level 2, and possibly their interactions (for slopes as outcomes models). $Z_{i j}$ is a subset of $X_{i j}$ that includes Level 1 covariates only because it specifies which Level 1 covariates’ slopes are cluster-specific. Usually, both $X_{i j}$ and $Z_{i j}$ include the intercept.

Let $θ$ be all parameters to be estimated. The probability density function of $(Y_{i j}, T_{i j})$ given covariates $X_{i j}$ is

\begin{array}{l} p (Y_{i j}, T_{i j} | X_{i j}, θ) & = \sum_{k = 1}^{K} p (Y_{i j}, T_{i j}, C_{i j} = k | X_{i j}, θ) \\ = \sum_{k = 1}^{K} p (Y_{i j} | T_{i j}, C_{i j} = k, X_{i j}, θ) Pr (T_{i j} | C_{i j} = k, X_{i j}, θ) Pr (C_{i j} = k | X_{i j}, θ) . \end{array}

That is, the probability density of each individual is a mixture of K components. Each component can be decomposed into three parts: the outcome model for $Y_{i j}$ , the selection model for $T_{i j}$ , and the mixing probability for $C_{i j}$ . They can be modeled as follows.

Selection model

This component models how the propensity score, the probability that the individual is assigned to the treatment group, is related to covariates given the latent class. Assuming $T_{i j}$ follows a logit mixed-effects model, we have

Pr (T_{i j} | C_{i j} = k, X_{i j}, θ) = \frac{exp [(X_{i j}' α_{k} + Z_{i j}' v_{k j}) T_{i j}]}{1 + exp (X_{i j}' α_{k} + Z_{i j}' v_{k j})},

where $α_{k}$ is the fixed effects of class k, and

(v_{1 j}', v_{2 j}', \dots, v_{K j}')^{'} \sim N (0, Σ_{v})

are the random effects of cluster j for all K latent classes. For each cluster, by stacking random effects across all latent classes together and specifying a multivariate normal prior, we allow random effects of different latent classes to be correlated. This is reasonable because the correlation captures cluster-level effects across latent classes. For example, we may expect that individuals within the same cluster all tend to have relatively high or low probabilities to be treated regardless of their latent classes, which implies that all the K random intercepts (corresponding to the K latent classes respectively) should be positively correlated.

Outcome model

This component models how the outcome is related to covariates when the treatment indicator and the latent class are given, so it can just be a commonly used regression model for multilevel data. Assuming $Y_{i j}$ follows a linear mixed-effects model, we have

Y_{i j} | T_{i j}, C_{i j} = k, X_{i j}, θ \sim N (τ_{k} T_{i j} + X_{i j}' β_{k} + Z_{i j}' u_{k j}, σ_{k}^{2}),

where $τ_{k}, β_{k}$ and $σ_{k}$ are the ATE, fixed effects, and standard deviation of individual-level error term of class k respectively, and

(u_{1 j}', u_{2 j}', \dots, u_{K j}')^{'} \sim N (0, Σ_{u})

are the random effects of cluster j for all K latent classes. Here, we again allow random effects from one latent class to be correlated with those from another within each cluster. The structure of this component is very similar to that of the selection model.

Mixing probability

This component models the probability that the individual belongs to each latent class given only its covariates. That is, the covariates of individuals from different latent classes may follow different distributions, so that the covariates alone give us information about latent class membership. Assuming $C_{i j}$ follows a multinomial logit model, we have

Pr (C_{i j} = k | X_{i j}, θ) = \frac{exp (X_{i j}' γ_{k})}{\sum_{ℓ = 1}^{K} exp (X_{i j}' γ_{ℓ})},

where $γ_{ℓ}$ is the coefficient vector of class $ℓ$ for $ℓ = 1, 2, \dots, K$ , and $γ_{1} = 0$ , so that the model is identified. This model simply tries to separate individuals from different latent classes with linear decision boundaries, so it works well for normally distributed covariates with the same covariance matrix across latent classes. Otherwise, nonlinear terms such as interactions help build better decision boundaries. Similar to the selection model and the outcome model, it is also possible to add hierarchical structure to (3) by considering random effects in decision boundaries:

Pr (C_{i j} = k | X_{i j}, θ) = \frac{exp (X_{i j}' γ_{k} + Z_{i j}' w_{k j})}{\sum_{ℓ = 1}^{K} exp (X_{i j}' γ_{ℓ} + Z_{i j}' w_{k ℓ})},

where

(w_{2 j}', w_{3 j}', \dots, w_{K j}')^{'} \sim N (0, Σ_{w}), w_{1 j} = 0,

are the random effects of cluster j for all K latent classes. Equation 4 allows different clusters to have different decision boundaries. A relevant practical issue about this extension is that unlike models (1) and (2), the outcome of (3) and (4) (i.e., $C_{i j}$ ) is unobserved, so the introduction of random effects will make the estimation much slower and more difficult. As will be described below, we can choose between (3) and (4) using information criteria. An example on how to compare these two will be shown in the case study section of this article.

Putting all the components above together, the log-likelihood function of parameter $θ = (τ, φ, β, σ, u, Σ_{u}, α, v, Σ_{v}, γ)$ is

\begin{array}{l} log L (θ) & = \sum_{j = 1}^{M} \sum_{i = 1}^{n_{j}} log p (Y_{i j}, T_{i j} | X_{i j}, θ) \\ = \sum_{j = 1}^{M} \sum_{i = 1}^{n_{j}} log [\sum_{k = 1}^{K} p (Y_{i j} | T_{i j}, C_{i j} = k, X_{i j}, θ) Pr (T_{i j} | C_{i j} = k, X_{i j}, θ) Pr (C_{i j} = k | X_{i j}, θ)] . \end{array}

With parameter estimates $\hat{θ}$ , the posterior probability of class membership can be computed through Bayes’s theorem:

Pr (C_{i j} = k | Y_{i j}, T_{i j}, X_{i j}, \hat{θ}) = \frac{p (Y_{i j} | T_{i j}, C_{i j} = k, X_{i j}, \hat{θ}) Pr (T_{i j} | C_{i j} = k, X_{i j}, \hat{θ}) Pr (C_{i j} = k | X_{i j}, \hat{θ})}{\sum_{ℓ = 1}^{K} p (Y_{i j} | T_{i j}, C_{i j} = ℓ, X_{i j}, \hat{θ}) Pr (T_{i j} | C_{i j} = ℓ, X_{i j}, \hat{θ}) Pr (C_{i j} = ℓ | X_{i j}, \hat{θ})} .

As in many model-based classification methods, our approach estimates which class each individual belongs to probabilistically, so we classify some individuals into a particular class more confidently than others. However, unlike most methods where the class membership is determined by one model (e.g., multinomial logit or cluster analysis), the proposed method fits the selection model, outcome model, and mixing probability jointly and simultaneously, and the three components collectively contribute to the estimation of class membership and possibly heterogeneous treatment effects.

Bayesian Estimation

Our goal is to estimate the latent class model with log-likelihood function (5). A Bayesian estimation procedure is applied here, given the complex structure of the model. First, we standardize all covariates in $X_{i j}$ and $Z_{i j}$ (except intercepts) and the outcome $Y_{i j}$ by using their z-scores before estimating $θ$ , so that the following estimation procedure becomes faster and convergence is achieved more quickly. Since it is a linear transformation, all coefficient estimates can be changed back to their original scales in the end. The following priors are used in this study:

\begin{array}{l} τ_{k}, φ_{k}, β_{k}, α_{k} \sim N (0, 1), & k = 1, 2, \dots, K, \\ γ_{1} = 0, γ_{k} \sim N (0, 1), & k = 2, 3, \dots, K, \\ (u_{1 j}', u_{2 j}', \dots, u_{K j}')^{'} \sim N (0, Σ_{u}), Σ_{u} = σ_{u}' Ω_{u} σ_{u}, & j = 1, 2, \dots, M, \\ (v_{1 j}', v_{2 j}', \dots, v_{K j}')^{'} \sim N (0, Σ_{v}), Σ_{v} = σ_{v}' Ω_{v} σ_{v}, & j = 1, 2, \dots, M, \\ σ_{k}, σ_{u}, σ_{v} \sim HalfCauchy (0, 1), & k = 1, 2, \dots, K, \\ Ω_{u}, Ω_{v} \sim LKJCorr (1) . \end{array}

These priors are mostly weakly informative because we hardly see very large standardized coefficients in practice. For specific research questions and data sets, we may want to use even weaker priors, although they may slow down the estimation procedure and make estimates less stable. We also recommend using more informative priors according to the specific research question. After the estimation procedure is done, the posterior mean of each parameter is used as the point estimate, the posterior standard deviation of each parameter measures the accuracy of the estimation, and specific quantiles of posterior samples can be used as bounds of credible intervals.

Akaike information criterion and Bayesian information criterion are widely used by frequentists for model comparison, but they are computed under the maximum likelihood framework (Vrieze, 2012). For Bayesian estimation, the leave-one-out cross-validation information criterion (LOOIC) can be used to select the best model (Vehtari et al., 2017). In this article, we fit the model with different Ks, compute their LOOIC, and test the differences between their expected predictive accuracy to decide the number of latent classes. More specifically, we consider two ways for selecting the number of latent classes K for the Bayesian method. The “minimum” way is to simply ignore the uncertainty in estimating LOOIC and choose the one with the smallest LOOIC. The other approach is a “sequential” process relying on testing the difference between expected log pointwise predictive densities (ELPDs) using $c = 2$ as the critical value. Suppose we have fit the data with $K = 1, 2, 3, 4$ . First, test whether $Δ ELPD$ is significantly different from zero for $K = 1$ and $K = 2$ . If $K = 1$ leads to smaller LOOIC or $Δ ELPD$ is not significant, then we choose $K = 1$ . Otherwise, we repeat this procedure for $K = 2$ and $K = 3$ , and then $K = 3$ and $K = 4$ . If $Δ ELPD$ is always significant, then we choose $K = 4$ .

There are some additional computational issues. First, in this study, we use Stan, which performs Hamiltonian Monte Carlo for Bayesian estimation (Gelman et al., 2015). We run four independent chains simultaneously. Second, the model is not identified because of label switching, that is, class labels can be arbitrarily permuted, and as a result, the four chains typically do not mix well. A common way to avoid this is to add a constraint to distinguish between different latent classes, such as $τ_{1} \leq τ_{2} \leq \dots \leq τ_{K}$ . This approach, however, does not work if, say, two latent classes have similar ATEs but different outcome models. One may also employ post hoc relabeling algorithms to the four chains, so that they have the same labels for each latent class (White & Murphy, 2014), but this approach is quite complex and requires additional computation. Our approach is to compare LOOIC for each chain and only keep samples from the chain with the lowest LOOIC. We lose information because other chains are discarded, but this approach is fast and easy to implement. Also, our simulation results show that this strategy leads to good performance.

In short, the model we propose uses three pieces of information to do classification: different distributions of covariates, different selection models, and different outcome models. The model here is very general because individuals from different latent classes are allowed to have not only different ATEs but also different coefficients, both fixed and random. Correspondingly, each cluster can have different random effects for different latent classes. Note that the formulation can be simplified or complicated depending on our assumptions. For example, we may assume that individuals from different latent classes actually follow the same outcome model except that they have different ATEs or, as described above, the mixing probability model also has a hierarchical structure. The corresponding estimation procedures are similar.

Doubly Robust Estimation for Model Checking

For each latent class k, ${\hat{τ}}_{k}$ , the posterior mean of $τ_{k}$ , is the point estimator of $τ_{k}$ based on the outcome model only. It is possible to incorporate the propensity score $e_{k i j} = Pr (T_{i j} = 1 | C_{i j} = k, X_{i j}, θ)$ from the selection model into the outcome model to obtain more robustness for estimating the ATE. Let

Y_{k i j} (0) = X_{i j}' β_{k} + Z_{i j}' u_{k j} and Y_{k i j} (1) = τ_{k} + X_{i j}' β_{k} + Z_{i j}' u_{k j}

be the two potential outcomes for class k, and let

Δ_{k i j} = [\frac{T_{i j}}{e_{k i j}} Y_{i j} + \frac{e_{k i j} - T_{i j}}{e_{k i j}} Y_{k i j} (1) - \frac{1 - T_{i j}}{1 - e_{k i j}} Y_{i j} - \frac{T_{i j} - e_{k i j}}{1 - e_{k i j}} Y_{k i j} (0)] .

Then, the weighted average

{\tilde{τ}}_{k} = \frac{\sum_{{i, j : 0.1 < e_{k i j} < 0.9}} Δ_{k i j} \cdot p_{k i j}}{\sum_{{i, j : 0.1 < e_{k i j} < 0.9}} p_{k i j}}

is the doubly robust estimator of the ATE of latent class k, where $p_{k i j} = Pr (C_{i j} = k | Y_{i j}, T_{i j}, X_{i j}, θ)$ is the posterior probability that individual i from cluster j belongs to class k. The double robustness property is desirable in practice because it is consistent when either the selection model or the outcome model is correctly specified (Scharfstein et al., 1999). Moreover, for numerical stability, we only keep observations whose propensity scores are between $0.1$ and $0.9$ for computing the doubly robust estimator (Crump et al., 2009). As postestimation quantities, one can easily change these cutoff values and check whether ${\tilde{τ}}_{k}$ is stable. We compute ${\tilde{τ}}_{k}$ for each posterior sample during the Bayesian estimation process, but these samples cannot be regarded as draws from the posterior distribution of $τ_{k}$ because ${\tilde{τ}}_{k}$ is an estimator of $τ_{k}$ rather than the parameter $τ_{k}$ itself. As a result, the probability that the true $τ_{k}$ lies in the $100 (1 - α) %$ credible interval of ${\tilde{τ}}_{k}$ is smaller than $100 (1 - α) %$ . This suggests that ${\tilde{τ}}_{k}$ should not be used for interpretation. Only posterior samples of $τ_{k}$ from model (2) give valid information about the posterior distribution of $τ_{k}$ .

Still, we propose using the posterior distribution of ${\tilde{τ}}_{k}$ for model checking on possibly underspecified K and misspecification of the selection or the outcome model. For each ${\tilde{τ}}_{k}$ , we check whether the outcome-only point estimate ${\hat{τ}}_{k}$ is within the $50 %$ credible interval of ${\tilde{τ}}_{k}$ . If it is not true, then we suspect misfit. The idea behind this approach is that if there is no model misspecification, then since ${\hat{τ}}_{k}$ and ${\tilde{τ}}_{k}$ are computed from the same data, they should both provide valid estimates of $τ_{k}$ and hence be extremely close. Then, we expect that ${\hat{τ}}_{k}$ is within the credible interval of ${\tilde{τ}}_{k}$ at a low level such as $50 %$ (Gelman et al., 2020; Park et al., 2006).

Other Estimation Approaches

We have already shown how to estimate the model (5) using Bayesian methods. Another popular way for estimating mixture models is the expectation–maximization (EM) algorithm. It is an iterative procedure, where during each iteration, we estimate the models based on the current class membership and then reassign individuals to new classes with the highest likelihood. More specifically, this procedure is as follows:

Randomly assign each individual to one of the K latent classes.

Estimate the mixing probability model based on the current class membership. Within each latent class, estimate its selection model and outcome model.

For each individual, compute their log-likelihood function for the K latent classes respectively. This log-likelihood function is the sum of the log-likelihood from the mixing probability model, the selection model, and the outcome model.

Assign each individual to the latent class with the highest log-likelihood.

If the total number of iterations exceeds our preset limit or the class membership stays unchanged, then stop and choose the iteration with the highest total log-likelihood as the output. Otherwise, go back to step 2.

The EM algorithm essentially estimates the same models as the Bayesian procedure except that we do not explicitly let random effects of different latent classes be correlated because models of each latent class are estimated separately. As will be shown in the simulation later, the EM algorithm is much less stable than the Bayesian method when the sample size is small due to the presence of singularities under the maximum likelihood framework (Bishop, 2006, pp. 433–434), while priors of the Bayesian method help restrict parameters from being too off. Another disadvantage of the EM algorithm is that obtaining standard errors and confidence intervals for parameters is not straightforward. Bootstrap is applicable, but it can be very computationally intensive.

For comparison, we also consider a sequential estimation procedure where we simply apply the k-means algorithm to covariates for classifying individuals and then estimate parameters based on the class membership. The procedure is similar to the EM algorithm described above except that in step 3, the log-likelihood function is replaced by the negative of the Euclidean distance between the individual and the mean of each latent class. This approach is fast and easy, but its classification only takes mixing probability into account while all information within the outcome model and the selection model is discarded. Therefore, we expect that it has good performance only when the distributions of covariates of different latent classes are so different that k-means alone separates latent classes well. Moreover, this sequential approach does not provide standard errors or confidence intervals for parameters directly, so bootstrap is needed.

Simulation Study

There are three important questions regarding our proposed method: How well does this method classify individuals into their latent classes? How accurate are the ATE estimates? Can this approach find the correct number of latent classes? These questions will be answered through simulation studies. We consider three different scenarios here, where for each case, there are N individuals from M clusters, and random effects consist of random intercepts only.

Data Generating Process and Estimation

The data generating process for the simulation study is as follows. First, each individual is assigned to one of the three classes according to $Pr (C_{n} = k) = p_{k}$ for $k = 1, 2, 3$ and generates Level 1 covariates by

x_{n}^{(1)} | C_{n} = k \sim N (μ_{k}, [\begin{matrix} 3 & 2 \\ 2 & 3 \end{matrix}]),

where

μ_{1} = [\begin{matrix} - 1 \\ - 2 \end{matrix}], μ_{2} = [\begin{matrix} 0 \\ 0 \end{matrix}], μ_{3} = [\begin{matrix} 2 \\ 1 \end{matrix}] .

Second, generate Level 2 covariates by

x_{j}^{(2)} \sim N ([\begin{matrix} 0 \\ 0 \end{matrix}], [\begin{matrix} 4 & - 1 \\ - 1 & 5 \end{matrix}]),

and generate random effects by

\begin{array}{l} [\begin{matrix} u_{1 j} \\ u_{2 j} \\ u_{3 j} \end{matrix}] \sim N ([\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}], [\begin{matrix} 3 \\ 4 \\ 2 \end{matrix}] [\begin{matrix} 1 & 0.6 & 0.6 \\ 0.6 & 1 & 0.6 \\ 0.6 & 0.6 & 1 \end{matrix}] [\begin{matrix} 3 \\ 4 \\ 2 \end{matrix}]), \\ [\begin{matrix} v_{1 j} \\ v_{2 j} \\ v_{3 j} \end{matrix}] \sim N ([\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}], [\begin{matrix} 4 \\ 3 \\ 5 \end{matrix}] [\begin{matrix} 1 & 0.4 & 0.4 \\ 0.4 & 1 & 0.4 \\ 0.4 & 0.4 & 1 \end{matrix}] [\begin{matrix} 4 \\ 3 \\ 5 \end{matrix}]) . \end{array}

Third, assign individuals to clusters by the multinomial logit model

Pr (J (n) = j) = \frac{exp (0.04 x_{n 1}^{(1)} x_{j 1}^{(2)} + 0.02 x_{n 1}^{(1)} x_{j 2}^{(2)} + 0.03 x_{n 2}^{(1)} x_{j 1}^{(2)} - 0.02 x_{n 2}^{(1)} x_{j 2}^{(2)})}{\sum_{ℓ = 1}^{M} exp (0.04 x_{n 1}^{(1)} x_{ℓ 1}^{(2)} + 0.02 x_{n 1}^{(1)} x_{ℓ 2}^{(2)} + 0.03 x_{n 2}^{(1)} x_{ℓ 1}^{(2)} - 0.02 x_{n 2}^{(1)} x_{ℓ 2}^{(2)})},

and then, the covariate for individual n who is the ith individual from cluster $j = J (n)$ is

X_{i j} = [\begin{matrix} 1 & x_{n}^{(1)}' & x_{j}^{(2)}' \end{matrix}]' .

Since the probability that each individual belongs to each cluster depends on the interaction between individual-level and cluster-level covariates, correlation is created between these two kinds of variables. Fourth, generate treatment indicator $T_{i j}$ by

Pr (T_{i j} = 1) = {\begin{matrix} {logit}^{- 1} (v_{1 j} + [\begin{matrix} - 0.5 & 0.9 & - 0.7 & 0.4 & 0.2 \end{matrix}] X_{i j}), & C_{i j} = 1, \\ {logit}^{- 1} (v_{2 j} + [\begin{matrix} - 0.1 & - 0.6 & - 0.2 & 0.1 & - 0.2 \end{matrix}] X_{i j}), & C_{i j} = 2, \\ {logit}^{- 1} (v_{3 j} + [\begin{matrix} 0.6 & - 0.8 & 0.7 & - 0.5 & 0.1 \end{matrix}] X_{i j}), & C_{i j} = 3, \end{matrix}

where $C_{i j}$ is the class membership of the ith individual from cluster j. Finally, generate outcome $Y_{i j}$ by

\begin{array}{l} Y_{i j} & \sim N (μ_{i j}, σ_{C_{i j}}^{2}), \\ μ_{i j} & = {\begin{matrix} u_{1 j} + τ_{1} T_{i j} + [\begin{matrix} 5 & 7 & - 5 & 1 & - 2 \end{matrix}] X_{i j}, & C_{i j} = 1, \\ u_{2 j} + τ_{2} T_{i j} + [\begin{matrix} 8 & - 5 & 0 & - 1 & 1 \end{matrix}] X_{i j}, & C_{i j} = 2, \\ u_{3 j} + τ_{3} T_{i j} + [\begin{matrix} - 1 & 2 & 3 & 1 & - 1 \end{matrix}] X_{i j}, & C_{i j} = 3, \end{matrix} \end{array}

where $τ = [\begin{matrix} τ_{1} & τ_{2} & τ_{3} \end{matrix}]' = [\begin{matrix} - 3 & 0 & 5 \end{matrix}]'$ are the true treatment effects and $σ = [\begin{matrix} σ_{1} & σ_{2} & σ_{3} \end{matrix}]' = [\begin{matrix} 4 & 3 & 5 \end{matrix}]'$ are the standard deviations of Level 1 error terms for the three latent classes.

We include all the five covariates (intercept, x ₁, x ₂, z ₁, and z ₂) as fixed effects in all the three models—mixing probability, selection model, and outcome model—and random intercepts are included in the latter two models. Since the covariance matrices of x for the three latent classes are the same, the corresponding best decision boundaries are indeed linear curves that do not depend on the cluster, so random effects are not needed in the mixing probability model. Therefore, all the three models are correctly specified.

Three different numbers of latent classes are considered: When $K = 1$ , we let $p_{1} = p_{2} = 0$ and $p_{3} = 1$ ; when $K = 2$ , we let $p_{1} = 0.4$ , $p_{2} = 0$ , and $p_{3} = 0.6$ ; when $K = 3$ , we let $p_{1} = 0.2$ , $p_{2} = 0.3$ , and $p_{3} = 0.5$ . For each scenario, $200$ replications of simulated data sets are generated. Table 1 provides a summary of cluster sizes across all the replications. The overall average cluster sizes are fixed to be $25$ . Our data generating process generates clusters of varying while reasonable sizes: More than half of the clusters have $20$ to $30$ individuals, and even the smallest cluster still has six individuals; each cluster has at least one individual from each latent class. Moreover, the selection models lead to good overlap of true propensity scores: Across all the replications, more than $50 %$ of individual propensity scores are between $0.2$ and $0.8$ , and more than $99 %$ are between $0.01$ and $0.99$ .

Table 1.

Summary of Cluster Sizes of Simulation Study

K	N	M	Class	Minimum	2.5%	25%	50%	Mean	75%	97.5%	Maximum
1	250	10	3	8	14	20	25	25.00	29	40	55
			All	8	14	20	25	25.00	29	40	55
	500	20	3	7	13	20	24	25.00	29	41	60
			All	7	13	20	24	25.00	29	41	60
	1,000	40	3	6	13	20	24	25.00	29	41	66
			All	6	13	20	24	25.00	29	41	66
2	250	10	1	1	4	8	10	10.01	12	18	24
			3	3	7	12	15	14.99	18	25	36
			All	11	16	21	25	25.00	28	35	45
	500	20	1	1	4	8	10	10.02	12	18	31
			3	2	7	11	15	14.98	18	26	40
			All	9	15	21	25	25.00	28	36	50
	1,000	40	1	1	4	7	10	9.97	12	18	28
			3	2	7	12	15	15.03	18	26	40
			All	9	15	21	25	25.00	28	36	52
3	250	10	1	1	1	3	5	5.01	6	10	14
			2	1	3	6	7	7.49	9	13	20
			3	3	6	10	12	12.56	15	22	30
			All	9	16	21	25	25.00	28	35	44
	500	20	1	1	1	3	5	4.99	6	10	17
			2	1	3	6	7	7.54	9	13	18
			3	1	5	9	12	12.52	15	22	36
			All	6	15	21	25	25.00	28	36	50
	1,000	40	1	1	1	3	5	5.07	6	10	19
			2	1	3	6	7	7.53	9	13	19
			3	1	5	9	12	12.45	15	22	34
			All	6	15	21	25	25.00	28	36	46

Note. K = number of latent classes; N = number of Level 1 units; M = number of Level 2 units.

Each data set is estimated using all the three methods: the Bayesian method, the EM algorithm, and the sequential procedure based on k-means. For the Bayesian estimation, we run $2, 000$ iterations for each of the four chains and use the mean value over the last $1, 000$ iterations from the chain with the lowest LOOIC as our point estimate. Our criterion for convergence is that maximum split $\hat{R}$ over all parameters for the chain chosen does not exceed $1.1$ (Vehtari et al., 2021). For each replication, the Bayesian estimation procedure is rerun repeatedly until convergence. For the EM algorithm, we run at most $1, 000$ iterations for $100$ random initial starting points and then pick the iteration with the highest log-likelihood as our point estimate. For k-means, we run at most $1, 000$ iterations for $100$ random initial starting points and then pick the iteration with the smallest sum of squared distances from each individual to its center. After estimation, we reorder the latent classes, so that ${\hat{τ}}_{k}$ , the point estimates of ATE, are increasing. The Stan code used in this study for the Bayesian estimation procedure is in Supplementary Material.

Results

Table 2 shows the summary statistics (means and standard deviations) of point estimates ${\hat{τ}}_{k}$ across replications and the coverage of the $95 %$ credible intervals. More specifically, for each scenario and each estimation method, we have $200$ point estimates of ATEs from $200$ replications, and Table 2 shows the mean and the standard deviation of these $200$ point estimates. In addition, since the Bayesian approach provides a natural way to find the credible intervals of each parameter directly from posterior samples, we report the mean coverage of $95 %$ credible intervals across all replications. There is no straightforward way to find the confidence intervals of parameters for EM and k-means except bootstrap, which is computationally intensive. Therefore, we do not compute the coverage of confidence intervals for EM and k-means here. It is an advantage of the Bayesian approach: We have more information about parameter estimates.

Table 2.

Summary of Average Treatment Effect Point Estimates ( ${\hat{τ}}_{k}$ ) and Coverage of Bayesian Credible Intervals Across Replications of Simulation Study

K	N	M	Parameter	True	Bayesian			EM		k-Means
K	N	M	Parameter	True	Mean	SD	Coverage	Mean	SD	Mean	SD
1	250	10	τ₁	5	4.98	0.73	.96	5.02	0.73	5.02	0.73
1	500	20	τ₁	5	4.93	0.54	.94	4.94	0.54	4.94	0.54
	1,000	40	τ₁	5	5.03	0.40	.94	5.04	0.40	5.04	0.40
2	250	10	τ₁	−3	−2.69	1.22	.95	−2.92	1.41	2.29	2.01
			τ₂	5	5.05	1.20	.94	5.12	1.29	6.20	1.46
	500	20	τ₁	−3	−2.86	0.88	.92	−2.93	1.01	2.31	1.41
			τ₂	5	5.02	0.82	.96	4.99	0.84	6.31	0.94
	1,000	40	τ₁	−3	−2.99	0.57	.94	−3.13	0.69	2.36	0.94
			τ₂	5	5.06	0.60	.94	4.97	0.64	6.26	0.70
3	250	10	τ₁	−3	−2.43	1.77	.94	−3.24	14.54	2.15	2.89
			τ₂	0	0.39	1.22	.98	2.08	3.14	5.18	1.75
			τ₃	5	5.10	1.29	.96	8.34	8.48	7.59	2.03
	500	20	τ₁	−3	−2.65	1.32	.94	−2.74	1.88	3.03	2.00
			τ₂	0	0.05	0.76	.98	0.24	1.20	5.30	1.16
			τ₃	5	5.08	0.93	.94	5.65	4.96	6.78	1.43
	1,000	40	τ₁	−3	−2.88	0.99	.95	−2.92	1.40	3.34	1.45
			τ₂	0	0.08	0.56	.96	0.20	0.72	5.23	0.92
			τ₃	5	5.06	0.64	.97	5.03	0.69	6.57	0.98

Note. K = number of latent classes; N = number of Level 1 units; M = number of Level 2 units; $τ_{k}$ = ATE of latent class k; EM = expectation–maximization.

The Bayesian method gives good results according to Table 2: The mean of ${\hat{τ}}_{k}$ is close to true $τ_{k}$ , and as the sample size N increases, the bias generally gets smaller; the coverage of the $95 %$ credible interval is close to $0.95$ . The EM algorithm gives good estimates for $K = 1$ and $K = 2$ as well, although the standard deviation is greater than that of the Bayesian method. When $K = 3$ , however, the EM algorithm has poor performance by showing large bias and huge variation unless the sample size is large enough. We find that EM frequently groups many outliers, such as treated cases with large outcomes and untreated cases with small outcomes, into a single latent class with extreme parameter estimates. The problem is that the EM algorithm tends to care too much about increasing the log-likelihood of these outliers, and consequently, the overall model estimation is “improved” by grouping them into an extreme latent class. These singularities make the maximization of log-likelihood an ill-posed problem (Bishop, 2006, pp. 433–434). In contrast, the Bayesian approach seldom faces this issue and is more stable even with small sample sizes because it relies on weakly informative priors that help avoid such extreme results. The estimates of the k-means approach are not reasonable, suggesting that the mixing probability model itself does not provide enough information for classification, at least for the simulated data.

Table 3 shows the mean correct classification rates across replications. Again, the Bayesian method shows the best performance by giving the highest mean correct classification rates. The EM algorithm is almost as good as the Bayesian method when $K = 2$ , but is much worse when $K = 3$ . The results of k-means do not look too bad when $K = 2$ , which is a bit surprising, given that the corresponding ATE estimates in Table 2 are very inaccurate. Actually, even if we take the much better classification results from the Bayesian procedure and then estimate the outcome model for each estimated latent class separately, we end up with much worse ATE estimates than those from the Bayesian joint estimation procedure shown in Table 2. Although only fewer than 10% individuals are classified incorrectly on average, the classification errors hurt ATE estimation in the next step. This supports our claim that joint estimation leads to better results, because unlike sequential estimation, all components of model (5) work together rather than transmit errors to each other. K-means looks like random guessing when $K = 3$ , again suggesting that solely relying on different distributions of covariates for classification is far from sufficient for the data generating process. Another finding is that it seems easier to classify the untreated than the treated, although this may be due to the specific data generating process. This phenomenon may be related to the uncertainty in estimating the ATE, which transfers to the difficulty in classifying the treated. In short, the Bayesian method is the best in the sense that its estimates are most accurate and stable, and it correctly classifies the greatest percentage of observations in most cases.

Table 3.

Mean Correct Classification Rates Across Replications of Simulation Study

K	N	M	Bayesian			Expectation–Maximization			K-Means
K	N	M	All	Untreated	Treated	All	Untreated	Treated	All	Untreated	Treated
2	250	10	.92	.95	.89	.91	.94	.87	.74	.80	.67
2	500	20	.93	.96	.90	.92	.95	.88	.80	.85	.73
	1,000	40	.93	.96	.90	.93	.96	.89	.80	.86	.74
3	250	10	.79	.81	.77	.53	.55	.49	.35	.36	.35
	500	20	.84	.85	.81	.79	.81	.76	.34	.35	.33
	1,000	40	.86	.88	.84	.83	.85	.81	.35	.37	.34

Note. K = number of latent classes; N = number of Level 1 units; M = number of Level 2 units.

Table 4 shows the results of the two ways for selecting the number of latent classes K for the Bayesian method. K is never underestimated: $k < K$ is never picked by either approach, which is desirable. Also, both approaches pick the correct K in most cases. The minimum approach sometimes overestimates K and can choose four latent classes even when there is only one. The sequential approach at most overestimates K by one, and even this kind of overestimation is very rare. Therefore, we recommend the sequential approach, which has satisfying performance.

Table 4.

Frequency Table of Estimated Number of Latent Classes Across Replications of Simulation Study

K	N	M	Approach	k = 1	k = 2	k = 3	k = 4
1	250	10	Minimum	178	15	5	2
			Sequential	200	0	0	0
	500	20	Minimum	181	14	4	1
			Sequential	199	1	0	0
	1,000	40	Minimum	162	27	10	1
			Sequential	196	4	0	0
2	250	10	Minimum	0	177	17	6
			Sequential	0	200	0	0
	500	20	Minimum	0	173	22	5
			Sequential	0	197	3	0
	1,000	40	Minimum	0	183	14	3
			Sequential	0	199	1	0
3	250	10	Minimum	0	0	167	33
			Sequential	0	0	195	5
	500	20	Minimum	0	0	158	42
			Sequential	0	0	195	5
	1,000	40	Minimum	0	0	166	34
			Sequential	0	0	196	4

Note. N = number of Level 1 units; M = number of Level 2 units; K = true number of latent classes; k = estimated number of latent classes.

As additional model checking, Table 5 shows whether the point estimate ${\hat{τ}}_{k}$ lies outside of the $50 %$ credible intervals of ${\tilde{τ}}_{k}$ . When $K = 2$ or $K = 3$ , there is very little misfit detected when K is at least the true number of classes, while more misfit is found when K is smaller than the true value. The power of detecting underspecification of K increases as the sample size increases. The results do not look very good when there is only one latent class. We emphasize that ${\tilde{τ}}_{k}$ is an estimator rather than the parameter $τ_{k}$ itself, so its posterior samples cannot be regarded as samples from the posterior distribution of parameter $τ_{k}$ : ${\tilde{τ}}_{k}$ has too small a posterior standard deviation and too narrow a credible interval for $τ_{k}$ , while ${\hat{τ}}_{k}$ provides correct information. In short, although ${\tilde{τ}}_{k}$ itself is not used as parameter estimate to report, we think ${\tilde{τ}}_{k}$ provides a promising way to check model fit, especially for large sample sizes, although future study is needed for finding the best way to take advantage of this information.

Table 5.

Percentages That Misfit Is Detected Across Replications of Simulation Study

K	N	M	k = 1	k = 2	k = 3	k = 4
1	250	10	.25	.15	.13	.17
1	500	20	.38	.11	.23	.28
	1,000	40	.29	.06	.18	.24
2	250	10	.42	.04	.03	.01
2	500	20	.61	.07	.05	.03
	1,000	40	.78	.09	.03	.03
3	250	10	.46	.10	.01	.05
	500	20	.48	.25	.02	.04
	1,000	40	.58	.60	.01	.04

Note. N = number of Level 1 units; M = number of Level 2 units; K = true number of latent classes; k = estimated number of latent classes.

The simulation study has shown that our proposed method is a good way for classifying observations and estimating ATEs: The correct classification rate is high, the ATE estimates are accurate and have the correct coverage, and in almost all cases, the correct number of latent classes is found. The performance of this method generally improves when the sample size increases.

Application With Trends in International Mathematics and Science Study (TIMSS) Data

Data and Variables

We use the 2015 TIMSS data of eighth graders in Korea to demonstrate our proposed approach. Briefly speaking, TIMSS is a large-scale educational assessment that investigates the progress of students’ performance in mathematics and science. Since 1995, TIMSS has been collected by the International Association for the Evaluation of Educational Achievement, and it has been conducted for Grades 4 and 8 every 4 years in more than 40 countries. The data are collected by a two-stage stratified cluster sampling, where schools are first selected based on school demographic variables and then at least one classroom is randomly selected within each school (Martin et al., 2016).

Our treatment variable is whether a student received a private science lesson or not. The outcome of interest is the first plausible value of the science achievement. For simplicity, we do not incorporate multiple plausible values of student science achievements or sampling weights in our data analysis. A set of covariates include six student-level covariates (student’s gender, father’s highest education level, the number of books at home, student's home study supports, student’s confidence in science, and student’s perceived value of science) and six school-level covariates (school’s gender type, the percentage of economically disadvantaged students, school location, science instruction affected by resource shortage, school’s emphasis on academic success, and school discipline problems). Table 6 summarizes all the covariates used in the analysis. Our analysis suggests nonlinearity of the variable on student’s confidence in science, so its quadratic term is included in the models as well. Note that some of the covariates, such as student’s confidence in science and perceived value of science, may be post-treatment confounders. If it is the case, then both controlling them and omitting them will lead to bias in the outcome regression model, so more complex treatment is needed (Moerkerke et al., 2015). We do not dive deeper into the rationality of variable selection shown in Table 6. In practice, researchers should rely on subject matter theory to determine which covariates to use and the time that each variable is measured is useful for judging whether each covariate is a pretreatment or a post-treatment confounder. After removing observations with missing values, our final sample consists of $4, 874$ students from $149$ schools, among which $seven$ have fewer than $20$ students and $10$ have more than $60$ students, and more than half of the schools have around $30$ students.

Table 6.

All Covariates Used in the TIMSS Analysis

Variable	Level	Description
male	Student	Student is male
dad_cll	Student	Student’s father is a college graduate
dad_q	Student	Student’s father’s education level is unknown
books25	Student	Student has more than 25 books at home
hspprt	Student	Student has both their own room and internet connection
conf	Student	Student’s confidence in science
value	Student	Student’s perceived value of science
girlsch	School	School has girls only
coedu	School	School has both boys and girls
disad_11	School	Percentage of economically disadvantaged students is between 11% and 25%
disad_26	School	Percentage of economically disadvantaged students is between 26% and 50%
disad_m50	School	Percentage of economically disadvantaged students is more than 50%
city_u	School	School is in urban area
city_sub	School	School is in suburban area
city_m	School	School is in a medium size city
resshort	School	School’s science instruction affected by resource shortage
academph	School	School’s emphasis on academic success
dscpn	School	School discipline problems

Note. $N = 4, 874$ students; $M = 149$ schools. TIMSS = Trends in International Mathematics and Science Study.

Procedure

As a demonstration of how to use our proposed method in practice, our data analysis procedure is as follows:

For each different postulated number of latent classes, fit the proposed model using the Bayesian method described above. We suggest running several independent chains, discarding the first part of each chain as warmup iterations, and keeping only the chain with the smallest LOOIC. Then, check whether the chain chosen converges. One popular criterion is to make sure that the maximum $\hat{R}$ across all the parameters is smaller than $1.05$ or $1.1$ . If convergence is not achieved, then rerun the estimation procedure with possibly more iterations.

Use the sequential approach to decide the number of latent classes: if $K = 1$ leads to smaller LOOIC than $K = 2$ does, or $Δ E L P D$ between these two is not significant, then choose $K = 1$ ; otherwise, we repeat the same procedure to compare $K = 2$ and $K = 3$ , then $K = 3$ and $K = 4$ , and so on, until the best number of latent classes is found.

For the best K, check whether all posterior means ${\hat{τ}}_{k}$ are within their corresponding $50 %$ credible intervals of ${\tilde{τ}}_{k}$ . If not, then there may be model misspecification.

Summarize the results: For each parameter, the posterior mean is the point estimate, the $95 %$ credible interval is the interval estimate, and the posterior standard deviation measures the uncertainty of the parameter. For each individual, its estimated class membership is the latent class with the highest posterior probability.

It is worth noting that we report ${\hat{τ}}_{k}$ , along with its corresponding posterior standard deviation and credible intervals, as our final estimates. The doubly robust estimate ${\tilde{τ}}_{k}$ is only for model checking because its posterior samples do not reflect the posterior distribution of $τ_{k}$ .

In this TIMSS analysis, both selection and outcome models include random intercepts just like the simulation study, and mixing probability models both with and without random intercepts are considered. We fit the model using the Bayesian procedure with $3, 000$ iterations for each of the four chains, where the last $1, 500$ iterations from the chain with the lowest LOOIC are used for computing posterior distributions of parameters. Our criterion for convergence is that maximum split $\hat{R}$ over all parameters for the chain chosen is smaller than $1.05$ . The codes for this TIMSS analysis are shown in Supplementary Material.

Results

Table 7 shows the model comparison results that help decide the number of latent classes and whether random intercepts are needed in the mixing probability model. Consider random effects in the mixing probability model first. The mixing probability model vanishes when $K = 1$ . When $K = 2$ or 3, random effects actually lead to greater LOOICs and hence worse model fit. Looking at the 95% credible intervals of the standard deviations of these random intercepts, we find that their lower bounds are very small and close to zero, which again suggests that these random effects are not necessary. Moreover, parameter estimates only change a little after random effects are added. Therefore, we exclude random effects from mixing probability. Then, consider the number of latent classes. Using critical value $c = 2$ , Table 7 suggests that $K = 2$ is significantly better than $K = 1$ , and the difference between $K = 3$ and $K = 2$ is not significant although $K = 3$ leads to smaller LOOIC, so the sequential process chooses $K = 2$ . In brief, we choose $K = 2$ latent classes and do not include random intercepts in the mixing probability model.

Table 7.

Model Comparison (Number of Latent Classes and Whether Mixing Probability Includes Random Intercepts [RIs]) for TIMSS Analysis

RI in Mixing	K	LOOIC	ELPD	ΔELPD	SE(ΔELPD)
No	1	17,868.9	−8,934.5
No	2	17,579.8	−8,789.9	162.5	21.2
	3	17,544.0	−8,772.0	17.9	11.1
Yes	1	17,868.9	−8,934.5
	2	17,581.5	−8,790.7	160.3	21.0
	3	17,548.3	−8,774.2	16.6	10.6

Note. K = number of latent classes; LOOIC = leave-one-out cross-validation information criterion; ELPD = expected log pointwise predictive density; TIMSS = Trends in International Mathematics and Science Study.

Tables 8 and 9 show the summary of ATE estimates and class memberships, respectively. Although we have already decided that $K = 2$ , results for $K = 1$ and $K = 3$ are also shown for comparison. Only when $K = 1$ does ${\hat{τ}}_{k}$ lie outside of the $50 %$ credible interval of ${\tilde{τ}}_{k}$ , suggesting that there may be model misspecification for the one class solution, but no model misspecification is detected for two or three classes, so we only consider $K = 2$ or $K = 3$ . When $K = 3$ , the first latent class has very unstable ATE estimates and very few students, suggesting either overfitting or that this class is too small to be reliably estimated. Therefore, two latent classes seem sufficient for the data. Among all the $149$ schools in the sample, $86$ have at least $40 %$ students in both latent classes, while there are only five schools with more than $80 %$ students in the same class. Since most schools have large proportions of students in both latent classes, it is necessary to do the classification at the student level rather than at the school level.

Table 8.

Summary of Posterior Distributions of Average Treatment Effects and Doubly Robust Estimates of TIMSS Analysis

K	Parameter	Mean	SD	2.5%	25%	50%	75%	97.5%
1	$τ_{1}$	4.34	1.96	0.66	2.97	4.31	5.66	8.21
	${\tilde{τ}}_{1}$	4.75	0.49	3.76	4.41	4.77	5.08	5.72
2	$τ_{1}$	−0.64	2.93	−6.14	−2.59	−0.70	1.41	5.14
	${\tilde{τ}}_{1}$	−0.81	1.38	−3.43	−1.74	−0.83	0.12	2.06
	$τ_{2}$	12.16	3.62	5.19	9.88	12.16	14.52	19.39
	${\tilde{τ}}_{2}$	11.97	2.00	8.15	10.66	11.95	13.25	15.94
3	$τ_{1}$	−16.91	54.86	−127.25	−55.33	−14.73	21.84	86.39
	${\tilde{τ}}_{1}$	−18.07	55.51	−126.04	−57.43	−16.91	20.31	88.45
	$τ_{2}$	0.33	3.03	−5.84	−1.73	0.39	2.42	6.34
	${\tilde{τ}}_{2}$	0.16	1.67	−3.10	−0.98	0.14	1.30	3.43
	$τ_{3}$	13.10	3.72	5.81	10.65	13.06	15.55	20.47
	${\tilde{τ}}_{3}$	12.85	2.11	8.66	11.43	12.83	14.23	17.22

Note. Note that ${\tilde{τ}}_{k}$ is only used for model checking. Only rows corresponding to $τ_{k}$ provide valid information about the posterior distribution of average treatment effect. $τ_{k}$ = average treatment effect of latent class k; ${\tilde{τ}}_{k}$ = doubly robust estimate of $τ_{k}$ ; K = number of latent classes; TIMSS = Trends in International Mathematics and Science Study.

Table 9.

Number of Students in Each Latent Class of TIMSS Analysis

K	Class 1	Class 2	Class 3
1	4,874
2	2,444	2,430
3	17	2,528	2,329

Note. K = number of latent classes; TIMSS = Trends in International Mathematics and Science Study.

As discussed in the simulation study, doubly robust estimates ${\tilde{τ}}_{k}$ are solely used for checking model misspecification. Only rows corresponding to $τ_{k}$ in Table 8 provide valid information about the posterior distribution of ATE. It is essential that only these rows are used when interpreting ATE estimates.

Table 10 shows some descriptive statistics of the two latent classes. Both latent classes have approximately half of the students. Class 1 has almost zero ATE, suggesting that on average, the treatment has little effect on these students’ science scores; the ATE estimate of Class 2 is around 12 and is significant, although this effect is still not large practically, given that the standard deviations of the science score in the whole sample and the subsample of Latent Class 2 are 76.15 and 64.58, respectively. Class 1 has a higher percentage of treated students than Class 2, which looks largely ineffective. It would be better if more students in Class 2 were treated while fewer students in Class 1 were treated. In addition, Class 1 has much higher mean science scores than Class 2, and the private science lesson is far from able to fill this gap. This might be related to the ceiling effect: It is more difficult for students already good at science to further increase their scores, so they are more likely to be in Class 1, which has almost zero ATE.

Table 10.

Descriptive Statistics of the Two Latent Classes of TIMSS Analysis

k	n_k	${\hat{τ}}_{k}$	${\bar{t}}_{k}$ (%)	${\bar{y}}_{k} (0)$	${\bar{y}}_{k} (1)$	${\bar{y}}_{k}$	$sd (y_{k})$
1	2,444	−0.64	36.58	597.91	602.91	599.74	60.66
2	2,430	12.16	27.94	508.44	525.30	513.15	64.58

Note. TIMSS = Trends in International Mathematics and Science Study; n_k = number of individuals in latent class k; ${\hat{τ}}_{k}$ = average treatment effect point estimate of latent class k; ${\bar{t}}_{k}$ = percentage of treated of latent class k; ${\bar{y}}_{k} (0)$ = mean outcome on the untreated of latent class k; ${\bar{y}}_{k} (1)$ = mean outcome on the treated of latent class k; ${\bar{y}}_{k}$ = mean outcome of latent class k; $s d (y_{k})$ = standard deviation of outcome of latent class k.

Table 11 shows the posterior means and standard deviations of fixed effects parameters for the two latent class solution. The two classes have very different selection models and outcome models, so it is reasonable to divide students into different latent classes with different model parameters. Generally, student-level variables are important in all the models: These variables are related to class membership and affect both the probability of taking private science lessons and the science score. School-level variables affect students’ science scores, but there is little evidence that they are also related to the class membership or whether students take private science lessons.

Table 11.

Posterior Means and Standard Deviations of Model Parameters for $K = 2$ Latent Classes of TIMSS Analysis

Variable	Mixing	Selection		Outcome
Variable	$γ_{1}$	$α_{1}$	$α_{2}$	$β_{1}$	$β_{2}$
Intercept	18.774	−2.455	−2.547	134.455	381.195
	(3.674)	(1.502)	(1.008)	(44.758)	(25.362)
male	0.789	0.209	−0.175	3.948	−5.744
	(0.404)	(0.126)	(0.156)	(3.383)	(4.260)
dad_cll	0.377	0.114	0.062	13.138	−0.970
	(0.422)	(0.124)	(0.163)	(3.349)	(4.484)
dad_q	1.084	−0.015	−0.099	1.177	−16.721
	(0.428)	(0.156)	(0.150)	(4.237)	(4.097)
books25	−0.881	−0.366	0.560	24.434	30.470
	(0.428)	(0.218)	(0.175)	(5.895)	(4.360)
hspprt	0.447	−0.043	0.246	−8.666	1.250
	(0.343)	(0.122)	(0.134)	(3.336)	(3.720)
conf	−0.728	−0.065	−0.290	66.946	15.112
	(0.315)	(0.221)	(0.166)	(6.787)	(6.189)
conf²	−0.048	0.008	0.031	−2.412	−1.126
	(0.019)	(0.009)	(0.014)	(0.293)	(0.548)
value	−0.849	0.159	0.141	1.282	9.310
	(0.207)	(0.037)	(0.044)	(1.135)	(1.135)
girlsch	−0.611	−0.119	−0.250	−9.390	−13.906
	(0.639)	(0.256)	(0.309)	(6.886)	(7.624)
coedu	−0.465	−0.012	−0.271	−11.273	−10.829
	(0.459)	(0.197)	(0.235)	(5.274)	(5.896)
disad_11	−0.059	0.130	0.141	−6.118	−1.957
	(0.418)	(0.181)	(0.235)	(4.339)	(5.434)
disad_26	0.170	0.047	0.191	−10.509	−8.031
	(0.453)	(0.204)	(0.246)	(4.757)	(5.742)
disad_m50	−0.403	0.058	−0.490	−25.731	−0.485
	(0.810)	(0.272)	(0.386)	(6.855)	(7.942)
city_u	0.097	−0.229	−0.045	15.936	12.428
	(0.413)	(0.181)	(0.228)	(4.533)	(5.277)
city_sub	0.532	−0.137	−0.315	11.028	2.238
	(0.637)	(0.272)	(0.363)	(6.702)	(8.314)
city_m	−0.067	−0.032	0.047	5.556	2.409
	(0.425)	(0.185)	(0.238)	(4.522)	(5.550)
resshort	0.093	0.003	0.085	0.446	−0.203
	(0.082)	(0.035)	(0.044)	(0.881)	(1.008)
academph	−0.181	0.081	−0.046	1.654	1.372
	(0.107)	(0.044)	(0.056)	(1.045)	(1.295)
dscpn	−0.053	−0.043	−0.016	−0.880	−1.161
	(0.109)	(0.038)	(0.051)	(0.966)	(1.179)

Note. Boldface values indicate that the $95 %$ credible interval does not contain zero. TIMSS = Trends in International Mathematics and Science Study.

The mixing probability model suggests that male students, those having more books at home, and those with more confidence or higher perceived value in science are more likely to be in Class 1. This agrees with the mean outcomes shown in Table 10: On average, these students tend to be good at science even without private lessons, so more science lessons are not very helpful. They may consider spending time in other ways. The descriptive statistics in Table 10 and the estimates of the mixing probability model in Table 11 together help us understand how students in the sample are classified. When new students who are not in the sample want to decide whether to take private science lessons, however, they only have access to their current covariates, so only the mixing probability model is applicable. By plugging in their covariates to the model, they can estimate the probability that they belong to each class, which helps make the decision. In addition, if future outcomes can be roughly predicted, then the descriptive statistics in Table 10 help estimate which latent class they might be in.

It is tricky to interpret the selection models because some covariates might be post-treatment confounders, so that it is hard to tell the direction of causality between the treatment indicator and the covariates. A carefully designed empirical study should take subject matter theory into account and pay more attention to the relationship among variables. Whether students in Class 1 take private science lessons seems only related to their perceived value of science but uncorrelated with any of the other covariates. For students in Class 2, more books at home, more confidence in science, and higher perceived value of science are all associated with higher probability of being treated. Students from schools whose science instruction suffers from resource shortage tend to take private science lessons, which makes much sense: Private lessons serve as a supplement to insufficient school teaching.

The outcome models of the two latent classes show some things in common. Students with more books at home or from schools in urban areas tend to have higher science scores. Students from schools with more economically disadvantaged students or schools with both boys and girls have lower science scores on average, although these covariates are significant for Class 1 only. For students in Class 1, their confidence in science and fathers’ education levels are positively associated with science scores; those having their own room and internet connection have lower science scores on average, which suggests that these conditions might be a distraction. For students in Class 2, their perceived values in science are positively associated with science scores; interestingly, the quadratic function of confidence in the outcome model suggests that more confidence leads to higher science scores for very few unconfident students but leads to lower scores for others.

In brief, we start from the descriptive statistics and the mixing probability model to relate covariates, treatment indicators, and outcomes to latent class membership and then interpret the selection model and the outcome model of each class. In this TIMSS analysis, we find that some variables have quite different coefficient estimates (different significance, signs, or magnitudes) in the selection or the outcome model across the two latent classes. Researchers should pay special attention to such heterogeneity.

Finally, the point estimates of Level 1 error term standard deviation and covariance matrices of random effects are

\begin{matrix} {\hat{σ}}_{1} = 54.374, {\hat{σ}}_{2} = 61.410, \\ {\hat{Σ}}_{u} = [\begin{matrix} 8.386 \\ 10.125 \end{matrix}] [\begin{matrix} 1 & 0.681 \\ 0.681 & 1 \end{matrix}] [\begin{matrix} 8.386 \\ 10.125 \end{matrix}], \\ {\hat{Σ}}_{v} = [\begin{matrix} 0.447 \\ 0.615 \end{matrix}] [\begin{matrix} 1 & 0.620 \\ 0.620 & 1 \end{matrix}] [\begin{matrix} 0.447 \\ 0.615 \end{matrix}] . \end{matrix}

The outcome model and the selection model show similar patterns on random effects: The two latent classes’ random intercepts are moderately positively correlated, and Class 2 has greater variation than Class 1. Class 2 also has greater variation in the error term. We conclude that students in Class 2 are more heterogeneous than those in Class 1.

Discussion

In this article, we propose a latent class model for estimating heterogenous treatment effects under the multilevel setting. The model relies on three pieces of information to classify Level 1 units into latent classes: different distributions of covariates, different selection models, and different outcome models across latent classes. Instead of proceeding with each piece of information sequentially, where errors from one model will be directly carried over to the next model and gradually accumulate, we adopt a Bayesian joint procedure that takes all the information into account simultaneously to estimate this model. The Bayesian approach also gives additional information like credible intervals for the parameters and LOOIC for model comparison. The simulation study results suggest that the Bayesian procedure gives better estimates than the EM algorithm, which is based on the same model, and LOOIC provides sufficient information for determining the correct number of latent classes. For model checking, we also compute the doubly robust ATE estimates to examine possible functional form misspecification and underspecification of the number of latent classes. Finally, we apply the procedure to a real-world data set to demonstrate how to use our proposed method and interpret the results. The empirical example clearly shows the necessity of using our latent class model instead of simply ignoring the heterogeneity among individuals.

A main advantage of our approach is that a large amount of information about the estimation, including the posterior standard deviations of parameters and the posterior probability that each individual belongs to each latent class, is provided by the estimation procedure automatically. With the three pieces of the estimated model, mixing probabilities, selection processes, and outcomes, researchers can easily interpret the results and relate them to existing subject matter theory. The posterior probability of each individual shows the uncertainty in classification and can be helpful for deciding the optimal treatment for each individual through benefit–cost analysis.

It is also worth noting that our proposed method is very flexible. In this article, we try to model the joint distribution of parameters conditional on the covariates and the cluster membership of Level 1 units. We only try to model the selection processes and the outcomes, while how the covariates are generated and how Level 1 units are assigned to Level 2 clusters are taken as given. It is straightforward to extend the model, so that these two processes can be incorporated if one has good models for them. Although we only consider subgroup-specific ATEs in this study, it is also possible to estimate treatment effects conditional on observed covariates by simply including interaction terms between the treatment indicator and covariates in the outcome models. The idea of subgroup-specific models for conditional treatment effects can give researchers even more detailed information on heterogeneous treatment effects.

There are some possible future directions that may help improve our current method. First, we found that the Bayesian procedure has a hard time converging in a few cases when we specify a too large number of latent classes. Nonconvergence itself might be an indicator of overfit, but it needs to be verified. Second, how to better use the doubly robust estimator to detect possible model misspecification is still not clear. In this article, we use the $50 %$ credible intervals of doubly robust estimates for model checking, and the results look promising. However, this approach is not perfect because the choice of $50 %$ is somewhat arbitrary and lacks theoretical support, and double robustness does not help improve estimation but only works as additional evidence on model fit. Incorporating the double robustness into the Bayesian estimation procedure is difficult because we cannot write it into a single log-likelihood function, but promising new approaches have been proposed to circumvent this problem. Luo et al. (2022) proposed a new approach where a loss function which involves the doubly robust estimator replaces the log-likelihood function. A possible future direction is to build mixture models based on this idea. Third, as Vehtari et al. (2019) have pointed out, LOOIC may fail under multilevel settings because conditional exchangeability does not hold. Our simulation shows that LOOIC, even without theoretical justification, still gives the correct number of latent classes in most cases. How to better compare models for clustered data is still under study.

Supplemental Material

Supplemental Material, sj-pdf-1-jeb-10.3102_10769986221115446 - Estimating Heterogeneous Treatment Effects Within Latent Class Multilevel Models: A Bayesian Approach

Supplemental Material, sj-pdf-1-jeb-10.3102_10769986221115446 for Estimating Heterogeneous Treatment Effects Within Latent Class Multilevel Models: A Bayesian Approach by Weicong Lyu, Jee-Seon Kim and Youmi Suk in Journal of Educational and Behavioral Statistics

Footnotes

Authors’ Note

This research was performed using the compute resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation, and is an active member of the OSG Consortium, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research and/or authorship of this article.

ORCID iDs

Weicong Lyu

Jee-Seon Kim

Youmi Suk

References

Athey

Imbens

(2016). Recursive partitioning for heterogeneous causal effects. Proceedings of the National Academy of Sciences, 113(27), 7353–7360. https://doi.org/10.1073/pnas.1510489113

Bishop

C. M.

(2006). Pattern recognition and machine learning. Springer.

Clogg

C. C.

(1995). Latent class models. In Arminger

Clogg

C. C.

Sobel

M. E.

(Eds.), Handbook of statistical modeling for the social and behavioral sciences (pp. 311–359). Springer.

Crump

R. K.

Hotz

V. J.

Imbens

G. W.

Mitnik

O. A.

(2009). Dealing with limited overlap in estimation of average treatment effects. Biometrika, 96(1), 187–199. https://doi.org/10.1093/biomet/asn055

Fuentes

Lüdtke

Robitzsch

(2021). Causal inference with multilevel data: A comparison of different propensity score weighting approaches. Multivariate Behavioral Research, 1–24. https://doi.org/10.1080/00273171.2021.1925521

Gelman

Hullman

Wlezien

Morris

G. E.

(2020). Information, incentives, and goals in election forecasts. Judgment and Decision Making, 15(5), 863–880.

Gelman

Lee

Guo

(2015). Stan: A probabilistic programming language for Bayesian inference and optimization. Journal of Educational and Behavioral Statistics, 40(5), 530–543. https://doi.org/10.3102/1076998615606113

Heinen

(1996). Latent class and discrete latent trait models: Similarities and differences. Sage Publications.

Hill

J. L.

(2011). Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 217–240. https://doi.org/10.1198/jcgs.2010.08162

10.

Holland

P. W.

(1986). Statistics and causal inference. Journal of the American Statistical Association, 81(396), 945–960. https://doi.org/10.2307/2289064

11.

Hong

Raudenbush

S. W.

(2006). Evaluating kindergarten retention policy: A case study of causal inference for multilevel observational data. Journal of the American Statistical Association, 101(475), 901–910. https://doi.org/10.1198/016214506000000447

12.

Hong

Raudenbush

S. W.

(2013). Heterogeneous agents, social interactions, and causal inference. In Morgan

S. L.

(Ed.), Handbook of causal analysis for social research (pp. 331–352). Springer. https://doi.org/10.1007/978-94-007-6094-3_16

13.

Imai

Ratkovic

(2013). Estimating treatment effect heterogeneity in randomized program evaluation. The Annals of Applied Statistics, 7(1), 443–470. https://doi.org/10.1214/12-aoas593

14.

Imai

Strauss

(2011). Estimation of heterogeneous treatment effects from randomized experiments, with application to the optimal planning of the get-out-the-vote campaign. Political Analysis, 19(1), 1–19. https://doi.org/10.1093/pan/mpq035

15.

Imbens

G. W

Rubin

D. B

. (2015). Causal inference in statistics, social, and biomedical sciences. Cambridge University Press. https://doi.org/10.1017/CBO9781139025751

16.

Kim

H. J.

Nehus

E. J.

Kim

M.-O.

(2019). Estimating heterogeneous treatment effects for latent subgroups in observational studies. Statistics in Medicine, 38(3), 339–353. https://doi.org/10.1002/sim.7970

17.

Kim

J.-S.

Lim

W.-C.

Steiner

P. M.

(2017). Causal inference with observational multilevel data: Investigating selection and outcome heterogeneity. In van der Ark

L. A.

Wiberg

Culpepper

S. A.

Douglas

J. A.

Wang

W.-C.

(Eds.), Quantitative psychology (pp. 287–308). Springer. https://doi.org/10.1007/978-3-319-56294-0_26

18.

Kim

J.-S.

Steiner

P. M.

(2015). Multilevel propensity score methods for estimating causal effects: A latent class modeling strategy. In van der Ark

L. A.

Bolt

D. M.

Wang

W.-C.

Douglas

J. A.

Chow

S.-M.

(Eds.), Quantitative psychology research (pp. 293–306). Springer. https://doi.org/10.1007/978-3-319-19977-1_21

19.

Kim

J.-S.

Steiner

P. M.

Lim

W.-C.

(2015). Mixture modeling methods for causal inference with multilevel data. In Harring

J. R.

Stapleton

L. M.

Beretvas

S. N.

(Eds.), Advances in multilevel modeling for educational research: Addressing practical issues found in real-world applications (pp. 335–359). IAP.

20.

Künzel

S. R.

Walter

S. J.

Sekhon

J. S.

(2019). Causaltoolbox—Estimator stability for heterogeneous treatment effects. Observational Studies, 5(2), 105–117. https://doi.org/10.1353/obs.2019.0005

21.

Leite

. (2016). Practical propensity score methods using R. Sage Publications. https://doi.org/10.4135/9781071802854

22.

Lindsay

Clogg

C. C.

Grego

(1991). Semiparametric estimation in the Rasch model and related exponential response models, including a simple latent class model for item analysis. Journal of the American Statistical Association, 86(413), 96–107. https://doi.org/10.1080/01621459.1991.10475008

23.

Loh

W. W.

Kim

J.-S

. (2020). Evaluating sensitivity to classification uncertainty in subgroup effect analyses. https://doi.org/10.31234/osf.io/nea94

24.

Luo

Stephens

D. A.

Graham

D. J.

McCoy

E. J

. (2022). Bayesian doubly robust causal inference via loss functions. https://doi.org/10.48550/arXiv.2103.04086

25.

Luque-Fernandez

M. A.

Schomaker

Rachet

Schnitzer

M. E.

(2018). Targeted maximum likelihood estimation for a binary treatment: A tutorial. Statistics in Medicine, 37(16), 2530–2546. https://doi.org/10.1002/sim.7628

26.

Martin

M. O.

Mullis

I. V. S.

Hooper

(2016). Methods and Procedures in TIMSS 2015. TIMSS & PIRLS International Study Center, Boston College.

27.

McLachlan

G. J.

Lee

S. X.

Rathnayake

S. I.

(2019). Finite mixture models. Annual Review of Statistics and Its Application, 6, 355–378. https://doi.org/10.1146/annurev-statistics-031017-100325

28.

Moerkerke

Loeys

Vansteelandt

(2015). Structural equation modeling versus marginal structural modeling for assessing mediation in the presence of posttreatment confounding. Psychological Methods, 20(2), 204–220. https://doi.org/10.1037/a0036368

29.

Neyman

Iwaszkiewicz

(1935). Statistical problems in agricultural experimentation. Supplement to the Journal of the Royal Statistical Society, 2(2), 107–180. https://doi.org/10.2307/2983637

30.

Nie

Wager

(2021). Quasi-oracle estimation of heterogeneous treatment effects. Biometrika, 108(2), 299–319. https://doi.org/10.1093/biomet/asaa076

31.

Park

D. K.

Gelman

Bafumi

(2006). State-level opinions from national surveys: Poststratification using multilevel logistic regression. In Cohen

J. E.

(Ed.), Public opinion in state politics (pp. 209–228). Stanford University Press. https://doi.org/10.11126/stanford/9780804753005.003.0011

32.

Pearl

(2017). Detecting latent heterogeneity. Sociological Methods & Research, 46(3), 370–389. https://doi.org/10.1177/0049124115600597

33.

Rosenbaum

P. R.

Rubin

D. B.

(1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41–55. https://doi.org/10.1093/biomet/70.1.41

34.

Rubin

D. B.

(1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5), 688. https://doi.org/10.1037/h0037350

35.

Rubin

D. B.

(1978). Bayesian inference for causal effects: The role of randomization. The Annals of Statistics, 6(1): 34–58. https://doi.org/10.1214/aos/1176344064

36.

Rubin

D. B.

(1986). Statistics and causal inference: Comment: Which ifs have causal answers. Journal of the American Statistical Association, 81(396), 961–962. https://doi.org/10.2307/2289065

37.

Scharfstein

D. O.

Rotnitzky

Robins

J. M.

(1999). Adjusting for nonignorable drop-out using semiparametric nonresponse models. Journal of the American Statistical Association, 94(448), 1096–1120. https://doi.org/10.2307/2669923

38.

Shahn

Madigan

(2017). Latent class mixture models of treatment effect heterogeneity. Bayesian Analysis, 12(3), 831–854. https://doi.org/10.1214/16-ba1022

39.

Spilt

J. L.

Koot

J. M.

van Lier

P. A.

(2013). For whom does it work? Subgroup differences in the effects of a school-based universal prevention program. Prevention Science, 14(5), 479–488. https://doi.org/10.1007/s11121-012-0329-7

40.

Suk

Kim

J.-S.

Kang

(2021). Hybridizing machine learning methods and finite mixture models for estimating heterogeneous treatment effects in latent classes. Journal of Educational and Behavioral Statistics, 46(3), 323–347. https://doi.org/10.3102/1076998620951983

41.

Vehtari

Gelman

Gabry

(2017). Practical Bayesian model evaluation using leave-one-out cross-validation and waic. Statistics and Computing, 27(5), 1413–1432. doi: https://doi.org/10.1007/s11222-016-9696-4

42.

Vehtari

Gelman

Simpson

Carpenter

Bürkner

P.-C.

(2021). Rank-normalization, folding, and localization: An improved b R for assessing convergence of mcmc (with discussion). Bayesian Analysis, 16(2), 667–718. https://doi.org/10.1214/20-BA1221

43.

Vehtari

Simpson

D. P.

Yao

Gelman

(2019). Limitations of “limitations of Bayesian leave-one-out cross-validation for model selection.” Computational Brain & Behavior, 2(1), 22–27. https://doi.org/10.1007/s42113-018-0020-6

44.

Vermunt

J. K.

Magidson

(2002). Latent class cluster analysis. In Hagenaars

J. A.

McCutcheon

A. L.

(Eds.), Applied latent class analysis (pp. 89–106). Cambridge University Press. https://doi.org/10.1017/cbo9780511499531.004

45.

Vrieze

S. I.

(2012). Model selection and psychological theory: A discussion of the differences between the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). Psychological Methods, 17(2), 228. https://doi.org/10.1037/a0027127

46.

Wager

Athey

(2018). Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523), 1228–1242. https://doi.org/10.1080/01621459.2017.1319839

47.

Wendling

Jung

Callahan

Schuler

Shah

Gallego

(2018). Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases. Statistics in Medicine, 37(23), 3309–3324. https://doi.org/10.1002/sim.7820

48.

White

Murphy

T. B.

(2014). BayesLCA: An R package for Bayesian latent class analysis. Journal of Statistical Software, 61(13), 1–28. https://doi.org/10.18637/jss.v061.i13

49.

Willke

R. J.

Zheng

Subedi

Althin

Mullins

C. D.

(2012). From concepts, theory, and evidence of heterogeneity of treatment effects to methodological approaches: A primer. BMC Medical Research Methodology, 12(1), 185. https://doi.org/10.1186/1471-2288-12-185

50.

Xie

Brand

J. E.

Jann

(2012). Estimating heterogeneous treatment effects with observational data. Sociological Methodology, 42(1), 314–347. https://doi.org/10.1177/0081175012452652

51.

Zeldow

III, V. L

Roy

. (2019). A semiparametric modeling approach using Bayesian additive regression trees with an application to evaluate heterogeneous treatment effects. The Annals of Applied Statistics, 13(3), 1989–2010. https://doi.org/10.1214/19-aoas1266

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.10 MB