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

Now for the pooled estimates. The pooled estimate of the regression coefficients is easy.

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

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

(betweenVar <- mean(mat[,2])) # mean of variances
(withinVar <- sd(mat[,1])^2) # variance of variances
(dfCorrection <- (nrow(mat)+1)/(nrow(mat))) # dfCorrection

(totVar <- betweenVar + withinVar*dfCorrection) # total variance
(pooledSE <- sqrt(totVar)) # standard error

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

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

Put them all together

(pooledEstimates <- round(c(pooledMean, pooledSE, pooledP),5))

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

fit <- with(data=nhimp,exp=lm(chl~bmi))

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$ – Peter Ellis Feb 7 '18 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 '18 at 2:09
  • 1
    $\begingroup$ This should help: github.com/stefvanbuuren/mice/blob/master/R/mice.df.r $\endgroup$ – Peter Ellis Feb 7 '18 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 '18 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 '20 at 14:16

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

  lambda <- (withinVar + (withinVar/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))

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. So we can pass it into degrees of freedom argument in the pt() function in order to obtain the p-value for the imputed model.

pt(q = pooledMean / pooledSE, df = nu_BR, lower.tail = FALSE) * 2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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