Empirical Likelihood Inference of Variance Components in Linear Mixed-Effects Models

Longitudinal and clustered data commonly arise from observational studies or clinical trials, where subjects are measured repeatedly over time or within a cluster. The repeated measures within a subject or a cluster are often correlated. To analyze such data, linear mixed-effects models that incorporate both fixed and random effects are widely used. Many statistical methods have been developed for such linear mixed-effects models, especially methods for inference of the fixed effects. However, inference on the variance components is less studied and often requires strong distributional assumptions on the random effects and the error terms. When the underlying distributions are known, classical inference methods, including the likelihood ratio tests and the score tests, can be applied. However, these parametric methods are often restrictive and not robust if the model assumptions are violated.

Empirical likelihood (el) method, as an alternative to parametric likelihood-based methods, was first proposed by Owen, (1988) and has been applied to many statistical inference problems. Without a prespecified distributional assumption on the data, el methods incorporate side information through constraints or prior distributions and have favorable statistical properties, including but not limited to Bartlett correctability, transformation invariance, better coverage accuracy of the corresponding confidence internals and greater power. A comprehensive introduction to el methods can be found in Owen, (2001). el methods have been applied to inferences of mixture models (Zou et al.,, 2002) and censored survival data (Chang and McKeague,, 2016), and have also been considered for longitudinal data modeling. For example, You et al., (2006) proposed a block el method for inference of the regression parameters assuming a working independence covariance, and Xue and Zhu, (2007) considered a semiparametric regression model, where the repeated within-subject measures are summarized as a function over time in order to address the dependence issue. Wang et al., (2010) proposed a generalized el method that takes into account the within-subject correlations. Li and Pan, (2013) defined an empirical likelihood ratio (elr) test by utilizing the extended score from quadratic inference functions for longitudinal data, which does not involve direct estimation of the correlation parameters.

The el methods mentioned above only focus on the inference of fixed effects in linear mixed effect models. In this paper, we consider a general setting of linear mixed-effects models and develop el methods for the inference of the variance components. Specifically, suppose there are $n$ subjects and denote by $n_{i}$ the number of repeated measures for the $i$th subject. For the $i$th subject, we observe a response vector $y_{i}\in R^{n_{i}}$, an $n_{i}\times p$ design matrix $X_{i}$ for the fixed effects $\beta^{*}\in R^{p}$, and $d$ $n_{i}\times n_{i}$ semi-positive design matrices $\Phi_{iq}~{}(q=1,\cdots,d)$ for the variance components $\theta^{*}\in(R_{+}\cup\{0\})^{d}$. The general linear mixed-effects model can be written as

$$y_{i}=X_{i}\beta^{*}+r_{i},\quad i=1,\cdots,n,$$

(1)

where $r_{i}\in R^{n_{i}}$ is a zero-mean random variable with variance-covariance $H_{i}(\theta^{*})$. We assume that $H_{i}(\theta^{*})$ has a linear structure,

$$H_{i}(\theta^{*})=\sum_{q=1}^{d}\theta^{*}_{q}\Phi_{iq},\quad\theta^{*}=(\theta^{*}_{1},\cdots,\theta^{*}_{d})^{T}=(\theta^{*}_{1},\theta^{*T}_{(1)})^{T},$$

where $\theta^{*}=(\theta^{*}_{1},\cdots,\theta^{*}_{d})^{T}$ is the vector of the variance components. We emphasize that this general setting does not require any assumptions on the distributions of the data or the distributions of the random effects.

In many real applications, we are interested in making statistical inference on the variance components $\theta^{*}$ in model (1). For example, in the study of heritability based on twin data, each monozygotic twin or dizygotic twin is treated as one cluster, and the linear variance structure can be constructed based on the twin type (see details in Section 6). In the heritability analysis, a key question is whether there exists an genetic effect, which motivates us to study the inference of one of the variance components, say, $\theta^{*}_{1}$. We propose to develop an el based inference method for $\theta^{*}_{1}$ without any assumptions on the random components. The method can effectively account for the nuisance parameters, including the unknown fixed effects $\beta^{*}$ and the variance components $\theta^{*T}_{(1)}$. The key difficulty when compared to the el inference of the fixed effects is to deal with the boundary value problem when $\theta^{*}_{1}=0$ in local testing problem $H_{0}:\theta^{*}_{1}=\theta_{1}^{0}$. To solve the issues, we propose a new empirical likelihood ratio test by utilizing an unbiased estimator of $\beta^{*}$ under very mild conditions, and prove that the asymptotic distribution of the test statistic is a mixture of $\chi^{2}$ distribution.

