Causal Inference for Multiple Treatments using Fractional Factorial Designs

TL;DR

This paper introduces a framework for designing and analyzing multi-factor experiments with fractional factorial and incomplete designs, especially useful when resources are limited or data is sparse, connecting design-based and regression methods.

Contribution

It develops a novel approach for applying fractional factorial designs to observational studies and links them with standard regression techniques.

Findings

Effective analysis of multi-factor experiments with limited data.

Application to real biomedical data demonstrating practical utility.

Insights into pesticide effects on body mass index.

Abstract

We consider the design and analysis of multi-factor experiments using fractional factorial and incomplete designs within the potential outcome framework. These designs are particularly useful when limited resources make running a full factorial design infeasible. We connect our design-based methods to standard regression methods. We further motivate the usefulness of these designs in multi-factor observational studies, where certain treatment combinations may be so rare that there are no measured outcomes in the observed data corresponding to them. Therefore, conceptualizing a hypothetical fractional factorial experiment instead of a full factorial experiment allows for appropriate analysis in those settings. We illustrate our approach using biomedical data from the 2003-2004 cycle of the National Health and Nutrition Examination Survey to examine the effects of four common pesticides on body mass index.

Tables

Table 1. Table 1: Example of a 2 3 superscript 2 3 2^{3} factorial design.

Treatment	Factor $1$	Factor $2$	Factor $3$	Outcomes
${𝒛}_1$	-1	-1	-1	$\overline{Y}\left({𝒛}_1\right)$
${𝒛}_2$	-1	-1	+1	$\overline{Y}\left({𝒛}_2\right)$
${𝒛}_3$	-1	+1	-1	$\overline{Y}\left({𝒛}_3\right)$
${𝒛}_4$	-1	+1	+1	$\overline{Y}\left({𝒛}_4\right)$
${𝒛}_5$	+1	-1	-1	$\overline{Y}\left({𝒛}_5\right)$
${𝒛}_6$	+1	-1	+1	$\overline{Y}\left({𝒛}_6\right)$
${𝒛}_7$	+1	+1	-1	$\overline{Y}\left({𝒛}_7\right)$
${𝒛}_8$	+1	+1	+1	$\overline{Y}\left({𝒛}_8\right)$
	${𝒈}_1$	${𝒈}_2$	${𝒈}_3$	$\overline{𝒀}$

Table 2. Table 2: Example of a 2 3 − 1 superscript 2 3 1 2^{3-1} factorial design.

Treatment	Factor $1$	Factor $2$	Factor $3$	Outcomes
${𝒛}_1^{\ast }$	-1	-1	+1	$\overline{Y}\left({𝒛}_1^{\ast }\right)$
${𝒛}_2^{\ast }$	-1	+1	-1	$\overline{Y}\left({𝒛}_2^{\ast }\right)$
${𝒛}_3^{\ast }$	+1	-1	-1	$\overline{Y}\left({𝒛}_3^{\ast }\right)$
${𝒛}_4^{\ast }$	+1	+1	+1	$\overline{Y}\left({𝒛}_4^{\ast }\right)$
	${𝒈}_1^{\ast }$	${𝒈}_2^{\ast }$	${𝒈}_3^{\ast }$	${\overline{𝒀}}_{\ast }$

Table 3. Table 3: Example of a 2 3 superscript 2 3 2^{3} factorial design with no observations for one treatment combination.

Treatment	$1$	$2$	$3$	Observed Outcomes
${𝒛}_1$	-1	-1	-1	${\overline{Y}}^{obs}\left({𝒛}_1\right)$
${𝒛}_2$	-1	-1	+1	${\overline{Y}}^{obs}\left({𝒛}_2\right)$
${𝒛}_3$	-1	+1	-1	${\overline{Y}}^{obs}\left({𝒛}_3\right)$
${𝒛}_4$	-1	+1	+1	${\overline{Y}}^{obs}\left({𝒛}_4\right)$
${𝒛}_5$	+1	-1	-1	${\overline{Y}}^{obs}\left({𝒛}_5\right)$
${𝒛}_6$	+1	-1	+1	${\overline{Y}}^{obs}\left({𝒛}_6\right)$
${𝒛}_7$	+1	+1	-1	?
${𝒛}_8$	+1	+1	+1	${\overline{Y}}^{obs}\left({𝒛}_8\right)$
	${𝒈}_1$	${𝒈}_2$	${𝒈}_3$	${\overline{𝒀}}^{obs}$

Table 4. Table 4: Design matrix and counts. +1 refers to exposure to the pesticide, -1 refers to no detectable exposure to the pesticide. Red treatments are used when recreating a fractional factorial design. “Original” column on the right gives the counts for each treatment for the full data set. “Trimmed” column on the right gives counts for those treatments used in the fractional factorial design, after trimming.

	Factor Levels				Number of Obs.
	beta-Hex	Hept Epox	Mirex	p,p’-DDT	Original	Trimmed
${𝒛}_1$	-+1	-+1	-+1	-+1	426	-
${𝒛}_2$	+-1	-+1	-+1	-+1	12	6
${𝒛}_3$	-+1	+-1	-+1	-+1	70	22
${𝒛}_4$	+-1	+-1	-+1	-+1	51	-
${𝒛}_5$	-+1	-+1	+-1	-+1	291	102
${𝒛}_6$	+-1	-+1	+-1	-+1	25	-
${𝒛}_7$	-+1	+-1	+-1	-+1	94	-
${𝒛}_8$	+-1	+-1	+-1	-+1	54	10
${𝒛}_9$	-+1	-+1	-+1	+-1	21	8
${𝒛}_{10}$	+-1	-+1	-+1	+-1	1	-
${𝒛}_{11}$	-+1	+-1	-+1	+-1	19	-
${𝒛}_{12}$	+-1	+-1	-+1	+-1	19	10
${𝒛}_{13}$	-+1	-+1	+-1	+-1	42	-
${𝒛}_{14}$	+-1	-+1	+-1	+-1	19	3
${𝒛}_{15}$	-+1	+-1	+-1	+-1	37	8
${𝒛}_{16}$	+-1	+-1	+-1	+-1	78	-
	${𝒈}_1$	${𝒈}_2$	${𝒈}_3$	${𝒈}_4$	1259	169

Equations

$$\tau (k)=\frac{1}{2^{K-1}}{\mathbf{g}}_k^T\bar{\mathbf{Y}}\text{and}\tau (k,k^{\prime })=\frac{1}{2^{K-1}}{\mathbf{g}}_{k,k^{\prime }}^T\bar{\mathbf{Y}},$$

$${\tau }_i(k)=\frac{1}{2^{K-1}}{\mathbf{g}}_k^T{\mathbf{Y}}_i\text{and}{\tau }_i(k,k^{\prime })=\frac{1}{2^{K-1}}{\mathbf{g}}_{k,k^{\prime }}^T{\mathbf{Y}}_i,$$

$$\mathbf{G}=[{\mathbf{g}}_0,{\mathbf{g}}_1,\dots ,{\mathbf{g}}_K,{\mathbf{g}}_{1,2},\dots ,{\mathbf{g}}_{1,\dots ,K}{]}^T$$

$$\mathbf{\tau }=[2\tau (0),\tau (1),\dots ,\tau (K),\tau (1,2),\dots ,\tau (K-1,K),\dots ,\tau (1,2,\dots ,K){]}^T.$$

$$2^{-(K-1)}\mathbf{G}\mathbf{Y}=\mathbf{\tau }.$$

$$P(\mathbf{W}=\mathbf{w})=\prod _{j=1}^J{n}_j!/n!$$

$${\bar{Y}}^{obs}({\mathbf{z}}_j)=\frac{1}{n_j}\sum _{i=1}^n{W}_i({\mathbf{z}}_j)Y_i({\mathbf{z}}_j)=\frac{1}{n_j}\sum _{i:W_i({\mathbf{z}}_j)=1}^n{Y}_i^{obs}.$$

$$\hat{\tau }(k)=\frac{1}{2^{K-1}}{\mathbf{g}}_k^T{\bar{\mathbf{Y}}}^{obs}\text{and}\hat{\tau }(k,k^{\prime })=\frac{1}{2^{K-1}}{\mathbf{g}}_{k,k^{\prime }}^T{\bar{\mathbf{Y}}}^{obs},$$

$$\text{Var}\left(\hat{\tau }(k)\right)=\frac{1}{2^{2(K-1)}}\sum _{j=1}^J\frac{1}{n_j}{S}^2({\mathbf{z}}_j)-\frac{1}{n}{S}_k^2,$$

$$S^2({\mathbf{z}}_j)=\frac{1}{n-1}\sum _{i=1}^n{\left(Y_i({\mathbf{z}}_j)-\bar{Y}({\mathbf{z}}_j)\right)}^2\text{and}{S}_k^2=\frac{1}{n-1}\sum _{i=1}^n{\left({\tau }_i(k)-\tau (k)\right)}^2.$$

$$\text{Cov}\left(\hat{\tau }(k),\hat{\tau }(k^{\prime })\right)=\frac{1}{2^{2(K-1)}}\left[\sum _{j:g_{kj}=g_{k^{\prime }j}}\frac{1}{n_j}{S}^2({\mathbf{z}}_j)-\sum _{j:g_{kj}\ne g_{k^{\prime }j}}\frac{1}{n_j}{S}^2({\mathbf{z}}_j)\right]-\frac{1}{n}{S}_{k,k^{\prime }}^2,$$

$$P({\mathbf{W}}^{\ast }={\mathbf{w}}^{\ast })=\prod _{j=1}^{J^{\prime }}{n}_j^{\ast }!/n!$$

$${\hat{\tau }}^{\ast }(k)=\frac{1}{2^{K-p-1}}{\mathbf{g}}_k^{\ast T}{\bar{\mathbf{Y}}}_{\ast }^{obs}.$$

$${\hat{\tau }}^{\ast }(k,k^{\prime })=\frac{1}{2^{K-p-1}}{\mathbf{g}}_{k,k^{\prime }}^{\ast T}{\bar{\mathbf{Y}}}_{\ast }^{obs}.$$

$$E\left[{\hat{\tau }}^{\ast }(k)\right]=\frac{1}{2^{K-p-1}}{\mathbf{g}}_k^{\ast T}{\bar{\mathbf{Y}}}_{\ast }=\sum _{j\in \mathcal{A}[k]}(-1{)}^{A_{k,j}}\tau (j)=\tau (k)+\sum _{j\in \mathcal{A}[k]\backslash \{k\}}(-1{)}^{A_{k,j}}\tau (j).$$

$${\tilde{S}}_k^2=\frac{1}{n-1}\sum _{i=1}^n{\left({\tilde{\tau }}_i(k)-\tilde{\tau }(k)\right)}^2,$$

$$\text{Var}\left({\hat{\tau }}^{\ast }(k)\right)$$

$$\text{Var}\left(\hat{\tau }(k)\right)=\frac{1}{2^{2(K-1)}}\frac{J^2}{n}{S}^2=\frac{4}{n}{S}^2$$

$$\text{Var}\left({\hat{\tau }}^{\ast }(k)\right)=\frac{1}{2^{2(K-p-1)}}\frac{(J^{\prime }{)}^2}{n}{S}^2=\frac{4}{n}{S}^2.$$

$${\tilde{S}}_{k,k^{\prime }}^2$$

