Testing the Correlation and Heterogeneity for Hierarchical Nonlinear Mixed-effects Models

Nonlinear mixed-effects models are very useful in analyzing repeated-measures data and have received a lot of attention in the field. It is of common interest to test for the correlation within clusters and the heterogeneity across different clusters. In this paper, we address these problems by proposing a class of score tests for the null hypothesis that all components of within-and between-subject variance are zeros in a kind of nonlinear mixed-effects model, and the asymptotic properties of the proposed tests are studied. The finite sample performance of this test is examined through simulation studies, and an illustrative example is presented.


Introduction
Repeated-measures data are frequent observations in different areas of investigation, such as economics and pharmacokinetics.For instance, in longitudinal studies, observations on the same subject are usually made at different times.Analysis of such data requires accounting for the within-cluster correlation and the between-subject heterogeneity of the data.Randomeffects models are commonly used for analyzing clustered and repeated-measures data.It is of no doubt that the linear mixed-effects models play an important role in evaluating and analyzing the repeated-measures and clustered data.For an example, see Laird and Ware 1 .However, many repeated-measures data, such as growth data, dose-response and pharmacokinetics data, are inherently nonlinear with respect to a given response regression function.Several different nonlinear mixed-effects model and various inference procedures have been proposed 2-5 .For models considered in these literatures, one of the interesting issues is whether there exists correlation within clusters or/and heterogeneity between subjects.It is well known that the misspecification of model can have serious impact

Nonlinear Mixed-Effects Model
In this paper, we consider the following nonlinear mixed-effects model as proposed by Pinheiro and Bates 5 .In the first stage the jth observation on the ith subject is modeled as ε ij , i 1, . . ., m, j 1, . . ., n i , 2.1 where f •, • is a twice-differentiable smooth nonlinear function of a subject-specific parameter vector φ ij and the predictor x ij , ε ij is a normally distributed noised term, m is the total number of subjects, and n i is the number of observations on the ith subject, and n m i 1 n i is the total number of observations.In the second stage, the subject-specific vector is modeled as where β is a p × 1 vector of unknown parameters, b i , i 1, . . ., m, are independent q i × 1 of random effects associated with the ith subject, and A ij and B ij are design matrices for fixed and random effects, respectively.The random effects b i account for the correlation with the same cluster.Furthermore, we make the following distributional assumptions.The ε ij , i 1, . . ., m, j 1, . . ., n i , are independently distributed as N 0, σ 2 m ij , where • is a twice differentiable function, γ is a r × 1 vector of unknown parameters, z ij are the known design vectors, and m z ij , γ 0 1 if and only if γ γ 0 , and T denote the vector obtained from stacking up the m cluster-specific entries and assume that b is generated from some distribution F with mean zero and covariance matrix D θ with D θ diag D 1 θ , . . ., D m θ , where θ is a s × 1 vector of unknown variance component.The magnitude of θ can be used to measure the degree of correlation and heterogeneity of the cluster-specific response vector within each subject.We postulate that each component of We further assume that the third-and higher-order moments of the random effects b i are of order o θ .These conditions are consistent with Lin 13 , Hall and Praestgaard 16 , and Zhu and Fung 14 .This model can be regarded as a hierarchical model that in some aspects generalized the linear mixed-effects model of Laird and Ware 1 , the usual nonlinear model for independent data of Bates and Watts 17 .
It should be pointed out that the distribution of the random error ε ij , i 1, . . ., m, j 1, . . ., n i , need not to be normal.In this paper, we still assume the ε ij , i 1, . . ., m, j 1, . . ., n i , to be normal.The main reason is to have the computation and deduction in mathematics become relatively simple.From the results presented in the below of this paper, it can be found that the derivation of the score test statistics is on the basis of the complete-data log likelihood, which requires making some specifications for the random error in model and further the random error is directly related to the complexity and difficulty of computation and deduction in mathematics.In nonlinear mixed-effects model, the random effects entering nonlinearly in model make the likelihood analysis of nonlinear mixed-effects model more difficult and complicated than that of their linear counterpart.It can be seen that even if the error distribution is assumed to be normal, the computation and deduction of the test score are cumbersome, let alone the error distribution is the others.In addition, the score test does not depend on the normality of the random error, but some good properties of the normal distribution help us obtain the test statistics.Therefore, we choose the normal distribution as that of the random error.
For model 2.1 and 2.2 , there is either the correlation within clusters or heterogeneity among observations.First, we investigate whether there exists correlation among observations in the same subject; it is equivalent to test D θ 0 or not.Thus we can use the hypothesis

