Embedded Files
Chan, K. W. & Meng, X.-L. (2022). Multiple Improvements of Multiple Imputation Likelihood Ratio Tests. Statistica Sinica, 32, 1489–1514.
Abstract
Multiple imputation (MI) inference handles missing data by first properly imputing the missing values m times, and then combining the m analysis results from applying a complete-data procedure to each of the completed datasets. However, the existing method for combining likelihood ratio tests has multiple defects: (i) the combined test statistic can be negative in practice when the reference null distribution is a standard F distribution; (ii) it is not invariant to re-parametrization; (iii) it fails to ensure monotonic power due to its use of an inconsistent estimator of the fraction of missing information (FMI) under the alternative hypothesis; and (iv) it requires non-trivial access to the likelihood ratio test statistic as a function of estimated parameters instead of datasets. This paper shows, via both theoretical derivations and empirical investigations, that essentially all of these problems can be straightforwardly addressed if we are willing to perform an additional likelihood ratio test by stacking the m completed datasets as one big completed dataset. A particularly intriguing finding is that the FMI itself can be estimated consistently by a likelihood ratio statistic for testing whether the m completed datasets produced by MI can be regarded effectively as samples coming from a common model. Practical guidelines are provided based on an extensive comparison of existing MI tests.
Main paper and Supplementary note (last update: 31 September 2021)
R-code (The key functions are copied below)
Figure: Idea of the proposed test.
Figure: Power curves of different multiple imputation (MI) tests. Null hypothesis: the means of a bivariate Normal random vector are equal. Level of significance is 0.5%. The tests were performed under non-standard parametrizations; see Figure 5 of the main paper for a detailed discription.
W-1,2,3,4: different MI Wald tests; see Rubin (2004, Wiley);
L-1: the best existing MI likelihood ratio test (LRT) by Meng and Rubin (1992, Biometrika);
L-2: a trivial modification of L-1 (to ensure positivity of the test statistics and the estimator of fraction of missing information (FMI));
L-3: first proposed MI LRT, which guarantees invariance and positivity of test statistics;
L-4: a trivial modification of L-3 (to ensure positivity of the estimator of FMI);
L-5: second proposed MI LRT, which guaranteed (i) positivity of both test statistic and estimate of FMI, (ii) invariance to parametrization, and (iii) consistency of estimator of FMI; and
L-0: benchmark, i.e., the LRT without using MI.
milrt.stack = function(data, lrt, k){ # data = MI dataset (list) # lrt = complete-data LRT statistic, which converges weakly to chi^2_k (function) # k = dim(H1)-dim(H0) (numeric) m = length(data) n = dim(data[[1]])[1] d = dim(data[[1]])[2] datastack = matrix(NA,nr=n*m,nc=d) for(l in 1:m){ datastack[(l-1)*n+(1:n),]=data[[l]] } dbar = mean(unlist(mapply(lrt, data=data, SIMPLIFY=FALSE))) dhat = lrt(datastack)/m rhat = max(0,(m+1)/(k*(m-1))*(dbar-dhat)) fhat = rhat/(1+rhat) Dhat = dhat/(k*(1+rhat)) df1 = k df2 = k*(m-1)/fhat^2 pval = c("p-value"=1-pf(Dhat,df1,df2)) pval}
milrt.robust = function(data, lrt, tml, k, h){ # data = MI dataset (list) # lrt = complete-data LRT statistic, which converges weakly to chi^2_k (function) # tml = twice maximized complete-data log-likelihood (function) # k = dim(H1)-dim(H0) (numeric) # h = dim(H1) (numeric) m = length(data) n = dim(data[[1]])[1] d = dim(data[[1]])[2] datastack = matrix(NA,nr=n*m,nc=d) for(l in 1:m){ datastack[(l-1)*n+(1:n),]=data[[l]] } deltabar = mean(unlist(mapply(tml, data=data, SIMPLIFY=FALSE))) deltahat = tml(datastack)/m dhat = lrt(datastack)/m rhatR = (m+1)/(h*(m-1))*(deltabar-deltahat) fhatR = rhatR/(1+rhatR) DhatR = dhat/(k*(1+rhatR)) df1 = k df2R = h*(m-1)/fhatR^2 pvalR = c("p-value"=1-pf(DhatR,df1,df2R)) pvalR}
milrt.robust = function(data, lrt, tml, k, h){ # data = MI dataset (list) # lrt = complete-data LRT statistic, which converges weakly to chi^2_k (function) # tml = twice maximized complete-data log-likelihood (function) # k = dim(H1)-dim(H0) (numeric) # h = dim(H1) (numeric) m = length(data) n = dim(data[[1]])[1] d = dim(data[[1]])[2] datastack = matrix(NA,nr=n*m,nc=d) for(l in 1:m){ datastack[(l-1)*n+(1:n),]=data[[l]] } deltabar = mean(unlist(mapply(tml, data=data, SIMPLIFY=FALSE))) deltahat = tml(datastack)/m dhat = lrt(datastack)/m rhatR = (m+1)/(h*(m-1))*(deltabar-deltahat) fhatR = rhatR/(1+rhatR) DhatR = dhat/(k*(1+rhatR)) df1 = k df2R = h*(m-1)/fhatR^2 pvalR = c("p-value"=1-pf(DhatR,df1,df2R)) pvalR}
Google Sites
Report abuse