$$\text{Cov}\left({\hat{\tau }}^{\ast }(k),{\hat{\tau }}^{\ast }(k^{\prime })\right)=\frac{1}{2^{2(K-p-1)}}\left[\sum _{j:g_{kj}^{\ast }=g_{k^{\prime }j}^{\ast }}\frac{1}{n_j^{\ast }}{S}^2({\mathbf{z}}_j^{\ast })-\sum _{j:g_{kj}^{\ast }\ne g_{k^{\prime }j}^{\ast }}\frac{1}{n_j^{\ast }}{S}^2({\mathbf{z}}_j^{\ast })\right]-\frac{1}{n}{\tilde{S}}_{k,k^{\prime }}^2.$$

$$\widehat{\text{Var}}\left({\hat{\tau }}^{\ast }(k)\right)$$

$$E\left[\widehat{\text{Var}}\left({\hat{\tau }}^{\ast }(k)\right)\right]=\frac{1}{2^{2(K-p-1)}}\sum _{j=1}^{J^{\prime }}\frac{1}{n_j^{\ast }}{S}^2({\mathbf{z}}_j^{\ast }).$$

$$2^{-(K-p-1)}\underset{{\mathbf{G}}^{\ast }}{\underbrace{\left[\begin{array}{c}{\mathbf{g}}_0^{\ast T}\\ {\mathbf{g}}_1^{\ast T}\\ {\mathbf{g}}_2^{\ast T}\\ {\mathbf{g}}_3^{\ast T}\end{array}\right]}}\underset{{\bar{\mathbf{Y}}}^{\ast }}{\underbrace{\left[\begin{array}{c}\bar{Y}({\mathbf{z}}_1^{\ast })\\ \bar{Y}({\mathbf{z}}_2^{\ast })\\ \bar{Y}({\mathbf{z}}_3^{\ast })\\ \bar{Y}({\mathbf{z}}_4^{\ast })\end{array}\right]}}=\underset{\tilde{\mathbf{\tau }}}{\underbrace{\left[\begin{array}{c}2\tilde{\tau }(0)\\ \tilde{\tau }(1)\\ \tilde{\tau }(2)\\ \tilde{\tau }(3)\end{array}\right]}}.$$

$$\sqrt{n}\left({\hat{\mathbf{\tau }}}^{\ast }-\tilde{\mathbf{\tau }}\right)\stackrel{d}{\to }\text{N}(0,\mathbf{V}).$$

$$\{\mathbf{\mu }:{\left({\hat{\mathbf{\tau }}}^{\ast }-\mathbf{\mu }\right)}^T\hat{\mathbf{V}}\left({\hat{\mathbf{\tau }}}^{\ast }-\mathbf{\mu }\right)\le q_{J^{\prime },1-\alpha }\},$$

$$\dot{\tau }(1)$$

$$=\frac{1}{3}\left[\tau (1|F_2=-1,F_3=-1)+\tau (1|F_2=-1,F_3=+1)+\tau (1|F_2=+1,F_3=+1)\right],$$

$$\dot{\tau }(1)$$

$$=\tau (1)+\frac{1}{3}\left[-\tau (13)+\tau (23)+\tau (123)\right].$$

$$\dot{\tau }(1)$$

Full text

Keywords: potential outcomes; interactions; joint effects; multiple treatments; Neymanian inference; observational studies

1 Introduction

What is the effect of exposing you to pesticide A, compared to no exposure, on your health? What about exposing you to pesticide B? Or exposing you to both pesticides at the same time? To answer these questions, we need to estimate the effects, including interactions, of multiple treatments. A randomized factorial experiment uses a random assignment of all possible treatment combinations to units to estimate these different effects. There is much interest in assessing the effects of multiple treatments, as reflected by the recent growth in the literature regarding the use of factorial designs in causal inference (e.g., Branson et al.,, 2016; Dasgupta et al.,, 2015; Dong,, 2015; Egami and Imai,, 2019; Espinosa et al.,, 2016; Lu, 2016a, ; Lu, 2016b, ; Lu and Deng,, 2017; Mukerjee et al.,, 2018; Zhao et al.,, 2018).

However, a full factorial design, in which all possible combinations of treatments are randomized, may not be possible due to constraints on resources such as units or cost. Instead, a fractional factorial or an incomplete factorial design, which uses a subset of all possible treatment combinations, may be more practical. We layout a Neyman-style randomization-based framework for the design and analysis of these designs using potential outcomes. Consideration of observational studies with multiple factors, in which certain treatment combinations may be missing from the data set, further motivates the usefulness of these designs.

With only a subset of all treatment combinations, saturated linear regression models have underlying assumptions resulting in estimators that are not always transparent, especially with respect to the implicit imputation of the missing potential outcomes. We discuss regression estimates in such setting in Section 3.2 and assess how they connect to design-based estimates. Analyses using regression models in observational studies are even trickier, as the use of regression without a careful design phase, in which one tries to uncover or approximate some underlying randomized experiment, can lead to incorrect conclusions (e.g., see Rubin,, 2008). Despite this difficulty, regression models with interaction terms are commonly used to estimate the effects of multiple treatments, in particular in observational studies (Bobb et al.,, 2015; Oulhote et al.,, 2017; Patel et al.,, 2010; Valeri et al.,, 2017). For instance, Bobb et al., (2015) considers a Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures. Therefore, exploration of design-based methods for multi-factor observational studies is an important missing piece in the literature.

The paper proceeds as follows: Section 2 reviews full factorial designs within the potential outcomes framework described in Dasgupta et al., (2015). Section 3 discusses extensions of this framework to fractional and incomplete factorial designs, expanding upon current inference results for variance and variance estimation. Section 4 examines how to embed an observational study into one of these experimental designs and illustrates our method and the challenges when working with observational data with multiple treatments with an application examining the effects of four pesticides on body mass index (BMI) using data from the National Health and Nutrition Examination Survey (NHANES), which is conducted through the Centers for Disease Control and Prevention (CDC). Section 5 concludes.

Throughout, we will use a running example of an experiment examining the impact of three pesticides (which we call pesticides A, B, and C) on weight gain in zebrafish. Our data example provides an illustration of the way we foresee these designs being useful when experimentation is not possible, such as when discussing the impact of pesticides on humans.

2 Set up

2.1 Notation and factorial estimands

We closely follow the potential outcomes (Splawa-Neyman,, 1990; Rubin,, 1974) framework for $2^{K}$ factorial designs in Dasgupta et al., (2015). We focus on designs with multiple two-level factors assigned in combination; that is, there are a number of distinct treatments (e.g., medications), each having two levels (e.g., receiving a certain medication or not). Let there be $K$ two-level factors, which creates $2^{K}=J$ total possible treatment combinations. The $j^{th}$ treatment combination is denoted $\bm{z}_{j}$. See Table 1 for an example with $K=3$, which illustrates the notation. For consistency, we list treatment combinations in lexicographic order, the standard ordering. Let $z_{j,k}\in\{-1,+1\}$ be the level of the $k^{th}$ factor in the $j^{th}$ treatment combination, so $\bm{z}_{j}=(z_{j,1},...,z_{j,K})$. Let there be $n$ units in the sample. We assume the Stable Unit Treatment Value Assumption (SUTVA): there is no interference and no different forms of each treatment level (Rubin,, 1980). Then the potential outcome for unit $i$ under treatment combination $\bm{z}_{j}$ can be written as $Y_{i}(\bm{z}_{j})$ and the sample average potential outcome under treatment $\bm{z}_{j}$ as $\bar{Y}(\bm{z}_{j})=\sum_{i=1}^{n}Y_{i}(\bm{z}_{j})/n$. $\bar{\bm{Y}}=(\bar{Y}(\bm{z}_{1}),\bar{Y}(\bm{z}_{2}),\ldots,\bar{Y}(\bm{z}_{J}))$ is the vector of mean potential outcomes for all $2^{K}$ treatment combinations.

Using the framework from Dasgupta et al., (2015), denote the contrast column $j$ in the design matrix by $\bm{g}_{j}$, as illustrated in Table 1. Following that paper, we can also define the contrast vector for the two-factor interaction between factors $k$ and $k^{\prime}$ as $\bm{g}_{k,k^{\prime}}=\bm{g}_{k}\circ\bm{g}_{k^{\prime}}$, where $\circ$ indicates the Hadamard (element-wise) product. Similarly, the contrast vector for three-factor interactions is $\bm{g}_{k,h,i}=\bm{g}_{k}\circ\bm{g}_{h}\circ\bm{g}_{i}$, and all higher order interaction contrast vectors can be calculated analogously.

Continuing to follow Dasgupta et al., (2015), define the finite population main causal effect for factor $k$ and the finite population interaction effect for factors $k$ and $k^{\prime}$ as

[TABLE]

respectively, where by finite population we mean that we are only interested in inference for the units we have in the experiment. Higher-level interaction terms are defined analogously. We can similarly define the individual-level effects as

[TABLE]

where $\bm{Y}_{i}=(Y_{i}(\bm{z}_{1}),\ldots,Y_{i}(\bm{z}_{J}))^{T}$. We also define the average potential outcome across treatments as $\tau(0)=2^{-K}\bm{g}^{T}_{0}\bar{\bm{Y}}$, where $\bm{g}_{0}$ is a vector of length $2^{K}$ of all $+1$’s.

Consider our example of estimating the effects of pesticides A, B, and C on weight gain in zebrafish. We take the three pesticides to be three factors, each with two levels: applying the pesticide (+1) or not (-1). In Table 1, zebrafish assigned to $\bm{z}_{6}$ would be those that get pesticides A and C, but not B. The main effect of pesticide A ($\tau(1)$) is the average effect of receiving pesticide A, compared to not receiving it, uniformly averaged over all combinations of receiving or not receiving pesticides B and C. The two-factor interaction between pesticides A and B, would be the difference in the effect of receiving pesticide A, compared to not receiving it, when pesticide B is applied versus when it is not, averaging uniformly over receiving or not receiving pesticide C.

Consider a full factorial model matrix, $\bm{G}$, that includes the mean and interactions and whose rows are comprised of $\bm{g}_{k}^{T}$, where $k\in\{0;1;...;K;(1,2);...;(1,2,\dots,K)\}$. Based on the definition of our estimands in Section 2.1, the matrix $\bm{G}$ relates the mean potential outcomes and the factorial effects as follows: Let $\bm{G}$ be a matrix with the $\bm{g}$ vectors as rows,

[TABLE]

and $\bm{\tau}$ be the column vector of all factorial effects,

[TABLE]

Then,

[TABLE]

Due to orthogonality, $\bm{G}^{-1}=\frac{1}{2^{K}}\bm{G}^{T}$, as argued by Dasgupta et al., (2015) in the context of imputing potential outcomes under Fisher’s sharp null hypothesis. The mean potential outcomes can be rewritten in terms of the factorial effects as ${\bar{\bm{Y}}=2^{-1}\bm{G}^{T}\bm{\tau}}$. Let the $j^{th}$ row of $\bm{G}^{T}$ be denoted by $\bm{h}_{j}$. The first entry of $\bm{h}_{j}$ is $+1$, corresponding to the mean, the next $K$ entries are equal to the entries of $\bm{z}_{j}$, and the remaining entries correspond to interactions, with the order given by the order of the rows of $\bm{G}$. In the $2^{3}$ factorial design shown in Table 1, $\bm{h}_{1}=(+1,-1,-1,-1,+1,+1,+1,-1)$. That is, it is the entries of the first row of Table 1, prepended with an entry for $2\tau(0)$ and appended with entries for the interactions. We have $\bar{Y}(\bm{z}_{j})=2^{-1}\bm{h}_{j}\bm{\tau}$. We will use this representation in Section 3.2 when considering causal estimands in settings lacking units for some treatment combinations.

