Bayesian semi-parametric G-computation for causal inference in a cohort study with MNAR dropout and death

TL;DR

This paper introduces a Bayesian semi-parametric G-computation method using Bayesian additive regression trees to estimate causal effects in longitudinal studies with MNAR dropout and death, addressing limitations of traditional parametric approaches.

Contribution

It develops a flexible Bayesian approach with sensitivity parameters for causal inference under MNAR, applicable to complex longitudinal observational data.

Findings

Successfully applied to a cohort study on cognition and aging.

Provided estimates of widowhood effect on memory accounting for MNAR dropout.

Compared favorably to standard methods in simulation and real data analysis.

Abstract

Causal inference with observational longitudinal data and time-varying exposures is often complicated by time-dependent confounding and attrition. The G-computation formula is one approach for estimating a causal effect in this setting. The parametric modeling approach typically used in practice relies on strong modeling assumptions for valid inference, and moreover depends on an assumption of missing at random, which is not appropriate when the missingness is missing not at random (MNAR) or due to death. In this work we develop a flexible Bayesian semi-parametric G-computation approach for assessing the causal effect on the subpopulation that would survive irrespective of exposure, in a setting with MNAR dropout. The approach is to specify models for the observed data using Bayesian additive regression trees, and then use assumptions with embedded sensitivity parameters to identify and estimate the causal effect. The proposed approach is motivated by a longitudinal cohort study on cognition, health, and aging, and we apply our approach to study the effect of becoming a widow on memory. We also compare our approach to several standard methods.

Tables

Table 1. Table 2: Algorithm for estimation of τ 𝜏 \tau in [ 2 ] using the G-computation formula. Details of the algorithm can be found in the Web Appendix section A.4.

1.	Sample the observed data posteriors as described in Section 5.
2.	For each posterior sample of the parameters sample pseudo data $({\overline{y}}_{j-1}^{\ast },{\overline{w}}_j^{\ast },{\overline{r}}_j^{\ast },{\overline{s}}_j^{\ast },x_0^{\ast })$
	and sensitivity parameters ${\gamma }_j$, ${\nu }_j$, $c\left(z_j\right)$, and $c\left(z_j^{\prime }\right)$ of size $N^{\ast }$. Additionally, sample
	one set of ${\Delta }_{{\overline{z}}_j^{\prime }}$.
3.	Implement G-computation for ${\overline{z}}_j$, and similarly for ${\overline{z}}_j^{\prime }$, using the pseudo data and
	sensitivity parameters from Step 2 by computing $E\left[Y_j\mid {\overline{y}}_{j-1},{\overline{z}}_j,{\overline{w}}_j,{\overline{r}}_j,{\overline{s}}_j=1,x_0\right]$
	and ${\prod }_{k=0}^j\mathrm{Pr}[S_k=1\mid {\overline{z}}_k,{\overline{w}}_k,{\overline{r}}_k,{\overline{y}}_{k-1},{\overline{S}}_{k-1}=1,x_0]$. Furthermore, implement Monte
	Carlo integration using the pseudo data to compute $\mathrm{Pr}\left[{\overline{S}}_j=1\mid {\overline{z}}_j\right]$ and
	$E\left[Y_j,{\overline{S}}_j=1\mid {\overline{z}}_j\right]$.
4.	Use the quantities in Step 3 to compute one posterior sample of $\tau$ as defined in [Result 1:]-[4].
5	Repeat step 2 - 4 for each of the posterior sample of the parameters.

Table 2. Table 3: Simulation results for causal effect estimation with n=1000 with the true causal effect τ = − 0.05 𝜏 0.05 \tau=-0.05 , using two settings for our proposed approach: error specified as normal (BSP-GC1) and error specified using a t-distribution (BSP-GC2), a parametric version of the proposed procedure (BP-GC), inverse probability of treatment weights using unstabilized weights (IPTW-W) and stabilized weights (IPTW-SW), and Targeted minimum loss-based estimation approach for longitudinal data structures (TMLE). Mean squared error (MSE) are multiplied by 100 for ease of presentation. ESD denotes empirical standard deviation and CP denotes coverage probability of 95% credible intervals.

	Bias	ESD	MSE	CP
BSP-GC1	-0.002	0.013	0.02	98.7
BSP-GC2	-0.003	0.013	0.02	98.4
BP-GC	-0.065	0.021	0.47	63.6
IPTW-W	0.024	0.040	0.22	64.3
IPTW-SW	-0.031	0.013	0.12	20.4
TMLE	-0.002	0.021	0.04	92.8

Table 3. Table 4: Simulation results for causal effect estimation with n=1000 with the true causal effect τ = − 0.05 𝜏 0.05 \tau=-0.05 . Two scenarios, a) where the error distribution for the outcome is misspecified using a t-distribution with 3 degrees of freedom and b) a setting with limited overlap. Comparing our proposed approach (BSP-GC1), a parametric version of the proposed procedure (BP-GC), inverse probability of treatment weights using unstabilized weights (IPTW-W) and stabilized weights (IPTW-SW), and Targeted minimum loss-based estimation approach for longitudinal data structures (TMLE). Mean squared error (MSE) are multiplied by 100 for ease of presentation. ESD denotes empirical standard deviation and CP denotes coverage probability of 95% credible intervals.

	a) $t_3$				b) non-overlap
	Bias	ESD	MSE	CP	Bias	ESD	MSE	CP
BSP-GC1	-0.002	0.037	0.13	97.2	-0.003	0.014	0.02	98.9
BP-GC	-0.074	0.038	0.69	69.5	-0.078	0.024	0.67	43.0
IPTW-W	0.028	0.068	0.55	59.7	-0.026	0.033	0.18	56.2
IPTW-SW	-0.031	0.022	0.15	43.5	-0.042	0.015	0.20	10.6
TMLE	-0.001	0.036	0.13	92.3	0.002	0.023	0.05	93.9

Table 4. Table 5: Comparison of methods used for causal effect estimation of the Betula data, setting Δ z ¯ j ′ = 0 subscript Δ subscript superscript ¯ 𝑧 ′ 𝑗 0 \Delta_{\bar{z}^{\prime}_{j}}=0 , γ j = 0 subscript 𝛾 𝑗 0 \gamma_{j}=0 , and c ​ ( z j ) = c ​ ( z j ′ ) = 0 𝑐 subscript 𝑧 𝑗 𝑐 subscript superscript 𝑧 ′ 𝑗 0 c(z_{j})=c(z^{\prime}_{j})=0 , using our proposed approach (BSP-GC), a parametric version of the proposed procedure (BP-GC), inverse probability of treatment weights using unstabilized weights (IPTW-W) and stabilized weights (IPTW-SW), and Targeted minimum loss-based estimation approach for longitudinal data structures (TMLE).

	Estimate [95% CI]
BSP-GC	-0.98 [-2.78, 0.73]
BP-GC	-0.53 [-1.73, 0.68]
IPTW-W	-1.67 [-5.96, 1.51]
IPTW-SW	-0.44 [-3.06, 1.39]
TMLE	-0.96 [-3.11, 0.99]

Equations

$$\mathrm{E}[Y_j({\bar{z}}_j)-Y_j({\bar{z}}_j^{\prime })\mid {\bar{S}}_j({\bar{z}}_j)={\bar{S}}_j({\bar{z}}_j^{\prime })=1].$$

