2.1 The General Framework

Researchers often rely on unconfoundedness assumptions to infer causality from observational data, where they do not directly control the treatment assignment mechanism. A canonical example of such an assumption, often called conditional ignorability, posits that the treatment assignment is independent of the potential outcomes given a set of observed pre-treatment confounders, i.e.,

(1)

where $D_i \in \{0,1\}$ is the treatment assigned to unit i, $ Y_{i}(d) $ is unit i’s potential outcome under treatment $d \in \{0,1\}$ , and $X_i$ is a vector of observed pre-treatment covariates. In many real-world applications, however, this assumption may not hold. A typical example is observational studies with non-randomized treatment assignment, where researchers cannot rule out the possibility of unobserved confounders correlated with both the treatment and potential outcomes. Another example, which we focus on in this article, is longitudinal data where researchers collect repeated observations for each unit. In such data, unobserved unit-level heterogeneity in the potential outcomes and time-varying common shocks to the outcomes may also affect the adoption of treatment.

A common approach to addressing the violation of an unconfoundedness assumption is sensitivity analysis. Specifically, instead of the relationship (1), it is assumed that there is an unobserved confounder $ U_i $ such that

(2)

which is also known as “latent ignorability” (Rubin Reference Rubin1976). Various sensitivity analysis methods have been developed on the basis of assumptions similar to (2) (e.g., Imbens Reference Imbens2003; Oster Reference Oster2019; VanderWeele and Ding Reference VanderWeele and Ding2017). In this article, we propose a fully Bayesian approach, which directly builds upon the Bayesian causal inference framework of Rubin (Reference Rubin1978) (see also Imbens and Rubin Reference Imbens and Rubin1997; Rubin et al. Reference Rubin, Wang, Yin and Zell2008). Specifically, we consider causal inference to be the problem of imputing missing potential outcomes ( $Y_{mis}$ ) given observed data (D, X and $Y_{obs}$ ) as well as the maintained assumptions. Under latent ignorability, we have

(3) $$ \begin{align} \begin{aligned} &Pr(Y_{mis} | D, X, Y_{obs}) \\& \quad =\frac{Pr(D, X, Y(0), Y(1))}{Pr(D, X, Y_{obs})} \\& \quad \propto Pr(D, X, Y(0), Y(1)) \\& \quad =\int Pr(D, X, U, Y(0), Y(1)) dU \\& \quad =\int Pr(D | X, U, Y(0), Y(1)) Pr(X, U, Y(0), Y(1)) dU \\& \quad = \int Pr(D | X, U) Pr(X, U, Y(0), Y(1)) dU, \end{aligned} \end{align} $$

where the last equality holds because of (2). Note that $ Pr(D | X, U) $ is the treatment assignment mechanism and $ Pr(X, U, Y(0), Y(1)) $ is what Rubin et al. (Reference Rubin, Wang, Yin and Zell2008) call the “science.” Thus, the posterior predictive distribution of the missing outcome $Pr(Y_{mis} | D, X, Y_{obs})$ is proportional to the treatment assignment mechanism multiplied by the underlying science, marginalizing over the unobserved confounder. Many existing sensitivity analysis techniques implicitly utilize this general relationship with additional parametric assumptions for $ Pr(D | X, U) $ and $ Pr(X, U, Y(0), Y(1))$ (e.g., Imbens Reference Imbens2003).

Here, we pursue an alternative factorization which establishes an explicit connection between Rubin’s Bayesian causal inference framework and the Bayesian, non-counterfactual framework for sensitivity analysis (McCandless et al. Reference McCandless, Gustafson and Levy2007). That is, we have

(4) $$ \begin{align} &Pr(Y_{mis} | D, X, Y_{obs}) \nonumber\\& \quad \propto \int Pr(D | X, U) Pr(X, U, Y(0), Y(1)) dU \nonumber\\& \quad = \int \frac{Pr(D, X, U) }{Pr(X, U) } Pr(X, U, Y(0), Y(1)) dU \\& \quad = \int Pr(U|D, X) Pr(D, X) \frac{Pr(X, U, Y(0), Y(1))}{Pr(X, U)} dU \nonumber\\& \quad \propto \int Pr(U|D, X) Pr(Y(0), Y(1)|X, U) dU, \nonumber \end{align} $$

suggesting an alternative procedure for inference based on the imputation of the latent variable U as a function of D and X, as well as a model for the potential outcomes. Indeed, the standard BSA framework utilizes this imputation-based strategy for inference (Gustafson et al. Reference Gustafson, McCandless, Levy and Richardson2010; McCandless et al. Reference McCandless, Gustafson and Levy2007), first imputing U using a parametric model for $Pr(U|D, X)$ , and then including the imputed values of U in an outcome model $f(Y_{obs}|D, U, X)$ . That is, a model for the observed data $(Y_{obs}, D, X)$ can be obtained by integrating U out:

(5) $$ \begin{align} \begin{aligned} Pr(Y_{obs}, D, X) &\propto Pr(Y_{obs}| D, X) \\&=\int Pr(Y_{obs}, U|D, X) dU \\&\propto \int Pr(U|D, X) Pr(Y_{obs}|D, X, U) dU. \end{aligned} \end{align} $$

Comparing (5) to the posterior predictive distribution (4) reveals that the standard BSA framework can be used for sensitivity analysis with respect to causal quantities involving missing potential outcomes. That is, a latent ignorability assumption, such as (2), implies sufficient constraints on the joint potential outcome distribution conditional on X and U so that a causal quantity of interest is identified by the distribution of the observed outcome $Y_{obs}$ conditional on D, $X,$ and U.

2.2 Setup and Notation

Now, we set up our notational framework. Suppose that we observe a panel data set that consists of N units and spans over T periods, with $i = 1, 2, \ldots , N$ and $t = 1, 2, \ldots , T$ indexing unit and time, respectively. We are interested in the effect of a binary treatment indicator, $D_{it} \in \{0, 1 \} $ , on some outcome variable $Y_{it} \in \mathbb {R}$ . For notational convenience, we assume that the panel data set is balanced, i.e., there are no missing observations. We adopt the potential outcomes framework (Rubin Reference Rubin1974) to define the causal quantities of interest. We denote by $\boldsymbol {D}_i = (D_{i1}, D_{i2}, \ldots , D_{iT}) $ the full treatment assignment vector for unit i and $\boldsymbol {D}_i^t = (D_{i1}, D_{i2}, \ldots , D_{it})$ the treatment assignment vector up to time t for unit i. $Y_{it} (\boldsymbol {d}_i) $ and $Y_{it}(\boldsymbol {d}^{\prime }_i)$ are two different potential outcomes for unit i in period t given treatment assignment vectors $\boldsymbol {d}_i = (d_{i1}, d_{i2}, \ldots , d_{iT}) $ and $\boldsymbol {d}^{\prime }_i = (d^{\prime }_{i1}, d^{\prime }_{i2}, \ldots , d^{\prime }_{iT}) $ . We make the following assumptions.

Assumption 1 No anticipation of future treatments

Treatment assignment in the future does not affect current potential outcomes, i.e.,