2.2 Estimation and inference in a $2^{K}$ factorial design

We now review inference for $2^{K}$ factorial designs. Let $W_{i}(\bm{z}_{j})=1$ if unit $i$ received treatment $\bm{z}_{j}$ and $W_{i}(\bm{z}_{j})=0$ otherwise. In a full $2^{K}$ factorial design, under complete randomization, all $2^{K}$ factors are assigned to units completely at random (possibly unbalanced). That is, there is a fixed number, $n_{j}\geq 2$, of units randomly assigned to treatment combination $\bm{z}_{j}$. Let $\bm{W}$ be the $n\times 2^{K}$ vector of $W_{i}(\bm{z}_{j})$ for all $i$ and $\bm{z}_{j}$. Let $\bm{w}$ be some $n\times 2^{K}$ matrix with entries in the set $\{0,1\}$. This design imposes

[TABLE]

if $\sum_{j=1}^{J}w_{i}(\bm{z}_{j})=1$ for all $i$ and $\sum_{i=1}^{n}w_{i}(\bm{z}_{j})=n_{j}$ for all $\bm{z}_{j}$, and $P(\bm{W}=\bm{w})=0$ otherwise. The observed potential outcome for unit $i$ is $Y_{i}^{obs}=\sum_{j=1}^{J}W_{i}(\bm{z}_{j})Y_{i}(\bm{z}_{j})$. An unbiased estimator of the average potential outcome for treatment $\bm{z}_{j}$ is

[TABLE]

Denote $\bm{\bar{Y}}^{obs}=(\bar{Y}^{obs}(\bm{z}_{1}),\bar{Y}^{obs}(\bm{z}_{2}),\ldots,\bar{Y}^{obs}(\bm{z}_{J}))$ as the vector of observed mean potential outcomes for all $2^{K}$ treatment combinations.

Unbiased estimators for $\tau(k)$ and $\tau(k,k^{\prime})$ are defined as

[TABLE]

respectively. Estimators for higher-level interaction terms and $\tau(0)$ are defined analogously.

The variance of the factorial effect estimators under an unbalanced design is given in Lu, 2016b (extended from balanced case in Dasgupta et al.,, 2015), as follows:

[TABLE]

where

[TABLE]

An expression for the covariance between two factorial effect estimators for unbalanced factorial designs is given by Lu, 2016b (extended from Dasgupta et al.,, 2015) as

[TABLE]

where $S^{2}_{k,k^{\prime}}=\sum_{i=1}^{n}\left(\tau_{i}(k)-\tau(k)\right)\left(\tau_{i}(k^{\prime})-\tau(k^{\prime})\right)/(n-1)$.

Further discussion of Neymanian and Fisherian inference for full factorial designs is reviewed in Supplementary Material A.

3 Alternate designs to full factorial: Fractional and incomplete

3.1 Fractional factorial designs

There are situations in which a full $2^{K}$ factorial design experiment cannot be conducted or is suboptimal. Limited resources may mean there are an insufficient number of units to randomly assign to each of the $J$ ($=2^{K}$) treatment groups. Or the full factorial design may not be the most efficient allocation of resources; for example, if the experimenter believes that higher-order interactions are negligible (Wu and Hamada,, 2000). Instead, an experimenter might implement a $2^{K-p}=J^{\prime}$ fractional factorial design in which only $J^{\prime}$ of the total $J$ treatment combinations are used. Here we will give a brief overview of this design, but we recommend Wu and Hamada, (2000) for a more detailed review. We follow notation in Dasgupta et al., (2015) and Dasgupta and Rubin, (2015).

To create a $2^{K-p}$ design, we can write out a full factorial design for the first $K-p$ factors, and fill in the $p$ columns for the remaining factors using multiplicative combinations of subsets of the previous contrast columns. We use a generator to choose the factors whose treatment levels we multiply together to get the treatment levels for the other factors (Wu and Hamada,, 2000). For example, to create a $2^{3-1}$ design using the generator “$3=12$”, we write out a $2^{2}$ design for factors 1 and 2, as in column 2 and 3 of Table 2. We generate the third column, corresponding to treatment levels for factor 3, by multiplying together the contrast vectors for factors 1 and 2, as shown in Table 2. Which two factors are used initially is irrelevant in this case because of symmetry. Because $3=12$, we also have $1=23$ and $2=13$. However, this may not hold in general and one should choose which factors to use based on the final aliasing structure, which tells us which effects are confounded.

Under this design, the contrast columns, $\bm{g}_{k}^{*}$, are now shortened versions or subsets of the contrast columns of a full $2^{3}$ factorial design, $\bm{g}_{k}$. The treatment combinations in the fractional design, $\bm{z}^{*}_{j}$, are a subset of the full set of treatment combinations, $\bm{z}_{j}$, for the $2^{3}$ design. Referring to Table 2, $\bm{g}_{3}^{*}=\bm{g}_{1}^{*}\circ\bm{g}_{2}^{*}$. The generator $3=12$ indicates that the main effect of factor 3 is aliased with the two-factor interaction 12, which means that we cannot distinguish these effects – they are confounded and combined in our estimators. The full aliasing in this design is as follows: $I=123$, $3=12$, $2=13$, $1=23$, where $I$ corresponds to a vector of all +1’s ($\bm{g}_{0}$). The relation $I=123$ is called the defining relation, which characterizes the aliasing and how to generate the rest of the columns. The main effects, as defined in Section 2.1, are aliased with the two-factor interactions. If the two-factor interactions are negligible, the 23-1 design is a parsimonious design to estimate the main effects and we can construct unbiased estimators for the main effects (see Section 3.1.1). The resulting design is also orthogonal (i.e., all pairs of columns are orthogonal) and balanced in the sense that each $\bm{g}^{*}_{k}$ has an equal number +1’s and -1’s (Wu and Hamada,, 2000). These properties simplify the aliasing structure.

We typically choose the generator or defining relation based on the maximum resolution criterion for the design. Resolution is defined as the word length (i.e., the number of factors) in the shortest word of the defining relation (see Wu and Hamada,, 2000, for more details) and indicates the aliasing structure and which levels of effects the main effects are aliased. In our example of a $2^{3-1}$ design, we only have one defining relation, $I=123$, the word is length 3, and main effects are aliased with two-factor interactions. Aliasing some main effects with other main effects is an alternate aliasing structure. The effect hierarchy principle assumes that lower-order effects have larger magnitude than higher-order interaction effects (Wu and Hamada,, 2000). Therefore, one generally chooses a fractional design where main effects and some other lower-order interaction effects are aliased with higher-order terms, i.e., the main effects and two-factor interactions are clear. This assumption goes along with the assumption of effect sparsity (the number of important effects is small), as a justification of the use of fractional designs over the full factorial (Wu and Hamada,, 2000). It is possible that there may be multiple defining relations that lead to the same resolution. In addition to the maximum resolution criterion, it is standard to consider minimum aberration, which says we want the design with the smallest number of short words in its defining relations.

3.1.1 Estimators

In this section we review estimators for the fractional design, which are similar to the full factorial case, and follow the framework laid out in Dasgupta et al., (2015) and Dasgupta and Rubin, (2015). We assume units are assigned using complete randomization to the $J^{\prime}$ treatments. In other words, if $\bm{W}^{*}$ is the $n\times 2^{K-p}$ matrix of $W_{i}(\bm{z}_{j}^{*})$ for all $i$ and $\bm{z}_{j}$, and $\bm{w}^{*}$ is some $n\times 2^{K-p}$ matrix with entries in the set $\{0,1\}$, then

[TABLE]

if $\sum_{j=1}^{J^{\prime}}w_{i}(\bm{z}_{j}^{*})=1$ for all $i$ and $\sum_{i=1}^{n}w_{i}(\bm{z}_{j}^{*})=n_{j}^{*}$ for all $\bm{z}_{j}^{*}$, and $P(\bm{W}^{*}=\bm{w}^{*})=0$ otherwise. Randomization is the same as in the full factorial case, but with fewer treatment groups. All of the following inference results are derived under complete randomization.

The estimator for $\tau(k)$ is defined as

[TABLE]

The estimator for $\tau(k,k^{\prime})$ is defined as

[TABLE]

These estimators are no longer unbiased. Let $\mathcal{A}[k]$ be the set of all effects aliased with factorial effect $k$ as well as $k$ itself. The number of factors aliased with factor $k$ is $2^{p}-1$ (Montgomery,, 2017, Chapter 8). So, $\mathcal{A}[k]$ has $2^{p}$ elements. Factor $k$ may be aliased with the negative of a main effect or interaction. Let $A_{k,j}$ be the indicator for whether factor $j$ is negatively aliased with factor $k$ ($A_{k,j}=1$) or positively aliased ($A_{k,j}=0$).

Result 1*.*

The expectation of the factorial estimator is

[TABLE]

Result 1 can be shown using the orthogonality of the $\bm{g}_{k}$ vectors. Hence, by aliasing these effects, they get combined in our estimator.

Returning to our running example, let us say we used the design in Table 2. In that design, the main effect of pesticide A is aliased with the interaction between pesticide B and C. The expectation of the main effect for pesticide A would be $E\left[\widehat{\tau}^{*}(1)\right]=\tau(1)+\tau(23)$. If the the interaction effect between pesticides B and C is small, our estimator for the main effect of pesticide A will be close to unbiased. Otherwise, our estimate of the effect of pesticide A will be biased by the size of the interaction. As an example, if the main effect of pesticide A is $\tau(1)=-2$ and the interaction between pesticides B and C is $\tau(23)=2$, the expectation of our estimator will be 0 even though the true average effect of pesticide A is nonzero.

More generally, the estimator for $\tau(k)$ is unbiased if the effects aliased with factorial effect $k$ are zero, and will be close to unbiased if the aliased effects are negligible, which may be justified by the effect hierarchy principle. That is, we expect higher-order effects to be smaller in magnitude. When aliasing occurs such that main effects are aliased with low-order interaction terms, such as two-factor interactions, this assumption may be unrealistic. However, if we have a large number of factors and are only aliasing main effects with higher-order effects, this may be reasonable. For example, in a $2^{6-1}$ fractional design, main effects can be aliased with five-factor interactions and two-factor interactions can be aliased with four-factor interactions, where both higher-order interactions may be assumed to be small.

We extend the variance and variance estimator expressions from Dasgupta et al., (2015) (see Section 2) to this setting. Recall $J^{\prime}=2^{K-p}$ and define $\tilde{\tau}_{i}(k)=\tau_{i}(k)+\sum_{j\in\mathcal{A}[k]\backslash\{k\}}(-1)^{A_{k,j}}\tau_{i}(j)$ and $\tilde{\tau}(k)=\frac{1}{N}\sum_{i=1}^{N}\tilde{\tau}_{i}(k)=E[\widehat{\tau}^{*}(k)]$. Further define the analog of $S^{2}_{k}$,

[TABLE]

which is the variation in our newly defined aliased effects, $\tilde{\tau}_{i}(k)$. Let $n^{*}_{j}$ be the number of units assigned to treatment combination $\bm{z}^{*}_{j}$.

Result 2*.*

