Power analysis, sample size calculation for testing the largest binomial probability

Abstract

A procedure is developed for power analysis and sample size calculation for a class of complex testing problems regarding the largest binomial probability under a combination of treatments. It is shown that the asymptotic null distribution of the likelihood-ratio statistic is not parameter-free, but ${\chi }_1^2$ is a conservative asymptotic null distribution. A nonlinear Gauss-Seidel algorithm is proposed to uniquely determine the alternative for the power and sample size calculation given the baseline binomial probability. An example from an animal clinical trial is discussed.

Keywords:

• Asymptotic null distribution
• binomial probability
• complex hypotheses
• Gauss-Seidel
• logistic regression
• power
• sample size
• tests

1. Introduction

In biological and medical research, it is often necessary to perform power analysis, or determine the sample size required for binomial trials involving multiple treatments. For example, such an analysis/calculation is required by the National Institutes of Health (NIH) for research grant applications. One particular type of questions that need to be answered is regarding the largest binomail probability under these treatments. Suppose that there are r treatments, denoted by $1,\dots ,r$. Here the ‘treatment’ can be a real treatment (e.g., drug), or a factor that has two levels (e.g., male/female). Let $x_q$ be the indicator for the qth treatment (0 or 1), $1\le q\le r$. Suppose that n independent trials are run under each combination of the treatments, $x_1,\dots ,x_r$, resulting a binary outcome for each trial (1 – success, 0 – failure). Let $N(x_1,\dots ,x_r)$ be the total number of successes under $x_1,\dots ,x_r$. It is assumed that $N(x_1,\dots ,x_r)$ has a Binomia$\{n,p(x_1,\dots ,x_r\left)\right\}$ distribution, where $p(x_1,\dots ,x_r)$ is the probability of success in a single trial under $x_1,\dots ,x_r$. It is believed that all the treatments have at least nonnegative effects, so $p(0,\dots ,0)$ is the smallest probability and $p(1,\dots ,1)$ is the largest. The question is to determine the sample size, n, so that one has at least the power $1-\gamma$ to prove the case that $p(1,\dots ,1)$ is higher than the rest of the probabilities, if it is indeed higher by at least as much as δ. In some cases, the researcher has a target sample size (e.g., the maximum under the budget constraint). The question then is to perform the power analysis for the given sample size. We illustrate with an example.

Example 1.1

Researchers at the Oregon Health & Science University (OHSU) were preparing to meet the deadline of a grant submission. One of the research aims of the grant proposal had to do with comparative transplantations via animal clinical trials (i.e., mice) to determine the best overall protocol, which can then be further optimised. Successful liver repopulation is defined as greater than 70% cell replacement and successful blood reconstitution is defined as greater than 50% human cells in the bone marrow. From previous studies, it was known that these levels of liver repopulation could be achieved in 20–80% of transplanted mice with a good hepatocyte donor, the average being 50–60%. For cord blood transplants into neonates, about 75% of the mice reach the desired human repopulation levels, again with a range of 50–90% between experiments. Using these numbers from the single cell type transplants, we made the assumption that a successful protocol would yield about 30% success in double-repopulation. Here the protocol involved a treatment high and low dose levels, and a control factor of youth and grown-up mice. In other words, there are two treatments, $x_1=1$ for high dose and $x_1=0$ for low dose; $x_2=1$ for youth mice and $x_2=0$ for grown-up one. The trials were to be carried out independently on n different mice, resulting a binomial proportion, under each protocol. It was determined that n should be no more than 30 due to the budget constraint. An initial request was made to perform a power analysis in order to detect a 10% difference that separates the success rates of the ‘optimal’ protocol (i.e., with high dose and youth mice) and the rest.

A naive approach to the power analysis/sample size calculation would be to compare $p(1,\dots ,1)$ with each of the other probabilities (actually, only those with exactly one of the treatment indicators equal to zero), and perform a two-sample t-test for the difference in proportions. However, this approach is low power, and often results in a sample size that the researcher cannot afford. The inefficiency of the naive approach is not surprising, because only a (small) portion of the data are involved in the two-sample t-test (for each comparison). One can do better by utilising the entire data in the analysis. To do so we need to assume a model for the binomial probability. Suppose that (1) $p(x_1,\dots ,x_r)=h({\text{β}}_0+{\text{β}}_1{x}_1+\cdots +{\text{β}}_r{x}_r),$(1) where $h(\cdot )$ satisfies $0<h\left(x\right)<1$ and is strictly increasing. A well-known example of h is the logistic function, $h\left(x\right)={\mathrm{e}}^x/(1+{\mathrm{e}}^x)$. Under the assumed model, the problem of interest can be expressed more precisely as testing the hypothesis (2) ${\mathrm{H}}_0:\text{at least one of the}{\text{β}}_j,1\le j\le r,\mathrm{i}\mathrm{s}\le 0$(2) versus ${\mathrm{H}}_1:{\text{β}}_j>0,1\le j\le r$. Natually, the likelihood-ratio test is considered. The latter is based on the test statistic (3) $\mathcal{L}=2(\hat{l}-{\hat{l}}_0),$(3) where $\hat{l}$ is the maximised log-likelihood and ${\hat{l}}_0$ the maximised log-likelihood under ${\mathrm{H}}_0$. A first step for the likelihood-ratio test (LRT) would be to determine the critical value, $c_{\alpha }$, corresponding to the given level of significance α, such that () $\underset{\text{β}\in {\mathrm{H}}_0}{sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>c_{\alpha })\le \alpha ,$() where $\text{β}=({\text{β}}_0,\dots ,{\text{β}}_r{)}^{\prime }$ and ${\mathrm{P}}_{\text{β}}$ is the probability distribution given that β is the true vector of parameters. If the log-likelihood function is log-concave, which is the case, for example, for the logistic regression, the event inside the probability of (4() $\underset{\text{β}\in {\mathrm{H}}_0}{sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>c_{\alpha })\le \alpha ,$() ) implies that the maximum for ${\hat{l}}_0$ must take place on the boundary of ${\mathrm{H}}_0$, provided that $c_{\alpha }>0$. Therefore, the critical value is computed by considering the boundary of ${\mathrm{H}}_0$, which is a subset of (5) ${\tilde{\mathrm{H}}}_0:{\text{β}}_j=0\text{for some}1\le j\le r.$(5) Unfortunately, even with the latest simplification, the asymptotic null distribution of the LRT is not parameter-free. This is shown in the next section. On the other hand, the arguments also show that a conservative asymptotic null distribution (CAND) for the LRT is ${\chi }_1^2$, which is parameter-free. Here the CAND is in the sense that (6) $\underset{\text{β}\in {\tilde{\mathrm{H}}}_0}{sup}\underset{n\to \mathrm{\infty }}{lim sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)=\alpha .$(6)

