For reasons I won't go into I need to calculate parameter estimates from several imputed datasets. Based on this CV post about Rubin's rules I have determined how to manually calculate both the pooled coefficient and standard error. However the method for the p-value eludes me. If it were up to me I would make do with the pooled coefficient and se, but I know my boss will want the p-value.

Here is the toy analysis.

First the imputation.

nhimp <- mice(nhanes)
v <- sapply(1:5, function(i) {
  fit <- lm(chl ~ bmi, data=complete(nhimp, i))
  print(c('coef'=coef(fit)[2], 'var'=vcov(fit)[2, 2], 'p'=summary(fit)$coefficients["bmi",4]))

Create a matrix of extracted estimates from the model applied to each of the five complete imputed datasets. 1st column are the coefficients, 2nd column are the variances, 3rd column are the p-values.

(mat <- t(v))

#      coef.bmi      var          p
# [1,] 2.195180 5.467482 0.35758491
# [2,] 4.231113 3.603410 0.03587456
# [3,] 3.470647 5.475586 0.15159999
# [4,] 2.937763 4.776171 0.19198054
# [5,] 2.208305 2.972943 0.21304488

Now for the pooled estimates. The pooled estimate of the regression coefficients is easy: it's just the mean of the coefficients (first column

(pooledMean <- mean(mat[,1]))

# [1] 3.008602

Calculating the pooled estimate of the standard error is a bit more tricky, but still relatively simple. Total variance is the sum of within-variance and between-variance*degrees-of-freedom-correction.

Within variance is the average of the imputation specific point estimate variances

(withinVar <- mean(mat[,2])) # mean of variances
# [1] 4.459118

Between variance is the variance of the coefficients (variance of first column, or sd of first column squared.

(betweenVar <- sd(mat[,1])^2) # variance of coefficients
# [1] 0.7537916

The degrees of freedom correction is (m+1)/m where m is the number of imputations

(dfCorrection <- (nrow(mat)+1)/(nrow(mat))) # dfCorrection
# [1] 1.2

Now we can calculate total variance

(totVar <- withinVar + betweenVar*dfCorrection) 
# [1] 5.363668

The pooled standard error is just the square root of the total variance

(pooledSE <- sqrt(totVar)) # standard error
# [1] 2.315959

Now is the part I don't know: how to get the pooled estimate of the p-value

(pooledP <- mean(mat[,3])) #??????
# [1] 0.190017

Put them all together

(pooledEstimates <- round(c(pooledMean, pooledSE, pooledP),5))
# [1] 3.00860 2.31596 0.19002

These should be exactly the same as the pooled values for these parameters returned by mice

fit <- with(data=nhimp,exp=lm(chl~bmi))
#          term   estimate std.error statistic       df   p.value
# 1 (Intercept) 111.958092 61.373512  1.824209 15.93028 0.0869345
# 2         bmi   3.008602  2.315959  1.299074 15.68225 0.2126945

The manually calculated pooled coefficient and se are the same as those yielded by the pool() function; but not the p-value. Can anyone explain simply the way mice calculates the pooled p-value? This post explains how to do it with software but I need to calculate it manually.

  • 1
    $\begingroup$ The trick is the degrees of freedom, right? Because it's just the upper tail probability of getting more than the absolute value of a t statistic pt(q = pooledMean / pooledSE, df = 18.21792, lower.tail = FALSE) * 2. So you just need to know where 18.21792 comes from? The help file for ?pool says it uses the "Barnard-Rubin adjusted degrees of freedom". $\endgroup$ Feb 7, 2018 at 1:46
  • $\begingroup$ Thank you @Peter Ellis. That was surprisingly simple. Every answer begs another question (as it should), so now I need to work out how to calculate the Barnard-Rubin adjusted degrees of freedom, and that should get me to my final destination. $\endgroup$
    – llewmills
    Feb 7, 2018 at 2:09
  • 1
    $\begingroup$ This should help: github.com/stefvanbuuren/mice/blob/master/R/mice.df.r $\endgroup$ Feb 7, 2018 at 2:13
  • $\begingroup$ Thanks @Peter Ellis. Just had my Eureka moment, and have posted it below. Even though I quit smoking years ago I feel like a cigarette. $\endgroup$
    – llewmills
    Feb 7, 2018 at 3:34
  • 1
    $\begingroup$ Might betweenVar and withinVar been switched in the opening post? Also asking because of the solution you mentioned. $\endgroup$
    – user297944
    Oct 2, 2020 at 14:16

1 Answer 1


This is for anyone who is interested, after reading pp. 37-43 in Flexible Imputation of Missing Data by Stef van Buuren. If we call the adjusted degrees of freedom nu

  m <- nrow(mat)
  lambda <- (betweenVar + (betweenVar/m))/totVar
  n <- nrow(nhimp$data)
  k <- length(coef(lm(chl~bmi,data = complete(nhimp,1))))
  nu_old <- (m-1)/lambda^2  
  nu_com <- n-k
  nu_obs <- (nu_com+1)/(nu_com+3)*nu_com*(1-lambda)
  (nu_BR <- (nu_old*nu_obs)/(nu_old+nu_obs))
  # [1] 15.68225

nu_BR, the Barnard_Rubin adjusted degrees of freedom, matches up with the degrees of freedom for the bmi variable yielded from the the summary(pool(fit)) call above: 15.68225. So we can pass this value into degrees of freedom argument in the pt() function in order to obtain the two-tailed p-value for the imputed model.

pt(q = pooledMean / pooledSE, df = nu_BR, lower.tail = FALSE) * 2
# [1] 0.2126945

And this manually calculated p-value now matches the p-value from the mice function output.


Your Answer

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

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