$$\tau =\frac{{\sum }_{j=1}^J\mathrm{E}[Y_j({\bar{z}}_j)-Y_j({\bar{z}}_j^{\prime })\mid {\bar{S}}_j({\bar{z}}_j)={\bar{S}}_j({\bar{z}}_j^{\prime })=1]\times \mathrm{Pr}[{\bar{S}}_j({\bar{z}}_j)={\bar{S}}_j({\bar{z}}_j^{\prime })=1]}{{\sum }_{k=1}^J\mathrm{Pr}[{\bar{S}}_k({\bar{z}}_k)={\bar{S}}_k({\bar{z}}_k^{\prime })=1]}.$$

$$Y_j({\bar{z}}_j)\perp \perp Z_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_j,{\bar{r}}_j,{\bar{s}}_j=1,x_0$$

$$S_j({\bar{z}}_j)\perp \perp Z_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_j,{\bar{r}}_j,{\bar{s}}_{j-1}=1,x_0.$$

$$p(y_j\mid {\bar{y}}_{j-1},{\bar{z}}_j,{\bar{w}}_j,r_j=0,{\bar{s}}_j=1,x_0)=p(y_j\mid {\bar{y}}_{j-1},{\bar{z}}_j,{\bar{w}}_j,{\bar{r}}_j=1,{\bar{s}}_j=1,x_0)$$

$$\mathrm{Pr}[z_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_j,r_j=0,{\bar{r}}_{j-1},s_j=0,{\bar{s}}_{j-1}=1,x_0]$$

$$=\mathrm{Pr}[z_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_j,{\bar{r}}_j,{\bar{s}}_j=1,x_0],$$

$$Pr[w_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_{j-1},r_j=0,{\bar{r}}_{j-1},s_j=0,{\bar{s}}_{j-1}=1,x_0]$$

$$=\mathrm{Pr}[w_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_{j-1},{\bar{r}}_j,{\bar{s}}_j=1,x_0],$$

$${\nu }_j=$$

$$\mathrm{Pr}[z_j\mid {\bar{y}}_{j-1},{\bar{z}}_{j-1},{\bar{w}}_j,{\bar{r}}_j,{\bar{s}}_j=1,x_0],$$

$$E[Y_j({\bar{z}}_j)-Y_j({\bar{z}}_j^{\prime })\mid \bar{S}({\bar{z}}_j)={\bar{S}}_j({\bar{z}}_j^{\prime })=1]=$$

$$\frac{E_{\mathcal{A}}[E(Y_j,{\bar{S}}_j=1\mid {\bar{z}}_j,c({\bar{z}}_j),\overline{{\gamma }_j},\mathcal{A})]}{E_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_j=1\mid {\bar{z}}_j,c({\bar{z}}_j),\overline{{\gamma }_j},{\nu }_j,\mathcal{A})]}-\frac{E_{\mathcal{A}}[E(Y_j,{\bar{S}}_j=1\mid {\bar{z}}_j^{\prime },c({\bar{z}}_j^{\prime }),\overline{{\gamma }_j},\mathcal{A})]}{E_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_j=1\mid {\bar{z}}_j^{\prime },c({\bar{z}}_j^{\prime }),\overline{{\gamma }_j},{\nu }_j,\mathcal{A})]}-$$

$${\Delta }_{{\bar{z}}_j^{\prime }}\times \left(1-\frac{E_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_j=1\mid {\bar{z}}_j,c({\bar{z}}_j),\overline{{\gamma }_j},{\nu }_j,\mathcal{A})]}{E_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_j=1\mid {\bar{z}}_j^{\prime },c({\bar{z}}_j^{\prime }),\overline{{\gamma }_j},{\nu }_j,\mathcal{A})]}\right),$$

$$\frac{\mathrm{Pr}[{\bar{S}}_j({\bar{z}}_j)={\bar{S}}_j({\bar{z}}_j^{\prime })=1]}{{\sum }_{k=1}^J\mathrm{Pr}[{\bar{S}}_k({\bar{z}}_k)={\bar{S}}_k({\bar{z}}_k^{\prime })=1]}=\frac{E_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_j=1\mid {\bar{z}}_j,c({\bar{z}}_j),\overline{{\gamma }_j},{\nu }_j,\mathcal{A})]}{{\sum }_{k=1}^J{E}_{\mathcal{A}}[\mathrm{Pr}({\bar{S}}_k=1\mid {\bar{z}}_k,c({\bar{z}}_k),\overline{{\gamma }_k},{\nu }_k,\mathcal{A})]}.$$

Full text

Bayesian semi-parametric G-computation for causal inference in a cohort study with MNAR dropout and death.

Josefsson

Centre for Demographic and Ageing Research,

Umeå University, Sweden.

Corresponding author: Maria Josefsson,

Umeå University, SE-901 87 Umeå, Sweden.

E-mail: maria.josefsson@umu.se

Maria

Centre for Demographic and Ageing Research,

Umeå University, Sweden.

Corresponding author: Maria Josefsson,

Umeå University, SE-901 87 Umeå, Sweden.

E-mail: maria.josefsson@umu.se

Daniels

Department of Statistics, University of Florida.

Michael J

Department of Statistics, University of Florida.

Abstract

Causal inference with observational longitudinal data and time-varying exposures is often complicated by time-dependent confounding and attrition. The G-computation formula is one approach for estimating a causal effect in this setting. The parametric modeling approach typically used in practice relies on strong modeling assumptions for valid inference, and moreover depends on an assumption of missing at random, which is not appropriate when the missingness is missing not at random (MNAR) or due to death. In this work we develop a flexible Bayesian semi-parametric G-computation approach for assessing the causal effect on the subpopulation that would survive irrespective of exposure, in a setting with MNAR dropout. The approach is to specify models for the observed data using Bayesian additive regression trees, and then use assumptions with embedded sensitivity parameters to identify and estimate the causal effect. The proposed approach is motivated by a longitudinal cohort study on cognition, health, and aging, and we apply our approach to study the effect of becoming a widow on memory. We also compare our approach to several standard methods.

Keywords

BART, Cognitive aging, Longitudinal data, Observational data, Non-ignorable missing, Sensitivity analysis, Survivor Average Causal Effect, Time-varying exposure, Time-varying confounding.

1 Introduction

Causal inference in non-randomized longitudinal studies with time-varying exposures is often complicated by time-dependent confounding and attrition. Attrition is inevitable especially if individuals in the studied population are older and followed over a long time period. Additionally, for cohort studies, an individual’s data is only recorded if that person completes follow-up testing. Hence, data for not only the outcome but also exposure level and confounders are missing at subsequent test waves.

The G-computation formula ([33]) is one approach for estimating a causal effect of time-varying exposures when time-varying confounding is present. The approach is completely nonparametric in its original form, although a parametric modeling approach based on maximum likelihood estimation is most typically used in practice (e.g. [39, 43]). Valid inference with the parametric G-formula requires correct model specification. This can be extremely difficult when there is a large set of regressors, the relationship is non-linear and/or includes interaction terms, and there are multiple observation times. Non- and semi-parametric estimation techniques that do not require prespecified distributional or functional forms of the data, have become popular in the causal inference literature (e.g. [13, 11, 19, 17, 40, 42]). One such modeling strategy is Bayesian Additive Regression Trees (BART, [3]). BART is a sum-of-trees model that adds together the predictions of a number of regression trees regularized by prior distributions. BART does not rely on strong modeling assumptions, and in contrast to other tree-based algorithms BART yields interval estimates for full posterior inference.