The main objective of the current paper is to determine the sample size, n, for the LRT so that the test will have the designated power, or to obtain the power of the LRT, under a given sample size. Although standard power and sample size problems in logistic regression are well studied (e.g., Alam, Rao, & Cheng, 2010; Borenstein, Rothstein, & Cohen, 2001; Demidenko, 2007; Hsieh, Bloch, & Larsen, 1998; Novikov, Fund, & Freedman, 2010; Whittemore, 1981), to the best of our knowledge, the kind of problems that we are dealing with have not been addressed. Clearly, the power of the test depends on the alternative, and there are infinitely many possible alternatives to (2(2) ${\mathrm{H}}_0:\text{at least one of the}{\text{β}}_j,1\le j\le r,\mathrm{i}\mathrm{s}\le 0$(2) ), that is, in ${\mathrm{H}}_1$. On the other hand, a practitioner would prefer a ‘short answer’, as opposed to something that is case-by-case. This issue is addressed in Section 3, where a unique alternative, based on a reasonable argument, is determined. A simple Gauss-Seidel type algorithm is proposed to compute the alternative. A Monte Carlo procedure is then proposed to compute the power or sample size. A real-life application is considered in Section 4. Technical details are deferred to Appendix.

2. Asymptotic null distribution

In the standard situation, the LRT is known to have an asymptotic ${\chi }^2$ distribution, under the null hypothesis, with a certain degrees of freedom that does not depend on the parameter under the null. For example, if, instead of (2(2) ${\mathrm{H}}_0:\text{at least one of the}{\text{β}}_j,1\le j\le r,\mathrm{i}\mathrm{s}\le 0$(2) ), one were to test (7) ${\mathrm{H}}_0^j:{\text{β}}_j=0$(7) versus ${\mathrm{H}}_1^j:{\text{β}}_j\ne 0$ for a fixed $1\le j\le r$, then the asymptotic distribution of the LRT is ${\chi }_1^2$, regardless of the value of the true ${\text{β}}_k,k\ne j$, as long as the true ${\text{β}}_j$ is zero. In other words, the asymptotic null distribution is parameter-free. However, this is not the case for testing (2(2) ${\mathrm{H}}_0:\text{at least one of the}{\text{β}}_j,1\le j\le r,\mathrm{i}\mathrm{s}\le 0$(2) ). We show this for the case of logistic regression with r=2.