The variance of the estimator $\hat{\tau}^{*}(k)$ is

[TABLE]

We can compare the variance under a factorial design to one under a fractional factorial design easily in a simplified setting. Take sample size to be $n$ for both designs, and both designs to be balanced ($n_{j}=n/J$, $n_{j}^{*}=n/J^{\prime}$). Further, assume we have additive effects (no effect heterogeneity), which implies $S^{2}(\bm{z}_{j})=S^{2}$ for all $\bm{z}_{j}$ and $S_{k}^{2}=\tilde{S}_{k}^{2}=0$. In this setting, the variance of the estimator for $\widehat{\tau}(k)$ for the full factorial design will be

[TABLE]

and for the fractional factorial design will be

[TABLE]

The variance of the estimator under the full factorial design will be equal to that under a fractional factorial design. As the fractional factorial design has a bias when the interaction effects being aliased are nonzero, the full factorial is preferred when possible. The general case without additive effects and a balanced design is more complicated, especially if the treatments included in the fractional factorial design happen to be those that produce less variable outcomes. However, this would be difficult to predict before running the experiment, and therefore a full factorial is preferred in general. Unbalanced designs will similarly change this result, potentially increasing or decreasing the weight on more variable groups in the two designs. A full factorial design may not always be feasible or practical, in which case inducing some (hopefully negligible) bias by using a fractional factorial design may be acceptable.

We can also obtain the covariance between two fractional factorial effect estimators. We define an altered version of $S^{2}_{k,k^{\prime}}$:

[TABLE]

Result 3*.*

The covariance of $\hat{\tau}^{*}(k)$ and $\hat{\tau}^{*}(k^{\prime})$ is

[TABLE]

These variance and covariance expressions result from an application of Theorem 3 in Li and Ding, (2017), as explained in Appendix D.1. The variance expressions for fractional factorial designs are similar to the full factorial case, but are defined over aliased effects. In particular, the treatment effect heterogeneity term has the interpretation of the heterogeneity of a linear combination of aliased effects.

3.1.2 Statistical inference

It is straightforward to extend the analysis for factorial designs in Dasgupta et al., (2015) and reviewed in Section 2 to fractional factorial designs. We can use a conservative Neyman-style variance estimator, similar to the full factorial setting (see Supplementary Material A):

[TABLE]

Result 4*.*

The expectation of the fractional factorial variance estimator is

[TABLE]

Thus $\widehat{\text{Var}}\left(\hat{\tau}^{*}(k)\right)$ is a conservative estimator and unbiased if and only if $\tilde{S}_{k}^{2}=0$, which would occur if the $\tilde{\tau}_{i}(k)$ are constant. This condition holds when all effects aliased with factor $k$ are constant additive effects.

We can build confidence regions and confidence intervals analogously to what was reviewed in Section A. To be more rigorous in this section, we show how to use the results of Li and Ding, (2017) to build asymptotic Wald-type confidence regions. There is a new matrix $\bm{G}^{*}$ that relates the treatment combinations in the fractional experiment to the unique treatment effect estimands and estimators defined for this experiment. For the $2^{3-1}$ example, recalling the aliasing structure which gives $\bm{g}^{*}_{0}=\bm{g}^{*}_{123}$ and $\tilde{\tau}(2)=\tilde{\tau}(13)$, we have

[TABLE]

Define $\hat{\bm{\tau}}^{*}$ as the $2^{K-p}$ vector of unique treatment effect estimators that corresponds to estimands in $\tilde{\bm{\tau}}$.

Result 5*.*

Assume that all $S^{2}(\bm{z}^{*}_{j})$ and $\tilde{S}_{k,k^{\prime}}^{2}$ have limiting values, that the $n^{*}_{j}/n$ have positive limiting values, and $\text{max}_{1\leq j\leq J^{\prime}}\text{max}_{1\leq i\leq n}\left(Y_{i}(z^{*}_{j})-\bar{Y}(z^{*}_{j})\right)^{2}/n\to 0$. According to Theorem 5 of Li and Ding, (2017), $n\text{var}(\hat{\bm{\tau}}^{*})$ has a limiting value, which, following their notation, we denote by $\bm{V}$ and

[TABLE]

Result 6*.*

Further, let $\hat{\bm{V}}=\sum_{j=1}^{J^{\prime}}n^{*}_{j}\bm{G}^{*}_{j}s^{2}(\bm{z}^{*}_{j})\bm{G}^{*T}_{j}$, where $\bm{G}^{*}_{j}$ is the $j$th column of $\bm{G}^{*}$. In addition to the assumptions required for Result 5, we require that the limit of $\sum_{j=1}^{J^{\prime}}n^{*}_{j}\bm{G}^{*}_{j}S^{2}(\bm{z}^{*}_{j})\bm{G}^{*T}_{j}$ is nonsingular. Under Proposition 3 of Li and Ding, (2017) and following their notation, the Wald-type confidence region,

[TABLE]

where $q_{J^{\prime},1-\alpha}$ corresponds to the $1-\alpha$ quantile of the $\chi^{2}_{J^{\prime}}$ distribution, has at least $1-\alpha$ asymptotic coverage.

We assume that the number of treatments is constant as $n\to\infty$, but Li and Ding, (2017) discuss the case where the number of treatment combinations grows as well.

As in the full factorial setting (see Lu, 2016b , Zhao and Ding, (2021), and Supplementary Material A), linear regression on the full vector of observed outcomes, using a saturated model with factor levels coded as -1 and +1, yields the same point estimates (divided by 2) and the HC2 variance estimator yields the same variance estimate (divided by 4) as the Neymanian estimates. For proof, see Appendix E. Zhao and Ding, (2021) show that weighted-least squares can also be used to obtain similar results in an unsaturated regression in the full factorial setting.

3.2 Incomplete factorial designs

In this section, we discuss an alternative experimental design called the incomplete factorial design (Byar et al.,, 1993), which we consider (novelly, to our knowledge) from the potential outcome and design-based perspective. This design uses a subset of data from a full factorial design but a different subset than the fractional design. In particular, the estimators we discuss will not use all of the available data. We can therefore consider each estimator to be associated with a particular “design” that corresponds to an experiment that randomizes units only to treatment combinations used in the estimator. Different estimators may give non-zero weight to different treatment combinations, and so the designs may be estimand-specific. Due to this inherent linking between the design and the estimators, we discuss the design and corresponding estimators together in the next section.

3.2.1 Design and estimators

Consider if we ran an experiment with $K$ binary treatments where, whether due to human error or lack of resources, some of treatment combinations from the full $2^{K}$ factorial design were not included in the experiment. We may have a factorial structure as in Table 3, but no outcome measurements for treatment $\bm{z}_{7}$. How should we analyze this?

Focus on estimation of $\tau(1)$ first. We could reduce the data to a fractional factorial design, which requires removing multiple treatment groups. If we recreate a fractional factorial design that aliases the main effects with the two-way interactions, we estimate $\tau(1)+\tau(23)$, which involves using outcome data from units assigned to treatment combinations $\bm{z}_{2}$, $\bm{z}_{3}$, $\bm{z}_{5}$, and $\bm{z}_{8}$, but not from units assigned to treatment combinations $\bm{z}_{1}$, $\bm{z}_{4}$, and $\bm{z}_{6}$.

Instead, we might consider building an estimator using all of the treatment combinations except for $\bm{z}_{3}$, which has the same levels for factors 2 and 3 but the opposite level for factor 1 as combination $\bm{z}_{7}$. Thus, in some sense, removing $\bm{z}_{3}$ “balances” the remaining treatment combinations and is essentially the “naïve estimator” discussed in Byar et al., (1993). This strategy creates a different hypothetical experimental design with a different aliasing structure. In this case, we would be estimating

[TABLE]

where we let $F_{k}$ denote the level of the $k^{th}$ factor so that ${\tau(k|F_{j}=x,F_{i}=y)}$ is the main effect of factor $k$ conditional on level $x$ of factor $j$ and level $y$ of factor $i$, as in Dasgupta et al., (2015). In words, we estimate the average of the conditional effects of factor 1 given the combinations (-1,-1), (-1,+1), and (+1,+1) for factors 2 and 3.

To find the aliasing structure, we can refer to the matrix $\bm{G}^{T}$ in Section 2.1 to rewrite the estimand above as

[TABLE]

We have partially aliased the main effect for factor 1 with the two-factor interactions for factors 1 and 3 and factors 2 and 3, as well as the three-factor interaction, all divided by three. This is partial aliasing because the factors are neither fully aliased nor completely clear of each other (Wu and Hamada,, 2000, Chapter 7). Return to our pesticide example. If $\tau(13)=2$, $\tau(23)=2.5$ and $\tau(123)=1$, the main effect of pesticide A will be positively biased by $1/2$. Whether this is preferable to aliasing the main effect with just the two-factor interaction between factors 2 and 3, as in the fractional design, depends on context. If we know that the three-factor interaction is negligible, (or at least smaller than the two-factor interactions, as in the numerical example just given) this new estimator may have lower bias, as we are dividing both two-factor interactions by 3. However, even if we had knowledge that the two-factor interactions were of the same sign, we would not know the direction of the bias in this case without knowing the relative magnitudes of the two-factor interactions.

When estimating the main effect for factor 2, we would naturally approximate a different design. Using the same logic as before, we would remove treatment combination $\bm{z}_{5}$.

As an alternative design, we may alias our main effects with the highest-order interaction possible, allowing for a different design for each main effect estimator. So, when estimating the main effect of factor 1, we would use a fractional factorial design that aliases factor 1 with the three-factor interaction. This leads to a design using treatment combinations $\bm{z}_{1}$, $\bm{z}_{4}$, $\bm{z}_{5}$, and $\bm{z}_{8}$ for which we can build the following estimand:

[TABLE]

Based on the hierarchy principle, an estimator with this aliasing structure should be a superior estimator to the original fractional factorial estimator mentioned because the three-factor interaction is more likely to be negligible than the two-factor interactions. In particular, using the previous numbers from our empirical example, the bias for the main effect of pesticide A would be $1/3$. In general, based on the hierarchy principle, we want to choose a design that aliases the effect we are looking at with the highest-order interaction possible (which may be limited by the missing data structure). For example, if we are only missing data for one treatment combination, when estimating main effects we can always choose the $2^{K-1}$ fractional factorial design that aliases the main effect with the (negative of) the highest-order interaction. This is not a design we would want to use to collect data, because it would alias our main effects with each other. However, differing from the fractional factorial setting, we have more data and would use a different design for estimating each factorial effect.

Denote by $\dot{\bm{g}}_{k}$ the analog of $\bm{g}_{k}$ with zeros where outcome data are missing or excluded in a given estimator for factor $k$. If there is a single treatment combination with no outcome measured, we can choose the aliasing such that factor $k$ is aliased with the $K$-factor interaction. When more rows are missing, the pattern of missingness will dictate what aliasing structure is possible. For example, if we are missing two treatment combinations but they are missing from the same design that aliases factor $k$ with the positive of the $K$-factor interaction, we can still recreate the design that aliases factor $k$ with the negative of the $K$-factor interaction. But if we are missing a row from each of these designs, neither option is possible and we must choose a different aliasing structure. If this method is continued for each factor, one ends up with a set of $\dot{\bm{g}}_{k}$, each with zeros for different treatment combinations.