A number of methodologies have been applied to address missing response or missing covariate data in causal effect estimation of longitudinal data under an assumption of missing at random (MAR; [2, 34]). These methods, however, are generally invalid when the missingness is missing not at random (MNAR) or due to death ([21]). Partly conditional models have been proposed to address the combination of dropout and truncation by death, where inference is conditioning on the sub-population being alive at a specific time-point ([20, 38, 26, 32, 46]) However, conditioning on survival may introduce bias due to the fact that survival is a post-randomization event. One estimand that has gained much attention to address this issue is the ”survivors average causal effect” (SACE), i.e. the causal effect on the subpopulation of those surviving irrespective of exposure ([7, 8]). Several approaches have been developed for estimation of the SACE in longitudinal randomized control studies (e.g. [23, 24, 44, 45]), or in context of semicompeting risks ([6, 48]). For observational data [41] (\citeyeartchetgen2014identification) developed a weighting estimator to identify the SACE without missingness, and [37] (\citeyearshardell2014doubly) identified the SACE with MAR missingness using also a weighting technique. Moreover, [15] (\citeyearjosefsson2016causal) proposed assumptions to identify the SACE of a baseline exposure on a longitudinal outcome under MNAR missingness for the outcome using parametric methods. These approaches however, do not appropriately account for MNAR data among survivors when the exposure and confounding are time-varying. [36] (\citeyearshardell2018joint) proposed a parametric shared parameter model with g-computation to identify a principal stratum causal effect for observational longitudinal data with time-dependent confounding. A drawback of their approach is that unbiased estimation depends on correct model specification and it does not appropriately account for MNAR data among survivors.

Widowhood has been identified as an important social factor associated with increased mortality ([12]) and cognitive impairment (e.g. [28]). Here, our goal is to develop a framework for assessing the impact of becoming a widow on memory, a monotone exposure, by estimating the SACE in a setting with MNAR dropout among survivors. The proposed approach is motivated by the Betula study ([29]), where individuals are followed over multiple test waves to study how cognitive functions potentially deteriorate with age and identify risk factors for dementia.

The remainder of the paper is organized as follows. In Section 2, we introduce the notation and the causal estimand. In Section 3, we propose identifying default assumptions and sensitivity parameters to allow deviations from these assumptions, followed by the identification of the SACE in Section 4. In Section 5, we propose a Bayesian semi-parametric (BSP) modeling approach for the observed data distributions and the algorithm for estimation of the SACE. In Section 6, we provide a simulation and in Section 7 an application to the Betula data. Finally, we conclude with a discussion and possible future work in Section 8.

2 Notation and the causal effect of interest

2.1 Data structure and notation

We begin with a formal description of the data. Let $i=1,2,\ldots,N$ denote individual and $j=0,1,\ldots,J$ denote time (the data used from the Betula study has $J=3$ follow-up test waves). We denote the vector of baseline confounders by $X_{i0}$ (gender, education, and age cohort) and the time-varying confounder by $W_{ij}$ (if the spouse has been seriously ill between the $j-1$th and $j$th test wave). The continuous memory outcome is denoted by $Y_{ij}$ and the binary exposure is denoted by $Z_{ij}$. We assume a monotone exposure where initially all subjects are unexposed ($Z_{10}=0$ for all $i$). If a subject is exposed (widowed) at test wave $j$ $Z_{ij}=1$ and if $Z_{ij}=1$, then $Z_{ik}=1$ for $k>j$. Let $S_{ij}$ denote survival, where $S_{ij}=1$ if an individual is alive at the time of the testing and 0 otherwise. Let $R_{ij}$ be a dropout indicator, where $R_{ij}=1$ if an individual has completed the cognitive testing or 0 otherwise. We have monotone missingness, so if $R_{ij}=0$, $R_{ik}=0$ for $k>j$. Note that vital status information is presumed to be available even after dropping out of the study. The history of the time-varying variables are denoted with an overbar. For example, the exposure history for individual $i$ through test wave $j$ is denoted by $\bar{Z}_{ij}=\{Z_{i0},Z_{i1},\ldots,Z_{ij}\}$. Furthermore, for individual $i$, $J^{r}_{i}$ denotes the number of test waves (s)he participates in the study, and $J^{s}_{i}\geq J^{r}_{i}$ denotes the number of test waves (s)he is alive. A simplified version of the study design restricted to two test waves is depicted in a causal diagram in Figure 1.

2.2 Causal estimand

The goal of the study is to estimate the causal effect of becoming a widow (within 5 years) on memory among those who would survive irrespective of being widowed or not. We consider two contrasting exposure regimes, $\bar{z}_{ij}=\{z_{i0}=0,\ldots,z_{ij-1}=0,z_{ij}=1\}$, i.e individuals exposed (widowed) between the $j-1$th and $j$th wave, and the contrasting regime $\bar{z}^{\prime}_{ij}=\{z_{i0}=0,\ldots,z_{ij}=0\}$, i.e. individuals unexposed through test wave $j$, for $j=1,2,3$. Below, we generally suppress the subscript $i$ to simplify notation. The potential memory outcome at wave $j$ is denoted by $Y_{j}(\bar{z}_{j})$ for an individual under exposure regime $\bar{z}_{j}$. Similarly, let $S_{j}(\bar{z}_{j})$ be the potential survival outcome at wave $j$, denoting survival under exposure regime $\bar{z}_{j}$.

We consider a principal stratum causal effect of a time-varying exposure on the outcome, at wave $j$, for those who would survive under either exposure regime,

[TABLE]

However, main interest is not the effect at a specific wave, but rather the effect aggregated over test waves, defined as

[TABLE]

3 Identifying assumptions and sensitivity parameters

To identify the causal effect in [2] from the observed data we first introduce a set of assumptions followed by a set of sensitivity parameters to assess the impact of violations to some of the assumptions. The sensitivity parameters (and their values) will be explained in relation to the Betula data in Section 7.2.

3.1 Assumptions

Assumptions $1-3$ are a set of standard assumptions for causal inference of longitudinal observational data:

Assumption 1 Consistency: For a given individual, if $\bar{Z}_{j}=\bar{z}_{j}$, then $Y_{j}=Y_{j}(\bar{z}_{j})$ and $S_{j}=S_{j}(\bar{z}_{j})$.

Assumption 2 Positivity for a monotone exposure: $\Pr[z_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1}=0,\bar{w}_{j},\bar{r}_{j-1},\bar{s}_{j-1}=1,x_{0}]>0$ for $z_{j}=0,1$ and for all individuals, such that all unexposed individuals have a nonzero probability of becoming exposed between test wave $j-1$ and $j$ if $p(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},\bar{r}_{j-1},\bar{s}_{j-1}=1,x_{0})\neq 0$.

