Likelihood ratio test for multiply imputed datasets?

Question

I have a set of generalized linear models fit to 5 multiply imputed datasets. I am interested in testing the statistical significance of a set of predictors, coded from dummy variables.

Rubin's Rules makes conducting 1 degree of freedom Wald tests very easy, since the general form of the test statistic $\hat{\beta}/SE_{RR}(\hat{\beta})$ is available from basic arithmetic provided the model summaries for each multiply imputed dataset. I understand this test statistic also has good asymptotic properties, relative to the general Wald test. But I do not understand how to go about 2 degree-of-freedom (or higher) tests. Putting aside the issue of finding higher dimensional Rubin's Rules, is it tractable to just use the likelihoods for each model, and conduct a test by doing:

$$-2 \sum_{s=1}^5 \sum_{i=1}^n \left( \log \mathcal{L}(X_{s,i}, y_{s,i}, \hat{\beta}_s) - \log \mathcal{L} (X_{s,i}, y_{s,i}, \hat{\beta}_{s,0}) \right) \sim_d \chi^2_k$$

Where $\mathcal{H}_0: \beta = [\beta_1, \beta_2, \ldots, \beta_{p-k-1}, 0, 0, \ldots, 0]$ is a $p$-length vector having the last $k \le p$ entries equal to 0

Answer 1

Apparently, the mice package has an option pool.compare. The help file cites the following publication by Meng, Rubin "Performing Likelihood Ratio Tests with Multiply-Imputed Data Sets" in Biometrika 1992. I cannot access this with my limited research library. However, looking at the source code for pool.compare suggests that the implementation is a very clever trick indeed.

The outcome, of any scale, is compared to a binomial likelihood by creating $Z = (Y - \min(Y))/\text{range}(Y)$ and taking a logistic inverse transformation of the linear predictors $E[Z|X] = \mbox{logit}^{-1} (\mathbf{X}\hat{\beta})$. Two log likelihood ratios are defined for each dataset, the $D_L$ which uses the pooled coefficient estimate, $\bar{Q}$ and the $D_M$ which uses the coefficient estimates from each dataset. The $D_L$ and $D_M$ are averaged together, and some algebra is done to create the test-statistic. Namely, $r_m$ is defined as

$$r_m = \frac{m+1}{k (m-1)} (D_M - D_L)$$

$$T = D_L/(k (1+r_m))$$

which has a $\chi^2_{df}$ distribution with $df = 4(mk-4)(1+(1-2/(mk))^2)$ if $mk>4$.