{'fraction_non_alphanumeric': 0.03370912209613066, 'fraction_numerical': 0.012044299156028986, 'mean_word_length': 4.471572361262242, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 49, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating Heterogeneous Causal Effects of High-Dimensional Treatments: Application to Conjoint Analysis

Estimation of heterogeneous treatment effects is an active area of research in causal inference. Most of the existing methods, however, focus on estimating the conditional average treatment effects of a single, binary treatment given a set of pre-treatment covariates. In this paper, we propose a method to estimate the heterogeneous causal effects of high-dimensional treatments, which poses unique challenges in terms of estimation and interpretation. The proposed approach is based on a Bayesian mixture of regularized logistic regressions to identify groups of units who exhibit similar patterns of treatment effects. By directly modeling cluster membership with covariates, the proposed methodology allows one to explore the unit characteristics that are associated with different patterns of treatment effects. Our motivating application is conjoint analysis, which is a popular survey experiment in social science and marketing research and is based on a high-dimensional factorial design. We apply the proposed methodology to the conjoint data, where survey respondents are asked to select one of two immigrant profiles with randomly selected attributes. We find that a group of respondents with a relatively high degree of prejudice appears to discriminate against immigrants from non-European countries like Iraq. An open-source software package is available for implementing the proposed methodology.


Introduction
Over the past decade, a number of researchers have exploited modern machine learning algorithms and proposed new methods to estimate heterogeneous treatment effects using experimental data. They include tree-based methods (e.g., Imai and Strauss, 2011;Athey and Imbens, 2016;Wager and Athey, 2018;Hahn, Murray and Carvalho, 2020), regularized regressions (e.g., Imai and Ratkovic, 2013;Tian et al., 2014;Künzel et al., 2019), ensemble methods (e.g., van der Laan and Rose, 2011;Grimmer, Messing and Westwood, 2017), and frameworks that allow for the use of generic machine learning methods (e.g., Chernozhukov et al., 2019;Imai and Li, 2021). This methodological development, however, has largely been confined to settings with a single, binary treatment variable. Indeed, there exists little prior research that models heterogeneous effects of high-dimensional treatments, which pose unique challenges in terms of computation and interpretation.
In this paper, we consider the estimation of heterogeneous causal effects based on the data from a randomized experiment, in which the treatments of interest are high-dimensional with the number of possible treatment combinations exceeding the sample size. We address the challenge of effectively summarizing the complex patterns of heterogeneous treatment effects that are induced by the interactions among the treatments themselves as well as the interactions between the treatments and unit characteristics. We use an interpretable machine learning algorithm based on a Bayesian mixture of (regularized) logistic experts model to identify groups of units who exhibit similar patterns of treatment effects. To further facilitate the interpreta-tion of the results, we directly model cluster membership using unit characteristics.
Thus, the proposed methodology allows one to explore the unit characteristics that are associated with different patterns of treatment effects.
Our motivating application is conjoint analysis, which is a popular survey experimental methodology in social sciences and marketing research (e.g., Hainmueller, Hopkins and Yamamoto, 2014;Rao, 2014). Conjoint analysis is a variant of factorial designs (Dasgupta, Pillai and Rubin, 2015) with a large number of factorial treatments-so large that typically not all possible treatments are seen. Under the most commonly used "forced-choice" design, respondents are asked to evaluate a pair of profiles whose attributes are randomly selected based on factorial variables with several levels. In the specific experiment we reanalyze, the original authors used a conjoint analysis to measure immigration preferences by presenting each survey respondent with several pairs of immigrant profiles with varying attributes including education, country of origin, and job experience (Hainmueller and Hopkins, 2015).
For each pair, the respondent was asked to choose one profile they prefer. The authors then analyzed the resulting response patterns to understand which immigrant characteristics play a critical role in forming the immigration preferences of American citizens.
In the methodological literature on factorial designs and conjoint analysis, researchers have almost exclusively focused upon the estimation of marginal effects, which represent the average effect of one factor level relative to another level of the same factor averaging over the randomization distribution of the remaining factors (Hainmueller, Hopkins and Yamamoto, 2014;Dasgupta, Pillai and Rubin, 2015). Although some have explored the estimation of interaction effects among the factorial treatments (e.g., Dasgupta, Pillai and Rubin, 2015;Egami and Imai, 2019;De la Cuesta, Egami and Imai, 2022), there exists little methodological research that in-vestigates how to estimate heterogeneous treatment effects of high-dimensional treatments. As a result, in conjoint analysis, a dominant approach to the estimation of heterogeneous treatment effects is based on subgroup analysis, in which researchers estimate the marginal effects of interest using a particular subset of respondents (e.g., Hainmueller and Hopkins, 2015;Newman and Malhotra, 2019). Unfortunately, this approach often results in low statistical power and may suffer from multiple testing problems. In addition, subgroup analysis that focuses on marginal effects fails to capture complex patterns of treatment effects based on interactions that may be of scientific interest.
To overcome this challenge, we propose a Bayesian mixture of (regularized) logistic experts model by building a model that draws together two distinct strands of methodological research. First, a number of scholars have developed a finite mixture of regularized regressions (e.g., Khalili and Chen, 2007;Städler, Bühlmann and Van De Geer, 2010), with some allowing covariates to predict cluster membership under a mixture of experts framework (Khalili, 2010). The marketing literature using conjoint analysis has long applied similar mixture models to analyzing heterogeneity (e.g., Gupta and Chintagunta, 1994;Andrews, Ainslie and Currim, 2002), although their treatments are typically low dimensional and do not require regularization. Using a mixture of experts gives a small number of treatment effect patterns that can be easily interpreted while also incorporating a rich set of moderators to characterize how they relate to different patterns of treatment effects.
Second, to analyze the data from a high-dimensional factorial experiment, we draw on a growing literature that shows how to regularize factor variables by grouping and fusing levels together rather than shrinking levels to zero (e.g., Bondell and Reich, 2009;Post and Bondell, 2013;Stokell, Shah and Tibshirani, 2021). This facilitates the interpretation of empirical findings by identifying a set of factor levels that character-ize treatment effect heterogeneity. Our proposed model accommodates interactions, respecting a hierarchical structure where the main effects may be fused only if the interactions themselves are also fused (Yan and Bien, 2017). For regularization, we use an 2 norm for the group overlapping lasso penalty for computational tractability although existing work has relied on a ∞ norm penalty to induce sparsity (Post and Bondell, 2013;Egami and Imai, 2019).
To fit this model, we use a Bayesian approach and apply an Expectation-Maximization (EM) algorithm through data augmentation to find a posterior mode (Dempster, Laird and Rubin, 1977;Meng and van Dyk, 1997). We exploit the representation of 1 and 2 penalties as a mixture of Gaussians and derive a tractable EM algorithm (see, e.g., Figueiredo 2003;Polson and Scott 2011;Ratkovic and Tingley 2017;Goplerud 2021, for some uses of EM algorithm for sparse models). This formulation allows for the application of commonly available methods to accelerate EM algorithms while maintaining monotonic convergence such as SQUAREM (Varadhan and Roland, 2008).
The rest of the paper is organized as follows. In Section 2, we discuss the motivating application, which is a conjoint analysis of American citizens' preferences regarding immigrant features. We also briefly describe a methodological challenge to be addressed. In Section 3 , we present our proposed methodology. In Section 4, we show our method performs well in a realistic numerical simulation. In Section 5, we apply this methodology and reanalyze the data from the motivating conjoint analysis.
Section 6 concludes with a discussion.

Motivating Application: Conjoint Analysis of Immigration Preferences
Our motivating application is a conjoint analysis of American immigration preferences. In this section, we introduce the experimental design and discuss the results of previous analyses that motivate our methodology for estimating heterogeneous treatment effects.

The Experimental Design
In an influential study, Hainmueller and Hopkins (2015) use conjoint analysis to estimate the effect of immigrant attributes on preferences for admission to the United States 1 . The authors conduct an online survey experiment using a sample of 1,407 American adults. Each survey respondent assessed five pairs of immigrant profiles with randomly selected attributes. For each pair, a respondent was asked to choose one of the two immigrant profiles they preferred to admit to the United States.
The attributes of immigrant profiles used in this factorial experiment are gender, education, employment plans, job experience, profession, language skills, country of origin, reasons for applying, and prior trips to the United States. For completeness, these factors and their levels are reproduced as Table 1 of Appendix A. In total, there exist over 1.4 million possible profiles, implying more than 2 × 10 12 possible comparisons of two profiles that are possible in the experiment. It is clear that with 1,407 respondents, even though each respondent performs five comparisons, not all possible profiles can be included.
The levels of each factor variable were independently randomized to yield one immigration profile. Randomization was subject to some restrictions such that pro- fession and education factors result in sensible pairings (e.g., ruling out doctors with with less than two-years of college education) and immigrants whose reason for applying is persecution are from Iraq, Sudan, Somalia, or China. The ordering of the attributes was also randomized for each respondent. The experiment additionally collected data on the respondents, including various demographic information, partisanship, attitudes towards immigration, and ethnocentrism. A rating for each immigrant profile was also recorded, but that metric is not the focus of our analysis.

Heterogeneous Treatment Effects
In the original article, Hainmueller and Hopkins (2015) conducted their primary analysis based on linear regression model where the unit of analysis is an immigrant profile (rather than a pair) and the outcome variable is an indicator for whether a given profile was chosen. The predictors of the model include the indicator variable for each immigrant attribute. The model also includes the interactions between education and profession, as well as between country of origin and reasons for applying, in order to account for the restricted randomization scheme mentioned above. Finally, the standard errors are clustered by respondent.
As formalized in Hainmueller, Hopkins and Yamamoto (2014), the regression coefficient represents the average marginal component effect (AMCE) of each attribute averaging over all the other attributes including those of the other profile in a given pair. The left plot of Figure 1 reproduces the estimated overall AMCEs of country of origin where the baseline category is Germany. There is little country effect with the exception of Iraq, which negatively affects the likelihood of being preferred by a respondent.
Beyond the AMCEs, these authors and others including Newman and Malhotra (2019) have explored the heterogeneous treatment effects among respondents by conducting many sub-group analyses for respondent characteristics including parti-  (2015), whereas the right two plots present the results of a subgroup analysis based on the level of respondent's Hispanic prejudice score.
sanship and level of education. Here, we follow Newman and Malhotra (2019) and focus on a measure of ethnic and racial prejudice by utilizing the Hispanic prejudice score, which represents a 0-100 point feeling thermometer towards Hispanics. We dichotomize this score in the same way as Newman and Malhotra (2019), labeling those with score below 50 as "high prejudice," and all others as "low prejudice." The results for high and low prejudice groups are shown in the right two plots of Figure 1 We now turn to our methodology which is designed to address this challenge.

The Proposed Methodology
In this section, we describe the proposed methodology for estimating heterogeneous effects of high-dimensional treatments. To simplify the exposition and notation, we focus on a general factorial design. This design corresponds to conjoint analysis with a single task per person, where there is only one treatment or profile assessed rather than a comparison to two profiles, and the complete and independent randomization of all levels within each factor. The extensions to more realistic conjoint analyses are immediate and will be discussed and applied in Section 5.

Set Up
Suppose that we have a simple random sample of N units. Consider a factorial design with J factors where each factor j ∈ {1, · · · , J} has L j ≥ 2 levels. The treatment variable for unit i, denoted by T i , is a J-dimensional vector of random variables, each of which represents the assigned level of the corresponding factor variable. For example, the jth element of this random vector T ij ∈ {0, 1, 2, . . . , L j − 1} represents the level of factor j which is assigned to unit i.
Following Dasgupta, Pillai and Rubin (2015), we define the potential outcome for unit i as Y i (t) where t ∈ T represents the realized treatment with T representing the support of the randomization distribution for T i . Then, the observed outcome is given by Y i = Y i (T i ). The notation implicitly assumes no interference between units (Rubin, 1980). In this paper, for the sake of concreteness, we focus on the binary outcome Y i ∈ {0, 1}. Extensions to non-binary outcomes are straightforward. Lastly, we observe a vector of p x pre-treatment covariates for each unit i and denote it by X i . All together, we observe (Y i , T i , X i ) for each unit i.
To illustrate the notation, consider a simplified version of our motivating example where each participant i observes a single immigrant profile and must decide whether to support that immigrant's admission or not. Then T i is a vector indicating the level participant i sees for each of the nine immigrant attributes. Y i is an indicator for whether participant i chooses to support admission for the immigrant they are presented with. X i is a vector of covariates for participant i that we believe might moderate the treatment effect (i.e., we believe there may be treatment effect heterogeneity based on the X i ). In our example, X i included political party, education, demographics of their ZIP code, ethnicity, and Hispanic prejudice score (see Section 5.1 for more details on how these variables were defined).
The randomness in our data, (Y i , T i , X i ) comes from two sources: random sampling of units into the study and random assignment of units to treatments. For simplicity, we assume units are sampled via simple random sampling (though this can be relaxed and survey weights incorporated into our method). We further assume that the design used to randomize the assignment of of treatment is such that for each i (this will hold for common designs where all units have the same probabilities of receiving each treatment). The exact mode of randomization will determine the distribution of T i . In many practical applications of conjoint analysis, researchers independently and uniformly randomize each factor. However, in some cases including our application, researchers may exclude certain unrealistic combinations of factor levels (e.g., doctor without a college degree), leading to the dependence between factors. In all cases, researchers have complete knowledge of the randomization distribution of the factorial treatment variables.
Based on random sampling and random treatment assignment alone, we can conduct valid inference for treatment effects of interest using basic regression or difference of means. However, if we wish to explore heterogeneity in responses to treatments, a model-based approach is useful. We next introduce our model, which will allow us to explore heterogeneous effects in a principled manner while also handling the high-dimensional nature of the data.

Modeling Heterogeneous Treatment Effects
The most basic causal quantity of interest is the average marginal component effect (AMCE; Hainmueller, Hopkins and Yamamoto, 2014), which is defined for any given factor j as, where l = l ∈ T j with T j representing the support of the randomization distribution for T j . The expectation in Equation (1)  . This means we will end up with a set of interpretable clusters. For our empirical example, we might recover two-cluster results similar to those shown in Figure 1, where, for example, cluster 1 corresponds to highly prejudiced individuals and cluster 2 corresponds to low prejudice individuals. The model consists of two parts. First, for each cluster, a regularized logistic regression is used to model the outcome variable as a function of treatments T i where an ANOVA-style sum-to-zero constraint is imposed separately for each factor. Regularization is used to facilitate merging of different levels within each factor. This modeling strategy identifies a relatively small number of treatment combinations while avoiding the specification of a baseline level for each factor (Egami and Imai, 2019). Second, we model the probability of unit i's assignment to each cluster using a set of p x unit characteristics X i . We refer to X i as "moderators." These moderators describe the characteristics of units that belong to each cluster. All together, the proposed modeling strategy allows researchers to identify interpretable patterns of heterogeneous effects across treatment combinations and units. Our model is conditional on the X i and T i .
Formally, let K be the number of clusters. Then, the model is defined as follows, where i = 1, 2, . . . , N and for k = 1, 2, . . . , K, we use the logistic regression to model the outcome variable and the multinomial logit model for the cluster membership, , and π k ( For identification, we set φ 1 = 0. We use the following linear equation for ψ k (T i ) where we include both main effects and two-way interaction effects with a common intercept µ shared across all clusters, for each k = 1, 2, . . . , K where t i is the vector of indicators, 1{T ij = l} and 1{T ij = l, T ij = l }, and β k is a stacked column vector containing all coefficients for cluster k. We use the following ANOVA-type sum-to-zero constraints, β j kl = 0, and for each j, j = 1, 2, . . . , J with j > j, and l = 0, 1, . . . , L j .
We write these constraints compactly as, Each row of C β k corresponds to one of the sum-to-zero constraints given in Equation (4). For example, in a model without interactions, C takes the following where 0 p and 1 p denote the column vectors of zeros and ones, respectively, with length p.

Sparsity-induced Prior
Given the high dimensionality of this model, we use a sparsity-inducing prior. For example, in our application, we have 315 β parameters for each cluster. In factorial experiments, it is desirable to regularize the model such that certain levels of each factor are fused together when their main effects and all interactions are similar (Post and Bondell, 2013;Egami and Imai, 2019). For example, we would like to fuse levels l 1 and l 2 of factor j if β j l 1 ≈ β j l 2 and β jj l 1 l ≈ β jj l 2 l for all other factors j and all of its levels l .
We encourage such fusion by applying a structured sparsity approach of Goplerud (2021) that generalizes the group and overlapping group LASSO (e.g., Yuan and Lin, 2006;Yan and Bien, 2017) while allowing positive semi-definite constraint matrices.
For computational tractability, we use 2 norm instead of the ∞ norm, which is used in GASH-ANOVA (Post and Bondell, 2013).
The first three terms encourages the pairwise fusion of the levels of factor one whereas the fourth encourages the fusion of the two levels of factor two. For compact notation, the penalty can also be written using the sum of Euclidean norms of quadratic forms, where F 1 , F 2 , F 3 are appropriate positive semi-definite matrices to encourage the fusion of the pairs of levels in factor one and F 4 encourages the fusion of the two factors in factor two, and β = [β 1 0 β 1 1 β 1 2 β 2 0 β 2 1 β 12 00 β 12 10 β 12 20 β 12 01 β 12 11 β 12 21 ] . Note that the sum to zero constraints make this type of fusion of factors together sensible for sparsity.
We generalize this formulation to an arbitrary number of factors and factor levels.
For each factor that contains L j levels, we have L j 2 penalty matrices to encourage pairwise fusion. Imposing additional constraints, e.g. for ordered factors, is a simple represent the total number of penalty matrices. For g = 1, 2, . . . , G, we use F g to denote a penalty matrix such that β F g β is equivalent to the 2 norm on the vector of differences between all main effects and interactions containing a main effect.
We cast this penalty as the following Bayesian prior distribution so that our inference is based on the posterior distribution given the conditional prior on β k , whereπ k = N i=1 π k (X i )/N and m = rank ([F 1 , · · · , F G ]). We follow existing work in allowing the penalty on the treatment effects β k to be scaled by the clustermembership sizeπ k when γ = 1 (Khalili and Chen, 2007;Städler, Bühlmann and Van De Geer, 2010). On the other hand, when γ = 0 theπ k disappears, implying no use of the X i in the prior. We note that the prior on p(β | {φ k } K k=2 ) is guaranteed to be proper when all pairwise fusions are encouraged by {F g } G g=1 , although in other circumstances it may be improper (Goplerud, 2021). Appendix B provides additional details. Following Zahid and Tutz (2013), we use a normal prior distribution for the coefficients for the moderators.
The resulting regularization is invariant to the choice of baseline category φ 1 = 0, which is the first row of the K × p x coefficient matrix φ. The prior distribution is given by, where Σ φ is a (K − 1) × (K − 1) matrix with [Σ φ ] kk = (K − 1)/K if k = k and [Σ φ ] kk = −1/K otherwise. We set σ 2 φ to 1/4 for a relatively diffuse prior.
As noted in a recent survey, "ensuring generic identifiability for general [mixture of expert] models remains a challenging issue" (Gormley and Frühwirth-Schnatter, 2019, p. 294). Although mixtures with a Bernoulli outcome variable are generally unidentifiable, several aspects of our methodology are expected to alleviate the identifiability problem. First, a typical conjoint analysis has repeated observations per unit i (Grün and Leisch, 2008). Second, our model is a mixture of experts rather than a mixture model (Jiang and Tanner, 1999). Third, our covariates are randomized.
Finally, our model regularizes the coefficients through an informative prior. Unfortunately, a formal identifiability analysis of our model is beyond the scope of this paper.
Our simulation analysis (Section 4), however, shows that our model can accurately recover the coefficients in a realistic setting for applied research. It is also possible to use a bootstrap-based procedure to examine the identifiability issue in a specific setting (Grün and Leisch, 2008). Equations (2), (6), and (7), where we collect all model parameters as θ.

Estimation and Inference
For now, we assume the value of regularization parameter λ is fixed, although we discuss this issue in Section 3.5. The linear constraints on β k given in Equation (A15) still hold but are suppressed for notational simplicity; we return to this point at the end of this subsection.
Algorithm 1 summarizes our approach to finding a maximum of Equation (8).
Each iteration of our AECM algorithm involves two "cycles" where the data augmentation scheme enables iterative updating of the treatment effect parameters β and moderators φ.
To update β, our data augmentation strategy requires three types of missing data.
and Z i ∼ Multinomial (1, π i ) with the kth element of π equal to π k (X i ). Note that β only enters Equation (9) via a quadratic form.
The first cycle of the AECM algorithm involves, therefore, maximizing or improving the following Q-function with respect to β given θ (t) .
where all expectations are taken over the conditional distribution of the missing data given the current parameter estimates. We note that the E-Step involves computing ) which factorizes into, respectively, a collection of Polya-Gamma, categorical, and Inverse-Gaussian random variables (Appendix C.2 provides the full conditional distributions and all relevant expectations). Thus, maximizing Q β reduces to a linear ridge regression for β.
The constraint implies that β k belongs to the null space of C (Lawson and Hanson, 1974). Let B C represent a basis of the null space of C . Rather than performing inference directly on β k in the constrained space, we consider unconstrained infer- Similarly, we re-define F g as B C F g B C . Once the algorithm converges, the constrained parameters β k are obtained noting β k = B C β k . Appendix C.1 provides a full derivation of the removal of the linear constraints.
To update the moderator parameters φ, we use the second cycle of the AECM algorithm where only the Z i are treated as missing data. The E-step involves recom- , given the updates in the first cycle. The implied Q-function is shown below, where π k (X i ) andπ k = N i=1 π k (X i )/N are functions of φ k . Note that if γ = 0, this simplifies to a (weighted) multinomial logistic regression with E(z ik ) as the outcome.
We perform the M -Step by using a few steps of a standard optimizer (e.g., L-BFGS-B) to increase Q φ and thus obtain φ (t+1) .

Additional Considerations
Since fitting the proposed model is computationally expensive, we use an information criterion approach based on BIC, rather than cross validation, to select the value of Lastly, our experience suggests that the proposed penalty function, which consists of overlapping groups, often finds highly sparse solutions. Appendix D.4 details the integration of the latent overlapping group formulation of Yan and Bien (2017) into our framework to address this issue.
Once the model parameters are estimated, we can compute quantities of interest such as the AMCEs, defined in Equation (1). We do this separately for each cluster, such that δ jk (l, l ) is the AMCE for factor j, changing from level l to l in cluster k. Our estimator is the average of the estimated difference in predicted responses when changing from level l to l of factor j, where the average is taken over the empirical distribution of the assignment on the other factors. This estimation is described in more detail in Appendix E: Appendix E.1 discusses factorial designs (with or without randomization restrictions) whereas Appendix E.2 discusses forcedchoice conjoint designs (with or without randomization restrictions). We can use the empirical distribution here because treatment is randomly assigned.
To quantify the uncertainty of the parameter estimates, we rely on a quadratic approximation to the log posterior distribution. To ensure its differentiability, we follow a standard approach in the regularized regression literature (e.g., Fan and Li, 2001) and fuse pairwise factor levels that are sufficiently close together. Appendix F describes this process, with Appendix F.1 and F.2 deriving the Hessian of the logposterior using Louis (1982)'s method; we then use the Delta method, as described in Appendix F.3, for inference on other of quantities of interest, e.g. the AMCE.

Simulations
We explore the performance of our method using a simple but realistic simulation study. Specifically, we consider the case of a conjoint experiment with ten factors (J = 10) each with three levels (L j = 3). To evaluate the performance of the proposed method, we consider two different settings; in the first, we assume there are 1,000 respondents who each perform five comparison tasks. In the second, we assume a larger experiment with 2,000 respondents who each perform ten tasks.
In all cases, we assume the true data generating process has three distinct clusters K = 3 and we calibrate the true β k such that the implied average marginal component effects (AMCE) are comparable in magnitude to the empirical effects presented in Section 5. We use a set of six correlated continuous moderators to again mimic a realistic empirical setting and choose {φ k } 3 k=2 to relatively clearly separate respondents into different clusters. Appendix G presents complete description of the simulation settings and the true parameter values used for the β k and marginal effects.
For each sample size, we independently generate 1,000 simulation data sets by drawing N observations of moderators, randomly assigning a true cluster membership to each observation based on the implied probabilities given their moderators, and generating the observed profiles completely at random. We fit our model to the  size; ρ = 0.999 for larger sample size). While the performance overall is reasonably strong, we see that even when the dataset is large there is some degree of attenuation bias due to shrinkage.
The right panel illustrates the frequentist evaluation of our Bayesian posterior standard deviations. We compare the average posterior standard deviation against the standard deviation of the estimated effects across the 1,000 Monte Carlo simulations.
The average posterior standard deviations are noticeably smaller than the standard deviation of the estimates when the sample size is small. For the large sample size, however, there is a much closer relationship. Thus, given sufficiently large data, our approximate Bayesian posterior standard deviations in this simulated example are roughly the same magnitude of the standard deviation of the sampling distribution of the estimator. Even though our method's frequentist coverage is somewhat below the nominal level (see Appendix G.2), this undercoverage appears to be primarily attributable to the shrinkage bias in our regularized estimation procedure rather than the large sample discrepancy between our posterior standard deviations and the corresponding standard deviation of sampling distribution.