Advances in Decision Sciences
for the test of the correlation within clusters in nonlinear mixed model.Second, we study whether there exists heterogeneity of between-subject variance and correlation within clusters at the same time.To address this problem, we can use the composite hypothesis The testing of random effects in the nonlinear mixed-effects models has been discussed in some literatures e.g., Jacqmin-Gadda and Commenges 12 , Hall and Praestgaard 16 , and Zhu and Fung 14 .However, the models they investigated are confined to the additive nonlinear mixed-effects model; that is, the random effects added to a nonlinear function or the rand noise is i.i.d random.Moreover, they generally studied the hypothesis 2.3 , which tested whether there exists the correlation within clusters.In our paper, we consider the testing of random effects under some considerably general conditions for hierarchical nonlinear mixed-effects models.We not only study the testing correlation within clusters, but also investigate whether there exists heterogeneity of between-subject variance at the same time. Let and M i diag m i1 , . . ., m in i .Denote the the vector obtained from stacking up the m cluster-specific entries of the same symbol by To derive score tests for null hypothesis in which θ θ 0 and θ θ 0 , γ γ 0 .We firstly study the properties of the log likelihood in model 2.1 and 2.2 .For a given b, the conditional log-likelihood of nonlinear mixed-effects model and its derivative are as follows: where where and fij is the second partial derivative of f φ ij , x ij with respect to φ ij evaluated at b 0, and then

Score Test for Correlation and Heterogeneity within Clusters
In this section, we firstly use the Laplace expansions to develop a score test for the null hypothesis H 0 : θ 0 under model 2.1 and 2.2 , which corresponds to having no correlation within clusters.Then using the same approach, we obtain a score test statistic for the null hypothesis H 0 : θ 0, γ 0 under model 2.1 and 2.2 , which corresponds to having no correlation within clusters and heterogeneity of between-subject variance at the same time.
Let ξ β T , σ 2 , γ T T , for an n × 1 data vector Y following model 2.1 and 2.2 .The likelihood function is Using the moment assumptions on the random effects b and the Laplace expansion similar to Lin 13 and Hall and Praestgaard 16 , we expand the integrated likelihood 3.1 and obtain the marginal log-likelihood for ξ, θ as follows: 0. For model 2.1 and 2.2 , from the log-quasilikelihood expansion 3.2 , some calculations give the efficient score U θ ξ 0 U θ 1 ξ 0 , . . ., U θ s ξ 0 under the null hypothesis H 0 : θ θ 0 as follows: where e is the magnitude of e Y − f evaluated at ξ, θ ξ 0 , θ 0 , Ḋk 0 ∂D/∂θ k | θ θ 0 k 1, . . ., s and ξ 0 is the maximum likelihood estimator of ξ under θ θ 0 .
To test for H 0 : θ θ 0 , we construct a score statistics as follows: where I θ ξ 0 I θθ −I θξ I −1 ξξ I ξθ is the efficient information matrix of θ evaluated under H 0 .Here, where l l ξ, θ , and the scores ∂l/∂ξ, ∂l/∂θ and the expectations are all calculated at θ θ 0 .Noting the properties of the normal distribution, after some computations, we obtain, for h, k 1, . . ., s, u, v 1, . . ., r, 3.6
One important feature of the proposed score test statistic ST is that a detailed specification of the distribution F b, θ of the random effects is not necessary.Therefore, the test is robust against arbitrary mixed model alternation where only the first two moments of the random effects are specified.The following gives the asymptotic properties of the proposed score test statistic.The "asymptotic" in the theorem refers to the number of clusters m → ∞ with cluster sizes n i bounded in model 2. 1  where l l α, η and the score ∂l/∂α, ∂l/∂η and the expectations are all evaluated at η η 0 , then, similar to 3.4 , one constructs a score statistic as follows: where αα I αη is the efficient information matrix of η.Similar to the derivation of ST , the expression form of the score statistics ST 1 can be obtained.Therefore, for the sake of space, the detailed derivation of ST 1 is omitted.The asymptotic properties of the score test statistic ST 1 are given as follows.
Corollary 3.2.For model 2.1 and 2.2 , under some regularity conditions, when H 0 in 2.4 is true, the asymptotic distribution of the score test statistic ST 1 is a χ 2 -distribution with s r degrees of freedom.Remark 3.3.It is obvious that obtaining MLEs of parameter under null hypothesis is very crucial in the score test.To obtain the MLEs of the unknown parameters in the hierarchical nonlinear mixed-effects models 2.1 and 2.2 , we can take the Newton-Raphson method or the method of score according to the score functions ∂l/∂β, ∂l/∂σ 2 and ∂l/∂γ k see Appendix 4.2 .
Denote by ξ β, σ 2 , γ T the unknown parameter vector, and let l • be the likelihood function defined in 3.2 under the null hypothesis θ θ 0 and l • , l • are the score function and Hessian matrix, respectively.If ξ is the value of parameter vector ξ at the previous iteration, the New-Raphson methods gets the new estimates ξ such that ξ ξ − l ξ −1 l ξ , 3.12 and the methods of score get