Motivated by heritability analysis of daily activity distribution as measured by wearable device such as actigraphy, we also consider the setting when linear mixed-effects models are fitted to a sequence of dependent outcomes. The wearable device data have been increasingly collected for continuous activity monitoring in large observational or experimental studies (Burton et al.,, 2013; Krane-Gartiser et al.,, 2014). In typical wearable activity tracking data, the activity is measured at one-minute resolution over several days for a given subject. Such wearable device data with repeated measures enable us to account for day-to-day variability of the activity. Instead of focusing on the activity counts at any minute of the day, daily activity distribution or the amount of time with the activity count above a given threshold provides a biologically meaningful measure of the activity traits. When the activity counts are summarized as distributions, we can consider the activity quantile profile as a phenotype measure. In analysis of such wearable device data, we fit a linear mixed-effects model for each of the activity level or quantile $y_{i}(t)$ at $t$. Denote by $\theta^{*}(t)$ the variance components for activity profile at level $t$. We are then interested in testing the global null $H_{0}:\theta^{*}_{1}(t)\equiv\theta_{1}^{0},~{}t\in[t_{1},t_{2}]$. We develop a max-type statistic for this global testing problem. Since the numerator of the proposed empirical likelihood ratio (elr) tests can be rewritten as the sum of approximately independent random variables over different subjects, a random perturbation method is developed to approximate the $p$-values of the proposed global test.

We first introduce some notation. Denote by $(A)_{-1}$ the submatrix of $A$ without the first column of $A$. For two vectors or matrices $A$ and $B$ of compatible dimension, define the inner product $\langle A,B\rangle=\text{tr}(A^{T}B)$. For a matrix $D_{m\times n}=(D_{1},\cdots,D_{n})$, where $D_{i}$ is the $i$th column of $D$, the vectorized $D$ is defined by $(D_{1}^{T},\cdots,D_{n}^{T})^{T}$. Let $E(x)$ and $\text{var}(x)$ be the expectation and variance of a random vector $x$, and let $\text{cov}(x,y)$ be the covariance of random vectors $x$ and $y$. When $x$ is a random matrix, $E(x)$ and $\text{var}(x)$ represent the expectation and variance of the vectorized $x$. When $x$ and $y$ are random matrices, $\text{cov}(x,y)$ denotes the covariance of the vectorized $x$ and vectorized $y$. We use $a=O(b)$ to denote that $a$ and $b$ are of the same order, and $a=o(b)$ to denote that $a$ is of a smaller order than $b$. We use $x=O_{p}(y)$ to denote that $x$ and $y$ are of the same order in probability, and $x=o_{p}(y)$ to denote that $x$ is of a smaller order than y in probability.

Statistical tests for the fixed effects in the linear mixed-effects model (1), $H_{0}:\beta^{*}=\beta_{0}$, has been well studied. We first briefly review the subject-wise el method proposed in Wang et al., (2010), where the covariance structure for each subject is considered.

Let $\hat{H}_{in}$ be an estimator of $H_{i}$, and assume that $\hat{H}_{in}$ converges to some $H_{i}^{*}$ in probability uniformly over all $i=1,\cdots,n$. One such a nonparametric sample covariance matrix $\hat{H}_{in}$ can be obtained using a simple two-step procedure, including estimating the residuals $\hat{r}_{i}=y_{i}-X_{i}\hat{\beta}$, where $\hat{\beta}$ is the least-squares estimator using working independence correlation matrices, and solving the constrained optimization problem $\min_{\theta\geq 0}\sum_{i=1}^{n}\|H_{i}(\theta)-\hat{r}_{i}\hat{r}_{i}^{T}\|_{F}^{2}$. Let