Assumption 3 Conditional exchangeability: If $X_{0}$ and \vbox{\hrule height=0.03pt\kern 1.50696pt\hbox{\kern-1.99997ptW\kern-1.99997pt}}_{j} contains all pre-exposure covariates related to exposure, potential outcomes and survival, then for all exposure regimes

[TABLE]

That is, at each test wave $j$, being exposed $z_{j}$ is as if randomized conditional on the set of the temporally preceding variables. The assumption of conditional exchangeability is likely to be violated in many settings and is impossible to assess from the data. Therefore, we introduce a sensitivity parameter to investigate sensitivity for unmeasured confounding in Section 3.2.

In cohort studies $Y_{j},Z_{j}$ and $W_{j}$ are not observed (but defined) for individuals who are alive but who drop out of the study. We make an MAR type assumption conditional on being survival at time $j$ (MAR-S) to identify the distribution of dropouts among survivors.

Assumption 4 Dropout among survivors For all $j\geq 1$ and all $t\leq j$

[TABLE]

That is, the outcome is distributed the same among dropouts and non-dropouts conditional on survival and the temporally preceding variables. Similarly, $p(w_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j-1},r_{j}=0,\bar{s}_{j}=1,x_{0})=p(w_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j-1},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0})$ and $p(z_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},r_{j}=0,\bar{s}_{j}=1,x_{0})=p(z_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0})$. Previous studies of the Betula data have shown that individuals who drop out have lower cognitive performance and steeper decline ([14]). In Section 3.2 we introduce sensitivity parameters to allow the dropout to deviate from this MAR type assumption.

We also need three further assumptions for identification of the potential outcomes for those individuals who would survive regardless of exposure history, i.e. the principal strata. We start with two standard assumptions.

Assumption 5 Monotonicity. $S_{j}(\bar{z}_{j})\leq S_{j}(\bar{z}^{\prime}_{j})$; if an individual were to be alive under exposure regime $\bar{z}_{j}$ then (s)he would also be alive under the contrasting regime $\bar{z}^{\prime}_{j}$. Deterministic monotonicity can be too strong in many settings and we discuss a weakening of this in Section 8.

Assumption 6 Differences in outcomes when comparing different strata. For the contrasting exposure regime $\bar{z}^{\prime}_{j}$ we assume, $E[Y_{j}(\bar{z}^{\prime}_{j})\mid\bar{S}_{j}(\bar{z}_{j})=\bar{S}_{j}(\bar{z}^{\prime}_{j})=1]=E[Y(\bar{z}^{\prime}_{j})\mid\bar{S}_{j}(\bar{z}^{\prime}_{j})=1,\bar{S}_{j}(\bar{z}_{j})\neq 1]$. That is, there is no difference in potential outcomes when comparing the ”always survivor” strata to the strata where individuals were to live under the contrasting regime $\bar{z}^{\prime}_{j}$ but not under exposure regime $\bar{z}_{j}$. In Section 3.2 we introduce a sensitivity parameter to investigate sensitivity to this assumption, due to the fact that individuals in the always survivor strata are likely healthier and have better cognitive performance.

A common problem encountered in longitudinal cohort studies is that an individual’s exposure level $z_{j}$, hence the exposure regime $\bar{z}_{j}$, and time-varying confounder $w_{j}$ is only observed if (s)he is alive and participates at the $j$th test wave. Hence we need to introduce a new assumption to be able to identify the probability of survival among exposed and non-exposed; this is necessary for the identification of the potential outcomes among always survivors.

Assumption 7 Exposure and confounding among non-survivors If $s_{j}=0$ and $\bar{s}_{j-1}=1$ for an individual, $z_{j}$ and $w_{j}$ may have occurred before the event of death, thus, $z_{j}$ and $w_{j}$ are not observed but could still be well-defined. We assume,

[TABLE]

and

[TABLE]

i.e. the exposure and confounder are distributed the same among survivors and non-survivors conditional on the temporally preceding variables.

This assumption is used for identification of the principal strata. In the Betula study the cognitive testing is performed at 5 year intervals. Since 5 years is a rather long time period it is likely that some of the participants who died before follow-up were also widowed before death. Thus, the number of widowed participants in the sample may be underestimated and must be accounted for.

3.2 Sensitivity parameters

To investigate sensitivity of Assumption 3 we follow the procedure of [1] (\citeyearbrumback2004sensitivity). The unmeasured confounding is quantified through a parameter which describes the outcome confounding. That is, for exposure regime $\bar{z}_{j}$, $c(\bar{z}_{j})=E[Y_{j}(\bar{z}_{j})\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]-E[Y_{j}(\bar{z}_{j})\mid\bar{y}_{j-1},\bar{z}^{\prime}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]$, where $c(\bar{z}_{j})$ is the average difference in potential outcomes because of unmeasured confounding. The conditional exchangeability assumption does not hold if $c(\bar{z}_{j})\neq 0$. Thus, estimating $E[Y_{j}(\bar{z}_{j})\mid\bar{y}_{j-1},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]$ using the naive estimator $E[Y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]$ leads to a bias of $c(\bar{z}_{j})\times\Pr[z^{\prime}_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]$. Further, since the two regimes only differ in $z_{j}$, for $\bar{z}^{\prime}_{j}$, the bias becomes $c(\bar{z}^{\prime}_{j})\times\Pr[z_{j}\mid\bar{y}_{j-1},\bar{z}^{\prime}_{j-1},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j}=1,x_{0}]$. Sensitivity to several types of unmeasured confounding can be assessed using this form. Here, we restrict to an unmeasured confounder independent of the history of the joint processes $(\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j},x_{0})$.

To investigate sensitivity of Assumption 4 we first make an assumption of non-future dependence (NFD) conditional on survival (NFD-S) for the outcome and then instroduce sensitivity parameters within this partial identifying restrictions ([27]). NFD is a special case of MNAR ([18]), and NFD-S is defined as, $p(y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\{r_{0}=1,\ldots,r_{t-1}=1,r_{t}=0,\ldots,r_{j}=0\},\bar{s}_{j}=1,x_{0})=p(y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0}),$ for all $j>1$ and all $t<j$. Here it is defined conditional on being alive at time $j$. The NFD-S assumption leaves one conditional distribution per incomplete dropout pattern unidentified, that is when $t=j$. To identify the unidentified conditional distribution left by the NFD-S assumption, we introduce a sensitivity parameter $\gamma_{j}$ such that $p(y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j}=\{1,\ldots,1,0\},\bar{s}_{j}=1,x_{0})=p(y_{j}+\gamma_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0})$, when $\gamma_{j}<0$ implies a negative location shift in the outcome at the first unobserved test wave. This assumption implies dropout at time $j$ depends on being alive at that time, the history up to that time, the exposure, time-varying confounder and the outcome at time $j$, but not outcomes or time-varying variables after time $j$. This assumption of dropout not depending on the ’future’ is often viewed as realistic and was proposed originally as a remedy to concerns about many pattern mixture models implicitly having future dependence. Table 1 displays a description of the possible mortality- and missing data patterns under the NFD-S assumption.