Empirical Analysis
In this section, we re-analyze the immigration conjoint data introduced in Section 2, using our methodology. We find evidence of effect heterogeneity for immigrant choice based on respondents traits. In particular, we consistently find that the immigrant's country of origin plays a greater role in forming the immigration preference of respondents with increased prejudice, as measured by a Hispanic prejudice score. However, outside of this cluster, which accounts for about one third or more of the respondents, the country of origin factor plays a much smaller role.

Data and Model
Following the original analysis, our model includes indicator variables for each factor and interactions between country and reason of application factors as well as those between education and job factors in order to account for the restricted randomization.
We additionally include interactions between country and job as well as those between country and education, in accordance with the skill premium theory of Newman and Malhotra (2019). This theory hypothesizes that prejudiced individuals prefer highly skilled immigrants only for certain immigrant countries. This results in a total of 41 AMCEs and 222 AMIEs for each cluster.
For modeling cluster membership, we include the respondents' political party, education, demographics of their ZIP code (we follow the original analysis and include the variables indicating whether respondents' ZIP code had few immigrants (< 5%) and for those from ZIPs with more than 5% foreign-born, whether the majority were from Latin America), ethnicity, and Hispanic prejudice score. The Hispanic prejudice score was used by Newman and Malhotra (2019), though we negate it to make lower values correspond to lower prejudice for easier interpretation. The score is based on a standardized (and negated) feeling thermometer for Hispanics. The score ranges from −1.61 to 2.11 for our sample, where higher scores indicate higher levels of prejudice.
We remove respondents who are Hispanic since the Hispanic prejudice score, which is our moderator, was not measured for these resopndents. After removing entries with missing data, we have a sample of 1,069 respondents. Most respondents evaluated five pairs of profiles, though five respondents have fewer responses in the data set used. The total number of observations is, therefore, 5,337 pairs of profiles. We do not incorporate the survey weights into our analysis to better demonstrate our methods without adding extra variability due to potentially noisy weights.
The original conjoint experiment was conducted using the forced choice design, in which a pair of immigrant profiles are presented and a respondent is asked to choose one of them. To accommodate this design, we follow Egami and Imai (2019) and slightly modify the model specification. In particular, we model the choice as a function of differences in treatments as follows, where T L i and T R i represent the factors for the left and right profiles. The outcome variable Y i is equal to 1 if the left profile is selected and is equal to 0 if the right profile is chosen. With this new linear predictor formulation, the estimation and inference proceed as explained in Section 3.
We conduct the two analyses, one with two clusters and the other with three clusters. These two models perform equally well in terms of out-of-sample classification, a measure that can be used to choose the number of clusters. Using more than three clusters does not appear to give significantly improved substantive insights and provides little improvement in model performance. As noted previously, each analysis optimizes the Bayesian information criterion (BIC) to calibrate the amount of regularization and employs standardization weights to account for factors with different number of levels (see Appendices C.3 and D.3, respectively, for details). We treat education and job experience as ordered factors and only penalize the differences between adjacent levels.