$$\phi_{i}(\beta)=X_{i}^{T}\hat{H}_{in}^{-1}(y_{i}-X_{i}\beta),$$

which satisfies $E\{\phi_{i}(\beta)\}=0$ when $\beta$ is the true value. Denote by $p_{i}$ the point mass at the $i$th subject. The nonparametric empirical likelihood is defined as

$$L_{0}(\beta)=\sup_{p_{i}}\Big{\{}\prod_{i=1}^{n}p_{i}:p_{i}\geq 0,\sum_{i=1}^{n}p_{i}=1,\sum_{i=1}^{n}p_{i}\phi_{i}(\beta)=0\Big{\}}.$$

Since it can be proved that $\max_{\beta}L_{0}(\beta)=1/n^{n}$ (Owen,, 2001), Wang et al., (2010) proposed the elr statistic

$$\textsc{elr}_{0}(\beta_{0})=\frac{L_{0}(\beta_{0})}{\max_{\beta}L_{0}(\beta)}=n^{n}L_{0}(\beta_{0}).$$

To obtain the asymptotic distribution of the elr statistic, the following regularity conditions are needed.

Conditions 1–3 are common conditions for the empirical likelihood methods (Owen,, 1991). Condition 4 assumes mild constraints on $\hat{H}_{in}^{-1}$ to ensure that the difference between the statistic $\textsc{elr}_{0}(\beta_{0})$ defined with $\hat{H}_{in}$ and the one using $H_{i}^{*}$ vanishes as $n\rightarrow\infty$. Under these regularity conditions, the following theorem provides the asymptotic distribution of the elr test $\textsc{elr}_{0}(\beta_{0})$ (Wang et al.,, 2010) under the null.

The asymptotic result only requires that the $\hat{H}_{in}$ converge uniformly to some $H_{i}^{*}$, which may not be the true $H_{i}$ (Wang et al.,, 2010). When the correlation structure is correctly specified, the estimator $\hat{H}_{in}$ is a consistent estimator of $H_{i}^{*}=H_{i}$. The statistic defined with the true $H_{i}$ is asymptotically locally most powerful among all the choices of the weight matrices.

We consider the local test $H_{0}:\theta^{*}_{1}=\theta_{1}^{0}$ in the framework of the empirical likelihood, including the null $H_{0}:\theta_{1}^{*}=0$, which is of the most interest. We define $r_{i}=y_{i}-X_{i}\beta^{*}$ and $R_{i}=r_{i}r_{i}^{T}$. Since $E(r_{i})=0$ and $\text{var}(r_{i})=H_{i}(\theta^{*})$, we have

$$R_{i}=H_{i}(\theta^{*})+\delta_{i}=\sum_{q=1}^{d}\theta^{*}_{q}\Phi_{iq}+\delta_{i},$$

where $E({\delta_{i}})=0$ and $\text{var}({\delta_{i}})$ exists. Since $\beta^{*}$ is unknown, we first need an estimator of $\beta^{*}$, denoted by $\hat{\beta}$. One simple choice is the least-squares estimator using the all data. Specifically, we stack the data from all subjects by denoting $X=(X_{1}^{T},\cdots,X_{n}^{T})^{T}$, $y=(y_{1}^{T},\cdots,y_{n}^{T})^{T}$, and $r=(r_{1}^{T},\cdots,r_{n}^{T})^{T}$. Model (1) can be rewritten as

$$y=X\beta^{*}+r.$$

Then the least-squares estimator is $\hat{\beta}=(X^{T}X)^{-1}X^{T}y$. For $i=1,\cdots,n$, let $\hat{r}_{i}=y_{i}-X_{i}\hat{\beta}=r_{i}+X_{i}(\beta^{*}-\hat{\beta}).$ We have

	$\displaystyle\hat{R}_{i}$	$\displaystyle=\hat{r}_{i}\hat{r}_{i}^{T}=r_{i}r_{i}^{T}+r_{i}(\beta^{*}-\hat{\beta})^{T}X_{i}^{T}+X_{i}(\beta^{*}-\hat{\beta})r_{i}^{T}+X_{i}(\beta^{*}-\hat{\beta})(\beta^{*}-\hat{\beta})^{T}X_{i}^{T}$
		$\displaystyle=R_{i}+\hat{\epsilon}_{i}=H_{i}(\theta^{*})+\delta_{i}+\hat{\epsilon}_{i},$

