0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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/… $\endgroup$ – Reinstate Monica Mar 28 at 20:59

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.