On the Correction of the Asymptotic Distribution of the Likelihood Ratio Statistic If Nuisance Parameters Are Estimated Based on an External Source

1 Introduction

Statistical models often contain multiple unknown parameters, some of interest, and some not, referred to as nuisance parameters. The presence of the latter parameters may make estimation, hypothesis testing and the construction of confidence intervals for the parameters of interest more complex, because more parameters must be estimated, possibly of infinite dimension.

Estimates of the nuisance parameters may be available in the literature or can sometimes be obtained from an independent dataset. Inserting these estimates into the model and next performing inference on the parameters of interest may simplify the statistical inference enormously. However, treating the nuisance parameters as known may affect the level of a test or the coverage of a confidence interval, making these higher and lower than nominal, respectively. In this paper, we investigate these effects for likelihood ratio procedures.

For standard statistical models, the likelihood ratio statistic for a finite-dimensional parameter is well known to be asymptotically distributed under the null hypothesis as a chi-squared random variable with degrees of freedom equal to the number of free parameters [1–3]. This result extends to likelihood ratio statistics for the parameter of interest in semiparametric models, as shown in various settings in Thomas and Grunkemeier [4] (right censoring), Murphy and van der Vaart [5] (semiparametric models), Banerjee and Wellner [6] (current status data), Banerjee [7] (monotone missingness), Kosorok [8], with empirical ratio likelihood statistics for functionals [9–11] as an important special case. The likelihood ratio statistic in the latter semiparametric models can be viewed as the likelihood ratio statistic for the parameter of interest in the model with the nuisance parameter replaced by its maximum likelihood estimator, based on the same data; equivalently as the likelihood ratio statistic based on the profile likelihood Murphy and van der Vaart [12]. With this estimate of the nuisance parameter, the asymptotic chi-squared distribution is retained, much as when testing only a part of the parameter in a finite-dimensional model. In contrast, it is shown below that the distribution changes if the nuisance is estimated using an independent sample of observations.

Naturally, this change will be small or even negligible if the additional sample is much bigger than the sample used to construct the likelihood ratio statistics. An example is estimating the disease-allele frequency in an inbred population, described in Example 4 below (see also Jonker et al. [13] and Teeuw et al. [14]). In this example, only 20 individuals were available for estimating the nuisance parameter, and 80 individuals to estimate the parameter of interest. It turns out that for two different approaches using a likelihood ratio, in the one case the distribution is well approximated by a chi-squared distribution as if the nuisance were known, whereas in the other case a correction is needed. One motivation of the present paper is to gain insight into the fact that not only the relative sample sizes determine the amount of correction needed but also the structure of the model. This is explained below in terms of a relevant Fisher information quantity.

The original motivation for the present paper was a study on the genetic determinants of the onset of migraine, described in Jonker et al. [15]. In this example, a baseline survival distribution was estimated using an additional, much larger dataset. As these data were subject to current status censoring, this estimate was relatively imprecise, of order comparable to the estimation error of the direct data, notwithstanding the size of the dataset. However, in a simulation study, it was found that substitution of the estimator in the likelihood ratio test caused only a minor error in the level of the test. This finding is explained by the detailed derivation of the likelihood ratio statistic in this paper, which reveals that the estimate for the baseline survival function enters (asymptotically) into the likelihood ratio through a smoothing operation, which repairs its inefficiency as an estimate for the survival function itself. See Example 3 below for details.

To our knowledge, the asymptotic distribution of the likelihood ratio statistic in these and other examples was not studied before. On the other hand, several authors have studied plug-in estimators within the context of empirical likelihood for parameters defined by estimating equations [9], among others Qin and Jing [16, 17], Li and Wang [18] and Wang and Jing [19]. In particular, Hjort et al. [20] study a similar problem in the setting of empirical likelihood. Besides the limit distribution, these authors also discuss a bootstrap approximation and allow an increasing number of estimating equations.

In the next section, we describe the theorem and give the proof. Thereafter we consider four examples; two theoretical examples, for a parametric and a semiparametric model, and two practical problems, for a large and a small dataset for estimating the nuisance parameter.

2 The likelihood ratio test

Suppose we observe a random sample X1,…,Xn from a distribution Pθ,η with corresponding density pθ,η indexed by two parameters: a vector θ∈Θ⊂ℝd and a second parameter η belonging to an arbitrary parameter set H. We let η˜ be an estimator for η based on a second sample Y1,…,Ym of observations that is independent of X1,…,Xn and possesses a distribution depending on η only.

