Sequential estimation for GEE with adaptive variables and subject selection

TL;DR

This paper introduces a sequential estimation method for GEE that adaptively selects subjects and variables, improving efficiency and variable identification in correlated data analysis.

Contribution

It proposes a novel GEE-based sequential estimation approach with adaptive subject recruitment and variable selection, including an adaptive shrinkage feature.

Findings

Method effectively identifies relevant variables during estimation.

Demonstrates improved estimation accuracy on simulated and real data.

Provides confidence sets with pre-specified precision.

Abstract

Modeling correlated or highly stratified multiple-response data becomes a common data analysis task due to modern data monitoring facilities and methods. Generalized estimating equations (GEE) is one of the popular statistical methods for analyzing this kind of data. In this paper, we present a sequential estimation procedure for obtaining GEE-based estimates. In addition to the conventional random sampling, the proposed method features adaptive subject recruiting and variable selection. Moreover, we equip our method with an adaptive shrinkage property so that it can decide the effective variables during the estimation procedure and build a confidence set with a pre-specified precision for the corresponding parameters. In addition to the statistical properties of the proposed procedure, we assess our method using both simulated data and real data sets.

Tables

Table 1. Table 1: Results of sequential sampling method, performance in variable identification, and estimation of number of nonzero components for identity link function with ar (1) correlation when p 0 : p k = 4 : 20 : subscript 𝑝 0 subscript 𝑝 𝑘 4 : 20 p_{0}:p_{k}=4:20 .

$\alpha$	$d$	Method	$N^{\ast }$	$C.P.$	$\kappa$	$Nu{m}_c$	$Nu{m}_{ic}$
0.3	0.4	Oracle	34.498 (4.931)	0.887	1.079	0.000	0.000
		ASE-D	36.452 (6.278)	0.669	1.126	19.714	0.000
		ASE-R	40.947 (11.430)	0.604	1.133	19.552	0.000
		GEE	156.375 (12.067)	0.690	1.009	0.033	0.000
	0.3	Oracle	55.091 (9.457)	0.901	1.021	0.000	0.000
		ASE-D	52.116 (10.537)	0.827	1.038	19.948	0.000
		ASE-R	63.010 (17.939)	0.779	1.078	19.834	0.000
		GEE	245.548 (16.103)	0.810	1.006	0.048	0.000
	0.2	Oracle	118.769 (14.745)	0.928	1.008	0.000	0.000
		ASE-D	113.290 (16.222)	0.904	1.010	19.997	0.000
		ASE-R	120.243 (16.195)	0.882	1.012	19.986	0.000
		GEE	488.373 (25.909)	0.892	1.002	0.066	0.000
0.7	0.4	Oracle	40.677 (8.140)	0.880	1.042	0.000	0.000
		ASE-D	40.920 (8.357)	0.769	1.052	19.981	0.000
		ASE-R	41.707 (8.640)	0.781	1.051	19.950	0.000
		GEE	170.199 (18.057)	0.695	1.005	0.053	0.000
	0.3	Oracle	76.714 (12.918)	0.918	1.013	0.000	0.000
		ASE-D	76.442 (12.987)	0.862	1.015	20.000	0.000
		ASE-R	77.342 (13.187)	0.882	1.015	19.999	0.000
		GEE	303.726 (23.976)	0.830	1.003	0.072	0.000
	0.2	Oracle	172.986 (18.759)	0.931	1.005	0.000	0.000
		ASE-D	173.148 (19.087)	0.919	1.005	20.000	0.000
		ASE-R	173.864 (19.172)	0.917	1.005	20.000	0.000
		GEE	676.856 (36.565)	0.905	1.001	0.108	0.000
* Standard deviations are given in parentheses

Table 2. Table 2: Performance in variable identification and estimation of nonzero components under sequential sampling method based on ASE and GEE for a logit link function and ar (1) correlation between clusters when p 0 : p k = 3 : 12 : subscript 𝑝 0 subscript 𝑝 𝑘 3 : 12 p_{0}:p_{k}=3:12 .

$\alpha$	$d$	Method	$N$	$C.P.$	$\kappa$	$Nu{m}_c$	$Nu{m}_{ic}$
0.3	0.5	Oracle	384.997 (41.880)	0.949	1.004	0.000	0.000
		ASE-D	297.593 (52.953)	0.873	1.022	11.675	0.434
		ASE-R	415.522 (114.443)	0.858	1.016	11.618	0.342
		GEE	1752.497 (77.876)	0.937	1.001	0.014	0.000
	0.4	Oracle	603.168 (56.069)	0.950	1.003	0.000	0.000
		ASE-D	381.430 (93.703)	0.923	1.014	11.872	0.402
		ASE-R	604.040 (167.815)	0.908	1.010	11.873	0.279
		GEE	2680.709 (93.850)	0.938	1.001	0.015	0.000
	0.3	Oracle	1070.389 (72.207)	0.956	1.002	0.000	0.000
		ASE-D	611.206 (175.545)	0.943	1.007	11.973	0.331
		ASE-R	1090.787 (248.695)	0.951	1.005	11.983	0.167
		GEE	4693.950 (123.752)	0.950	1.000	0.021	0.000
0.7	0.5	Oracle	384.904 (43.847)	0.941	1.004	0.000	0.000
		ASE-D	296.432 (54.766)	0.902	1.020	11.660	0.454
		ASE-R	412.007 (108.277)	0.855	1.017	11.615	0.345
		GEE	1748.791 (73.967)	0.948	1.001	0.012	0.000
	0.4	Oracle	602.695 (54.798)	0.951	1.003	0.000	0.000
		ASE-D	383.120 (92.708)	0.916	1.014	11.879	0.383
		ASE-R	615.738 (162.033)	0.903	1.010	11.861	0.265
		GEE	2685.861 (92.699)	0.955	1.001	0.015	0.000
	0.3	Oracle	1071.102 (74.639)	0.953	1.002	0.000	0.000
		ASE-D	616.704 (173.345)	0.921	1.007	11.969	0.314
		ASE-R	1087.811 (250.271)	0.941	1.006	11.974	0.169
		GEE	4696.926 (129.885)	0.961	1.000	0.019	0.000

Table 3. Table 3: Sample size, number of TFs selected, and estimates of the correlation parameter in the yeast cell-cycle process with ASE-R and ASE-D.

Correlation	Method	$N$	${\hat{p}}_0$	$\hat{\alpha }$
ar(1)	ASE-R	240(26.20)	22.7(5.60)	0.498(0.021)
	ASE-D	133 (8.37)	21.0(5.96)	0.379(0.045)
Exchangeable	ASE-R	240(26.20)	21.3(5.63)	0.339(0.023)
	ASE-D	133 (7.88)	20.2(5.63)	0.215(0.050)
Independence	ASE-R	231(26.80)	15.0(3.31)	–
	ASE-D	129 (6.52)	17.1(4.32)	–
* Empirical standard deviations are given in parentheses

Table 4. Table 4: Sample size, number of simulation runs for which the entire sample has been used, number of variables selected, and correlation parameter estimates in the multiple sclerosis data process with ASE-R and ASE-D.

Correlation	Method	$N$	$N^{+}$	${\hat{p}}_0$	$\hat{\alpha }$
ar(1)	ASE-R	40.8(9.00)	310	3.31(0.493)	-0.080(0.013)
	ASE-D	36.2(6.72)	68	3.07(0.253)	-0.072(0.012)
Exchangeable	ASE-R	42.0(9.67)	368	3.38(0.489)	0.018(0.013)
	ASE-D	37.0(7.50)	110	3.12(0.325)	0.026(0.014)
Independence	ASE-R	41.6(9.11)	360	3.55(0.512)	–
	ASE-D	36.2(6.87)	89	3.09(0.290)	–
* Empirical standard deviations are given in parentheses

Table 5. Table 5: Number of coefficients that shrunk to zero over 1000 simulations of each variable.

