probability of treatment regimen $(Z_1^{(b)}=z_1^{(b)},Z_2^{(b)}=z_2^{(b)})$ is denoted as $P^*_{z_2^{(b)}}=P_{z_2^{(b)}}P^*_{z_1^{(b)}}$.

8.  Estimate the model $$Y_{z_1^{(b)},z_2^{(b)}}^{(b)}|P^*_{z_2^{(b)}},Z_1^{(b)}=z_1^{(b)},Z_2^{(b)}=z_2^{(b)},Y_{z_1^{(b)}}^{(b)},X_{z_1^{(b)}}^{(b)},X_{z_1^{(b)},z_2^{(b)}}^{(b)},W_1^{(b)},W_2^{(b)},V^{(b)}$$ again restricting to subjects that are observed and predicted to survive under the treatment regimes of interest at $t=2$. Use the respective models to impute the counterfactual of $Y_{z_1^{(b)},z_2^{(b)}}^{(b)}$.

9.  Using a similar procedure for steps 5-8 with the restriction determined by $S_{z_1^{(b)},\ldots,z_{t-1}^{(b)}}^{(b)}=S_{z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}=1$ for time $t$ where at least one $z_t^{(b)}\neq z_t^{'(b)}$, extend the estimation until the desired time point $t=T$.

10. Repeat Steps 1-9 to obtain $B$ bootstrap values for $$\hat{\Delta}_{z_1^{(b)},\ldots,z_{t-1}^{(b)},z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}=E[Y_{z_1^{(b)},\ldots,z_{t-1}^{(b)}}^{(b)}-Y_{z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}|S_{z_1^{(b)},\ldots,z_{t-1}^{(b)}}^{(b)}=S_{z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}=1].$$ with associated pooled variance $W_{z_1^{(b)},\ldots,z_{t-1}^{(b)},z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}$.

11. The estimate of $$\Delta_{Z_1,\ldots,Z_t,Z'_1,\ldots,Z'_t}=E[Y_{Z_1,\ldots,Z_t}-Y_{Z'_1,\ldots,Z'_t}|S_{Z_1,\ldots,Z_{t-1}}=S_{Z'_1,\ldots,Z'_{t-1}}=1]$$ is then $$\bar{\Delta}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}=\sum_{b=1}^B(\hat{\Delta}_{z_1^{(b)},\ldots,z_{t-1}^{(b)},z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)})/B,$$ and the estimate of the variance of $\bar{\Delta}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}$ is $$T_B=\bar{W}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}+(1+1/B)D_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B},$$ where $$\bar{W}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}=\sum_{b=1}^B(W_{z_1^{(b)},\ldots,z_{t-1}^{(b)},z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)})/B$$ and $$D_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}=\sum_{b=1}^B\frac{(\hat{\Delta}_{z_1^{(b)},\ldots,z_{t-1}^{(b)},z_1^{'(b)},\ldots,z_{t-1}^{'(b)}}^{(b)}-\bar{\Delta}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B})^2}{B-1}.$$ The estimate $\Delta_{Z_1,\ldots,Z_t,Z'_1,\ldots,Z'_t}$ follows a $t$ distribution with degree of freedom $\nu$, $$\frac{\Delta_{Z_1,\ldots,Z_t,Z'_1,\ldots,Z'_t}-\bar{\Delta}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}}{\sqrt{T_B}}\sim t_{\nu},$$ where $\nu=(B-1)(1+\frac{\bar{W}_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}}{D_{z_1,\ldots,z_t,z'_1,\ldots,z'_t,B}(B+1)})^2$.

*Remark*. The idea of including the BART estimated propensity score within BART as a predictor in Steps 4 and 8 is not new. [@hahn] showed that including a BART estimated propensity score as a predictor within BART improved the estimation of heterogenous treatment effects for observational studies. [@tan] also reported that the inclusion of the BART estimated propensity score as a predictor within BART to impute missing data, under the missing 