In the following, we will test the two-sided hypothesis H0:θ∈Θ0 which is of the form H0:Σθ∈{(c1,…,cd−k,xd−k+1,…,xd)T:xd−k+1,…,xd∈ℝ} for c1,…,cd−k∈ℝ given and Σ a given d×d invertible matrix. If Σ equals the identity matrix and c1=⋯=cd−k=0 it is tested whether the first d−k coordinates of θ equal zero. We assume that, n(Θ−θ0) and n(Θ0−θ0) tend to d and k dimensional spaces. The convergence of the sets n(Θ−θ0)={n(θ−θ0):θ∈Θ} and n(Θ0−θ0)={n(θ−θ0):θ∈Θ0} is defined as follows (see van der Vaart [2], chapter 16). A set Hn converges to a set H if H is the set of all limits of all converging sequences hn with hn∈Hn for all n and, if h=limihni for n1<n2<… with hni∈Hni for every i, then h is contained in H.

To test the null hypothesis, we use the likelihood ratio statistic in which we replace the parameter η by its estimator η˜:

In the theorem below, we will determine the asymptotic distribution of the test-statistic under θ=θ0 with θ0∈Θ0⊂Θ.

For k=0 (the dimension of the null space is zero), the asymptotic distribution has a nice form and is given directly after the assumptions. For k>0 the values di can be found in the proof of the theorem.

If n/m→0, λ=0 and Λm,n converges in distribution to a chi-squared distribution with d−k degrees of freedom. So, if m>>n the distribution of the test-statistic is hardly affected if the estimated nuisance parameters are seen as the true values and the level of the test is close to the supposed level under the assumption of known nuisance parameters. However, also if m<n and the values of di are close to zero, the distribution of the likelihood ratio statistic may be close to the chi-squared distribution with d−k degrees of freedom (see Example 4).

In the theorem, the exact asymptotic distribution is not given. However, if the dimension of the null space Θ0 is zero (Θ0 consists of a single point, k=0) the asymptotic distribution of the likelihood ratio statistic has a nice form: the likelihood ratio statistic Λm,n is asymptotically distributed as ∑i=1d(1+λdi)Zi2 with λ=limm,nn/m and d1,…,dd the eigenvalues of the covariance matrix Iθ0,η0−1/2Jθ0,η0Iθ0,η0−1/2 and Z12,…,Zd2 independent random variables with a chi-squared one distribution (the matrices Iθ0,η0 and Jθ0,η0 were defined in the assumptions 3 and 5). If d=1 then, Iθ0,η0−1/2Jθ0,η0Iθ0,η0−1/2=Jθ0,η0Iθ0,η0−1 and the asymptotic distribution simplifies to (1+λJθ0,η0Iθ0,η0−1)Z2; a chi-squared one distribution multiplied with a correction factor bigger than one.

For a statistical model probably the most difficult assumption to check is the asymptotic normality of mPθ0,η˜l˙θ0,η0 in case η˜ is a non-parametric estimator of η which converges at a lower rate than m1/2 (assumption 5). For example, if η˜ is the non-parametric maximum likelihood estimator (NPMLE) of a distribution function based on interval censored data, we know that η˜ converges pointwise at rate n1/3, but, under some assumptions, some functionals of it, for instance ∫x(η˜(x)−η(x))dx, converge at rate n1/2 to normal distributions (theorem 5.1 in Huang and Wellner [21]). So, even if the convergence rate of η˜ is lower than m1/2, it may be possible to prove that mPθ0,η˜l˙θ0,η0 is asymptotically normal.

In order to test the hypothesis H0:θ∈Θ0 with test-statistic Λm,n, we have to compute the (asymptotic) distribution of Λm,n under the null hypothesis. Remember that θ0 is a point in Θ and Θ0 so that n(Θ−θ0) and n(Θ0−θ0) tend to d and k dimensional spaces. Now write

for θˆ the maximum likelihood estimator of θ and θˆ0 the restricted maximum likelihood; restricted to the parameter space Θ0 after inserting η˜ into the likelihood. Under assumption 2 and under θ0∈Θ0 the maxima over θ∈Θ and θ∈Θ0 are in 1/n-shrinking neighborhoods of θ0. That means that θˆ=θ0+hˆ/n for some value hˆ∈n(Θ−θ0) and θˆ0=θ0+hˆ/n for hˆ∈n(Θ0−θ0). Instead of maximizing the likelihood ∏i=1npθ,η˜ in the test-statistic with respect to θ∈Θ or θ∈Θ0, the product ∏i=1npθ0+h/n,η˜ can be maximized with respect to h∈Hn=n(Θ−θ0) or h∈Hn,0=n(Θ0−θ0). This yields the likelihood ratio statistic

By a Taylor expansion and by assumption 7

with

The term in eq. (2) should be maximized with respect to h.