Coefficients	ar(1)		Exchangeable		Independence
	ASE-R	ASE-D	ASE-R	ASE-D	ASE-R	ASE-D
T	0	0	0	0	0	0
T2	994	998	994	997	992	996
Gender	571	900	514	845	539	867
Durr	1000	1000	1000	1000	999	999
EDSS	2	2	4	3	4	1
Ltrt	1000	923	951	994	943	997
Htrt	230	57	166	61	188	60

Equations

$$\mathbf{S}_{n}(\beta)=\sum_{i=1}^{n}\mathbf{X}_{i}^{\mathrm{\scriptscriptstyle T}}\mathbf{A}_{i}(\beta)\mathbf{V}_{i}^{-1}(\beta)\{\mathbf{y}_{i}-\mathbf{h}_{i}(\beta)\bigr{\}}=0,$$

$${\mathbf{H}}_n(\beta )=\sum _{i=1}^n{\mathbf{X}}_i^{\mathrm{T}}{\mathbf{A}}_i(\beta ){\mathbf{V}}_i^{-1}(\beta ){\mathbf{A}}_i^{\mathrm{T}}(\beta ){\mathbf{X}}_i.$$

$${\mathbf{M}}_n(\beta )=\sum _{i=1}^n{\mathbf{X}}_i^{\mathrm{T}}{\mathbf{A}}_i(\beta ){\mathbf{V}}_i^{-1}(\beta ){\mathbf{ϵ}}_i(\beta ){\mathbf{ϵ}}_i^{\mathrm{T}}(\beta ){\mathbf{V}}_i^{-1}(\beta ){\mathbf{A}}_i^{\mathrm{T}}(\beta ){\mathbf{X}}_i.$$

$$\underset{i\ge 1}{\sup }E\{\Vert {\mathbf{e}}_i{\Vert }^2\mid {\mathcal{F}}_{i-1}\}<\infty \text{almost surely}.$$

$${\mathbf{M}}_n^{-1/2}{\mathbf{H}}_n/{\rho (n)}^{1/2}\to {\mathbf{\Sigma }}^{1/2}\text{ almost surely.}$$

$$\Vert {\tilde{\beta }}_n-{\beta }_0\Vert =o(\{({\overline{\lambda }}_n\log {\overline{\lambda }}_n{)}^{1/2}(\log \log {\overline{\lambda }}_n{)}^{1/2+\alpha }\}/{\underline{\lambda }}_n).$$

$${L_r}^{1/2}\kappa \to 0\text{and}{L_r}^{1/2+\gamma \delta }\kappa \to \infty ,\text{as}n\to \infty .$$

$$(O_{n}\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})(O_{n}\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})^{\mathrm{\scriptscriptstyle T}}=\left(\begin{array}[]{cc}{\Sigma_{11}(n)}_{\hat{p}_{0}\times\hat{p}_{0}}&{\Sigma_{12}(n)}_{\hat{p}_{0}\times(p_{0}-\hat{p}_{0})}\\ {\Sigma_{21}(n)}_{(p-\hat{p}_{0})\times\hat{p}_{0}}&{\Sigma_{22}(n)}_{(p-\hat{p}_{0})\times(p-\hat{p}_{0})}\\ \end{array}\right).$$

$$\displaystyle\;O_{n}I_{n}(\epsilon)\bigl{(}(\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})(\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})^{\mathrm{\scriptscriptstyle T}}\bigr{)}^{-1}I_{n}(\epsilon){O_{n}}^{\mathrm{\scriptscriptstyle T}}$$

$$=$$

$$\displaystyle\;(\widehat{\beta}_{n}-\beta_{0})^{\mathrm{\scriptscriptstyle T}}\bigl{(}I_{n}(\epsilon)((\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})(\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})^{\mathrm{\scriptscriptstyle T}})^{-1}I_{n}(\epsilon)\bigr{)}^{-}(\widehat{\beta}_{n}-\beta_{0})$$

$$=$$

$$({\hat{\beta }}_{N1}-{\beta }_{01}{)}^{\mathrm{T}}{\tilde{\Sigma }}_{11}(N)({\hat{\beta }}_{N1}-{\beta }_{01})\to {\chi }^2(p_0)\text{in distribution.}$$

$$N=N_d=\inf \left\{k:k\ge n_0\text{and}{\nu }_k\le \frac{d^2}{a_k^2}\right\},$$

$$R_N=\left\{\beta \in R^p:S_N\le \frac{d^2}{{\nu }_N}\text{mbox}and{\beta }_j=0\text{mbox}for{I}_{Nj}(ϵ)=0,1\le j\le p\right\},$$

$$\mathbf{X}^{*}=\arg\max_{\mathbf{X}^{*}\in\mathcal{U}}\det\Bigl{(}\mathbf{H}_{(n+1)1}(\widehat{\beta}_{n})\mathbf{M}_{(n+1)1}^{-1}(\widehat{\beta}_{n})\mathbf{H}_{(n+1)1}(\widehat{\beta}_{n})\Bigr{)}.$$

$${\mathbf{X}}^{\ast }$$

$${\mathbf{G}}_n=\sum _{i=1}^n{\mathbf{g}}_i=\sum _{i=1}^n{\mathbf{X}}_i{\mathbf{A}}_i({\hat{\beta }}_n){\hat{\mathbf{V}}}_i^{-1}({\hat{\beta }}_n){\mathbf{A}}_i^{\mathrm{T}}({\hat{\beta }}_n){\mathbf{X}}_i^{\mathrm{T}}=\sum _{i=1}^n{\mathbf{X}}_i{\mathbf{A}}_i^{1/2}({\hat{\beta }}_n){\bar{\mathbf{R}}}^{-1}{\mathbf{A}}_i^{1/2}({\hat{\beta }}_n){\mathbf{X}}_i^{\mathrm{T}}$$

$${\mathbf{y}}_i={\mathbf{X}}_i^{\mathrm{T}}{\beta }_0+{\mathbf{ϵ}}_i,$$

$$\text{logit}({\mu }_{ij})={\mathbf{x}}_{ij}^{\mathrm{T}}{\beta }_0,i=1,\dots ,5000,j=1,2,3,$$

$$y_{ij}={\alpha }_0+{\alpha }_1{t}_{ij}+\sum _{k=1}^{60}{\beta }_k{\mathbf{x}}_{ik}+{ϵ}_{ij},$$