Write the log-likelihood function as $l({\text{β}}_0,{\text{β}}_1,{\text{β}}_2)$. Let $({\widehat{\text{β}}}_0,{\widehat{\text{β}}}_1,{\widehat{\text{β}}}_2{)}^{\prime }$, $({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0{)}^{\prime }$, and $({\widehat{\text{β}}}_0^{\left[1\right]},0,{\widehat{\text{β}}}_2^{\left[1\right]}{)}^{\prime }$ be the maximiser of l without constraint, over ${\mathrm{H}}_0^2=\left\{\right({\text{β}}_0,{\text{β}}_1,0):{\text{β}}_0,{\text{β}}_1\in R\}$, and over ${\mathrm{H}}_0^1=\left\{\right({\text{β}}_0,0,{\text{β}}_2):{\text{β}}_0,{\text{β}}_2\in R\}$, respectively. Suppose that the true parameter vector is $({\text{β}}_0,{\text{β}}_1,0)\in {\mathrm{H}}_0^2$, where ${\text{β}}_1\ne 0$. By the standard asymptotic theory, we have ${\widehat{\text{β}}}_j^{\left[2\right]}\stackrel{\mathrm{P}}{⟶}{\text{β}}_j,j=0,1$, as $n\to \mathrm{\infty }$. On the other hand, by White (1982), we have ${\widehat{\text{β}}}_j^{\left[1\right]}\stackrel{\mathrm{P}}{⟶}{\text{β}}_j^{\left[1\right]},j=0,2$, as $n\to \mathrm{\infty }$ for some ${\text{β}}_0^{\left[1\right]}$ and ${\text{β}}_2^{\left[1\right]}$. Thus, by the Taylor expansion, we have (8) $\begin{array}{rl}& \frac{l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)}{n}=\frac{l({\text{β}}_0,{\text{β}}_1,0)}{n}\\ & +\frac{1}{n}\sum _{j=0,1}\left\{{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }{\text{β}}_j}\right|}_{{\widetilde{\text{β}}}^{\left(2\right)}}\right\}\left({\widehat{\text{β}}}_j^{\left[2\right]}-{\text{β}}_j\right),\end{array}$(8) (9) $\begin{array}{rl}& \frac{l({\widehat{\text{β}}}_0^{\left[1\right]},0,{\widehat{\text{β}}}_2^{\left[1\right]})}{n}=\frac{l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & +\frac{1}{n}\sum _{j=0,2}\left\{{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }{\text{β}}_j}\right|}_{{\widetilde{\text{β}}}^{\left(1\right)}}\right\}\left({\widehat{\text{β}}}_j^{\left[1\right]}-{\text{β}}_j^{\left[1\right]}\right),\end{array}$(9) for some ${\widetilde{\text{β}}}^{\left(j\right)},j=1,2$. It is easy to show that the partial derivatives are uniformly bounded when divided by n, so the second terms on the right sides of (8(8) $\begin{array}{rl}& \frac{l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)}{n}=\frac{l({\text{β}}_0,{\text{β}}_1,0)}{n}\\ & +\frac{1}{n}\sum _{j=0,1}\left\{{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }{\text{β}}_j}\right|}_{{\widetilde{\text{β}}}^{\left(2\right)}}\right\}\left({\widehat{\text{β}}}_j^{\left[2\right]}-{\text{β}}_j\right),\end{array}$(8) ) and (9(9) $\begin{array}{rl}& \frac{l({\widehat{\text{β}}}_0^{\left[1\right]},0,{\widehat{\text{β}}}_2^{\left[1\right]})}{n}=\frac{l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & +\frac{1}{n}\sum _{j=0,2}\left\{{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }{\text{β}}_j}\right|}_{{\widetilde{\text{β}}}^{\left(1\right)}}\right\}\left({\widehat{\text{β}}}_j^{\left[1\right]}-{\text{β}}_j^{\left[1\right]}\right),\end{array}$(9) ) are $o_{\mathrm{P}}\left(1\right)$. As for the first terms, it is easy to show that, for any β, we have $\begin{array}{rl}l\left(\text{β}\right)& =c+\sum _{x_1,x_2=0,1}\left\{N(x_1,x_2)({\text{β}}_0+{\text{β}}_1{x}_1+{\text{β}}_2{x}_2)\right.\\ & \left.-n\log (1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1+{\text{β}}_2{x}_2})\right\},\end{array}$ where c does not depend on the parameter. Thus, we have $\begin{array}{rl}& \frac{l({\text{β}}_0,{\text{β}}_1,0)-l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & =\sum _{x_1,x_2=0,1}\left[\frac{N(x_1,x_2)}{n}\left\{\right({\text{β}}_0+{\text{β}}_1{x}_1)-({\text{β}}_0^{\left[1\right]}\right.\\ & \left.+{\text{β}}_2^{\left[1\right]}{x}_2\left)\right\}-\log \left(\frac{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0^1+{\text{β}}_2^1{x}_2}}\right)\right].\end{array}$ By the weak law of large numbers, we have $n^{-1}N(x_1,x_2)\stackrel{\mathrm{P}}{⟶}p(x_1,x_2)=h({\text{β}}_0+{\text{β}}_1{x}_1)$, with $h\left(x\right)={\mathrm{e}}^x/(1+{\mathrm{e}}^x)$, $x_1,x_2=0,1$. Thus, we have (10) $\begin{array}{rl}& \frac{l({\text{β}}_0,{\text{β}}_1,0)-l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & \stackrel{\mathrm{P}}{⟶}\sum _{x_1,x_2=0,1}\left[\frac{{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}\left\{\right({\text{β}}_0+{\text{β}}_1{x}_1)\right.\\ & \left.-({\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2)\}-\log \left(\frac{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2}}\right)\right].\end{array}$(10) For fixed $x_1,x_2\in \{0,1\}$, write $p_0=h({\text{β}}_0+{\text{β}}_1{x}_1)$ and $p_1=h\left({\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2\right)$. Then, by the inequality in Appendix A.1, we have $\begin{array}{rl}& \frac{{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}\left\{\right({\text{β}}_0+{\text{β}}_1{x}_1)-({\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2\left)\right\}\\ & -\log \left(\frac{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2}}\right)\\ & =p_0\log \left(\frac{p_0}{p_1}\right)+(1-p_0)\log \left(\frac{1-p_0}{1-p_1}\right)\ge 0,\end{array}$ with the equality holding if and only if $p_0=p_1$. It follows that the right side of (10(10) $\begin{array}{rl}& \frac{l({\text{β}}_0,{\text{β}}_1,0)-l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & \stackrel{\mathrm{P}}{⟶}\sum _{x_1,x_2=0,1}\left[\frac{{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}\left\{\right({\text{β}}_0+{\text{β}}_1{x}_1)\right.\\ & \left.-({\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2)\}-\log \left(\frac{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2}}\right)\right].\end{array}$(10) ) is positive unless $p_0=p_1$ for all $x_1,x_2=0,1$. Because the latter implies ${\text{β}}_1=0$, a contradiction, the right side of (10(10) $\begin{array}{rl}& \frac{l({\text{β}}_0,{\text{β}}_1,0)-l({\text{β}}_0^{\left[1\right]},0,{\text{β}}_2^{\left[1\right]})}{n}\\ & \stackrel{\mathrm{P}}{⟶}\sum _{x_1,x_2=0,1}\left[\frac{{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}\left\{\right({\text{β}}_0+{\text{β}}_1{x}_1)\right.\\ & \left.-({\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2)\}-\log \left(\frac{1+{\mathrm{e}}^{{\text{β}}_0+{\text{β}}_1{x}_1}}{1+{\mathrm{e}}^{{\text{β}}_0^{\left[1\right]}+{\text{β}}_2^{\left[1\right]}{x}_2}}\right)\right].\end{array}$(10) ) must be positive. Therefore, in conclusion, we have with probability tending to one that $l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)>l({\widehat{\text{β}}}_0^{\left[1\right]},0,{\widehat{\text{β}}}_2^{\left[1\right]})$, hence $\mathcal{L}=2\left\{l\right({\widehat{\text{β}}}_0,{\widehat{\text{β}}}_1,{\widehat{\text{β}}}_2)-l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0\left)\right\}\stackrel{\mathrm{d}}{⟶}{\chi }_1^2$, by the standard asymptotic result [see the note below (7(7) ${\mathrm{H}}_0^j:{\text{β}}_j=0$(7) )].

Now suppose that the true parameter vector is $({\text{β}}_0,0,0)$. Then, it can be shown (see Appendix A.2) that $\mathcal{L}\stackrel{\mathrm{d}}{⟶}\eta -{\eta }_1\vee {\eta }_2$, as $n\to \mathrm{\infty }$, where $\begin{array}{rl}\eta & ={\left(\begin{array}{c}{\xi }_0\\ {\xi }_1\\ {\xi }_2\end{array}\right)}^{\prime }{\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right)}^{-1}\left(\begin{array}{c}{\xi }_0\\ {\xi }_1\\ {\xi }_2\end{array}\right),\\ \left(\begin{array}{c}{\xi }_0\\ {\xi }_1\\ {\xi }_2\end{array}\right)& \sim N\left[\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right),\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right)\right],\\ {\eta }_j& ={\left(\begin{array}{c}{\xi }_0\\ {\xi }_j\end{array}\right)}^{\prime }{\left(\begin{array}{cc}4& 2\\ 2& 2\end{array}\right)}^{-1}\left(\begin{array}{c}{\xi }_0\\ {\xi }_j\end{array}\right),j=1,2.\end{array}$ Note that $\eta -{\eta }_1\vee {\eta }_2$ is not distributed as ${\chi }_1^2$. To see this, note that ${\eta }_1<{\eta }_2$ if and only if $({\xi }_1-{\xi }_0/2{)}^2<({\xi }_2-{\xi }_0/2{)}^2$. Because the pdf of ξ is positive everywhere, there is a positive probability that ${\eta }_1<{\eta }_2$, hence ${\eta }_1\vee {\eta }_2>{\eta }_1$. On the other hand, one always has ${\eta }_1\vee {\eta }_2\ge {\eta }_1$. It follows that $\mathrm{E}({\eta }_1\vee {\eta }_2)>\mathrm{E}\left({\eta }_1\right)$, hence $\mathrm{E}(\eta -{\eta }_1\vee {\eta }_2)<\mathrm{E}\left(\eta \right)-\mathrm{E}\left({\eta }_1\right)=3-2=1$.

The results so far in this section have shown that the asymptotic null distribution of $\mathcal{L}$ depends on the values of the true parameters, namely, if $\text{β}=({\text{β}}_0,{\text{β}}_1,0)$, where ${\text{β}}_1\ne 0$ [or $\text{β}=({\text{β}}_0,0,{\text{β}}_2)$, where ${\text{β}}_2\ne 0$, by a similar argument], the asymptotic null distribution is ${\chi }_1^2$; if $\text{β}=({\text{β}}_0,0,0)$, the asymptotic null distribution is not ${\chi }_1^2$. Thus, the asymptotic null distribution is not parameter-free.

Nevertheless, ${\chi }_1^2$ is, in general (not restricted to logistic and r=2), a CAND in the sense of (6(6) $\underset{\text{β}\in {\tilde{\mathrm{H}}}_0}{sup}\underset{n\to \mathrm{\infty }}{lim sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)=\alpha .$(6) ), provided that the log-likelihood is log-concave. This can be shown with a simple argument. For any $\text{β}\in {\tilde{\mathrm{H}}}_0$, there is a $1\le j\le r$ such that ${\text{β}}_j=0$. Due to the log-concavity, the event $\mathcal{L}>{\chi }_{1,\alpha }^2$ is the same as the event $2(\hat{l}-{\tilde{l}}_0)>{\chi }_{1,\alpha }^2$, where ${\tilde{l}}_0$ is the maximised log-likelihood over ${\tilde{\mathrm{H}}}_0$ [see the note above (5(5) ${\tilde{\mathrm{H}}}_0:{\text{β}}_j=0\text{for some}1\le j\le r.$(5) )]. On the other hand, we have $2(\hat{l}-{\tilde{l}}_0)\le 2(\hat{l}-{\hat{l}}_{0j})$, where ${\hat{l}}_{0j}$ is the maximised log-likelihood over ${\tilde{\mathrm{H}}}_{0j}=\{\text{β}:{\text{β}}_j=0\}$. Therefore, we have ${\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)={\mathrm{P}}_{\text{β}}\left\{2\right(\hat{l}-{\tilde{l}}_0)>{\chi }_{1,\alpha }^2\}\le {\mathrm{P}}_{\text{β}}\left\{2\right(\hat{l}-{\hat{l}}_{0j})>{\chi }_{1,\alpha }^2\}\to \alpha$, by the standard asymptotic result. Therefore, we have $\underset{n\to \mathrm{\infty }}{lim sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)\le \alpha$.

On the other hand, if $\text{β}\in {\tilde{\mathrm{H}}}_0$ such that ${\text{β}}_j=0$ while ${\text{β}}_k\ne 0,k\ne j$, by a similar argument as the above for the special case of logistic regression with r=2, it can be shown that, with ${\mathrm{P}}_{\text{β}}$ tending to one, we have $\mathcal{L}=2(\hat{l}-{\hat{l}}_{0j})\stackrel{\mathrm{d}}{⟶}{\chi }_1^2$, hence $\underset{n\to \mathrm{\infty }}{lim}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)=\alpha$. Because $\underset{n\to \mathrm{\infty }}{lim}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)$ achieves its supremum at ${\text{β}}_j=0$ while ${\text{β}}_k\ne 0,k\ne j$, (6(6) $\underset{\text{β}\in {\tilde{\mathrm{H}}}_0}{sup}\underset{n\to \mathrm{\infty }}{lim sup}{\mathrm{P}}_{\text{β}}(\mathcal{L}>{\chi }_{1,\alpha }^2)=\alpha .$(6) ) must hold. Note that, by the definition of ${\tilde{H}}_0$, which has $j\ge 1$, all the indexes j,k mentioned in this paragraph are assumed to be $\ge 1$; in other words, ${\text{β}}_0$ is not involved.

3. Power and sample size calculation

By the results of the previous section, we can use ${\chi }_{1,\alpha }^2$ as the critical value of the LRT. As for the power calculation, although there are infinitely many alternatives that influence the power, it is often reasonable to assume, in practice, that the baseline probability is known. In other words, $h\left({\text{β}}_0\right)$ is known according to (1(1) $p(x_1,\dots ,x_r)=h({\text{β}}_0+{\text{β}}_1{x}_1+\cdots +{\text{β}}_r{x}_r),$(1) ); therefore, ${\text{β}}_0$ is known.

In addition, the minimum probabilistic increase of the largest probability over the other probabilities, δ, is given. In other words, we consider all the alternatives such that (11) $p(1,\dots ,1)\ge p(x_1,\dots ,x_r)+\delta$(11) for all $(x_1,\dots ,x_r)\ne (1,\dots ,1)$. It follows, under (1(1) $p(x_1,\dots ,x_r)=h({\text{β}}_0+{\text{β}}_1{x}_1+\cdots +{\text{β}}_r{x}_r),$(1) ), that all of the ${\text{β}}_j,1\le j\le r$ must be positive, and that (11(11) $p(1,\dots ,1)\ge p(x_1,\dots ,x_r)+\delta$(11) ) is equivalent to (12) $\begin{array}{rl}& h\left({\text{β}}_0+\sum _{j=1}^r{\text{β}}_j\right)\ge h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta ,\\ & 1\le k\le r.\end{array}$(12) The minimum amount of increase of the left sides of (12(12) $\begin{array}{rl}& h\left({\text{β}}_0+\sum _{j=1}^r{\text{β}}_j\right)\ge h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta ,\\ & 1\le k\le r.\end{array}$(12) ) over the right sides takes place when the equalities hold in all of the inequalities, that is, when (13) $\begin{array}{rl}& h\left({\text{β}}_0+\sum _{j=1}^r{\text{β}}_j\right)=h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta ,\\ & 1\le k\le r.\end{array}$(13) This results in r equations, from which we can uniquely determine the alternative. Note that (13(13) $\begin{array}{rl}& h\left({\text{β}}_0+\sum _{j=1}^r{\text{β}}_j\right)=h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta ,\\ & 1\le k\le r.\end{array}$(13) ) is a nonlinear equation system; however, it can be solved conveniently by utilising a Gauss-Seidel type algorithm (e.g., Jiang, 2000). Namely, from (13(13) $\begin{array}{rl}& h\left({\text{β}}_0+\sum _{j=1}^r{\text{β}}_j\right)=h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta ,\\ & 1\le k\le r.\end{array}$(13) ) we have (14) $\begin{array}{rl}{\text{β}}_k& =h^{-1}\left\{h\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right)+\delta \right\}\\ & -{\text{β}}_0-\sum _{1\le j\le r,j\ne k}{\text{β}}_j\\ & =g\left({\text{β}}_0+\sum _{1\le j\le r,j\ne k}{\text{β}}_j\right),\end{array}$(14) where $g\left(x\right)=h^{-1}\left\{h\right(x)+\delta \}-x$ (the inverse of h exists because h is assumed to be strictly increasing). Thus, given the initial values ${\text{β}}_j^{\left(0\right)},1\le j\le r-1$ (e.g., all zero), we have (15) $\begin{array}{rl}{\text{β}}_k^{\left(l\right)}& =g\left\{{\text{β}}_0+\sum _{j=1}^{k-1}{\text{β}}_j^{(l-1)}+\sum _{j=k+1}^r{\text{β}}_j^{\left(l\right)}\right\},\\ & k=r,\dots ,1\end{array}$(15) for $l=1,2,\dots$. The convergence is guaranteed, and fast. We illustrate with an example.

