Using Amelia in R, I obtained multiple imputed datasets. After that, I performed a repeated measures test in SPSS. Now, I want to pool test results. I know that I can use Rubin's rules (implemented through any multiple imputation package in R) to pool means and standard errors, but how do I pool p-values? Is it possible? Is there a function in R to do so? Thanks in advance.

  • $\begingroup$ You might want to check out information about p-value meta-analysis. One good starting point: en.wikipedia.org/wiki/Fisher%27s_method $\endgroup$ – user29889 Sep 4 '13 at 7:00
  • $\begingroup$ You might refer to miceadds package pool_mi function. It provide more and complete functionality including p-value and confidence interval estimation. The readme page of that page provides lots of example. $\endgroup$ – J.D. Jul 8 '20 at 19:42

Yes, it is possible and, yes, there are R functions that do it. Instead of computing the p-values of the repeated analyses by hand, you can use the package Zelig, which is also referred to in the vignette of the Amelia-package (for a more informative method see my update below). I'll use an example from the Amelia-vignette to demonstrate this:

amelia.out <- amelia(freetrade, m = 15, ts = "year", cs = "country")

zelig.fit <- zelig(tariff ~ pop + gdp.pc + year + polity, data = amelia.out$imputations, model = "ls", cite = FALSE)

This is the corresponding output including $p$-values:

  Model: ls
  Number of multiply imputed data sets: 15 

Combined results:

lm(formula = formula, weights = weights, model = F, data = data)

                Value Std. Error t-stat  p-value
(Intercept)  3.18e+03   7.22e+02   4.41 6.20e-05
pop          3.13e-08   5.59e-09   5.59 4.21e-08
gdp.pc      -2.11e-03   5.53e-04  -3.81 1.64e-04
year        -1.58e+00   3.63e-01  -4.37 7.11e-05
polity       5.52e-01   3.16e-01   1.75 8.41e-02

For combined results from datasets i to j, use summary(x, subset = i:j).
For separate results, use print(summary(x), subset = i:j).

zelig can fit a host of models other than least squares.

To get confidence intervals and degrees of freedom for your estimates you can use mitools:

imp.data <- imputationList(amelia.out$imputations)
mitools.fit <- MIcombine(with(imp.data, lm(tariff ~ polity + pop + gdp.pc + year)))
mitools.res <- summary(mitools.fit)
mitools.res <- cbind(mitools.res, df = mitools.fit$df)

This will give you confidence intervals and proportion of the total variance that is attributable to the missing data:

              results       se    (lower    upper) missInfo    df
(Intercept)  3.18e+03 7.22e+02  1.73e+03  4.63e+03     57 %  45.9
pop          3.13e-08 5.59e-09  2.03e-08  4.23e-08     19 % 392.1
gdp.pc      -2.11e-03 5.53e-04 -3.20e-03 -1.02e-03     21 % 329.4
year        -1.58e+00 3.63e-01 -2.31e+00 -8.54e-01     57 %  45.9
polity       5.52e-01 3.16e-01 -7.58e-02  1.18e+00     41 %  90.8

Of course you can just combine the interesting results into one object:

combined.results <- merge(mitools.res, zelig.res$coefficients[, c("t-stat", "p-value")], by = "row.names", all.x = TRUE)


After some playing around, I have found a more flexible way to get all necessary information using the mice-package. For this to work, you'll need to modify the package's as.mids()-function. Use Gerko's version posted in my follow-up question:

as.mids2 <- function(data2, .imp=1, .id=2){
  ini <- mice(data2[data2[, .imp] == 0, -c(.imp, .id)], m = max(as.numeric(data2[, .imp])), maxit=0)
  names  <- names(ini$imp)
  if (!is.null(.id)){
    rownames(ini$data) <- data2[data2[, .imp] == 0, .id]
  for (i in 1:length(names)){
    for(m in 1:(max(as.numeric(data2[, .imp])))){
        indic <- data2[, .imp] == m & is.na(data2[data2[, .imp]==0, names[i]])
        ini$imp[[names[i]]][m] <- data2[indic, names[i]]

With this defined, you can go on to analyze the imputed data sets:

imp.data <- do.call("rbind", amelia.out$imputations)
imp.data <- rbind(freetrade, imp.data)
imp.data$.imp <- as.numeric(rep(c(0:15), each = nrow(freetrade)))
mice.data <- as.mids2(imp.data, .imp = ncol(imp.data), .id = NULL)

mice.fit <- with(mice.data, lm(tariff ~ polity + pop + gdp.pc + year))
mice.res <- summary(pool(mice.fit, method = "rubin1987"))

This will give you all results you get using Zelig and mitools and more:

                  est       se     t    df Pr(>|t|)     lo 95     hi 95 nmis   fmi lambda
(Intercept)  3.18e+03 7.22e+02  4.41  45.9 6.20e-05  1.73e+03  4.63e+03   NA 0.571  0.552
pop          3.13e-08 5.59e-09  5.59 392.1 4.21e-08  2.03e-08  4.23e-08    0 0.193  0.189
gdp.pc      -2.11e-03 5.53e-04 -3.81 329.4 1.64e-04 -3.20e-03 -1.02e-03    0 0.211  0.206
year        -1.58e+00 3.63e-01 -4.37  45.9 7.11e-05 -2.31e+00 -8.54e-01    0 0.570  0.552
polity       5.52e-01 3.16e-01  1.75  90.8 8.41e-02 -7.58e-02  1.18e+00    2 0.406  0.393

Note, using pool() you can also calculate $p$-values with $df$ adjusted for small samples by omitting the method-parameter. What is even better, you can now also calculate $R^2$ and compare nested models:


mice.fit2 <- with(mice.data, lm(tariff ~ polity + pop + gdp.pc))
pool.compare(mice.fit, mice.fit2, method = "Wald")$pvalue
  • 1
    $\begingroup$ Great answer, just wanted to point out a slight typo, I think you meant: mice.res <- summary(pool(mice.fit, method = "rubin1987")). $\endgroup$ – FrankD Oct 14 '15 at 10:39
  • $\begingroup$ Good catch. I have corrected the typo. $\endgroup$ – crsh Oct 14 '15 at 12:41

Normally you would take the p-value by applying Rubin's rules on conventional statistical parameters like regression weights. Thus, there is often no need to pool p-values directly. Also, the likelihood ratio statistic can be pooled to compare models. Pooling procedures for other statistics can be found in my book Flexible Imputation of Missing Data, chapter 6.

In cases where there is no known distribution or method, there is an unpublished procedure by Licht and Rubin for one-sided tests. I used this procedure to pool p-values from the wilcoxon() procedure, but it is general and straightforward to adapt to other uses.

Use procedure below ONLY if all else fails, as for now, we know little about its statistical properties.

lichtrubin <- function(fit){
    ## pools the p-values of a one-sided test according to the Licht-Rubin method
    ## this method pools p-values in the z-score scale, and then transforms back 
    ## the result to the 0-1 scale
    ## Licht C, Rubin DB (2011) unpublished
    if (!is.mira(fit)) stop("Argument 'fit' is not an object of class 'mira'.")
    fitlist <- fit$analyses
        if (!inherits(fitlist[[1]], "htest")) stop("Object fit$analyses[[1]] is not an object of class 'htest'.")
    m <- length(fitlist)
    p <- rep(NA, length = m)
    for (i in 1:m) p[i] <- fitlist[[i]]$p.value
    z <- qnorm(p)  # transform to z-scale
    num <- mean(z)
    den <- sqrt(1 + var(z))
    pnorm( num / den) # average and transform back
  • 1
    $\begingroup$ @ Stef van Buuren what do you mean by 'take the p-value by applying Rubin's rules on conventional statistical parameters like regression weights'? How does the pool() function in your package (which is excellent by the way) arrive at the pooled p-value? $\endgroup$ – llewmills Feb 7 '18 at 1:17
  • $\begingroup$ I used this for computing pooled p-values on simulated Fisher exact tests; is there any more recent literature on this? $\endgroup$ – gregmacfarlane Dec 21 '20 at 20:35
  • $\begingroup$ Found this dissertation which seems to describe the coded methods in Chapter 4. d-nb.info/101104966X/34 $\endgroup$ – gregmacfarlane Dec 21 '20 at 20:40
  • $\begingroup$ Apart from the work from Licht, there's also a suggestion by Eekhout et al (2017) to take the median P-value if you're doing logistic regression <stefvanbuuren.name/fimd/sec-multiparameter.html#sec:chi> $\endgroup$ – Stef van Buuren Dec 22 '20 at 21:22

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