$$\text{Ltrt}=\{\begin{array}{ll}1,& \text{Low Dose}\\ 0,& \text{Otherwise},\end{array}\text{Htrt}=\{\begin{array}{ll}1,& \text{High Dose}\\ 0,& \text{Otherwise}.\end{array}$$

$$\text{logit}({\mu }_{ij})={\beta }_1{\text{T}}_j+{\beta }_2{\text{T}}_j^2+{\beta }_3{\text{Gender}}_i+{\beta }_4{\text{Dur}}_i+{\beta }_5{\text{EDSS}}_i+{\beta }_6{\text{Ltrt}}_i+{\beta }_7{\text{Htrt}}_i,$$

$${\mathbf{M}}_n^{-1/2}{\mathbf{H}}_n({\hat{\beta }}_n-{\beta }_0)=$$

$$\displaystyle\times\bigl{\{}\mathbf{H}_{n}^{-1/2}\mathbf{M}_{n}^{1/2}\bigr{\}}^{-1}\mathbf{M}_{n}^{-1/2}\mathbf{S}_{n}(\beta_{0}),$$

$${\Vert {\mathbf{X}}_i{\mathbf{A}}_i{\mathbf{V}}_i^{-1}{\mathbf{e}}_i\Vert }^2$$

$$\le {\lambda }_{\max }\left({\mathbf{X}}_i{\mathbf{A}}_i{\mathbf{V}}_i^{-1}\mathrm{cov}({\mathbf{y}}_i\mid {\mathcal{F}}_{i-1}){\mathbf{V}}_i^{-1}{\mathbf{A}}_i^{\mathrm{T}}{\mathbf{X}}_i\right)\le c.$$

$$\mathrm{pr}\left\{\underset{k\le n\delta }{\max }|{\mathbf{S}}_{n+k}-{\mathbf{S}}_n|\ge ϵ\sqrt{n}\right\}\le \left(\frac{c}{n{ϵ}^2}\right)n\delta =\frac{c\delta }{{ϵ}^2}.$$

$$\sqrt{\rho (n)}({\hat{\beta }}_n-{\beta }_0)$$

$$\sqrt{\rho (n)}({\tilde{\beta }}_n-{\beta }_0)=$$

$$+\sqrt{\rho (n)}{\mathbf{\Sigma }}^{-1/2}{\mathbf{M}}_n^{-1/2}{\mathbf{H}}_n({\tilde{\beta }}_n-{\beta }_0).$$

$$\biggl{\{}\rho(N)(\widehat{\beta}_{N}-\beta_{0})^{\mathrm{\scriptscriptstyle T}}\{I_{N}(\epsilon){\mathcal{M}_{N}}^{-1}I_{N}(\epsilon)\}^{-}(\widehat{\beta}_{N}-\beta_{0})\leq\frac{d^{2}\rho(N)}{\nu_{N}}\text{ and }\beta_{j}=0\text{ for }I_{Nj}(\epsilon)=0\biggr{\}}$$

Full text

Adaptive design,

Adaptive sampling,

Maximum quasi-likelihood estimate,

Generalized estimating equations,

Stopping time,

and

1 Introduction

Correlated or highly stratified response data are common in studies where subjects are observed at multiple time points (Diggle et al., 2002), such as in some medical, epidemiological, and financial studies and other likewise longitudinal studies. Liang and Zeger (1986) proposed the generalized estimating equations (GEE) method, as an extension of generalized linear models (Wedderburn, 1974), for analyzing this kind of correlated responses. Due to modern computational/informational techniques, it is now easier in those scenarios to collect a long list of variables with highly clustered observations. However, from both a theoretical and computational perspective, statisticians are still hard-pressed to efficiently decide the effective variables when building a model using GEE methods.

In addition to variable determination issues, there are usually many redundant data points in a modern, automatically collected, disorganized data set, which becomes a data analysis burden, especially when the number of such data points is large and only limited computational capacity is available. Using a smaller, representative subset from a large data set would be an easy way to overcome such a computational obstacle. Nevertheless, to achieve this, we need to address new challenges, such as subset size and how to obtain a representative subset from such a data set through a statistically legitimate way. This situation becomes more complicated when there is a lengthy list of variables in the data set. Thus, it is important to have an efficient way to simultaneously determine important variables and detect the informative data points for analyzing this kind of data. Moreover, model interpretation would suffer if a large number of variables are included in the model. Thus, how to decide the effective variables to use in a model is essential for improving our ability to interpret it.

To balance the statistical inference needs and computational cost due to data set size, our goal is to use a sufficient amount of data to build a statistically interpretable model with effective selected variables. Sequential methods can help to identify and select the informative data points from a data pool for our analysis goal without using all of them at a time. This motivated us to study sequential methods with both adaptive variable and sample selection features for GEE-based estimation problems in highly stratified multiple response data.

Sequential methods are commonly used when there is no fixed sample size solution(see Siegmund, 1985). The idea of sequential methods emerged in the 1940s (Wald, 1945) and many researchers apply this concept, under different statistical setups, to many kinds of applications, such as clinical trials in medical studies, quality control in industry applications, computerized educational testing, and so on (Lai, 2001; Chang and Lu, 2010; Bartroff, Lai and Shih, 2013; Wainer, 2014; Tartakovsky, Nikiforove and Basseville, 2015; Park and Chang, 2016).

Conventional sequential methods analyze data when they become available so that we can end an experiment as soon as we have a satisfactory result with the desired statistical properties. We adopt the idea of sequential analysis, but recruit data from an existing, large data pool instead of collecting new observations, which differs from conventional sequential analysis. For a given precision in the parameter estimates, we sequentially find the most informative data points from an existing data set, for a GEE procedure, and simultaneously select a high-impact subset of variables during the estimation process. This kind of method allows us to retrench the sample size used without diminishing study quality. Taking advantage of the ability of sequential analysis to deal with adaptive/random subject recruiting and using a statistical experimental design criterion, we adaptively recruit new subjects, into the analysis to diminish the computational obstacles caused by large sample sizes. This adaptive feature of our procedure makes it closely related to stochastic regression, sequential experimental design, and nonlinear optimal design (Lai and Wei, 1982; Wu, 1985; Fedorve and Leonov, 2014)..

In the rest of this paper, we first briefly review the method of obtaining a maximum quasi-likelihood estimate (MQLE) and then present a sequential estimation procedure based on it for correlated multiple response data. This procedure features adaptive data selection and impact variable detection and builds a fixed-size confidence set for these effective variables when the data recruiting procedure is stopped. In Section 3, we use both simulated data and real examples to illustrate the proposed methods for various models and under different estimation strategies. This is then followed by the Conclusion section. We present technical details and supplemental numerical results in the Appendix and Supplementary materials.

2 Methodology

2.1 Model

Let $(y_{ij},\mathbf{x}_{ij})$ denote a pair representing the $j$th measurement on the $i$th subject, $j=1,\ldots,m_{i}$ and $i=1,\ldots,n$, where $y_{ij}$ is a scalar response and $\mathbf{x}_{ij}$ is a $p\times 1$ covariate vector. Observations from the same subject are assumed to be correlated; otherwise, they are independent. Let $\mathbf{y}_{i}=(y_{i1},\ldots,y_{im_{i}})^{\mathrm{\scriptscriptstyle T}}$, $i=1,\ldots,n_{i}$, be the vector of responses for the $i$th cluster and $\mathbf{X}_{i}=(\mathbf{x}_{i1},\ldots,\mathbf{x}_{im_{i}})^{\mathrm{\scriptscriptstyle T}}$ be the associated $m_{i}\times p$ matrix of covariates. We assume $m_{i}=m<\infty$ for simplicity. Let $\mathbf{h}_{i}(\beta)=\mathbf{h}(\mathbf{X}_{i}\beta)=E\{\mathbf{y}_{i}\mid\mathbf{X}_{i},\beta\}$ be the conditional expectation of $\mathbf{y}_{i}$ given $\mathbf{X}_{i}$ and $\beta$, where $\mathbf{h(\cdot)}$ is an $m$-dimensional sufficiently smooth one-to-one real-valued function and $\beta$ is the unknown parameter vector to be estimated. Wedderburn (1974) proposed an estimate for $\beta$, denoted by $\tilde{\beta}_{n}$, through solving the following equations:

[TABLE]

where $\mathbf{A}_{i}(\beta)=\mathbf{A}(\mathbf{X}_{i}\beta)=\partial\mathbf{h}^{\mathrm{\scriptscriptstyle T}}(t)/\partial t|_{t=\mathbf{X}_{i}^{\mathrm{\scriptscriptstyle T}}\beta}$ and $\mathbf{V}_{i}(\beta)=\mathbf{V}(\mathbf{X}_{i}^{\mathrm{\scriptscriptstyle T}}\beta)$ is a suitably-chosen known function of $\operatorname{cov}\{\mathbf{y}_{i}|\mathbf{X}_{i},\beta\}$. $\tilde{\beta}_{n}$ is commonly called the MQLE. If we specify the covariance of $\mathbf{y}_{i}$ as $\mathbf{V}_{i}(\beta,\alpha)=\mathbf{A}_{i}^{1/2}(\beta)\mathbf{R}_{i}(\alpha)\mathbf{A}_{i}^{1/2}(\beta)$ through a working matrix $\mathbf{R}_{i}(\alpha)$, then (1) is a generalized estimating equation (GEE), as the one in Liang and Zeger (1986), with equal-sized clusters. Note that the $\tilde{\beta}_{n}$ derived from (1) requires no distributional assumption.

Let $\boldsymbol{\epsilon}_{i}(\beta)=\mathbf{y}_{i}-\mathbf{h}_{i}(\beta)$. Then, following notations similar to those used in Xie and Yang (2003) and Balan and Schiopu-Kratina (2005), let

[TABLE]

Then, replacing the true (unknown) correlation matrix of the response variable by $\boldsymbol{\epsilon}_{i}(\beta)\boldsymbol{\epsilon}_{i}^{\mathrm{\scriptscriptstyle T}}(\beta)$ in the covariance matrix of $\mathbf{S}_{n}(\beta)$, we have

[TABLE]

Let $\beta_{0}\in R^{p}$ be the unknown true parameter vector to be estimated. (To simplify our notations, we will suppress $\beta_{0}$ (or $\widehat{\beta}_{n}$) in the sequel when it is clear that a term is a function of $\beta$ and evaluated at $\beta_{0}$ (or $\widehat{\beta}_{n}$), that is, $\mathbf{H}_{n}=\mathbf{H}_{n}(\beta_{0})$ and $\mathbf{M}_{n}=\mathbf{M}_{n}(\beta_{0})$ or, likewise, $\widehat{\mathbf{H}}_{n}=\widehat{\mathbf{H}}_{n}(\widehat{\beta}_{n})$ and $\widehat{\mathbf{M}}_{n}=\widehat{\mathbf{M}}_{n}(\widehat{\beta}_{n})$.)

2.2 Sequential method for MQLE

Using the notations defined before, let $\mathcal{F}_{n}=\sigma\{(\mathbf{y}_{j},\mathbf{X}_{j}):j=1,\ldots,n\}$ for $n\geq 1$ and let $\mathcal{F}_{0}=\sigma\{\emptyset,\Omega\}$; then, $\{\mathcal{F}_{n}:n\geq 0\}$ is an increasing sequence of $\sigma$-fields. Assume that for each $i$, the observed data satisfy $E\{\mathbf{y}_{i}\mid\mathcal{F}_{i-1}\}=\mathbf{h}_{i}(\beta_{0})$. Define $\mathbf{e}_{i}=\mathbf{y}_{i}-E\{\mathbf{y}_{i}\mid\mathcal{F}_{i-1}\}$ for each $i\geq 1$; then, $\{\mathbf{e}_{n}\}$ is a martingale difference sequence with respect to the $\sigma$-fields $\mathcal{F}_{n}$. Following Lai and Wei (1982), we also call Equation (1) a stochastic quasi-likelihood estimating equation.

If there are only a small number $p_{0}$ out of the $p$ components of $\beta_{0}$ that truly have an impact on the response and we know which they are in advance, then fitting a model with these $p_{0}$ variables, without using all $p$ variables, will certainly be more efficient – requiring fewer observations and less computational time. In practice, we usually have no such information, nor the number $p_{0}$, and fitting a model with all of the variables would increase the computational cost and instability of the model. The situation will be worse when $p$ is much larger than $p_{0}$. (Throughout the rest of this paper, we will refer to the $p_{0}$ variables as ‘effective variables’.) The conventional sequential estimation procedure, with no ability to detect the effective variables, cannot work well in this situation, which motivates us to incorporate the idea of adaptive shrinkage estimate (ASE) (Wang and Chang, 2013) to the current procedure.

There is a great deal of discussion about the asymptotic properties of the estimate for stochastic linear regression. Here, without linearity assumptions, we report some asymptotic results for MQLEs under conditions similar to those used in Lai, Robbins and Wei (1979); Lai and Wei (1982) for adaptive covariate/design vectors:

1. (C1)

For any $t\in R^{q}$, $\mathbf{V}(t)>0$, $\det\mathbf{A}(t)\neq 0$, each element of $\mathbf{V}(t)$ is continuously differentiable, and $\mathbf{h}(t)$ is twice continuously differentiable. 2. (C2)

The matrix $\mathbf{X}_{i}$, $i\geq 1$, satisfies $\sup_{i\geq 1}\|\mathbf{X}_{i}\|<\infty$ almost surely, where $\|\cdot\|$ stands for the Euclidean norm. In addition, the maximum and minimum eigenvalues of $\sum_{i=1}^{n}\mathbf{X}_{i}\mathbf{X}_{i}^{\mathrm{\scriptscriptstyle T}}$, denoted respectively as $\bar{\lambda}_{n}$ and $\underline{\lambda}_{n}$, satisfy $\underline{\lambda}_{n}\rightarrow\infty$ and {\color[rgb]{0,0,0}\lim\inf_{n\rightarrow\infty}\underline{\lambda}_{n}/\{(\bar{\lambda}_{n}\log\bar{\lambda}_{n})^{\frac{1}{2}}(\log\log\bar{\lambda}_{n})^{\frac{1}{2}+\alpha}\}>0} almost surely for some $\alpha>0$. 3. (C3)

The martingale difference sequence $\{\mathbf{e}_{n},n\geq 1\}$ with respect to an increasing sequence of $\sigma$-fields $\{\mathcal{F}_{n},n\geq 1\}$ satisfies

[TABLE] 4. (C4)

$\operatorname{cov}\{\mathbf{e}_{i}\mid\mathcal{F}_{i-1}\}>c{I}_{q}$ for all $i\geq 1$, and $\sup_{i\geq 1}E\{\|\mathbf{e}_{i}\|^{r}\mid\mathcal{F}_{i-1}\}<\infty$ almost surely for some $r>2$, where $c$ is a positive constant independent of $n$ and $I_{q}$ is a $q\times q$ identity matrix. 5. (C5)

There exists a non-random positive definite symmetric matrix $\tilde{\mathbf{M}}_{n}$, for which

$\tilde{\mathbf{M}}_{n}^{-1/2}\mathbf{M}_{n}\tilde{\mathbf{M}}_{n}^{-1/2}\rightarrow{I}_{p}$ in probability, where $\tilde{\mathbf{M}}_{n}^{-1/2}$ is the positive definite symmetric square root of $\tilde{\mathbf{M}}_{n}$. 6. (C6)

There exists a continuously increasing function $\rho(\cdot)$ and a non-random positive definite symmetric matrix $\boldsymbol{\Sigma}$ such that

[TABLE]

Assume that conditions (C1)–(C5) are satisfied. Then, we have

Lemma 1**.**

The sequence of random variable \bigl{\{}\widehat{\mathbf{M}}_{n}^{-1/2}\widehat{\mathbf{H}}_{n}(\widehat{\beta}_{n}-\beta_{0}),n=1,2,\ldots\bigr{\}} is uniformly continuous in probability (u.c.i.p.).

2.3 Adaptive shrinkage estimate

Under Conditions (C1)–(C3), Yin, Zhao and Wei (2006); Yin, Zhang and Zhao (2008) showed that the convergence rate of the MQLE $\tilde{\beta}_{n}$ is:

[TABLE]

Hence, if we let $L_{r}=\{(\overline{\lambda}_{n}\log\overline{\lambda}_{n})^{1/2}(\log\log\overline{\lambda}_{n})^{1/2+\alpha}\}/\underline{\lambda}_{n}$, then, as $n\rightarrow\infty$, with probability one, ${L_{r}}^{1/2}\kappa|\tilde{\beta}_{nj}|^{-\gamma}\longrightarrow 0\times I(\beta_{0j}\neq 0)+\infty\times I(\beta_{0j}=0)$, where $I(\cdot)$ is the indicator function (presuming $\infty\times 0=0$) and $\kappa=\kappa(n)$ is a non-random function of $n$ such that for $0<\delta<\frac{1}{2}$ and $\gamma>0$,

[TABLE]

(Note that (4) can be easily satisfied; for example, we can take $\kappa={L_{r}}^{-\theta}$ with $\theta\in(\frac{1}{2},\frac{1}{2}+\gamma\delta)$.) Let $I_{nj}(\epsilon)=I\{{L_{r}}^{1/2}\kappa|\tilde{\beta}_{nj}|^{-\gamma}<\epsilon\}$ for each $j$; then, by (3), we can use $I_{nj}(\epsilon)$ to detect whether the absolute difference between the $j$th component of $\beta_{0}$ and zero is asymptotically significantly larger than a predetermined constant $\epsilon(>0)$. Let $I_{n}(\epsilon)=\operatorname{diag}\{I_{n1}(\epsilon),\ldots,I_{np}(\epsilon)\}$; then, $\widehat{\beta}_{n}\equiv I_{n}(\epsilon)\tilde{\beta}_{n}$ defines an ASE of $\beta_{0}$. Thus, we have the following theorem:

Theorem 1**.**

Let $N(t)$ be a positive integer-valued random variable for which $N(t)/t$ converges to $1$ in probability as $t\rightarrow\infty$. Assume that conditions in Lemma 1 and (C6) are satisfied, then for any small $\epsilon>0$, as $n\rightarrow\infty$ we have $\sqrt{\rho(N(t))}(\widehat{\beta}_{N(t)}-\beta_{0})\rightarrow N(0,I_{0}\boldsymbol{\Sigma}^{-1}I_{0})$ in distribution, where I_{0}=\operatorname{diag}\bigl{\{}I(\beta_{01}\neq 0),\ldots,I(\beta_{0p}\neq 0)\bigr{\}} is a $p\times p$ diagonal matrix.

For convenience, we rearrange the matrix according to the effective variables detected at the current stage. Let $O_{n}$ be an orthonormal matrix satisfying $O_{n}^{\mathrm{\scriptscriptstyle T}}O_{n}=I_{p}$ so that $(\widehat{\beta}_{n1}^{\mathrm{\scriptscriptstyle T}},\widehat{\beta}_{n2}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}=O_{n}\widehat{\beta}_{n}$, where $(\widehat{\beta}_{n1}^{\mathrm{\scriptscriptstyle T}},\widehat{\beta}_{n2}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle T}}$ are the order-rearranged components of $\widehat{\beta}_{n}$ in which the indicators ${I_{nj}}$ corresponding to $\widehat{\beta}_{n1}$ are equal to 1, and the remaining ones, corresponding to $\widehat{\beta}_{n2}$, are equal to 0. We partition the matrix $(O_{n}\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})(O_{n}\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1/2})^{\mathrm{\scriptscriptstyle T}}$ according to the first $\hat{p}_{0}$ nonzero components of $O_{n}\widehat{\beta}_{n}$ so that

