Submitted to the Annals of Statistics
arXiv: arXiv:0000.0000

MULTIPLE IMPROVEMENTS OF MULTIPLE
IMPUTATION LIKELIHOOD RATIO TESTS

arXiv:1711.08822v1 [] 23 Nov 2017

By Kin Wai Chan and Xiao-Li Meng
Harvard University
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.

1. Historical Successes and Failures. Missing-data problems are
ubiquitous in practice, to the extent that the absence of any missingness
often is a strong indication that the data have been pre-processed or manipulated in some way (Blocker and Meng, 2013). Multiple imputation (MI)
(Rubin, 1978, 2004) has been a preferred method by many practitioners,
especially those who are ill-equipped to handle missingness on their own,
due to lack of information or skills or resources. MI relies on the data collector (e.g., a census bureau) to build a reliable imputation model to fill
in the missing data m times, so the users of the data can apply their favorite software or procedures that are designed to handle complete data,
AMS 2000 subject classifications: Primary 62D05; Secondary 62F03, 62E20.
Keywords and phrases: Fraction of missing information; missing data; invariant test;
monotonic power; robust estimation.
Meng thanks NSF and JTF for partial financial support. He is also embarrassed by
the research immaturity displayed in his thesis, and very thankful to Keith’s (Kin Wai)
creativity and diligence which lead to the beautiful remedies presented here.

1

2

K.W. CHAN AND X.-L. MENG

and do so m times. MI inference, e.g., hypothesis testing, is then performed
by appropriately combining these m complete-data results.
Although MI was designed initially for public-use datasets, over the years
it has become a method of choice for handling missing data in general, because it separates the handling missingness from conducting analysis (e.g.,
Tu et al., 1993; Rubin, 2004). This convenient separation, however, comes
with the thorny issue of uncongeniality, which refers to the incompatibility between the imputation model and the subsequent analysis procedures
(Meng, 1994a). This issue is examined in greater detail recently by Xie and
Meng (2017), which shows that the problem of uncongeniality is easiest
to address when the imputation model is more saturated than the user’s
model/procedure, and when the user is conducting most efficient analysis,
such as likelihood inference. The current paper, therefore, tackles the basic question of how to conduct MI likelihood ratio tests (LRTs), assuming
the imputation model is sufficiently saturated to render the validity of the
common assumptions, which we will review, made in the literature about
conducting LRTs with MI.
Like many hypothesis testing procedures in common practice, the exact
null distributions of various MI test statistics, LRTs or not, are intractable.
This intractability is not computational, but rather statistical due to the
well-known issue of nuisance parameter, that is, the lack of pivotal quantity,
as highlighted, historically, by the Behrens-Fisher problem (Wallace, 1980).
Indeed, the nuisance parameter in the MI context is the so-called “the fraction of missing information” (FMI), which is determined by the ratio of the
between-imputation variance to within-imputation variance (and its multivariate counterparts), and hence the challenge we face is almost identical to
the one faced by the Behrens-Fisher problem (Meng, 1994b).
An added challenge in the MI context is that the user’s complete-data
procedures can be very restrictive, especially when a user is using a routine/software as a black box. What is available to the user could vary from
the entire likelihood function, to point estimators such as MLE and Fisher
information, to a single p-value. As it should be intuitive, the more which
is available to the user, the more reliable and effective MI tests we can construct. Therefore, there have been a variety of procedures proposed in the
literature, depending on what quantities we assume the user has access to,
as we shall review shortly.
Among them, a statistically and computationally more promising idea is
to directly combine individual LRTs. However, the execution of this idea as
presented in Meng and Rubin (1992) relied too heavily on the usual asymptotic equivalence (in terms of the data size, not the number of imputations,

MULTIPLE IMPUTATION LR TESTS

3

m) between the LRT and Wald test under the null. Whereas this provides
an asymptotically valid test in the sense that it has the correct type-I error
asymptotically, its performance in practice can deteriorate rather quickly
when the data size becomes small, capable of delivering negative “F test
statistic” or FMI. Worse of all, the test can have essentially zero power
because of its badly inconsistent estimate of FMI under the alternative hypothesis. In addition, the combining rule of Meng and Rubin (1992) requires
the user to have access to the LRT as a function of parameter values, not
just as a function of the data, a requirement which typically is not met when
the user relies on a software routine.
The main purpose of this paper is to derive MI LRTs that are free of the
defects as listed in the abstract and detailed in § 1.3 below. We achieve this
by switching the order of two main operators in the combining rule of Meng
and Rubin (1992): Maximizing the average of the m log-likelihoods instead of
averaging the maximums of them. This automatically retrieves the standard
properties of LRTs, namely, positivity, invariance and monotonic power, and
nearly eliminates the problem of negative estimates of the FMI. To fully
resolve the remaining problem, we propose a new estimator of FMI that is
guaranteed to be non-negative and consistent regardless of the validity of
the null hypothesis. A computationally feasible algorithm is also derived. In
particular, one of the proposed test statistics requires users to have access
to the complete data LRT as a function of data only.
The remainder of the current section provides necessary background and
notation. § 2 then discusses the defects of the existing MI LRT and our remedies. § 3 investigates computational requirements for our proposals, including theoretical considerations and comparisons. § 4 is devoted to empirical
investigations including simulation experiments and a real-world example.
§ 5 concludes the article with a call for future work. Supplementary results
can be found in Appendix A. All technical proofs are relegated to Appendix
B. Extra simulation results are placed in Appendix C.
1.1. Notation and Complete-data Tests. Let Xobs and Xmis be, respectively, the observed and missing parts of an intended complete dataset
X “ Xcom “ tXobs , Xmis u, which consists of n observations. Denote the
sampling model — probability or density, depending on the data type —
of X by f p¨ | ψq, where ψ P Ψ Ă Rh is a vector of parameters. The goal
is to test H0 : θ “ θ0 against H1 : θ ­“ θ0 , when only Xobs is available,
where θ “ θpψq P Θ Ă Rk is a function of ψ, and θ0 is a specified vector. (In
this paper, we focus on the standard two-sided alternative, but the results
here generalize to more specific forms of H1 .) Denote the true values of ψ

4

K.W. CHAN AND X.-L. MENG

