3
$\begingroup$

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')

enter image description here


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)
$\endgroup$

1 Answer 1

1
$\begingroup$

To start with

> 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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – jay.sf
    Commented Jun 7 at 2:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.