Section 4.5 of Wu and Hamada, (2000) gives a general strategy to design experiments while attempting to reduce aliasing for certain main effects. There are other designs we can construct, not considered here, such as nonregular designs that have partial aliasing (see Wu and Hamada,, 2000, Chapter 7 for more details on nonregular designs).

3.2.2 Statistical inference

Once we pick which subsets of treatments we will use (based on the aliasing structure and practical constraints), we assume units are assigned using complete randomization to the $\tilde{J}$ treatments in the incomplete design. Let $\bm{W}$ be the $n\times 2^{K}$ matrix of $W_{i}(\bm{z}_{j})$ for all $i$ and $\bm{z}_{j}$, and $\bm{w}$ is some $n\times 2^{K}$ matrix with entries in the set $\{0,1\}$. If $\tilde{\mathcal{J}}$ is the set of indices for the $\tilde{J}$ treatments included in the design, then

[TABLE]

if $\sum_{j\in\tilde{\mathcal{J}}}w_{i}(\bm{z}_{j})=1$ for all $i$ and $\sum_{i=1}^{n}w_{i}(\bm{z}_{j})=n_{j}$ for all $\bm{z}_{j}$, where $n_{j}=0$ for any treatments $\bm{z}_{j}$ not included in this design, and $P(\bm{W}=\bm{w})=0$ otherwise. Assignment is similar to the full factorial or fractional factorial, just for a different subset of treatments. The major differences between a fractional factorial design and incomplete design are therefore not in the way that the treatments are randomized but in (i) that the treatment subset used is not (in general) a fractional subset of the original $2^{K}$ treatments and (ii) the way in which the analysis is proposed to proceed. All of the following inference results are derived under the assumption of complete randomization.

More generally, denote our estimator for $k^{th}$ factor under one of these alternative incomplete factorial designs as $\hat{\dot{\tau}}(k)=\dot{\bm{g}}_{k}^{T}\bar{\bm{Y}}^{obs}$. It is straightforward to extend the variance expression $\text{Var}\left(\hat{\dot{\tau}}(k)\right)$ and variance estimator $\widehat{\text{Var}}\left(\hat{\dot{\tau}}(k)\right)$, as well as the covariance of $\hat{\dot{\tau}}(k)$ and $\hat{\dot{\tau}}(k^{\prime})$, from Section 3.1.2; in fact, it is easy to similarly extend these results to any linear combinations of interest on the $2^{K}$ treatment combinations. However, terms for some treatment combinations will be set to zero in these expressions because they are present in one estimator but not the other. See Appendix F for the specification and derivation.

If we run a regression with all interactions on a dataset with missing treatment combinations, it is ambiguous what design the resulting estimators correspond to and, as discussed above, multiple designs may be plausible to estimate effects well. If we specify a regression of the outcome on the fully interacted treatments, but some treatment combinations are not present in the data, the regression will not be able to estimate all interactions. Not including a set of interactions implies that the linear model assumes that those interaction effects are zero, and we can assess how reasonable this assumption is. However, the specific aliasing structure between the effects included in the model and those dropped by the regression is not obvious from the usual computer output alone, though can be discerned from the design matrix. Byar et al., (1993) and Byar et al., (1995) give more discussion of these types of estimators and Appendix F.2 has further discussion.

3.3 Fractional factorial designs vs incomplete factorial designs

Sometimes an experimenter may have a choice between running a fractional factorial design or an incomplete factorial design. However, in many cases this decision will be dictated in part by the ability to assign at least 2 units to each treatment in the design. The fractional factorial design strictly uses $2^{K-p}$ of the original treatments, whereas an incomplete design may use more treatment groups.

When either design may be run, we expect there to be an advantage in terms of bias by using an incomplete factorial design. In the incomplete factorial example given in Table 3, we were able to show that we can alias the main effect of factor 1 with the three factor interaction. In our numerical example with the three pesticides, this implied a bias for our main effect estimator of pesticide A of $1/3$. On the other hand, the fractional factorial design of Table 2 resulted in factor 1 being aliased by the interaction between factors 2 and 3. In our example, this meant the estimator for main effect of pesticide A was biased by 2. We expect this larger bias by aliasing with two-factor interactions over three-factor interactions in general due to the hierarchy principle.

In terms of variance, consider the case where we have the same total sample size for each design, a balanced design, and additive effects. The fractional factorial design uses all $2^{K-p}$ treatment groups in its estimators of effects, and thus all $n$ units. The incomplete design will only use a subset of the treatment combinations for estimating each effect, and therefore less than the total $n$ units. It should be clear that we expect lower variance under the fractional design, in which more units are used in estimation. Therefore, there is a bias-variance trade-off comparing the fractional factorial to the incomplete factorial design.

If instead we are comparing using different estimators (that correspond to different hypothetical design) after an incomplete design, such as in Table 3, was run, it should be clear that if the design is approximately balanced, we should expect similar variances for estimators that use the same number of treatments.

4 Discussion: Embedding observational studies in fractional factorial designs

4.1 General issues

To address causality in a non-randomized study, it has been argued that one needs to conceptualize a hypothetical randomized experiment that could correspond to the observational data (Bind and Rubin,, 2019; Rosenbaum,, 2002; Rubin,, 2008; Stuart,, 2010). The hypothetical randomization is plausible if the treatment groups are “similar” with respect to confounding variables (Rubin,, 2007, 2008). With multiple treatments of interest, we would aim to recreate a hypothetical factorial randomized experiment. However, in observational studies, there may be no units that received certain treatment combinations. Therefore, we propose to recreate a hypothetical fractional factorial or incomplete factorial experiment, which further motivates the usefulness of these designs beyond randomized experiments.

For the fractional factorial design, we must decide which fraction to use. This decision should be based on criteria such as the maximum resolution criteria. To reduce the amount of aliasing a small $p$, i.e., a small fraction of total design being removed, is desirable. In practice, we may not be able to control the aliasing structure. We must choose a $2^{K-p}$ design such that the treatment combination(s) with no observations is not used. However, this strategy usually results in the removal of units assigned to treatment combinations that were present in the observational data set but not used in the design. If only one treatment combination is not present in the data set, we could use a $2^{K-1}$ design, but in that case we are not using $2^{K-1}-1$ treatment combinations for which we have data.

An alternative strategy to using a subset of the data would be to use an unsaturated regression model. Zhao and Ding, (2021) discusses this model from the potential outcome perspective, though not in the case where there are entire treatment groups missing from the data or containing only a single unit. Those authors find that the unsaturated model will result in bias, in general settings, even with all treatment groups present.

After a design is chosen, a strategy to balance covariates should be used to ensure that the units across treatment groups are similar. We assume strong ignorability; i.e., conditional on measured covariates, the assignment mechanism is individualistic, probabilistic, and unconfounded (Rosenbaum and Rubin,, 1983). Then we can obtain unbiased causal estimates by analyzing the data as if it arose from a hypothetical randomized experiment after balancing.

Importantly, these assumptions apply to all treatment combinations in the final experimental design used. We must also assume that the treatment assignment is, at least hypothetically, manipulable such that all potential outcomes are well-defined. This ensures that our estimands of interest are defined. See Supplementary Material G for more discussion. A similar argument must hold for the unconfoundedness assumption.

4.2 Covariate balance

An important stage when estimating the causal effect of non-randomized treatments is the design phase (Rubin,, 2007, 2008). At this stage, we attempt to obtain a subset of units for which we can assume unconfoundedness of the treatment assignment. That is, units for which $P(W_{i}|Y_{i}(z),\bm{X}_{i})=P(W_{i}|\bm{X}_{i})$, where $\bm{X}_{i}$ is an $m$-dimensional vector of covariates (Imbens and Rubin,, 2015). Matching strategies are often used to ensure no evidence of covariate imbalance between treatment groups, as reviewed by Stuart, (2010).

With multiple treatments, matching can be difficult and methods for obtaining covariate balance specifically for factorial type designs require further exploration. There are extensions of propensity score balancing to multiple treatments, notably the generalized propensity score (GPS) (Hirano and Imbens,, 2004) and the covariate balancing propensity score (CBPS) (Imai and Ratkovic,, 2014). Lopez and Gutman, (2017) review techniques, including matching, for observational studies with multiple treatments and Bennett et al., (2020) uses template matching, which matches units to a “template” population of units, for multiple treatments. Although neither of those works explore factorial (full or fractional) designs, it may be straightforward to extend their methods to this setting. Nilsson, (2013) discusses matching in the $2^{2}$ design. In our data illustration in Section 4.4, we employ sequential trimming and checks on covariate balance. Tests for covariate balance are reviewed in Supplementary Material B.

One issue likely to occur with matching is the difficulty in achieving good balance across many groups. This could result in dropping a large number of units, leading to lower power to detect effects. Weighting techniques can avoid the issue of completely dropping units, but the weights may be unstable if there is little overlap between all treatment groups.

Due to challenges in obtaining full balance, a first step might be to only obtain covariate balance between two treatment groups, a task commonly done in the causal inference literature, and compare outcomes for these groups. To start, we could estimate the difference in outcomes between the units that were assigned level $+1$ for all factors and the units that were assigned $-1$ for all factors. Under certain assumptions, testing the difference in these groups can act as a global test for whether any effects of interest are significant. We discuss the test and assumptions that are needed for this simple comparison to be meaningful in Supplementary Material C.

4.3 Comparing designs

There may be complications in the observational setting that make attempting to recreate one design easier or more desirable than another. If we use a different design for each estimator, as in Section 3.2, we are able to use more of the data than a fractional factorial design, but this can also incur a cost. If we have no observations for only one treatment combination, we would use all $2^{K}-1$ treatment combinations if we did an analysis for all effects and used a different design for each. This approach would require us to either first obtain balance among all $2^{K}-1$ treatment groups or to obtain balance among the treatment groups within each design separately. The former option may be difficult; as the number of treatment combinations grows, obtaining covariate balance across all treatment groups becomes increasingly difficult and may result in smaller and smaller sample sizes, especially if trimming is used. The latter option will make joint inferences more challenging because different units would be used in each analysis. Therefore, although these incomplete factorial designs may improve the bias of our estimators, the fractional factorial design in which we are using the same $2^{K-p}$ rows may be more attractive in terms of obtaining covariate balance.

Further, estimators for different designs result in different types of bias due to aliasing of effects. In order to believe that the biases of main effect estimators are negligible, we must believe that the interactions aliased with that main effect are negligible. Researchers need to carefully unpack the aliasing structure for the design they choose in the observational setting in order to understand the assumptions necessary for unbiased estimation. This requires even more care and thought than a standard observational study.

4.4 Data illustration

4.4.1 Data description

We now give an illustration of the implementation of our methods using data on pesticide exposure and body mass index (BMI), the ratio between weight and height-squared. Previous studies have found an association between pesticide exposures and body mass index (BMI) (Buser et al.,, 2014; Ranjbar et al.,, 2015). We use the 2003-CDC cycle of the National Health and Nutrition Examination Survey (NHANES) collected by the Centers for Disease Control and Prevention (CDC). We access the data via the R (R Core Team,, 2017) package RNHANES (Susmann,, 2016). We focus on four organochlorine pesticides, measured via a blood serum test and then dichotomized based on whether they were above (+1) or below (-1) the detection limit, as given in the NHANES dataset, as factors. Organochlorine pesticides are persistent in the environment and adverse health effects have been reported by the CDC (CDC), making them an interesting group of pesticides to study.