where $\hat{\epsilon}_{i}=r_{i}(\beta^{*}-\hat{\beta})^{T}X_{i}^{T}+X_{i}(\beta^{*}-\hat{\beta})r_{i}^{T}+X_{i}(\beta^{*}-\hat{\beta})(\beta^{*}-\hat{\beta})^{T}X_{i}^{T}$.

To control the rates of $E({\hat{\epsilon}_{i}})$, $\text{cov}({r_{i}r_{i}^{T}},{\hat{\epsilon}_{j}})$, and $\text{cov}({\hat{\epsilon}_{i}},{\hat{\epsilon}_{j}})$, we need the following condition, which is also commonly used for empirical likelihood methods.

Under Condition 5, we see that the least-squares estimator $\hat{\beta}$ is good enough.

Let $\Xi=(\Xi_{kl})_{d\times d}$ with $\Xi_{kl}=\sum_{i=1}^{n}\text{tr}(\Phi_{ik}\Phi_{il})$, and $\hat{\Upsilon}=(\hat{\Upsilon}_{1},\cdots,\hat{\Upsilon}_{d})^{T}$ with $\hat{\Upsilon}_{k}=\sum_{i=1}^{n}\text{tr}(\Phi_{ik}\hat{R}_{i})$. We define

$$\hat{Z}_{i}(\theta_{1})=\text{tr}\Bigg{\{}\Phi_{i1}\Bigg{(}\hat{R}_{i}-\Phi_{i1}\theta_{1}-\sum_{q=2}^{d}\hat{\theta}_{q}\Phi_{iq}\Bigg{)}\Bigg{\}},~{}i=1,\cdots,n,$$

where

$$\hat{\theta}_{(1)}=(\hat{\theta}_{2},\cdots,\hat{\theta}_{q})^{T}=(\Xi^{-1})_{-1}^{T}\hat{\Upsilon}.$$

(2)

Since Proposition 1 implies $E\hat{Z}_{i}(\theta_{1})=O(n^{-1})$ if $\theta_{1}$ is the true value (see (A18) in the appendix), we define the nonparametric likelihood as

$$L(\theta_{1})=\max_{p_{i}}\Bigg{\{}\prod_{i=1}^{n}p_{i}|p_{i}\geq 0,\sum_{i=1}^{n}p_{i}=1,\sum_{i=1}^{n}p_{i}\hat{Z}_{i}(\theta_{1})=0\Bigg{\}}$$

and the corresponding elr statistic as

$$\textsc{elr}(\theta_{1}^{0})=\frac{L(\theta_{1}^{0})}{\max_{\theta_{1}\geq 0}L(\theta_{1})}.$$

(3)

If the true value $\theta_{1}^{*}=0$ (i.e., the null hypothesis under the case $\theta_{1}^{0}=0$), the denominator in (3) would not be $1/n^{n}$ as usual owing to the boundary value issue, and thus the existing results are inapplicable. To derive the asymptotic distribution of the proposed test $\textsc{elr}(\theta_{1}^{0})$, we assume the following condition similar to Condition 1.

Under Conditions 5 and 6, we have the following theorem on the asymptotic distribution of the elr test under the null.

Although the elr statistic $\hat{c}_{n}(\theta_{1}^{0})\left(-2\log\textsc{elr}(\theta_{1}^{0})\right)$ in Theorem 2 involves optimizations in the numerator and denominator, the following lemma shows that the statistic has an asymptotically equivalent expression that can be used to calculate the statistic efficiently.

We next provide an estimator of the asymptotic variance of $n^{-1/2}\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0})$. We rewrite $\Xi$ as

$$\Xi=\left(\begin{smallmatrix}E_{11}&E_{12}\\ E_{21}&E_{22}\end{smallmatrix}\right)$$

with $E_{11}$ being a scalar. Let $F=E_{22}^{-1}E_{21}=(F_{1},\cdots,F_{d-1})^{T}$ and $\alpha=1-E_{12}F/E_{11}\in(0,1]$. It can be verified that