[TABLE]

With simple matrix algebra, we have that

[TABLE]

where $\tilde{\Sigma}_{11}^{-1}(n)={\Sigma^{-1}_{11}}(n)+{\Sigma^{-1}_{11}}(n)\Sigma_{12}(n){\Sigma^{-1}_{22.1}}(n)\Sigma_{21}(n){\Sigma^{-1}_{11}}(n)$ and ${\Sigma^{-1}_{22.1}}(n)=\Sigma_{22}(n)-\Sigma_{21}(n){\Sigma^{-1}_{11}}(n)$ $\Sigma_{12}(n).$ Let $M^{-}$ denote a general inverse matrix $M$. It follows that

[TABLE]

where $\beta_{01}$ is sub-vector of $\beta_{0}$ corresponding to $\widehat{\beta}_{n1}$. Using Theorem 1, it implies that as $t\rightarrow\infty$,

[TABLE]

2.4 Sequential estimate with adaptive shrinkage estimate

Our goal now is to have a sequential estimation procedure that can detect the effective variables and then construct a fixed-sized confidence set for them with a pre-specified accuracy. Theorem 1 paves the way for models with stochastic regressors. However, due to the unknown true value of $p_{0}$, which should also be estimated using the observations, we cannot use (9) directly. Conventional sequential estimation focuses on the sample size required for the parameter estimates; here, we need to decide the effective variables sequentially based on the recruited observations as well.