3.13
Applying 3.12 or 3.13 iteratively, we can obtain the MLEs ξ.For the asymptotic distribution of the estimate ξ of ξ, we can use the Taylor expansion to show that where ξ 0 is the true value of ξ under the null hypothesis and I −1 • is the Fisher information matrix.
It should be pointed out that the estimating parameters in the nonlinear model are a challenge.In general, the iterating procedure of obtaining the estimates of unknown parameters in a nonlinear model is a popular way.For example, Wong et al. 19 used the Newton-Raphson iteration approach to get the maximum likelihood estimation of ARMA model with error process for replicated observation.In this paper, we adopt a similar method to obtain the MLEs of hierarchical nonlinear mixed-effects model.We do not show in detail the convergence of the algorithm and only give the idea of proof, but we explore the behavior of this algorithm both in the simulation study and the analysis of a real data set in Section 4 and find that the convergence of this algorithm is guaranteed and the precision of convergence is well.The practical results show that the algorithm is reasonable and feasible.

Simulation Studies
We perform some simulation studies to evaluate the sizes and the powers of the score tests proposed in Section 3. We first draw data for m subjects with n i 10 measurements on each unit from the model where x ij are independent random variates from N 0, 1/4 , the ε ij are random noise having the N 0, σ 2 m ij , and m ij exp z ij γ with z ij as independent random drawing from N 0, 1/4 , and the random effects b i i 1, . . ., m are independent random from distribution N 2 0, θI 2 , where I 2 is an 2 × 2 identity matrix.We vary θ from 0 to 0.2, to study the sizes and powers of the score test ST and ST1.
We consider four different sample sizes, m 10, 20, 40, 60.The experiment is replicated 2000 times for each parameter configuration.The nominal sizes of the tests are set to be 0.05.The empirical sizes and powers of test statistics ST and ST1 are presented in Tables 1 and 2 respectively.
The results in Tables 1 and 2 show that the empirical sizes of the tests are very close to 0.05.As the θ, γ, and m increase, the power of the test increases quickly and approaches 1.Furthermore, we notice that, for the score test ST1, when θ 0 but γ / 0, the test has lower power.We speculate that the higher power is obtained because the discrepancy between the model and the postulated model is larger when there exist the random effects in the model; otherwise, the lower power is obtained.These findings are consistent with the theoretical results and also show that random effects in the models may be main factors influencing the inference performance.To further confirm the performance of the test ST and ST1, we present the Q-Q plots of ST and ST1 with m 20, 40 in Figures 1 and 2, respectively.The others are omitted for the sake of space.The Q-Q plots also confirm that the asymptotic distribution of the test ST and ST1 is χ 2 -distribution, which is consistent with the theoretical results.

