Bootstrap, Rubin's rules, and uncertainty of sub-estimates? Can someone provide an intuition for why, when using bootstrap to calculate the variability of an estimate (say a regression coefficient $\beta$) we don't need to incorporate the uncertainty of each bootstrapped estimate ($\hat\beta_i$)?
On the other hand, why is it that when using Rubin's rules we do need to incorporate the uncertainty of the parameter estimate on each imputed dataset?
What is the intuition for this difference between these procedures?
(Edit:) In other words, why is the uncertainty of the MI sub-estimates an underestimate of the uncertainty of the  MI estimate? Also, what would happen if we tried to incorporate the uncertainty of the bootstrap sub-estimates (e.g. in linear regression, by sampling from each $\beta_i$ t-distribution) into the bootstrap procedure? Would it change anything?
 A: (This is a fairly long answer, there is a summary at the end)
The theoretical justifications given below, based mainly on this, this, this, this, and this articles, can help you get the intuition you are after. Otherwise, let me know.
I. The theoretical justification on MI is large sample bayesian approximation, as follows:


*

*Using iterative procedures, we create draws from the posterior distribution of $θ$ (or $\beta$ to keep with your notation).

*In that case, a large number of draws $D$ are needed.

*If we assume normality of the observed-data posterior distribution, we need to estimate only the mean and variance–much less $D$ are needed.

*In that case, a very limited $D$ are required to estimate reliably the distribution mean.

*MI is based on this idea.

*If we assume $p(θ | y_{obs})$ is approximately normal, the observed-data posterior can be effectively determined by the posterior mean and variance, $ E(θ | y_{obs})$ and $Var(θ | Y_{obs}) $

*Note that $E(θ | y_{obs}) = E[E(θ | y_{mis}, y_{obs}) | Y_{obs}] =
\int E(θ | y_{mis}, y_{obs})p(y_{mis} | y_{obs})dy_{mis}$
where the outer expectation is taken with respect to the
posterior predictive distribution, $p(y_{mis} | y_{obs})$


*And that,
$Var(θ | Y_{obs}) = E[Var(θ | Y_{mis}, Y_{obs}) | Y_{obs}]+
Var[E(θ | Y_{mis}, Y_{obs}) | Y_{obs}]$
where the outer expectation is taken with respect to the
posterior predictive distribution, $p(y_{mis} | y_{obs})$


*In that case, $E[Var(θ | Y_{mis}, Y_{obs}) | Y_{obs}] =
\int Var(θ | Y_{mis}, Y_{obs})p(Y_{mis} | Y_{obs})dY_{mis}$ and, 


$Var[E(θ | Y_{mis}, Y_{obs}) | Y_{obs}] = \int E^2(θ | Y_{mis}, Y_{obs})p(Y_{mis} | Y_{obs})dY_{mis}−(\int  E(θ | Y_{mis}, Y_{obs})p(Y_{mis} | Y_{obs}))^2$


*For large $D$,