Example 1.1

continued

In the animal clinical trial example, the researchers suggested a baseline probability of 0.3. Using the logistic regression, we have ${\text{β}}_0=\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}\left(0.3\right)=-0.8472979$. Furthermore, the minimum probability increase was set (again by the researchers) as $\delta =0.1$. Recall that r=2 in this case. The Gauss-Seidel algorithm converged within three iterations. Table 1 shows the R outputs of the first five iterations.

Table 1. Convergence of Gauss-Seidel algorithm.

	Initial value: ${\text{β}}_1=0$		Initial value: ${\text{β}}_1=1$
Iteration	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_1$	${\text{β}}_2$
0–1	0.0000000	0.4418328	1.0000000	0.4144315
1–2	0.4054651	0.4070466	0.4066332	0.4069919
2–3	0.4069726	0.4069760	0.4069751	0.4069759
3–4	0.4069759	0.4069759	0.4069759	0.4069759
4–5	0.4069759	0.4069759	0.4069759	0.4069759

Given the target sample size, the power of the LRT at the alternative can be computed by a Monte Carlo method. Consider, for example, the case of logistic regression. Under the alternative $\text{β}=({\text{β}}_0,{\text{β}}_1,\dots ,{\text{β}}_r{)}^{\prime }$, one can simulate data, $N(x_1,\dots ,x_r),x_1,\dots ,x_r\in \{0,1\}$, under the logistic regression. For each simulated data set, r+1 logistic regressions are fit. The first one is under the full model, the next one under the model without $x_1$,…, and the last one under the model without $x_r$. Let $\hat{l},{\hat{l}}_{01},\dots ,{\hat{l}}_{0r}$ denote the maximised log-likelihoods as results of these logistic regressions. We compute $\mathcal{L}=2(\hat{l}-\underset{1\le j\le r}{max}{\hat{l}}_{0j})$ for the simulated data set. This is repeated B times, resulting ${\mathcal{L}}^{\left(b\right)},b=1,\dots ,B$. The power of the LRT is then approximated by $B^{-1}\mathrm{\#}\{1\le b\le B:{\mathcal{L}}^{\left(b\right)}>{\chi }_{1,\alpha }^2\}$.

If, instead, the goal is to determine the sample size so that the LRT has a designated power, say, γ, we can use the following bisection procedure to speed up the search for the minimum sample size. First pick a couple of initial sample sizes, $n_0$ and $n_1$, and compute the power under $n_0$ and $n_1$ using the above procedure. Suppose that the power under $n_0$ is less than γ, and the power under $n_1$ is greater than γ. We then let $n_2=(n_0+n_1)/2$ (take the integer part, if necessary), and compute the power under $n_2$ using the above procedure. If the power under $n_2$ is greater than γ, let $n_3=(n_0+n_2)/2$; otherwise, let $n_3=(n_2+n_1)/2$, and so on. The procedure should converge quickly to a single integer, $n_{\ast }$, so that either $n_{\ast }$ or $n_{\ast }+1$ is the minimum sample size to have the power greater than or equal to γ.

4. Animal clinical trial revisited

Let us go back to Example 1.1 of Section 1. Recall the initial request was to make a power analysis based on the sample size n=30 for detecting a 10% difference between $p(1,1)$ and the rest of the binomial probabilities, and the baseline probability was set as 30%. We considered a logistic regression with r=2 for this case. Here $p_0=0.3$ and $\delta =0.1$. The alternative was computed by the Gauss-Seidel algorithm [see Example 1.1 (continued) in Section 3] as ${\text{β}}_0=-0.8472979$ and ${\text{β}}_1={\text{β}}_2=0.4069759$. The corresponding probabilities are $p(0,0)=0.3$, $p(0,1)=p(1,0)=0.3916643$, and $p(1,1)=0.4916643$. As we can see, the minimum difference between $p(1,1)$ and the rest of the p's is 0.1. For this alternative, the power of the LRT at 5% level of significance was computed, using the Monte-Carlo method described in Section 3 with B=1000, as approximately 88%.

As the 80% power was considered satisfactory by the researchers, it appeared that the sample size might be reduced a little. However, when the same procedure was applied to n=20, the power was computed as approximately 78%. In communicating with the lead researcher, the researcher suggested that he would rather sacrifice a little regarding the minimum probabilistic difference in exchange for a reduced sample size (i.e., n=20). Thus, the new δ was set as 0.15. The new alternative was computed by the Gauss-Seidel algorithm (which, again, converged in three iterations) as ${\text{β}}_0=-0.8472979$ and ${\text{β}}_1={\text{β}}_2=0.6051085$. The corresponding probabilities are $p(0,0)=0.3$, $p(0,1)=p(1,0)=0.4397469$, and $p(1,1)=0.5897469$. As we can see, the minimum probabilistic increase of $p(1,1)$ over the rest of the p's is 0.15. For the new alternative, the power of the LRT at 5% level of significance was computed, again using the Monte-Carlo method with B=1000, as approximately 86%. The researcher was satisfied with the result.

Acknowledgments

The authors are grateful to Dr. Markus Grompe of the Doernbecher Children's Hospital of the Oregon Health & Science University for presenting the problem from their research, and for information and helpful discussions. The authors also wish to thank a reviewer for helpful comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The authors' research is partially supported by the National Institutes of Health (NIH) grant R01-GM085205A1. In addition, Thuan Nguyen's research is partially supported by the National Science Foundation (NSF) grant SES-1118469; Jiming Jiang's research is partially supported by the National Science Foundation (NSF) grant SES-1121794.

Notes on contributors

Thuan Nguyen

Thuan Nguyen is Associate Professor, Department of Public Health and Preventive Medicine, Oregon Health and Science University, USA.

Jiming Jiang

Jiming Jiang is Professor, Department of Statistics, University of California, Davis, USA.

Appendix

A.1. An inequality

Consider the function $g\left(p\right)=p^{p_0}(1-p{)}^{1-p_0}$ and $h\left(p\right)=\log \left\{g\right(p\left)\right\},0<p<1$. We have $h^{\prime }\left(p\right)=\left\{p\right(1-p\left){\}}^{-1}\right(p_0-p)$. Thus, $h^{\prime }\left(p\right)>0$, $=0$, or $<0$ depending on $p<p_0$, $p=p_0$, or $p>p_0$. It follows that $h(\cdot )$, hence $g(\cdot )$, has a unique maximum at $p=p_0$, and $g\left(p\right)<g\left(p_0\right)$ for any $p\ne p_0$. Thus, for any $0<p_1<1,p_1\ne p_0$, we have ${\left(\frac{p_0}{p_1}\right)}^{p_0}{\left(\frac{1-p_0}{1-p_1}\right)}^{1-p_0}=\frac{g\left(p_0\right)}{g\left(p_1\right)}>1.$

A.2. Some derivation in Section 2

We show that $\mathcal{L}\stackrel{\mathrm{d}}{⟶}\eta -{\eta }_1\vee {\eta }_2$ as $n\to \mathrm{\infty }$, where the η's are defined in Section 2, if $({\text{β}}_0,0,0)$ is the true parameter vector. With the notation introduced in Section 2, we have, by the Taylor expansion, $\begin{array}{rl}& l({\text{β}}_0,0,0)=l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)\\ & +\frac{1}{2}{\left(\begin{array}{c}{\text{β}}_0-{\widehat{\text{β}}}_0^{\left[2\right]}\\ -{\widehat{\text{β}}}_1^{\left[2\right]}\end{array}\right)}^{\prime }\\ & \times {\left.\left(\begin{array}{cc}{\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0^2& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0\mathrm{\partial }{\text{β}}_1\\ {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1\mathrm{\partial }{\text{β}}_0& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1^2\end{array}\right)\right|}_{({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)}\\ & \times \left(\begin{array}{c}{\text{β}}_0-{\widehat{\text{β}}}_0^{\left[2\right]}\\ -{\widehat{\text{β}}}_1^{\left[2\right]}\end{array}\right)\end{array}$$+o_{\mathrm{P}}\left(1\right)$, implying $\begin{array}{rl}& l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)=l({\text{β}}_0,0,0)\\ & -\frac{1}{2}{\left(\begin{array}{c}{\widehat{\text{β}}}_0^{\left[2\right]}-{\text{β}}_0\\ {\widehat{\text{β}}}_1^{\left[2\right]}\end{array}\right)}^{\prime }\\ & \mathrm{E}\left\{{\left.\left(\begin{array}{cc}{\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0^2& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0\mathrm{\partial }{\text{β}}_1\\ {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1\mathrm{\partial }{\text{β}}_0& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1^2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}\right\}\left(\begin{array}{c}{\widehat{\text{β}}}_0^{\left[2\right]}-{\text{β}}_0\\ {\widehat{\text{β}}}_1^{\left[2\right]}\end{array}\right)\end{array}$$+o_{\mathrm{P}}\left(1\right)$. Also, by the standard asymptotic expansion (e.g., Jiang, 2010, Ch. 4), we have $\begin{array}{rl}& \left(\begin{array}{c}{\widehat{\text{β}}}_0^{\left[2\right]}-{\text{β}}_0\\ {\widehat{\text{β}}}_1^{\left[2\right]}\end{array}\right)={\left[\mathrm{E}\left\{{\left.\left(\begin{array}{cc}{\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0^2& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0\mathrm{\partial }{\text{β}}_1\\ {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1\mathrm{\partial }{\text{β}}_0& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1^2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}\right\}\right]}^{-1}\\ & \times {\left.\left(\begin{array}{c}\mathrm{\partial }l/\mathrm{\partial }{\text{β}}_0\\ \mathrm{\partial }l/\mathrm{\partial }{\text{β}}_1\end{array}\right)\right|}_{({\text{β}}_0,0,0)}+O_{\mathrm{P}}\left(n^{-1}\right).\end{array}$ It follows that $\begin{array}{rl}& l({\widehat{\text{β}}}_0^{\left[2\right]},{\widehat{\text{β}}}_1^{\left[2\right]},0)=l({\text{β}}_0,0,0)\\ & -\frac{1}{2}{\left.\left(\begin{array}{c}\mathrm{\partial }l/\mathrm{\partial }{\text{β}}_0\\ \mathrm{\partial }l/\mathrm{\partial }{\text{β}}_1\end{array}\right)\right|}_{({\text{β}}_0,0,0)}^{\prime }\\ & {\left[\mathrm{E}\left\{{\left.\left(\begin{array}{cc}{\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0^2& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0\mathrm{\partial }{\text{β}}_1\\ {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1\mathrm{\partial }{\text{β}}_0& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_1^2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}\right\}\right]}^{-1}\\ & \times {\left.\left(\begin{array}{c}\mathrm{\partial }l/\mathrm{\partial }{\text{β}}_0\\ \mathrm{\partial }l/\mathrm{\partial }{\text{β}}_1\end{array}\right)\right|}_{({\text{β}}_0,0,0)}+o_{\mathrm{P}}\left(1\right).\end{array}$ Similarly, we have $\begin{array}{rl}& l({\widehat{\text{β}}}_0^{\left[1\right]},0,{\widehat{\text{β}}}_2^{\left[1\right]})=l({\text{β}}_0,0,0)\\ & -\frac{1}{2}{\left.\left(\begin{array}{c}\mathrm{\partial }l/\mathrm{\partial }{\text{β}}_0\\ \mathrm{\partial }l/\mathrm{\partial }{\text{β}}_2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}^{\prime }\\ & {\left[\mathrm{E}\left\{{\left.\left(\begin{array}{cc}{\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0^2& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_0\mathrm{\partial }{\text{β}}_2\\ {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_2\mathrm{\partial }{\text{β}}_0& {\mathrm{\partial }}^{2}l/\mathrm{\partial }{\text{β}}_2^2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}\right\}\right]}^{-1}\\ & \times {\left.\left(\begin{array}{c}\mathrm{\partial }l/\mathrm{\partial }{\text{β}}_0\\ \mathrm{\partial }l/\mathrm{\partial }{\text{β}}_2\end{array}\right)\right|}_{({\text{β}}_0,0,0)}+o_{\mathrm{P}}\left(1\right),\end{array}$

and $\begin{array}{rl}& l({\widehat{\text{β}}}_0,{\widehat{\text{β}}}_1,{\widehat{\text{β}}}_2)=l({\text{β}}_0,0,0)\\ & -\frac{1}{2}{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }\text{β}}\right|}_{({\text{β}}_0,0,0)}^{\prime }{\left\{\mathrm{E}\left({\left.\frac{{\mathrm{\partial }}^{2}l}{\mathrm{\partial }\text{β}\mathrm{\partial }{\text{β}}^{\prime }}\right|}_{({\text{β}}_0,0,0)}\right)\right\}}^{-1}{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }\text{β}}\right|}_{({\text{β}}_0,0,0)}\\ & +o_{\mathrm{P}}\left(1\right).\end{array}$ Furthermore, we have (16) $\mathrm{E}\left({\left.\frac{{\mathrm{\partial }}^{2}l}{\mathrm{\partial }\text{β}\mathrm{\partial }{\text{β}}^{\prime }}\right|}_{({\text{β}}_0,0,0)}\right)=-n{h}^{\prime }\left({\text{β}}_0\right)\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right),$(16) and, by the standard asymptotic result, (17) $\frac{1}{\sqrt{n}}{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }\text{β}}\right|}_{({\text{β}}_0,0,0)}\stackrel{\mathrm{d}}{⟶}N\left[\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right),h^{\prime }\left({\text{β}}_0\right)\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right)\right],$(17) Let ${\xi }_n=({\xi }_{n,0},{\xi }_{n,1},{\xi }_{n,2}{)}^{\prime }$ denote the left side of (A2(17) $\frac{1}{\sqrt{n}}{\left.\frac{\mathrm{\partial }l}{\mathrm{\partial }\text{β}}\right|}_{({\text{β}}_0,0,0)}\stackrel{\mathrm{d}}{⟶}N\left[\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right),h^{\prime }\left({\text{β}}_0\right)\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right)\right],$(17) ) divided by $\sqrt{h^{\prime }\left({\text{β}}_0\right)}$, and A denote the $3\times 3$ matrix in (A1(16) $\mathrm{E}\left({\left.\frac{{\mathrm{\partial }}^{2}l}{\mathrm{\partial }\text{β}\mathrm{\partial }{\text{β}}^{\prime }}\right|}_{({\text{β}}_0,0,0)}\right)=-n{h}^{\prime }\left({\text{β}}_0\right)\left(\begin{array}{ccc}4& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right),$(16) ). For $a=(a_0,a_1,a_2{)}^{\prime }$ and $A=(a_{st}{)}_{s,t=0,1,2}$, let $a[0,j]$ and $A[0,j]$ denote the subvector $(a_0,a_j{)}^{\prime }$ and submatrix $(a_{st}{)}_{s,t=0,j}$, respectively, j=1,2. Then, by the above expressions, it is easy to show that $\mathcal{L}={\xi }_n^{\prime }{A}^{-1}{\xi }_n-\underset{j=1,2}{max}\left\{{\xi }_n\right[0,j{]}^{\prime }A[0,j{]}^{-1}{\xi }_n[0,j\left]\right\}+o_{\mathrm{P}}\left(1\right).$ Because ${\xi }_n\stackrel{\mathrm{d}}{⟶}\xi =({\xi }_0,{\xi }_1,{\xi }_2{)}^{\prime }\sim N(0,A)$, as $n\to \mathrm{\infty }$, by the continuous mapping theorem (e.g., Jiang, 2010, p.30), we have $\mathcal{L}\stackrel{\mathrm{d}}{⟶}\eta -{\eta }_1\vee {\eta }_2$, where $\eta ={\xi }^{\prime }A\xi$ and ${\eta }_j=\xi [0,j{]}^{\prime }A[0,j{]}^{-1}\xi [0,j],j=1,2$.