$$\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0})=\sum_{i=1}^{n}\hat{D}_{i}(\theta_{1}^{0})=\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0}),$$

(4)

where

	$\displaystyle\hat{D}_{i}(\theta_{1}^{0})$	$\displaystyle=\alpha^{-1}\langle\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1},\hat{R}_{i}-\theta_{1}^{0}\Phi_{i1}\rangle,$
	$\displaystyle\hat{M}_{i}(\theta_{1}^{0})$	$\displaystyle=\alpha^{-1}\langle\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1},\hat{R}_{i}-H_{i}((\theta_{1}^{0},\hat{\theta}_{(1)}^{T})^{T})\rangle.$

In addition, for $i\neq j$,

$\displaystyle\text{cov}({\hat{R}_{i}},{\hat{R}_{j}})=$

$\displaystyle\text{cov}({\delta_{i}}+{\hat{\epsilon}_{i}},{\delta_{j}}+{\hat{\epsilon}_{j}})=\text{cov}({\delta_{i}},{\hat{\epsilon}_{j}})+\text{cov}({\hat{\epsilon}_{i}},{\delta_{j}})+\text{cov}({\hat{\epsilon}_{i}},{\hat{\epsilon}_{j}})=O(n^{-2})$

based on proposition 1. Therefore, ${\hat{D}_{i}(\theta_{1}^{0})}~{}(i=1,\cdots,n)$ are asymptotically independent with expectation $E(\hat{D}_{i}(\theta_{1}^{0}))=\alpha^{-1}\langle\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1},\sum_{q=2}^{d}\theta_{q}^{*}\Phi_{iq}\rangle+O(n^{-1})$, while the expectations of $\hat{Z}_{i}(\theta_{1}^{0})$ and $\hat{M}_{i}(\theta_{1}^{0})$ are $O(n^{-1})$ (see (A18) in the appendix). We have

	$\displaystyle\hat{D}_{i}(\theta_{1}^{0})-E(\hat{D}_{i}(\theta_{1}^{0}))$	$\displaystyle=\alpha^{-1}\langle\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1},\hat{R}_{i}-H_{i}((\theta_{1}^{0},(\theta_{(1)}^{*})^{T})^{T})\rangle+O(n^{-1})$
		$\displaystyle=\hat{M}_{i}(\theta_{1}^{0})+o_{p}(1).$

Therefore,

	$\displaystyle\text{var}\big{\{}n^{-1/2}\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0})\big{\}}=$	$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\{\hat{D}_{i}(\theta_{1}^{0})-E(\hat{D}_{i}(\theta_{1}^{0}))\}^{2}+o_{p}(1)$
	$\displaystyle=$	$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0})^{2}+o_{p}(1),$		(5)

which leads a consistent estimator of the variance of $n^{-1/2}\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0})$ as

$$\hat{\nu}_{1n}^{2}(\theta_{1}^{0})=n^{-1}\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0})^{2}.$$

In some applications, we are interested in testing whether the variance components are all zero over a set of possibly correlated outcomes. One example of such applications is to test the variance components for the activity distribution based on wearable device data where we are interested in testing the variance component at each of the quantiles $t$ of the activity distribution. Extending model (1), we assume the following outcome model at level $t$,

$$y_{i}(t)=X_{i}\beta^{*}(t)+r_{i}(t),\quad i=1,\cdots,n,$$

(6)

where $r_{i}(t)\in R^{n_{i}}$ is a zero-mean random variable with variance $H_{i}(\theta^{*}(t))$. We assume that $H_{i}(\theta^{*}(t))$ has the same linear structure for each $t$,

$$H_{i}(\theta^{*}(t))=\sum_{q=1}^{d}\theta^{*}_{q}(t)\Phi_{iq},\quad\theta^{*}(t)=\{\theta^{*}_{1}(t),\cdots,\theta^{*}_{d}(t)\}^{T}=\{\theta^{*}_{1}(t),\theta^{*T}_{(1)}(t)\}^{T}.$$

We are interested in testing the null $H_{0}:\theta^{*}_{1}(t)\equiv\theta_{1}^{0},~{}t\in[t_{1},t_{2}]$, where $[t_{1},t_{2}]$ is a pre-defined interval. We propose the following maximally selected empirical likelihood ratio statistic (gelr),