Suppose that $C_{k}$ is a set of $k$ observations $\{(\mathbf{y}_{i},\mathbf{X}_{i}):i=1,\ldots,k\}$ at the $k$th stage. For an $\epsilon>0$, let $\hat{p}_{0}(k)=\sum_{j=1}^{p}I_{kj}(\epsilon)$ based on $C_{k}$. Because $I_{nj}(\epsilon)$ converges to $I(\beta_{0j}\neq 0)$ almost surely as $n\rightarrow\infty$, this implies that $\hat{p}_{0}=\hat{p}_{0}(n)$ converges to $p_{0}$ almost surely; in addition, we have that $\lim_{n\rightarrow\infty}E\{\hat{p}_{0}(n)\}=p_{0}$ (see also Wang and Chang, 2013).

Let $a_{k}^{2}\in R$ be a constant satisfying the conditional probability $\mathrm{pr}(\chi^{2}_{\hat{p}_{0}(k)}\leq a_{k}^{2}|C_{k})=1-\alpha$ for a given $\alpha$. For a given $d>0$, define the stopping time $N_{d}$ as follows:

[TABLE]

where $\nu_{k}$ is the maximum eigenvalue of $I_{k}(\epsilon)(\widehat{\mathbf{H}}_{k}\widehat{\mathbf{M}}_{k}^{-1}\widehat{\mathbf{H}}_{k})^{-}I_{k}(\epsilon)$. We then conduct a sequential estimation procedure, collecting one new observation at a time until the stopping criterion $N_{d}$ is satisfied. If the inequality in (10) is fulfilled, then using all the $N_{d}$ observations, we construct a confidence ellipsoid for $\beta_{0}$:

[TABLE]

where $S_{N}=(\beta_{N1}-\widehat{\beta}_{N1})^{\mathrm{\scriptscriptstyle T}}\tilde{\Sigma}_{11}(N)(\beta_{N1}-\widehat{\beta}_{N1})$. It follows that the maximum axis of ${(\beta_{N1}-\widehat{\beta}_{N1})}^{\mathrm{\scriptscriptstyle T}}{\tilde{\Sigma}}_{11}(N){(\beta_{N1}-\widehat{\beta}_{N1})}={d^{2}}/{\nu_{N}}$, say $R_{N}^{\hat{p}_{0}}$, is $2d$. Please note that this $R_{N}^{\hat{p}_{0}}$ is only for the $\hat{p}_{0}$ detected effective variables, which are sequentially decided based on the observations.

The proposed method allows the sequential estimation procedure to focus on the effective variables so that the confidence ellipsoid $R_{N}^{\hat{p}_{0}}$ is a projection of $R_{N}$ onto the $\hat{p}_{0}$-dimensional space spanned by the axes with nonzero components of $\widehat{\beta}_{N}$. This is the reason why the proposed method uses fewer observations to achieve the required properties compared to those used in procedures with no variable detection feature. Theorem 2 summarises these properties; its proof is in the Appendix.

Theorem 2**.**

Let $N$ be the stopping time as defined in (10). If the assumptions of Theorem 1 hold, then (i) \lim_{d\rightarrow 0}{d^{2}{\color[rgb]{0,0,0}\rho(N)}}/{a^{2}\nu}=1 almost surely, (ii) $\lim_{d\rightarrow 0}\mathrm{pr}(\beta_{0}\in R_{N})=1-\alpha$, (iii) $\lim_{d\rightarrow 0}\hat{p}_{0}(N)=p_{0}$ almost surely, and (iv) $\lim_{d\rightarrow 0}E\{\hat{p}_{0}(N)\}=p_{0}$, where $\nu$ is the maximum eigenvalue of matrix $I_{0}\boldsymbol{\Sigma}^{-1}I_{0}$.

