# Calculating Barnard-Rubin degrees of freedom from Sums of Squares in an F-test

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

require(mice)
require(car)
set.seed(123)
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) mat 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?

• I swapped p and F for $p$ and $F$ because $F$ renders in titles while F appears as *F*. If you don't agree with this change, you can roll back the edit with my apologies. More information can be found in math.meta.stackexchange.com/questions/5020/… – Reinstate Monica Mar 28 at 20:59