Findings
We focus on the AMCE for each factor as the primary quantity of interest and separately estimate it for each cluster. Under our modeling strategy for the forced choice design, the AMCE of level l versus level l of factor j within cluster k can be written as, The expectation is over individuals and the distribution of the factors not involved in this AMCE. That is, we compute the AMCE separately for the left and right profiles and then average them to obtain the overall AMCE. We estimate this quantity using the fitted model and averaging over the empirical distribution of the factorial treatments. Figure 3 presents the estimated AMCEs and their 95% Bayesian credible intervals for the two-cluster and three-cluster analyses in the left and right panels, respectively.
Cluster 2 in the two-cluster analysis and Cluster 3 in the three-cluster analysis display stronger impacts of country of origin than the other clusters. The respondents in these clusters give the most preference to immigrants from Germany and the least preference to immigrants from Iraq (followed by Sudan). The significant negative effects of Iraq in Cluster 2 of the two-cluster analysis and Clusters 2 and 3 of the three-cluster analysis are consistent with the significant negative effect for Iraq found     Figure 3: Estimated average marginal component effects using a two-cluster (left) and three-cluster (right) analysis. The point estimates and 95% Bayesian credible intervals are shown. A solid circle represents the baseline level of each factor. Number after colon give average posterior predictive probability for each cluster.
by Hainmueller and Hopkins (2015). The patterns we observe for the other factors are also similar for these two clusters in the two analyses.
Across all clusters, respondents prefer educated and experienced immigrants who already have contracts (over those who have no contracts or plans). Respondents also prefer immigrants who have better language skills, although this feature matters less for respondents in Cluster 1 of the three cluster analyses.
For both analyses, the respondents in Cluster 1 do not care much about immigrant's country of origin. Instead, they place a greater emphasis on education and reason for immigration when compared to those in the other clusters. For the threecluster analysis, Clusters 1 and 2 together correspond roughly to Cluster 1 of the two-cluster analysis. In fact, about 81% of the respondents who belong to Clus-ter 1 of the two-cluster analysis are the members of either Cluster 1 or 2 in the three-cluster analysis using a weighted average of their estimated cluster membership posterior probabilities. The differences between Clusters 1 and 2 in the three-cluster analysis are substantively small, but those in Cluster 2 appear to place more emphasis on education and prior entry without legal authorization. Those in Cluster 1, on the other hand, give a slight benefit to immigrants whose reason for immigration is persecution.
It is worth noting that concerns have been raised about comparison of AMCEs across subgroups as they are inherently dependent on the choice of baseline (Leeper, Hobolt and Tilley, 2020). As an alternative, we can compare (and visualize) the marginal means or the values of β, which have sum-to-zero constraints, across clusters.
These alternative avoids issues of baseline dependency in comparisons. A plot of the estimated marginal means is provided in Appendix H. In this case, the results are generally similar to the analysis examined here using AMCEs.
Who belongs to each cluster? The left panel of Figure 4 shows the distribution of Hispanic prejudice score for each cluster weighted by the corresponding posterior cluster membership probability for each individual respondent. The plot shows that for the two-cluster analysis, those with high prejudice score are much more likely to be part of Cluster 2. For the three-cluster analyses, those with high prejudice are more likely to be in Cluster 3. This is consistent with the finding above that the respondents in those clusters put more emphasis on immigrant's country of origin.
The right panel of the figure shows the distribution of other respondent characteristics. In general, Cluster 2 in the two-cluster analysis and Cluster 3 in the three-cluster analysis consist of those who live in ZIP codes with few immigrants and have lower educational achievements. For the three-cluster analysis, those in Cluster 2 tend to be Republican, where-as those in Cluster 1 are more likely Democrats. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2 cluster analysis 3 cluster analysis −1 0 1 2 −1 0 1 2  Figure 4: Distribution of respondent characteristics for each cluster. Left set of plots shows weighted box plots of the Hispanic prejudice moderator within each cluster over the posterior predictive distribution using a two-cluster (left) and three-cluster (right) analysis. Right set of plots shows the distribution of categorical moderators within each cluster over the posterior predictive distribution using a two-cluster (left) and three-cluster (right) analysis.
This is consistent with the finding of a larger penalty for entry without legal authorization in Cluster 2 of the three-cluster analysis. Cluster 3 contains a mix of political ideologies, though with somewhat more respondents who identify as Undecided/Independent/Other or not strong Republican than the other two clusters.
What respondent characteristics are predictive of the cluster membership? We examine how the predicted probabilities of cluster memberships change across respondents with different characteristics. Specifically, we estimate where x 0 and x 1 are different values of covariate of interest X ij . If X ij is a categorical variable, we set x 0 to the baseline level and x 1 to the level indicated on the vertical axis. If X ij is a continuous variable as in the case of the Hispanic prejudice score, then x 0 and x 1 represent the 25th and 75th percentile values. The solid arrows represent  Figure 5: The impact of moderator values on likelihood of being assigned to clusters, for two-cluster (left two plots) and three-cluster (right three plots) analysis. Dark arrows indicate that there is a significant effect of the moderator on cluster membership, i.e. that the corresponding quantity defined in Equation (12) is statistically significant.
whether the corresponding 95% Bayesian credible interval covers zero or not.
Consistent with the earlier findings, Figure 5 shows that those with high Hispanic prejudice scores are predicted to belong to Cluster 2 in the two-cluster analysis and Cluster 3 in the three-cluster analysis even after controlling for other factors. These respondents are also less likely to belong to Cluster 1 in either analysis. For the three-cluster analysis, party ID appears to play a statistically significant role (indicated by dark arrow) with Republicans tending to belong to Cluster 2 or 3 and Democrats tending to belong to Cluster 1. Respondents in Cluster 1 tend to have higher education in both analyses.
Finally, we estimate the average marginal interaction effects (AMIE) between two factors (Egami and Imai, 2019), which can be computed by subtracting the two AMCEs from the average effect of changing the two factors of interest at the same time. Thus, the AMIE represents the additional effect of changing the two factors beyond the sum of the average effects of changing one of the factors alone. Formally, we can define the AMIE of changing factors j and j from levels l j and l j to levels l j and l j , respectively, as follows, All of the AMIE effects found are quite small, so we do not include those results here. According to the skill-premium theory of Newman and Malhotra (2019), we expect to find an interaction between job and country or education in country, in at least some clusters. Unfortunately, our analysis does not find support for this hypothesis.