To investigate sensitivity of Assumption 6 we let, $\Delta_{\bar{z}^{\prime}_{j}}=E[Y_{j}(\bar{z}^{\prime}_{j})\mid\bar{S}_{j}(\bar{z}_{j})=\bar{S}_{j}(\bar{z}^{\prime}_{j})=1]-E[Y(\bar{z}^{\prime}_{j})\mid\bar{S}_{j}(\bar{z}^{\prime}_{j})=1,\bar{S}_{j}(\bar{z}_{j})\neq 1]$, for the contrasting exposure regime $\bar{z}^{\prime}_{j}$. That is, the mean difference in potential outcomes when comparing the ”always survivor” strata to the strata where individuals were to live under the contrasting regime $\bar{z}^{\prime}_{j}$ but not under exposure regime $\bar{z}_{j}$. In our analysis we assume $\Delta_{\bar{z}^{\prime}_{j}}\geq 0$ which implies that memory performance is on average higher in the ”always survivors”-strata (the always survivors-strata is healthier). We further assume this difference is independent of the preceding variables.

To investigate sensitivity of Assumption 7, we introduce a sensitivity parameter $\nu_{j}$ for the exposure such that,

[TABLE]

representing the mean difference in the proportion exposed between non-survivors and survivors. The first probability on the right-hand side of each expression is not identified. However, bounds can be derived for $\nu_{j}$; see the Web Appendix section A.2 for details. In particular, the upper bound for $\nu_{j}$, $U_{\nu_{j}}$, is obtained when $\Pr[z_{j}\mid\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},r_{j}=0,\bar{r}_{j-1},s_{j}=0,\bar{s}_{j-1}=1,x_{0}]=1$. This reflects that among non-survivors, all subjects were exposed before the event of death between the $j-1$th and $j$th wave. Further, by using Assumption 1 and 5, the lower bound for $\nu_{j}$ is obtained when $\Pr[S_{j}=1\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j-1}=1,x_{0}]=\Pr[S_{j}=1\mid\bar{y}_{j-1},\bar{z}^{\prime}_{j},\bar{w}_{j},\bar{r}_{j},\bar{s}_{j-1}=1,x_{0}]$. This reflects an equal survival probability among those exposed or unexposed at wave $j$. Here, by using the law of total probability and Bayes theorem, the lower bound $L_{\nu_{j}}$ becomes 0.

4 Identification

Identification of the SACE in [2] follows from two results.

Result 1:

The causal contrasts in [1] can be identified as follows

[TABLE]

where $\mathcal{A}$ denotes the set of temporally preceding variables $(\bar{y}_{j-1},\bar{w}_{j},\bar{r}_{j},x_{0})$.

Result 2:

$\tau$ in [2] can further be identified using Assumption 5 by weighting the contrasts in [Result 1:] with

[TABLE]

The proofs of the results can be found in the Web Appendix section A.3. The causal effect is identifiable based on the observed data and Assumptions 1-7, conditional on the fixed values for the sensitivity parameters $c(\bar{z}_{j})$, $c(\bar{z}^{\prime}_{j})$, $\Delta_{\bar{z}^{\prime}_{j}}$, $\nu_{j}$ and, $\gamma_{j}$. For a Bayesian analysis, the sensitivity parameters can be given informative priors. In Section 5.3 and Table 2 we describe the estimation algorithm where the sensitivity parameters are given informative, non-degenerate, priors.

5 Modeling of the observed data distributions and computation of the causal effect

The joint distribution of the observed data is specified as a marginal model for the baseline confounders and a set of sequential conditional models for the time-varying variables, given the history of the joint process (the outcome, exposure, confounders, and missingness). Details of the joint distribution are given in Web Appendix section A.1. The baseline confounders $x_{i0}$ are all observed before an individual enters the study. For each visit $j$ we postulate the time-varying variables in the following order: $(s_{ij},r_{ij},w_{ij},z_{ij},y_{ij})$, even though the exposure, the time-varying confounder, and survival all occurred between $(j-1)$st and $j$th test wave. Of course, $y_{ij}$, $w_{ij}$ and $z_{ij}$, are only observed if $\bar{r}_{ij}=1$ and $\bar{s}_{ij}=1$. It is further allowed that $w_{ij}$ and $z_{ij}$ may have occurred before $s_{ij}$.

5.1 Bayesian semi-parametric modeling

We propose a Bayesian semi-parametric modeling approach based on Bayesian Additive Regression Trees (BART, [3]) for the observed data distribution.

For the time varying components, we specify BART models for the responses as a function of prior histories for all individuals alive and not dropped out at a given test wave. The model consists of two parts: a sum-of-trees model and a regularization prior on the parameters of that model. The model for the continuous response $Y_{j}$ is conditioned on the history of the joint process $(\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},x_{0})$ for the subset that satisfies $\bar{r}_{j}=1$ and $\bar{s}_{j}=1$, and can be expressed as $Y_{j}=\sum_{k=1}^{K_{Y_{j}}}g_{Y_{j}}\left((\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},x_{0});T_{Y_{j}}^{k},M_{Y_{j}}^{k}\right)+\varepsilon_{j}.$ The model consists of $K_{Y_{j}}$ distinct binary regression trees denoted by $T_{Y_{j}}^{k}$. Each tree constitute a set of interior node decision rules leading down to $b_{Y_{j}}^{k}$ terminal nodes, and for a given $T_{Y_{j}}^{k}$, $M_{Y_{j}}^{k}=(\rho_{Y_{j}}^{k,1},\ldots,\rho_{Y_{j}}^{k,b^{k}})$ is the associated terminal node parameters. The conditional distribution of the continuous outcome is specified as normal, $Y_{j}\sim N\left(\mu_{Y_{j}}(\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},x_{0}),\sigma_{j}^{2}\right),$ where the mean function, $\mu_{Y_{j}}(\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},x_{0})$, is given by the sum-of-trees.

The BART models for our binary responses $Z_{j},W_{j},R_{j}$, and $S_{j}$ are specified as probit models. For example the model for the exposure can be expressed as: $\pi_{Z_{j}}(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},x_{0})$

$=\Phi\left(\sum_{k=1}^{K_{Z_{j}}}g_{Z_{j}}\left((\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},x_{0});T_{Z_{j}}^{k},M_{Z_{j}}^{k}\right)\right),$ where $\Phi$ denotes the cumulative density function of the standard normal distribution and $\pi_{Z_{j}}(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j-1},x_{0})$ is the probability of being exposed at wave $j$ given $(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j},x_{0})$ for the subset that satisfies $\bar{r}_{j}=1$ and $\bar{s}_{j}=1$. The BART model for $S_{j}$ is fitted for the subset that satisfies $\bar{r}_{j-1}=1$ and $\bar{s}_{j-1}=1$, and for $R_{j}$ the subset that satisfies $\bar{r}_{j-1}=1$ and $\bar{s}_{j}=1$. The predicted probabilities of $r_{j}=1$ and $s_{j}=1$ are: $\pi_{R_{j}}(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j-1},x_{0})$ and $\pi_{S_{j}}(\bar{y}_{j-1},\bar{z}_{j-1},\bar{w}_{j-1},x_{0}).$ Note that, $s_{0}=1$ and $r_{0}=1$ for all individuals, $\pi_{R_{j}}=0$ if $r_{j-1}=0$, and $\pi_{S_{j}}=0$ if $s_{j-1}=0$.