(6) $$ \begin{align} Y_{it} (\boldsymbol{d}_i) = Y_{it} (\boldsymbol{d}^{\prime}_i) \quad \text{if} \quad \boldsymbol{d}_i^t = \boldsymbol{d}_i^{'t} \end{align} $$

for any $i\in \{1,\ldots ,N\}$ and $t \in \{1, \ldots , T\}$ .

Assumption 2 No carryover effects

Past treatment assignments do not affect current potential outcomes, i.e.,

(7) $$ \begin{align} Y_{it} (\boldsymbol{d}_i^t) = Y_{it} (\boldsymbol{d}_i^{'t}) \quad \text{if} \quad d_{it} = d^{\prime}_{it} \end{align} $$

for any $i\in \{1,\ldots ,N\}$ and $t\in \{1,\ldots ,T\}$ . Footnote ¹

Under these assumptions, we can re-write the potential outcome as $ Y_{it} (\boldsymbol {d}_i) = Y_{it} (d_{it}) $ , and there are two potential outcomes $Y_{it}(1)$ and $Y_{it}(0)$ for each unit in each period. The individual treatment effect is defined as the difference in potential outcomes:

(8) $$ \begin{align} \delta_{it} = Y_{it}(1) - Y_{it}(0). \end{align} $$

We focus on the identification and estimation of the average treatment effect on the treated (ATT), defined as $ ATT = \frac {\sum _{i,t} D_{it} \delta _{it} }{\sum _{i,t} D_{it} } $ . The results also apply to other convex combinations of individual effects, like ATT in a specific period relative to the occurrence of the treatment.

We consider three types of unobserved confounders: unit-specific time-invariant confounders $\gamma _i$ , time-specific unit-invariant confounders $f_t$ , and confounders that vary both across units and time $U_{it}$ . Below, we will develop a method for sensitivity analysis based on the following latent ignorability assumption with respect to these unobserved confounders.

Assumption 3 Panel latent ignorability

(9)

for $i\in \{1,\ldots ,N\}$ and $t\in \{1,\ldots ,T\}$ , where $X_{it}$ is a vector of observed confounders.

2.3 Causal Panel Data Models

The panel latent ignorability assumption implies that we need to control for observed confounders $X_{it}$ as well as the unobserved confounders $\gamma _i$ , $f_t$ , and $U_{it}$ to identify the effect of $D_{it}$ on $Y_{it}$ . Allowing $U_{it}$ to be time-varying enables the model to capture both unobserved time-varying confounders and the time-varying influence of time-invariant confounders, making it a more flexible specification. Suppose the observed confounder $X_{it}$ is a vector of length P, i.e., $ X_{it} = (X_{it,1}, X_{it,2}, \ldots , X_{it,P}) $ , where each element $X_{it,p}$ is normalized to have zero mean and unit variance. Existing literature has proposed both design-based approaches like doubly-robust estimator (Arkhangelsky and Imbens Reference Arkhangelsky and Imbens2022) and model-based approach to counterfactual outcome prediction. The latter is widely adopted in empirical studies, where scholars assume a parametric model for the potential outcomes. A general functional form for the observed outcome model is written as:

(10) $$ \begin{align} Y_{it} = g(D_{it}, X_{it}, \gamma_i, f_t, U_{it}, \varepsilon_{it}), \end{align} $$

where $\varepsilon _{it}$ is an error term with zero mean. In empirical studies, scholars add more constraints on the functional form to identify causal parameters despite the presence of unobserved confounding. For example, Pang, Liu, and Xu (Reference Pang, Liu and Xu2022) assume that the observed terms in $g(\cdot )$ are additive and that unobserved confounding can be represented by functions of $\gamma _i$ and $f_t$ , such as the TWFE and the latent factor term or IFE in the SC methods. With a statistical model for potential outcomes, scholars can predict the counterfactual outcomes under control for the treated observations and then average the predicted differences between the observed outcomes and predicted counterfactuals as an estimate for the ATT. Such counterfactual prediction-based approaches are called “causal panel data models” in Athey et al. (Reference Athey, Bayati, Doudchenko, Imbens and Khosravi2021).

Here, we focus on a linear model incorporating $\gamma _i$ and $f_t$ by way of a two-way interaction term or a latent factor term:

(11) $$ \begin{align} \begin{gathered} Y_{it} = \delta_{it} D_{it} + \beta_0 + X_{it}' \beta + \gamma_i^{\prime}f_t + \varepsilon_{it}, \\ \gamma_i | \boldsymbol{\Omega}_0 \sim \mathcal{N}(0, \boldsymbol{\Omega}_0), \quad f_t \sim \mathcal{N}(0, I_r), \quad \varepsilon_{it} \sim \mathcal{N}(0, \sigma^2), \end{gathered} \end{align} $$

where $\beta _0$ is the global intercept and $\beta = (\beta _1, \beta _2, \ldots , \beta _P) $ is the slope coefficients for $X_{it}$ . In this setup, $\gamma _i$ can be regarded as an ( $r \times 1$ ) vector of factor loadings and $f_t$ an ( $r \times 1$ ) vector of common factors. The number of factors r is unknown by the researcher a priori. The functional form above implies that $ Y_{it}(0) = \beta _0 + X_{it}' \beta + \gamma _i^{\prime }f_t + \varepsilon _{it} $ , i.e., we assume a latent factor model (LFM) for counterfactual outcomes under control. From a Bayesian perspective, Model (11) is a static variant of the state space model (Klinenberg Reference Klinenberg2022; Pang et al. Reference Pang, Liu and Xu2022). Following the literature of Bayesian factor analysis, we assume that the factor loadings are drawn from an i.i.d. multivariate normal distribution $\mathcal {N}(0, \boldsymbol {\Omega }_0)$ , where $ \boldsymbol {\Omega }_0 = Diag\{\omega _1^2, \ldots , \omega _r^2 \} $ is an ( $r \times r$ ) diagonal matrix, and that common factors are drawn from an i.i.d. multivariate normal distribution $\mathcal {N}(0, I_r)$ , where $I_r$ is the identity matrix. The error term $\varepsilon _{it}$ is drawn from i.i.d. normal distribution with mean 0 and variance $\sigma ^2$ . Model (11) entails the assumption that the panel latent ignorability (9) holds without conditioning on $U_{it}$ , i.e., unobserved confounding that varies both over units and time can be fully captured by the r-dimensional factor term $\gamma _i'f_t$ .

It is worth noting that there is an extensive literature on estimating Model (11) using both frequentist and Bayesian approaches.Footnote ² While the Bayesian approach offers greater flexibility in estimation, it relies on additional distributional assumptions for the factor term and the error term. In the frequentist approach, the factor term $\gamma _i{\prime }f_t$ , also known as “interactive fixed effects” (Bai Reference Bai2009), can be used to address endogeneity for the identification of coefficients $\beta $ , when some covariates in $X_{it}$ are endogenous. Additionally, instead of assuming that error terms are drawn from an independent standard normal distribution, they can be serially correlated within a unit.

In this article, we focus on BSA for treatment effects and thus assume independence between the factor term and the observed covariates, with i.i.d. error terms. However, Model (11) can be extended to incorporate endogenous covariates by directly modeling the relationship between $X_{it}$ and $\gamma _i{\prime }f_t$ . For the error terms, while we assume that $\varepsilon _{it}$ are i.i.d., the factor term, which includes unit random effects as a special case, provides a way to account for within-unit correlation (Cameron and Miller Reference Cameron and Miller2015). More complex structures for the error term, such as stochastic volatility (Belmonte, Koop, and Korobilis Reference Belmonte, Koop and Korobilis2014), can also be incorporated.

Determining the number of factors is a frequently discussed issue in the literature on estimating LFMs like (11). The prior variances for factor loadings, $\{ \omega _1^2, \ldots , \omega _r^2 \}$ , which control the effective number of factors, play a key role in the Bayesian analysis of factor models. For example, if $\omega _1^2 = 0$ , then the factor loading that corresponds to the first factor equals 0 for each unit, and thus it has no effect on the potential outcomes and is essentially excluded from the model. However, a well-known challenge in specifying prior for variance is that the conjugate inverse gamma distribution performs poorly for the purpose of selection of factor numbers. Since the support of the inverse gamma distribution includes only positive real numbers, even a conventional vague inverse gamma prior exerts a substantial influence on the posterior of $\omega ^2$ when its true value is close to 0 (Frühwirth-Schnatter and Wagner Reference Frühwirth-Schnatter and Wagner2010). There are two main approaches to this problem, both based on Bayesian shrinkage: spike-and-slab prior and continuous shrinkage prior. In this article, we adopt the continuous shrinkage approach. We follow Frühwirth-Schnatter and Wagner (Reference Frühwirth-Schnatter and Wagner2010) and re-parameterize the factor loadings as follows:

(12) $$ \begin{align} \begin{gathered} y_{it} = \beta_0 + \delta_{it} D_{it} + X_{it}^{\prime}\beta + (\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \varepsilon_{it}, \\ \tilde{\gamma}_i \sim \mathcal{N}(0, I_r), \quad f_t \sim \mathcal{N}(0, I_r), \quad \varepsilon_{it} \sim \mathcal{N}(0, \sigma^2), \end{gathered} \end{align} $$

where the factor loadings $ \gamma _i $ in model (11) are re-parameterized as the product of $ \omega _{\gamma } $ and $\tilde {\gamma }_i$ . In model (12), $\tilde {\gamma }_i$ is a vector of “standardized” factor loadings with each entry drawn from i.i.d. standard normal distribution and $\omega _\gamma $ is now interpreted as a vector of coefficients controlling the influence of each factor. If the $k^{th}$ element in $\omega _\gamma $ equals $0$ , then the $k^{th}$ factor is excluded from the model. The re-parameterization allows the use of a prior distribution for $\omega _\gamma $ that contains 0 in the interior of its support, enabling model selection by shrinkage.Footnote ³

2.4 Bayesian Sensitivity Analysis

Although the factor term in model (11) is highly flexible, it may not fully capture the unobserved confounding that varies across units and time. To address this common concern, we develop a sensitivity analysis method for treatment effect estimation under Assumption 3, where the existence of general confounding varying across units and time is permitted. As shown above in terms of a general setting (4) and (5), Assumption 3 enables imputation-based sensitivity analysis for causal effects by jointly modeling $Y_{it}$ conditional on ( $D_{it}$ , $X_{it}$ , $\gamma _i$ , $f_t$ , $U_{it}$ ) and $U_{it}$ conditional on ( $D_{it}, X_{it}, \gamma _i, f_t$ ). Specifically, we expand Model (12) to incorporate $U_{it}$ in the model:Footnote ⁴

(13) $$ \begin{align} \begin{gathered} Y_{it} = \beta_0 + \delta_{it} D_{it} + X_{it}^{\prime}\beta + \beta_u U_{it} + (\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \varepsilon_{it}, \\ U_{it} = \lambda_0 + \lambda_d D_{it} + X_{it}^\prime \lambda_x + \lambda_f (\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \nu_{it}, \\ \tilde{\gamma}_i \sim \mathcal{N}(0, I_r), \quad f_t \sim \mathcal{N}(0, I_r), \\ \varepsilon_{it} \sim \mathcal{N}(0, \sigma^2), \quad \nu_{it} \sim \mathcal{N}(0, c_1^2), \end{gathered} \end{align} $$

where the coefficient $\beta _u$ measures the effect of $U_{it}$ on the outcome of interest. We follow common practice (e.g., Gustafson et al. Reference Gustafson, McCandless, Levy and Richardson2010) to model $U_{it}$ as a linear function of the treatment indicator $D_{it}$ , observed confounders $X_{it}$ , as well as the factor terms $(\omega _{\gamma } \cdot \tilde {\gamma }_i)^{\prime } f_t$ , with their associated coefficients denoted by $\boldsymbol \lambda \equiv [\lambda _0, \lambda _d, \lambda _x^\prime , \lambda _f]^\prime $ . The variance of the normal error term $\nu _{it}$ equals $c_1^2$ , a fixed constant specified by the researcher. We call Model (13) “Bayesian sensitivity analysis for the linear factor model” (BSA-LFM).

$\boldsymbol \lambda $ and $\beta _u$ represent our sensitivity parameters and deserve several additional remarks. First, as discussed in the next section, we specify the prior for each element in $\boldsymbol \lambda $ to be symmetric and centered at zero, indicating that the grand mean of the unobserved confounder is $0$ a priori. Second, existing methods for sensitivity analysis often assume that the unobserved confounder is only determined by the treatment indicator to simplify the analysis (e.g., Lin, Psaty, and Kronmal Reference Lin, Psaty and Kronmal1998). However, the assumption that $U_{it}$ is independent of the observed covariates and the factor term conditional on treatment is often unrealistic and has been challenged by VanderWeele (Reference VanderWeele2008).Footnote ⁵ In addition, the common correlations of the outcome and the unobserved confounder imply additional learning (Gustafson et al. Reference Gustafson, McCandless, Levy and Richardson2010). That is, while $\beta _u$ and $\lambda _d$ are the only two sensitivity parameters that affect the estimation of the treatment effect,Footnote ⁶ allowing non-zero correlation of $U_{it}$ with observed confounders as well as the factor terms can help narrow the inference. If at least one element in $\lambda _x$ or $\lambda _f$ is not equal to 0, we can achieve substantive updating in $\beta _u$ given the data. Furthermore, allowing for correlation between observed covariates and unobserved confounders implies that some covariates may be endogenous (Diegert, Masten, and Poirier Reference Diegert, Masten and Poirier2022), which is a weaker condition than treating all covariates as exogenous.

Finally, it is important to note that we must assume some correlation between $\lambda _0$ and $\lambda _d$ . Since Model (13) allows for arbitrary heterogeneous treatment effects, like most counterfactual imputation methods, only observations under control ( $D_{it} = 0$ ) are used for model fitting. For these observations, Model (13) reduces to:

(14) $$ \begin{align} \begin{gathered} Y_{it|D_{it} = 0} = \beta_0 + X_{it}^{\prime}\beta + \beta_u U_{it} + (\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \varepsilon_{it}, \\ U_{it|D_{it} = 0} = \lambda_0 + X_{it}^\prime \lambda_x + \lambda_f (\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \nu_{it}. \end{gathered} \end{align} $$

Thus, we can only infer $\lambda _0$ directly from the data and indirectly estimate $\lambda _d$ through its assumed correlation with $\lambda _0$ . Details of the prior specifications are provided below.

2.5 Prior Specification

Prior specification plays a key role in BSA. Because data contain no direct information about unobserved confounders, posteriors of sensitivity parameters ( $\boldsymbol \lambda $ and $\beta _u$ ) are dominated by priors assigned to them. In fact, priors for the sensitivity parameters encode the researcher’s a priori belief about the direction and magnitude of unobserved confounding. Such beliefs are typically incorporated in frequentist sensitivity analysis only partially and informally by fixing the sensitivity parameters to particular values. In contrast, our Bayesian approach provides a more formal, principled way to represent researcher’s belief about confounding in the form of a probability distribution.

From a practical standpoint, the choice of a prior distribution for sensitivity parameters entails a bias–variance trade-off. An overly diffuse prior would result in a credible interval that has low power to detect a non-zero treatment effect, while a prior with high precision risks bias that never disappears asymptotically. Here, we adopt the approach proposed by McCandless, Gustafson, and Levy (Reference McCandless, Gustafson and Levy2008), which makes use of observed confounders to benchmark our prior on the effects of unobserved confounders. That is, we assume that the confounding effects of $U_{it}$ and $X_{it, j}$ ( $j \in \{1, \ldots , p\}$ ) are comparable a priori after standardization, in the sense that $\beta _u$ and $\beta $ have exchangeable prior distributions. Formally, we have

(15)

where the second line in (15) is by de Finetti’s theorem, and the first line adds several constraints on the prior distributions of sensitivity parameters $\boldsymbol \lambda $ as well as the variance of $\nu _{it}$ . We assign conjugate normal priors to $\boldsymbol {\lambda }$ with constraints on prior variancesFootnote ⁷ and exchangeable priors to $\beta _u$ and $\beta $ . More specifically, our prior on $\boldsymbol {\lambda }$ is specified as follows:

• • $\lambda _0$ and $\lambda _d$ :

  (16) $$ \begin{align} \begin{gathered} (\lambda_0, \lambda_d)' \sim \mathcal{N}(0, \Lambda), \\ \Lambda = \begin{pmatrix} c_2^2, -c_3^2/2\\ -c_3^2/2, c_3^2 \end{pmatrix}, \end{gathered} \end{align} $$
  where the hyperparameters are constrained such that $0<c_2\leq 1$ , $c_3> 0,$ and $ 2c_2> c_3 $ to guarantee positive definiteness of the covariance matrix.

• • $\lambda _x$ :

  (17) $$ \begin{align} \lambda_x \sim \mathcal{N}\left(0, \frac{c_4^2}{p} I_p\right). \end{align} $$

• • $\lambda _f$ :

  (18) $$ \begin{align} \lambda_f \sim \mathcal{N} \left(0, \frac{c_5^2(k_1-1)}{2rk_2}\right), \end{align} $$

where $k_1$ and $k_2$ are hyperparameters for $\omega _\gamma $ defined below. Taken together, our prior specification for $\boldsymbol {\lambda }$ implies the following proposition.Footnote ⁸

Proposition 1 Prior variance of unobserved confounder

Priors (16)–(18) imply that the prior variance of $U_{it}$ is independent of treatment status $D_{it}$ and given by,

(19) $$ \begin{align} Var(U_{it}|D_{it}) = Var(U_{it}) = c_1^2 + c_2^2 + c_4^2 + c_5^2. \end{align} $$

A formal proof of Proposition 1 can be found in Appendix 1 of the Supplementary Material. Since we standardize the observed and unobserved confounders to have unit variance, Proposition 1 implies the following constraint on the hyperparameters:

(20) $$ \begin{align} c_1^2 + c_2^2 + c_4^2 + c_5^2 = 1. \end{align} $$

Now, we turn to parameters in the outcome model and prior variances of factor loadings. There are several sensible options for exchangeable priors assigned to $\beta $ and $\beta _u$ . Of particular note, while an independent, weakly informative normal prior is perhaps the most intuitive choice, Gustafson et al. (Reference Gustafson, McCandless, Levy and Richardson2010) use the i.i.d. t distribution with a certain degree of freedom, noting that the latter often yields substantially narrower inferences than a normal independence prior with the same marginal variability. This prior has a hierarchical form: the priors for $ \beta _0, \beta _1, \ldots , \beta _p, \beta _u $ are i.i.d. $N(0, \tau ^2 S)$ given S, and the hyperparameter S is drawn from inverted-gamma distribution $S \sim Inv-Gamma(d/2, d/2)$ .

Here, we propose a more flexible specification, allowing each $\beta _k$ and $\beta _u$ to have heterogeneous prior variances drawn from a multi-level model. Specifically, we use the Laplacian prior (Park and Casella Reference Park and Casella2008), which is a continuous shrinkage prior that have heavier tails than the t distribution and more densities near $0$ to “let the data speak for itself.”

• • $ \beta _j $ for $ j \in \{0, 1, \ldots , p \} $ :Footnote ⁹

  (21) $$ \begin{align} \begin{gathered} \beta_j | \tau_{\beta_j}^2 \sim \mathcal{N}(0, \tau_{\beta_j}^2), \\ \tau_{\beta_j}^2 | \xi_{\beta} \sim Exp(\frac{\xi_{\beta}^2}{2}), \\ \xi_{\beta}^2 \sim Gamma(a_1, a_2). \end{gathered} \end{align} $$

• • $ \beta _u $ :

  (22) $$ \begin{align} \begin{gathered} \beta_u | \tau_{\beta_u}^2 \sim \mathcal{N}(0, \tau_{\beta_u}^2), \\ \tau_{\beta_u}^2 | \xi_{\beta} \sim Exp(\frac{\xi_{\beta}^2}{2}), \\ \xi_{\beta}^2 \sim Gamma(a_1, a_2). \end{gathered} \end{align} $$

We also assign the Laplacian prior to $\omega _\gamma $ for the purpose of factor number selection.

• • $ \omega _{\gamma _j} $ for $ j \in \{1, \ldots , r \} $ :

  (23) $$ \begin{align} \begin{gathered} \omega_{\gamma_j} | \tau_{\gamma_j}^2 \sim N(0, \tau_{\gamma_j}^2), \\ \tau_{\gamma_j}^2 | \xi_{\gamma} \sim Exp(\frac{\xi_{\gamma}^2}{2}), \\ \xi_{\gamma}^2 \sim Gamma(k_1, k_2). \end{gathered} \end{align} $$

Finally, we use the conventional inverted-gamma prior for the error variance $\sigma ^2$ .

• • $ \sigma _\varepsilon ^2 $ :

  (24) $$ \begin{align} \sigma^{2} \sim Inv-Gamma(e_1, e_2). \end{align} $$

Taken together, our model for sensitivity analysis requires 11 hyperparameter values to be set by the researcher, subject to the constraints as specified above. Appropriately setting hyperparameter values is often a challenging task in implementing Bayesian methods in practice, particularly for complex models. Fortunately, sensible default values that approximate a lack of prior knowledge are available for most of the hyperparameters in our model, including those related to the error variance, coefficients in the outcome model, and factor selection. This substantially simplifies the task of specifying priors in applied research.

Table 1 summarizes such default values for these hyperparameters, along with their substantive interpretations. For the prior variances of the sensitivity parameters $\lambda $ s, we impose constraints to keep their scales small while allowing $\beta _u$ to be potentially large. In practice, selecting these values may require some prior knowledge. For instance, if researchers believe that $U_{it}$ is informative given $X_{it}$ or the factor terms, they may choose relatively large values for $c_4^2$ or $c_5^2$ . Otherwise, they may favor a larger random component, reflected in a larger $c_1^2$ .

Table 1 Suggested default specifications for the hyperparameters.

3.1 Transparent Parameterization of the Nonidentified Model

We begin by briefly summarizing the general idea of transparent parameterization proposed by Gustafson (Reference Gustafson2005). Suppose the original parameter vector in the nonidentified model is $\theta $ . we reparameterize $\theta $ to new parameter vector $\phi \equiv \{ \phi _I, \phi _N \} $ such that $ f(\text {data}|\phi ) = f(\text {data}|\phi _I) $ . Priors assigned to $\phi $ can be factorized as $ f(\phi ) = f(\phi _N|\phi _I)f(\phi _I) $ , and the corresponding posteriors are:

(25) $$ \begin{align} \begin{gathered} f(\phi_I|\text{data}) \propto f(\text{data}|\phi_I) f(\phi_I), \\ f(\phi_N|\text{data}) = \int f(\phi_N|\phi_I) f(\phi_I|\text{data}) d\phi_I. \end{gathered} \end{align} $$

The posteriors (25) imply qualitatively different asymptotic behaviors for $\phi _I$ and $\phi _N$ . That is, direct Bayesian learning occurs for $\phi _I$ , in that under regularity conditions, $f(\phi _I|\text {data})$ converges to the point probability mass at the true value of $\phi _I$ as sample size grows, and indirect Bayesian learning for $\phi _N$ , which converges to the conditional distribution $f(\phi _N|\phi _I)$ at the true value of $\phi _I$ . This, in turn, implies that $\phi _I$ and $\phi _N$ will also exhibit distinct behaviors in an MCMC sampler.

Gustafson et al. (Reference Gustafson, McCandless, Levy and Richardson2010) suggest that transparent parameterization could be exploited for the purpose of building a computationally efficient posterior sampler, making jumps in spaces for those parameters informed by the data while moving fast through spaces for those parameters with little information from the data. We adopt a similar strategy for our proposed BSA-LFM. Specifically, we substitute $U_{it}$ in the first line for potential outcome with the second line in model (13) that describe the DGP of $U_{it}$ , essentially integrating out the unobserved confounder:

(26) $$ \begin{align} Y_{it} = (\delta_{it} + \beta_u \lambda_d) D_{it} + (\beta_0 + \beta_u \lambda_0) + X_{it}^{\prime}(\beta + \beta_u \lambda_x) + (1 + \beta_u \lambda_f)(\omega_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \beta_u \nu_{it} + \varepsilon_{it}. \end{align} $$

Define $\tilde {\delta }_{it} = \delta _{it} + \beta _u \lambda _d$ , $\tilde {\beta }_0 = \beta _0 + \beta _u \lambda _0$ , $\tilde {\beta } = \beta + \beta _u \lambda _x$ , $\tilde {\omega }_{\gamma } = (1 + \beta _u \lambda _f)\omega _{\gamma }$ , and $\tilde {\epsilon }_{it} = \beta _u \nu _{it} + \varepsilon _{it}$ . Model (13) can then be re-written as:

(27) $$ \begin{align} \begin{gathered} Y_{it} = \tilde{\delta}_{it} D_{it} + \tilde{\beta}_0 + X_{it}^{\prime}\tilde{\beta}+ (\tilde{\omega}_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t + \tilde{\epsilon}_{it}, \\ \tilde{\epsilon}_{it} \sim \mathcal{N} (0, \tilde{\sigma}^2), \\ \tilde{\sigma}^2 = \beta_u^2 c_1^2 + \sigma_Y^2. \end{gathered} \end{align} $$

The original parameter vector is $\theta = \{ \beta _0, \beta , \boldsymbol {\omega }, \Gamma , F, \sigma _Y^2, \beta _u, \lambda _0, \lambda _d, \lambda _x, \lambda _f \}$ , and the re-parameterized vector is:

(28) $$ \begin{align} \phi = \left\{ \tilde{\beta}_0, \tilde{\beta}, \tilde{\boldsymbol{\omega}}, \Gamma, F, \tilde{\sigma}^2, \beta_u, \lambda_0, \lambda_d, \lambda_x, \lambda_f \right\}, \end{align} $$

where $\phi _I = \left \{\tilde {\beta }_0, \tilde {\beta }, \tilde {\boldsymbol {\omega }}, \Gamma , F, \tilde {\sigma }^2\right \}$ and $\phi _N = \left \{\beta _u, \lambda _0, \lambda _d, \lambda _x, \lambda _f\right \}$ . One can verify that the data likelihood only depends on $ \phi _I $ :Footnote ¹⁰

(29) $$ \begin{align} f(Y^{obs} | \phi ) &\propto (\sqrt{\sigma^2 + \beta_u^2 c_1^2}) ^{-(NT - \sum_{i,t} D_{it})} \nonumber\\& \quad \exp [-\frac{1}{2(\sigma^2 + \beta_u^2 c_1^2)} \sum_{D_{it} = 0} (y_{it} - (\beta_0 + \beta_u \gamma_0) - X_{it}^{\prime}(\beta + \beta_u \gamma_x) - (\omega_{\gamma} (1 + \beta_u \gamma_f) \cdot \tilde{\gamma}_i)^{\prime} f_t)^2 ] \nonumber\\&= \tilde{\sigma}^{-(NT - \sum_{i,t} D_{it})} \exp [-\frac{1}{2\tilde{\sigma}^2} \sum_{D_{it} = 0} (y_{it} - \tilde{\beta}_0 - X_{it}^{\prime}\tilde{\beta} - (\tilde{\omega}_{\gamma} \cdot \tilde{\gamma}_i)^{\prime} f_t)^2 ]\nonumber\\&\propto f(Y^{obs} | \phi_I ). \\[-40pt]\nonumber \end{align} $$

Now, we derive the posteriors for parameters in the re-parameterized model (27). Denote the one-to-one mapping of the transparent re-parameterization by $ \phi = g(\theta ) $ . Prior distributions assigned to the original parameter vector $\theta $ imply priors assigned to $\phi $ are:

(30) $$ \begin{align} f_\phi(\phi) = f_\theta(g^{-1}(\phi))|J|, \end{align} $$

where $|J|$ is the determinant of the Jacobian $ \theta = g^{-1}(\phi ) $ . Thus, posterior distributions of $\phi $ can be written as:

(31) $$ \begin{align} \begin{aligned} \pi(\phi | Y^{obs}) &\propto f(Y^{obs} | \phi ) f_\phi(\phi) \\& = f(Y^{obs} | \phi ) f_\theta(g^{-1}(\phi))|J| \\& = f(Y^{obs} | \phi ) |1 + \beta_u\lambda_f|^{-r} \\&\quad \prod_{j = 0}^{p} \pi(\tilde{\beta}_j - \beta_u \lambda_x | \tau_{\beta_j}^2) \pi(\tau_{\beta_j}^2|\xi_\beta^2) \prod_{k = 1}^{r}\pi(\frac{\tilde{\omega}_{\gamma_j}}{1 + \beta_u\lambda_f} | \tau_{\gamma_j}^2) \pi(\tau_{\gamma_j}^2|\xi_\gamma^2) \pi(\xi_\gamma^2) \\&\quad \pi(\beta_u | \tau_{\beta_u}^2) \pi(\tau_{\beta_u}^2|\xi_\beta^2) \pi(\tilde{\Gamma})\pi(F) \pi(\tilde{\sigma}^{2} - \beta_u^2 c_1^2) \pi(\lambda_0, \lambda_d) \pi(\lambda_x) \pi(\lambda_f), \end{aligned} \end{align} $$

where $\pi (\cdot )$ denotes the prior distribution of the corresponding original parameter. Then, we can recover the posteriors of the original parameter vector $\theta $ from the posteriors (31) of $\phi $ . Since prior distributions are assigned to $\theta $ , the implied priors are not conjugate for some parameters in $\phi $ . Therefore, we develop a Metropolis–Hastings step within Gibbs sampling procedure to simulate the posteriors summarized as Algorithm 1. Details of the MCMC algorithm can be found in Appendix 2 of the Supplementary Material.

Remark. In Step 10 of Algorithm 1, we estimate the individual treatment effect for observations under treatment ( $D_{it} = 1$ ) using $\delta _{it} = \tilde {\delta }_{it} - \beta _u \cdot \lambda _d$ , where $\tilde {\delta }_{it} = Y_{it} - Y_{it}(0)$ , with $Y_{it}(0)$ representing the imputed counterfactual outcome under control. This step demonstrates the “frequentist analogy” of the proposed BSA. To see this, suppose the treatment effects are homogeneous ( $\delta _{it} \equiv \delta $ ). In that case, the identified parameter is $\tilde {\delta } = \delta + \beta _u \cdot \lambda _d$ . Researchers can obtain an interval estimate for $\delta $ by fixing the sensitivity parameters at specific values and varying these values. The BSA approach, rather than selecting specific values for the sensitivity parameters, marginalizes the treatment effect estimate over the posterior distributions of $\beta _u$ and $\lambda _d$ .