Calculating pooled p-values manually

Question

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.

require(mice)
set.seed(123)
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))
summary(pool(fit))
#          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.

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.