and θ by ψ ‹ and θ‹ . Here we assume either Xobs is rich enough that the
missing data mechanism is ignorable (Rubin, 1976), or more generally the
imputation model has incorporated the missing data mechanism (since the
imputer may have additional confidential data that are not a part of Xobs ).
p
Let θp “ θpXq
and U “ Uθ “ Uθ pXq be respectively the complete-data
p (e.g., the inverse of the observed
MLE of θ and an efficient estimator of Varpθq
p
p
p
Fisher information). Also, let ψ0 “ ψ0 pXq and ψp “ ψpXq
be respectively
the H0 -constrained and unconstrained complete data MLEs of ψ, and Uψ “
p For testing H0 against H1 , the
Uψ pXq be an efficient estimator of Varpψq.
p U q{k and the LRT
common choices include the Wald statistic DW “ dW pθ,
p ψp0 | Xq{k, where
statistic DL “ dL pψ,
p U q “ pθp ´ θ0 q| U ´1 pθp ´ θ0 q,
dW pθ,

p ψp0 | Xq “ 2 log
dL pψ,

p
f pX | ψq
.
f pX | ψp0 q

Under regularity conditions (RCs), such as those in § 4.2.2 and § 4.4.2 of Serfling (2001) when the rows of X are independent and identically distributed,
we have the following classical results.
Property 1.1. Under H0 , (i) DW ñ χ2k {k and DL ñ χ2k {k; and (ii)
pr
pr
DW ´ DL Ñ 0 as n Ñ 8, where “ñ” and “Ñ” denote convergence in
distribution and in probability, respectively.
When only Xobs is available, testing θ “ θ0 is more involved. For MI,
` u, ` “ 1, . . . , m, be the m completed datasets, where
let X ` “ tXobs , Xmis
` are drawn conditionally independently from a proper imputation model
Xmis
given Xobs ; see Rubin (2004). We can then carry out a complete-data estimation or testing procedure on X ` , ` “ 1, . . . , m, resulting in a set of m
quantities. The so-called MI inference is to appropriately combine them to
obtain a single answer, such as a p-value.
1.2. MI Wald’s Test and Fraction of Missing Information. Let d`W “
p ` q and U ` “ U pX ` q be the imputed counterparts of
dW pθp` , U ` q, θp` “ θpX
p U q, θp and U , respectively, for each `. Also, write their averages as
dpθ,
(1.1)

dW “

m
1 ÿ `
d ,
m `“1 W

θ“

m
1 ÿ p`
θ,
m `“1

U“

m
1 ÿ `
U .
m `“1

Under congeniality (Meng, 1994a), one can show that asymptotically (Rubin
and Schenker, 1986) Varpθq can be consistently estimated by
(1.2)

T “ U ` p1 ` 1{mqB,

where

B“

m
1 ÿ p`
pθ ´ θqpθp` ´ θq|
m ´ 1 `“1

5

MULTIPLE IMPUTATION LR TESTS

is known as the between-imputation variance, in contrast to Ū in (1.1), which
measures within-imputation variance. Intriguingly, 2T serves as a universal
(estimated) upper bound of Varpθq under uncongeniality (Xie and Meng,
2017). Under RCs, we have that, as m, n Ñ 8,
pr

pr

npU ´ Uθ q Ñ 0,

npT ´ Tθ q Ñ 0,

pr

npB ´ Bθ q Ñ 0

for some deterministic matrices Uθ , Tθ and Bθ “ Tθ ´ Uθ , where 0 denotes
a matrix or vector of zeros having appropriate dimensions, and the subscript
θ is used to highlight that these matrices are for the estimation of θ, because
there are also corresponding Tψ , Bψ , Uψ for the entire parameter ψ. Similar
to U , T and B, we also define U ψ , Tψ and Bψ for the component ψ. Note that
if θpcom and θpobs are the MLEs of θ based on Xcom and Xobs (upon assuming
congeniality), respectively, then Uθ l Varpθpcom q and Tθ l Varpθpobs q as n Ñ
8, where An l Bn means that An ´ Bn “ op pAn ` Bn q.
The straightforward MI Wald test DW pT q “ dW pθ, T q{k is not practical
because T is singular when m ă k (usually 3 ď m ď 10). Even when it is
not singular, it is usually not a very stable estimator of Tθ because m is
small. To resolve it, Rubin (1978) made the following assumption of equal
fraction of missing information (EFMI).
Assumption 1 (EFMI of θ).

There is r ě 0 such that Tθ “ p1 ` rqUθ .

This assumption implies that we believe the missing data have caused
an equal amount of increase in variance for all components of θ, clearly a
strong assumption. In this case, the quantity f “ r{p1 ` rq is the FMI,
and r is the odds of the missing information, though r itself sometimes
is referred to as FMI in the literature. Under Assumption 1, Rubin (2004)
1 qU , where r
1
proposed to replace T by the full-rank matrix p1 ` rrW
rW
“
´1
p1`1{mq trpBU q{k is an estimator of rm “ p1`1{mqr (a finite imputation
1 can be computed via
corrected version of r), and shown that rrW
1
rrW
“

(1.3)
(1.4)

1

dW “

m`1
1
pd ´ dr1W q,
kpm ´ 1q W

m
1 ÿ p`
dpθ , U q
m `“1

and

dr1W “ dpθ, U q,

where the prime “1” emphasizes U is used instead of individual tU ` um
`“1 .
Then, Rubin (2004) proposed a very simple MI Wald test statistic:
(1.5)

r1 “
D
W

dr1W
1 q.
kp1 ` rrW

6

K.W. CHAN AND X.-L. MENG

It is worthwhile to emphasize the intuition behind (1.3)–(1.5) because the
forms here are the fundamental building blocks for virtually all the subsequent developments. The “obvious” Wald test form dr1W {k is too large as a
test statistic because it uses the inverse of the complete-data variance U ,
which is larger than the inverse of the correct (observed-data) variance T .
1 q factor attempts to correct this, and the amount of correction
The p1 ` rrW
is determined by the (average) amount of between-imputation variance relative to the within-imputation variance. Expression (1.3) shows that this
relative amount can be estimated by contrasting the average of individual
Wald test statistics and the Wald test statistic based on an average of individual estimates. This theme of using the difference between “average of
functions” and “function of average”, namely,
(1.6)

AvetGpxqu ´ GpAvetxuq

is a common practice, including the usual variance by taking Gpxq “ x2 ; see
Meng (2002) for other examples.
r 1 is intractable, Li et al. (1991b)
Since the exact null distribution of D
W
proposed to approximate it by Fk,dfpr
r r1 ,kq , the F distribution with degrees
W
r r1 , kq, where Km “ kpm ´ 1q and
of freedom k and dfpr
W

"
(1.7)

r m , kq “
dfpr

4 ` pKm ´ 4qt1 ` p1 ´ 2{Km q{rm u2 , if Km ą 4;
pm ´ 1qp1 ` 1{rm q2 pk ` 1q{2,
otherwise.

This approximation assumes n is sufficiently large so that the standard
asymptotic χ2 distribution in Property 1.1 can be used. If n is small, the
small sample degree of freedom in Barnard and Rubin (1999) can be used.
1.3. The Current MI Likelihood Ratio Test and Its Defect. Let d`L “
p ` q be the imputed counterparts
dL pψp0` , ψp` | X ` q, ψp0` “ ψp0 pX ` q and ψp` “ ψpX
p respectively, for each `. Let their averages be
of dL pψp0 , ψp | Xq, ψp0 and ψ,
(1.8)

dL “

m
1 ÿ `
d ,
m `“1 L

ψ0 “

m
1 ÿ p`
ψ ,
m `“1 0

ψ“

m
1 ÿ p`
ψ.
m `“1

1 , Meng and Rubin (1992) proposed estimating r by
Similar to rrW
m

(1.9) rrL “

m`1
pdL ´ drL q,
kpm ´ 1q

where

m
1 ÿ
drL “
dL pψ 0 , ψ | X ` q,
m `“1

and hence it is again in the form of (1.6). Computation of rrL requires users
to have access to (i) a subroutine for pX, ψ0 , ψq ÞÑ dL pψ0 , ψ | Xq, and (ii) the

MULTIPLE IMPUTATION LR TESTS

7

estimates ψp0` and ψp` , rather than the matrices U and B. Therefore computing
1 . The resulting MI LRT is
rrL is easier than computing rrW
(1.10)

rL “
D

drL
,
kp1 ` rrL q

whose null distribution can be approximated by Fk,dfpr
r rL ,kq .
The main theoretical justification (and motivation) was the asymptotic
equivalence between the complete-data Wald test statistic and LRT statistic, under the null, as stated in Property 1.1. This equivalence permitted
1
the replacement of dW and dr1W in (1.3) respectively by dL and drL in (1.9).
However, this is also where the problems lie.
First, with finite samples, 0 ď drL ď dL is not guaranteed, consequently
r L ě 0 or rrL ě 0. Since D
r L is referred to an F distribution and rrL
nor are D
r L or rrL will cause trouble.
estimates rm ě 0, clearly negative values of D
r
Second, DL is not invariant to re-parameterization of ψ because it relies
on the average of the m estimates of ψ. For each individual LRT statistic
d`L and bijective map g such that ϕ “ gpψq, we have
(1.11)

p`0 q, g ´1 pϕ
p` q | X ` q,
d`L “ dL pψp0` , ψp` | X ` q “ dL pg ´1 pϕ

where ϕ
p`0 and ϕ
p` are the constrained and unconstrained MLEs of ϕ based on
X ` . However, the MI LRT statistic drL no longer has this property because
m
m
1 ÿ
1 ÿ
`
r
dL “
dL pψ 0 , ψ | X q ­“
dL pg ´1 pϕ0 q, g ´1 pϕq | X ` q
m `“1
m `“1

ř
ř
p` .
p`0 and ϕ “ m´1 m
for most bijective maps g, where ϕ0 “ m´1 m
`“1 ϕ
`“1 ϕ
r L may vary dramatically among different
In some cases, the values of D
parametrizations in finite samples; see § 4.
Third, the estimator rrL involves the estimators of ψ under H0 , i.e., ψp0` and
ψ 0 . When H0 is not true, they are not consistent for ψ. As a result, rrL is
no longer consistent for rm . A serious consequence is that the power of the
r L is not guaranteed to monotonically increase as H1 moves
test statistic D
away from H0 . Indeed our simulations (see § 3.2) show that under certain
parametrizations, the power may nearly vanish for obviously false H0 .
Fourth, in order to compute drL in (1.9), users need to have access to
r L : pX, ψ0 , ψq ÞÑ dL pψ0 , ψ | Xq. But, in most software
the LRT function D
p
packages, the LRT function is built as DL : X ÞÑ dL pψp0 pXq, ψpXq
| Xq.
r L.
Hence users would need to write themselves a subroutine for evaluating D
This may not be feasible for users because of lack of information or skills.

8

K.W. CHAN AND X.-L. MENG

In short, four problems need to be resolved: (i) lack of non-negativity,
(ii) lack of invariance, (iii) lack of consistency and power, and (iv) lack of
a computationally feasible algorithm. Problems (i)–(iii) are resolved in § 2
below, where § 2.1 presents an invariant combining rule, which fully resolves
(ii). Next, we propose two estimators of rm (or equivalently r) in § 2.2 and
§ 2.4. The first one is an always non-negative estimator, which partially
resolves (i); whereas the second one is a robust and always non-negative
estimator, which completely resolves (i) and (iii). Finally, in § 3, we derive
a very handy algorithm, which resolves (iv).
2. Improved MI Likelihood Ratio Tests.
2.1. An Invariant Combining Rule. To derive a MI LRT that is invariant
to re-parametrization, we replace drL by an asymptotically equivalent version
that behaves like a standard LRT statistic, that is, the average of the m
completed-data log-likelihoods. Specifically, let
(2.1)

Lpψq “

m
1 ÿ `
L pψq,
m `“1

where

L` pψq “ log f pX ` | ψq.

We emphasize that Lpψq is not a real log-likelihood (even if we drop the divider m), because it is not derived from an authentic sampling or imputation
model about the collection of all completed data sets: X “ tX 1 , . . . , X m u.
For one thing, all X ` s share the same Xobs , and hence if we want to write
an authentic likelihood function for X, we need to take into account the
(strong) dependence among the X ` s. Nevertheless, Lpψq can be treated as
a log-likelihood for computational purposes. In particular, we can obtain its
constrained and unconstrained MLEs of ψ via
(2.2) ψp0˚ “ ψp0˚ pXq “

arg max Lpψq,

ψp˚ “ ψp˚ pXq “ arg max Lpψq,

ψPΨ : θpψq“θ0

ψPΨ

where the asterisks “˚” both symbolize that they are different from ψp0 and
p and that they are MLEs from an artificially constructed “likelihood”. The
ψ,
corresponding log-likelihood ratio test statistic is given by
(2.3)

m
!
)
1 ÿ
dpL “ 2 Lpψp˚ q ´ Lpψp0˚ q “
dL pψp0˚ , ψp˚ | X ` q.
m `“1

Thus, in contrast to drL of (1.9), dpL aggregates MI datasets through averaging
MI LRT functions as in (2.1), rather than averaging MI test statistics and
pr
pr
moments, as in (1.8). Although ψp0˚ ´ ψ 0 Ñ 0 and ψp˚ ´ ψ Ñ 0 as n Ñ 8

MULTIPLE IMPUTATION LR TESTS

9

for each m, only dpL , not drL , is guaranteed to be non-negative and invariant
to parametrization of ψ for all m, n. Using these properties, we can derive
some important properties of dpL . But to do so, we need the usual RCs on
MLE as well as on MI.
Assumption 2.

The sampling model f pX | ψq satisfies the following:

(a) ψ ÞÑ Lpψq “ n´1 log f pX | ψq is twice continuously differentiable;
p
(b) the complete data MLE ψpXq
is the unique solution of BLpψq{Bψ “ 0;
2
(c) if Ipψq “ ´B Lpψq{BψBψ | , then for each ψ, there exists a positive
pr
definite matrix Ipψq “ Uψ´1 such that Ipψq Ñ Ipψq as n Ñ 8; and
(d) the observed MLE ψpobs of ψ obeys
(2.4)

”
´
¯
ı
´1{2 p
Tψ
ψobs ´ ψ ‹ | ψ ‹ ñ Nh p0, Ih q

as n Ñ 8, where Ih is the h ˆ h identity matrix.
Assumption 3. The imputation model is proper (Rubin, 2004):
”
´
¯
ı
´1{2 p`
(2.5)
Bψ
ψ ´ ψpobs | Xobs ñ Nh p0, Ih q,
”
¯
ı
”
´
ı
pr
pr
Tψ´1 Uψ` ´ Uψ | Xobs Ñ 0,
Tψ´1 pBψ ´ Bψ q | Xobs Ñ 0
(2.6)
independently for ` “ 1, . . . , m, as n Ñ 8, provided that Bψ´1 is well-defined.
Assumption 2 holds under the usual RCs that guarantee normality and
1 , . . . , X`
consistency of MLEs. Note that, if the imputed missing data Xmis
mis
are drawn independently from the posterior predictive distribution f pXmis |
Xobs q with a suitably chosen prior f pψq, then Assumption 3 is satisfied. Usually, the condition (2.6) is known as the lower order variability assumption of
Uψ` and Bψ . Clearly, we can replace ψ by its sub-vector θ in Assumptions 2
and 3. These θ-version assumptions are sufficient to guarantee the validity of
the following Theorem 2.4 and Corollary 2.3. For simplicity, Assumption 1,
the θ-version of Assumptions 2 and 3, and conditions that are strong enough
to guarantee Property 1.1 are collectively written as RCθ . These conditions
are commonly assumed for MI inference.
Theorem 2.1. Assume RCθ . Under H0 , we have (i) dpL ě 0 for all m, n;
(ii) dpL is invariant to parametrization of ψ for all m, n; and (iii) dpL l drL
as n Ñ 8 for each m.

10

K.W. CHAN AND X.-L. MENG

Consequently, an improved combining rule is defined as
(2.7)

p L prm q “
D

dpL
,
kp1 ` rm q

for a given value of rm . It follows the forms of (1.5) and of (1.10). The issue
then is how to estimate rm that avoids the defects of rrL of (1.9).
2.2. An Improved Estimator of rm . Using dpL in (2.3), we can modify rrL
in (1.9) to a potentially better estimator:
(2.8)

rpL “

m`1
pdL ´ dpL q.
kpm ´ 1q

Although dpL ě 0 is guaranteed by our construction, rpL ě 0 does not hold in
general for a finite m. However, it is guaranteed in the following situation.
Proposition 2.2. Write ψ “ pθ| , η | q| , where η represents a nuisance
parameter that is distinct from θ. If there exist functions L: and L; such
that, for all X, the log-likelihood function Lpψ | Xq “ log f pX | ψq is of the
form Lpψ | Xq “ L: pθ | Xq ` L; pη | Xq, then rpL ě 0 for all m, n.
The condition in Proposition 2.2 means that the likelihood function of ψ
is separable. Then, the profile likelihood estimator of η given θ, i.e., ηpθ “
arg maxη Lpθ, η | Xq, does not depend on θ. Trivially, if there is no nuisance
parameter η, the separation condition is satisfied. More generally, we have
pr

Corollary 2.3. Assume RCθ . We have (i) under H0 , rpL Ñ r as
pr
m, n Ñ 8; and (ii) under H1 , rpL Ñ r0 as m, n Ñ 8, where r0 ě 0 is
some finite value depending on θ0 and θ‹ .
Corollary 2.3 ensures that, under H0 , rpL is non-negative asymptotically
and converges in probability to the true r. But it also reveals another fundamental defect of rpL : under H1 , the limit r0 may not equal to r, a problem we
pr
will address in § 2.2. Fortunately, since dpL Ñ 8 under H1 , the LRT statistic
p L pp
D
rL q is still powerful, albeit the power may be somewhat reduced. Similarly, rrL of (1.9) has the same asymptotic properties and defects, but rpL
behaves more nicely than rrL for finite m. This hinges closely on the high
sensitivity of rrL to the parametrization of ψ for small m, e.g., in some cases,
rrL becomes more negative as H1 moves away from H1 ; see § 4.1.
Whereas practically we can fix the occasional negativeness of rpL by using
rpL` “ maxp0, rpL q, theoretically such an ad hoc fix misses the opportunity

11

MULTIPLE IMPUTATION LR TESTS
Table 1
The definitions of hypotheses H00 , H01 , H10 , H11 .
C0 : ψ1 “ ¨ ¨ ¨ “ ψm P Ψ
(i.e., r “ 0)

C1 : ψ1 , . . . , ψm P Ψ
(i.e., r ě 0)

C0 : θ1 “ ¨ ¨ ¨ “ θm “ θ0 P Θ
(i.e., H0 -constrained)

H00 “ C0 X C 0

H01 “ C0 X C 1

C1 : θ1 , . . . , θm P Θ
(i.e., not H0 -constrained)

H10 “ C1 X C 0

H11 “ C1 X C 1

H00

H01

H10

H11

Fig 1. The relationships between the four hypotheses H00 , H01 , H10 , H11 . Each arrow denotes
an implication, e.g., H00 ñ H01 means that H00 implies H01 .

to improve upon rpL , which cannot possibly be the best estimator precisely
because it permits negative values, a sign of a structural defect. More importantly, the inconsistency of rpL under H1 cannot be resolved trivially.
2.3. A Complication Caused by Nuisance Parameter. To better understand the source of the negativity of rpL , we extend Lpψq in (2.1) to allow it
take m different arguments:
(2.9)

Lpψ 1 , . . . , ψ m q “

m
1 ÿ ` `
L pψ q.
m `“1

Using the “log-likelihood” Lpψ 1 , . . . , ψ m q, we can construct, at least conceptually, four hypotheses H00 , H01 , H10 , H11 defined in Table 1. Each of them
consists of zero, one or two of the constraints
C0 : θ1 “ ¨ ¨ ¨ “ θm “ θ0

and

C0 : ψ1 “ ¨ ¨ ¨ “ ψm.

The constraint C0 is equivalent to H0 , and the constraint C 0 means that all
ψ ` s are equal, and hence it is effectively equivalent to r “ 0, i.e., no missing
information. The relationships among H00 , H01 , H10 , H11 can be visualized
in Figure 1. Define the maximized value of Lpψ 1 , . . . , ψ m q under hypothesis
H P tH00 , H01 , H10 , H11 u by LpHq. Then we can re-express pdL ´ dpL q{2 as
(
(
(2.10)
pdL ´ dpL q{2 “ LpH11 q ´ LpH10 q ´ LpH01 q ´ LpH00 q .

12

K.W. CHAN AND X.-L. MENG

IL (H 00)

IL (H 01)

IL (H 10)

IL (H 11)

Constraint

0

Constraint

0

Fig 2. A schematic illustration of the sign of (2.10). The contour lines of Lpψ 1 , . . . , ψ m q
are plotted. The two straight lines refer to constraints C0 and C 0 . Since LpH11 q( “
0.082, LpH01 q “ LpH10 q “ 0.08, and LpH00 q “ 0.01, we have LpH11 q ´ LpH10 q ´
(
LpH01 q ´ LpH00 q “ 0.002 ´ 0.007 ă 0. Note that the function Lpψ 1 , . . . , ψ m q in (2.9)
is at least 4-dimensional (i.e., θ1 , θ2 , η 1 , η 2 ) generally, so the above illustration in a 2dimension space is just conceptual.

Whereas the two bracketed terms in (2.10) are non-negative because they
correspond to two LRT statistics, the difference between these two terms is
not guaranteed to be non-negative. A simple example illustrates this well.
For the regression model rY | X1 , X2 s „ Npβ0 ` β1 X1 ` β2 X2 , σ 2 q, the
LRT statistic for testing H10 : β1 “ 0, β2 P R against H11 : β1 , β2 P R is not
necessarily larger (or smaller) than that for testing H00 : β1 “ β2 “ 0 against
H01 : β1 P R, β2 “ 0. A schematic illustration is provided in Figure 2.
The decomposition (2.10) provides another interpretation of rpL . The test
statistic LpH11 q ´ LpH10 q seeks evidence for detecting the falsity of r “ 0 in
both θ and η, whereas LpH01 q ´ LpH00 q seeks evidence only in η. For cases
where θ and η are orthogonal (at least locally), the left-hand side of (2.10)
can be viewed as a measure of evidence against r “ 0 solely from θ; Proposition 2.2 already hinted this possibility. However, the “test statistic” (2.10)
has a very serious problem apart from being possibly negative. Because C0
requires all θ` s to coincide with a specific θ0 , C0 is not nested within C 0 , i.e.,
C 0 œ C0 . Hence rpL is guaranteed to consistently estimate rm only under
H0 . This explains Corollary 2.3, and leads to an improvement below.
2.4. A Consistent and Non-negative Estimator of rm . Our new estimator simply drops the second term in (2.10), that is, we estimate rm by
(2.11)
(2.12)

m`1
pδ L ´ δpL q, where
hpm ´ 1q
δ L “ 2Lpψp1 , . . . , ψpm q,
δpL “ 2Lpψp˚ , . . . , ψp˚ q,
rpL♦ “

MULTIPLE IMPUTATION LR TESTS

13

where h is the dimension of ψ, and the rhombus “♦” symbolizes a robust
estimator. It is robust, because it is consistent under either H0 or H1 , as long
as we are willing to impose the EFMI assumption on the entire parameter
ψ, a stronger requirement than Assumption 1. This is necessary because the
Wald test involves only θ, whereas LRT has to deal with the entire ψ. The
collection of Assumptions 2–4 will be referred to as RCψ .
Assumption 4 (EFMI of ψ).

There is r ě 0 such that Tψ “ p1 ` rqUψ .

Theorem 2.4. Assume RCψ . Under both H0 and H1 , we have (i) rpL♦ ě 0
for all m, n; (ii) rpL♦ is invariant to parametrization of ψ for all m, n; and
pr
(iii) rpL♦ Ñ r as m, n Ñ 8, where r is given in Assumption 4.
p L prm q of (2.7) and improved estimaWith the improved combining rule D
tors for rm , we are ready to propose two MI LRT statistics:
p` “ D
p L pp
D
rL` q
L

(2.13)

and

p L pp
p♦ “ D
rL♦ q.
D
L

pL “ D
p L pp
rL q.
For comparison, we also study the test statistic D
2.5. Reference Null Distributions. The estimators rpL` and rrL have the
same functional form asymptotically (n Ñ 8) and rely on the same set of
assumptions, hence they have the same asymptotic distribution.
Lemma 2.5.

Suppose RCθ and m ą 1. Under H0 , we have, jointly,
rpL`
ñ M2
rm

(2.14)

and

p ` ñ p1 ` rm q M1
D
L
1 ` rm M2

as n Ñ 8, where M1 „ χ2k {k and M2 „ χ2kpm´1q {tkpm´1qu are independent.
p` “ D
p L pp
Consequently, the null distribution of D
rL` q can be approximated
L
by Fk,dfpp
r r` ,kq , but a better approximation will be provided shortly.
L

pr

For the other proposal, although rpL` ´ rpL♦ Ñ 0 as n Ñ 8 under H0 , their
non-degenerated distributions (after proper scaling) are different because rpL♦
relies on an average FMI in ψ, but rpL` only on an average FMI in θ.
Theorem 2.6.
(2.15)

Suppose RCψ and m ą 1. Under H0 and H1 , we have
χ2hpm´1q
rpL♦
ñ M3 „
rm
hpm ´ 1q

as n Ñ 8, where M3 is independent of the M1 defined in (2.14).

14

K.W. CHAN AND X.-L. MENG

Theorem 2.6 implies that, if n can be regarded as infinity and rpL♦ is uniformly integrable in L 2 , then
rL♦ q “ Epp
rL♦ q ´ rm “ 0
Biaspp

and

rL♦ q “
Varpp

2
2rm
“ Op1{mq
hpm ´ 1q

?
as m Ñ 8. Therefore rpL♦ is a m-consistent estimator of rm (or r) in L 2 .
Moreover, for each m ą 1 and as n Ñ 8, we have
rL` q
Biaspp
rL♦ q
Biaspp

Ñ1

and

Varpp
rL` q
rL♦ q
Varpp

Ñ

h
ě 1,
k

which implies that rpL♦ is no less efficient than rpL` when RCψ holds. This is of
no surprise because of the extra information in rpL♦ brought in by the stronger
Assumption 4. Result (2.15) also give us the exact (i.e., for any m ą 1, but
p ♦ , as given below.
assuming n Ñ 8) reference null distribution of D
L
Theorem 2.7.
(2.16)

Assume RCψ and m ą 1. Under H0 , we have
p ♦ ñ p1 ` rm q M1 ” D
D
L
1 ` rm M 3

as n Ñ 8, where M1 „ χ2k {k and M3 „ χ2hpm´1q {thpm´1qu are independent.
For a given rm , simulating the above reference null distribution is fast and
straightforward, so the quantile (i.e., the critical value for rejecting H0 ) can
be approximated by Monte Carlo methods. Moreover, according to (2.16)
p ♦ and D
p ` converge in distribution to M1 “ χ2 {k as
and (2.14), both D
k
L
L
pr
m, n Ñ 8, which is a pivotal distribution, i.e., free of rm . Because M3 Ñ 1
pr
p ♦ is expected to be more robust with respect
at least as fast as M2 Ñ 1, D
L
p ♦ over D
p `.
to rm . This is another advantage of using D
L
L
Nevertheless, a simple approximate null distribution might be handy in
practice. The right-hand side of (2.16) can be approximated by Fk,dfpr
r m ,hq
using the same rationale of Li et al. (1991b); see (1.7). But this approximation is poor (see below). A better approximation is by matching the first
two moments of the denominator of (2.16), 1 ` rm M3 , with that of a scaled
´1 q2 hpm ´ 1q, and hence
χ2 : aχ2b {b. This yields a “ 1 ` rm and b “ p1 ` rm
the right-hand side of (2.16) can be approximated by Fk,dfpr
p m ,hq , where
"
(2.17)

p m , hq “
dfpr

1 ` rm
rm

*2
hpm ´ 1q “

hpm ´ 1q
.
2
fm

15

MULTIPLE IMPUTATION LR TESTS

This degrees of freedom is very appealing because it can be viewed as the
degrees of freedom for estimating the nuisance parameter rm after being
inflated by the square of fm “ rm {p1 ` rm q, i.e., the finite imputation
2
corrected FMI. The smaller the fm , the closer Fd,dfpr
p m ,hq is to χk {k, which
is equivalent to the usual large-n asymptotic χ2 test; as mentioned before,
for small n, see Barnard and Rubin (1999).
To check how Fk,dfpr
p m ,hq compares with Fk,dfpr
r m ,hq as an approximation
to the distribution of D given in (2.16), we compute via simulations
!
)
!
)
r “ P D ą F ´1r
α
p1 ´ αq
and α
p “ P D ą F ´1p
p1 ´ αq ,
k,dfprm ,hq

k,dfprm ,hq

´1
where Fk,df
pqq denotes the q-quantile of Fk,df . We draw N “ 218 independent
copies D for each of the following possible combinations: m P t3, 5, 7u, k P
t1, 2, 4, 8u, τ “ h{k P t1, 2, 3u, fm P t0, 0.1, . . . , 0.9u, and following Benjamin
et al. (2017)’s recommendation, we use both α P t0.5%, 5%u. The result for
α “ 0.5% is shown in Figure 3, while that for α “ 5% is deferred to the
p approximates α much better than α
r, especially
Appendix. In general, α
when m, k, h are small. When m, h are larger, their performances are closer
to each other. It makes sense because both Fk,dfpr
r m ,hq and Fk,dfpr
p m ,hq are
2
essentially χk {k when hpm ´ 1q is large. Note also that the performances
of α
r and α
p are not monotonic in fm . Generally speaking, the performance
of Fk,dfpr
p m ,hq is particularly good for 0% À fm À 30%. Consequently, we
p ♦ , and
recommend using F p ♦ as an approximate null distribution for D
rL ,hq
k,dfpp

L

p`
Fk,dfpp
p r` ,kq for DL , as employed in the rest of this paper.
L

3. Computational Considerations and Comparisons. The statistic dL of (1.8) is easy to compute because only the standard complete-data
p
procedure DL : X ÞÑ dL pψpXq,
ψp0 pXq | Xq is needed. However, dpL of (2.3)
♦
and rpL of (2.8) in general cannot be computed solely by DL , e.g., dpL requires
ř
`
p˚ p˚
tX 1 , . . . , X m u ÞÑ m
`“1 dL pψ , ψ0 | X q. Creating a subroutine for this computation requires additional effort and information that may beyond a user’s
capacity. Here we show how to compute or approximate dpL and rpL♦ solely by
DL or an extremely trivial modification of DL .
3.1. Computationally Feasible Combining Rule. We begin with precise
notation for our complete data X and its sampling model f pX|ψq. For the
vast majority of real-world datasets, especially those in public domain, X is
of the form of an n ˆ p matrix, with rows indicating subjects and columns
variables/attributes. We then write X “ pX1 , . . . , Xn q| , and its sampling

16

K.W. CHAN AND X.-L. MENG
Nominal size α = 0.5 %
k=1

k=8

k=4

k=2

FMI
10% 30% 50% 70% 90% 10% 30% 50% 70% 90% 10% 30% 50% 70% 90% 10% 30% 50% 70% 90%

m=3

Empirical size

0.5%
0.4%
0.3%
0.2%
0.1%
0%
0.5%

m=5

0.4%
0.3%
0.2%

3
3
1 2
2
2
3 3
1 2
2 2 3
3
3 1
2
3 3 3
2
2
1
2
2
3
2 2
2
3
1
3
2 3
3 3 3
2
1
2
1 1
1

3 3
2
2
3
1 3
2
1
3
1 1 2 3
2 3
2
1 3
1
3
3
2
1
1
2
3
2
2
1

2
3
3
1 1
2
3 2
3 3 3 3 3
2 2
2 1
3 3
3
2 2 2
1
3
2
3
3
2
1 3 3 3
2
1
2
1
2
2
2 2 1
1

2 2
3
3
1
1 3
3
2
1 2
2
3 2
3 3
2
1
3 3
3 3
2
1
2
1
3
2
2
1
2
1

2
1

1

3 3 3 3
2 2
2
2 3
3 3 3
1

2 1
1 2
2 2 1 1
1
1
1 1 1

1

3
1 3
2
3 3 3 3 3 3
3 3
1
1 2
2
2
3 3
3
2
2 3
3
2 3
1 2 1
2 1
2 3
1
2
2 2
2 1
2 1
2 2
1
1 1
2 1
1
1
1
1 1 1
1
3
2
3
2
1

1

1
1 1 1 1

3 3 3 3 3
1 3
1 1 1 2
2
3
1
2 2
3
1
1
2
3 3
3
3
1 3
2 2
2 3
3 3
2 2
2 2
2
2
2
1
2 1
1
2
1 1
2 1
1
1

1
1

2 2
1 2 2
2
3 3
1 3
2
3
2
3
1 3
3 3 3
3 3 3
2 2
2 2
1
3
2
2 3 3 3 3 3
1
1
2
2
2
1 2 2
2
1 1 1

2
2
2
2 3
1 3
2
3 3
2
2
2
2 3
2
2
1 3
1
2
1 3
2 3
3
3
2 3
1 3
3
2
1 3
2 2
2
3
3 3
1 1
1
1
1
1 1
1 2
1
2 2
1 2
1 1 1 1 1
1
1
1
1 1

Approximations
LRR (91)
New proposal

1
1 1 1

0.1%
0%
0.5%

m=7

0.4%
0.3%
0.2%
0.1%

2
2
1 3
2 3
2
1 3
2 3
2
1
1
3
1 3
2
2
1 3
2
2
3
2
2
1
1 3
2
2 3
3 2
1 2
2 3
2 2
3 3
3 3
3
2 3
1
2 3
1
2 3
3
1 1
1
3 3 3 3 2
2 3
2 1
1
2
2 2
1 2
1
1
1 1
1
1

2 2
3 3 3
3
3
3 2
2 3
3 3
1
3
1 3 2 1 3
2
1 3
2
3
1 3
1
3
2
3
2
2 2
2
2
1
1 3
2
3
2
1 3
3 2
3
2 3
1 3
1
1 2
2
1 2
3
2
2
1
1
1
3
1 1
2
1 3
1 2
2
1 2
1 2
2
2
2
1 1
2
1 3
1 2
1 1
1 1 1
1 2
1 2
1
1
1 1 1 1
1
1
1 1

1
1
1

1
1 1

1

0%

Fig 3. The performances of two different approximate null distributions when the nomp or α
r, and the horizontal axis denotes
inal size is α “ 0.5%. The vertical axis denotes α
the value of fm . The number attached to each line denotes the value of τ “ h{k. The
p is denoted by thick solid lines with triangles, and the existing
proposed approximation α
r is denoted by thin dashed lines with circles.
approximation α

model by fn pX | ψq. Correspondingly, the `th completed-dataset by MI is
X ` “ pX1` , . . . , Xn` q| . Define the stacked dataset by X S “ rpX 1 q| , . . . , pX m q| s| ,
a matrix having mn rows, which is conceptually different from the collection
of datasets X “ tX 1 , . . . , X m u. Treating X S as a dataset with size mn, we
can define
S

(3.1)

L pψq “

1
log fmn pX S | ψq,
m

which, other than the scaling factor 1{m, is just the ordinary log-likelihood
function of ψ based on the dataset X S . Since the scaling factor 1{m does
not alter the maximizers, as long as the user’s complete-data procedure can
handle size mn instead of just n, the user can apply it to X S to obtain
(3.2)

ψp0S “

S

arg max L pψq
ψPΨ : θpψq“θ0

and

S
ψpS “ arg max L pψq.
ψPΨ

MULTIPLE IMPUTATION LR TESTS

17

Consequently, the quantities
(3.3)

S
δp0,S “ 2L pψp0S q

S
δpS “ 2L pψpS q

and

are readily available from the user’s complete-data procedure. It is then
S
desirable if we can replace Lpψq by L pψq in the proposed test statistics.
Precisely, in parallel to (2.7), (2.8) and (2.11), we define
(3.4)
(3.5)
(3.6)

dpS
,
kp1 ` rm q
m`1
rpS “
pdS ´ dpS q,
kpm ´ 1q
m`1
rpS♦ “
pδ S ´ δpS q,
hpm ´ 1q

p S prm q “
D

with dpS “ δpS ´ δp0,S of (3.3);
with dS “ dL of (1.8);
with δ S “ δ L of (2.12);

p ♦ and its existing
and rpS` “ maxp0, rpS q. The “stacked” counterparts of D
L
p L and D
p ` (see (2.13)) then are given by
counterparts D
L
(3.7)

p♦ “ D
p S pp
rS♦ q,
D
S

pS “ D
p S pp
rS q,
D

p` “ D
p S pp
rS` q.
D
S

Proposition 3.1. IfśX “ pX1 , . . . , Xn q| is row-independent for arbitrary n, i.e., f pX | ψq “ ni“1 f pXi | ψq, then (2.1) and (3.1) are the same:
S
pS ” D
p L and D
p♦ ” D
p ♦.
Lpψq ” L pψq. Consequently, D
S
L
Since for many applications, the rows correspond to individual subjects,
the row-independence assumption typically holds for arbitrary n (hence we
can extend from n to mn, assuming the user’s complete-data procedure is
not size-limited). Even if it does not hold, we can still have dpL l dpS under
S
some RCs that guarantee Lpψq and L pψq are close; see Appendix A, where
we also reveal a subtle but important difference between rpL♦ and rpS♦ .
Now define complete data functions (with data X being the only input)
(3.8)
(3.9)

p
DL,0 pXq “ 2 log f pX | ψp0 pXqq, DL,1 pXq “ 2 log f pX | ψpXqq,
p
DL pXq “ dL pψp0 pXq, ψpXq
| Xq.

The subroutine for evaluating the complete data LRT function X ÞÑ DL pXq
is usually available, as is the subroutine for the maximum of complete data
log-likelihood X ÞÑ DL,1 pXq; for example, the function logLik in R extracts
the maximum of complete data log-likelihood for objects belonging to classes
"glm", "lm", "nls" and "Arima". If the user is able to edit the subroutine for
X ÞÑ DL pXq “ DL,1 pXq ´ DL,0 pXq, a subroutine for X ÞÑ DL,1 pXq can be

18

K.W. CHAN AND X.-L. MENG

p`
Algorithm 1: MI LRT statistic D
S
Input: MI datasets X 1 , . . . , X m ; dimension k; function DL in (3.9).
begin
Compute dL “ m´1 tDL pX 1 q ` ¨ ¨ ¨ ` DL pX m qu .
(i) Stack the datasets to form X S “ rpX 1 q| , . . . , pX m q| s| .
(ii) Compute dpS “ m´1 DL pX S q.
p ` according to (3.4) and (3.7).
Calculate rpS` according to (3.5), and D
S
`
p r , kq according to (2.17).
Calculate dfpp
S
p`
The p-value is Sk,dfpp
x r ` ,kq pDS q, where Su,v p¨q is the survival function of Fu,v .
S

p♦
Algorithm 2: (Robust) MI LRT statistic D
S
Input: MI datasets X 1 , . . . , X m ; dimensions h, k; functions DL,1 , DL in (3.8), (3.9).
begin
Compute δ S “ m´1 tDL,1 pX 1 q ` ¨ ¨ ¨ ` DL,1 pX m qu .
(i) Stack the datasets to form X S “ rpX 1 q| , . . . , pX m q| s| .
(ii) Compute dpS “ m´1 DL pX S q and δpS “ m´1 DL,1 pX S q.
p ♦ according to (3.4) and (3.7).
Calculate rpS♦ according to (3.6), and D
S
♦
p
Calculate dfpp
rS , hq according to (2.17).
p♦
The p-value is Sk,dfpp
x r ♦ ,hq pDS q, where Su,v p¨q is the survival function of Fu,v .
S

obtained easily by dropping the second term. However, editing a subroutine
may not be possible in general.
Now, let the survival function of Fu,v be Su,v pxq “ PpFu,v ą xq. Then, we
p ` and D
p ♦ , respectively.
can use Algorithms 1 and 2 to compute D
S
S
3.2. Computational Comparison with Existing Tests. First, we list some
existing estimators of rm and their computation. Let s2W,1 and s2W,1{2 be the
‘ ` m
sample variances of td`W um
`“1 and t dW u`“1 , respectively. By the method of
moments concerning s2W,1 and s2W,1{2 , Rubin (2004) and Li et al. (1991a)
respectively proposed estimating rm by
p1 ` 1{mqs2W,1
c!

(3.10) rrW,1 “
2dW `

2
4dW

´ 2ks2W,1

)` ,

rrW,1{2 “ p1 ` 1{mqs2W,1{2 ,

where tau` “ maxp0, aq. Note that when k is large and m is small, using
(3.10) may lead to power loss, although the users have no choice when they
rL of (1.9), i.e., rrL` “
are given only td`W um
`“1 . A trivial modification of r
maxp0, rrL q, is a better alternative if the user is able to compute rrL .

19

MULTIPLE IMPUTATION LR TESTS

Second, we list some alternative MI combining rules. Having the above
estimators of rm , we can insert them into the following combining rules:
r 1 prm q “
D
W

dr1W
,
kp1 ` rm q

r L prm q “
D

drL
,
kp1 ` rm q

!
)`
r ` prm q “ D
r L prm q .
D
L

Using (1.3) and (1.9), we can also define the following combining rules:
1

(3.11)

1
DW prm q

dW ´ kpm´1q
m`1 rm
“
,
kp1 ` rm q

DL prm q “

dL ´ kpm´1q
m`1 rm
.
kp1 ` rm q

1

The combining rule DW prm q is useful when computing dr1W is difficult but
1
computing dW and estimating rm are simple. However, the resulting power
may deteriorate if the estimator of rm is inefficient or inaccurate. This type
of test statistics is also mentioned in Li et al. (1991a). Indeed, there are
infinitely many asymptotically equivalent test statistics, e.g., any convex
r W prm q ` p1 ´ φqDW prm q, for
r W prm q and DW prm q, i.e., φD
combination of D
φ P r0, 1s. When rrW,1 or rrW,1{2 is used for estimating rm , the null distributions of the resulting MI test statistics can be approximated by Fk,df
r 1 pr ,kq ,
m

r 1 prm , kq “ pm ´ 1qp1 ` r ´1 q2 k ´3{m ; see Li et al. (1991a).
where df
m
Next, we introduce some notation to facilitate the comparison. Define the
complete data moments estimation procedures (MW and ML ), and complete
data testing procedures (DW and DL ) by
!
)
!
)
p
p
MW pXq “ θpXq,
U pXq ,
ML pXq “ ψpXq,
ψp0 pXq ,
p
DW pXq “ dW pθpXq,
U pXqq,

p
DL pXq “ dL pψp0 pXq, ψpXq
| Xq.

Also define the following non-standard complete data procedures:
p
r L pX, ψ0 , ψq “ dL pψ0 , ψ | Xq, DL,1 pXq “ 2 log f pX | ψpXqq,
D
m
m
1 ÿ
2 ÿ
DL pXq “
dL pψp0˚ pXq, ψp˚ pXq | X ` q, DL,1 pXq “
log f pX ` | ψp˚ pXqq,
m `“1
m `“1
r L u prowhere X “ tX 1 , . . . , X m u. Clearly, MW produces DW , and tML , D
r L produces tML , DL , DL,1 u.
duces DL . If users can perform optimization, D
Note that an un-normalized density can be used in DL,1 and DL,1 pXq.
Now, we investigate the computational and statistical properties for differr 1 p¨q, D1 p¨q, D
r L p¨q, D
r ` p¨q, D
p S p¨qu
ent MI test statistics Dprq, where Dp¨q P tD
W
W
L
`
` ♦
1
1
1
rW , rrW,1 , rrW,1{2 , rrL , rrL , rpS , rpS , rpS u are generic notations for MI comand r P tr
bining rules (functions of rm ) and estimator of rm , respectively. The required
computer routines and the following properties are summarized in Table 2.

Combining Rule
Test

No.

Formula

Approx. null distribution

Estimator of rm

Properties

Original

Proposed

MW
MW

a
Fk,dfpr
Ă m ,kq
Fk,df
Ă 1 prm ,kq

Fk,dfpr
x m ,kq
NA

1
rrW
1
rrW,1

MW
DW

R04
R04

´
´

`b
´

`
`

´
´

´
´

´
´

`
`

θ
θ

MW
MW

Fk,df
Ă 1 prm ,kq
Fk,dfpr
Ă m ,kq

NA
Fk,dfpr
x m ,kq

1
rrW,1{2
1
rrW

DW
DW

LMRR91
R04

´
´

´
`

`
`

´
´

´
´

´
´

`
`

θ
θe

W-5 DW prm q

DW

Fk,df
Ă 1 prm ,kq

NA

1
rrW,1

DW

R04

´

´

´

´

´

`

`

θ

1
DW prm q

DW

Fk,df
Ă 1 prm ,kq

NA

1
rrW,1{2

DW

LMRR91

´

´

´

´

´

`

`

θ

LRT L-1
L-2

r L prm q
D
r ` prm q
D

rL
ML , D
rL
ML , D

Fk,dfpr
Ă m ,kq
Fk,dfpr
Ă m ,kq

Fk,dfpr
x m ,kq
Fk,dfpr
x m ,kq

rrL
rrL`

rL
ML , D
rL
ML , D

MR92
MR92g

´
´

´
´

´
`

´
´

`
`

´f
´

`
`

θ
θ

L-3
L-4
L-5

p S prm q
D
p S prm q
D
p S prm q
D

DL
DL
DL

Fk,dfpr
Ă m ,kq
Fk,dfpr
Ă m ,kq
Fk,dfpr
Ă m ,hq

Fk,dfpr
x m ,kq
Fk,dfpr
x m ,kq
Fk,dfpr
x m ,hq

rpS
rpS`
rpS♦

DL
DL
DL,1

Proposal
Proposal
Proposal

`
`
`

´
´
`

´
`
`

´
´
`

`
`
`

`
`
`

´h
´
´

θ
θ
ψ

W-4

DW pT q

d

1

W-6

a

L

Formula Routine Reference Inv Rob ě 0 Pow Def Sca Dep EFMI

The original approximate null distribution in Rubin (2004) was modified by Li et al. (1991a). This footnote also applies to W-2,4,5.
1
The estimator rrW
does not depend on θ0 , but its MSE may be inflated under H1 if a bad parametrization of θ is used.
1
1
c
The originally proposed combining rule is DW prm q; see (3.11). Although DW prm q is more computational feasible, the power loss
1
1
r
is more significant than DW prm q after inserting an inefficient estimator rrW,1 for rm . This footnote also applies to W-3.
d
Computing the test statistic DW pT q “ dW pθ, T q{k does not require estimating rm .
e
EFMI is not required for the test statistic DW pT q, but it is required for its approximate null distribution.
f
Averaging and processing vector estimators of ψ, but not their covariance matrixes, is needed. This footnote also applies to L-2.
g
It is a trivial modification of the original proposal in MR92 by replacing rrL with rrL` “ maxt0, rrL u.
h
p L p¨q and rpL are used instead of D
p S p¨q and rpS , then the resulting test can be applied to dependent data. Similarly, if we replace
If D
the subscript “S ” by “L ” in L-4,5, then the resulting tests can also be applied to dependent data.
b

K.W. CHAN AND X.-L. MENG

Routine

1
rW
prm q
WT W-1 D
1
r
W-2 DW prm qc
1
rW
prm q
W-3 D

20

Table 2
Computational requirements and statistical properties of MI test statistics, their associated combining rules and estimators of FMI rm .
The symbols “`” and “´” mean that the test statistic (or estimator) is equipped and not equipped with the indicated property,
respectively; see the end of § 3.2 for heading descriptions. The reference papers/book are abbreviated as follows: Rubin (2004) (R04), Li
et al. (1991a) (LMRR91) and Meng and Rubin (1992) (MR92).

21

MULTIPLE IMPUTATION LR TESTS
Table 3
The values of parameters used in the simulation experiment in § 4.1.
Parameters
Experiment

Fixed

Variable

No.

Varying

ρ

p

f

Case 1

Case 2

Case 3

Case 4

Case 5

I
II
III

Correlation ρ
Dimension p
FMI f

–
0.4
0.4

2
–
2

0.5
0.5
–

´0.8
2
0.1

´0.4
3
0.3

0
4
0.5

0.4
5
0.7

0.8
6
0.9

•
•
•
•
•
•
•
•

(Inv) Dprq and r are invariant to re-parametrization of ψ;
(Rob) r is robust against θ0 , i.e., consistent under both H0 and H1 ;
(ě 0) Dprq and r are non-negative for all m and n;
(Pow) the test has high power to reject H0 under H1 ;
(Def) Dprq and r are always well-defined and numerically well-conditioned;
(Sca) the MI procedure requires users only to deal with scalars;
(Dep) X1 , . . . , Xn can be dependent; and
(EFMI) whether EFMI is assumed for θ or for ψ.

p S pp
rS` q is the most computationally attractive test statistic.
In summary, D
p S pp
If the user is willing to make stronger assumptions, D
rS♦ q has better statisp L pp
tical properties, and is still computationally feasible. Nevertheless, D
rL♦ q
is the most general test statistic and has the best statistical properties.
4. Empirical Investigation and Findings.
4.1. Monte Carlo Simulation Experiments. Suppose that X1 , . . . , Xn „
Np pµ, Σq independently, where Σ “ σ 2 tp1 ´ ρqIp ` ρ1p 1|p u, and 1p is the
p-vector of ones. The values of p, σ 2 , ρ and µ are specified below. Further
assume that only nobs “ tp1´fqnu data are observed, where f P p0, 1q is the
FMI. Let Xobs “ tXi : i “ 1, . . . , nobs u and Xmis “ tXi : i “ nobs `1, . . . , nu.
Suppose that we want to test whether the means of all components are equal,
i.e., H0 : µ “ µ0 1p , where µ0 P R is an unknown constant.
Obviously, one may directly use the observed dataset to construct the
LRT statistic DL without MI. The test DL (denoted by L-0) is regarded as
a benchmark for comparison. Throughout this subsection, W-1,2,3,4 and L1,2,3,4,5 stated in Table 2 are compared. In the imputation step, a Bayesian
model is employed for imputation. Assume a multivariate Jeffreys prior on
pµ, Σq, i.e., f pµ, Σq9|Σ|´pp`1q{2 . Then, the `th imputed missing dataset can
be produced by the following procedure, for ` “ 1, . . . , m.
Step 1: Draw a posterior sample Σ` from the inverse-Wishart distribution

22

K.W. CHAN AND X.-L. MENG

´1
withřn ´ 1 degree of freedom and scale matrix Sobs
, where Sobs “
nobs
´1
´1 řnobs
|
nobs i“1 pXi ´ X obs qpXi ´ X obs q and X obs “ nobs i“1 Xi .
Step 2: Draw one posterior sample µ` from Np pX obs , Σ` {nobs q.
Step 3: Draw pn ´ nobs q imputed missing values tXi` : i “ nobs ` 1, . . . , nu from
Np pµ` , Σ` q independently. Also, denote Xi` “ Xi for i “ 1, . . . , nobs .

With the `th completed dataset, the unconstrained MLEs for µ and Σ are
µ
p` “

n
1ÿ `
X ,
n i“1 i

n ´
¯|
¯´
ÿ
p` “ 1
p` .
p` Xi` ´ µ
Σ
Xi` ´ µ
n i“1

Whereas we generate data using a covariance matrix with common variance
and correlation, our model does not assume any structure for Σ. The only
restriction we can impose is the common-mean assumption under the null,
for which the constrained MLEs are
#
+
| p ` ´1 `
´
¯´
¯|
µ
p
1
p
Σ
q
p
`
`
`
`
p` “ Σ
p` ` µ
µ
p`0 “
p
µ
p
µ
p
µ
p
1
,
Σ
´
´
.
p
0
0
0
p ` q´1 1p
1|p pΣ
In the experiment, we study the impact of parametrization on different
test statistics. For the Wald tests, three parametrizations of θ are examined:
(i) θ “ pµ2 ´ µ1 , . . . , µp ´ µp´1 q| — differences of means,
(ii) θ “ `pµ2 {µ1 , . . . , µp {µp´1 q| ´˘1p´1 — relative differences of means, and
|
(iii) θ “ µ32 ´ µ31 , . . . , µ3p ´ µ3p´1 — differences of cubic means.
For any case above, H0 can be re-expressed as θ0 “ 0p´1 , an pp ´ 1q-vector
of zeros. For LRTs, the following parametrizations of ψ are used:
(i) ψ “ tµ; Σu — means and covariances,
?
(ii) ψ “ t σii {µi , 1 ď i(ď p; Σu — noise-to-signal and covariances, and
(iii) ψ “ µ| Σ´1{2 ; Σ´1 — standardized means and precisions,
where Σ “ pσij q and Σ1{2 is the symmetric square root of Σ. The dimension
of ψ is h “ pp2 ` 3pq{2.
In the first part of the experiment, we study the distribution of p-values
derived from each test under H0 . In particular, we use n “ 100, m “ 3,
σ 2 “ 5 and µ “ 1p , with various values of ρ, p and f specified in Table 3. All
simulations are repeated 212 times. The comparison under parametrization
(ii) is shown in Figure 4; whereas those under parametrizations (i) and (iii)
are deferred to Appendix C. Note that, for Wald tests under parametrization (ii), the matrix U ` is numerically singular in some of the replications
(less than 0.25%). Those ill-conditioned cases are removed from the analysis
(which should favor the Wald tests).

23

MULTIPLE IMPUTATION LR TESTS

Table 4
r L . The results under parametrizations (ii)
The empirical proportions of negative rrL and D
r L ě 0 in the experiments.
and (iii) are shown. For parametrizations (i), rrL ě 0 and D
Case
1

2

3

4

5

1

% of rrL ă 0

2

3

4

5

Experiment

Parametrization

rL ă 0
% of D

I

(ii)
(iii)

1
6

2
6

3
7

4
7

5
7

26
1

16
1

13
1

12
1

12
2

II

(ii)
(iii)

4
7

1
3

0
1

0
1

0
1

12
1

5
0

3
0

4
0

3
0

III

(ii)
(iii)

13
18

6
9

4
7

4
5

3
4

55
20

25
5

12
1

5
1

2
0

The empirical sizes of the MI Wald tests generally deviate from the nominal size α under parametrization (ii). In contrast, the sizes of all LRTs
are closer to α. However, the original L-1 and its trivial modification L-2
do not have accurate sizes when |ρ| or f is large. They can be over-sized
or under-sized depending on which parametrization is used. Moreover, the
trivial modification L-2 does not help to correct the size, and it may even
worsen the test. For our first proposed test statistics L-3 and L-4, they are
invariant to parametrizations and have quite accurate sizes, although they
are under-sized in challenging cases where both p and f are large. Moreover,
they are exactly the same throughout our simulation experiments, i.e., we
never observed rpL ă 0. For our second proposed test statistic L-5, it gives
the most satisfactory overall results. It generally has very accurate size, except that it is slightly over-sized for large p, a problem that should diminish
when we use m beyond the smallest recommended m “ 3.
Interestingly, as seen clearly in Figure 4, the benchmark L-0 performs very
badly for large p and f. The sample size per parameter, n{h, is small; for
p ě 4, n{h ď 100{14 ă 8. The asymptotic null distribution χ2k {k then can
fail badly under arbitrary or even all parametrizations; (ii) apparently falls
into this category. An F approximation would be more appropriate. But this
is exactly what is being used for MI tests, albeit with different choices of the
r L and
denominator degrees of freedom. Table 4 documents how often rrL , D
rpS are negative. In some cases, nearly half of the simulated values of rrL and
r L are negative. In contrast, rpS is always non-negative in our simulation,
D
despite the fact that it can be negative in theory.
To study the power of each test, we set f “ 0.5, p “ 2, ρ “ 0.8,

24

K.W. CHAN AND X.-L. MENG

Case 2

Case 3

Case 4

Nominal size α
1% 3% 5% 7% 9%

1% 3% 5% 7% 9%

1% 3% 5% 7% 9%

1% 3% 5% 7% 9%

9%
Empirical size

Experiment I ( Vary ρ )

Parametrization (ii)
Case 1 ( Small ρ , p , f )

7%
5%
3%

m
M

1%

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

M
m

m
M
M
m

m
M
M
m
m
M

M
m

M
m

M
m
M
m

M
m

m
M
m
M

m
M

Case 5 ( Large ρ , p , f )

M
m

M
m

M
m

M
m

M
m

M
m

M
m

M
m

m
M

M
M m
M m
M m
M m
m

M
M m
m

1% 3% 5% 7% 9%

M
M m
m

M
m

m
M

M
M m
M m
M m
M m
M m
M m
M m
M m
m

m
M

Experiment II ( Vary p )
Experiment III ( Vary f )

Test statistics
9%

9%

7%
5%
3%

M
m

1%

M
M m
M m
M m
m
M
m

M
M m
m

M
M m
m
M
m

M
m

M
m
M
m

M
m

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

M
m

m
M

m
M
m
M

m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M

3%

m
M

M
m

m
M

m
M

m
M

m
M

M
m

M
m

m
M

5%

m
M

M
m

M
m

7%

M
m

m
M

M
m
m
M
m
M
m
M

m
M

m
M

M
m

M
m

M
m

m
M

1%

M
m

M
M m
M m
M m
M m
m

M
M m
m

m
M
m
M
m
M
M M
m
M
M M
m
M
M
m
M M M
m M
M

m
M

M
M m
m

M

M

M

m

m

m

M

M

m m

M

m

M

m

M

m

M

m

M

m
M

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m

m
m
M

m
M

m

m m
m m
m m
m m

m
M

Fig 4. The comparison between empirical size and nominal size α under parametrization
(ii) for α P p0, 10%s. The Wald tests (W-1,2,3,4) and LRTs (L-0,1,2,3,4,5) are represented
r L (L-1: MR92)
by grey dashed and black solid lines, respectively. The LRT statistic D
r ` (L-2: MR92(+)) are the tests that greatly improved upon by our
and its modification D
L
p S (L-3: Proposal), D
p ` (L-4: Proposal(+)) and D
p ♦ (L-5: Proposal(R)).
proposals D
S
S

σ 2 “ 5 and µ “ p´2 ` δ, ´2 ` 2δq| for different values of m P t3, 10, 30u,
n P t100, 400, 1600u and δ “ µ2 ´ µ1 P r0, 4s. The empirical power functions
for size 0.5% tests under parametrizations (i), (ii) and (iii) are plotted in Figure 5. The results for size 5% tests are deferred to Table 11 of the Appendix.
Generally, none of the Wald tests exhibits monotonically increasing power as
δ increases, and their performance is affected significantly by parametrization. In particular, the powers can be as low as zero when 1 À δ À 2 under
parametrizations (ii) and (iii). Under parametrization (ii), L-1 is not powerful even for large δ. Moreover, its trivial modifications L-2 cannot retrieve
all the power it should have. On the other hand, our first proposed test
statistics L-3 and L-4 perform better than L-1 and L-2 at least for large m,
however, they also lose a significant amount of power when m is small.
Compared with all these, our second proposed test statistic L-5 performs
extremely well for all m and n, with power very close to the benchmark L-0
even for small m. To ensure the comparisons of power are fair, we also investigate the empirical (actual) size, α
p, in comparison to the nominal type-I
error α. Table 5 shows the minimum and maximum of the empirical sizes
over the three parametrizations considered in each test — and only one value

25

MULTIPLE IMPUTATION LR TESTS

Parametrization (i)
δ = µ2 − µ1
1
2
3

m = 30 , n = 100 m = 10 , n = 100 m = 3 , n = 100

Power

m = 3 , n = 400 m = 3 , n = 1600

0
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%

m
M

Nominal size α = 0.5 %
Parametrization (ii)

4 0

1

mm
mM
m
mm
Mm
Mm
MM
Mm
Mm
mm
MM
MM
Mm
Mm
Mm
Mm
mM
mM
M

MM
M
mmm

m
M

M
m

2
MMMM

3

Parametrization (iii)
4 0

M
Mmmmm
m
mmmmm

mmmmmm

m
M

MM
MMM
MM M
mmM
M
mmmm
m
m
mmmmmmmmmmm

m
M

m
mM
M

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M
mm
mM
MM
mm
mM
M
m
mM
M

m
M

MMMMMM
MM
MM
MMMMM
MMM
m
m
m
m
mmmmm mmmmmm
Mm mm
m

mM
mM
mM
mM
m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M

m
M
M
m

m
M
m
mM
M

m
M

m
mM
M

M
m

M

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
M

MMMMMMM
MM
MM
MMMmmmmmmmmmmm
Mmmm
M
M
mmm

m
M
m
M
m
M

m
M
m
mM
M

MMMMMMMM
mM
mM
M
mM
mM
mM
mM
mmmmmmmmm

m
M

m

m
M

M
m

m
M

m
mM
M

M
m

MMMMMMMMMMMMMM
MM
M
mmmmmmmmmmmmmmmm
Mm
m

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m
mM
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
M

m
M

Test statistics

MMMMM
MMM
MM
MMM
M
mmmmmmmmmm
m
m
M m
M m
M
mmm

m
M

m
M

M
mm
M

M
m

M MMMMMMM
MMMMMM M
MM
Mmmmmmmmmmmmmmmmmm
m

m
M

m
M

m
M

4

M
m

M
m

m
M

3

m
M

MMMMMMMMM

mM
mM
mM
mM
mM
mM
mM
M
mM
mM
mM
mM
mm
Mm
mM
MM
mM
m

m
M

2

m
M

m
M

m
M

1

MMMMMMMMMMMMMMM
MMM
M
mmmmmmmmmmmmmmmmmmm

MMMMMMMMMM

m
mM
M

Fig 5. The power curves under different parametrizations. The nominal size is α “ 0.5%.
In each plot, the vertical axis denotes the power, whereas the horizontal axis denotes the
value of δ “ µ2 ´ µ1 . The legend in Figure 4 also applies here.

is needed for those tests that are invariant to parametrization — when the
nominal size α “ 0.5% (results for α “ 5% are in Table 8 of the Appendix).
We see the deviations from the nominal α can be noticeable, especially
when m “ 3. To take that into account, we report the empirical size adα, which also has the interpretation as (an
justed power, that is, O “ power{p
approximated) posterior odds of H1 to H0 (Bayarri et al., 2016). Figure 6
plots the result for size 0.5% tests (the results for size 5% tests are in Figure
12 of the Appendix). Compared with the benchmark L-0, the odds O of the
proposed robust MI test (L-5) is closer to the nominal value 1{α as δ Ñ 8.
Nevertheless, the finite sample performances of all size 0.5%-tests are less
satisfactory than size 5%-tests because a very large sample size n is required
in order to approximate the tail behavior of test statistics satisfactorily.
We also investigate the performance of different estimators of rm for different values of δ and under different parametrizations. In our experiment,
we have rm “ 1 ` 1{m because we have set r “ 1. The mean square error
(MSE) of estimators fp “ rp{p1 ` rpq of fm “ rm {p1 ` rm q are computed.
The results are shown in Figure 7. Clearly, the only estimator that is con-

26

K.W. CHAN AND X.-L. MENG

Table 5
The range of empirical size rmin α
p, max α
ps in percentage, where max and min are taken
over the three parametrizations. Only one value is recorded for those tests that are
invariant to parametrization. The nominal size is α “ 0.5%.
p, max α
ps{%
Range of empirical size: rmin α
pn, mq

p1600, 3q

p400, 3q

p100, 3q

p100, 10q

p100, 30q

W-1
W-2
W-3
W-4

r0.90, 1.05s
r0.90, 1.05s
r0.98, 1.05s
r0.90, 1.05s

r0.76, 1.05s
r0.98, 1.22s
r0.98, 1.25s
r0.76, 1.05s

r0.20, 1.22s
r0.93, 1.25s
r0.90, 1.29s
r0.20, 1.22s

r0.07, 0.56s
r0.32, 0.73s
r0.34, 0.71s
r0.07, 0.56s

r0.02, 0.49s
r0.20, 0.85s
r0.22, 0.73s
r0.02, 0.49s

L-1
L-2
L-3
L-4
L-5
L-0

r0.90, 1.03s
r0.90, 1.05s
0.90
0.90
0.46
0.39

r1.10, 1.64s
r1.10, 1.76s
1.10
1.10
0.44
0.66

r1.15, 1.49s
r1.15, 2.37s
0.83
0.83
0.68
0.66

r0.37, 1.05s
r0.37, 0.98s
0.24
0.24
0.46
0.66

r0.10, 0.46s
r0.10, 0.49s
0.07
0.07
0.42
0.66

sistent, invariant to parametrization and robust against δ is our proposal
fpL♦ “ rpL♦ {p1 ` rpL♦ q. It concentrates at the true value fm quite closely even for
small m and n. Since fpL♦ is the only reliable estimator of fm , it verifies why L5 has the greatest power. On the other hand, the estimator frL “ rrL {p1 ` rrL q
has very large MSE when δ ­“ 0. It also explains why L-1 is not powerful.
4.2. Applications to a Care-Survival Data. Meng and Rubin (1992) applied their test to the data given in Table 6, where i, j and k index, respectively, clinic (A or B), amount of parental care (more or less) and survival
status (died or survived). However, the clinic label k is missing for some of
the observations (and the missing-data mechanism was assumed to be ignorable). Two hypotheses were tested in Meng and Rubin (1992). The first is
whether the clinic and parental care are conditionally independent given the
survival status, and the second is whether all three variables are mutually
independent. The MI datasets are generated from a Bayesian model in § 4.2
of Meng and Rubin (1992). Our aim here is to investigate the impact on the
r S, D
p ` and D
p ♦ by different parametrizations of tπijk u; and
test statistics D
S
S
the impact on the estimators rrL , rpS` and rpS♦ under different null hypotheses.
the `th imputed log-likelihood function is log f pX ` | πq “
ř Specifically,
`
`
`
c nc log πc , where X are the cell counts nc in the `th imputed dataset.
ř
`
pc “ n`c {n`` , where n`` “ c n`c .
Hence the unconstrained MLE of πc is π

27

MULTIPLE IMPUTATION LR TESTS

Parametrization (i)
δ = µ2 − µ1
1
2
3

Parametrization (ii)

4 0

1

2

3

Parametrization (iii)
4 0

1

2

3

4

Power/ Size

400

m = 30 , n = 100 m = 10 , n = 100 m = 3 , n = 100

m = 3 , n = 400 m = 3 , n = 1600

0

Nominal size α = 0.5 %

300
200
100
0
400

m
M

Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
Mm
M
Mm
Mm
mm
Mm
mM
M

m
M

m
M

m
M

MMMM MMMMMMMMMM
mM
mM
mM
M
mmmmmM
mmmmmmmmmmm
m
M

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
mM
mM
mM
mM
M
mM
mM

300
200
100
0
400

m
mM
M

mM
mm
Mm
Mm
mM
Mm
mM
mM
MM
Mm
mM
mM
mM
mM
mM
mM
mM
m
mM
M

m
mM
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M
mM
mM
mM
mM
mM
mM
mM
m

m
mM
M

Test statistics

mM
mM
mM
m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M

300
200
100
0
400

m
M
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
M
m
M
mM
mm
mM
Mm
mM
MM
mm
mM
mM
Mm
MM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m M
mM
mM
MM
mm
mM
mM
mM
mM
mM
m
M
M
M

300

m
M

200

mM
mM
mM
mM
m
mM
mM
mM
mM
mM
mM
MM
Mm
mm
mM
M

m
M

100
0
400

m
M

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m
M
m
mM
M

m
M

m
mM
mM
M

MMMMMMMM
mmmmmmmmm
mM
mM
mM
mM
mM
mM
mM
mM
mM
M

m
mM
mM
M

m
M

m
M
mM

M
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM

m
M

300
200
m
M

100
0

m
mM
M

M
mm
M

m
M

MMMMMM
MMMM
mmmmmmmmmmm
MMMM
Mmmm
mM
mm
mM
M

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
M

m
M
M
mm
M

m
M

Fig 6. The ratios of empirical power to empirical size under different parametrizations.
The nominal size is α “ 0.5%. In each plot, the vertical axis denotes the ratio, whereas
the horizontal axis denotes δ “ µ2 ´ µ1 . The legend in Figure 4 also applies here.

Table 6
Data from Meng and Rubin (1992). The notation “?” indicates missing label.
Survival Status (j)
Clinic (k)

Parental care (i)

Died

Survived

A

Less
More

3
4

176
293

B

Less
More

17
2

197
23

?

Less
More

10
5

150
90

Consequently, the joint log-likelihood based on the stacked data is
(4.1)

log f pX S | πq “

m ÿ
ÿ
`“1 c

where n`
c “

řm

`
`“1 nc .

n`c log πc “

ÿ

n`
c log πc ,

c

Thus the unconstrained MLE with respect to (4.1)

28

K.W. CHAN AND X.-L. MENG
Mean squared error of estimators of FMI
Parametrization (i)
δ = µ2 − µ1
1
2
3

m = 30 , n = 100 m = 10 , n = 100 m = 3 , n = 100

log(MSE)

m = 3 , n = 400 m = 3 , n = 1600

0

Parametrization (ii)
4 0

1

10

-5

m

m

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m M
M
m

1

2

M

m

M

M

m

m

M

m

m

m

m

M m
M

M

M

M

M

m

m

m

m

m
M m
M

M

M

M

M

m

m

m

m

m

M

M

M

M

M

m

m

m

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m M
M
m

M

M

M

M

M

M

M

m
M

m
M

M

m
M

4
m

m

m

m
M

3

m

m

5

-5

4 0

m

10

0

Parametrization (iii)

3

m

5
0

2

M

m

m

M

M

M

m

m

M

M

M

m

m
M

m

M

Estimators in

m

10
5
0
-5

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
m M
M

m

m

m

M

M

M

m

m

m

M

M

M

m

m

M

M

m

m
m

m

M

M

M

m

m

10

m

m

m

M

M

M

m

0
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

M

m
M

m

M

M

M

m

m

m
M M

M

M

m

5

-5

m

m

M

m

m

M

M

M

M

m

m

m

m

M

M

M

M

m

m
M

m
M

W1
W2
W3
W4
L1-MR92
L2-MR92(+)
L3-Proposal
L4-Proposal(+)
L5-Proposal(R)

m
m

m
M

M

m

m

m
m
M

m
M

M

M

m

m

M

M

m
M

m
M

M

10
5

-5

m

m
m

0
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

m
M

M

M

m
m
M

M

M

m

m
m
M

M

M

m

M
m
M

m
M

m
M

m
M

M

m
M

M

M

Fig 7. The MSEs of estimators of fm used in the test statistics. In each plot, the vertical
axis denotes the log of MSE, whereas the horizontal axis denotes the value of δ “ µ2 ´ µ1 .
The legend in Figure 4 also applies here.

ř `
`
`
is π
pcS “ n`
c {n` , where n` “
c nc . Similarly, we can find the constrained
MLEs under a given null. We consider the following parametrizations:
(i) ψijk “ πijk — the identity map,
(ii) ψijk “ eπijk — the exponential transformation, and
(iii) ψijk “ logtπijk {p1 ´ πijk qu — the logit transformation.
r L, D
p ` and D
p ♦ and the associated
p♦S of the tests D
The p-values prL , pp`
S and p
S
S
♦
`
estimates of rm , i.e., rrL , rpS and rpS , are shown in Table 7. A more detailed
comparison is deferred to Table 9 of the Appendix.
p ` and D
p ♦ are invariant to
The simulation outputs demonstrate that D
S
S
r L is not. Moreover, the impact on D
r L is large
parametrizations, whereas D
under the parametrization (iii). In particular, the value of rrL is inflated; and
r L are negative. We remark that rpS ě 0 always in our
some of the values of D
p` « D
p ♦ for testing conditional independence, a
simulations. In addition, D
S
S
hypothesis that is not rejected by either test. In contrast, for testing full
p ` and D
p ♦ are not very close to each other, but they both
independence, D
S
S
lead to essentially zero p-value, and hence both reject the null hypothesis.

29

MULTIPLE IMPUTATION LR TESTS

Table 7
r L, D
p ` and D
p ♦ under different parametrizations in § 4.2. The p-values
The LRTs using D
S
S
prL , pp`
p♦
S and p
S for testing full independence are essentially all zero, and hence omitted.
Parametrization (i)
H0 : Conditional independence
m

rrL , rpS` , rpS♦

3 0.54, 0.54, 0.38
10 0.50, 0.50, 0.70
50 0.31, 0.31, 0.45

r L, D
p `, D
p♦
D
S
S
0.08, 0.08, 0.09
0.14, 0.14, 0.12
0.11, 0.11, 0.10

prL , pp`
p♦
S ,p
S

H0 : Full independence
rrL , rpS` , rpS♦

r L, D
p `, D
p♦
D
S
S

0.93, 0.93, 0.92 0.31, 0.31, 0.38 54.2, 54.2, 51.4
0.87, 0.87, 0.88 0.56, 0.56, 0.70 45.4, 45.4, 41.7
0.90, 0.90, 0.91 0.33, 0.33, 0.45 51.5, 51.5, 47.3

Parametrization (ii)
H0 : Conditional independence
m

rrL , rpS` , rpS♦

3 0.54, 0.54, 0.38
10 0.50, 0.50, 0.70
50 0.30, 0.31, 0.45

p♦
r L, D
p `, D
D
S
S
0.08, 0.08, 0.09
0.14, 0.14, 0.12
0.11, 0.11, 0.10

prL , pp`
p♦
S ,p
S

H0 : Full independence
rrL , rpS` , rpS♦

p♦
r L, D
p `, D
D
S
S

0.92, 0.93, 0.92 0.30, 0.31, 0.38 54.4, 54.2, 51.4
0.87, 0.87, 0.88 0.54, 0.56, 0.70 46.1, 45.4, 41.7
0.90, 0.90, 0.91 0.32, 0.33, 0.45 52.1, 51.5, 47.3

Parametrization (iii)
H0 : Conditional independence
m

rrL , rpS` , rpS♦

p♦
r L, D
p `, D
D
S
S

p♦
prL , pp`
S
S ,p

H0 : Full independence
rrL , rpS` , rpS♦

p♦
r L, D
p `, D
D
S
S

3 1.08, 0.54, 0.38 ´0.07, 0.08, 0.09 1.00, 0.93, 0.92 0.61, 0.31, 0.38 43.9, 54.2, 51.4
10 0.99, 0.50, 0.70 ´0.10, 0.14, 0.12 1.00, 0.87, 0.88 1.09, 0.56, 0.70 33.7, 45.4, 41.7
50 0.63, 0.31, 0.45 ´0.10, 0.11, 0.10 1.00, 0.90, 0.91 0.65, 0.33, 0.45 41.3, 51.5, 47.3

These results reconfirm the conclusions in Meng and Rubin (1992). Last but
not least, the estimator rpS♦ does not change under different null hypotheses,
however it is not true for the estimators rrL and rpS` .
5. Conclusions, Limitations and Future Work. In addition to conducting a general comparative study of MI tests, we proposed two particup` “ D
p S pp
p♦ “ D
p S pp
rS` q and D
rS♦ q. Both
larly promising MI LRT based on D
S
S
test statistics are non-negative, invariant to re-parametrizations, powerful
to reject a false null hypothesis (at least for large enough m), and easy to
compute. Some additional features are discussed below.
p ` is very easy to compute, especially for users who have access only
• D
S
to standard complete-data computer subroutines for likelihood ratio
tests. One drawback is that rpS` “ maxp0, rpS q is an ad hoc way of ensuring non-negativity, and may be inconsistent in general. Fortunately,
the inconsistency does not appear to significantly affect the asymptotic
power, and the test performs rather well in our simulation studies.

30

K.W. CHAN AND X.-L. MENG

p ♦ is obtained in a principled way, with non-negativity and consis• D
S
tency of rpS♦ without ad hoc fixing. The resulting test has the desirable
monotonically increased power as H1 departs from H0 . However, the
consistency of rpS♦ requires a stronger assumption of EFMI for ψ, not
just for θ; and row independence of Xcom is needed for the ease of
computation. Nevertheless, when the rows are not independent, we
p L pp
rL♦ q.
can resort to the most general test statistic D
p ` and D
p ♦ are significant improvements over
In short, the tests based on D
S
S
their existing counterparts. Nevertheless, more study remains to be done,
and the following is only a partial list of the open problems.
• One grand challenge is that when the missing data mechanism is not
ignorable and the imputers fail to fully take that into account, there
will be an issue of uncongeniality (Meng, 1994a). The most recent work
by Xie and Meng (2017) provides theoretical tools for addressing such
an issue in the context of estimation, and research is needed to extend
their findings to the setting of hypothesis testing.
• All results in this article assume EFMI. It is a strong assumption
that is unlikely to hold in practice. Although the violation of this
assumption may not invalidate a test, it will affect its power. It is
therefore desirable to explore MI tests without this assumption.
p ♦ relies on a stronger assumption of EFMI on ψ. We can
• The robust D
S
modify it so only EFMI on θ is required, but the modification may
be very difficult to compute and may require users to have access to
non-trivial complete-data procedures. Hence a computational feasible
robust test that only assumes EFMI on θ needs to be developed.
• Because the FMI is a fundamental nuisance parameter here and there
is no (known) pivotal quantity, all MI tests are approximate in nature.
In particular, they all have the potential of doing poorly when FMI is
large and/or m is small. It is therefore of both theoretical and practical
interest to seek powerful MI tests that are least affected by FMI.
References.
Barnard, J. and Rubin, D. B. (1999) Small-sample degrees of freedom with multiple imputation. Biometrika, 86, 948–955.
Bayarri, M. J., Benjamin, D. J., Berger, J. O. and Sellke, T. M. (2016) Rejection odds
and rejection ratios: A proposal for statistical practice in testing hypotheses. Journal
of Mathematical Psychology, 72, 90–103.
Benjamin, D. J., Berger, J. O., Johannesson, M., Nosek, B. A., Wagenmakers, E.-J., Berk,
R., Bollen, K. A., Brembs, B., Brown, L., Camerer, C. et al. (2017) Redefine statistical
significance. Nature Human Behaviour, 1–5.
Blocker, A. W. and Meng, X.-L. (2013) The potential and perils of preprocessing: Building
new foundations. Bernoulli, 19, 1176–1211.

MULTIPLE IMPUTATION LR TESTS

31

Cox, D. R. and Reid, N. (1987) Parameter orthogonality and approximate conditional
inference. Journal of the Royal Statistical Society B, 49, 1–39.
Li, K. H., Meng, X.-L., Raghunathan, T. E. and Rubin, D. B. (1991a) Significance levels
from repeated p-values with multiply-imputed data. Statistica Sinica, 1, 65–92.
Li, K. H., Raghunathan, T. E. and Rubin, D. B. (1991b) Large-sample significance levels
from multiply imputed data using moment-based statistics and an F reference distribution. Journal of the American Statistical Association, 86, 1065–1073.
Meng, X.-L. (1994a) Multiple-imputation inferences with uncongenial sources of input.
Statistical Science, 9, 538–573.
Meng, X.-L. (1994b) Posterior predictive p-values. The Annals of Statistics, 22, 1142–1160.
Meng, X.-L. (2002) Discussion of “Bayesian measures of model complexity and fit” by
Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and Van Der Linde, A. Journal of the
Royal Statistical Society B, 64, 633.
Meng, X.-L. and Rubin, D. B. (1992) Performing likelihood ratio tests with multiplyimputed data sets. Biometrika, 79, 103–111.
Rubin, D. B. (1976) Inference and missing data. Biometrika, 63, 581–592.
Rubin, D. B. (1978) Multiple imputations in sample surveys - a phenomenological Bayesian
approach to nonresponse. Proceedings of the Survey Research Methods Section of the
American Statistical Association, 20–34.
Rubin, D. B. (2004) Multiple Imputation for Nonresponse in Surveys. John Wiley & Sons.
Rubin, D. B. and Schenker, N. (1986) Multiple imputation for interval estimation from
simple random samples with ignorable nonresponse. Journal of the American Statistical
Association, 81, 366–374.
Serfling, R. J. (2001) Approximation Theorems of Mathematical Statistics. WileyInterscience.
Tu, X. M., Meng, X.-L. and Pagano, M. (1993) The AIDS epidemic: Estimating survival
after AIDS diagnosis from surveillance data. Journal of the American Statistical Association, 88, 26–36.
van der Vaart, A. W. (2000) Asymptotic Statistics. Cambridge University Press.
Wallace, D. L. (1980) The Behrens-Fisher and Fieller-Creasy problems. In R. A. Fisher:
An Appreciation (eds. S. E. Fienberg and D. V. Hinkley), 119–147. Springer New York.
Wang, N. and Robins, J. M. (1998) Large-sample theory for parametric multiple imputation procedures. Biometrika, 85, 935–948.
Xie, X. and Meng, X.-L. (2017) Dissecting multiple imputation from a multi-phase inference perspective: What happens when God’s, imputer’s and analyst’s models are
uncongenial? (with discussion). Statistica Sinica, 27, 1485–1594.

32

K.W. CHAN AND X.-L. MENG

APPENDIX A: SUPPLEMENTARY RESULTS
♦
A.1. Another Motivation for rpL
. The definition of rpL♦ is also motivated by the following observation. First, observe that one simple method to
construct an always non-negative estimator of rm is to perturb ψp0˚ and ψp0`
by a suitable amount, say ∆, so that the perturbed version of rpL is always
non-negative, and is still asymptotically equivalent to the original rpL . We
show, in Theorem A.1 below, that the right amount of ∆ is ∆ “ ψp˚ ´ ψp0˚ .
Using the perturbed version of rpL , we obtain

rpL4 “

m ` 1 p4
δ ,
kpm ´ 1q L

where
δpL4

m
2 ÿ
“
log
m `“1

#

f pX ` | ψp` q f pX ` | ψp0˚ ` ∆q
f pX ` | ψp˚ q f pX ` | ψp` ` ∆q

+

0

m
1 ÿ
“
dL pψp0` ` ∆, ψp` | X ` q.
m `“1

Then we have the following result.
Theorem A.1. Suppose RCθ . Under H0 , we have (i) rpL4 ě 0 for all
m, n; and (ii) rpL4 l rpL as n Ñ 8 for each m.
Although rpL4 ě 0, it is only invariant to affine transformations, and not
robust against θ0 , and less computational feasible than rpL ; see § 3. However,
it gives us some insights on how to construct a potentially better estimator.
Note that, in (A.1), the constrained MLE is not used in dL p¨, ¨ | X ` q, but it is
still always non-negative. We call this a “pseudo” LRT statistics. Then, δpL4
is just a multiple of an average of many “pseudo” LRT statistics. In order
to find a good estimator of rm , we may seek for an estimator which admits
this form. Indeed, our proposed estimator rpL♦ also takes the same form:
rpL♦ “

m
m`1 1 ÿ
dL pψp˚ , ψp` | X ` q.
hpm ´ 1q m `“1

A.2. Results for Dependent Data. This is a supplement for § 3.1.
If the data are not independent, dpL l dpS is still true under the following
conditions.
S

Assumption 5.
(a) Define Rpψq “ L pψq ´ Lpψq, where Lpψq “
řm
S
`
´1 log f pX S | ψq. For
´1
pmnq
`“1 log f pX | ψq and L pψq “ pmnq
each m, as n Ñ 8,
ˇ
ˇ
ˇ B
ˇ
sup |Rpψq| “ Op p1{nq,
sup ˇˇ Rpψqˇˇ “ Op p1{nq.
Bψ
ψPΨ

ψPΨ

MULTIPLE IMPUTATION LR TESTS

33

(b) For each m, there exists a continuous function ψ ÞÑ Lpψq, which is
free of n but may depend on m, such that, as n Ñ 8,
ˇ
ˇ
sup ˇLpψq ´ Lpψqˇ “ op p1q.
ψPΨ

(c) Let ψ0˚ “ arg maxψPΨ : ψpθq“θ0 Lpψq and ψ ˚ “ arg maxψPΨ Lpψq. For
any fixed m, and for all ε ą 0, there exists δ ą 0 such that
(
(
sup
Lpψ0˚ q ´ Lpψq ě δ,
sup
Lpψ ˚ q ´ Lpψq ě δ.
ψPΨ : |ψ0˚ ´ψ|ąε
θpψq“θ0

ψPΨ : |ψ ˚ ´ψ|ąε

Conditions (b) and (c) in Assumption 5 are standard RCs that are usually assumed for M-estimators (see § 5 of van der Vaart (2000)); whereas
condition (a) is satisfied by many models (see Example A.1 below).
Theorem A.2. Suppose RCθ and Assumption 5. Under both H0 and
H1 , we have (i) dpS , rpS ě 0 for all m, n; (ii) dpS , rpS are invariant to the
parametrization of ψ for all m, n; and (iii) dpL l dpS and rpL l rpS as n Ñ 8
for each m.
p S and D
p ` approxiTheorem A.2 implies that the handy test statistics D
S
p L and D
p ` for dependent data, provided that Assumption 5 holds.
mate D
L
Example A.1. Consider a stationary autoregressive model of order one.
Suppose the complete data X “ pX1 , . . . , Xn q| is generated as following:
X1 „ Np0, v 2 q and rXi |Xi´1 s „ NpφXi´1 , σ 2 q for i ě 2, where v 2 “
σ 2 p1 ` φq{p1 ´ φq. Then ψ “ pφ, σ 2 q| , and
m
1
1
1 ÿ X1`
n´1
Lpψq “ ´ logp2πq ´
log v 2 ´
´
log σ 2
2
2
2n
mn `“1 2v
2n

´

m n
` q2
1 ÿ ÿ pXi` ´ φXi´1
,
mn `“1 i“2
2σ 2

1
1
pX11 q2
mn ´ 1
S
L pψq “ ´ logp2πq ´
log v 2 ´
log σ 2
´
2
2
2mn
2mnv
2mn
m n
m
` q2
1 ÿ ÿ pXi` ´ φXi´1
1 ÿ pX1` ´ φXn`´1 q2
´
´
.
mn `“1 i“2
2σ 2
mn `“2
2σ 2
Then, it is easy to see that condition (a) of Assumption 5 is satisfied.

34

K.W. CHAN AND X.-L. MENG

APPENDIX B: PROOFS
Proof of Theorem 2.1. (i, ii) From (2.3), we know dpL ě 0 is invariant
to parametrization ψ. (iii) Since dpL is invariant to transformation of ψ, we
assume, without loss of generality, that ψ admits a parameterization such
that Covpθp` , ηp` q l 0 by taking suitable linear transformation of ψ. Also write
η q based on X ` ; and recall that Uθ` “ U `
Uη` as an efficient estimator of Varpp
p based on X ` .
is an efficient estimator of Varpθq
ř
`
Using Taylor’s expansion on ψ ÞÑ Lpψq “ m´1 m
`“1 log f pX | ψq around
ψp˚ “ ppθp˚ q| , pp
η ˚ q| q| , we know that for ψ l ψp˚ ,
¯|
¯
´
1´
(B.1)
ψ ´ ψp˚ Ipψp˚ q ψ ´ ψp˚ ,
Lpψq l Lpψp˚ q ´
2
where Ipψq “ ´B 2 Lpψq{BψBψ | , which satisfies
˜ ´1
¸
0
U
θ
(B.2)
Ipψp˚ q l
´1
0
Uη
ř
`
p˚
p˚
with U η “ m´1 m
i“1 Uη . Under the null, ψ l ψ0 . So, using (B.1), we have
´
¯|
´
¯
dpL l
ψp0˚ ´ ψp˚ Ipψp˚ q ψp0˚ ´ ψp˚ ,
¸| ˜ ´1
¸
˜
¸˜
0
Uθ
θ0 ´ θp˚
θ0 ´ θp˚
l
´1
ηppθ0 q ´ ηppθp˚ q
ηppθ0 q ´ ηppθp˚ q
0
Uη
(B.3)

|

´1

|

l pθ ´ θ0 qU θ pθ ´ θ0 q “ dr1W ,

where we have used (a) θp˚ l θ; see, e.g., Lemma 1 of Wang and Robins
?
(1998), and (b) ηppθ0 q ´ ηppθp˚ q “ Op p1{nq if θ0 ´ θp˚ “ Op p1{ nq; see Cox and
Reid (1987). Since dr1W l drL (Meng and Rubin, 1992), we have dpL l drL .
Proof of Proposition 2.2. The given condition implies that
ψp` “ ppθp` q| , pp
η ` q| q| ,
ψp˚ “ ppθp˚ q| , pp
η ˚ q| q| ,

ψp0` “ pθ0| , pp
η ` q| q| ,
ψp˚ “ pθ| , pp
η ˚ q| q| .
0

0

Clearly, we also have the decomposition: L` pψq “ L`: pθq ` L`; pηq for all `,
where L`: pθq “ L: pθ | X ` q and L`; pηq “ L; pη | X ` q. Then,
dL ´ dpL “
“

m
)
2 ÿ ! ` p`
L pψ q ´ L` pψp0` q ´ L` pψp˚ q ` L` pψp0˚ q
m `“1
m
)
2 ÿ ! ` p`
L: pθ q ´ L`: pθp˚ q ě 0
m `“1

MULTIPLE IMPUTATION LR TESTS

35

since L`: pθp` q ě L`: pθp˚ q for all `.
Proof ofp Corollary 2.3. Applying Taylor’s expansion on ψ ÞÑ L` pψq,
we can find ψ ` lying on the line segment joining ψp` and ψp0` such that
p
1 ´ p` p` ¯| ` ` ´ p` p` ¯
` p`
` p`
I pψ q ψ0 ´ ψ ,
ψ0 ´ ψ
L pψ0 q “ L pψ q ´
2
p
`
2
`
|
`
p ´B L pψq{BψBψ
p . By thep lower order variability of I pψ ` q,
where I pψq “
we can find ψ ˚ such that I ` pψ ` q l I ` pψ ˚ q for all `. Then, using similar
techniques as in (B.2) and (B.3), we have
p
1 ´ p` p` ¯| ` ˚ ´ p` p` ¯
` p`
` p`
L pψ q ´ L pψ0 q l
I pψ q ψ0 ´ ψ
ψ0 ´ ψ
2
p
¯
¯| ´1 ´
1´
(B.4)
l
θ0 ´ θp`
θ0 ´ θp` U
2
p
for some matrix U . Similarly, we have
¯| p ´1 ´
¯
1´
˚
` p˚
` p˚
p
θ0 ´ θ
(B.5)
L pψ q ´ L pψ0 q l
U
θ0 ´ θp˚ .
2
Write Ab2 “ AA| for any appropriate matrix A. Using (B.4), (B.5) and the
cyclic property of trace, we have
m "
¯| p ´1 ´
¯ ´
¯| p ´1 ´
¯*
1 ÿ ´
`
`
˚
˚
p
p
p
p
p
dL ´ dL l
θ0 ´ θ
U
θ0 ´ θ ´ θ0 ´ θ
U
θ0 ´ θ
m `“1
+ff
«p #
m
¯b2
¯b2 ´
´1
1 ÿ´
˚
`
p
p
“ tr U
´ θ0 ´ θ
θ0 ´ θ
m `“1
«p
ff
ˆp
˙
m
)
´
¯
´1 1 ÿ !
´1
b2
´1
`
b2
l tr U
pθp q ´ θ
l tr U B l tr Uθ,0
Bθ
m `“1
p is a deterministic matrix that depends on both θ0
as m, n Ñ 8, where Uθ,0
pr
´1
and θ‹ , and satisfies npU ´ Uθ,0 q Ñ 0. Note that trpUθ,0
Bθ q “ kr0 , for some
pr

´1
Bθ q{k, proving (ii). (But
finite r0 by Assumption 2. Then rpL Ñ r0 “ trpUθ,0
Uθ,0 may not equal to Uθ , and hence rpLpmay not be consistent for rm .)
pr
pr
If H0 is true, then θ Ñ θ0 “ θ‹ and U l U l Uθ “ Uθ,0 . Then, rpL Ñ r
as m, n Ñ 8. So, (i) follows.

Proof of Theorem 2.4. (i, ii) It is trivial by the definition of prpL♦ . (iii)
Applying Taylor’s expansion to ψ ÞÑ L` pψq again, we know there is ψ ` lying

36

K.W. CHAN AND X.-L. MENG

on the line segment joining ψp` and ψp˚ such that
p
1 ´ p˚ p` ¯| ` ` ´ p˚ p` ¯
` p˚
` p`
I pψ q ψ ´ ψ .
ψ ´ψ
(B.6)
L pψ q “ L pψ q ´
2
p
p
By the lower order variability
of I ` pψ ` q, we know that I ` pψ ` q l Ipψp˚ q for
ř
m
all `, where Ipψq “ m´1 `“1 I ` pψq. We also know that ψp˚ l ψ. Thus
m
1 ÿ ´ p˚ p` ¯| p˚ ´ p˚ p` ¯
Ipψ q ψ ´ ψ
ψ ´ψ
m `“1
#
+
m ´
¯b2
ÿ
1
“ tr Ipψp˚ q
ψp˚ ´ ψp`
m `“1
#
+
m ´
¯b2
´
¯
ÿ
1
l tr Ipψp˚ q
ψp` ´ ψ
l tr Uψ´1 Bψ
m `“1

δ L ´ δpL l

(B.7)

pr

as m, n Ñ 8. By the assumption of EFMI of ψ, we have rpL♦ Ñ r.
Proof of Lemma 2.5. First, recall that, as n Ñ 8, the observed data
D
MLE θpobs of θ satisfies (2.4), which can be written as rθpobs | θ‹ s « Nk pθ‹ , Tθ q,
D

where A1,n « A2,n means that A1,n and A2,n have the same asymptotic
distribution, i.e., there exist deterministic sequences µn and Σn such that
´1{2
´1{2
ñ A for some non-degenerate
ñ A and pA2,n ´ µn qΣn
pA1,n ´ µn qΣn
random variable A. From Assumption 3, a proper imputation model is used.
So, we have (2.5), which is equivalent to say that, as n Ñ 8,
”
ıD
(B.8)
θp` | Xobs « Nk pθpobs , Bθ q,
independently for for ` “ 1, . . . , m. Therefore we can represent
(B.9)

D
1{2
θpobs « θ‹ ` Tθ W,
D
1{2
θp` « θpobs ` Bθ Z` ,

(B.10)

` “ 1, . . . , m

iid

where Z1 , . . . , Zm , W „ Nk p0, Ik q. Also write Z` “ pZ1` , . . . , Zk` q| , for ` “
D

1, 2, . . . , m, and W “ pW1 , . . . , Wk q| . Averaging (B.10) over `, we have θ «
ř
1{2
θpobs ` Bθ Z ‚ , where Z ‚ “ m´1 m
`“1 Z` . Since Bθ “ rUθ , we have
´1{2

Uθ

´1{2

Uθ

D

pθp` ´ θ‹ q « p1 ` rq1{2 W ` r 1{2 Z` ,
D

pθ ´ θ‹ q « p1 ` rq1{2 W ` r 1{2 Z ‚ .

37

MULTIPLE IMPUTATION LR TESTS

Note that (2.6) implies Uθ l U . Under H0 , we have θ‹ “ θ0 and
D

1

dL l dW «

k !
)2
ÿ
p1 ` rq1{2 Wi ` r 1{2 Zi` ,
i“1

D
dpL l drL l dr1W «

k !
)2
ÿ
p1 ` rq1{2 Wi ` r 1{2 Z i .
i“1

After some simple algebra, we obtain
D
rpL` «

k
pm ` 1qr ÿ 2
sZi
mk
i“1

D
p` «
D
L

and

m

řk

p1 ` rq1{2 Wi ` r 1{2 Z i‚
ř
mk ` pm ` 1qr ki“1 s2Zi

(2

i“1

,

ř
2
m
where s2Zi “ pm ´ 1q´1 m
`“1 pZi` ´ Z i‚ q is the sample variance of tZi` u`“1 .
Since Wi , Z i‚ and s2Zi are mutually independent for each fixed i, we can
p ` to
simplify the representation of D
L
D
rpL` «

k
pm ` 1qr ÿ 2
H
mpm ´ 1qk i“1 i
iid

and

řk
2
` D pm ´ 1qtm ` pm ` 1qru i“1 Gi
p
DL «
,
řk
mpm ´ 1qk ` pm ` 1qr i“1 Hi2

iid

where G2i „ χ21 and Hi2 „ χ2m´1 , for i “ 1, . . . , k, are all mutually independent. Clearly, they can be further simplified to (2.14).
Proof of Theorem 2.6. Similar to (B.9) and (B.10), we can have a
more general representation:
D
1{2
ψpobs « ψ ‹ ` Tψ W ;

D
1{2
ψp` « ψpobs ` Bψ Z` ,

` “ 1, . . . , m,

iid

where Z1 , . . . , Zh , W „ Nh p0, Ih q. Also write Z` “ pZ1` , . . . , Zh` q| , for ` “
1, 2, . . . , m, and W “ pW1 , . . . , Wh q| . Using (B.7), we have
#
+
m ´
¯´
¯|
ÿ
1
ψp` ´ ψ ψp` ´ ψ
δ L ´ δpL l tr Ipψp˚ q
m `“1
«
ff
m !
)b2
ÿ
`
˘
1
D
pTψ ´ Uψ q1{2 Z` ´ Z ‚
« tr Uψ´1
m `“1
“

m
m h
`
˘b2 )
r ÿÿ
1 ÿ !
tr rIh Z` ´ Z ‚
“
pZi` ´ Z i‚ q2 .
m `“1
m `“1 i“1

38

K.W. CHAN AND X.-L. MENG

Equivalently, we can say δ L ´ δpL ñ rχ2hpm´1q {m as n Ñ 8. Hence
rpL♦ ñ r ¨

m`1
¨ χ2
,
hmpm ´ 1q hpm´1q

which is equivalent to (2.15). Note that it is true under both H0 and H1 .
Proof of Theorem 2.7. From the representations of dp♦L and rpL♦ in Lemma
2.5 and Theorem 2.6, we know that they are asymptotically (n Ñ 8) independent. The proof then follows the derivation for Lemma 2.5.
Proof of Proposition 3.1. It is trivial.
Proof of Theorem A.1. (i) Using the representation (A.1), we can
easily see that rpL4 ě 0. (ii) It suffices to show
´1

m

m
ÿ

dL pψp0` ` ∆m , ψp` | X ` q l dL ´ drL ,

`“1

where ∆m “ ψp˚ ´ ψp0˚ . Under H0 , ∆m l 0 and ψp0` l ψp` , so ψp0` ` ∆m l ψp` .
Using Taylor’s expansion on ψ ÞÑ L` pψq around its maximizer ψp` , we have
for ψ l ψp` that
¯|
´
¯
1´
L` pψq l L` pψp` q ´
ψ ´ ψp` I ` pψp` q ψ ´ ψp` .
2
Under the parametrization of ψ in the proof of Theorem 2.1, we know that
` ˘´1
the upper k ˆ k sub-matrix of I ` pψp` q is U `
. Using the lower order vari` ` ˘´1
´1
`
ability of U , we have U
l U and
m
1 ÿ
dL pψp0` ` ∆m , ψp` | X ` q
m `“1
m
¯|
´
¯
1 ÿ ´ p`
l
ψ0 ` ∆m ´ ψp` I ` pψp` q ψp0` ` ∆m ´ ψp`
m `“1

l

m
1 ÿ p`
1
´1
pθ ´ θq| U pθp` ´ θq “ dW ´ dr1W l dL ´ dpL .
m `“1

Therefore, the desired result follows.
Proof of Theorem A.2. Throughout this proof, conditions (a), (b)
and (c) refer to the list given in Assumption 5. (i, ii) It trivially follows

MULTIPLE IMPUTATION LR TESTS

39

from the definitions of dpS and rpS . (iii) First, by the definition of maximizer
and condition (a), we have
S
S
Lpψp˚ q ´ LpψpS q “ Lpψp˚ q ´ L pψpS q ` L pψpS q ´ LpψpS q
S
S
ď Lpψp˚ q ´ L pψp˚ q ` L pψpS q ´ LpψpS q
ˇ
ˇ
S
ˇ
ˇ
ď 2 sup ˇLpψq ´ L pψqˇ “ Op p1{nq,
ψPΨ

which, together with condition (b), imply that
Lpψ ˚ q ´ LpψpS q

(B.11)

) !
)
( !
“ Lpψ ˚ q ´ Lpψ ˚ q ` Lpψ ˚ q ´ LpψpS q ` LpψpS q ´ LpψpS q
)
ˇ !
ˇ
ď 2 sup ˇLpψq ´ Lpψqˇ ` Lpψp˚ q ´ LpψpS q “ op p1q.
ψPΨ

Using (B.11)
pr
ψp˚ Ñ ψ ˚ . So,
(B.12)

pS pr ˚
ˇand (c), ˇwe have ψ Ñ ψ . By (b) and (c), we also have
ˇ pS p˚ ˇ pr
ˇψ ´ ψ ˇ Ñ 0 as n Ñ 8. By the definition of maximizer,
S
0 “ ∇L pψpS q “ ∇LpψpS q ` ∇RpψpS q,

where ∇gpψq “ Bgpψq{Bψ is the gradient of ψ ÞÑ gpψq. By condition (a), we
know ∇RpψpS q “ Op p1{nq. Thus, together with (B.12), we have ∇LpψpS q “
p we have ∇Lpψp˚ q “ 0.
Op p1{nq. Also, by the definition of MLE,
By Taylor’s expansion, there exits ψ such that
p ) ´
!
¯
|
˚
S
(B.13)
Lpψp q ´ Lpψp q “
∇Lpψq
ψp˚ ´ ψpS “ op p1{nq,
p
where we have used the continuity of ψ ÞÑ ∇Lpψq to yield ∇Lpψq “ Op p1{nq.
Rewriting (B.13), we have
(B.14)

S
Lpψp˚ q ´ L pψpS q “ RpψpS q ` op p1{nq.

Similar to (B.14), we have
(B.15)

S
Lpψp0˚ q ´ L pψp0S q “ Rpψp0S q ` op p1{nq.

Then, using (B.14) and (B.15), we have
ˇ
ˇ
ˇ!
) !
)ˇ
S
S
ˇp
ˇ
ˇ
ˇ
ˇdL ´ dpS ˇ “ 2n ˇ Lpψp˚ q ´ L pψpS q ´ Lpψp0˚ q ´ L pψp0S q ˇ
ˇ
ˇ
ˇ
ˇ
“ 2n ˇRpψpS q ´ Rpψp0S q ` op p1{nqˇ .
Now consider two cases.

40

K.W. CHAN AND X.-L. MENG

(i) Under H0 , we have dpL “ Op p1q and ψp0S l ψpS .
implies RpψpS q ´ RpψpS q “ op p1{nq. Then, we have
0

(a)
ˇThus condition
ˇ
ˇp
ˇ
p
p
ˇdL ´ dS ˇ “ op pdL q.

pr
(ii) Under H1 , we have dpL Ñ 8. Condition (a) and (B.11) imply that
S
S pS
p˚
Lpψp˚ q ´ L pψpS q ˇ“ Op p1{nq.
ˇ Similarly, we also haveˇ Lpψ0 q ˇ´ L pψ0 q “
ˇ
ˇ
ˇ
ˇ
Op p1{nq. Hence ˇdpL ´ dpS ˇ “ Op p1q. Thus we have ˇdpL ´ dpS ˇ “ op pdpL q.

ˇ
ˇ
ˇ
ˇ
Therefore, under either H0 or H1 , we also have ˇdpL ´ dpS ˇ “ op pdpL q. Since
dpL l dpS and dL “ dS , we know rpL l rpS .
Note that, even under the assumption of this theorem, rpS and rpS♦ are not
equivalent. From (3.5) and (3.6), rpS and rpS♦ are a “difference of difference”
estimator and a “difference” estimator, respectively. So, the “bias” of using
S
L pψq cannot be canceled out in rpS♦ .
APPENDIX C: ADDITIONAL FIGURES AND TABLE
This section presents additional figures and tables in § 2.5 and § 4
• Figure 8: the performance of different approximations to the reference
null distribution when α “ 5%; see § 2.5.
• Figures 9 and 10: the empirical distributions of the p-values under H0
and parametrizations (i) and (iii), respectively; see § 4.1.
• Figures 11 and 12: the empirical power functions and the empirical
ratio of power-to-size for size 5% tests, respectively; see § 4.1.
• Table 8: the ranges of empirical sizes over different parametrizations
for size 5% tests; see § 4.1.
• Table 9: detailed results of the care-survival example in § 4.2.
Department of Statistics
Harvard University
Cambridge, Massachusetts 02138, U.S.A.
E-mail: kinwaichan@g.harvard.edu; meng@stat.harvard.edu

41

MULTIPLE IMPUTATION LR TESTS

Nominal size α = 5 %
k=1

k=8

k=4

k=2

FMI
10% 30% 50% 70% 90% 10% 30% 50% 70% 90% 10% 30% 50% 70% 90% 10% 30% 50% 70% 90%

m=3

Empirical size

5%
4.5%
4%
3.5%
3%
2.5%
5%

m=5

4.5%
4%

2
3
2 3
3
1 1
2
3 3
2 2
1
3 3
2
3 3
2
3
2
1
2
2 2 2
2
1
3
3
2
1
3
2 3
3 3 3
1 2
2
1 1
3 2
3 2
3 3
1 2
2
3
1 3
2
2 3
1
2
2 3
1
3 3
1 3
2
1 3
1
2
2
2 2

3.5%

3
2
3

1
2

3
1 3
2
1
2
3 3
3 3 3
1 2
3
1
2 2
2 2
2 3
1 1
3 1
3
3
2 1
1 3
1
2
1
2
2 2
1
1

3 3 3
2
2 2 3
3 3 2
1
1 1
2
2

3
2
1
3
1 1
2
3 2
3 2
2 2
3
1 3
2 3
3
2 1
1
2

3
2
3
1
2

3
2
3
1
2

1

3 3
3 3
2
3 3
1 3
2
3 3
1
1 3
2
2
2
3
1
3
2 3
2
3
2
2 3
2
3
1 3
3
3 3
2
1
1
2
2 2
1
1 3
1 2
1
2 2
2 2
2 2
2
3 3 1
1
1
1
1 2
1
2
1
1 1
1
1 1
1

1
1

1 1

1

1
1 1
2
3 3 3
2
2
2
2
2
3 3
1 3
1
2
2
2
3 3
2
2
1 3
2 3
3
3
1 3
2
1 3
2
1 3
2
1 3
2
2
1
2
2 3
1
1 3
2
2
1
1 3
2
3
3
1
2
1
3
2
1 3
1 3
1 3
1
1 3
2
1 3
2
2 2 2
2 3
2
1 3
2 3
2 3
2 3
2
2 3
1 3
1 3
3 3
1 2 2
3
1
1 3
1
2
3
2 3
1 2
2 3
3 3
1 1 1 1 1
2
1 3
2
2
2
1 1
1
1
3
1 1 1
1
1
3
2 2 2
2
2 2
1
1 1
1
1
1
1 2
1
1
1
2
1
1
1
1 1 1

Approximations
LRR (91)
New proposal

3%
2.5%
5%

m=7

4.5%

2 2
2 2
2 3
2 3
2
2
1 3
3
3
1
2
2
2 3
2 2
2
1 3
2
2
1
3
2
1 3
3
2
3
1 1
1 3
2
1 3
2
1 3
2
1 3
2
3
3 3
2
1
1 3
2
2
1 3
2
1
1 3
3
3
1 2
2
2 2
1
3
1
2
3
2 3
2
1 1
2 1
1 2
1 3
3
2
1 3
1
2
2
2
3 3
2
2
2
1
3 3
1 3
1 3
1 2
1 3
2
2
3
1 2
1 2
1 3
1 3
1 3
1
1 3
1 3
1 3
2 1
3
1 1 2
2
3 3
1 3
1 3
1
1 3
2 3
1 1 1 1 1
1 2 2 2 2 2
1
2 1
2 3
2 1
2 1
1
1
1 1
1 1
1
1
1

4%

1

1

3.5%

1

1 1

1

3%
2.5%

Fig 8. The performances of two different approximate null distributions when the nomp or α
r, and the horizontal axis denotes
inal size is α “ 5%. The vertical axis denotes α
the value of fm . The number attached to each line denotes the value of τ “ h{k. The
proposed approximation α
p is denoted by thick solid lines with triangles, and the existing
approximation α
r is denoted by thin dashed lines with circles.

42

K.W. CHAN AND X.-L. MENG
Parametrization (i)
Case 1 ( Small ρ , p , f )
Nominal size α
1% 3% 5% 7% 9%

Empirical size

1% 3% 5% 7% 9%

m
M

m
M

m
M

7%

m
M

5%

m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M

Experiment II ( Vary p )

9%

m
M

m
M
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M

m
M
m
M

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

m
M

m
M
m
M

1%

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M

3%

m
M

m
M

m
M

m
M

5%

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M
m
M

m
M

m
M

m
M

m
M

m
M

m
M

Test statistics

m
M
m
M

m
M

m
M

m
M

m
M
m
M

7%

m
M

m
M

m
M

m
M

m
M

m
M
m
M
m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

1%

m
M
m
M

m
M

m
M

3%

m
M
m
M
m
M

m
M

m
M

m
M

5%

m
M
m
M

m
M

m
M

m
M
m
M

m
M

m
M

m
M

7%

m
M

m
M

m
M
m
M

m
M
m
M

m
M
m
M
m
M

m
M

m
M

m
M
m
M

m
M

m
M
m
M

m
M

1% 3% 5% 7% 9%

m
M

m
M
m
M

m
M

m
M

m
M

9%

1% 3% 5% 7% 9%

m
M

m
M

m
M

m
M

m
M

m
M

1%

m
M

m
M

m
M

3%

1% 3% 5% 7% 9%

m
M

Case 5 ( Large ρ , p , f )

Case 4

m
M

m
M

9%

Experiment III ( Vary f )

Experiment I ( Vary ρ )

m
M

Case 3

Case 2

m
M

m
M

Fig 9. The comparison between empirical size and nominal size α under parametrization
(i) for α P p0, 10%s. The Wald tests (W-1,2,3,4) and LRTs (L-0,1,2,3,4,5) are represented
r L (L-1: MR92) (Meng
by grey dashed and black solid lines, respectively. The LRT statistic D
r ` (L-2: MR92(+)) are the existing counterparts
and Rubin, 1992) and its modification D
L
p ♦ (L-5: Proposal(R)).
p S (L-3: Proposal), D
p ` (L-4: Proposal(+)) and D
of our proposals D
S
S

Parametrization (iii)
Case 1 ( Small ρ , p , f )

Case 3

Case 2
M

m
Nominal size α
M
m
1% 3% 5% 7% M 9%

M
m
1% 3% 5% 7% 9%

m

M
m

9%

M
m

7%
5%

M
m

3%

M
m

M
m

M
m

M
m

M
m

M
m

M
m

m
M

m
M

m
M

m
M

m
M

m
M

m
M

m
M
m
M

7%

m
M

m
M

m
M

5%

m
M
m
M

m
M

m
M
M
m

m
M

m
M

m
M

M
m

m
M

m
M

1%

m
M
M

m
M

m
M
m
M

M
m

M
m

M
m

M
m

M
m
m
M

m
M

M

M

m
m

M

M

M

M

M

M

Test statistics

m
M

M

m

m

m
M

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m

m
m

M
m
m

m

m

M
m

m
M

M

M

M
M

M
M m
m

M

M

M

M
m
M

m
M

m
M

m
M
m
M
m
M

m
M
m
M

M

m
M

m

m
M

m
M

m
M

M
m

m
M

m
M

m
M

M
m

m
M

m
M

m
M

m
M

M
m

M
M m
m

M
m

m
M

m
M

m
M

M
m

M
M m
m

m
M

1%

3%

m
M
m
M

m
M

m
M

m
M

M

m
1% 3% 5% 7% 9%
M

M
m

m
M

m
M

m
M

m
M

M
m

M
M m
m

M
m

m
M
m
M

m
M

m
M
m
M

M
m

M
m

M
m

m
M

m
M

m
M

M
m

m
M

m
M

M
m

M
M m
m

M
m

M
m
M
m
M
m

M
m

m
M

1% 3% 5% 7% 9%
M
m

M
m

M
m

M
m

1%

9%

M
m

M
m

M
m
M
m

3%

Experiment II ( Vary p )

Empirical size

M
m

M
m

5%

Experiment III ( Vary f )

Experiment I ( Vary ρ )

M
m
M
m

7%

M
m

M
m

M
m

Case 5 ( Large ρ , p , f )

Case 4

1% 3% 5% 7% 9%
M
m

M
m

M
m

9%

M
m

m

m

m

m

m

m

m

m

m
m
M

m
M

Fig 10. The comparison between empirical size and nominal size α under parametrization
(iii) for α P p0, 10%s. The legend in Figure 9 also applies here.

43

MULTIPLE IMPUTATION LR TESTS

Parametrization (i)
δ = µ2 − µ1
1
2
3

m = 30 , n = 100 m = 10 , n = 100 m = 3 , n = 100

Power

m = 3 , n = 400 m = 3 , n = 1600

0
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%
100%
80%
60%
40%
20%
0%

Nominal size α = 5 %
Parametrization (ii)

4 0

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
M

m
M

M
m

m
M
M
m

m
M

MMMMM
M M
M M mmmmmmmmmmmmmm
mmmm

M
m

m
M
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
M

MMMMMMM
MMM
MMMM
M
mmmmmmmmmm
M
mM
m
mMmmmm
m

m
M
m
M

m
M

m
M

M

MMMM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M
mmmmm

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
M

m

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

m
M

m
M

m
M

MMMMMMMMMM
MMMM
MM
M mmmmmmmmmmmmmmm
mm

m
M

m
M

m
M

mmmmmmmmmmmmmm

M
m

MMMMMMM
MM
MMMMM
M
M
mMmmmmmmmmmmmmmmm
mM
m

m
M

m
M

MMMMMMMMMMMMMM

M
m
m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
M

m
M

MM
MM
mmmm

MMMMMMM

M
m

m
M

m
M

4

Test statistics

MM

m
M

3

m
M

mM
MM
mm
mm
mM
M
mM
mm
MM
MM
Mm
mm
mM
mM
MM
mm
M

m
M

2

m
M

m
M

m
M

1

MMMMMMMMMMMMMMMM
M
mM
mM
mm
mmm
m
M
mmmmmmmmmmmm

m
M

m
M

m
M

m
M

Parametrization (iii)
4 0

MMMMMMMM
MM MMMMMMM
MmmM
m
mmmmmmmm
mmmmmmmm

m
M

m
M

3

m
M
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
mM
Mm
MM
Mm
Mm
mm
mM
M

m
M

2

M MMMMMMMMMMMMMMM
M
mM
mmM
m
M
m
mmmmm
mmm
mmmmmmm

m
M

m
M

1

m
M

m
M

m
M

Fig 11. The power curves under different parametrizations. The nominal size is α “ 5%.
In each plot, the vertical axis denotes the power, whereas the horizontal axis denotes the
value of δ “ µ2 ´ µ1 . The legend in Figure 9 also applies here.

Table 8
The range of empirical size rmin α
p, max α
ps in percentage, where max and min are taken
over the three parametrizations. Only one value is recorded for those tests that are
invariant to parametrization. The nominal size is α “ 5%.
Range of empirical size: rmin α
p, max α
ps{%
pn, mq

p1600, 3q

p400, 3q

p100, 3q

p100, 10q

p100, 30q

W-1
W-2
W-3
W-4

r5.62, 5.71s
r5.93, 6.05s
r5.81, 6.03s
r5.62, 5.71s

r5.30, 6.03s
r6.08, 7.18s
r6.01, 6.98s
r5.30, 6.03s

r3.22, 6.20s
r5.52, 8.69s
r5.37, 8.28s
r3.22, 6.20s

r1.64, 4.81s
r4.42, 8.47s
r4.20, 7.67s
r1.64, 4.81s

r1.37, 5.00s
r4.20, 8.50s
r4.10, 7.50s
r1.37, 5.00s

L-1
L-2
L-3
L-4
L-5
L-0

r5.57, 6.15s
r5.52, 6.10s
5.76
5.76
4.96
5.03

r6.37, 6.57s
r6.37, 6.52s
6.37
6.37
5.32
5.03

r5.88, 6.47s
r5.88, 7.47s
5.42
5.42
4.93
5.57

r4.39, 5.66s
r4.39, 5.66s
3.78
3.78
4.79
5.57

r4.22, 5.32s
r4.22, 5.32s
3.71
3.71
4.54
5.57

44

K.W. CHAN AND X.-L. MENG

Parametrization (i)
δ = µ2 − µ1
1
2
3

m = 30 , n = 100 m = 10 , n = 100 m = 3 , n = 100

Parametrization (ii)

4 0

1

2

3

Parametrization (iii)
4 0

1

2

3

4

40
Power/ Size

m = 3 , n = 400 m = 3 , n = 1600

0

Nominal size α = 5 %

30
20

m
M

10
0
40

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
M

m
M

m
M

M
mM
mM
MMMMMMMMMMMMMMMM
mm
mmmmmmmm
mmmmmmm

m
M

m
M

MMMMMMMMMMMMMM
M
mM
mM
mM
mM
mmmmmmmmmmmmmmm

m
M

30
20

mM
mM
mM
mM
mM
mM
mM
m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M

10
0
40

m
M

m
M

m
M

m
M

MM MMMMMMMMMMMMMMM
M
mmmM
mmmmmmmmmmmmmmmm

MMMMMMMMMMMMMMMMMM
M
mmmmmmmmmmmmmmmmmmm
m
M

Test statistics

m
M

30
20

m

10
0
40

m
mM
M

m
M

mM
m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M
m
mM
M

MMMMMMMM
MMMMMM
mM
mmmmmmmmmmmmmmm
mM
mM
mM
M

m
mM
M

m
M

m
mM
M

m
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
M
mM

30
20

m
M

10
0
40

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
m
MM

m
M
m
mM
M

M
mm
M

m
M

MMMMMMMMM
MMMMMM
mM
mM
M
mmmmmmmmmmmmmmmm

m
M
m
mM
M

M

W-1
W-2
W-3
W-4
L-1: MR92
L-2: MR92(+)
L-3: Proposal
L-4: Proposal(+)
L-5: Proposal(R)
L-0: Benchmark

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
M

m
M

30
20

m
M

0

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
M

m
M

10
m
mM
M

m
M
m
mM
M

MM
MMMMMMMMM
MMMM
mmmmmmmmmmmm
mM
M
mM
mmmm

m
M

mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
mM
m
mM
mM
M

m
M
m
mM
M

Fig 12. The ratios of empirical power to empirical size under different parametrizations.
The nominal size is α “ 5%. In each plot, the vertical axis denotes the ratio, whereas the
horizontal axis denotes δ “ µ2 ´ µ1 . The legend in Figure 9 also applies here.

45

MULTIPLE IMPUTATION LR TESTS

Table 9
r L, D
p ` and D
p ♦ under different parametrizations in § 4.2. The p-values
The LRTs using D
S
S
prL , pp`
p♦
S and p
S for testing conditional independence are essentially zero.
Parametrization (i)
H0 : Conditional independence

H0 : Full independence

m

rrL , rpS` , rpS♦

r L, D
p `, D
p♦
D
S
S

prL , pp`
p♦
S ,p
S

rrL , rpS` , rpS♦

r L, D
p `, D
p♦
D
S
S

2
3
5
7
10
25
50

0.63, 0.64, 0.83
0.54, 0.54, 0.38
0.49, 0.48, 0.89
0.23, 0.23, 0.47
0.50, 0.50, 0.70
0.35, 0.35, 0.47
0.31, 0.31, 0.45

0.14, 0.14, 0.12
0.08, 0.08, 0.09
0.12, 0.12, 0.10
0.06, 0.06, 0.05
0.14, 0.14, 0.12
0.06, 0.06, 0.06
0.11, 0.11, 0.10

0.87, 0.87, 0.89
0.93, 0.93, 0.92
0.89, 0.89, 0.91
0.94, 0.94, 0.95
0.87, 0.87, 0.88
0.94, 0.94, 0.95
0.90, 0.90, 0.91

0.53, 0.53, 0.83
0.31, 0.31, 0.38
0.72, 0.72, 0.89
0.31, 0.31, 0.47
0.56, 0.56, 0.70
0.35, 0.35, 0.47
0.33, 0.33, 0.45

44.4, 44.4, 37.1
54.2, 54.2, 51.4
40.8, 40.8, 37.1
53.2, 53.2, 47.6
45.4, 45.4, 41.7
51.4, 51.4, 47.0
51.5, 51.5, 47.3

Parametrization (ii)
H0 : Conditional independence

H0 : Full independence

m

rrL , rpS` , rpS♦

p♦
r L, D
p `, D
D
S
S

p♦
prL , pp`
S
S ,p

rrL , rpS` , rpS♦

p♦
r L, D
p `, D
D
S
S

2
3
5
7
10
25
50

0.59, 0.64, 0.83
0.54, 0.54, 0.38
0.51, 0.48, 0.89
0.22, 0.23, 0.47
0.50, 0.50, 0.70
0.35, 0.35, 0.47
0.30, 0.31, 0.45

0.15, 0.14, 0.12
0.08, 0.08, 0.09
0.11, 0.12, 0.10
0.06, 0.06, 0.05
0.14, 0.14, 0.12
0.06, 0.06, 0.06
0.11, 0.11, 0.10

0.86, 0.87, 0.89
0.92, 0.93, 0.92
0.90, 0.89, 0.91
0.94, 0.94, 0.95
0.87, 0.87, 0.88
0.94, 0.94, 0.95
0.90, 0.90, 0.91

0.47, 0.53, 0.83
0.30, 0.31, 0.38
0.68, 0.72, 0.89
0.28, 0.31, 0.47
0.54, 0.56, 0.70
0.33, 0.35, 0.47
0.32, 0.33, 0.45

46.3, 44.4, 37.1
54.4, 54.2, 51.4
41.8, 40.8, 37.1
54.3, 53.2, 47.6
46.1, 45.4, 41.7
51.9, 51.4, 47.0
52.1, 51.5, 47.3

Parametrization (iii)
H0 : Conditional independence

H0 : Full independence

m

rrL , rpS` , rpS♦

r L, D
p `, D
p♦
D
S
S

prL , pp`
p♦
S ,p
S

rrL , rpS` , rpS♦

r L, D
p `, D
p♦
D
S
S

2
3
5
7
10
25
50

1.23, 0.64, 0.83
1.08, 0.54, 0.38
1.02, 0.48, 0.89
0.45, 0.23, 0.47
0.99, 0.50, 0.70
0.71, 0.35, 0.47
0.63, 0.31, 0.45

0.01, 0.14, 0.12
´0.07, 0.08, 0.09
´0.09, 0.12, 0.10
´0.07, 0.06, 0.05
´0.10, 0.14, 0.12
´0.14, 0.06, 0.06
´0.10, 0.11, 0.10

0.99, 0.87, 0.89
1.00, 0.93, 0.92
1.00, 0.89, 0.91
1.00, 0.94, 0.95
1.00, 0.87, 0.88
1.00, 0.94, 0.95
1.00, 0.90, 0.91

0.98, 0.53, 0.83
0.61, 0.31, 0.38
1.40, 0.72, 0.89
0.58, 0.31, 0.47
1.09, 0.56, 0.70
0.68, 0.35, 0.47
0.65, 0.33, 0.45

34.2, 44.4, 37.1
43.9, 54.2, 51.4
29.0, 40.8, 37.1
43.9, 53.2, 47.6
33.7, 45.4, 41.7
41.0, 51.4, 47.0
41.3, 51.5, 47.3