The baseline confounders are all categorical (age cohort, sex, and education level). We create a saturated multinomial random variable, $x_{0}\sim Multi(N,\pi_{x_{0}}^{1},\pi_{x_{0}}^{2},\ldots,\pi_{x_{0}}^{L}),$ based on these categorical variables. $L$ is the number of categories and each category corresponds to a unique combination of the categorical variables. $\pi_{x_{0}}=(\pi_{x_{0}}^{1},\pi_{x_{0}}^{2},\ldots,\pi_{x_{0}}^{L})$ is given a Dirichlet prior with parameters equal to one.

5.2 Posterior

Draws from the posterior distribution of the sum-of-trees models are generated using Markov chain Monte Carlo (MCMC). The parameters of the conditional distributions for $Y_{j},Z_{j},W_{j},R_{j}$, and $S_{j}$ are assumed independent and thus their posteriors can be sampled simultaneously. BART is implemented in the R package bartMachine ([16]) for continuous and binary responses. We use default priors on all of the parameters of the sum-of-trees model, that is, on the tree structure, the terminal node parameters, and the error variance. For details see [16] (\citeyearkapelner2013bartmachine).

5.3 Computation of the SACE

The algorithm for generating samples from the posterior distribution of $\tau$ in [2] using the G-computation formula is given in Table 2. Details can be found in the Web Appendix section A.4. The algorithm provides the details of generating posterior samples of the causal quantities in Results 1 and 2 (from Section 4) using the posterior distribution of the observed data model parameters (Section 5.1) and the identifying restrictions with sensitivity parameters (Sections 3.1 and 3.2). Recall the expressions in Results 1 and 2 are a function of the observed data distribution and the sensitivity parameters.

For implementation of the algorithm in practice, a number of the initial posterior samples are discarded as burn-in. Parallel computation can be implemented to speed up computations. For example, instead of running one long chain in Step $1$, it is possible to run multiple shorter chains in parallel, although each chain still needs to converge. Also, Step $2$ may be divided into $k$ blocks of size $N^{*}/k$, and in Steps $3-4$ the parameters of interest are computed by combining the pseudo data from the $k$ blocks. We give further details on computation with Betula data in Section 7.3.

6 Simulation study

We performed a simulation study to evaluate the performance of the BSP G-computation algorithm. For simplicity of comparison to other appropriate methods we estimate $\mathrm{E}[Y_{j}(\bar{z}_{j})-Y_{j}(\bar{z}^{\prime}_{j})\mid\bar{s}_{j}=1]$ and set $\Delta_{\bar{z}^{\prime}_{j}}=0$, $\gamma_{j}=0$, and $c(z_{j})=c(z^{\prime}_{j})=0$, i.e. a setting with MAR missingness and and no deaths. Details are found in the Web Appendix section A.5.

We consider two settings for our BSP approach. First, where we specify a normal distribution for the outcome as described in the algorithm (BSP-GC1), and second (BSP-GC2), when specifying a t-distribution with 3 degrees of freedom ($t_{3}$). We compare our approach with three other methods used for causal effect estimation of longitudinal data with time-varying confounding. The three other methods implemented are: (i) A parametric version of the proposed procedure (BP-GC). Here we specified Bayesian linear and logistic additive regression models instead of the BART models described in Section 5.1. (ii) Inverse probability of treatment weights (IPTW; [5]). Here, the mean $E[Y_{j}\mid\bar{s}_{j}=1,\bar{z}_{j}]$ is estimated by averaging the memory outcome for the subset with $\bar{Z}_{j}=\bar{z}_{j}$ in a pseudo-population constructed by weighting each individual using both unstabilized weights (IPTW-W) and stabilized weights (IPTW-SW), to adjust for confounding and for attrition among survivors. The IPTW-W and IPTW-SW were implemented using the ipw and survey packages in R. (iii) Targeted minimum loss-based estimation approach for longitudinal data structures (TMLE; [22]). We implemented the TMLE using the ltmle package using default settings ([25]). Confidence intervals were calculated using nonparametric bootstrap. We used 5000 bootstrap samples. The bootstrap confidence intervals were calculated using the 2.5th and 97.5th percentiles of the resulting estimates.

Data were generated based on a simplified version of the Betula data. We simulated 1000 datasets of size n = 1000. We considered $J_{i}=2$ follow-up test waves, a continuous baseline covariate, $X_{i0}$, generated as $X_{0}\sim Unif\left(0,1\right)$. The outcome, $Y_{ij}$, was considered a continuous time-varying variable. The binary variable $Z_{ij}$ indicated if the subject was widowed or not, and $W_{ij}$ indicated if the spouse been severely sick. Widowhood was an absorbing state, such that, if $Z_{ij}=1$ then $Z_{ik}=1$ for $k\geq j$. Note, that $Z_{i0}=0$ for all subjects. As in the Betula data, all time-varying variables had a highly nonlinear relationship with the baseline covariate and the time-varying confounder interacted with the baseline covariate in the exposure model. Data for the simulation study was generated as $X_{0}\sim Unif\left(0,1\right)$, $W_{j}\sim$ $Ber\left(expit(-2+0.5X_{0}-2X_{0}^{2}+5X_{0}^{3}+0.25W_{j-1})\right)$,

$Z_{j}\sim$ $Ber\left(expit(-5+X_{0}-4X_{0}^{2}+6X_{0}^{3}+0.6W_{j}+0.3W_{j-1}+0.5X_{0}W_{j}-X_{0}^{2}W_{j}+2X_{0}^{3}W_{j})\right)$, and $Y_{j}=0.5-0.05Z_{j}-0.05W_{j}-0.25Y_{j-1}-0.1X_{0}+0.25X_{0}^{2}-0.25X_{0}^{3}+\epsilon_{j}$, where $\epsilon_{j}\sim N\left(0,0.1^{2}\right)$. R code for the data generation is provided in the Web Appendix section A.5.

Table 3 shows the bias, empirical standard deviation (ESD), mean squared error (MSE), and coverage of $95\%$ confidence intervals from the simulation study for BSP-GC1, BSP-GC2, BP-GC, IPTW-W, IPTW-SW, and TMLE. The causal effect estimates for BSP-GC1, BSP-GC2 and TMLE are nearly unbiased. BSP-GC1 and BSP-GC2 are however more efficient (smaller MSE and ESD) and have higher coverage (larger than $95\%$) than TMLE (lower than $95\%$). The simulation results for BSP-GC1 and BSP-GC2 are very similar. As expected, the three other methods; BP-GC, IPTW-W, and IPTW-SW are all biased. Additionally, of all methods, BP-GC has the highest bias and MSE, IPTW-W is the least efficient, and IPTW-SW has the lowest coverage.

To see how the proposed approach performs when the error distribution is misspecified data were instead generated from a $t_{3}$-distribution for the error of the outcome $Y_{j}$ in [LABEL:data_generation_simulation]. Bias, ESD, MSE, and coverage from the simulation are found in Table 4. The results are similar in terms of bias and coverage compared to the previous simulation with correctly specified error. However, ESD and MSE are higher and are now comparable to TMLE.

To see how the proposed approach performs when there is lack of overlap, data were generated as $Ber\left(expit(a_{Z_{j}})\right)I_{X_{0}>0.5}$ for the exposure, mimicking the Betula study where only individuals at older ages were exposed (widowed). The results from the simulation are found in Table 4. The results are similar to the first simulation with correctly specified error and non-linear effects for BSP-GC1 and TMLE. But for the other three approaches bias was higher and CP was lower.