Theorem 2 parts (i) and (ii) were termed ‘asymptotic consistency and efficiency,’ respectively (Chow and Robbins, 1965), which mean that (i) the coverage probability of the proposed sequential procedure will converge to the prescribed one and (ii) the proposed sequential procedure is, asymptotically, as efficient as the (unknown) fixed sample procedure in terms of the ratio of the sample sizes. Theorem 2 parts (iii) and (iv) state that when the sequential estimation procedure is stopped, the estimate of the number of effective variables converges to the true number of effective variables in probability and its expectation will also converge to the true number of effective variables.

By (10), it is clear that the procedure requires a larger sample size for a smaller $d$. From Equation (11), we can see that the constant $d$ confines the maximum axis of the confidence ellipsoid, which specifies the precision of the final parameter estimates of a sequential procedure. Together, these guarantee that the proposed sequential procedure can use a ‘minimum’ required sample size such that the confidence set will approximately have the pre-specified coverage probability for the effective components of $\beta_{0}$. (Note that to sample a batch of new observations at a time is also feasible with slight modifications of the current arguments.)

Remark 1**.**

We can treat the conventional sequential procedure, without ASE, as a simplified case of the current procedure. Hence, we also summarize the results for the fixed-sized confidence set estimation without the variable detection feature in Supplementary materials.

Modified Sequential D-criterion

In the theorem above, the regularity conditions on the design are not for any specific design scheme. Using the D-optimality criterion in statistical experimental design to select points will fulfill those conditions. However, to find the new observations based on conventional D-optimality may be computational intensive. Therefore, we use a modified sequential D-criterion described below to select new observations adaptively and sequentially.

Let $\mathcal{D}$ and $\mathcal{U}$ be the recruited sample set and the inactive sample set, respectively. The regular D-criterion is to choose a data point $\mathbf{X}^{*}$ from the data pool, which satisfies the following equation:

[TABLE]

Because this procedure is usually computationally intensive, we propose the following modified sequential D-criterion instead:

[TABLE]

where

[TABLE]

