Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
98,769 characters · 25 sections · 43 citation commands
Synthetic Combinations: A Causal Inference Framework for Combinatorial Interventions
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Matrix Completion, Fourier Analysis of Boolean Functions, Latent Factor Models, Synthetic Controls, Learning to Rank
\spacingset{1.75}
Modern-day decision-makers, in settings from e-commerce to public policy to medicine, often encounter settings where they have to pick a combination of interventions and ideally would like to do so in a highly personalized manner. Examples include recommending a curated basket of items to customers on a commerce platform, deciding on a combination of therapies for a medical patient, enacting a collection of socio-economic policies for a specific geographic location, feature selection in a machine learning model, doing a conjoint analysis in surveys, etc. Despite the ubiquity of this setting, it comes with significant empirical challenges: with $p$ interventions and $N$ units, a decision maker must evaluate $N \times 2^p$ potential combinations in order to confirm the optimal personalized policy. With large $N$ and even with relatively small $p$ (due to its exponential dependence), it becomes infeasible to run that many experiments. Of course, in observational data there is the additional challenge of potential unobserved confounding. Current methods tackle this problem by following one of two approaches: (i) they impose structure on how combinations of interventions interact, or (ii) they assume latent similarity in potential outcomes across units. However, as we discuss in detail below, these approaches require a large number of observations to estimate all $N \times 2^p$ potential outcomes because they do not exploit structure across both units and combinations. Hence the question naturally arises: {\em how can one effectively share information across both units and combinations of interventions?}
{\bf Contributions.} Our contributions may be summarized as follows.
(1) For a given unit $n \in [N]$, we represent its potential outcomes over the $2^p$ combinations as a Boolean function from $\{-1, 1\}^p$ to $\mathbb{R}$, expressed in the Fourier basis. To impose structure across combinations, we assume that for a unit $n$, the Fourier coefficients $\boldsymbol{\alpha}_n \in \mathbb{R}^{2^p}$ induced by this basis are sparse, i.e., have at most $s$ non-zero entries. To impose structure across units, we assume that this matrix of Fourier coefficients across units $\mathcal{A} = [\boldsymbol{\alpha}_n]_{n \in [N]} \in \mathbb{R}^{N \times 2^p}$ has rank $r$. This simultaneous sparsity and low-rank assumption is indeed what allows one to share information across both units and combinations.
(2) We establish identification for the $N \times 2^p$ potential outcomes of interest, which requires that any confounding is mediated by the (unobserved) matrix of Fourier coefficients $\mathcal{A}$.
(3) We design a two-step algorithm “Synthetic Combinations” and prove it consistently estimates the various causal parameters, despite potential unobserved confounding. The first step of Synthetic Combinations, termed “horizontal regression”, learns the structure across combinations of interventions via the Lasso. The second step, termed “vertical regression”, learns the structure across units via principal component regression (PCR).
(4) Our results imply that Synthetic Combinations is able to consistently estimate unit-specific potential outcomes given a total of $\text{poly}(r) \times \left( N + s^2p\right)$ observations (ignoring logarithmic factors). This improves over previous methods that do not exploit structure across both units and combinations, which have sample complexity scaling as $\min(N \times s^2p, \ \ \text{poly}(r) \times (N + 2^p))$. A summary of the sample complexities required for different methods can be found in Table (ref). We provide an example showing how replacing Lasso with “classification and regression trees” (CART) leads to a sample complexity of $\text{poly}(r) \times \left( N + s^2\right)$ under additional regularity assumptions on the Fourier coefficients. A key technical challenge in our theoretical analysis is studying how the error induced in the first step of Synthetic Combinations percolates through to the second step. To tackle it, we reduce this problem to that of high-dimensional error-in-variables regression with linear model misspecification, and do a novel analysis of this statistical setting. Under additional conditions, we establish asymptotic normality of the estimator.
(5) We show how Synthetic Combinations can be used to inform experiment design for combinatorial interventions. In particular, using ideas from compressed sensing, we propose and analyze an experimental design mechanism that ensures the key assumptions required for consistency of Synthetic Combinations are satisfied.
(6) To empirically evaluate Synthetic Combinations, we apply it to a real-world dataset for user ratings on sets of movies sharma2019sets, and find it outperforms both Lasso and matrix completion methods. Further, we show two of our key modeling assumptions---low-rank and sparse Fourier coefficients---hold in this dataset. We perform additional numerical simulations that corroborate our theoretical findings and show the robustness of Synthetic Combinations to unobserved confounding.
(7) We discuss how to extend Synthetic Combinations to estimate counterfactual outcomes when the intervention is a permutation over $p$ items, i.e., rankings. Learning to rank is a widely studied problem and has many applications such as search engines and matching markets.
{\bf Learning over combinations.} To place structure on the space of combinations, we use tools from the theory of learning (sparse) Boolean functions, in particular, the Fourier transform. Sparsity of the Fourier transform as a complexity measure was proposed by karpovsky1976finite, and was used by brandman1990spectral and others to characterize and design learning algorithms for low-depth trees, low-degree polynomials, and small circuits. Learning Boolean functions is now a central topic in learning theory, and is closely related to many important questions in ML more broadly; see e.g. mossel2003learning for discussion of the $k$-Junta problem and its relation to relevant feature selection. We refer to o2008some for further background on this area. In this paper, our focus is on general-purpose statistical procedures for learning sparse Boolean functions with noise. We build on the work of negahban2012learning on variants of the Lasso procedure, and the works of breiman2017classification,syrgkanis2020estimation,klusowski2020sparse,klusowski2021universal and others on CART. In particular, we highlight the works of negahban2012learning and syrgkanis2020estimation which respectively show that Lasso and CART can be used to efficiently learn sparse Boolean functions. Recent work bravohermsdorff2023intervention has also explored the use of graphical models to learn across combinations.
{\bf Matrix completion.} Previous works have shown that imputing counterfactual outcomes with latent factor structure and potential unobserved confounding can be equivalently expressed as low-rank matrix completion with missing not at random data bai2019matrix, athey2021matrix, agarwal2020synthetic, agarwal2021causalmatrix. The observation that low-rank matrices may typically be recovered from a small fraction of the entries by nuclear-norm minimization has had a major impact on modern statistics candes2010power,recht2011simpler,candes2012exact. In the noisy setting, proposed estimators have generally proceeded by minimizing risk subject to a nuclear-norm penalty, such as in the SoftImpute algorithm of mazumder2010spectral, or minimizing risk subject to a rank constraint as in the hard singular-value thresholding (HSVT) algorithms analyzed by keshavan2009matrix,keshavan2010matrix,gavish2014optimal, chatterjee2015matrix. We refer the reader to ieee_matrix_completion_overview, nguyen2019low for a comprehensive overview of this vast literature.
{\bf Causal latent factor models.} There is a rich literature on how to learn personalized treatment effects for heterogeneous units. This problem is of particular importance in the social sciences and in medicine, where experimental data is limited, and has led to several promising approaches including instrumental variables, difference-in-differences, regression discontinuity, and others; see angrist1996identification,imbens2015causal for an overview. Of particular interest is the synthetic control method abadie2010synthetic, which exploits an underlying factor structure to effectively impute outcomes under control for treated units. Building on this, recent works agarwal2020synthetic, agarwal2021causalmatrix have shown how to use latent factor representations to efficiently share information across heterogeneous units and treatments. This paper generalizes these previous works to where a treatment is a combination of interventions and {\em most treatments have no units that receive it}. To do so, we impose latent structure across combinations of interventions.
In doing so, we aim to bridge such latent factor models with highly practical settings where multiple interventions are delivered simultaneously, such as slate recommendations, medicine (e.g. basket trials, combination therapies), conjoint analysis in surveys egami2019causal,haimmueller2014causal,goplerud2023estimating, and factorial design experiments (e.g., multivariate tests in digital experimentation) duflo2007using,dasgupta2015causal,wu2011experiments,zhao2022regression . For example, structure across combinations is often assumed in the analysis of factorial design experiments, where potential outcomes are generally assumed to only depend on main effects and pairwise interactions between interventions bertrand2004emily,delaCuesta2021ImprovingTE,eriksson2014employers,george2005statistics. We discuss factorial design experiments in detail in Section (ref).
In this section, we first describe requisite notation, background on the Fourier expansion of real-valued functions over Booleans, and how it relates to potential outcomes over combinations. Next, we introduce the key modeling assumptions and data generating process (DGP), along with the target causal parameter of interest. We also provide a brief summary of how the formalism can be extended to potential outcomes over permutations in Section (ref), and a more detailed discussion in Appendix (ref).
{\bf Representing combinations as binary vectors}. Let $[p] = \{1,\ldots p\}$ denote the set of $p$ interventions. Denote by $\Pi$ the power set of $[p]$, the set of all possible combinations of $p$ interventions, where we note $|\Pi| = 2^p$. Then, any given combination $\pi \in \Pi$ induces the following binary representation $\mathbf{v}(\pi) \in \{-1,1\}^p$ defined as follows: $\mathbf{v}(\pi)_{i} = 2~\mathbbm{1}\{i \in \pi\} -1.$
{\bf Fourier expansion of Boolean functions.} Let $\mathcal{F}_{\text{bool}} = \{f : \{-1,1\}^p \to \mathbb{R}\}$ be the set of all real-valued functions defined on the hypercube $\{-1,1\}^p$. Then $\mathcal{F}_{\text{bool}}$ forms a Hilbert space defined by the following inner product: for any $f,g \in \mathcal{F}_{\text{bool}}$, $\langle f, g \rangle_{B} = \frac{1}{2^p}\sum_{\mathbf{x} \in \{-1,1\}^p}f(\mathbf{x})g(\mathbf{x}).$ This inner product induces the norm $\langle f, f \rangle_{B} \coloneqq \lVert f \rVert_{B}^2 = \frac{1}{2^p}\sum_{\mathbf{x} \in \{-1,1\}^p}f^2(\mathbf{x}).$ We construct an orthonormal basis for $\mathcal{F}_{\text{bool}}$ as follows: for each subset $S \subset [p]$, define a basis function $\chi_{S}(\mathbf{x}) = \prod_{i \in S} x_i$ where $x_i$ is the $i^\text{th}$ coefficient of $\mathbf{x} \in \{-1,1\}^p$. One can verify that for any $S \subset [p]$ that $\lVert \chi_{S} \rVert_B = 1$, and that $\langle \chi_S, \chi_{S'} \rangle_B = 0$ for any $S' \neq S$. Since $|\{\chi_S: S \subset [p] \}| = 2^p$, the functions $\chi_S$ are an orthonormal basis of $\mathcal{F}_{\text{bool}}$. We will refer to $\chi_S$ as the Fourier character for the subset $S$.
Hence, any $f \in \mathcal{F}_{\text{bool}}$ can be expressed via the following Fourier decomposition: $f(\mathbf{x}) = \sum_{S \subset[p]}\alpha_{S}\chi_{S}(\mathbf{x}),$ where the Fourier coefficient $\alpha_S$ is given by computing $ \alpha_S = \langle f, \chi_S \rangle_B $. For a function $f$, we refer to $\boldsymbol{\alpha}_f = [\alpha_S]_{S \in [p]} \in \mathbb{R}^{2^p}$ as the vector of Fourier coefficients associated with it. For a given binary vector $\mathbf{x} \in \{-1, 1\}^{p}$, we refer to $\boldsymbol{\chi}(x) = [\chi_{S}(\mathbf{x})]_{S \in [p]} \in \{-1, 1\}^{2^p}$ as the vector of Fourier character outputs associated with it. Hence any function $f : \{-1,1\}^p \to \mathbb{R}$ can be re-expressed as follows: $f(\mathbf{x}) = \langle \boldsymbol{\alpha}_f, \boldsymbol{\chi}(x) \rangle.$ For $\pi \in \Pi$, abbreviate $\chi_S(\mathbf{v}(\pi))$ and $\boldsymbol{\chi}(\mathbf{v}(\pi))$ as $\chi^{\pi}_S$ and $\boldsymbol{\chi}^\pi$ respectively.
{\bf Observed and potential outcomes.} Let $Y^{(\pi)}_n \in \mathbb{R}$ denote the {\em potential outcome} for unit $n$ under combination $\pi$ and $Y_{n\pi} \in \{\mathbb{R} \cup \star \}$ as the {\em observed outcome}, where $\star$ indicates a missing value, i.e., the outcome associated with the unit-combination pair $(n, \pi)$ was not observed. Let $\mathbf{Y} = [Y_{n\pi}] \in \{\mathbb{R} \cup \star \}^{N \times 2^p}.$ Let $\mathcal{D} \subset [N] \times [2^p],$ refer to the subset of observed unit-combination pairs, i.e.,
Note that (ref) implies stable unit treatment value assignment (SUTVA) holds. Let $\Pi_S \subseteq \Pi$ denote a subset of combinations. For a given unit $n$, let $\mathbf{Y}_{\Pi_S,n} = [Y_{n\pi_i} : \pi_i \in \Pi_S] \in \{\mathbb{R} \cup \star \}^{|\Pi_S|}$ represent the vector of observed outcomes for all $\pi \in \Pi_S$. Similarly, let $\mathbf{Y}_n^{(\Pi_S)} = [Y^{(\pi_i)}_{n} : \pi_i \in \Pi_S] \in \mathbb{R}^{|\Pi_S|}$ represent the vector of potential outcomes. Denote $\boldsymbol{\chi}(\Pi_S) = [\boldsymbol{\chi}^{\pi_i} : \pi_i \in \Pi_S] \in \{-1, 1\}^{|\Pi^S| \times 2^p }$.
Define $Y_n^{(\cdot)} : \pi \to \mathbb{R}$ as a real-valued function over the hypercube $\{-1,1\}^p$ associated with unit $n$. It takes as input a combination $\pi$, converts it to a $p$-dimensional binary vector $\mathbf{v}(\pi)$, and outputs the real-valued potential outcome $Y_n^{(\pi)}$. Given the discussion in Section (ref) on the Fourier expansion of Boolean functions, it follows that $Y_n^{(\pi)}$ always has the representation $\langle \boldsymbol{\alpha}_n, \boldsymbol{\chi}^\pi \rangle$ for some $\boldsymbol{\alpha}_n \in \mathbb{R}^{2^p}$. Thus, without any loss of generality, the $\boldsymbol{\alpha}_n$ are unit-specific latent variables, i.e., the Fourier coefficients, encoding the treatment response function. Below, we state our key assumption on these induced Fourier coefficients.
The assumption is then that each $\boldsymbol{\alpha}_n$ is $s$-sparse and $\mathcal{A}$ is rank-$r$; $\epsilon_n^{\pi}$ is the residual from this sparse and low-rank approximation, and it serves as the source of uncertainty in our model.
Given $\mathbb{E}[\epsilon_n^{\pi} ~ | ~ \mathcal{A}] = 0$ and that $\mathcal{A}$ is rank-$r$, it implies the matrix $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}] = [\mathbb{E}[\mathbf{Y}_n^{(\Pi)}]: n \in [N]] \in \mathbb{R}^{2^p \times N}$ is also rank $r$, where the expectation is defined with respect to $\epsilon_n^{\pi}$. This is because $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}]$ can be written as $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}] = \boldsymbol{\chi}(\Pi) \mathcal{A}^T$, and since $\boldsymbol{\chi}(\Pi)$ is an invertible matrix, $\text{rank}(\mathbb{E}[\mathbf{Y}_N^{(\Pi)}]) = \text{rank}(\mathcal{A})$. The low-rank property places structure across units; that is, there is sufficient similarity across units so that $\mathbb{E}[\mathbf{Y}_n^{(\Pi)}]$ for any unit $n$ can be written as a linear combination of $r$ other rows of $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}]$. This is a standard assumption used to encode latent similarity across units in matrix completion and its related applications (e.g., recommendation engines).
Sparsity establishes unit-specific structure; that is, the potential outcomes for a given user only depend on a small subset of the Fourier characters $\{\chi_S: S \subset [p] \}$. This subset of Fourier characters can be different across units. As discussed in Section (ref), sparsity is commonly employed when studying the learnability of Boolean functions. In the context of recommendation engines, sparsity is implied if the ratings for a set of goods only depend on a small number of combinations of items within that set---if say only $k < p$ items mattered, then $s \le 2^k$. Sparsity is also often assumed implicitly in factorial design experiments, where analysts typically only include pairwise interaction effects between interventions and ignore higher-order interactions george2005statistics---here $s \le p^2$. These various models of imposing combinatorial structure in greater detail in Section (ref).
Next, we present an assumption that formalizes the dependence between the missingness pattern induced by the treatment assignments $\mathcal{D}$ and the potential outcomes $Y^{(\pi)}_n$, i.e., the type of confounding we can handle, and provide an interpretation of the induced data generating process (DGP) for potential outcomes.
Data Generating Process. Given Assumptions (ref) and (ref), one sufficient DGP that is line with the assumptions made can be summarized as follows: (i) unit-specific latent Fourier coefficients $\mathcal{A}$ are either deterministic or sampled from an unknown distribution; we will condition on this quantity throughout. (ii) Given $\mathcal{A}$, we sample mean-zero random variables $\epsilon_n^{\pi}$, and generate potential outcomes according to our model $Y_n^{(\pi)} = \langle \boldsymbol{\alpha}_n, \chi^{\pi} \rangle + \epsilon_n^{\pi}$. (iii) $\mathcal{D}$ is allowed to depend on unit-specific latent Fourier coefficients $\mathcal{A}$ (i.e., $\mathcal{D}=f(\mathcal{A})$). We define all expectations w.r.t. noise, $\epsilon_n^{\pi}$. This DGP introduces unobserved confounding since $Y_n^{(\pi)} \not\protect\mathpalette{\protect\independenT}{\perp} \mathcal{D}$. However, this DGP does imply that Assumption (ref) holds, i.e., conditional on the Fourier coefficients $\mathcal{A}$, the potential outcomes are independent of the treatment assignments $\mathcal{D}$. This conditional independence condition can be thought of as “selection on latent Fourier coefficients”, which is analogous to the widely made assumption of “selection on observables.” The latter requires that potential outcomes are independent of treatment assignments conditional on {\em observed} covariates---we reemphasize that $\mathcal{A}$ is unobserved.
Target parameter. For any unit-combination pair $(n,\pi)$, we aim to estimate $\mathbb{E}[Y^{(\pi)}_{n} \mid \mathcal{A}]$, where the expectation is w.r.t. $\epsilon_n^{\pi}$, and we condition on the set of Fourier coefficients $\mathcal{A}$.
We provide a brief summary of how the formalism established above for combinations can be extended to permuations (i.e., rankings), and provide a detailed discussion in Appendix (ref).
Binary Representation of Permutations. Let $\tau: [p] \rightarrow [p]$ denote a permutation on a set of $p$ items such that $\tau(i)$ denotes the rank of item $i \in [p]$. There are a total of $p!$ different permutations, and denote the set of permutations by $\mathbb{S}_p$. Every permutation $\tau \in \mathbb{S}_p$ induces a binary representation $\mathbf{v}(\tau) \in \{-1,1\}^{\binom{p}{2}}$, which can be constructed as follows. For an item $i$, define $\mathbf{v}^{i}(\tau) \in \{-1,1\}^{p-i}$ as follows: $\mathbf{v}_j^{i}(\tau)= \mathbbm{1}\{ \tau(i) > \tau(j)\} - \mathbbm{1}\{ \tau(i) < \tau(j)\}$ for items $1 \leq i < j \leq p$. That is, each coordinate $\mathbf{v}_j^{i}(\tau)$ indicates whether items $i$ and $j$ have been swapped for items $j > i$. Then, $\mathbf{v}(\tau) = [\mathbf{v}^{i}(\tau): i \in [p]] \in \{-1,1\}^{\binom{p}{2}}$. For example, with $p = 4$, the permutation $\tau([1,2,3,4]) = [1,3,4,2]$ has the binary representation $\mathbf{v}(\tau) = (-1,-1,-1,1,1,-1)$.
Fourier Expansion of Functions of Permutations. Since a permutation $\tau$ can be expressed as a binary vector, any function $f: \mathbb{S}_p \rightarrow \mathbb{R}$ can be thought of as a Boolean function. Then, given the discussion on Fourier expansions of Boolean functions in Section (ref), any function $f: \mathbb{S}_p \rightarrow \mathbb{R}$ admits the following Fourier decomposition: $f(\tau) = \sum_{S \subset [\binom{p}{2}]} \alpha_S \chi_{S}(\mathbf{v}(\tau)) \coloneqq \langle \boldsymbol{\alpha}_f, \boldsymbol{\chi}^{\tau} \rangle$, where $\boldsymbol{\alpha}_f = [\alpha_S]_{S \in [\binom{p}{2}]} \in \mathbb{R}^{\binom{p}{2}}$, and $\boldsymbol{\chi}^{\tau} = [\chi_S(\mathbf{v}(\tau)]_{S \in [\binom{p}{2}]} \in \{-1,1\}^{\binom{p}{2}}$ for $\tau \in \mathbb{S}_p$.
Model and DGP. We propose a similar model to the one discussed for combination, that is, model the potential outcome for a unit-permutation pair $(n,\tau)$ as $Y_n^{(\tau)} = \langle \boldsymbol{\alpha}_n, \boldsymbol{\chi}^{\tau} \rangle + \epsilon_n^{\tau}$, where we assume sparsity ($||\boldsymbol{\alpha}||_0 = s \leq 2^{\binom{p}{2}}$), low-rank structure ($\text{rank}(\mathcal{A}) = r \in \{\min\{N,2^{\binom{p}{2}}\}\}$), and $\mathbb{E}[\epsilon_n^{\tau} ~|~ \mathcal{A}] = 0$. Additionally, one sufficient DGP for permutations is one that is analogous to what is stated above for combinations.
Expressing potential outcomes over permutations $Y_n^{(\tau)}$ as Boolean functions allows us to easily adapt the proposed estimator, and our theoretical results to rankings. For simplicity, we focus on combinations for the rest of the paper and provide a detailed discussion of our results for permutations in Appendix (ref).
In this section, we discuss how classical models for learning functions over Booleans relates to the proposed potential outcome model. In particular, we discuss two well-studied models for functions over combinations: low-degree boolean polynomials and $k$-Juntas. We also discuss applications such as factorial design experiments and recommendation systems.
Low-degree boolean Polynomials. A special instance of sparse boolean functions is a low-degree polynomial which is defined as follows. For a positive integer $d \leq p$, $f : \{-1,1\}^p \rightarrow \mathbb{R}$ is a $d$-degree polynomial if its Fourier transform satisfies the following for any input $\mathbf{x} \in \{-1,1\}^p$
In this setting, $s \leq \sum^d_{i=0} \binom{p}{i} \approx p^d$. That is, degree $d$-polynomials impose sparsity on the potential outcome by limiting the degree of interaction between interventions. Next, we discuss how potential outcomes are typically modeled as low-degree polynomials in factorial design experiments.
Factorial Design Experiments. Factorial design experiments consist of $p$ treatments where each treatment arm can take on a discrete set of values, and units are assigned different combinations of these treatments. In the special case that each treatment arm only takes on two possible values, this experimental design mechanism is referred to as a $2^p$ factorial experiment. Factorial design experiments are widely employed in the social sciences, agricultural, and industrial applications duflo2007using,dasgupta2015causal,wu2011experiments, amongst others. A common strategy to determine treatment effects in this setting is to assume that the potential outcome only depends on main effects and pairwise interactions, with higher-order interactions being negligible bertrand2004emily,eriksson2014employers,george2005statistics. By only considering pairwise interactions, it is easy to see that $d = 2$, i.e., $s = p^2$. The reader can refer to zhao2022regression for a detailed discussion of various estimation strategies and their validity in different settings.
The model and algorithm proposed (see Assumption (ref)) capture these various modeling choices commonly used in factorial design experiments by imposing sparsity on $\boldsymbol{\alpha}_n$ and learning which of the coefficients are non-zero in a data-driven manner. That is, the Synthetic Combinations algorithm can adapt to low-degree polynomials without pre-specifying the degree, $d$, i.e., it automatically adapts to the inherent level of interactions in the data. However, the additional assumption made compared to the literature is that there is structure across units, i.e., the matrix $\mathcal{A} = [\boldsymbol{\alpha}_n]_{n \in [N]}$ is low-rank. Further, we require that units are observed under multiple combinations, which is the case in recommendation engines where users are exposed to multiple combinations of goods or in crossover design studies in which units receive a sequence of treatments vonesh1996linear,jones2003design.
$k$-Juntas. Another special case of sparse Boolean functions are $k$-Juntas which only depend on $k < p$ input variables. More formally, a function $f : \{-1,1\}^p \to \mathbb{R}$ is a $k$-junta if there exists a set $K \subset [p]$ with $|K| = k$ such that the fourier expansion of $f$ can be represented as follows
Therefore, in the setting of the $k$-Junta, the sparsity index $s \leq 2^k$. In contrast to low-degree polynomials where the function depends on all $p$ variables but limits the degree of interaction between variables, $k$-Juntas only depend on $k$ variables but allow for arbitrary interactions amongst them. We discuss how two important applications can be modeled as a $k$-Junta.
Recommendation Systems. Recommendation platforms such as Netflix are often interested in recommending a combination of movies that maximizes a user's engagement with a platform. Here, units $n$ can be individual users, and $\pi$ represents combinations of different movies. The potential outcomes are the engagement levels for a unit $n$ when presented with a combination of movies $\pi$, with $\boldsymbol{\alpha}_n$ representing the latent preferences for that user. In this setting, a user's engagement with the platform may only depend on a small subset of movies. For example, a user who is only interested in fantasy will only remain on the platform if they are recommended movies such as Harry Potter. Under this behavioral model, the potential outcomes (i.e., engagement levels) can be modeled as a $k$-Junta with the non-zero coefficients of $\boldsymbol{\alpha}_n$ representing the combinations of the $k$ movies that affect engagement level for a user $n$. The potential outcome observational model (Assumption (ref)) captures this form of sparsity while also reflecting the low-rank structure commonly assumed when studying recommendation systems. Once again, Synthetic Combinations can adapt to $k-$Juntas without pre-specifying the subset of features $K$, i.e., it automatically learns the important features for a given user. This additional structure in $k-$Juntas is well-captured by the CART estimator, which leads to tighter finite-sample bounds, as detailed in Section (ref).
Knock-down Experiments in Genomics. A key task in genomics is to identify which set of genes are responsible for a phenotype (i.e., physical trait of interest such as blood pressure) in a given patient. To do so, geneticists use knock-down experiments which measure the difference in the phenotype after eliminating the effect of a set of genes in an individual. To encode this process in the language of combinatorial causal inference, we can think of units as different individual patients or sub-populations of patients, an intervention as knocking a particular gene out, and $\pi$ as a combination of genes that are knocked out. The potential outcome $Y_n^{(\pi)}$ is the expression of the phenotype for a unit $n$ when the combination of genes $\pi$ are eliminated via knock-down experiments, and the coefficients of $\boldsymbol{\alpha}_n$ represent the effect of different combination of genes on the phenotype for unit $n$. A typical assumption in genomics is that phenotypes only depend on a small set of genes and their interactions. In this setting, one can model this form of sparsity by thinking of the potential outcome function as a $k$-junta, as well as capturing the similarity between the effect of genes on different patients via our low-rank assumption.
We show $\mathbb{E}[Y^{(\pi)}_{n} \mid \mathcal{A}]$ can be written as a function of observed outcomes, i.e., we establish identification of our target causal parameter. As discussed earlier, our model allows for unobserved confounding: whether or not a unit is seen under a combination may be correlated with its potential outcome under that combination due to unobserved factors, as long as certain conditions are met. We introduce necessary notation and assumption required for our result. For a unit $n \in [N]$, denote the subset of combinations we observe them under as $\Pi_n \subseteq \Pi$. For $\pi \in \Pi$, let $\tilde{\boldsymbol{\chi}}^{\pi}_n \in \mathbb{R}^{2^p}$ denote the restricted Fourier characteristic vector where we zero out all coordinates of $\boldsymbol{\chi}^{\pi} \in \mathbb{R}^{2^p}$ that correspond to the coefficients of $\boldsymbol{\alpha}_n$ which are zero. For example, if $\boldsymbol{\alpha}_n = (1,1,0,0,\ldots 0)$ and $\boldsymbol{\chi}^{\pi} = (1,1,\ldots 1)$, then $\tilde{\boldsymbol{\chi}}^{\pi}_n = (1,1,0, \ldots 0)$. We then make the following assumption.
Horizontal span inclusion requires that the set of observed combinations for any “donor” unit is “diverse” enough that the projection of the Fourier characteristic for a target intervention is in the span of characteristics of observed interventions for that unit. Linear span inclusion requires that the donor unit set is diverse enough such that the Fourier coefficient of any “non-donor” unit is in the span of the Fourier coefficients of the donor set.
Motivating example. The following serves as motivation for existence of a donor set $\mathcal{I}$ satisfying Assumption (ref). Suppose Assumption (ref) holds and that each unit belongs to one of $r$ types. That is, for all $n \in [N]$, $\boldsymbol{\alpha}_n \in \{\boldsymbol{\alpha}(1),\boldsymbol{\alpha}(2),\ldots,\boldsymbol{\alpha}(r)\} \subset \mathbb{R}^{2^p},$ where each $\boldsymbol{\alpha}(i)$ is $s$-sparse. Having only $r$ distinct sets of Fourier coefficients implies the low-rank property in Assumption (ref) (a). Suppose the donor set $\mathcal{I} \subset [N]$ is chosen by sampling a subset of units independently and uniformly at random with size satisfying $\Omega(r\log(r/\gamma))$. Similarly, sample $\Pi_\mathcal{I} \subset \Pi$ combinations independently and uniformly at random with size satisfying $\Omega(s\log(r|\mathcal{I}|/\gamma))$, and assign all donor units this set of combinations. Then, the following result shows this sampling scheme ensures both horizontal and linear span inclusion are satisfied with high probability as long as there are $\Omega(N)$ units of each type.
Proposition (ref) shows that (ignoring logarithmic factors) sampling $r$ units and assigning these units $s$ combinations at random is sufficient to induce a valid donor set. In Section (ref), we show that Assumption (ref) holds even when we relax our assumption that there are only $r$ types of units. Specifically, it is established that a randomly selected donor set of size $\omega(r\log(rs))$ will satisfy linear span inclusion as long as the $r$ non-zero singular values of the matrix $\mathcal{A}$ are of similar magnitude (see Assumption (ref) below). Note that the motivating example discussed above does not induce unobserved confounding since the treatment mechanism is completely randomized (i.e., $\mathcal{D}$ does not depend on $\mathcal{A}$). Hence, next we present a simple but representative example that both induces unobserved confounding and satisfies Assumption (ref).
Natural model of unobserved confounding. For any unit $n$, suppose $\boldsymbol{\alpha}_n \in \{\boldsymbol{\alpha}(1) = (1, 1, 4, 0 \ldots, 0), \ \boldsymbol{\alpha}(2) =(1, 1, 5, 0 \ldots, 0), \ \boldsymbol{\alpha}(3) = (0.5, 0.5, 2, 0 \ldots, 0), \ \boldsymbol{\alpha}(4) = (0.5, 0.5, 2.5, 0 \ldots, 0)\}$. One can verify $\text{rank}(\mathcal{A}) = 2 $. Define the treatment assignment to be such that combination $\pi$ for unit $n$ is observed only if $| \langle \boldsymbol{\alpha}_n, \boldsymbol{\chi}^{\pi} \rangle | \geq 2$ (i.e., $| \mathbb{E}[Y^{(\pi)}_{n}] | \ge 2$). Missingness patterns where only outcomes with large absolute values are observed is common in applications such as recommendation engines, where one is only likely to observe ratings for combinations that users either strongly like or dislike. For units of type $\{\boldsymbol{\alpha}(3),\boldsymbol{\alpha}(4)\}$, we verify in Appendix (ref) that they do not satisfy horizontal span inclusion (i.e., Assumption (ref) (a) does not hold). For units of type $\{\boldsymbol{\alpha}(1),\boldsymbol{\alpha}(2)\}$, combinations $\{\pi_1,\pi_2,\pi_3,\pi_4\}$ with the following binary representations $\{v(\pi_1) = (1,1,1,1,\ldots,1), v(\pi_2) = (1,-1,-1,1,\ldots,1)\},v(\pi_3) = (1,-1,1,1,\ldots,1),v(\pi_4) = (1,1,-1,1,\ldots,1)$ are observed. These combinations have associated restricted Fourier characteristics $\Tilde{\boldsymbol{\chi}}^\pi$ as follows: $\{\Tilde{\boldsymbol{\chi}}^{\pi_1} = (1,1,1,0,\ldots,0), \Tilde{\boldsymbol{\chi}}^{\pi_4} = (1,-1,-1,0,\ldots,0)\}, \Tilde{\boldsymbol{\chi}}^{\pi_3} = (1,-1,1,0,\ldots,0),\Tilde{\boldsymbol{\chi}}^{\pi_2} = (1,1,-1,0,\ldots,0)$. One can verify that these observed combinations ensure that horizontal span inclusion holds for units with type $\{\boldsymbol{\alpha}(1),\boldsymbol{\alpha}(2)\}$. Since $\text{rank}\{\boldsymbol{\alpha}(1),\boldsymbol{\alpha}(2)\} =2$, vertical span inclusion also holds. It is straightforward to generalize this example to setting with $r$ types.
Given these assumptions, we now present our identification theorem.
Theorem (ref) gives conditions under which the donor set $\mathcal{I}$ and the treatment assignments $\mathcal{D}$ are sufficient to recover the full set of unit specific potential outcomes $\mathbb{E}[Y_n^{(\pi)}|\mathcal{A}]$ in the noise-free limit. Part (a) establishes that for every donor unit $u \in \mathcal{I}$, the causal estimand can be written as a function of its own observed outcomes $\mathbb{E}[\mathbf{Y}_{\Pi_u}]$, given knowledge of $\boldsymbol{\beta}^{\pi}_{\Pi_u}$. Part (b) states that the target causal estimand $\mathbb{E}[Y_n^{(\pi)}]$ for a non-donor unit and combination $\pi$ can be written as a linear combination of the outcomes of the donor set $\mathcal{I}$, given knowledge of $\mathbf{w}_n^n$. Previous work that establishes identification under a latent factor model requires a growing number of donor units to be observed under all treatments agarwal2021causalmatrix. This is infeasible in our setting because the vast majority of combinations have no units that receive it. As a result, one has to first identify the outcomes of donor units under all combinations (part (a)), before transferring them to non-donor units (part (b)). In order to do so, Theorem (ref) suggests that the key quantities in estimating $\mathbb{E}[Y_n^{(\pi)}]$ for any unit-combination pair $(n,\pi)$ are $\boldsymbol{\beta}_{\Pi_u}^{\pi}$ and $\mathbf{w}_u^n$. In the following section, we propose an algorithm to estimate both $\boldsymbol{\beta}_{\Pi_u}^{\pi}$ and $\mathbf{w}_u^n$, as well as concrete ways of determining the donor set $\mathcal{I}$.
We now describe the Synthetic Combinations estimator, a simple and flexible two-step procedure for estimating our target causal parameter. A pictorial representation of the estimator is presented in Figure (ref).
Step 1: Horizontal Regression. For notational simplicity, denote the vector of observed responses $\mathbf{Y}_{n,\Pi_n} = [Y_{n\pi} : \pi \in \Pi_n] \in \mathbb{R}^{|\Pi_n|}$ for any unit $n$ as $\mathbf{Y}_{\Pi_n}$. Then, for every unit $u$ in the donor set $\mathcal{I}$, $\mathbb{E}[Y_u^{(\pi)}]$ is estimated via the Lasso, i.e., by solving the following convex program with penalty parameter $\lambda_u$:
where recall that $\boldsymbol{\chi}(\Pi_u) = [\boldsymbol{\chi}^\pi: \pi \in \Pi_u] \in \mathbb{R}^{|\Pi_u| \times 2^p}$. Then, for any donor unit-combination pair $(u,\pi)$, let $\hat{\mathbb{E}}[Y_u^{(\pi)}] = \langle \hat{\boldsymbol{\alpha}}_u, \boldsymbol{\chi}^{\pi} \rangle$ denote the estimate of the potential outcome $\mathbb{E}[Y_u^{(\pi)}]$.
Step 2: Vertical Regression. Next, estimate potential outcomes for all units $n \in [N] \setminus \mathcal{I}$. To do so, some additional notation is required. For $\Pi_S \subseteq \Pi$, define the vector of estimated potential outcomes $\hat{\mathbb{E}}[\mathbf{Y}^{(\Pi_S)}_{u}] = [ \hat{\mathbb{E}}[Y_u^{(\pi)}]: \pi \in \Pi^S] \in \mathbb{R}^{|\Pi_S|}$. Additionally, let $\hat{\mathbb{E}}[\mathbf{Y}^{(\Pi_S)}_{\mathcal{I}}] = [\hat{\mathbb{E}}[\mathbf{Y}^{(\Pi_S)}_{u}]: u \in \mathcal{I}] \in \mathbb{R}^{|\Pi_S| \times |\mathcal{I}|}$.
Step 2(a): Principal Component Regression. Perform a singular value decomposition (SVD) of $\hat{\mathbb{E}}[\mathbf{Y}^{(\Pi_n)}_{\mathcal{I}}]$ to get $\hat{\mathbb{E}}[\mathbf{Y}^{(\Pi_n)}_{\mathcal{I}}] = \sum^{\min(|\Pi_n|,|\mathcal{I}|)}_{l = 1} \hat{s_l}\hat{\boldsymbol{\mu}}_{l}\hat{\boldsymbol{\nu}}^T_{l}$. Using a hyper-parameter $\kappa_n \leq \min(|\Pi_n|,|\mathcal{I}|)$ \footnote{ Both $\lambda_u$ and $\kappa_n$ can be chosen in a data-driven manner (e.g., via cross-validation).}, compute $\hat{\mathbf{w}}^{n} \in \mathbb{R}^{|\mathcal{I}|}$ as follows:
Step 2(b): Estimation. Using $\hat{\mathbf{w}}^{n} = [\hat{w}_u^n : u \in \mathcal{I}]$, we have the following estimate for any intervention $\pi \in \Pi$
Suitability of Lasso and PCR. Lasso is appropriate for the horizontal regression step because it adapts to the sparsity of $\boldsymbol{\alpha}_u$. However, it can be replaced with other algorithms that adapt to sparsity (see below for a larger discussion). For vertical regression, PCR is appropriate because $\mathcal{A}$ is low rank. As agarwal2019robustness,agarwal2020principal show, PCR implicitly regularizes the regression by adapting to the rank of the covariates, i.e., $\mathbb{E}[\mathbf{Y}^{(\Pi_n)}_{\mathcal{I}}]$. As a result, the out-of-sample error of PCR scales with $r$ rather than the ambient covariate dimension, which is given by $\mathcal{I}$.
Determining Donor Set $\mathcal{I}$. Synthetic Combinations requires the existence of a subset of units $\mathcal{I} \subset [N]$ such that one can (i) accurately estimate their potential outcomes under all possible combinations, and (ii) transfer these estimated outcomes to a unit $n \in [N] \setminus \mathcal{I}$. Theoretically, we give sufficient conditions on the observation pattern such that we can perform (i) and (ii) accurately via the Lasso and PCR, respectively. In practice, the following practical guidance to determine the donor set $\mathcal{I}$ is recommended. For every unit $n \in [N]$, learn a separate Lasso model $\hat{\boldsymbol{\alpha}}_n$ and assess its performance through $k$-fold cross-validation (CV). Assign units with low CV error (with a pre-determined threshold) into the donor set $\mathcal{I}$, and estimate outcomes $\hat{\mathbb{E}}[Y_u^{(\pi)}]$ for every unit $u \in \mathcal{I}$ and $\pi \in \Pi$. For non-donor units, PCR performance can also be assessed via k-fold CV. For units with low PCR error, linear span inclusion (Assumption (ref)(b)) and the assumptions required for the generalization for PCR likely hold, and hence we estimate their potential outcomes as in (ref). For units with large PCR error, it is either unlikely that this set of assumptions holds or that $|\Pi_n|$ is not large enough (i.e., additional experiments need to be run for this unit), and hence we do not recommend estimating their counterfactuals. In our real-world empirics in Section (ref), we choose donor units via thie approach, and find Synthetic Combinations outperforms other methods.
Horizontal Regression Model Selection. Synthetic Combinations allows for any ML algorithm (e.g., random forests, neural networks, ensemble methods) to be used in the first step. We provide an example of this flexibility by showing how the horizontal regression can also be done via CART in Section (ref). Theorem (ref) shows CART leads to better finite-sample rates when $\mathbb{E}[Y_n^{(\pi)}]$ is a $k$-Junta. This model-agnostic approach allows the analyst to tailor the horizontal learning procedure to the data at hand and include additional structural information for better performance.
In this section, finite-sample consistency of Synthetic Combinations is established, starting with a discussion of the additional assumptions required for the results.
To define the next set of assumptions,necessary notation is introduced. For any subset of combinations $\Pi_S \subset \Pi$, let $\mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_S)}] = [\mathbb{E}[\mathbf{Y}_u^{(\Pi_S)}]: u \in \mathcal{I}] \in \mathbb{R}^{|\Pi_S| \times |\mathcal{I}|}$.
Assumption (ref) is necessary for finite-sample consistency when estimating $\boldsymbol{\alpha}_n$ via the Lasso estimator (ref), and is commonly made when studying the Lasso rigollet2015high. Incoherence can also be seen as an inclusion criteria for a unit $n$ to be included in the donor set $\mathcal{I}$. Lemma 2 in negahban2012learning shows that the Assumption (ref) is satisfied with high probability if $\Pi_u$ is chosen uniformly at random and grows as $\omega(s^2p)$. If Assumption (ref) does not hold, then the Lasso may not estimate $\boldsymbol{\alpha}_u$ accurately, and alternative horizontal regression algorithms may be required instead. Assumption (ref) requires that the non-zero singular values of $\mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)} ~ | ~ \mathcal{A}]$ are well-balanced. This assumption is standard when studying PCR agarwal2019robustness,agarwal2020synthetic, and within the econometrics literature bai2021matrix,fan2018eigenvector. It can also be empirically validated by plotting the spectrum of $\hat{\mathbb{E}}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$; if the singular spectrum of $\hat{\mathbb{E}}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$ displays a natural elbow point, then Assumption (ref) is likely to hold. Assumption (ref) is also commonly made when analyzing PCR agarwal2020principal,agarwal2020synthetic,agarwal2021causal. It can be thought of as a “causal transportability” condition from the model learned using $\Pi_n$ to the interventions $\pi \in \Pi \setminus \Pi_n$. That is, subspace inclusion enables accurate estimation of $\mathbb{E}[Y_n^{(\pi)} ~ | ~ \mathcal{A} ]$ using $\langle \hat{\mathbb{E}}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}], \hat{\mathbf{w}}^n \rangle$. In Section (ref), we propose a simple experimental design mechanism that ensures that Assumptions (ref), (ref), (ref) (and Assumption (ref))---the key conditions for Synthetic Combinations to work---hold with high probability.
The following result establishes finite-sample consistency of Synthetic Combinations. Without loss of generality, we will focus on estimating the pair of quantities $(\mathbb{E}[Y_u^{(\pi)}],\mathbb{E}[Y_n^{(\pi)}])$ for a given donor unit $u \in \mathcal{I}$, and non-donor unit $n \in [N] \setminus \mathcal{I}$ for combination $\pi \in \Pi$. To simplify notation, we utilize $O_p$ notation: for any sequence of random vectors $X_n$, $X_n$ = $O_p(\gamma_n)$ if, for any $\epsilon > 0$, there exists constants $c_\epsilon$ and $n_\epsilon$ such that $\mathbb{P}(\lVert X_n \rVert_2 \geq c_\epsilon\gamma_n) \leq \epsilon$ for every $n \geq n_\epsilon$. Additionally, let $O_p$ absorb any dependencies on $\sigma$. For notational simplicity, define $\tilde{O_p}(\gamma_n)$ which suppresses logarithmic terms.
Establishing Theorem (ref) requires a novel analysis of error-in-variables (EIV) linear regression. Specifically, the general EIV linear model is as follows: $Y = \mathbf{X} \beta + \epsilon$, $\mathbf{Z} = \mathbf{X} + \mathbf{H}$, where $Y$ and $\mathbf{Z}$ are observed. In our case, $Y = \mathbf{Y}_{\Pi_n}$, $\mathbf{X} = \mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$, $\beta = \mathbf{w}^n$ , $\mathbf{Z} = \hat{\mathbb{E}}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$, and $\mathbf{H}$ is the error arising in estimating $\mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$ via the Lasso. Typically one assumes that $\mathbf{H}$ is is a matrix of independent sub-gaussian noise. Our analysis requires a novel worst-case analysis of $\mathbf{H}$ (due to the 2-step regression of Lasso and then PCR), in which each entry of is $\mathbf{H}$ simply bounded.
We describe the conditions required on $|\Pi_n|$, $M$, for Synthetic Combinations to consistently estimate $(\mathbb{E}[Y_u^{(\pi)}], \ \mathbb{E}[Y_n^{(\pi)}])$. Recall $|\Pi_n|$ are the number of observations for the non-donor unit $n$ of interest, and $M$ is the minimum number of observations for a donor unit. If $|\Pi_n| = \omega(r^4_n)$ and $M = \omega(r^4_ns^2p)$, then $\max\left(|\hat{\mathbb{E}}[Y_u^{(\pi)}] - \mathbb{E}[Y_{u}^{(\pi)}]|, \\ |\hat{\mathbb{E}}[Y_n^{(\pi)}] - \mathbb{E}[Y_{n}^{(\pi)}]|\right) =\tilde{o}_p(1)$. Conversely, the corollary below quantifies how quickly the parameters $r_n$, $s$, and $p$ can grow with the number of observations $|\Pi_n|$ and $M$.
Corollary (ref) quantifies how $s, r_n$ can scale with the number of observations to achieve consistency. That is, it reflects the maximum “complexity” allowed for a given sample size.
As discussed in Section (ref), Synthetic Combinations allows for any ML algorithm to be used for horizontal regression. In this section, we present theoretical results when the horizontal regression is done via CART. When the potential outcome function is a $k$-Junta (see equation (ref) for the definition), CART is able to achieve stronger sample complexity guarantees than the Lasso. See Appendix (ref) for a detailed description of the CART algorithm. Below, a brief overview of the CART-based horizontal regression procedure is given.
Step 1: Feature Selection. For each donor unit $u \in \mathcal{I}$, divide the observed interventions $\Pi_u$ equally into two sets $\Pi^{a}_u$ and $\Pi^{b}_u$. Fit a CART model $\mathcal{\hat{T}}$ on the dataset $\mathcal{D}^a_u = [\{\mathbf{v}(\pi),Y_{u\pi}\}: \pi \in \Pi^a_u]$. Let $\hat{K}_u = \{\mathbf{v}(\pi)_j : \mathbf{v}(\pi)_j$ $\in \mathcal{\hat{T}}\}$ denote every feature that the CART splits on.
Step 2: Estimation. For every subset $S \subset \hat{K}_u$, compute $\hat{\alpha}_{u,S} = \frac{1}{|\Pi^b_u|}\sum_{\pi \in \Pi^b_u }Y_{u\pi}\chi^{\pi}_{S}$. The target causal parameter is then estimated as follows: $\hat{\mathbb{E}}[Y_u^{(\pi)}] = \sum_{S \in \hat{K}_u}\hat{\alpha}_{u,S}\chi^{\pi}_{S}$.
CART is able to take advantage of the $k$-Junta structure of the potential outcomes by first performing feature selection (i.e., learning the relevant feature set $\hat{K}_u$), before estimating the Fourier coefficient for each relevant subset. In contrast, Lasso estimates the Fourier coefficient of all $2^p$ subsets simultaneously.
Next, we present an informal theorem for the horizontal regression finite sample error when it is done via CART. The formal result (Theorem (ref)) can be found in Appendix (ref).
Compared to the Lasso horizontal regression error (see Theorem (ref) (a)), CART removes an additional factor of $p$ by exploiting the $k$-Junta structure. It can be verified that CART also removes a factor of $p$ for the vertical regression error (see Theorem (ref) (b)). As a result, if $|\Pi_n| = \omega(r^4_n)$, and $M = \omega(r^4_ns^2)$, then $\max\left(|\hat{\mathbb{E}}[Y_u^{(\pi)}] - \mathbb{E}[Y_{u}^{(\pi)}]|, \ |\hat{\mathbb{E}}[Y_n^{(\pi)}] - \mathbb{E}[Y_{n}^{(\pi)}]|\right) =\tilde{o}_p(1)$. That is, CART reduces the number of required observations for donor units by a factor of $p$.
We discuss the sample complexity of Synthetic Combinations to estimate all $N \times 2^p$ causal parameters, and compare it to that of other methods. To ease our discussion, we ignore dependence on logarithmic factors, and assume that the horizontal regression is done via the Lasso.
Even if potential outcomes $Y_n^{(\pi)}$ were observed for all unit-combination pairs, consistently estimating $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}]$ is not trivial. This is because, we only get to observe a single and noisy version $Y_{n\pi} = \langle \boldsymbol{\alpha}_n, \boldsymbol{\chi}^{\pi} \rangle + \epsilon_n^{\pi}$. Hypothetically, if $K$ independent samples of $Y_{n\pi}$ for a given $(n,\pi)$, denoted by $Y^{1}_{n\pi}, \ldots, Y^{K}_{n\pi}$ are observed, then the maximum likelihood estimator would be the empirical average $\frac{1}{K}\sum^K_{i=1}Y^{i}_{n\pi}$. The empirical average would concentrate around $\mathbb{E}[Y_n^{(\pi)}]$ at a rate $O(1/\sqrt{K})$ and hence would require $K = \Omega(\delta^{-2})$ samples to estimate $\mathbb{E}[Y_n^{(\pi)}]$ within error $O(\delta)$. Therefore, this naive (unimplementable) solution would require $N \times 2^p \times \delta^{-2}$ observations to estimate $\mathbb{E}[\mathbf{Y}_N^{(\Pi)}]$.
On the other hand, Synthetic Combinations produces consistent estimates of the potential outcome despite being given {\em at most only a single noisy sample} of each potential outcome. As the discussion after Theorem (ref) shows, Synthetic Combinations requires $|\mathcal{I}| \times r^4s^2p/\delta^2$ observations for the donor set, and $(N - |\mathcal{I}|) \times r^4/\delta^4$ observations for the non-donor units to achieve an estimation error of $O(\delta)$ for all $N \times 2^p$ causal parameters. To satisfy linear span inclusion (Assumption (ref) (b)), and the necessary conditions for PCR (i.e., Assumption (ref)), the donor set needs to have size $|\mathcal{I}| = \omega(r)$. Our experimental design mechanism in Section (ref) shows that uniformly choosing a donor set of size $|\mathcal{I}| = O(r)$ satisfies these assumptions. Hence, the number of observations required to achieve an estimation error of $O(\delta)$ for all pairs $(n,\pi)$ scales as $O\left(\text{poly}(r/\delta) \times \left(N + s^2p \right) \right)$.
Sample Complexity Comparison to Other Methods.
Horizontal regression: An alternative algorithm would be to learn an individual model for each unit $n \in [N]$. That is, run a separate horizontal regression via the Lasso for every unit. This alternative algorithm has sample complexity that scales at least as $O(N \times s^2p/\delta^2)$ rather than $O\left(\text{poly}(r)/\delta^4 \times \left(N + s^2p \right) \right)$ required by Synthetic Combinations. It suffers because it does not utilize any structure across units (i.e., the low-rank property of $\mathcal{A}$), whereas Synthetic Combinations captures the similarity between units via PCR.
Matrix completion: Synthetic Combinations can be thought of as a matrix completion method; estimating $\mathbb{E}[Y_n^{(\pi)}]$ is equivalent to imputing $(n,\pi)$-th entry of the observation matrix $\mathbf{Y} \in \{\mathbb{R} \cup \star\}^{N \times 2^p}$, where recall $\star$ denotes an unobserved unit-combination outcome. Under the low-rank property (Assumption (ref)(b)) and various models of missingness (i.e., observation patterns), recent works on matrix completion candes2010matrix,ma2019missing,agarwal2021causalmatrix (see Section (ref) for an overview) have established that estimating $\mathbb{E}[Y_n^{(\pi)}]$ to an accuracy $O(\delta)$ requires at least $O\left(\text{poly}(r/\delta) \times \left(N + 2^p \right)\right)$ samples. This is because matrix completion techniques do not leverage the sparsity of $\boldsymbol{\alpha}_n$. Moreover, matrix completion results typically report error in the Frobenius norm, whereas we give entry-wise guarantees. This leads to an extra factor of $s$ in our analysis as it requires proving convergence for $\lVert \hat{\boldsymbol{\alpha}_n} - \boldsymbol{\alpha}_n \rVert_{1}$ rather than $\lVert \hat{\boldsymbol{\alpha}_n} - \boldsymbol{\alpha}_n \rVert_{2}$.
Natural Lower Bound on Sample Complexity. We provide an informal discussion on the lower bound sample-complexity to estimate all $N \times 2^p$ potential outcomes. As established in Lemma (ref), $\mathcal{A}$ has at most $rs$ non-zero columns. Counting the parameters in the singular value decomposition of $\mathcal{A}$, only $r\times (N + rs)$ free parameters are required to be estimated. A natural lower bound on the sample complexity scales as $O(Nr + r^2s)$. Hence, Synthetic Combinations is only sub-optimal by a factor (ignoring logarithmic factors) of $sp$ and $\text{poly}(r)$. As discussed earlier, an additional factor of $s$ can be removed if we focus on deriving Frobenius norm error bounds. Further, the use of non-linear methods such as CART can remove the factor of $p$ (Theorem (ref)) under stronger assumptions on the potential outcomes, e.g., $\mathbb{E}[Y_n^{(\pi)}]$ is a $k$-Junta. It remains as interesting future work to derive estimation procedures that are able to achieve this lower bound.
In this section, we show how Synthetic Combinations can be used to design experiments that allow for combinatorial inference (i.e., learning all $N \times 2^p$ causal parameters).
Key Assumptions Behind Finite Sample Consistency of Synthetic Combinations. Synthetic Combinations requires the existence of a donor set $\mathcal{I}$ such that we are able to perform accurate horizontal regression for all donor units, and then transfer these estimated outcomes to non-donor units via PCR. The enabling conditions for accurate horizontal regression are (i) horizontal span inclusion (Assumption (ref) (a)), and (ii) incoherence of the Fourier characteristics (Assumption (ref)). Similarly, the critical assumptions required for consistency of PCR are (i) linear span inclusion (Assumption (ref) (b)), (ii) well-balanced spectrum (Assumption (ref)), and (iii) subspace inclusion (Assumption (ref)). In terms of experiment design, the donor set and treatment assignments must be carefully chosen such that these key assumptions hold. To that end, we introduce the following design.
Experimental Design Mechanism. Fix a probability threshold $\gamma \in (0,1)$ and estimation error threshold $\delta \in (0,1)$. Our design mechanism (see Figure (ref) for a a visual description) then proceeds as follows.
Step 1: Donor set selection. Choose the donor set $\mathcal{I} \subset [N]$ by sampling a subset of units independently and uniformly at random with size satisfying $\Omega\left(r\log(rs/\gamma) \right)$.
Step 2: Donor set treatment assignment. Sample $\Pi_{\mathcal{I}} \subset \Pi$ combinations independently and uniformly at random with size satisfying $|\Pi_{\mathcal{I}}| = \Omega\left(\frac{r^3s^2\log(|\mathcal{I}|2^p/\gamma)}{\delta^2}\right)$. Assign all donor units $u \in \mathcal{I}$ this set of combinations.
Step 3: Non-donor unit treatment assignment. Randomly sample $\Pi_N \subset \Pi$ combinations independently and uniformly at random of size $|\Pi_N| = \Omega\left(r\log(|\mathcal{I}|/\gamma) \vee r^4/\delta^4 \right)$. Assign all non-donor units $n \in [N] \setminus \mathcal{I}$ combinations $\Pi_N$.
Potential outcome estimation. Given the observation pattern described above, estimate outcomes for each unit-combination pair via Synthetic Combinations as described in Section (ref).
It turns out that this simple design mechanism satisfies the key assumptions presented above with high probability. In fact, the proposed mechanism ensures these key conditions hold under very limited assumptions. The only required assumptions are that: (i) the potential outcome model is satisfied (Assumption (ref)), (ii) $\mathbb{E}[Y_n^{(\pi)}]$ is bounded (Assumption (ref)). Additionally, only a weakened version of the balanced spectrum condition is needed (Assumption (ref)):
As compared to the original balanced spectrum condition, Assumption (ref) only requires that the singular values are balanced for the entire potential outcome matrix as opposed to a collection of submatrices. We then have the following result.
Theorem (ref) implies that the key enabling conditions for Synthetic Combinations are satisfied for every unit-combination pair $(n,\pi)$. Corollary (ref) further establishes that with this design, $O(\delta)$ error is achievable for all $N \times 2^p$ causal parameters. In contrast, the results in the observational setting do not guarantee accurate estimation of {\em all} $N \times 2^p$ parameters. Instead, they establish that $\mathbb{E}[Y_n^{(\pi)}]$ can be learned for any specific unit-combination pair $(n,\pi)$ that satisfies the required assumptions on the observation pattern. Additionally, given the discussion in Section (ref), one can verify that the number of observations required by this experiment design mechanism scales as $\tilde{O}\left(\text{poly}(r/\delta)\times \left(N + s^2p \right) \right)$. In practice, this design requires knowledge of $r,$ and $s$. To overcome this, one can sequentially sample donor units and their observations until the rank and lasso error stabilizes. This provides an estimate of $r,$ and $s$; a formal analysis of this procedure is left as future work.
We now establish asymptotic normality of Synthetic Combinations under additional conditions. Since Synthetic Combinations is agnostic to the learning algorithm used in the horizontal regression, as demonstrated by the discussion on CART, the Lasso can be replaced by any regression technique that achieves asymptotic normality. For example, previous works have studied variants of the Lasso liu2013asymptotic,zou2006adaptive, where it has been established that the predictions are asymptotically normal. Here, we propose replacing the Lasso with the “Select + Ridge” (SR) procedure of liu2013asymptotic, and specify the conditions on the donor units' observation pattern under which SR achieves asymptotic normality. Next, using asymptotic normality of the horizontal regression predictions, asymptotic normality of the vertical regression step of Synthetic Combinations is established.
The SR estimator consists of the following two steps.
Step 1: Subset Selection. For every unit $u$ in the donor set $\mathcal{I}$, fit a Lasso model with regularization parameter $\lambda_u$ using $\boldsymbol{\chi}(\Pi_u)$ as described in Step 1 of Synthetic Combinations to obtain $\hat{\boldsymbol{\alpha}}_u$. Let $\hat{\mathcal{S}}_u = [S \subset [p]: \hat{\alpha}_{u,S} \neq 0]$ denote the non-zero coefficients of $\hat{\boldsymbol{\alpha}}_u$.
Some additional notation needs to be defined for the next step. For a unit $u \in \mathcal{I}$, and combination $\pi$, let $\boldsymbol{\chi}_{\hat{\mathcal{S}}_u}^{\pi} = [\chi^{\pi}_{S} : S \in \hat{\mathcal{S}}_u] \in \{-1,1\}^{|\hat{\mathcal{S}}_u|}$. Let $\boldsymbol{\chi}_{\hat{\mathcal{S}}_u}(\Pi_u) = [\boldsymbol{\chi}_{\hat{\mathcal{S}}_u}^{\pi}: \pi \in \Pi_{u}] \in \{-1,1\}^{|\Pi_{u}| \times |\hat{\mathcal{S}}_u|}$.
Step 2: Ridge Regression. In this step, SR forms predictions using the selected subsets in Step 1, and ridge regression. Specifically, for every $u \in \mathcal{I}$, compute the Fourier coefficient as follows,
where $\mathbf{I}_{|\hat{\mathcal{S}}_u|} \in \mathbb{R}^{|\hat{\mathcal{S}}_u| \times |\hat{\mathcal{S}}_u|}$ is the identity matrix of size $|\hat{\mathcal{S}}_u|$. Then, let $\hat{\mathbb{E}}_{SR}[Y_u^{(\pi)}] = \langle \hat{\boldsymbol{\alpha}}^{SR}_{u}, \boldsymbol{\chi}^{\pi} \rangle$ denote the estimated potential outcome $\mathbb{E}[Y_u^{(\pi)}]$ for a combination $\pi$.
Next, we present our result establishing asymptotic normality of SR predictions. Some additional notation is required for our result. For a given unit $u \in \mathcal{I}$, let $\mathcal{S}_u = \{S \subset [p] ~ | ~ \alpha_{u,S} \neq 0 \}$. Like above, define $\boldsymbol{\chi}_{\mathcal{S}_u}^{\pi} = [\chi^{\pi}_{S} : S \in \mathcal{S}_u] \in \{-1,1\}^{|\mathcal{S}_u|}$ and $\boldsymbol{\chi}_{\mathcal{S}_u}(\Pi_u) = [\boldsymbol{\chi}_{\mathcal{S}_u}^{\pi}: \pi \in \Pi_{u}] \in \{-1,1\}^{|\Pi_{u}| \times |\mathcal{S}_u|}$. Denote $\mathbf{K}_{u} = \frac{1}{|\Pi_u|}((\boldsymbol{\chi}_{\mathcal{S}_u}(\Pi_u))^T\boldsymbol{\chi}_{\mathcal{S}_u}(\Pi_u))$ as the covariance matrix of the observed Fourier characteristics for a donor unit $u$. For a square matrix $\mathbf{X}$, let $\lambda_{\min}(\mathbf{X})$ denote its smallest eigenvalue.
Proposition (ref) builds upon Theorem 3 of liu2013asymptotic to establish asymptotic normality and as a result allows for construction of valid confidence intervals for the donor units. Condition (b) can be enforced by simply normalizing each row of $\boldsymbol{\chi}(\Pi_u)$ to have mean zero. Condition (c) is mild, and is required in order to ensure that the model is identifiable even if $\mathcal{S}_u$ is known wainwright2019high. Condition (d) limits the maximum sparsity and number of interventions allowed in terms of the sample size $|\Pi_u|$. Equivalently, condition (d) states that $|\Pi_u| = \Omega(s^{1/c_2})$ and $|\Pi_u| = \Omega(p^{1/c_3})$, i.e., at least a polynomial number of samples is required in both $s$ and $p$. Condition (e) states the probability of recovering the wrong support set (i.e., $\hat{\mathcal{S}}_u \neq \mathcal{S}_u$) decays exponentially quickly with $|\Pi_u|$. Consistency of recovering the support is necessary for asymptotic normality, since $\hat{\boldsymbol{\alpha}}^{SR}_u$ needs to be unbiased, which requires $\hat{\mathcal{S}}_u = \mathcal{S}_u$. Previous works wainwright2019high,zhao2006model have established that $\mathbb{P}(\hat{\mathcal{S}_u} \neq \mathcal{S}_u)$ decays exponentially quickly for the Lasso under certain regularity conditions. Next, we use normality of the horizontal regression predictions to establish asymptotic normality for non-donor units.
This section establishes asymptotic normality for non-donor units assuming that the horizontal regression predictions are asymptotically normal. This can be achieved, for example, by using the SR procedure described above. However, the result in this section does not assume the use of any particular horizontal regression method since step 1 of Synthetic Combinations is method agnostic. Rather, it is just assumed that the donor unit predictions are normal, allowing for the use of any horizontal regression method with asymptotically normal predictions---the variant of Lasso in Section (ref) being one such concrete example. Intuitively, since non-donor unit predictions are a linear combination of donor set estimates, normality of the horizontal regression is necessary to ensure asymptotic normality of the vertical regression.
We define some necessary notation for our result. Let $\Tilde{\mathbf{w}}^n = \mathbf{V}_{\mathcal{I}}^{(\Pi_n)}(\mathbf{V}_{\mathcal{I}}^{(\Pi_n)})^T\mathbf{w}^n \in \mathbb{R}^{|\mathcal{I}|}$, where $\mathbf{V}_{\mathcal{I}}^{(\Pi_n)}$ are the right singular vectors of $\mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$. That is, $\Tilde{\mathbf{w}}^n$ denotes the orthogonal projection of $\mathbf{w}^n$ onto the rowspace of $\mathbb{E}[\mathbf{Y}_{\mathcal{I}}^{(\Pi_n)}]$. Let $\Tilde{w}_u^n \in \mathbb{R}$ denote a coordinate of $\Tilde{\mathbf{w}}^n$ corresponding to a given donor unit $u \in \mathcal{I}$. Additionally, recall that $M \coloneqq \min_{u \in \mathcal{I}} |\Pi_u| $.
Theorem (ref) establishes asymptotic normality for non-donor unit and combination pair $(n,\pi)$, allowing for construction of valid confidence intervals. For (ref) to hold, $\lVert \Tilde{\mathbf{w}}_n \rVert_2$ needs to be sufficiently large (e.g., $\lVert \Tilde{\mathbf{w}}_n \rVert_2 \geq c$, for a positive constant $c$), and (ignoring $\log$ factors) if $M = \omega(r^{4}_n s^2p)$, and $|\Pi_n| = \omega(r_n^4)$. Recall from the discussion in Section (ref), $M = \omega(r^{4}_n s^2p)$, and $|\Pi_n| = \omega(r_n^4)$ are precisely the conditions required for consistency as well. As shown in Section (ref), (ref) holds when using the SR estimator for horizontal regression.
In this section, we corroborate our theoretical findings with numerical simulations in both the observational and experiment design setting. Synthetic Combinations is compared to the Lasso (i.e., running a separate horizontal regression for every unit) and two benchmark matrix completion algorithms: SoftImpute mazumder2010spectral, and IterativeSVD troyanskaya2001missing . We also try running our experiments using the nuclear-norm minimization method introduced in candes2012exact, but the method was unable to be run within a reasonable time frame (6 hours).\footnote{Code for Synthetic Combinations can be found at \url{https://github.com/aagarwal1996/synth_combo}.}
For this numerical experiment, we simulate an observation pattern commonly found in applications such as recommendation engines, where most users tend to provide ratings for combinations of goods they either particularly liked or disliked, i.e., there is confounding. The precise experimental set-up for this setting is as follows.
Experimental Set-up. We consider $N = 100$ units, and vary the number of interventions $p \in \{10,11, \ldots 15\}$. Further, let $r = 3$, and $s = rp^{3/2}$. Next, we describe how we generate the potential outcomes and observation.
Generating potential outcomes. For every $i \in [r],$ $\boldsymbol{\alpha}_i$ is generated by sampling $p^{3/2}$ non-zero coefficients at random. Every non-zero coefficient of $\boldsymbol{\alpha}_i$ is then sampled from a standard normal distribution. Denote $\mathcal{A}_r = [\boldsymbol{\alpha}_i : i \in [r]] \in \mathbb{R}^{r \times 2^p}$. Next, let $\mathcal{A}_{N - r}= \mathbf{B}\mathcal{A}_r \in \mathbb{R}^{(N - r) \times 2^p}$, where $\mathbf{B} \in \mathbb{R}^{(N - r) \times r}$ is sampled i.i.d from a standard Dirichlet distribution. Let $\mathcal{A} = [\mathcal{A}_r, \mathcal{A}_{N-r}]$; by construction, $\text{rank}(\mathcal{A}) = r$, and $\lVert \boldsymbol{\alpha}_n \rVert_0 \leq rp^{3/2} = s$ for all $n \in [N]$. Normally distributed noise $\epsilon \sim N(0,\sigma^2)$ is added to the potential outcomes $\mathbb{E}[Y_n^{(\pi)}]$, where $\sigma^2$ is chosen such that the signal-to-noise ratio is 1.0.
Observation pattern. Denote $p_{n,\pi}$ as the probability that $\mathbb{E}[Y_n^{(\pi)}]$ is observed. Let $p_{n,\pi} = |\mathbb{E}[Y_n^{(\pi)}]|/\sum_{\pi \in \Pi} |\mathbb{E}[Y_n^{(\pi)}]|$; this definition induces the missingness pattern where outcomes with larger absolute values are more likely to be seen. $|\mathcal{I}| = 2r$ donor units are chosen uniformly at random. For each donor unit $u \in \mathcal{I}$, randomly sample $|\Pi_u| = 2p^{5/2}$ combinations without replacement, where each combination $\pi$ is chosen with probability $p_{u,\pi}$. That is, each donor unit $u$ is assigned $2p^{5/2}$ treatments, where the outcome $Y_u^{(\pi)}$ is observed with probability $p_{u,\pi}$. For each non-donor unit $n$, generate $2r^4$ observation according to the same procedure. That is, $|\Pi_n| = 2r^{4}$ combinations are randomly sampled without replacement, where each combination $\pi$ is chosen with probability $p_{n,\pi}$.
Hyper-parameter choices. The hyper-parameters of the sub-procedures of Synthetic Combinations, i.e., Lasso and PCR, are tuned via $5-$fold CV. SoftImpute and IterativeSVD require that the rank $r$ of the underlying matrix to be recovered is provided as a hyper-parameter. We provide the true rank $r = 3$ to both algorithms.
Results. We measure mean squared error (MSE) averaged over 5 repetitions between the estimated potential outcome matrix and the true potential outcome matrix for each method. The results are displayed in Figure (ref) (a). Synthetic Combinations outperforms other approaches as $p$ grows. Further, the gap in performance between Synthetic Combinations and the Lasso enforces the utility of using PCR for non-donor units that do not have sufficient measurements.
Experimental Set-up. Outcomes are generated as in the observational setting. The observation pattern is generated by the experimental design mechanism described in Section (ref).
Results. We plot the MSE (averaged over 5 repetitions) for the different methods as the number of interventions is increased $p \in \{10,11,\ldots 15\}$ in Figure (ref) (b). Synthetic Combinations significantly outperforms other methods (and itself in the observational setting), which corroborates our theoretical findings that this experimental design utilizes the strengths of the estimator effectively.
This section details a real-world data experiment on recommendation systems for sets of movie ratings, comparing Synthetic Combinations to the algorithms listed above. Further, we empirically validate that our key modeling assumptions (i.e., low-rank condition of $\mathcal{A}$ and sparsity of donor unit Fourier coefficients) hold in this dataset. For all methods, hyper-parameters are chosen via $5-$fold CV. Additionally, the donor set $\mathcal{I}$ is chosen via the procedure described in Section (ref).
Data and Experimental Set-up. We use data collected in sharma2019sets which consists of user ratings of sets of movies. Specifically, users were asked to provide a rating of $1$-$5$ on a set of $5$ movies chosen at random. This resulted in a total of ratings from $854$ users over $29,516$ sets containing $12,549$ movies. $80\%$ of each user's ratings are chosen as the training set, and the other $20\%$ as the test set.
Results. We measure the RMSE averaged over 5 repetitions for all methods, and display the results in the Table below. with Synthetic Combinations outperforming all other approaches. The performance gap between Synthetic Combinations and other methods demonstrates the benefit of only performing the horizontal regression on the units with sufficient observations (i.e., the donor units), and using PCR for the non-donor units that have an insufficient number of measurements.
Key Assumptions of Synthetic Combinations hold. We also verify that our two key modeling assumptions (i.e., low-rankness and sparsity of $\mathcal{A}$) hold in this dataset. For the low-rank condition,the singular value spectrum of movies rated by all users is plotted on a log-scale in Figure (ref). The plot shows the outcomes (and hence the Fourier coefficients) are low-rank. To investigate sparsity, we examine $\hat{\boldsymbol{\alpha}}_u$ for the donor set. The MSE averaged across all donor units on the test set was 0.22, indicating that the estimated Fourier coefficient is an accurate representation of the true underlying Fourier coefficient. Further, the estimated donor unit Fourier coefficients are indeed sparse, and on average have 8.7% non-zero coefficients.
This paper introduces a causal inference framework for combinatorial interventions, a setting that is ubiquitous in practice. We propose a model that imposes both unit-specific combinatorial structure and latent similarity across units. Under this model, Synthetic Combinations, an estimation procedure, is introduced. Synthetic Combinations leverages the sparsity and low-rankness of the Fourier coefficients to efficiently estimate all $N \times 2^p$ causal parameters while implicitly allowing for unobserved confounding. Theoretically, finite-sample consistency and asymptotic normality of Synthetic Combinations is established. A novel experiment design mechanism is proposed, which ensures that the key assumptions required for the estimator to accurately recover all $N \times 2^p$ mean potential outcomes of interest hold. The empirical effectiveness of Synthetic Combinations is demonstrated through numerical simulations and a real-world case study on movie ratings. We discuss how Synthetic Combinations can be adapted to estimate counterfactuals under different permutations of items, i.e., rankings. This work suggests future directions for research such as providing an analysis of Synthetic Combinations that is agnostic to the horizontal regression algorithm used and deriving estimators that can achieve the sample complexity lower bound discussed. More broadly, we hope this work serves as a bridge between causal inference and the Fourier analysis of Boolean functions.
\if11{
We thank Alberto Abadie, Peng Ding, Giles Hooker, Devavrat Shah, Vasilis Syrgkanis, and Bin Yu for useful discussions and feedback. We also thank Austin Serif for his help in implementing Synthetic Combinations. } \fi
\bibliographystyleappendix{agsm}