7 Analysis of the Betula data

7.1 The Betula data

The goal is to estimate the causal effect of becoming a widow on memory among those who would survive irrespective of being widowed or not. As such, we limit our data set to those individuals who were married at enrollment. Of approximately 2000 participants $N=1059$ were married at study enrolment, and data were recorded at 4 fixed test waves ($j=0,\ldots,3$) with 5 years interval. The memory outcome was assessed at each wave using a composite of three episodic memory tasks. The score can range between 0 and 76, with a higher score indicating better memory (for details see [14]). We consider two contrasting exposure regimes, subjects who became a widow between the $j-1$th and $j$th wave, $\bar{z}_{j}=\{z_{0}=0,\ldots,z_{j-1}=0,z_{j}=1\}$, and subjects married through test wave $j$, $\bar{z}^{\prime}_{j}=\{z_{0}=0,\ldots,z_{j}=0\}$, for $j=1,2,3$. Baseline demographic characteristics included age-cohorts: $45,50,\ldots,80$ years of age at enrollment, gender, and education, categorized into low: 6-7 years of education (29%), intermediate: 8-9 years (31%), or high: >9 years (40%). We also measured a time-varying confounder; an indicator if the spouse has been sick within the last 5 years. We note that baseline confounders are always recorded.

7.2 Sensitivity parameters

Our approach allows uncertainty about untestable assumptions by specifying priors for the sensitivity parameters described in Section 3.2. We restrict the parameters to a plausible range of values, reflecting the authors’ beliefs about the unknown quantities.

In Section 3.2, the sensitivity parameter $c(\bar{z}_{j})$ reflects the average difference in potential outcomes due to unmeasured confounding (violation of Assumption 3). For the Betula data, when studying the effect of widowhood on cognition, one concern may be that the association is confounded by a healthy lifestyle, such as a healthy diet and/or exercise, something that is often shared within couples. Couples with a healthy lifestyle live longer and may have better cognitive performance than couples with a less healthy lifestyle. This information is not available from the database. Hence, it is a potential unmeasured confounder. Here, we assume $c(z_{j})<0$ and $c(z^{\prime}_{j})>0$, reflecting that exposed (widowed) individuals are less healthy compared to unexposed (married) individuals. We further assume the effect is equal for exposed and unexposed. That is, we assume $c(z_{j})=-\xi_{j}$ and $c(z^{\prime}_{j})=\xi_{j}$. Here, we specify a uniform prior on the sensitivity parameters, $\xi_{j}\sim\mathrm{Unif}(0,U_{\xi_{j}}),$ with upper bound $U_{\xi_{j}}=\frac{1}{2}\times SD(Y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0}).$ That is, we expect the sensitivity parameter not to be bigger than one-half standard deviation of the outcome conditional on the history of the joint process. This approximately corresponds to an effect size similar to that found in previous literature on the effect of Mediterranean diet on memory ([31]).

Departures from a MAR mechanism (Assumption 4) for the missingness among survivors can be investigated by varying $\gamma_{j}$ in Section 3.2. Our prior belief is that $\gamma_{j}<0$, reflecting a negative shift in memory performance occur immediately after the first unobserved test wave. Here, the prior is specified as $\gamma_{j}\sim\mathrm{Unif}(-L_{\gamma_{j}},0),$ where we assume the lower bound is one observed conditional standard deviation, $L_{\gamma_{j}}=1\times SD(Y_{j}\mid\bar{y}_{j-1},\bar{z}_{j},\bar{w}_{j},\bar{r}_{j}=1,\bar{s}_{j}=1,x_{0}).$ The effect is similar to what has been found in previous work examining differences in cognition between completers and those who withdraw, at the last cognitive testing visit before dropping out ([30]).

Sensitivity to Assumption 6, uses $\Delta_{\bar{z}^{\prime}_{j}}$, which reflects the difference in outcomes when comparing the ”always survivor” strata to the strata where individuals were to live under the contrasting regime $\bar{z}^{\prime}_{j}$ but not under exposure regime $\bar{z}_{j}$. We again specify a uniform prior $\Delta_{\bar{z}^{\prime}_{j}}\sim\mathrm{Unif}(0,U_{\Delta_{\bar{z}^{\prime}_{j}}}),$ with upper bound $U_{\Delta_{\bar{z}^{\prime}_{j}}}=1\times SD(Y_{j}\mid\bar{s}_{j}=1).$

Finally, sensitivity to Assumption 7 uses the sensitivity parameter $\nu_{j}$, which represents the difference in the probability of being exposed at wave $j$ for non-survivors and survivors conditioning on the history of the joint process. As shown in Section 3.2, $\nu_{j}$ is restricted to $[0,U_{\nu_{j}}]$. We assume the prior for $\nu_{j}$ is uniform over this range, $\nu_{j}\sim\mathrm{Unif}(0,U_{\nu_{j}}).$ The upper bound reflects that, between the $j-1$th and $j$th wave, all subjects were exposed before death.

7.3 Results and comparison with other methods

We estimated $\tau$ using the proposed BSP method and embedded sensitivity parameters. For each chain the first 1000 iterations were discarded as burn-in, and 2240 posterior samples of $\tau$ were obtained. We sampled pseudo data of size $N^{*}=25000$ at each iteration. Convergence of the posterior samples was monitored using trace plots of the samples. To reduce computation time we used 448 parallel chains. Total computation time was 1 hour and 18 minutes.

For longitudinal exposure regimes limited overlap is not uncommon. To avoid extrapolation of the outcome model outside the range of estimated propensities we restrict the overlap region for the longitudinal exposure regimes. Specifically, we restrict data to the set of individuals that have an estimated propensity score that lies within the range of the observed propensities for the two contrasting regimes (similar to the procedure used in [49]).

We consider two settings for our BSP approach. First, we specify a normal distribution for the residual of the outcome as described in the algorithm (BSP-GC1); second (BSP-GC2), we replace the normal distribution with a t-distribution with 3 degrees of freedom ($t_{3}$). For BSP-GC1, the posterior sampling results revealed a mean episodic memory score of 38.2 (95% CI; $35.4,40.8$) for exposed and 38.1 (95% CI; $35.6,40.1$) for unexposed individuals, and an estimate of $\tau$ of 0.18 (95% CI; -$1.43,1.86$), suggesting that there is no effect of becoming a widow on memory among those who would survive irrespective of exposure. For BSP-GC2, the posterior sampling results revealed a mean episodic memory score of 38.2 (95% CI; $35.5,40.9$) for exposed and 38.0 (95% CI; $35.6,40.1$) for unexposed individuals, and an estimate of $\tau$ of 0.21 (95% CI; -$1.42,1.82$). The conclusions are insensitive to the two choices of outcome residual distribution here.