We consider the three terms in eq. (2) separately. By assumption 4 the sequence Gnl˙θ0,η˜ converges in distribution to N(0,Iθ0,η0) and −Pnl¨θ0,η˜(Xi) converges in probability to Iθ0,η0 by assumption 3. Left is the second term in eq. (2). For any η in its parameter space H, Pθ0,ηl˙θ0,η=0. Write

By assumption 6 the last term at the right-hand side converges in probability to zero. Thus

what is, by assumption 5, asymptotically normal with mean zero and covariance matrix λJθ0,η0 with λ equal to the limit of n/m if m,n→∞.

Combining the previous results, by assumption 4 and by independence of the samples X1,…,Xn and Y1,…,Ym the sum

converges to the normal distribution N(0,Iθ0,η0+λJθ0,η0). The same computations hold for the second term of eq. (1).

The following is similar to the derivation of the asymptotic distribution of the likelihood ratio statistic, if there is no nuisance parameter (van der Vaart [2], Theorem 16.7). The differentiability of logpθ0+h/n,η˜(Xi) with respect to h implies that

converges in probability to zero for every h. By assumption 7 (the remainder term of the Taylor expansion converges uniformly in probability to zero) and the fact that

there is also uniform convergence of Z on a bounded set:

The convergence in the previous display also holds if M=Mn increases to infinity sufficiently slowly. Since θˆ and θˆ0 are both n consistent, both estimators are within a ball of radius Mn/n around θ0 with probability tending to 1. Then hˆ is within a ball of radius Mn around zero with probability tending to 1. Consequently, the limit distribution of Λm,n does not change if Hn or Hn,0 are replaced by the intersection of these sets and the ball around zero with radius Mn. These intersected sets also converge to H and H0, respectively. By uniform convergence we conclude that for Gn,η˜=Gnl˙θ0,η˜+nPθ0,η0l˙θ0,η˜

where the last equality follows from Lemma 7.13 in van der Vaart [2] with H and H0 the limit sets of Hn and Hn,0 as explained at the beginning of this section and G distributed as N(0,Iθ0,η0+λJθ0,η0), the limit distribution of Gn,η˜ (as was seen before).

Define Z:=Iθ0,η0−1/2G, then Z possesses a normal distribution with mean zero and covariance matrix

with I the d×d identity matrix. Furthermore, since θ0 is an inner point of Θ, the set H is the full space in ℝk and the second term in eq. (3) equals zero, and

for m,n to infinity.

Suppose that k=0, then H0 equals the singleton {0} in a d-dimensional space and the asymptotic distribution of Λm,n is equal to the law of ∥Z∥2 with Z normally distributed with mean zero and covariance matrix Qθ0,η0. The covariance matrix Qθ0,η0 is a symmetric positive definite matrix and can therefore be factorized as Qθ0,η01/2=UTD1/2U with U an orthonormal matrix and D a diagonal matrix with the eigenvalues of Qθ0,η0 on the diagonal. Then Qθ0,η0=Qθ0,η01/2Qθ0,η01/2 with Qθ0,η01/2=UTD1/2U and

Now, Qθ0,η0−1/2Z has a d dimensional normal distribution with mean zero and covariance matrix the identity matrix. Since U is an orthonormal matrix and only rotates the vector and thus ∥Ux∥=∥x∥, the last term in the previous display equals ∥D1/2U(Qθ0,η0−1/2Z)∥2. Furthermore, U(Qθ0,η0−1/2Z) and Qθ0,η0−1/2Z are equally distributed by symmetry of the standard normal distribution. Consequently

with e1,…,ed the diagonal elements of the matrix D and Z1,…,Zd independent standard normal distributed random variables. So, Zi2 has a chi-squared distribution with one degree of freedom and Λm,n converges in distribution to a linear combination of independent chi-squared random variables (all with one degree of freedom): ∑i=1deiZi2. The eigenvalues e1,…,ed of the matrix Qθ0,η0 equal ei=1+λdi,i=1,…,d with d1,…,dd the eigenvalues of the covariance matrix Iθ0,η0−1/2Jθ0,η0Iθ0,η0−1/2.

Now, let k>0. Remember that the asymptotic distribution of Λm,n equals the law of ∥Z−Iθ0,η01/2H0∥2 for H0 a subspace of dimension k. If P is the orthogonal projection onto the orthocomplement of Iθ0,η01/2H0, which is a d−k dimensional linear subspace of ℝd, then