$$\Gamma=\sup_{t\in[t_{1},t_{2}]}\hat{c}_{n}(\theta_{1}^{0},t)\left\{-2\log\textsc{elr}(\theta_{1}^{0},t)\right\},$$

(7)

where $\hat{c}_{n}(\theta_{1}^{0},t)\left\{-2\log\textsc{elr}(\theta_{1}^{0},t)\right\}$ is the elrstatistic for the outcome at $t$. It can be shown that $\Gamma=\sup_{t\in[t_{1},t_{2}]}S(t)+o_{p}(1),$ with

$$S(t)=\begin{cases}\frac{\{n^{-1/2}\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0},t)\}^{2}}{\hat{\nu}_{1n}^{2}(\theta_{1}^{0},t)},&\text{ if }\theta_{1}^{0}>0,\\ \frac{\{n^{-1/2}\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0},t)\}^{2}}{\hat{\nu}_{1n}^{2}(\theta_{1}^{0},t)}I\{\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0},t)\geq 0\},&\text{ if }\theta_{1}^{0}=0,\end{cases}$$

where

	$\displaystyle\hat{Z}_{i}(\theta_{1}^{0},t)$	$\displaystyle=\text{tr}\Bigg{\{}\Phi_{i1}\Bigg{(}\hat{R}_{i}(t)-\Phi_{i1}\theta_{1}^{0}-\sum_{q=2}^{d}\hat{\theta}_{q}(t)\Phi_{iq}\Bigg{)}\Bigg{\}},$
	$\displaystyle\hat{\nu}_{1n}^{2}(\theta_{1}^{0},t)$	$\displaystyle=n^{-1}\alpha^{-2}\sum_{i=1}^{n}\Bigg{\langle}\hat{R}_{i}(t)-H_{i}((\theta_{1}^{0},\hat{\theta}_{(1)}(t)^{T})^{T}),\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1}\Bigg{\rangle}^{2}.$

Assessment of the statistical significance of the statistic $\Gamma$ defined in (7) is challenging because of the dependence of $\hat{Z}_{i}(\theta_{1}^{0},t)$. We propose a simple way of evaluating its significance by perturbing the el statistic. Specifically, we apply (4) to rewrite $\sum_{i=1}^{n}\hat{Z}_{i}(\theta_{1}^{0},t)$ as $\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0},t)$, where $\hat{M}_{i}(\theta_{1}^{0},t)=\alpha^{-1}\left\langle\Phi_{i1}-\sum_{q=1}^{d-1}F_{q}\Phi_{iq+1},\hat{R}_{i}(t)-H_{i}((\theta_{1}^{0},\hat{\theta}_{(1)}^{T}(t))^{T})\right\rangle$. We can generate the null distribution of $\Gamma$ by perturbing the test statistic $\Gamma^{(g)}$. Specifically, for each $g~{}(g=1,\cdots,G)$, we generate $\xi_{i}^{(g)}$ from i.i.d. standard normal distribution and define

$$S^{(g)}(t)=\begin{cases}\frac{\{n^{-1/2}\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0},t)\xi_{i}^{(g)}\}^{2}}{\hat{\nu}_{1n}^{2}(\theta_{1}^{0},t)},&\text{ if }\theta_{1}^{0}>0,\\ \frac{\{n^{-1/2}\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0},t)\xi_{i}^{(g)}\}^{2}}{\hat{\nu}_{1n}^{2}(\theta_{1}^{0},t)}I\{\sum_{i=1}^{n}\hat{M}_{i}(\theta_{1}^{0},t)\xi_{i}^{(g)}\geq 0\},&\text{ if }\theta_{1}^{0}=0.\end{cases}$$

Define the corresponding perturbed test statistic as $\Gamma^{(g)}=\sup_{t\in[t_{1},t_{2}]}S^{(g)}(t)$. The following Proportion 2 shows that the perturbed test statistics have the same distribution as the original test statistic under the null. Therefore, the $p$-value of $\Gamma$ can be approximated by $\sum_{g=1}^{G}I(\Gamma^{(g)}>\Gamma)/G$.