Concluding Remarks
We have shown that a Bayesian mixture of (regularized) logistic experts can be effectively used to estimate heterogeneous treatment effects of high-dimensional treatments. The proposed approach yields interpretable results, illuminating how different sets of treatments have heterogeneous impacts on distinct groups of units. We apply our methodology to conjoint analysis, which is a popular survey experiment in marketing research and social sciences. Our analysis shows that individuals with high prejudice score tend to discriminate against immigrants from certain non-European countries. These individuals tend to be less educated and live in areas with few immigrants. Future research should consider the derivation of optimal treatment rules in this setting as well as the empirical evaluation of such rules. Another important research agenda is the estimation of heterogeneous effects of high-dimensional treatments in observational studies. Has a contract with a U.S. employer; Does not have a contract with a U.S. employer, but has done job interviews; Will look for work after arriving in the U.S.; Has no plans to look for work at this time Prior Trips to the U.S. 5 Never been to the U.S.; Entered the U.S. once before on a tourist visa; Entered the U.S. once before without legal authorization; Has visited the U.S. many times before on tourist visas; Spent six months with family members in the U.S.

B Propriety of the Structured Sparse Prior
The proof of propriety for the structured sparse prior used in our paper is an application of Theorem 1 established in Goplerud (2021) and is reproduced here.
Theorem B.1 (Goplerud (2021)). Consider the following structured sparse prior on β ∈ R p with regularization strength λ > 0 penalizes K linear constraints d k and L quadratic constraints F on the parameters where F is symmetric and positive semi-definite. The kernel of the prior is shown below.
Then, for λ > 0, the prior above is proper if and only ifD is full column rank.
In our specific case, we note that K = 0, L = G, and λ = λπ γ k . Prior propriety of p(β k | {φ k } K k=2 , λ), therefore, can be determined by empirically investigating whetherD, i.e. the vertically stacked F , is full column rank.
It is also possible to analytically show the propriety of the prior distribution in all cases considered in this paper. We focus on the case of K = 1 and arbitrary λ > 0 as the result follows automatically for the case in our paper.
Result B.1. Assume a structured sparse prior for a factorial or conjoint design with J factors each with L j levels where all pairwise interactions are included and levels of each factor are encouraged to be fused together (i.e. the model in the main text). The kernel of the prior is shown below where F g are as defined in the main text.
Assume that the linear sum-to-zero constraints C β = 0 hold. Then, the structured sparse prior on the unconstrainedβ such thatβ ∈ N (C ) is proper. Or, equivalently, the following result holds.
Proof. Let B C represent a basis for the linear constraints C . The integral for evaluating propriety can be written as, Note that F g can be expressed as a sum of N g outer products of |β|-length vectors of the form l i ∈ {−1, 0, 1} where −1 and 1 correspond to the two terms that are fused together and all other elements are 0, i.e., F g = Ng g =1 l g l g . Thus, one can define a matrix Q g = l 1 , · · · , l Ng such that Q g Q g = F g , which allows us to rewritek(β) as: By applying Theorem B.1 and noting that the nullspaces of A T A and A are identical, the integral ofk(β) is finite if and only if QB C is full column rank, where Q = [Q 1 , · · · , Q G ]. We demonstrate this fact in two steps. First, there exists a permutation matrix P Q such that P Q Q has a block diagonal structure with J + 1 diagonal blocks. The first J blocks corresponding to the main terms for each factor j and the last block corresponds to all interaction terms. The nullspace of each block is spanned by the vector 1 as the corresponding block of P Q Q is a (transposed) orientated incidence matrix of a fully connected graph. Thus, the nullspace of P Q Q, and hence Q, is spanned by the J + 1 columns of a block diagonal matrix with 1 on each block. Second, consider the linear constraints C β = 0. The only vector to satisfy this constraint and lie in the nullspace of Q must be 0 as, for each block, the only vector proportional to 1 and satisfying the corresponding sum-to-zero constraints must be 0. Thus, QB C is full column rank and the prior is proper.

