Does MCMC error mean you need to increase your iterations for reporting?

Question

I am working with multiple imputation and read almost universally how the number of imputations needed is relatively small. So it is recommended to use $M=5$ or thereabouts, that is the default is Stef Van Buuren's mice package in R. The rationale for this is not explained, but if we are to accept mice as being efficient and unbiased, it seems reasonable that even with MCMC error, the resulting estimates could be considered reasonable.

However, I find that when I set different seeds, create imputed datasets, then calculate the results, the results will differ out at the second or third decimal place, and I am finding that primary inference may even differ.

For reporting (such as in academic manuscripts), should estimates be free of MCMC error, or is it okay as long as the results are replicable and the Seed is prespecified rather than chosen to provide "ideal" results.

Ian White has a paper in Stats in Medicine 2010 377-399 which is a good summary of more recent thinking (well, up to 2010 anyway) — mdewey, Jul 18 '17 at 12:59

score 4 · Accepted Answer · answered Jul 18 '17 at 09:05

The answer to this question can boil down to one's personal opinion, so here is my personal opinion: as a statistician, I must, whenever possible, report standard errors. MCMC standard errors are no exception. If only the results are replicable and the seed is provided, the authors are assuming that the standard errors serve the purpose of reproducibility, which is only a secondary benefit. The primary purpose of standard errors in MCMC is to ensure enough samples were obtained to estimate the quantity of interest. That is, there were enough samples to ensure that the standard error was small.

The paper Markov Chain Monte Carlo: Can We Trust the Third Significant Figure? almost entirely addresses this question. They discuss how very few academic journals report MCMC standard errors, and that this must change. They then go on to explain how these errors can be calculated and reported. Since in MCMC, standard errors can theoretically be made arbitrarily small, they ask you to simulate samples until the standard error is less than some pre-specified number. In this case you don't have to directly report standard errors, since it is evident that the standard errors are small. The theory for the above paper comes from Jones et.al (2006).

Much work has been done since the above paper(s). New work has indicated that for multivariate problems, since standard errors are matrices, it doesn't make sense in reporting them. Instead, one can simulate until the multivariate effective sample size is smaller than a pre-specified lower bound (check my answer here).

To conclude, it is important to report standard errors in one form or the other. Since in MCMC, standard errors can be made arbitrarily small, one can replace reporting standard errors with reporting effective sample size as well. (see Vats et. al).

score 1 · Answer 2 · answered Jul 26 '17 at 20:09

A challenge with multiple imputation is that it usually makes use of Rubin's Rules to pool regression models fit on simulated datasets, or rather datasets with randomly generated imputations. The "between" model error summarizes both the uncertainty of the imputation model as well as the MCMC error. Since it is computationally inexpensive to fit potentially 100 or 1,000 fixed effects regression models, it's possible to consider a cumulative series of imputations and use a logarithmic plot to project the extent of MCMC error.

In the mice package in R, the b slot from a returned mipo (multiple imputation pooled outcome) object supposedly summarizes the between-imputation error. Interestingly, such a logarithmic plot tends to show an increase in such error over time, but converging to values at 30 imputations or more, which befuddles me slightly. If I took anything away from this, it suggests that more imputations may be safer bet. In this example, increasing the imputations led to a pooled fit which had less precise inference relative to the complete case analysis.

set.seed(123)
library(mice)
imp <- mice(nhanes, m=50, print=F)
simData <- lapply(1:50, complete, x=imp)
fits <- lapply(simData, lm, formula=chl ~ hyp + bmi + factor(age))

pool.x <- function(fits, ...) { # trick MICE to pooling custom fits, just needs one slot
  fits <- list(analyses=fits)
  class(fits) <- 'mira'
  pool(fits, ...)
}

vars <- sapply(2:50, function(i) sqrt(diag(pool.x(fits[1:i], method='RR')$b)))
matplot(t(vars), type='l', xlim=c(-20, 50))
legend('left', lty=1:5, col=1:5, rownames(vars), )

Does MCMC error mean you need to increase your iterations for reporting?

2 Answers2