In the situation that k=0 the covariance matrix Qθ0,η01/2 of Z was factorized in order to rewrite the quadratic norm ∥Z∥2 into a sum. The same can now be done for the covariance matrix PQθ0,η0PT of W. It follows that ∥W∥2 equals the sum ∑i=1deiZi2 with Z12,…,Zd2 chi-squared one distributed and e1,…,ed equal to the eigenvalues of PQθ0,η0PT. We still have to compute these eigenvalues; we will show that k eigenvalues equal zero and that the remaining d−k eigenvalues are strictly positive and equal to ei=1+λdi,i=1,…,d−k with d1,…,dd−k the strictly positive eigenvalues of PSθ0,η0PT with Sθ0,η0:=Iθ0,η0−1/2Jθ0,η0Iθ0,η0−1/2. Then, the sum reduces to ∑i=1d−k(1+λdi)Zi2.

The derivation of the eigenvalues is as follows. The covariance matrix PQθ0,η0PT=PPT+λPSθ0,η0PT=P+λPSθ0,η0PT. It vanishes on Iθ0,η01/2H0 and acts as the strictly positive-definite operator I+λPSθ0,η0PT on the orthocomplement of Iθ0,η01/2H0. It follows that PQθ0,η0PT has 0 as an eigenvalue of multiplicity k (the dimension of Iθ0,η01/2H0) and has d−k strictly positive eigenvalues of the form ei=1+λdi, where di is a positive eigenvalue of PSθ0,η0PT. □

If d=1 then Iθ0,η0−1/2Jθ0,η0Iθ0,η0−1/2=Jθ0,η0Iθ0,η0−1 and the asymptotic distribution of the likelihood ratio statistic simplifies to (1+λJθ0,η0Iθ0,η0−1)Z2 where Z2 follows the chi-squared distribution with one degree of freedom. The matrix Jθ0,η0 is the asymptotic covariance matrix of mPθ0,η˜l˙θ0,η0 with Pθ0,η˜l˙θ0,η0 the “plug in” estimator for the parameter η↦Pθ0,ηl˙θ0,η0 based on the sample Y1,…,Ym. The Fisher information matrix Iθ0,η0 is the asymptotic covariance matrix of the alternative estimator Pnl˙θ0,η0 based on X1,…,Xn. So, the term Jθ0,η0Iθ0,η0−1 measures the relative efficiency at η=η0 for estimating Pθ0,ηl˙θ0,η0 using one observation Y compared to one observation X. Then Jθ0,η0/m and Iθ0,η0/n are the variances of Pθ0,η˜l˙θ0,η0 and Pnl˙θ0,η0 based on the samples Y1,…,Ym and X1,…,Xn, respectively. The fraction (Jθ0,η0/m)(Iθ0,η0/n)−1=(n/m)(Jθ0,η0/Iθ0,η0) measures the relative efficiency of the two samples for estimating Pθ0,ηl˙θ0,η0. This fraction is approximately equal to λJθ0,η0Iθ0,η0−1 as in the asymptotic distribution of the likelihood ratio statistic. So, not only the sample sizes of the two samples are important but also the estimation precision of the estimators (the asymptotic variances). This is illustrated in Example 4.

Next we consider the asymptotic distribution of the likelihood ratio test-statistic for the situation that Θ is a half-space and Θ0 its boundary. Suppose we want to test the null hypothesis H0:θ∈Θ0 against H1:θ∈Θ/Θ0. Let θ0∈Θ0 and define, as before, Hn=n(Θ−θ0) and Hn,0=n(Θ0−θ0). Assume that Hn and Hn,0 tend to a half space H in ℝd and its boundary H0, respectively, like Θ0 and Θ.

Under the assumptions mentioned earlier, Λm,n converges, under θ=θ0, in distribution to a mixture of a point mass at 0 and the distribution of a linear combination of d−k independent chi-squared distributed random variables, each with one degree of freedom. The coefficients of the linear combination are as described before.

The proof is very similar up to the following. Like in the previous proof

but now the term ∥Iθ0,η0−1/2G−Iθ0,η01/2H∥ is not equal to zero with probability one. The sets H′:=Iθ0,η01/2H and H0′:=Iθ0,η01/2H0 are a halfspace and its boundary, because this is true for H and H0 by assumption, and the boundary contains 0, because θ0∈Θ0. If Z:=Iθ0,η0−1/2G falls in the halfspace H′, the second term in Λ vanishes, and Λ equals the square distance of Z to the boundary of the halfspace H′. If Z falls outside H′, then Λ equals zero, because in that case the distance of Z to H′ is equal to the distance to its boundary, as the projection of Z onto H′ is on the boundary of H′. Thus Λ is equal to zero with probability that Z falls outside H′. If Z falls into H′, the right term vanishes and Λ=∥Iθ0,η0−1/2G−Iθ0,η01/2H0∥2. This term can be considered as before.