C Derivations for the Basic Model
This section derives a number of results for the basic model. It first restates the main results concerning the elimination of the linear constraints C β k = 0. Then, it derives the Expectation Maximization algorithm, our measure of degrees of freedom, and some additional computational improvements used to accelerate estimation. In the following, with a slight abuse of notation, we use t i to denote the corresponding vector of indicators for whether certain treatments or interactions are present (i.e. stacking all I(T ij = 1), etc. from Equation A14). In addition, we use ψ ik to indicate the linear predictor for observation i and cluster k.

C.1 Removing the Linear Constraints
The inference problem in the main text is presented as an optimization problem subject to linear constraints on the coefficients β k . We note that inference is noticeably easier if these are eliminated via a transformation of the problem to a lower-dimensional one by noting that β k must lie in the null space of the constraint matrix C .
where B C is a basis for the null space of C . Given this, the problem can be solved in terms of the unconstrainedβ k ∈ R p−rank(C ) . The treatment design vectors (t i = B C t i ) and penalty matrices (F g as B C F g B C ) are adjusted and then the augmented problem can be written in Equations (9a) and (9b). The new linear predictor is ψ i,k = t i β Z i + µ.
Note that the only difference between the unconstrained and constrained problem is the choice of treatment design matrix. Given this similarity and to reduce the burden of notation, we thus present all results herein dropping the "tilde" notation and note that, once estimated,β k is projected back into the original space for the reported coefficients, average marginal component effects, etc. Given Appendix F's results on approximatingβ k as normal, β k will have a (singular) normal distribution distribution.