We focus on four pesticides (beta-Hexachlorocyclohexane (beta-Hex), heptachlor epoxide (Hept Epox), mirex, and p.p’-DDT), chosen primarily based on data availability and exposure rates. We removed 271 units with missing values of pesticide and BMI, noting that because this is an illustration and not intended to draw causal conclusions we simply drop those units. We decided to study a non-farmer population as farmers are more likely to be exposed to pesticides than the general population and may also differ on other unobserved covariates that affect health outcomes. This is our first step to achieving covariate balance, leaving a dataset with 1,259 observations (see Figure 5 in the Appendix). Importantly, as shown in Table 4, we have one treatment group with only one unit. Thus, using the full dataset can lead to unstable estimates and an inability to calculate Neyman variance.

4.4.2 Design stage

To show the process of how estimands change as we adjust our design stage, we consider three different hypothetical “experiment” designs:

1. a 24 factorial design, without adjusting the sample to balance for covariates
2. a fractional factorial 24-1 design, without adjusting the sample to balance for covariates
3. a fractional factorial 24-1, with trimming to obtain covariate overlap and no evidence of covariate imbalance with respect to gender (recorded as male vs. female), smoking status (smoker vs. non-smoker), and age at the time of survey (in years). Comparing designs (1) and (2) will show how estimates change between a full factorial and a fractional factorial design. Comparing designs (2) and (3) will show how estimates change when a design phase with the aim of obtaining covariate balance is implemented. We use simple trimming as the focus of this paper is not to compare different methods of balancing for covariates, though exploration of different balancing methods in this setting would be worthwhile. We consider gender, age, and smoking status as our first tier of most important covariates. For smoking status, we use the question “Have you smoked at least 100 cigarettes in your life time?” and categorize “Yes” and “Don’t know” as smokers. Only one observation with value “Don’t know” was recorded.

4.4.3 Statistical analysis

We analyzed the three datasets described in the design stage using:

1. a multiple linear regression that regresses BMI on treatment factors;
2. a multiple linear regression that regresses BMI on treatment factors, as well as the following covariates, as factors: race and ethnicity, income, gender, and smoking status;
3. Fisher-randomization tests of the sharp null hypothesis of no treatment effects. For the first analysis, recall that regression estimates when including all factors and interactions, but not covariates, correspond to the Neymanian estimates divided by two. The second analysis is motivated by Figures 6 and 7 in Appendix H.3, which show that further adjustment for ethnicity and income are needed, even after trimming. The second analysis also adjusts for our second tier of covariates: race and ethnicity and income. Income is defined as annual household income. From now on, we will refer to the “race and ethnicity” covariate as simply ethnicity, as a short hand and to make clear that this is one categorical variable in the NHANES dataset. In the second analysis for each design, units with missing covariate values were removed. Due to the right-skewed nature of the weight variable, BMI exhibits some degree of right-skewness. Therefore, we use log-transformed BMI as the outcome in all of these analyses. Additional data descriptions and the full analyses are available in Appendix H.

4.4.4 Resulting designs

Table 4 provides the counts of observations for each treatment ($\bm{z}_{j}$). Factor combination 10 ($\bm{z}_{10}$) has only one observation. Relying on only one observation for a treatment will lead to unstable estimates and we will not be able to use Neymanian variance estimates. Hence, the fractional factorial design aims to avoid this issue by embedding the observational study in a fractional factorial hypothetical experiment.

Recalling that $I$ corresponds to the intercept, in the $2^{4-1}$ fractional factorial design, instead of using $I=1234$, we chose $I=-1234$ to exclude row 10 in Table 4 that has only a single observation. The dataset in this hypothetical experiment consists of 523 observations. In this design, aliasing is as follows: $I=-1234$, $4=-123$, $3=-124$, $2=-134$, $1=-234$, $12=-34$, $13=-24$, $14=-23$. The main effects are aliased with the negative of the three-factor interactions and the two-factor interactions are aliased with each other with reversed signs. In order to identify main effects, we will assume that the three-factor interactions are negligible. In practice, researchers should assess whether this aliasing assumption is realistic.

For the third design we attempted to get balance on the covariates of gender, smoking status, and age. Figure 1 shows covariate imbalance with respect to gender, smoking, and age in the fractional factorial design before trimming. Recall that these form our first tier of covariates and ethnicity and income constitutes the second tier, which is adjusted for via linear regression. In practice, balancing a large set of covariates is not always possible and so choosing the most important (first tier) covariates should be done using subject matter knowledge. We used a rejection approach that sequentially pruned observations from the fractional factorial dataset until we found no evidence of imbalance across exposure groups with respect to our first tier of covariates, based on a MANOVA using the Wilks’ statistic (Wilks,, 1932, see also Supplementary Material B). Figure 2 shows the covariate distribution for our first tier after trimming. Gender and smoking status are still imbalanced even after trimming, and hence were adjusted for in the linear regression that includes covariates (this is true for all “designs”). The first dataset that resulted in no evidence of covariate imbalance consisted of 169 observations, and the new treatment counts are presented in Table 4. The number of units has been drastically reduced in our attempt to achieve covariate balance, a major challenge in this setting. In fact, one treatment group only has three observations.

4.4.5 Results comparison across different conceptualized experiments and statistical approaches

Figure 3 shows a comparison of the regression estimates of the main effects across designs and statistical analyses. To compare the different methods we present univariate analyses. That is, we utilize individual tests for each main effect rather than joint tests, for better illustration of the different methods. In practice, adjustment for multiple comparisons should be considered. All methods and designs estimated a positive “effect” of heptachlor epoxide and a negative “effect” of mirex on BMI. Although the full factorial estimates generally agree with the estimates of the two fractional factorial designs, differences in estimates of beta-Hexachlorocyclohexane (Beta-Hex) and p,p’-DDT may be due to the aliasing of the three-factor interactions with the main effects. However, it is also plausible that we have reduced our data in the fractional design to a subset of individuals with different average main effects than in the full data set.

Figure 4 shows a comparison of the significance of these estimates. The Fisher p-value is the p-value for the test of no effects of any pesticides, based on effect estimates for a given pesticide, which is suggested as a screening stage in Espinosa et al., (2016). We obtained low p-values for the main effect of mirex on BMI across all methods and designs. However, the p-values disagree for the other pesticides, especially the p-values testing the effects of p,p’-DDT and beta-Hexachlorocyclohexane. The HC2/Neyman p-value is the significance based on the HC2 variance estimate (or Neyman variance estimate as we have shown this estimator to be equivalent in settings with no covariates) and the normal approximation. This p-value was only calculated for the regression analysis without covariates and was unavailable for the full factorial model due to limited data.

4.4.6 Discussion of data illustration

We have performed a data illustration to show the benefits and challenges of using our method and working with observational data with multiple treatments in general.

First, it is important to note that simplifications were made in the statistical analyses to focus on illustrating how researchers can capitalize on using fractional factorial designs to estimate the main and interactive effects of multiple treatments in observational studies. We did not adjust for all important hypothetical covariates, such as diet. Another consideration is that we log-transformed BMI to address the fact that BMI is a ratio and so its distribution tends to have heavy tails. However, heavy-tailed distributions could have been considered.

There were also some major challenges in working with this data and our method related to sample size. Trimming helps mitigate bias but greatly reduced our sample size, potentially leading to decreased power and precision. There was also a great reduction in sample size by using the fractional factorial design which entailed dropping treatment combinations. We could have used an incomplete factorial type design, but this likely would have made the covariate balancing even more challenging, as discussed previously.

For p,p’-DDT, the full factorial design resulted in a lower p-value than the other designs, which could be an issue of aliased interactions watering down the effect in the fractional design, as a result of the bias noted in Result 1. It could also have occurred because the populations are different in these two designs. For beta-Hexachlorocyclohexane, the covariate balance adjusted fractional factorial design has a lower p-value compared to the other designs. We are more inclined to trust the balanced design as this should have reduced bias. This result may indicate that there was some confounding that made the effect appear less significant before trimming. In fact, the stark contrast between the unbalanced and balanced fractional design suggests that confounding may be to blame. Alternatively, by trimming we may have reduced our sample to a subpopulation where beta-Hexachlorocyclohexane has a larger effect on BMI than the rest of the population.

It is important to acknowledge that using the small subset of the NHANES dataset, we do not intend to provide policy recommendations on pesticide use. In the general population, organochlorine pesticide exposure primarily occurs through diet (excluding those with farm-related jobs), particularly eating foods such as dairy products and fatty fish (Centers for Disease Control and Prevention , 2009, CDC). Without further adjustments for diet, we are not be able to disentangle the causal effect of diet and pesticides. For instance, in our study individuals are likely to have been exposed to mirex largely through fish consumption (Agency for Toxic Substances and Disease Registry [ATSDR], ATSDR). Further studies could investigate BMI differences in a group of fish consumers with high levels of mirex and a “similar” group of fish consumers with low levels of mirex, where similar is with respect of important confounding variables. Indeed, it could be that eating fish causes individuals to have both high levels of mirex and also lower BMI.

5 Epilogue

In this paper, we have laid out how to design and analyze fractional factorial and incomplete factorial designs using the potential outcomes framework. These designs are particularly useful when running a full factorial design is infeasible or impractical, but the factorial estimands are still of interest. Our work includes extensions of some of the known factorial design results for variance and regression to the fractional factorial and incomplete design setting. In doing so, we have expanded the tool kit available to researchers wishing to perform causal inference.

Further, we have proposed to embed observational studies with multiple treatments in fractional factorial hypothetical experiments. This type of design is useful in settings with many treatments, especially when some treatment combinations have few or no observations and the aliasing assumptions are plausible. Once we recreate a factorial or fractional factorial experiment in the design phase, we can use standard methods, extended as in Sections 2 and 3, to estimate causal effects of interest. We first reviewed the basic setup for factorial and fractional factorial designs. We have also discussed covariate balance complications that may arise when dealing with multiple non-randomized treatments in practice and illustrated these methods on a data set with pesticide exposure and BMI.

We have given a short overview of factorial, fractional factorial, and incomplete factorial designs, in the potential outcomes framework. However, there are many related designs that we did not cover (see Wu and Hamada,, 2000), such as nonregular design types (e.g., Plackett-Burman designs). Practitioners may also use variable selection and the principle of effect heredity to select their model for estimating factorial effects.

We see many avenues of future exploration connected to our approach. For instance, coupling fractional factorial designs with a Bayesian framework would provide more statistical tools and would potentially offer different methodology for dealing with missing data. Additionally, development of observational tools, especially covariate balancing techniques for factorial designs with many treatment combinations should be an area of future exploration. A particular challenge is that as we increase the number of treatment combinations, matching becomes more and more difficult due to the increased dimensionality and weighting methods may produce unstable estimators. One direction could also be to choose the design of the observational study based upon the ability to balance different treatment groups. Random allocation designs (Dempster,, 1960, 1961), in which randomness of the design is incorporated, could also be utilized in this framework. Finally, we could explore other causal estimands that may be of interest in observational studies with multiple treatments, such as those in Egami and Imai, (2019) and De la Cuesta et al., (2019).

Figures

Appendix A Statistical inference for full factorial

In this section we provide some further review of full factorial designs under Neyman-style randomization-based inference. A conservative Neyman-style estimator for the variance was proposed by Dasgupta et al., (2015) and extended to the unbalanced case by Lu, 2016b :

[TABLE]

where

[TABLE]