As a sensitivity analysis we compare how the point estimates and uncertainty varied when setting one sensitivity parameter at a time to zero, while the remaining sensitivity parameters are given the priors described in Section 7.2. Setting $\gamma_{j}$ to zero resulted in a estimate of $\tau$ of $0.18$ (95% CI; $-1.42,1.91$); for $\nu_{j}=0$, $0.20$ (95% CI; $-1.36,1.83$); for $\delta=0$, $0.21$ (95% CI; $-1.38,1.94$); and for $\xi_{j}=0$, $-0.83$ (95% CI; $-2.43,0.75$). The largest differences was found for the analysis setting $\xi_{j}$ to zero (i.e. no unmeasured confounding); however the CI still cover zero and we expect this assumption to not hold. Fixing the other sensitivity parameters at zero had minimal impact.

We also compare our approach, BSP-GC1 with BP-GC, IPTW-W, IPTW-SW, and TMLE (described in the simulation study). For simplicity of comparison we estimate the causal contrasts described in the simulation study. Further to avoid limited overlap, we restrict our data to those age-cohorts where we observe both married and widowed participants over the study period, instead of restricting to the region as for the main analyses. For IPTW-W, IPTW-SW, and TMLE, confidence intervals were calculated using nonparametric bootstrap. We used 5000 bootstrap samples. The bootstrap confidence intervals were calculated using the 2.5th and 97.5th percentiles of the resulting estimates.

The results from all the methods are given in Table 5. First, all of the methods display a negative widowhood effect on memory, although all confidence/credible intervals (CI) cover zero. There is a large discrepancy between our semi-parametric approach, BSP-GC1, and the parametric counterpart, BP-GC. In the latter, the effect was attenuated and the CI was narrower. A likely explanation of the discrepancy in effect estimates is that BP-GC is more susceptible to bias caused by model misspecification. BP-GC and IPTW-SW yielded most similar results, although the weighting approach had much wider CI. Further, the effect estimate appeared most negative using IPTW-W and the CI was much wider than for any of the other methods. Weighting methods are known to be unstable and to have problems with large variance estimates in finite samples if the values of the weights are extreme. In our analysis the range of the weights was 0.06-14.3 for IPTW-W, compared to 0.06-5.4 for IPTW-SW. The large weights using IPTW-W may explain the deviating result using this method. Our BSP-GC1 approach yielded an estimate of $\tau$ most similar to TMLE, although TMLE had slightly wider CI. This is consistent with the results of the simulation study.

8 Concluding remarks

This paper has proposed a Bayesian semi-parametric (BSP) framework for estimating the SACE with longitudinal cohort data. Our approach allows for Bayesian inference under MNAR missingness and truncation by death, as well as the ability to characterize uncertainty about unverifiable assumptions. The proposed approach has several advantages compared to existing approaches: (i) the flexible modeling of the observed data as compared to parametric methods, while maintaining computational ease, (ii) interval estimates for full posterior inference, (iii) easy to introduce sensitivity parameters.

The simulation study, although simplified, mirrored the Betula data. All time-varying variables had a highly nonlinear relationship with the baseline covariate and interaction effects were included. The models for BP-GC, IPTW-W, IPTW-SW, and TMLE were specified using additive effects, and thus, were misspecified. The results showed that BSP-GC1, BSP-GC2 and TMLE were nearly unbiased. BSP-GC1 and BSP-GC2 were however more efficient and had better coverage than TMLE (though a bit conservative). The results are in line with previous research ([35]), suggesting that TMLE is less efficient than Bayesian semi-parametric and non-parametric modeling. This, however, must be explored more thoroughly in future work.

The three other methods; BP-GC, IPTW-W, and IPTW-SW, were all biased. This is expected since these methods make stronger distributional assumptions and thus are more sensitive to model misspecification. Similar to TMLE our approach does not rely on strong modeling assumptions, but unlike TMLE, it is quite easy to modify assumptions and incorporate sensitivity parameters. Recall we could not easily make direct comparisons of the proposed approach with the other approaches under our assumptions that include sensitivity parameters. We attempted to implement Super learner, implemented in the R package SuperLearner, but observed highly variable results for the Betula data (causal effect estimates varied between -0.34 and -1.33). This may be a result of the cross-validation step and the fact that the exposure is a rather rare event. Using our BSP approach these problems are avoided by increasing the size of the pseudo data and running longer chains. Although, computation time can be demanding for large pseudo sample sizes, the algorithm can be fully parallelized as discussed in Section 5.3 and Section 7.3, which would vastly reduced the total computation time.

For the Betula data we did not find an effect of widowhood on memory. The results were not sensitive to different specification of the errors as normal- or t-distributed, and changing the sensitivity parameters one at a time did not change the results significantly either. The difference in findings from previous studies may partly be explained by different estimands being used; ours is the only analysis using a SACE. Additionally, in this study we considered the immediate effect of widowhood (within 5 years) rather than a long term effect; it may take longer for degeneration to become apparent.

Our approach can be generalized in various ways. For example, it is possible to allow for multiple time-varying confounders and/or continuous baseline confounders using a sequential approach as proposed by ([47]). This would involve first ordering the confounders into sequential conditionals and then applying BART to model each of these univariate conditionals. Additionally, although widowhood status is thought of as a monotone exposure pattern and an absorbing state in this study, this is not essential for the proposed approach and other (non-monotone) exposure regimes, such as the effect of widowhood duration on memory at the last visit, might be of interest and are possible to study with a few modifications (for example, the positivity assumption).

Violations of the consistency assumption can be problematic when using observational data ([4]). For example, the effect of widowhood can affect memory via different pathways, e.g. for some subjects via stress or depression and others via reduced physical health due to poorer lifestyle choices ([10]). This is a limitation with the current study where the exposure is defined homogeneously, and should be explored more thoroughly in future work.

Several of our assumptions can be (further) relaxed. For example, Assumption 5 can be weakened to a stochastic Monotonicity, by following the procedure described in [24] (\citeyearlee2010causal). Also, in this study we have considered unmeasured outcome confounding; this assumption can easily be extended to allow unmeasured mortality confounding. Assumption 6 can be weakened by conditioning on the history of the joint process. However, a drawback with relaxing these assumptions is increasing the number of sensitivity parameters.

One limitation with BART is the restrictive, and sometimes unrealistic, assumption of IID normal errors, (although they can easily be replaced with heavier tail errors as in here). A fully non-parametric modeling approach could be obtained by extending BART to model the error distribution using the Dirichlet process mixtures ([9]). An additional limitation of the proposed approach is that we used existing R-functions for BART that are not most efficient for our setting. We will explore these limitations in future work, as well as, other choices for priors of the sensitivity parameters.

Supplementary materials

Web Appendices referenced in Sections 3, 4, 5, 6, and 7, as well as R code are available as Supplementary materials.

Acknowledgments

The authors would like to thank Dr Anna Sundström for helpful discussions on the interpretation of the results. This work is partially funded by The Swedish Foundation for Humanities and Social Sciences P17-0196:1 and Paths to Healthy and Active Ageing, funded by the Swedish Research Council for Health, Working Life and Welfare, (Dnr 2013 – 2056) to MJ. This work is partially funded by US NIH grants CA183854 and GM112327 to MJD. This publication is based on data collected in the Betula prospective cohort study, Umeå University, Sweden. The Betula Project is supported by Knut and Alice Wallenberg foundation (KAW) and the Swedish Research Council (K2010-61X-21446-01). The simulations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Umea University partially funded by the Swedish Research Council through grant agreement no. 2016-07213.

9 Tables and Figure