C.2 Expectation Maximization Algorithm
Assume that we have the model after the linear constraints have been removed as discussed above. We begin with the cycle of the AECM algorithm for updating {β k } K k=1 given {φ K k=2 }. The model with triple augmentation is shown below, mirroring Equation (9) from the main text.
We treat {ω i , Z i } N i=1 and {τ 2 gk } as missing data to be used in the E-Step and β k as parameters to be optimized in an M -Step. For this cycle of the AECM algorithm, we hold φ constant. The E-Step can be computed by noting that p(ω i , Z i | Y i , θ) and p(τ 2 gk | θ) are conditionally independent across i and gk. The distributions are shown below, The relevant expectations are shown below, Note that as β k F g β k approaches zero, E(τ −2 gk ) approaches infinity. To prevent numerical instability, we rely on the strategy in Goplerud (2021) (inspired by Polson and Scott 2011) where once it is sufficiently small, e.g. below 10 −4 , and thus the restriction is almost binding, we ensure that restriction holds in all future iterations. We do so by adding a quadratic constraint β k F g β k = 0. Note that this again implies that β k lies in the nullspace of F g and thus with an additional transformation, it can be removed and the problem be solved in an unconstrained space with a modified design.
With these results, the Q β -function can be thus expressed as the expectation of the logarithm of the augmented posterior where we suppress all terms that do not depend on β k and all expectations depend on prior value θ (t) . Recall that β = {{β k } K k=1 , µ}.
From the above, define a vector that concatenatesβ = [µ, β 1 , · · · , β K ] . We can create a corresponding designŤ = [1 N , With this in hand, we can rewrite the Q β function as proportional to the following ridge regression problem and the corresponding update in the M -Step.
Rβ + const., Note that sinceΩ,Y , and R depend on θ (t) , we can think about this algorithm as follows where we denote their dependence on prior values with the superscript (t).
In practice, we rely on a generalized EM algorithm where Q β is improved versus maximized for computational reasons; we do so using a conjugate gradient solver initialized ať β (t) . The update for {φ k } K k=2 is shown in the main text; as is required for AECM algorithms, it requires a re-computation of the E-Step for E(z ik ) that depends on β (t+1) . We use L-BFGS-B for a few steps to improve the objective.

C.3 Degrees of Freedom
We note that our procedure for estimatingβ (t) appears similar to the results in Oelker and Tutz (2017) where complex regularization and non-linear models can be recast as a (weighted) ridge regression. Using that logic, we take the trace of the "hat matrix" implied by our algorithm at stationarity to estimate our degrees of freedom. We also adjust upwards the degrees of freedom by the number of moderator coefficients (e.g., Khalili 2010;Chamroukhi and Huynh 2019). Equation (A5) shows our procedure where R andΩ contain expectations calculated at convergence. p x denotes the number of moderators, i.e. the dimensionality of φ k . Before evaluating Equation (A5), for any two factor levels that are sufficiently close (e.g., β k F g β k < 10 −4 ), we assume they are fused together and consider it as an additional linear constraint on the parameter vector β k .
From this, we can calculate a BIC criterion. We seek to find the regularization parameter λ that minimizes this criterion. To avoid the problems of a naive grid-search, we use Bayesian model-based optimization that attempts to minimize the number of function evaluations while searching for the value of λ that minimizes the BIC (mlrMBO; Bischl et al. 2018). We find that with around fifteen model evaluations, the optimizer can usually find a near optimal value of λ.

C.4 Computational Improvements
While the algorithm above provides a valid way to locate a posterior mode, our estimation problem is complex and high-dimensional. Furthermore, given the complex posterior implied by mixture of experts models, we derived a number of computational strategies to improve convergence. We use the SQUAREM algorithm (Varadhan and Roland 2008) and a generalized EM algorithm to update β using a conjugate gradient approach and φ using a few steps of L-BFGS-B.
We also outline a way to deterministically initialize the model to provide stability and, again, speed up estimation on large problems. To do this, we adapt the procedure from Murphy and Murphy (2020) for initializing mixture of experts: (i) initialize the clusters using some (deterministic) procedure (e.g. spectral clustering on the moderators); (ii) using only the main effects, estimate an EM algorithm-possibly with hard assignment at the E-Step (CEM; Celeux and Govaert 1992); (iii) iterate until the memberships have stabilized. Use those memberships to initialize the model. This has the benefit of having a deterministic initialization procedure where the cluster members are based on the moderators but guided by which grouping seem to have sensible treatment effects, at least for the main effects. Given the memberships, update β using a ridge regression and φ using a ridge regression and take those values as β (0) and φ (0) .

D Extensions to the Basic Model
As noted in the main text, there are five major extensions to the basic model that applied users might wish to include: 1. Repeated tasks (observations) for a single individual 2. A forced-choice conjoint experiment 3. Survey to weight the sample estimate to the broader population 4. Adaptive weights for each penalty 5. Latent overlapping groups All can be easily incorporated into the proposed framework above. This section outlines the changes to the underlying model.

D.1 Repeated Observations
This modification notes that in factorial and conjoint experiments it is common for individuals to perform multiple tasks. Typically, the number of tasks N i is similar across individuals. The updated likelihood for a single observation i is shown below; we show both the observed and complete case. y it represents the choice of person i on task t ∈ {1, · · · , N i }; p itk is the probability of Y it = 1 if person i was in cluster k.
Note that because of the conditional independence of (y it , ω it ) given Z i and the parameters, the major modifications to the EM algorithm is that the E-Step must account for all t observations, i.e. the terms summed in Equation (A6). Some additional book-keeping is required in the code as the design of the treatments has N i=1 N i rows whereas the design of the moderators has N rows. Repeated observations can be easily integrated into the uncertainty estimation procedure outlined below.

D.2 Forced Choice Conjoint Design
A popular design of a conjoint experiment is the forced choice design where the respondents are required to choose between two profiles. Therefore, the researcher does not observe an outcome for each profile separately, but rather a single outcome is observed for each pair indicating which is preferred. Egami and Imai (2019) show that this can be easily fit into the above framework with some adjustment. Specifically, the model is modified to difference the indicators of the treatment levels for the pair of profiles (subtracting, e.g., the levels of the profile presented on the left from those of the profile presented on the right). The intercept for this model can be interpreted as a preference for picking a profile presented in a particular location. With this modification, estimation proceeds as before.

D.3 Standardization Weights
An additional modification to the problem is to weight the penalty. This could be done for two reasons. First, there is an issue of the columns having different variances/Euclidean norms because of the different number of factor levels L j . Second, it is popular to weight the penalty based on some consistent estimator (e.g. ridge regression) to improve performance and, in simpler models, can be shown to imply various oracle properties (e.g. Zou 2006). We leave the latter to future exploration.
Define ξ gk as a positive weight for the g-th penalty and the k-th cluster. The kernel of the penalty is modified to include them.
This has no implication on the rank of the stacked F g (and thus the results in Appendix B) as they are all positive and thus only slightly modify the E-Step.
We employ weights in all of our analyses to account for the fact that different factors j may have different number of levels L j . We use a generalization of the weights in Bondell and Reich (2009) to the case of penalized differences. Specifically, consider the over-parameterized model in Appendix D.4 where the penalty can be written entirely on the differences δ Main , δ Int , δ Main−Copy . Note that each of those penalties has a simple (group) LASSO form and thus we adopt the approach in Lim and Hastie (2015) of weighting by the Frobenius norm of the associated columns in T LOG , i.e. the over-parameterized design matrix. At slight abuse of notation, define [T LOG ] g as the columns of T LOG corresponding to the differences penalized in the (group) lasso g, the weight can be expressed as follows: Ignoring the factor of √ N , this exactly recovers the weight proposed in Bondell and Reich (2009) in the non-latent-overlapping non-interactive model of (L j + 1) −1 N j l + N j l where N j l , N j l are the number of observations for factor j in level l and level l that are being encouraged to fuse together by the penalty in group g.

D.4 Latent Overlapping Groups
One feature of the above approach is that our groups are highly overlapping. Yan and Bien (2017) suggest that, in this setting, a different formulation of the problem may result in superior performance (see also Lim and Hastie 2015). Existing work on the topic has focused on group LASSO penalties (e.g. F g = I) and thus some modifications are needed for our purposes. To address this, we note that we can again recast our model in an equivalent fashion. Instead of penalizing β k F g β k , we can penalize the vector of differences between levels as long as we also impose linear constraints to ensure that the original model is maintained.
Consider a simple example with two factors each with two levels {1, 2} and {A, B}. The relevant differences are defined such that δ j 1−2 = β j 1 − β j 2 and δ jj (lm)−(l m ) = β j l,m − β j l ,m . The equivalent penalty can be imposed as follows: The latent overlapping group suggests a slight modification. In addition to the above penalization of the 2 norm of the main and interactive differences, 2 it duplicates the main effect and penalizes it separately while ensuring that all effects maintain the accounting identities between the "latent" groups and the overall effect. Specifically, it modifies the above penalty to duplicate the column corresponding to δ j 1−2 and adds a new parameter δ j (1−2)−Copy .
Scoping out to the full problem, define δ Main as the main effect differences, e.g. δ j 1−2 , and δ Int as the interaction differences and D Main as the matrix such that D Main β = δ Main , and D Int as the corresponding matrix to create the vector of interactions. Define δ Main−g as the sub-vector of δ Main−g that corresponds to the (main) effect differences between levels l and l of factor j penalized by F g in the original notation. Similarly define δ Int−g and δ Main−Copy−g .
This also requires a modification of the design matrix T to ensure that (i) its dimensionality conforms with the expanded parameter vector and (ii) that for any choice of the expanded parameter that satisfies the constraints, the linear predictor for all observation (and thus the likelihood) is unchanged. Consider first the simple case without latent-overlapping groups. In this case, following Bondell and Reich (2009), note that the expanded design can be expressed asT = TM † whereM = [I, D Main , D Int ] andM † is a left-inverse of M . The latent-overlapping group formulation is a simple extension; we copy the columns ofT that correspond to δ Main and append them to get T LOG .
With this new design and parameterization in hand, we can again use the above results on projecting out the linear constraints to turn the problem into inference on an unconstrained vector β k with a set of positive semi-definite constraints {F g } 2G g=1 and inference proceeds identically to before.

E Estimators
Here we provide further details on the estimators. In particular, we discuss estimation of Average Marginal Component Effects (AMCEs) and Average Marginal Interaction Effects (AMIEs) based on our model. We consider a traditional factorial design, where each unit receives one treatment (profile), and a conjoint design in which each unit compares two treatments (profiles). We also discuss the impact of randomization restrictions on estimators and implied changes in interpretation of estimands.

E.1.1 Without restrictions on randomization
For a unit in cluster k we have where i = 1, 2, . . . , N and for k = 1, 2, . . . , K, . (A13) We model ψ k (T i ) as for each k = 1, 2, . . . , K, with constraints where β k is a stacked column vector containing all coefficients for cluster k. We can rewrite this to aid in the interpretation of β k as follows: Thus, β j kl − β j kf is the AMCE going from level f to level l of factor j on the logit probability of Y i = 1 scale.
Let t be some combination of the J factors, where t j is the jth factor's level and t −j is the levels for all factors except j. This allows us to easily write, taking expectation over units in cluster k, where T ij = l indicates for unit i forcing factor j to be assigned level l and T i,−j = t −j indicates forcing the assignment on all factors except for j to be assigned levels as in t −j . The causal effects of interest (on the original Y scale) are defined as contrasts of these expectations. Without additional weighting (i.e., using traditional uniform weights for marginalization), the AMCE for level l vs f of factor j in cluster k is, where M is the number of possible combinations of the other J − 1 factors (e.g., if we had J 2-level factors, M = 2 J−1 ). We can estimate this by plugging in the coefficients directly. Note that, because of the nonlinear nature of the estimator, this approach is consistent (under model assumptions) but not unbiased. Alternatively, instead of summing over all possible t −j , we can use the empirical distribution of t −j in the sample. This potentially changes the estimand. Define estimators Then we can use the following overall estimator for the AMCE: This is a consistent estimator (under model assumptions) of conditioning on the treatments we actually observed. Now, we turn to examination of the AMIEs. Without additional weighting (i.e., using traditional uniform weights for marginalization), the AMIE for level l of factor j and level q of factor s vs f of factor j and level r of factor s in cluster k is where ACE * (l, f, q, r) where M * is the number of possible combinations of the other J − 2 factors (e.g., if we had J two-level factors, M * = 2 J−2 ). We can use the following overall estimator for the ACE: This is then combined with the estimators for the AMCEs to get

E.1.2 With restrictions on randomization
In this section we consider restricted randomization conditions. Let us assume that factor j and factor h are such that some levels of j are not well defined and hence excluded in combination with some levels of factor h under the randomization set up. Let S(j, h) ⊂ {1, . . . , L j } be the set of levels of factor j that are not defined for some levels of factor h. Similarly, let S(h, j) ⊂ {1, . . . , L h } be the set of levels of factor h that are not defined for some levels of factor j. In our example, if j is education and h is profession, we have S(j, h) = {No formal, 4th grade, 8th grade, High school} and S(h, j) = {Financial analyst, Research scientist, Doctor, Computer programmer}. When estimating the AMCE for level l vs f of factor J − 1 in cluster k, using the model rather than the empirical distribution, we consider, where M def (j,h) is the number of possible combinations of the other factors, restricted such that t h / ∈ S(h, j) (e.g., if we had J 3-level factors, and some of the levels of factor j were not defined for one level of factor h, this would be 2 × 3 J−2 ).
To use empirical distribution, we need a way to deal with profiles that are not well defined. We can accomplish this by only aggregating over those profiles that are sensible for all levels of factor j. That is, we use the following estimator, Consider the case where we are estimating the AMCE for "doctor" vs "gardener" for profession. Because of the randomization restriction between certain professions and level of education, we will remove any profiles that have "4th grade" as level of education. Although "gardener" with "4th grade" education is allowable under the randomization, we must remove such profiles to have an "apples-to-apples" comparison with profession of doctor, which is not allowed to have "4th grade" education. Note that we do this dropping of profiles even if we are comparing "waiter" vs "gardener" for profession, which are both allowed to have "4th grade" as level of education, to ensure that all AMCEs for profession comparable.
Similarly for the AMIEs, we restrict the profiles we marginalize over to be only those that are defined for both factors in the interactions. Let factor j be restricted by some other factor h and let factor s be restricted by some other factor w. Then we have the following estimator, ACE * (l, f, q, r) The relevant AMCEs should be similarly restricted within the AMIE estimator, with restrictions applied based on the restrictions for all levels both factors in the interaction.

E.2.1 Without restrictions on randomization
Consider a conjoint experiment in which each unit i only compares two profiles. The response Y i indicates a choice between two profiles. Let T L i be the levels for the left profile and T R i be the levels for the right profile that unit i sees. Here, we modify how we model ψ k to If we use Y i = 1 to indicate that unit i picks the left profile, then we have, .
We can use the symmetry assumption that choice order does not affect the appeal of individual attributes. That is, there may be some overall preference for left or right accounted for by µ, but this preference is not affected by profile attributes. Then, we can define our effects, on the original Y scale, as contrasts of these expectations. Without additional weighting, the AMCE for level l vs l of factor j in cluster k is, To save space, the outer expectation is over the random assignment, which corresponds to the expectation over theM possible combinations of the two profiles on the other J − 1 factors (e.g., if we had J two-level factors, this would be 4 J−1 ). We can again estimate this by plugging in our coefficient estimates directly. Alternatively, instead of summing over all possible t L −j and t R −j , we can use the empirical distribution of t L −j and t R −j in the sample. Define ) .
Then we can use the estimator Now we turn to examination of the AMIEs. Without additional weighting (i.e., using traditional uniform weights for marginalization), the AMIE for level l of factor j and level q of factor s vs m of factor j and level r of factor s in cluster k is Here we can use the estimator This gives us AMIE jsk (l, f, q, r) = ACE(l, f, q, r) − δ jk (l, f ) − δ sk (q, r).