is the estimated sampling variance of potential outcomes under treatment combination $\bm{z}_{j}$. Dasgupta et al., (2015) discussed other variance estimators and their performance under different assumptions (e.g. strict additivity, compound symmetry). In that paper, the authors also provided a Neyman-style estimator for the covariance by again substituting $s^{2}(\bm{z}_{j})$ for $S^{2}(\bm{z}_{j})$. Unfortunately, this estimator is not guaranteed to be conservative because $S^{2}_{k,k^{\prime}}$ may be positive or negative. The same authors provided a Neyman confidence region for $\bm{\tau}$, the vector of all $J-1$ factorial effects. First, they define, in Equation 26,

[TABLE]

where $\hat{\bm{\tau}}$ is the vector of Neyman estimators of $\bm{\tau}$ that we defined earlier and $\hat{\Sigma}_{\hat{\bm{\tau}}}$ is the estimator of the covariance matrix of $\hat{\bm{\tau}}$, $\Sigma_{\hat{\bm{\tau}}}$. Then the $100(1-\alpha)\%$ confidence region for $\bm{\tau}$ (Equation 27 of Dasgupta et al., (2015)), is

[TABLE]

where $0<\alpha<1$ and $p_{\alpha}$ is the is the $\alpha$ quantile of the asymptotic distribution of $T_{n}$. In that paper, a Neyman-style confidence interval for individual effects was also provided as

[TABLE]

where the authors rely on a normal approximation for the distribution of $\hat{\tau}(k)$, typically assumed to hold asymptotically. See Proposition 3 and Example 11 in Li et al., (2020) for conditions under which the asymptotic normality holds in this setting.

Note that a model that linearly regresses an outcome against the factors coded using contrast values $-1$ and $+1$ and all interactions between factors (but no other covariates) will result in the same point estimates for the factorial effects as presented here, divided by 2 (Dasgupta et al.,, 2015; Lu, 2016b, ). For a balanced design or when treatment effects are assumed to be additive, so that the variances $S^{2}(\bm{z}_{j})$ are all the same, the standard linear regression variance estimate, relying on homoskedasticity, will be the same as the Neymanian variance estimate (Dasgupta et al.,, 2015; Imbens and Rubin,, 2015). However, this is not true for an unbalanced design. Samii and Aronow, (2012) showed that the HC2 heteroskedasticity robust variance estimator (see MacKinnon and White,, 1985, for more details) is the same as the Neymanian variance estimator presented here for a single treatment experiment. Lu, 2016b extended this finding to factorial designs, showing that this estimator is also the same as the Neymanian variance estimator in that case. Note that because the regression estimator is different by a factor of 2, this regression variance is different by a factor of 4.

Fisherian and Bayesian types of analyses are also possible. See Dasgupta et al., (2015) for additional discussion. In particular, they explore creating “Fisherian Fiducial” intervals (Fisher,, 1930; Wang,, 2000) for factorial effects. The basic idea is to invert our estimands so that we write our potential outcomes in terms of the estimands. Then, under a Fisher sharp null hypothesis, we can impute all missing potential outcomes and generate the randomization distribution.

Appendix B Tests for covariate balance

Testing for covariate imbalance across multiple treatment groups can be done using multivariate analyses of variance (MANOVA), which uses the covariance between variables to test for mean differences across treatment groups, as used in Branson et al., (2016). Recall that the factorial design has $J$ treatment combinations. Define the following H and E matrices (Coombs and Algina,, 1996):

[TABLE]

[TABLE]

where $\overline{\bm{X}}_{k}$ is the $m$-dimensional vector of mean covariate values for treatment group $k$, $\overline{\bm{X}}$ is the average $m$-dimensional vector of mean covariate values for all units, that is $\overline{\bm{X}}=\sum_{k=1}^{J}\frac{n_{k}}{n}\overline{\bm{X}}_{k}$, and $\bm{X}_{kj}$ is the $m$-dimensional vector of covariates for the $j$th unit in treatment group $k$. Denote by $\theta_{k}$ the ordered eigenvalues of $\bm{H}\bm{E}^{-1}$, where $k\in\{1,...,s\}$ and $s=\min(m,K-1)$. Standard MANOVA statistics, which can be used to test covariate balance, are typically functions of the eigenvalues of $\bm{H}\bm{E}^{-1}$ (Coombs and Algina,, 1996). We chose the Wilks’ statistic (Wilks,, 1932),

[TABLE]

where $\theta_{k}$ corresponds to the $k^{th}$ eigenvalue of $\bm{H}\bm{E}^{-1}$.

As discussed in Imai et al., (2008), a potential drawback of testing for evidence against covariates imbalance is that as we drop units we lose power to detect deviations from the null hypothesis of no difference in covariates between the treatment groups. Another diagnostic for covariate balance is checking covariate overlap via plots and other visual summaries of the data. So called “Love plots”, which show standardized differences in covariate means between two treatment groups before and after adjustment (Ahmed et al.,, 2006), are difficult to generalize directly because of the multitude of treatment groups and comparisons. However, plots of standardized means or of distributions may be helpful to detect imbalance.

Appendix C Initial test for significance of effects

As discussed in the previous section, achieving covariate balance for more than two treatment groups can be a challenge. Therefore, instead of attempting to achieve balance among all treatment groups, a simple first step might be to examine two carefully chosen treatment groups and attempt to balance these two groups only. Obtaining balance between two treatment groups has been well-studied in causal inference; see Stuart, (2010) for a review of common matching methods. Once significant covariate imbalance can be ruled out, we can test whether the mean difference between these two groups is significantly different. But what can we learn from this comparison about our factorial effects?

We have from Section 2.1 that

[TABLE]

So when we subtract two observed means, assuming that the observed means are unbiased estimates of the true means (i.e. we have randomization or strong ignorability), we are estimating

[TABLE]

which is the sum of terms that are signed differently in $\bm{h}_{j}$ and $\bm{h}_{j^{\prime}}$.

As an example for a $2^{3}$ design, we have the following matrix for $\bm{G}^{T}$:

[TABLE]

Now consider taking the difference between $\bar{Y}(\bm{z}_{8})-\bar{Y}(\bm{z}_{1})=\bar{Y}(+1,+1,+1)-\bar{Y}(-1,-1,-1)$, which yields

[TABLE]

Hence, testing whether the difference between $\bar{Y}(\bm{z}_{8})$ and $\bar{Y}(\bm{z}_{1})$ is zero is the same as testing whether $\tau(1)+\tau(2)+\tau(3)+\tau(123)$ is zero. If it is reasonable to assume that all main effects are of the same sign based on subject matter knowledge and that the three-factor interaction is also of the same sign or negligible, then this global test of whether there are any treatment effects is relevant to the estimand of interest. According to the effect hierarchy principle (Wu and Hamada,, 2000), the main effects should dominate the three-factor interaction. Thus, even if the interaction differs in sign, we would expect to see an effect under this assumption if there is one. If the global test is not rejected, then we would move on to recreating the entire factorial design. If it is unclear whether the signs of the effects are the same (also referred as antagonistic effects), then this global test would be inappropriate because effects could still be different from zero but their sum could be (close to) zero.

If we are particularly interested in the causal effects involving the first factor, we can create a test specifically for those effects. In the the $2^{3}$ setting, we might use the estimand $\bar{Y}(\bm{z}_{8})-\bar{Y}(\bm{z}_{4})$ as a proxy for the effect of the first factor, as follows:

[TABLE]

In Equation 8, if all terms have the same sign and we find that the difference is significantly different than zero, then we would conclude that factor 1 has an effect, either on its own or through interactions with the other factors. A different choice of levels for the other factors, for instance comparing the mean potential outcome when one factor is high vs. low when all other factors are at the low level, would result in different signs for the interactions. The choice of which levels to compare should be based on subject matter knowledge and reasonable assumptions.

To reiterate, throughout this section we focused on the meaning of several estimands, which are easily estimated under a randomized experiment. In an observational study, we would first need to obtain balance for any treatment groups we would be using in our estimator.

Appendix D Variance derivations

D.1 Variance and covariance for fractional factorial design

This section gives details on the derivation of the variances and covariances of our estimators of factorial effects under a fractional factorial design, as given in Equation 3. The proofs are direct from Theorem 3 in Li and Ding, (2017) as we can express our estimands and estimators as linear combinations of potential outcomes. In particular, let $\bm{A}_{\bm{z}^{*}_{j}}$ be the vector with entries being the $j$th entries in the $\bm{g}^{*}$ for all factorial effects:

[TABLE]

Then $\left(\tau(1),\dots,\tau(K),\tau(1,2),\dots,\tau(1,2,\dots,K)\right)^{T}=\sum_{j=1}^{J^{\prime}}\bm{A}_{\bm{z}^{*}_{j}}\bm{Y}(\bm{z}^{*}_{j})$

Then by Theorem 3 of Li and Ding, (2017),

[TABLE]

Similarly, by Theorem 3 of Li and Ding, (2017)

[TABLE]

Appendix E Relating linear regression estimators to Neyman estimators in the fractional factorial design

This section gives a brief overview of a proof that the linear regression point estimates are the same as the Neymanian estimates for the fractional factorial design. For proofs of these results for the full factorial design, see Lu, 2016b . The linear regression coefficient estimate is

[TABLE]

where $\bm{X}$ is a $n\times 2^{K-p}$ matrix whose columns correspond to first an intercept, then include the main effects, then the second order interactions not aliased with the main effects, and so on such that no two columns are aliased. For the $2^{3-1}$ design laid out in Section 3, $\bm{X}$ would be a design matrix with a first column of 1’s and the rest of the columns corresponding to levels of the first, second, and third factor. No interactions would be included in this example because each of the two-factor interactions is aliased with a main effect and the three-factor interaction is aliased with the intercept, which are all already included in the model. Thus the design looks like the columns $\bm{g}^{*}_{1}$, $\bm{g}^{*}_{2}$, and $\bm{g}^{*}_{3}$ of Table 2, but with repeated rows for each of the units assigned to the same treatment combination. We need $(\bm{X}^{T}\bm{X})^{-1}\bm{X}^{T}\bm{X}=\bm{I}$. Denote $\bm{B}=(\bm{X}^{T}\bm{X})^{-1}\bm{X}^{T}$.

[TABLE]

where $\bm{b}_{i\cdot}$ is the $i$th row of $\bm{B}$ and $\bm{f}_{j}$ is the $j$th column of $\bm{X}^{T}$. Using notation from Section C, let $\bm{h}^{*}_{j}$ be the $j$th row of $\bm{G}^{*T}$ which is an expanded version of $\bm{z}^{*}_{j}$ which includes elements for interactions and the intercept. Then $f_{ji}$, the $i$th element of $\bm{f}_{j}$ is $h^{*}_{kj}$, the $j$th element of $\bm{h}^{*}_{k}$, where $k$ is the treatment combination for the $i$th individual. It must be true that $(\bm{B}\bm{X})_{ii}=1$ and $(\bm{B}\bm{X})_{ij}=0$ for $i\neq j$. Consider letting the $i$th row of $\bm{B}$ be $\bm{b}_{i\cdot}=\frac{1}{2^{K-p}}(\bm{\tilde{n}}^{-1}\circ\bm{f}_{i})^{T}$ where $\bm{\tilde{n}}^{-1}$ is the vector whose $i$th entry is $\frac{1}{n_{j}^{*}}$ where $j$ is number of the treatment combination that the $i$th unit is assigned to ($j=\sum_{k=1}^{J^{\prime}}kW_{i}(z_{k})$). Then we have