and $\bar{\mathbf{R}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{A}_{i}^{-1/2}(\widehat{\beta}_{n})\{\mathbf{y}_{i}-\mathbf{h}_{i}(\widehat{\beta}_{n})\}\{\mathbf{y}_{i}-\mathbf{h}_{i}(\widehat{\beta}_{n})\}^{\mathrm{\scriptscriptstyle T}}\mathbf{A}_{i}^{-1/2}(\widehat{\beta}_{n})$, as suggested in Balan and Schiopu-Kratina (2005). (For a matrix $\mathbf{W}$, the notation $\mathbf{W}_{n}$ represents the matrix arranged according to the detected effective variables at stage $n$ and $\mathbf{W}_{(n)1}$ denotes its sub-matrix. After $\mathbf{X}^{*}$ is selected, we delete it from $\mathcal{U}$ and add $(\mathbf{y}^{*},\mathbf{X}^{*})$ to $\mathcal{D}$, then update $\widehat{\beta}_{n}$ to $\widehat{\beta}_{n+1}$. Then, we continue this process until the stopping criterion (10) is fulfilled.

3 Numerical Studies

3.1 Simulated Data

We first assess the proposed method with simulated data sets by comparing the average sample sizes (stopping times, $N^{*}$), empirical coverage probabilities (C.P.), and consistency of variable detection of the GEE-based procedures with ASE under random sampling (ASE-R), adaptive sampling with D-optimality criterion (ASE-D), and the procedure with all $p$ variables (GEE). Additionally, we report the results using only the true $p_{0}$ variables (the Oracle) as a reference. We also report the average numbers of incorrectly identified zero variables ($\text{Num}_{ic}$) and the average number of correctly identified zero variables ($\text{Num}_{c}$) for the performance in identifying effective variables. The targeted coverage probability is equal to $95\%$ and the empirical coverage probability should approach the nominal $95\%$ as $d$ decreases. Let $\kappa=d^{2}N/(a^{2}_{N}\nu_{N})$; then, as stated in Theorem 2, we expect that $\kappa$ will be close to 1 as $d$ decreases. We consider both continuous and discrete responses in our simulation study.

Continuous Response Case.

The highlight of this example is about the ability to detect effective variables and we generate data with a linear link function for illustration purposes. Let $\mathbf{X}_{i}=(\mathbf{x}_{i1},\ldots,\mathbf{x}_{im})$ be the covariate matrices with $m=5$ and $p=8,10$, and $24$, and for different $p$, the parameter $\beta_{0}=(1,-1.1,1.5,-2,\mathbf{0}_{k})^{\mathrm{\scriptscriptstyle T}}$ with $k=4,8$, and $20$, where $\mathbf{0}_{k}$ denotes the $k$-dimensional vector of [math]s. Thus, the ratio of the number of nonzero coefficients ($p_{0}$) to the number of zero coefficients ($p_{k}$) are $4:4$, $4:8$ and $4:20$. We generate this data set sequentially so that successors depend on the information of their predecessors in the following adaptive manner: Let the first $\mathbf{X}_{1}$ be from a standard multivariate normal distribution with zero mean vector and identity covariance matrix. We then adaptively generate $\mathbf{X}_{n}$s, $n>1$, from a multivariate normal distribution with mean vector $\sum_{i=1}^{n-1}\sum_{j=1}^{m}\mathbf{x}_{ij}/m(n-1)$ and identity covariance matrix. We then generate the response $\mathbf{y}_{i}$ as follows:

[TABLE]

where the random error vector $\boldsymbol{\epsilon}_{i}$ follows a normal distribution with mean $\mathbf{0}$ and three different covariance structures with corresponding dimensional numbers. These three covariance matrices are the identity matrix, the exchangeable, and the ar(1) autoregressive correlation structure with autocorrelations $\alpha=0.3$ and $0.5$. We only report the results of the ar(1) case with $p_{0}:p_{k}=4:20$ and present the rest of the results in Supplementary materials. In this example, we set the parameters for ASE, as described in (4), as follows: $\gamma=1$, $\delta=0.45$, and $\theta=0.65$ and use the qic (Pan, 2001) as its information criterion.

We start with a data pool of $1000$ data points each time and obtain an initial estimate of $\beta_{0}$ using a random sample of initial data of size $25$ from this data pool. Table 1 summarizes the results based on 1000 replications with $d=0.2,0.3$, and $0.5$. The results show that the empirical coverage probabilities of ASE-D, ASE-R, and GEE are all less than that of the Oracle, whereas the sample sizes used in ASE-D and ASE-R are very close to that of the Oracle. In this example, the estimate from the conventional GEE uses around 4 or 5 times larger sample sizes than those of the other methods. The incorrectly identified zero variables ($\text{Num}_{ic}$) are all equal to 0. For ASE-D and ASE-R, the average number of correctly identified zero variables ($\text{Num}_{c}$) is close to the true number – 20 and all GEE cases have $\text{Num}_{c}<0.1$, which is reasonable, since conventional GEE do not consider the sparseness situation. These results confirm that the proposed sequential procedure with ASE works well for the multiple correlated responses data in the current simulation setup. Both ASE-D and ASE-R perform similarly here, because this model is simple. The advantage of ASE-D will be clearly revealed in the discrete response case next.

Discrete Responses Case

In this example, we use a logistic model to illustrate discrete response situations. We generate the covariates vector $\mathbf{x}_{ij}$ from a multivariate normal distribution with mean zero and an ar(1) correlation matrix with autocorrelation coefficient $0.5$ and marginal variance equal to $0.2$. The binary response vector for each cluster has an ar(1) correlation structure with correlation coefficient $\alpha$ and the marginal expectation $\mu_{ij}$ satisfies the following equation:

[TABLE]

We use the R program in Oman (2009) to generate the correlated binary response data under several different setups: two values of $\alpha=0.3$ and $0.7$, and three different length regression coefficient vectors: $\beta_{0}^{\mathrm{\scriptscriptstyle T}}=(0.6,-0.5,0.4,\mathbf{0}_{k}^{\mathrm{\scriptscriptstyle T}}$) with $k=3,6$ and $12$.

At the beginning, 5000 observations are generated as the data pool and we randomly choose $200$ from this pool as our initial set. We consider two situations in this nonlinear case: (1) the structure of $\mathbf{R}(\alpha)$ is known and (2) the correlation structure $\mathbf{R}(\alpha)$ is mis-specified. In both situations, we estimate $\alpha$ based on the suggestion of Liang and Zeger (1986).

Table 2 reports the results of ar(1) with $p_{0}:p_{k}=3:12$. The GEE method using all $p$ variables has better coverage probabilities than the other two procedures at the cost of the sample size used. In fact, GEE usually uses samples approximately 4 to 5 times larger than those used by the Oracle procedure in this study. Both ASE-D and ASE-R have similar values of $\text{Num}_{c}$, which are close to the targeted $p_{k}=12$, whereas the $\text{Num}_{ic}$ values are slightly higher than those in the previous case due to the nonlinear model situation. In the current model, we clearly see that the sample sizes used in ASE-D and ASE-R are very different: ASE-D uses fewer observations by fully taking the advantage of the D-optimality criterion in selecting data points. The results for the other settings are in Supplementary materials.

3.2 Real Data Examples

Yeast Cell-Cycle Gene Expression Data Analysis

We apply the proposed MQLE-based sequential procedure with ASE to the yeast cell-cycle gene expression data set collected in CDC15 (Spellman et al., 1998, see also http://genome-www.stanford.edu/cellcycle/data/rawdata/). In this experiment, genome-wide mRNA levels were recorded for 6178 yeast ORFs (an abbreviation for open reading frames, DNA sequences that can determine which amino acids will be encoded by a gene) at 7 minute intervals for 119 minutes, which covers two cell-cycle periods, for a total of 18 time points. The cell-cycle is a tightly regulated life process, where cells grow, replicate their DNA, segregate their chromosomes, and divide into as many daughter cells as the environment allows, and this process is commonly divided into M/G1-G1-S-G2-M stages. The M stage stands for ‘mitosis’, during which nuclear (chromosome separation) and cytoplasmic (cytokinesis) divisions occur. The G1 stage stands for ‘GAP 1’;the S stage stands for ‘synthesis’, during which DNA replication occurs; and, the G2 stage stands for ‘GAP 2’.

A sequence-specific DNA-binding factor, sometimes referred to as TF, is a protein that binds to specific DNA sequences, thereby controlling the flow (or transcription) of genetic information from DNA to mRNA. In their experiment, Spellman et al. (1998) identified approximately 800 genes that vary in a periodic fashion during the yeast cell-cycle. The regulation of most of these genes was not clear; however, TFs have been observed to play critical roles in gene expression regulation (see also Simon et al., 2001). We apply the proposed method to identify the TFs that influence the gene expression level at each stage of the cell-cycle process. This is essential for understanding how the cell-cycle is regulated and also how cell-cycles regulate other biological processes.

We analyze a subset of 297 cell-cycle-regularized genes obtained from the R package PGEE (Inan and Wang, 2017) as in Luan and Li (2003) and Wang, Li and Huang (2008). The response variable $y_{ij}$ is the log-transformed gene expression level of gene $i$, measured at time point $j$; the covariate $\mathbf{x}_{ik}$, $k=1,\ldots,96$, is the matching score of the binding probability of the $k$th transcription factor on the promoter region of the $i$th gene. We calculated the binding probability with a mixture modeling approach based on data from a ChIP binding experiment (see Wang, Chen and Li, 2007, for further details) and, to compare the performance of ASE-D, ASE-R, and GEE, we removed 20 less correlated covariates due to the sample size available.

We apply the sequential estimation procedure to the G1 stage (which contains a few time points from the cycle) of the cell-cycle process using the following model

[TABLE]

where $\mathbf{x}_{ik},k=1,\ldots,60$, is standardized to have mean zero and variance 1; $t_{ij}$ denotes time and imposes shrinkage on the $\beta_{k}$’s. Table 3 summarizes the number of TFs identified based on 1000 replications, using three different working correlation structures for $\epsilon_{ij}$: independence, ar(1), and exchangeable. In each replication, the initial regression coefficient estimates are based on $100$ random observations.. Our analysis reveals that at the G1 stage, the selected TFs, in terms of numbers and specific TFs, are not sensitive to the choice of the working correlation structure. Due to the space limitation, in Table 3, we only report the average sample sizes, the average number of selected TFs and the average estimates of the correlation parameter. Some of the TFs selected have already been confirmed by biological experiments using the genome-wide binding method. For example, MBP1, SWI6, and SWI4 are three TFs that have been proved important in stage G1 in the aforementioned biological experiments and they have been selected by the proposed sequential procedure for stage G1. We can see from Table 3 that ASE-D uses fewer observations ASE-R based on their averages; in fact, the standard deviation of the sample sizes used by ASE-D is also smaller than that of the sample sizes used by ASE-R. The numerical results here clearly show that the proposed method is beneficial and the sequential procedure with D-optimality for adaptive data recruiting is promising.

Multiple Sclerosis Data Analysis

This is a longitudinal clinical trial data set for assessing the effects of neutralising antibodies on interferon beta-1b (IFNB) in relapsing–remitting multiple sclerosis (MS), a disease that destroys the myelin sheath surrounding the nerves. This data set is from a Magnetic Resonance Imaging (MRI) sub-study of the Betaseron clinical trial conducted at the University of British Columbia in relapsing–remitting multiple sclerosis, which involved 50 patients and each one visited the university every six weeks. The patients were randomly allocated into three treatment groups: 17 patients treated by placebo, 17 by a low dose, and 16 by a high dose. This data set was previously analyzed in Petkau (2003); Petkau et al. (2004) and included in a book by Song (2007).

Exacerbation is used as the binary response variable, which indicates whether an exacerbation appeared since the previous MRI scan—1 for ‘yes’ and 0 for ‘no’. Seven explanatory variables were recorded: Treatment (Trt), Time (T) in weeks, Squared time (T2), Age, Gender, Duration of disease (Dur) in years, and an additional baseline covariate—initial EDSS (Expanded Disability Status Scale) scores. Due to the results in Song (2007) and Li et al. (2013) and the available sample size, we also delete Age from the original data as in the analysis conducted in Li et al. (2013). We recode Ltrt and Htrt as

[TABLE]

and consider the marginal logistic model

[TABLE]

where $\mu_{ij}$ is the probability of exacerbation at visit $j$ for subject $i$. For illustration, we use three correlation structures: ar(1), exchangeable, and independence. Table 4 reports the average sample size ($N$), the number of simulation runs that use the entire sample ($N^{+}$), the average number of variables selected ($\hat{p}_{0}$) and the estimate of $\alpha$ ($\widehat{\alpha}$) based on 1000 runs. We start with an initial estimate of regression coefficients using $25$ randomly selected observations in each run. In Table 5, for each variable, we also summarize how many times the corresponding coefficient shrinks to zero during the 1000 runs. Our numerical results show that ASE-D tends to use fewer observations and can complete the estimation procedure with the available sample sizes. Although the differences may not be statistically significant, they clearly show that ASE-D tends to select fewer variables and this information is usually beneficial for practitioners to design future studies.

4 Discussion

The collection of correlated or highly stratified response data is common due to modern methods of data collection and to the monitoring methods and facilities. Data set sizes are larger than before, and, hence, using an entire data set at once may not be convenient when our computational power and skills are not adequate. Selecting a random subset of such a large data set is a common and easy solution; however, it is well-known that this method usually suffers from sampling variations and cannot always provide a stable and consistent result. Moreover, to decide the size of such a random subset has never been an easy task when complicated models are involved. Thus, how to efficiently use the effective observations is an important issue and the idea of a sequential estimation method seem appropriate for our needs. We propose a method providing sequential stochastic MQLEs so that we can adaptively select effective observations and adopt the ASE method to detect the effective variables simultaneously. The proposed method recruits new observations from a data pool into the analysis until a pre-specified stopping criterion is fulfilled. The asymptotic properties of the sequential procedure and numerical studies using both simulation data and real data examples show that the proposed method can perform well and appears to be a promising method for practical uses. Our method is not limited to selecting one new observation at a time. Selecting a batch of observations at each stage is possible with only slight modifications.

Acknowledgment

A part of this research was supported by the Ministry of Science and Technology, Taiwan (MOST 106-2118-M-001-007-MY2). The first author also gratefully acknowledges the financial support from the China Scholarship Council(Grant No. 201706340175).

Supplementary materials

We state the results of the sequential procedure with no ASE feature in Supplementary materials, together with some extra simulation results, including the cases of the logit link with the ASE feature and different correlation structures. Readers interested in the results under different setups may refer to Supplementary materials.

Appendix

Technical details

Proof of Lemma 1.

Let $\mathbf{S}_{n}(\beta_{0})=\sum_{i=1}^{n}\mathbf{X}_{i}\mathbf{A}_{i}\mathbf{V}_{i}^{-1}\mathbf{e}_{i}$. Because $\mathbf{X}_{i}$ is $\mathcal{F}_{i-1}$ measurable, $\{\mathbf{e}_{i}:i\geq 1\}$ is a sequence of martingale differences with respect to $\{\mathcal{F}_{i},i\geq 1\}$, and $\mathbf{S}_{n}(\beta_{0})$ is a sum of martingale differences. By Taylor series expansion, there exists $\bar{\beta}_{n}$, which lies between $\widehat{\beta}$ and $\beta_{0}$, such that

[TABLE]

where $\mathbf{S}_{n}(\beta_{0})=\mathcal{D}_{n}(\bar{\beta}_{n})(\widehat{\beta}_{n}-\beta_{0})$. Yin, Zhao and Wei (2006) showed that $\mathbf{H}_{n}^{-1/2}\mathcal{D}_{n}(\bar{\beta}_{n})\mathbf{H}_{n}^{-1/2}\rightarrow I_{p}$ almost surely under conditions (C1)–(C3). Note that each element of $\mathbf{M}_{n}^{-1/2}\mathbf{H}_{n}^{-1/2}$ and its inverse matrix are bounded. Hence, by Lemma 1.4 of Woodroofe (1982), to show $\{\mathbf{M}_{n}^{-1/2}\mathbf{H}_{n}(\widehat{\beta}_{n}-\beta_{0}):n\geq 1\}$ is u.c.i.p., it suffices to show that $\{\mathbf{S}_{n}(\beta_{0}):n\geq 1\}$ is u.c.i.p.

By the assumption that $\sup_{i\geq 1}\|\mathbf{X}_{i}\|<\infty$ almost surely, all $\mathbf{X}_{i}^{\mathrm{\scriptscriptstyle T}}\beta_{0}$ fall in a compact set $T$ of $\mathbf{R}^{q}$ with probability one. A result of Yin, Zhao and Wei (2006, Lemma 1) implies that for each $i\geq 1$ and $t\in T$, $c_{1}I_{q}<\mathbf{A}_{i}(t)<c_{2}I_{q}$ and $c_{3}I_{q}<\mathbf{V}_{i}^{-1}(t)<c_{4}I_{q}$ with probability one. Let $\lambda_{\max}(\mathbf{T})(\lambda_{\min}(\mathbf{T}))$ denote the largest(smallest) eigenvalue of the matrix $\mathbf{T}$. Then,

[TABLE]

Let $U_{n}=\mathbf{S}_{n}(\beta_{0})$ and $b_{n}=1$ in the Hájek-Rényi inequality (see Chow and Teicher, 1988, Theorem 7.4.8 (iii)). Following the arguments of Example 1.8 in Woodroofe (1982) and replacing the Kolmogorov’s inequality there with the Hájek-Rényi inequality for martingale differences, we have, for $\epsilon,\delta>0$ and $k\leq n\delta$,

[TABLE]

Because ${c\delta}/{\epsilon^{2}}$ is independent of $n$ and goes to zero as $\delta\rightarrow 0$, this implies that $\{\mathbf{S}_{n}(\beta_{0}),n\geq 1\}$ is u.c.i.p. ∎

Proof of Theorem 1.

Let us rewrite

[TABLE]

Because $I_{n}(\epsilon)$ converges to $I_{0}$ almost surely, this implies that $\Delta_{2}(n)$ converges almost surely to [math] as $n$ goes to infinity. It has been shown that $\mathbf{M}_{n}^{-1/2}\mathbf{H}_{n}(\tilde{\beta}_{n}-\beta_{0})\rightarrow N(0,I_{p})$ in distribution (Yin, Zhao and Wei, 2006). Then, by Slusky’s theorem and (C6), $\sqrt{\rho}(n)(\widehat{\beta}_{n}-\beta_{0})\rightarrow N(0,I_{0}\boldsymbol{\Sigma}^{-1}I_{0})$ in distribution. Because $I_{n}(\epsilon)$ is bounded, to replace the sample size with a random variable it is sufficient to show that $\{\sqrt{\rho(n)}(\tilde{\beta}_{n}-\beta_{0})\}$ is u.c.i.p. (see Anscombe, 1952).

[TABLE]

Because $\sqrt{\rho(n)}\mathbf{M}_{n}^{-1/2}\mathbf{H}_{n}(\tilde{\beta}_{n}-\beta_{0})$ is stochastic bounded, it follows from condition (C6) that the first term of (Proof of Theorem 1.) goes to [math] almost surely as $n$ goes to infinity. Using arguments similar to those in Lemma 1, the sequence $\{\sqrt{\rho(n)}\boldsymbol{\Sigma}^{-1/2}\mathbf{M}_{n}^{-1/2}\mathbf{H}_{n}(\tilde{\beta}_{n}-\beta_{0})\}$ is u.c.i.p, and the remain proof, which is omitted here, follows from the results of Anscombe (1952) (see also Woodroofe (1982)). ∎

Proof of Theorem 2.

We know that with probability one,

$\rho(k)I_{k}(\epsilon)(\widehat{\mathbf{H}}_{k}\widehat{\mathbf{M}}_{k}^{-1}\widehat{\mathbf{H}}_{k})^{-}I_{k}(\epsilon)$ converges to $I_{0}\boldsymbol{\Sigma}^{-1}I_{0}$ and $\rho(k)\nu_{k}$ converges to $\nu$ as $k\rightarrow\infty$. Let $y_{k}={(\rho(k)\nu_{k})}/{\nu}$, $f(k)={(\rho(k)a^{2})}/{a_{n}^{2}}$ and $t={(\nu a^{2})}/{d^{2}}$. Then, (10) becomes $N=\min\{k:k\geq n\ \text{and}\ y_{k}\leq f(k)/t\}.$ Hence, $1=\lim_{t\rightarrow\infty}{f(N)}/{t}=\lim_{d\rightarrow 0}{(d^{2}\rho(N)}){(a^{2}\nu)}$ almost surely. This implies that, as $t\rightarrow\infty$, ${d^{2}\rho(N)}/{\nu_{N}}\rightarrow a^{2}$ and $\rho(N)/t\rightarrow 1$ almost surely. Hence, (i) follows from Chow and Robbins (1965, Lemma 1).

Let $\mathcal{M}_{n}=\widehat{\mathbf{H}}_{n}\widehat{\mathbf{M}}_{n}^{-1}\widehat{\mathbf{H}}_{n}$. Hence, when $t$ goes to infinity, the event $\{\mathbf{\beta}_{0}\in R_{N}\}$ is equivalent to

[TABLE]

Thus, Theorem 2 (ii) follows from Theorem 1. For the proofs of (iii) and (iv), we can use arguments similar to those of Wang and Chang (2013, Theorem 4), so we are omitting them here.

∎