E.2.2 With restrictions on randomization
Similar to Appendix E.1.2, adjustments to estimation need to be made when we have restricted randomizations. We again will do this by dropping profiles that have levels of factors not allowable for all levels of the factor(s) whose effects we are estimating (e.g., profiles with "4th grade" for education when estimating an effect for profession). However, now we estimate the effect for the right profile and the effect for the left profile, and then average the two (they should be equal under symmetry). When estimating the effect for the right profile, therefore, we will only drop pairings if the right profile has a level that is not allowed for some level of the factor we are estimating an effect of. For example, dropping pairings where the right profile has "4th grade" as level of education when estimating main effects of profession because "doctor" cannot have level "4th grade." Again, this will drop more profiles than those that are not allowed under randomization to ensure an "apples-toapples" comparison across levels of profession.
In this calculation, we use the empirical distribution for the levels of the left profile (which represents the "opponent"). Thus, the distribution of other factors for the profile we are calculating the effect of may differ than that distribution for its opponents. Similarly, when estimating the effect for the left profile, we only drop pairings in which the left profile has a restricted level for some level of the factor of interest. Estimation for the AMIE under randomization restrictions follows similarly.

F Quantification of Uncertainty
We quantify uncertainty in our parameter estimates by inverting the negative Hessian of the log-posterior at the estimatesθ, i.e. − ∂ ∂θθ T ln p(θ|Y i ) θ=θ or I(θ). This can be stably and easily computed using terms from the AECM algorithm following Louis (1982)'s method. Specifically, consider the model from the main text augmented with Z i , i.e. the cluster memberships. Recall that z ik = 1(Z i = k) for notational simplicity.
(A16) Louis (1982) notes that equation can be used to compute I L (θ), where the subscript L denotes its computation via this method.
To address the issue with the non-differentiablity of the penalty on β (and thus L c (θ)), we follow the existing research in two ways. First, for restrictions that are sufficiently close to binding, we assume them to bind and estimate the uncertainty given those restrictions.
That is, we identify the binding restrictions such that β k F g β k is sufficiently small (say 10 −4 ) and note that if these are binding, we can use the nullspace projection technique to transform β k such that it lies in an unconstrained space.
To further ensure stability, we modify the penalty with a small positive constant ≈ 10 −4 to ensure that the entire objective is (twice) differentiable. For notational simplicity, we derive the results below assuming β k represent the parameter vector after projecting into a space with no linear constraints. The approximated log-posterior is shown below and denoted with a tilde. We thus evaluate I L (θ) usingL c in place of L c .
This procedure has some pleasing properties that mirror existing results on approximate standard errors after sparse estimation; consider a simple three-level case: β j 1 , β j 2 , β j 3 . If β j 1 and β j 2 are fused, then their approximate point estimates and standard errors will be identical but crucially not zero. This is because while their difference is zero and assumed to bind with no uncertainty, this does not imply that the effects, themselves, have no uncertainty: β j 1 − β j 2 will have a standard error of zero in our method. This thus mirrors the results from Fan and Li (2001) where effects that are shrunken to zero by the LASSO are not estimated with any uncertainty. One might relax this with fully Bayesian approaches in future research.
Second, note that if all levels are fused together, i.e. β j 1 = β j 2 = β j 3 , then all point estimates must be zero by the ANOVA sum-to-zero constraint and all will have an uncertainty of zero. Thus, when an entire factor is removed from the model, the approximate standard errors return a result consist with existing research.