[TABLE]

For $k\neq j$

[TABLE]

So the $k$th row of $\bm{B}$ is $\frac{1}{2^{K-p}}(\bm{\tilde{n}}^{-1}\circ\bm{f}_{k})^{T}$. This means that

[TABLE]

So, indeed, the linear regression estimates are one half of the factorial effects.

Now let us consider the HC2 variance estimator. It has the form (MacKinnon and White,, 1985)

[TABLE]

where $\widehat{\Omega}=\text{diag}\left(\frac{\hat{e}_{i}}{1-h_{ii}}\right)$ with $\hat{e}_{i}$ being the residual for observation $i$ and $h_{ii}$ being the $ii$ value of the hat matrix, $\bm{X}\left(\bm{X}^{T}\bm{X}\right)^{-1}\bm{X}^{T}$. For discussion of this estimators in the single treatment case and the factorial case, see Samii and Aronow, (2012) and Lu, 2016b , respectively. We use similar ideas to those papers in the arguments below.

If unit $i$ was assigned to treatment $\bm{z}^{*}_{k}$ then the $i$th column of $\left(\bm{X}^{T}\bm{X}\right)^{-1}\bm{X}^{T}$ is $\bm{b}_{\cdot i}=\frac{1}{2^{K-p}}\frac{1}{n^{*}_{k}}\bm{h}^{*T}_{k}$. We have

[TABLE]

So then $1-h_{ii}=1-\frac{1}{n_{k}^{*}}=\frac{n_{k}^{*}-1}{n_{k}^{*}}$. This in turn means that

[TABLE]

Now we can solve for the whole expression of $\widehat{Var}_{HC2}(\widehat{\bm{\beta}})$. We focus on the diagonal entries.

[TABLE]

where $\bm{s}_{1}$ is the vector of whose $i$th entry, given entry $i$ is assigned treatment combination $\bm{z}^{*}_{k}$, is $\frac{\left(Y_{i}(\bm{z}_{k}^{*})-Y_{i}(\bm{z}_{k}^{*})\right)^{2}}{n_{k}^{*}-1}$. Thus, we have that $\left(\widehat{Var}_{HC2}(\widehat{\bm{\beta}})\right)_{kk}$ is $1/4$ times the Neyman style variance estimator.

Appendix F Incomplete factorial designs

F.1 Inference

This section discusses alternative incomplete factorial designs that one might use. See Byar et al., (1993) for more discussion of incomplete factorial designs. This discussion also applies more generally to recreating fractional designs in observational settings where we are using a subset of treatment combinations present in the data. We follow the same outline as proofs from Section D.1. In these incomplete factorial designs, for each main effect or other estimand of interest, we define a new experimental design tailored to estimating that effect. When this is done, we then analyze the data as if the treatment groups within the new design are the only possible treatment groups.

First we need to introduce some notation. Let $\dot{\bm{g}}_{j}$ be the same as $\bm{g}_{j}$ but with zero elements corresponding to treatment combinations that are not included in the particular design in use. Let $2^{m}$ treatment groups be used in the design, with half assigned to the $+1$ level of factor one and half assigned to the $-1$ level of factor one. Let $\dot{\tau}(k)=2^{-(m-1)}\dot{\bm{g}}^{T}_{k}\bar{\bm{Y}}$. We can write this using matrix form to apply Theorem 3 of Li and Ding, (2017). Define

[TABLE]

Then $\left(\dot{\tau}(1),\dots,\dot{\tau}(K),\dot{\tau}(1,2),\dots,\dot{\tau}(1,2,\dots,K)\right)^{T}=\sum_{j=1}^{J}\dot{\bm{A}}_{\bm{z}_{j}}\bm{Y}(\bm{z}^{*}_{j})$

Then we have

[TABLE]

An important note is that in terms of estimation of this variance, whether we analyze the data as if the treatment levels in this design are the only possible treatment combinations or if we keep the assumption that units can be assigned to any possible treatment combinations (which aids in the interpretation and inference), the variance estimator will be the same. This is because we can only estimate the first term in this expression which only involves the specific treatment levels in this design.

F.2 Regression with missing levels

This section discusses what will result from a standard regression interacting all factors when not all treatment combinations are observed in the data.

If the dataset is missing $m$ treatment combinations, then the regression will be able to estimate the first $2^{K}-m$ effects (including interactions) that are not aliased and the rest will be removed due to collinearity of the matrix columns. Then for each effect, there will be some aliasing structure imposed but the same “design” will not necessarily be used for each factor.

To explore this scenario more, let’s take the specific example of three factors where we only observe five of the eight treatment combinations. For simplicity let there be one observation for each treatment combination and let the model matrix be as follows (we expect similar results to hold in more realistic scenarios given we use a saturated model which allows a separate mean for each treatment group):

[TABLE]

The first column corresponds to the intercept, the second through fourth columns give the levels for the three factors, and the fifth column corresponds to the interaction between the first and second factor. Note that the first four rows correspond to a $2^{3-1}$ design, so it is possible to recreate that design here.

We have

[TABLE]

Recalling that $\hat{\bm{\beta}}=\left(\bm{X}^{T}\bm{X}\right)^{-1}\bm{X}^{T}\bm{Y}^{obs}$, the first three columns correspond to estimates we would get using the fractional factorial design using the defining relation $I=123$. The last two estimates, for factor 3 and the interaction between factors 1 and 2, have a different aliasing structure. In particular, the aliasing on factor 3 is similar to aliasing structures in Section C and we can find that factor 3 will be aliased with the two-factor interactions 13 and 23 as well as the three-factor interaction.

Appendix G Further comments on observational studies

As noted in the main text, we must assume that the treatment assignment is, at least hypothetically, manipulable such that all potential outcomes are well-defined. Although our fractional factorial estimators do not use all treatment combinations, the estimands and aliasing structure depend on potential outcomes for the unobserved treatment combinations being well defined. If this assumption on manipulation does not hold, we would need to consider a causal estimand that does not depend on unmeasurable or undefined potential outcomes. Thus, it is important that subject-matter knowledge guides us in deciding which covariates are informative about which units have well-defined potential outcomes under all treatment combinations. It may be that only a sub-population has potential outcomes for all treatments well-defined, similar to analyses for noncompliance in which the estimand is only well-defined for compliers. For other groups where not all treatment combinations could theoretically be received, there may be other interesting estimands based on the well-defined potential outcomes, but those are not our focus.

Remark*.*

Under some assumptions, the estimand based on only well-defined potential outcomes may be the same as the original. If we have a $2^{2}$ experiment and there is no interaction between factor 1 and factor 2, then $\tau(1)=\bar{Y}(+1,+1)-\bar{Y}(-1,+1)=\bar{Y}(+1,-1)-\bar{Y}(-1,-1)$. Hence, even if we only observed one level of factor 2 but both levels of factor 1, we could recover $\tau(1)$. However, $\tau(1)$ as previously defined would not be of scientific interest if all potential outcomes are not defined.

Appendix H Data illustration

This section gives some additional descriptors for the data.

Figure 5 shows that the proportion of farmers is different across the eight treatments of the hypothetical fractional factorial experiment defined in Section 4.4.2.

We can also examine the correlations in occurrence of the different pesticides. beta-Hexachlorocyclohexane (beta-Hex) and heptachlor epoxide have the highest correlation of 0.48. Mirex and heptachlor epoxide have the lowest correlation of 0.16.

For each design, we explore regression (both saturated, i.e. with all interactions among factors, and unsaturated, i.e. without interactions among factors), regression with covariates (both saturated and unsaturated for the factors), and Fisher tests for significance of effects. We also include the HC2 standard error estimate for the saturated model without covariates, where it is equivalent to the Neyman estimator, when possible. Note that all regression outputs are using log BMI as output. An important note is that individuals with missing values, either for factor levels, treatment levels, or covariates where they are used, were removed. This action likely changes the types of individuals within the analysis and therefore the generalizability of the results. However, as this is intended as simply an illustration of the methods and not as a full analysis to draw substantive conclusions from, this simple model suffices to allow us to continue with the analysis.

Code is available for those interested in seeing the outputs from all of the analyses. Some notes are provided for the various outputs below.

H.1 Full factorial design

We start by ignoring our limited data and use a full $2^{4}$ factorial approach. Note that the individual who received treatment combination $(-1,1,1,-1)$ has leverage 1 because they are the only individual with that treatment combination. Because of this data limitation, estimating variance using the HC2/Neyman variance estimator is not possible for the saturated case. Also note that in the saturated model, the variance estimates given in the linear model summary are all the same. This will be true of Neyman variance estimates too, since they are the same for each factorial effect estimators. We see changes in significance of the estimator and even a change in sign for the estimate of beta-Hexachlorocyclohexane (beta-Hex) going from the model with just main effects to the saturated model.

Next we analyze the full factorial design, adjusting for the covariates of income, ethnicity, gender and smoking status as linear factors in the regression model. For simplicity, we remove all individuals who had missing values for income or ethnicity, and assume that those values are missing at random. This resulted in 75 units being removed. In practice one should instead use multiple imputation to account for the missing values. One individual refused to give income (the only such individual) and so was removed. This did not affect the analysis. The unit assigned to the unique treatment combination, necessarily had leverage one in the saturated model.

Next we turn to the Fisherian analysis We assume the sharp null hypothesis of zero individual factorial effects, for all factors and interactions. This means that imputed missing potential outcomes for different assignments are just the observed potential outcomes. We do the imputation, or effectively rearrange the assignment vector, 1000 times. We find that only heptachlor epoxide (Hept Epox) and mirex appear to be significantly different from zero at the 0.05 level. We only examine the main effects here but could further consider interactions.

H.2 Fractional factorial design

This section gives the analysis of the fractional factorial design as is, using regression. Note that we remove farmers again, leaving 523 observations. The results without any covariate adjustment generally align with the full factorial analysis, in terms of sign and significance of terms. Again it makes sense that the standard errors estimates will be the same for all estimates in the saturated model because the same groups are being used to calculate them.

Next we turn to the analysis of the fractional factorial design, adjusting for the covariates of income, ethnicity, gender and smoking status as linear factors in the regression model. We again removed units who were missing values for income or race, reducing the sample size by 34. One unit replied “Don’t know” for income so this unit was removed. This did not affect the analysis.

With the Fisherian analysis, once again, only heptachlor epoxide (Hept Ex) and mirex appear to be significantly different than zero at the 0.05 level.

H.3 Fractional factorial with covariate adjustment

Finally, we have the analysis of the fractional factorial design after trimming to attain balance, using regression. 169 units are retained after trimming. We started with regressions without covariate adjustment.

We next performed the analysis on the trimmed data set adjusting for ethnicity, income, gender, and smoking status as linear factors in the regression. Figures 6 and 7 shows the balance across treatment groups for income and ethnicity after trimming, indicating further need to adjust. We see from Figure 2 that gender and smoking also require further adjustment, even after trimming.

Units who had missing values for income or ethnicity were removed, which reduced the sample size by 9. We found that two units had unique income values of “Over $20,000” and “0 to$4,999”. These units were removed. This changes the baseline level for income in the analysis from “$0 to$4,999” to “$5,000 to$9,999.”

In the Fisherian analysis, mirex is the only pesticide that appears to have a significant effect on BMI at the 0.05 level.