# Bootstrap standard errors after multiple imputation

Following Rubin's rules for multiple imputation, I've calculated pooled estimates, group means in this case, with pooled standard errors.

I checked this with a bootstrap and, assuming pooled standard errors are correct, the bootstrap standard errors are much smaller.

Since bootstrap standard errors are rather conservative, I suspect that we cannot simply use the standard deviations of the bootstrap estimates as standard errors. How can we bootstrap the standard errors correctly?

I'm happy to provide a minimal example in R below.

> ## imputation
> m. <- 20
> imp <- mice::mice(iris2, m=m., print=FALSE, seed=42)  ## impute
Warning message:
Number of logged events: 1
> long <- mice::complete(imp, action='long')  ## combine
> r0 <- do.call('rbind', by(long, ~ .imp, \(x) tapply(x$$Sepal.Length, x$$Species, mean)))
> colMeans(r0)  ## estimates
setosa versicolor
4.981452   5.971290

> ## compute estimates and standard errors analytically
> MV <- by(long, ~ .imp, \(x)
+          sapply(c(mean=mean, var=var), \(f)
+                 tapply(x$$Sepal.Length, x$$Species, f))) |>
+   simplify2array()
> Q <- rowMeans(MV[, 1, ])  ## calculate mean estimates
> U <- rowMeans(MV[, 2, ])  ## calculate within variances
> B <- rowSums(((MV[, 1, ] - Q)^2))/(m. - 1)  ## calculate between variances
> T <- U + (1 + 1/m.)*B  ## calculate total variances
> cbind(Estimate=Q, 'Std. Error'=sqrt(T))
Estimate Std. Error
setosa     4.981452  0.3923879
versicolor 5.971290  0.5775910


> ## bootstrap standard errors
> R <- 999
> b <- pbmcapply::pbmclapply(seq_len(R), \(i) {
+   tryCatch({
+     m. <- 20
+     imp <- mice::mice(iris2[sample.int(nrow(iris2), replace=TRUE), ],
+                       m=m., print=FALSE)
+     long <- mice::complete(imp, action='long')
+     r0 <- do.call('rbind', by(long, ~ .imp, \(x)
+                               tapply(x$$Sepal.Length, x$$Species,
mean)))
+     colMeans(r0)
+   }, error=\(e) array(dim=c(1, 2)))
+ }, mc.cores=7L)
> B <- b |> do.call(what='rbind')
> matrixStats::colSds(B, na.rm=TRUE) |>  t() |>
rownames<-('Bootstrap Std. Error')
setosa versicolor
Bootstrap Std. Error 0.06895931  0.1079427


> ## bootstrap distributions
> par(mfrow=c(1, 2))
> hist(B[, 'setosa'], 'fd'); hist(B[, 'versicolor'], 'fd')


Data:

> iris2 <- subset(iris, Species %in% c('setosa', 'versicolor')) |>
as.matrix()
> set.seed(42)
> iris2[sample.int(length(iris2), length(iris2)*.4)] <- NA
> iris2 <- type.convert(as.data.frame(iris2), as.is=TRUE)


> with(iris, by(Sepal.Length,Species,sd))
Species: setosa
[1] 0.3524897
------------------------------------------------------------
Species: versicolor
[1] 0.5161711
------------------------------------------------------------
Species: virginica
[1] 0.6358796
> sem<-function(x) sd(x)/sqrt(length(x))
> with(iris, by(Sepal.Length,Species,sem))
Species: setosa
[1] 0.04984957
------------------------------------------------------------
Species: versicolor
[1] 0.07299762
------------------------------------------------------------
Species: virginica
[1] 0.08992695


So, the standard error estimates you are getting analytically are close to the standard deviations of Sepal.Length by Species in the complete data. The bootstrap standard errors are a bit larger than the complete-data standard errors of the means. This suggests the analytic estimates are wrong

Using mice for the analytic standard errors gives

> mice::pool(with(imp, lm(Sepal.Length~1, subset=Species=="versicolor")))
Class: mipo    m = 20
term  m estimate       ubar           b          t dfcom       df
1 (Intercept) 20  5.97129 0.01068082 0.002386659 0.01318681    30 21.87703
riv    lambda      fmi
1 0.2346254 0.1900377 0.255155
> mice::pool(with(imp, lm(Sepal.Length~1, subset=Species=="setosa")))
Class: mipo    m = 20
term  m estimate        ubar           b           t dfcom       df
1 (Intercept) 20 4.981452 0.004847555 0.003518128 0.008541589    30 13.81826
riv    lambda       fmi
1 0.7620408 0.4324763 0.4999653
> sqrt(0.008541589)
[1] 0.09242072
> sqrt(0.01318681)
[1] 0.1148338


which is a lot closer to the bootstrap results.

Ok, so where is the problem with the analytic results? I note that B is much smaller than U, and that B is of the same order of magnitude as the bootstrap variances. And I think you're calculating U as the mean of the estimated variance of the data, not of the means.

So: the bootstrap seems about right.

• Thanks very much, Thomas. If I use c(mean=mean, var=\(x) var(x)/length(x)), i.e. use the variance of the sample mean for each imputation, I get the same results as mice for analytical standard errors. Commented Jun 7 at 2:41