F.1 Derivation of Hessian
To calculate the above terms, the score and gradient ofL c are required. They are reported below:S The above results use the following intermediate derivations: Second, the variance ofS c (θ) over p({z ik } | θ). This is derived blockwise below.
This provides all terms needed to compute I L (θ).

F.2 Repeated Observations
Now consider the case of repeated observations per individual i. In this scenario, each individual i performs N i tasks. Note, after augmentation, the score has exactly the same form and thus the complete ScoreS c and HessianH c are identical where the sum merely now runs over N i=1 N i t=1 . The average forπ k is similarly a weighted average by N i , although note that often each respondent answers an identical number of tasks so it is, effectively, the same as before. The covariance ofS c is adjusted as shown below.

F.3 Standard Errors on Other Quantities of Interest
Given the above results, we derive an approximate covariance matrix onθ. We calculate uncertainty on other quantities of interest, e.g. AMCE and marginal effects, using the multivariate delta method. As almost all of our quantities of interest can be expressed as (weighted) sums or averages over individuals i ∈ {1, · · · , N }, calculating the requisite gradient for the multivariate delta method simply requires calculating the relevant derivative for each observation. For example, all derivatives needed in the AMCE are of the following form; see Appendix E for more details.

G Simulations
We detail our simulations and provide additional results in this section.

G.1 Setup
We generate one set of true coefficients β k that are calibrated to generate average marginal component effects that are broadly comparable to the results we found in our empirical example, i.e. ranging between around -0.30 and 0.30. Our β k and {φ k } 3 k=2 are one draw from the following distributions: Simulating β k : 1. For each factor j and cluster k, draw the number of unique levels u with equal probability from {1, 2, 3}.
2. Draw u normal random variables independently from N (0, 1/3); call these b j ku . 3. For u = 1, set β j kl = 0 4. For u = 3, de-mean {b j ku } 3 u=1 drawn in (2) and set all β j kl equal to the corresponding value. 5. For u = 2, assign b j k3 equal to one of the two b j ku with equal probability. De-mean the {b j ku } 3 u=1 and set β j kl equal to the corresponding values.
Simulating φ k : {φ k } K k=1 ∼ N (0, 2 · I) We calculate the true average marginal component effects using a Monte Carlo based approach where we sample 1,000,000 pairs of treatment profiles for the other attributes to simulate marginalizing over the other factors analytically. The true distribution of our β k and average marginal component effects (with a baseline level of '1') are shown below.
For each simulation, we draw N individuals who rate T profiles where (N , T ) ∈ {(1000, 5), (2000, 10)}. For each individual i , we draw its moderators x i from a correlated 17 Figure   We see that the members are well-separated; the clusters are somewhat unbalanced, i.e. π = [0.217, 0.261, 0.522]. If we consider the maximum probability for each person i, i.e. π * i = max k∈{1,2,3} π ik , this distribution has a median of 0.93, a 25th percentile of 0.75 and a 75th percentile of 0.99.
In terms of simulating the treatment profiles and outcome, for each individual i , we draw their true cluster membership z i using their π i . For each task t , we then draw a pair of treatments at random and then, given z i , draw the outcome y i using the generative model outlined in the main text. After estimating our model with K = 3, we resolve the problem of label switching by permuting our estimate cluster labels to minimize the absolute error between the estimated posterior membership probabilities {E[z ik |θ]} K k=1 and z i (the one-hot assignment of cluster membership).

G.2 Additional Results
We first provide the results on the β k to complement the results on the average marginal component effects shown in the main text. The same broad patterns exist: There is a close correlation with some shrinkage of the estimates for both large and small sample sizes. For the large sample size, the average standard errors are broadly comparable to the standard deviation of the estimated coefficients.  Tasks) 2,000 People (10 Tasks) Figure A3: The performance of the estimator on simulated data. The black squares indicate the effects estimated with the smaller sample size (1,000 people completing 5 tasks); the red crosses indicate effects estimated with the smaller sample size (2,000 people completing 10 tasks).
We finally report the coverage. We do this in two ways; first, for each of the reported AMCEs or elements of β k , we calculate the coverage across the 1,000 simulations. Given that our method gives a standard error of zero if the AMCE or β k are all fused to zero, we also show post-selection coverage (i.e., the coverage conditional on a non-zero estimate). Figure A4 shows a boxplot of the coverage of the AMCEs and β k , respectively, at both sample sizes. The solid horizontal line indicates nominal frequentist coverage at 0.95; the dashed line is nominal frequentist coverage at 0.90. We see that coverage is somewhat below nominal in the case of the AME, although post-selection coverage is close to 0.90 on median. The coverage for the β k is noticeably more sluggish. Our results above, however, suggest that this is due to non-vanishing bias in the coefficient estimates rather than our standard errors being systematically too small-at least in the case of large sample sizes.

H Additional Results for Immigration Conjoint Experiment
We provide some additional results for the three-cluster model in Figure A5. We do this by an analogue to the "marginal means" estimator in Leeper, Hobolt and Tilley (2020). We compute the probability of a profile being chosen without specifying a baseline category. The equation is shown below for the forced choice case; note it consists of two of the terms used for the AMCE.
The below plot ignores randomization restrictions when estimating this quantity to center the estimate around 0.50 as in Leeper, Hobolt and Tilley (2020). The results are substantively similar to the analysis in shown in the main paper using AMCEs.