$E[E(θ| Y_{mis}, Y_{obs})|Y_{obs}] \approx \frac{1}{D} \sum_{d=1}^D \widehat θ_d$
$E[Var(θ|Y_{mis}^{(d)},y_{obs}) \approx 
\bar V = \frac{1}{D} \sum_{d=1}^D Var(θ|Y_{mis}^{(d)},y_{obs})$
$Var[E(θ|y_{mis}, y_{obs})|y_{obs}] \approx 
B = \frac{1}{D - 1}\sum_{d=1}^D (\widehat θ - \bar θ)^2$
where,
$Y_{mis}^{(d)}$ are independent draws of $Y_{mis} $ from the posterior
predictive distribution $p(y_{mis} | y_{obs})$, $\widehat θ_d = E(θ|Y_{mis}^{(d)},y_{obs})$ the complete-data posterior mean of $θ$ calculated for the $d$th imputed data set $(Y_{mis}^{(d)}, Y_{obs})$, 
$Var(θ|Y_{mis}^{(d)},y_{obs})$ the complete-data posterior
variance of $θ$ calculated for the calculated for the $d$th imputed data set $(Y_{mis}^{(d)}, Y_{obs})$ 
and,
$\bar θ = \frac{1}{D}\sum_{d=1}^D (\widehat θ)$
$\bar V$: within-imputation variance
$B$: between-imputation variance
Therefore:

Total variance $T$ can be estimated using $T = \bar V + (1 + D^{-1})B$
and approximated using $T \approx \bar V + B$
while MI point estimate for $E(θ|y_{obs}$), (that is, for $θ$) is obtained using $\bar θ = \frac{1}{D}\sum_{d=1}^D (\widehat θ)$

II. A bootstrap algorithm for estimating the variance of $θ_n$ (or $\beta$ to keep with your notation) and for constructing confidence intervals works as follows:


*

*We would like to find the variance of $θ_n$. 
Let, $Var_P(θ_n) = Var_P (g(X1, . . . , X_n)) ≡ S_n(P)$

*Note that $Var_P (θ_n)$ is some function of P (and n) so I have written $Var_P(θ_n) = S_n(P)$
If we knew $P$, we could approximate $ S_n(P)$ by simulation as follows:
(a) draw $X1, . . . , X_n ∼ P$
(b) compute $\widehat θ_n^{(1)} = g(X1, . . . , X_n)$
(c) draw $X1, . . . , X_n ∼ P$
(d) compute $\widehat θ_n^{(2)} = g(X1, . . . , X_n)$
(e) draw $X_1, . . . , X_n ∼ P$
(f) compute $\widehat θ_n^{(B)} = g(X1, . . . , Xn)$
Let $s^2$ be the sample variance of $\widehat θ_n^{(1)} , . . . , \widehat θ_n^{(B)} $
So $s^2 = \frac{1}{B} \sum_{i=1}^B (\widehat θ)^2 - (\frac{1}{B} \sum_{i=1}^B \widehat θ)^2$
By the law of large numbers, $s^2 \rightarrow E[(\widehat θ)^2] - (E[(\widehat θ)])^2 = Var(\widehat θ) = S_n(P)$
Since we can take $B$ as large as we want, we have that $s^2 \approx Var_P(\widehat θ_n)$
In other words, we can approximate $S_n(P)$ by repeatedly simulating $n$ observations from $P$.
But we don’t know $P$. So we estimate $S_n(P)$ with $S_n(P_n)$ where $P_n$ is the empirical distribution. Since $P_n$ is a consistent estimator, we expect that $S_n(P_n) ≈ S_n(P)$
In other words:
Bootstrap approximation of the variance: estimate $S_n(P)$ with $S_n(P_n)$ 
Or in other words $\widehat {Var_P(θ_n)} = Var_{Pn}(\widehat θ_n)$
Question: But how do we compute $S_n(P_n)$
Answer: We use the simulation method above, except that we simulate
from $P_n$ instead of $P$. This leads to the following algorithm:

Bootstrap Variance Estimator:

(a) Draw a bootstrap sample $X_1^*, ..., X_n^*$
(b) Compute $ \widehat θ_n^* = g(X_1^*, ..., X_n^*)$
(c) Repeat the previous step B times, yielding estimator $ \widehat θ_{n,1}^*, ..., \widehat θ_{n,B}^*$
(d) Compute:
$\widehat s = \sqrt{\frac{1}{B} \sum_{j=1}^B (\widehat θ_{n,1}^* - \bar θ)^2}$
where $\barθ = \frac{1}{B} \sum_{j=1}^B \widehat θ_{n,j}^*$
(e) Output $\widehat s$
Therefore, you can think about it like this:
$\frac{1}{B} \sum_{j=1}^B (\widehat θ_{n,1}^* - \bar θ)^2 \approx S_n(P_n) \approx
S_n(P) $
These last two approximations are our two sources of error in this
estimation of variance. The first is due to the fact that n is finite
and the second is due to the fact that B is finite. However, we can
make B as large as we like. (In practice, it usually suffices to take
B = 10, 000.) So we ignore the error due to finite B.

Summary:
Question 1: Why is the uncertainty of the MI sub-estimates an
  underestimate of the uncertainty of the MI estimate? 
Answer: Because of the law of total variance as illustrated in step
  I.8 above, $Var(θ | Y_{obs}) = E[Var(θ | Y_{mis}, Y_{obs}) |
 Y_{obs}]+ Var[E(θ | Y_{mis}, Y_{obs}) | Y_{obs}]$. Otherwise, you will
  be violating this law. See article 3 for more details about the
  justification of the law of total variance. However, if you do $D$ multiple imputations to get $D$ datasets, and then you generate $B$ bootstrap samples from each of the $D$ imputed datasets yielding $D * B$ datasets $D_{d,b}^*$, and if in each of these datasets the quantity of interest is estimated, that is, $ \widehat θ_{d,b}^*$. Then, the pooled sample of the ordered estimates $θ_{MIBP}^* = \{\widehat θ_{d,b}^*; b = 1, ..., B; d = 1, ..., D\}$  can be used to construct a correct confidence interval for $θ$ by using $[\widehat θ_{lower}, \widehat θ_{upper}]_{MIBP} = [\widehat θ_{MIBP}^{*, \alpha} \widehat θ_{MIBP}^{*,1 - \alpha}]$, instead of the law of total variance, where $\widehat θ_{MIBP}^{*, \alpha}$ is the $\alpha$-percentile of the ordered bootstrap estimates $θ_{MIBP}^*$. For more details see article 6.
Question 2: What would happen if we tried to incorporate the
  uncertainty of the bootstrap sub-estimates (e.g. in linear regression,
  by sampling from each  t-distribution) into the bootstrap
  procedure? Would it change anything?
Answer: Because of the law of large numbers as illustrated in step II.2.f) above, $s^2 \rightarrow
E[(\widehat θ)^2] - (E[(\widehat θ)])^2 = Var(\widehat θ) = S_n(P)$ is
  a valid approximation of the true variance. Therefore, doing what you
  are suggesting would simply add unnecessary noise to the variance. See
  article 4 for more details about the justification of the law of
  large numbers.

