I am wondering if it is possible to manually calculate a pooled $p$-value from an omnibus $F$-test on several multiply-imputed datasets. It is possible with a regression, where there are unique estimates and standard errors for every comparison (see here) but I am not sure how to do it when the test is an $F$-test of an omnibus effect and the only available variance information is sums of squares and degrees of freedom.

Here is a working example

nhanes$age <- factor(nhanes$age) # make age a factor
nhimp <- mice(nhanes) # multiply impute
anPool <- with(nhimp, aov(chl ~ bmi*age)) # anova on each of the mutiply imputed datasets
summary(pool(anPool)) # pool the results

Note that the pool() function automatically converts the aov() object back to an lm() format to get the pooled estimates.

I am going to try and do this manually with the Anova() function from car, since it returns type-3 sums of squares. I will try to derive a pooled $p$-value for the (factorised) age variable to illustrate, but the procedure could be applied to any of the omnibus effects.

First we run the Anova() on each of the 5 imputed datasets and extract the omnibus $F$-value for the age variable, and calculate mean squared error for the age variable (between-observation variance) and mean squared error for the residual term (within-observation variance)

v <- sapply(1:5, function(i) {
  an <- Anova(lm(chl ~ bmi*age, data=complete(nhimp, i)), type = 3)
  print(c(F_age = an$`F value`[3],
          mse_age = an$`Sum Sq`[3]/an$Df[3],
          mse_within = an$`Sum Sq`[5]/an$`Df`[5]))

Transpose this so we get a matrix of $F$, mse-between, and mse-within for each of the five imputed datasets

mat <- t(v)

        F_age  mse_age mse_within
[1,] 2.633596 2078.889   789.3726
[2,] 4.920537 4038.557   820.7554
[3,] 2.304795 2298.430   997.2384
[4,] 4.648544 2343.350   504.1040
[5,] 6.802849 4162.495   611.8753

Now we need to apply Rubin's rules (see here) to get the pooled standard deviation.

First we need the number of imputations

m <- nrow(mat)

Now we calculate the pooled estimates of the $F$ statistic, which is simply the average of the $F$s across all the imputed datasets

(f_pooled_age <- mean(mat[,"F_age"]))
[1] 4.262064

Now for the pooled estimate of the variance. This is the part I am not sure about because I am applying Rubin's rules to the mean square error for the omnibus $F$-test, whereas in regression you use the standard error for each coefficient. Anyhow total variance is sum of between-term and the within-term variance. We also need to calculate a correction term

# between
(betweenVar_age <- mean(mat[,"mse_age"])) # mean of F-values
[1] 2984.344

# within
(withinVar_age <- sd(mat[,"mse_age"])^2) # variance of F-values
[1] 1050152

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

Now the total variance is between + (within x dfCorrection)

(totVar_age <- betweenVar_age + withinVar_age*dfCorrection)
[1] 1263166

And the standard deviation we get from the square root of the variance

(pooledSD_age <- sqrt(totVar_age))
[1] 1123.907

Now we use these values to obtain the Barnard-Rubin adjusted degrees of freedom. First we calculate the lambda

(lambda_age <- (withinVar_age + (withinVar_age/m))/totVar_age)
[1] 0.9976374

(Which already looks wrong)

And then apply this to calculating the degrees of freedom

n <- nrow(nhimp$data)
k <- 3
nu_old_age <- (m-1)/lambda_age^2 
nu_com <- n-k
nu_obs_age <- (nu_com+1)/(nu_com+3)*nu_com*(1-lambda_age)
(nu_BR_age <- (nu_old_age*nu_obs_age)/(nu_old_age+nu_obs_age))
[1] 0.04725655

Now we pass these degrees of freedom into the upper tail probability of getting more than the absolute value of a $t$ statistic and use the square root of the pooled $F$ statistic (because t = sqrt(F) to get the pooled $p$-value

(pooledP_age <- pt(q = abs(sqrt(f_pooled_age) / pooledSD_age), df = nu_BR_age, lower.tail = FALSE) * 2)
[1] 0.9872084

Now I know this is totally wrong because when I run the Anova() on one of the imputed datasets...

Anova(lm(chl ~ bmi*age, data = complete(nhimp, 1)), type = 3)

             Sum Sq Df F value   Pr(>F)   
(Intercept)   200.3  1  0.2537 0.620269   
bmi         11690.3  1 14.8096 0.001083 **
age          4157.8  2  2.6336 0.097830 . 
bmi:age      3531.2  2  2.2367 0.134180   
Residuals   14998.1 19  

...they are way off. Does anyone have any idea if calculating pooled $p$-values manually from multiply-imputed omnibus $F$-tests is possible, and, if it is, how to do it?


Your Answer

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

Browse other questions tagged or ask your own question.