An Illustrative Example
We illustrate the use of the test in analysis of longitudinal study.The data of the example are taken from the guinea pig data in Johansen 20 .In the experiment, 50 tissue samples were taken from the intestine of each of eight guinea pigs.For each guinea pig, five tissue samples were assigned randomly to each different concentration of B-methyl-glucoside and the uptake volume was measured in micromoles per milligram of fresh tissue per 2 min.Only the means of the five tissue samples at each concentration for each animal were used.The data is previously analyzed by Lee and Xu 21 , to investigate the diagnostic measures through  104.1653P 0.0000 .The degree of freedom of ST is 3.We also calculate the test statistic ST1 and obtain β 1 0.2023 0.0041 , β 2 2.4175 0.0795 , β 3 0.0043 0.0002 , σ 2 0.0017, and ST1 91.4559 P 0.0000 under θ 0, γ 0. The degree of freedom of ST1 is 4. The small P -values of the test statistics suggest strongly rejecting the null hypothesis of independence and homogeneity θ 0 and γ 0. These results demonstrate that there exists heterogeneity of the between-subject variance and correlation within clusters at the same time for the guinea pig data.It should be pointed out that our results are consistent with those of Lee and Xu 21 .They also illustrated the dependence and heterogeneity for the guinea pig data.These results may suggest that the original model having both correlation within clusters and the heterogeneity across different clusters for the guinea pig data should be taken.The illustrative example also demonstrates that the score test can efficiently detect the random effects among the outcomes in practice.
It is worth noticing that the estimate of the random effects is an interesting topic when the existence of random effects is proved.However, in this paper, our interests focus on testing whether the random effects exist.So the estimate of random effects is out of our study scope.Lee and Xu 21 gave an estimate of θ 1 , θ 2 , θ 3 through SA-MCMC algorithm.

Appendices Appendix A
In what follows, we assume that the expectations are taken under θ θ 0 .According to 3.3 , the h, k th component h, k 1, . . .., s of I θθ is, where G, H, Ḋi 0 , Q j , Ω j , i 1, . . ., s, j 1, . . ., m, all being the same as that of in Section 3, respectively.Note the following properties of the normal distribution: if ε ∼ N 0, Σ , then for any nonnegative definite matrices G and H.After some calculations, we obtain A.
and U θ is given in 3.6 .Assume ξ 0 is the true value of ξ.For obtaining the asymptotic properties of the score test statistic ST.We assume the following regularity conditions under θ θ 0 .These assumptions are similar to those given in 13 by Lin.
Condition 1.The size of cluster is a finite sequence of positive integers, the first-and secondorder partial derivatives of f

Advances in Decision Sciences
Condition 3. The log-quasilikelihood l ξ, θ 0 of ξ β T , σ 2 , γ T T has the usual asymptotic properties, including consistency of ξ and the linear expansion Proof.For any given p 1 r s × 1 constant vector λ λ T 1 , λ 2 , λ T 3 , λ T 4 T , where λ 1 is an p × 1 vector, λ 2 is a constant, λ 3 is an r × 1 vector, and λ 4 is an s × 1 vector, we have

B.5
Because Ḋk 0 is a block diagonal matrix, according to the definition of m-dependent sequence, then λ T U is the summation of a sequence of m-dependent random variables; that is, λ T U can be written as where {U λ,i } is an m-dependent sequence.From Conditions 1 and 2, we have that {U λ,i } are uniformly bounded in O ξ 0 for any given λ.It can be shown that E U λ,i 0. By applying Theorem 7. converges in distribution to a chi-square distribution with s degrees of freedom as n → ∞, that is, m → ∞.

Table 1 :
Empirical sizes and powers of the score test statistics ST in 2000 replications.

Table 2 :
Empirical sizes and powers of the score test statistics ST1 in 2000 replications.
Proof of the Asymptotic Distribution of ST.Here we study the asymptotic distribution of ST under H 0 : θ θ 0 .Let •, • with respect to parameter are bounded.The components of Ḋk 0 , k 1, . . ., s, are uniformly bounded.