再探討參考曲線相等性之假設檢定

全文

(1)國立交通大學統計學研究所碩士論文. 再探討參考曲線相等性之假設檢定. Re-visit the Hypothesis Testing for Equality of Reference Charts. 研究生：陳怡頻指導教授：陳鄰安. 博士. 中華民國九十九年十二月.

(2) 再探討參考曲線相等性之假設檢定 Re-visit the Hypothesis Testing for Equality of Reference Charts 研究生：陳怡頻指導教授：陳鄰安. Student：Yi-Pin Chen 博士. Advisor：Dr. Lin-An Chen. 國立交通大學統計學研究所碩士論文. A Thesis Submitted to Institute of Statistics College of Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master in Statistics December 2010 Hsinchu, Taiwan, Republic of China. 中華民國九十九年十二月.

(3) 再探討參考曲線相等性之假設檢定學生：陳怡頻. 指導教授：陳鄰安博士. 國立交通大學理學院統計學研究所. 摘. 要. 探討兩個團體之個體之成長曲線是否相同在文獻上已有些許討論。但，這些討論大都侷限在代表兩個成長模型之迴歸線參數是否相等的問題上打轉。近來，張（2010）提出詳細之迴歸參數關係用以表現兩個成長曲線相等的問題。這個關係清楚說明如何借由參數之檢定達到是否兩團體參數曲線相等之探討。張（2010）也介紹了最大概似比之檢定。這個方法為有母數且較繁複。我們介紹廣義最小平方法之檢定並做了許多模擬分析。. 關鍵字：假設檢定；線性迴歸；檢定力；參考曲線；參考區間；迴歸分位數 i.

(4) Re-visit the Hypothesis Testing for Equality of Reference Charts Student：Yi-Pin Chen. Advisor：Dr. Lin-An Chen. Institute of Statistics National Chiao Tung University. Abstract Comparisons of reference charts for verifying if two populations of subjects have the same growth pattern have received some attention in literature. However, the proposals of comparison are restricted on equalities of regression parameters or regression functions. Recently, Zhang ( 2010 ) has a detailed description for relationships between model parameters and equalities of reference charts that provides a precise indication for verifying if reference charts in two models are identical. She also proposed the likelihood ratio technique as a test for this verification. We re-visit this problem with the generalized least squares estimation and use it to construct the classical F test. Through simulation study, we see that this proposal is very competitive in reference charts comparison.. Key words: Hypothesis testing; linear regression; power; reference chart; reference interval; regression quantile ii.

(5) 誌. 謝. 時光飛逝，兩年的研究生生涯也即將劃下美麗的句點。所上的每一位老師都很親切，不只是學問高深、教學豐富，對學生也充滿了關心與鼓勵。感謝老師們的指導與照顧，在這溫暖的環境中學習，將成為我最珍貴的經驗與回憶。感謝師長、家人以及同學，讓我在研究生涯中充滿了信心與歡笑，謹將此論文獻給你們，並致上我最誠摯的謝意。. 陳怡頻. 謹誌于. 國立交通大學統計學研究所中華民國九十九年十二月 iii.

(6) Contents 中文摘要………………………………………………...………..……..………i Abstract…………………………………….…………………………....………ii 誌謝………………………………………………...……………...……………iii Contents…………………...…………………………………………………….iv. 1. Introduction………………………………………………………..….......…..1 2. Characterization of Reference Charts…………………………..………....…..4 3. Comparison of Two Unknown Reference Charts………………………..........7 4. Reference Charts Comparison for Random Intercept Model…………….….12 5. Reference Charts Comparisons for Random Intercept with Covariate Model………………………..…………………………………………….....20 References……………………………………………………………….…..…26. iv.

(7) Re-visit the Hypothesis Testing for Equality of Reference Charts Abstract. Comparisons of reference charts for verifying if two populations of subjects have the same growth pattern have received some attention in literature. However, the proposals of comparison are restricted on equalities of regression parameters or regression functions. Recently, Zhang (2010) has a detailed description for relationships between model parameters and equalities of reference charts that provides a precise indication for verifying if reference charts in two models are identical. She also proposed the likelihood ratio technique as a test for this veri cation. We re-visit this problem with the generalized least squares etsimation and use it to construct the classical F test. Through simulation study, we see that this proposal is very competitive in reference charts comparison. Key words: Hypothesis testing

(8) linear regression

(9) power

(10) reference chart

(11) reference interval

(12) regression quantile.. 1. Introduction. Growth is a fundamental property of biological systems, occurring at the level of populations, individual animals and plants, as well as within organisms while the growth of a subject depends on nutritional, health, and environmental conditions. Typically the growth pattern for a population group depicts a family of symmetric quantile curves, called reference charts, as a function of some covariates (age or time). One di culty in reference charts problem is that the measurement variables taken over time are generally not independent. Much research has been devoted to modelling growth function and constructing reference charts in parametric or nonparametric way. For overview of parametric methodology, linear or nonlinear growth models, see Cole and Green (1992) and Laird and Ware (1982). When the measurements can be formulated as parametric regression model, the reference charts may be expressed as simple functions of parameters involved in the regression model Typeset by. 1. AMS-TEX.

(13) 2. so that its estimation may be done through estimations of these parameters. For example, the reference charts of a regression with normal errors model are linear functions of the mean and standard deviation. For growth characteristics that are approximately normal, proposals are available for transformations to normal where, among them, the most successful proposal is the LMS by Cole (1988). However, the Exponential-Normal distribution method by Wright and Royston (1997) has the advantage of being parametric with explicit expressions for estimating parameters and quantiles. Verifying the similarity of two growth patterns through comparing the reference charts is an important topic in application. Basically the use of reference charts tries to summarize individual dierences in the growth pattern and it is commonly known that the comparison of reference charts is done by studying the determinants of these dierences. The most common method of comparison considers parametric growth model that the determinals of growth pattern can be represented by a few model parameters so that the job can be done by comparison of these parameters. However, the reference charts comparison considered in literature mainly restrcited on the comparison of growth regression functions. For example, it is seen that most parametric comparison methods consider only those parameters involved in regression function such as testing equality of two or several regression parameter vectors (see, Hoel (1964), Chi and Weerahandi (1998) and Pan and Cole (2004)) or comparing relations between regression slope parameters and (or) intercept parameters (see Zucker, Zerbe and Wu (1995)). Instead of parametric reference charts comparison, there are nonparametric methods comparing the unknown regression functions (see, for examples, Scheike and Zhang (1998), Scheike, Zhang and Juul (1999), Richard, et al. (1989) and Gri ths, Iles et al. (2004)). Hoel (1964) showed that such methods are less e cient than those to compare values of regression parameters. For any comparison exercise, there needs to clarity its precise objectives. Assessment of growth pattern by charts is the most popular tool for de ning health and nutritional status in both individual and population (country) level. Hence, there needs more general study for public health purpose in.

(14) 3. verifying if two or several populations display in the same or similar growth pattern. This is an objective important to be answered in public health, espectially, for studying the developing countries. However, little research has been performed in reference charts comparison in this purpose. It can be seen that comparisons of mean regression functions or few regression parameters can not achieve this public health problem (see Henry (1992)). One exception with a closely related study is that Heckman and Zamar (2000) discussed the concepts of similarity and grouping in growth pattern based on rank correlation coe cient between regression functions. However, besides this is an estimation procedure that it is di cult to extend to hypothesis testing of comparison, regression function comparison is not enough to interprete the similarity or equality of growth patterns characterized by the reference charts. Recently, Zhang (2010) developed explicit relationships between two population reference charts. This developed the analytic relationships between model parameters of growth models that achieved equality of population reference charts. This relationships provides exact test for comparison of reference charts and this observation showed that testing equalities of regression parameters or regression mean functions often provides only a crude approximation to reality so that the conclusions for growth pattern comparison are very questionable. This approach by Zhang is heading in a right direction in a general investigation if two growth models are with the same growth pattern. We concerned that Zhang's likelihood ratio test is too complicated in parameters estimation. The main interest of this paper is to develop generalized least squares estimates for parameters estimation and study the F test for this comparison and study its e ciencies in all kinds of situations of reference charts between two models. In Section 2, we develop conditions on parameters between models to guarantee equality of reference charts and apply this concept to build up test for hypothesis of equal reference charts when the response variable follows a simple linear regression model. In Section 3, we select one interesting longitudinal linear model- the random intercept model - as example to dis-.

(15) 4. play these relations. Tests for this random intercept model are developed for models of without and with covariate, respectively, in Sections 3 and 4. All tests are evaluated with simulations.. 2. Characterization of Reference Charts. We consider that the response variable has a population type linear regression model y(t) = x(t) y + y (t) t 2 S = (0 1) (2.1) 0. where x(t) is vector of independent variables indexed in t, and y (t) is error variable with mean zero and S is the set of age. The th reference chart is the plot of the function Fy 1( jt) against t in S that can be reprented as ;. Cy ( ) = fFy 1 ( jt) : t 2 S g ;. where Fy 1 ( jt) is the conditional quantile of y given age t. The class of reference charts for a population of variable y is ;. fCy ( ) : 2 (0 1)g:. Suppose that there is another population of subjects with a response variable z(t) that follows the same linear regression model with possible dierent parameters as. z(t) = x(t) z + z (t) 0. (2.2). where z (t) is also error variable independent of y (t) with mean zero. Using the same explanatory variables x(t) indicates the balanced design that all the subjects in two groups are measured on the same set of time points. This design is for simplicity of discussion while the theory and method developed in this paper are valid for the unbalance design. For response variable z(t), the th reference chart may be analogously represented as. Cz ( ) = fFz 1 ( jt) : t 2 S g ;. where Fz 1 ( jt) is the th quantile of z at time t and the class of reference charts for the population of variable z is fCz ( ) : 2 (0 1)g. ;.

(16) 5. The interest of the comparison of reference charts is that the two sets of reference charts, respectively, constructed by these two regression models are identical. The parametric approaches of reference charts comparison consider to test equality of regression parameters as. H : y = z :. (2.3). However, the general hypothesis for comparison of reference charts from our formulation then is. H0 : Cy ( ) = Cz ( ) 2 (0 1). (2.4). that generally is not the problem of (2.3). With linear model assumption (2.1), it is seen that the th reference chart may be writtenas Fy 1( jt) = x(t) y + Fy 1 ( ) = x(t) y where 1 y = y + F0y ( ) is called the regression quantile (see Koenker and p 1 Bassett (1978)). The 100 %th reference chart then is ;. 0. ;. 0. ;. ;. Cy ( ) = fx(t) y : t 2 S g: 0. This reects the obaservation by Hoel (1964) that the estimation of reference charts is reduced to estimating the regression quantile . The th regression quantile for model (2.2) is Fz 1 ( jt) = x(t) z + 1 Fz 1 ( ) = x(t) z with z = z + F0z ( ) . Then the th reference p 1 chart for response variable z is ;. 0. ;. ;. 0. ;. Cz ( ) = fx(t) z : t 2 S g: 0. It is agreed, as investigated by Hoel (1964), that comparison of reference charts is more e cient conducted by comparing model parameters. It is then desired to verify when equality of reference charts in hypothesis (2.4) can be re-written into equations in terms of model parameters. The following theorem was provided by Zhang (2010) that gives a rule for testing hypothesis of equal reference charts..

(17) 6. Theorem 2.1. (a) The hypothesis of equal reference charts may be formulated as. Href : y = z Fy 1( ) = Fz 1 ( ) 2 (0 1): (b) If we further assume that Fy 1 ( ) = y F0 1 ( ) and Fz 1 ( ) = z F0 1( ) where y and z are two unknown constants not dependent of time t. Then the hypothesis reduces to ;. ;. ;. ;. ;. ;. Href : y = z y = z :. (2.5). Result of (b) in Theorem 2.1 tells us that solving a comparison of reference charts is valid to be treated as a problem of testing hypothesis for equalities of some model parameters. However, dierent growth models lead to varying hypothesis testing problems. We know that each individual (subject) practically is repeated measured with n-observations y1 ::: yn and x1 ::: xn available from model (2.1). Let us de ne vectors y = (y1 ::: yn) , X = (x1 ::: xn) and y = (y (t1) ::: y (tn )). A matrix form of this regression model for this individual is 0. 0. 0. y = Xy + y. (2.6). where we consider that y (ti ) s are not independent with means 0's and y has covariance matrix as y . Suppose that each individual from another population also has n;observations from model (2.2) so that a matrix form of observations for one individual is 0. z = Xz + z. (2.7). where we consider that z (ti ) s are not independent with means 0 s and z has covariance matrix as z . The di culty in reference chart problem of estimation and hypothesis testing is that the measurement variables taken over time are not statistically independent. Hence, generally the variables in fy (t) : t 2 S g have a complicated structure including correlation. In this consideration, we may test equalities of all model parameters as 0. H : y = z y = z :. 0. (2.8).

(18) 7. 3. Comparison of Two Unknown Reference Charts. When the reference chart is used for public health purposes, it is to compare general health and nutrition of two or more populations (developed and developing world). In this situation, exact test for reference charts comparison for populations is desired to be proposed and evaluated. For statistical inferences, we have m individuals and there are n observations for each individual. For j th individual, there are yj and j follow model of (2.6) as yj = Xy + yj for j = 1 ::: m. By setting vertical joinings y with y = (y1 y2 ::: ym) and y with y = (y1 ::: ym) , vector y has linear regression model of matrix form as 0. 0. 0. 0. 0. 0. 0. y = (1m X )y + y E (y ) = 1m 0n cov(y ) = Im y. (3.1). where represents the Kronecker product, 1m is m-vector of values 1 s and Im is m m identity matrix. Models of this type is interesting since the covariance matrices for various subjects are identical. Suppose that for reference population there are k subjects and n-obsevations follow model (2.7) as zj = Xz + zj j = 1 ::: k 0. where zj s are iid with mean 0n and common covariance matrix z . By setting vertical joinings z with z = (z1 z2 ::: zk ) and with z = (z1 ::: zk ) , vector z has linear regression model 0. 0. 0. 0. 0. 0. z = (1k X )z + z E (z ) = 1k 0n cov(z ) = Ik z :. 0. 0. (3.2). We here consider in this section the classical simple regression with y = y2In and z = z2In . This is the linear regression with repeated measures. First we consider hypothesis H0 : y = z with assumption that y = z . The LSEs of y and z are, respectively,. ^y = m1 (X X ) 0. ;. 1. m X j =1. X yj and ^z = k1 (X X ) 0. 0. 1. ;. k X j =1. X zj : 0. (3.3).

(19) 8. With normality assumption, the F statistic for this hypothesis is. Fbeta =. 0:5( m1 + k1 ) 1 (^y ; ^z ) X X (^y ; ^z ) P (y ; X ^ ) (y ; X ^ ) + Pk (z ; X ^ ) (z ; X ^ )] (m + k)(n ; 2)] 1 m y j y z j z j =1 j j =1 j ;. ;. 0. 0. 0. 0. which follows the distribution f (2 (m + k)(n ; 2)) when H0 is true and the rule for testing H0 is rejecting H0 if Fbeta > f(2 (m + k)(n ; 2)):. (3.4). For all experiments, we conduct a simulation with replications 10 000 under some settings of numbers of subjects m and k ans sample size n. The parameter values under H0 is y = z = (1 1) and y = z = 1. Table 1 displays the resulted sizes of the test. 0. Table 1. Type I error probabilities unedr F test m=k=5 m = k = 10 m = k = 15 m = k = 20 m = k = 25 m = k = 30. n = 50. n = 100. n = 200. n = 1 000. 0:0452 0:0477 0:0442 0:0425 0:0481 0:0416. 0:0466 0:0484 0:0465 0:0461 0:0456 0:0486. 0:0504 0:0514 0:0471 0:0516 0:0496 0:044. 0:0484 0:0511 0:0504 0:0496 0:05 0:0508. The results showing in Table 1 indicates that this test is approximately with sizes very close to the signi cance level = 0:05. We then further conduct simulation with various design of deviation from H0 to verify the power performance. The simulated results are displayed in Table 2.. Table 2. Power performance unedr F test for hypothesis of equal regression parameters.

(20) 9. Parameters y2 = z2 = 1 0y = 1:14 0y = 1:16 0y = 1:2 1y = 1:008 1y = 1:01 1y = 1:02. power. y2 = z2 = 2:5 0y = 1:25 0y = 1:3 0y = 1:4 1y = 1:014 1y = 1:03. 0:7351 0:8503 0:9679 0:7524 0:9147 1. Parameters y2 = z2 = 5 0y = 1:32 0y = 1:5 0y = 1:65 1y = 1:018 1y = 1:025 1y = 1:038. power 0:7508 0:9887 1 0:7551 0:9639 1. 0:8423 0:9503 0:998 0:8375 1. y2 = z2 = 10 0y = 1:5 0y = 1:6 0y = 1:7 1y = 1:03 1y = 1:05. 0:8374 0:9481 0:9894 0:8852 0:9999. We have two comments drawn from the results showing in Table 2: 1. The power is increasing when the parameter 0y or 1y moves away from the value under H0 . 2. The power is very sensitive to the change in slope parameter 1y and it is relatively less sensitive for a departure in location parameter 0y . However, the performance seems to be satisfactory. We next consider a test for hypothesis of dierences in variances as H0 : y = z vs H1 : y = 6 z . We consider the following test statistic. Pm (y ; X ^ ) (y ; X ^ )=(mn ; 2) y j y Fsigma = Pjk=1 j (z ; X ^ ) (z ; X ^ )=(kn ; 2) 0. j =1 j. z. 0. j. z. and the rule for testing H0 is rejecting H0 if Fsigma < f=2(mn;2 kn;2) or Fsigma > f1 =2(mn;2 kn;2): ;. We rst, again, verify the size of this test.. Table 3. Type I error probabilities unedr F test for hypothesis of equal variances.

(21) 10. n = 20. n = 30. n = 50. = 0:05 m = k = 10 0:0487 0:0477 0:0491 m = k = 20 0:0496 0:0522 0:0522 m = k = 30 0:05 0:0492 0:053 = 0:1 m = k = 10 0:1003 0:1005 0:1028 m = k = 20 0:0987 0:1045 0:0947 m = k = 30 0:0994 0:101 0:1001 = 0:15 m = k = 10 0:1454 0:1499 0:1568 m = k = 20 0:1513 0:1457 0:1451 m = k = 30 0:1461 0:1516 0:1437 The sizes in most cases are very close to the set signi cance levels. Hence, this test is approximately level test. For power comparison, we consider the following design in cases that H0 is not true: y = (1 1) z = (a b) y2 = 1 z2 = c: 0. This allows us to verify both regression parameters and error variances are varied in two growth models.. Table 4. Powers under F test for hypopthesis of equal variances (a b c). n = 20. n = 30. n = 50. (1 1 1:2) 0:6056 0:7818 0:9401 (1 1 1:3) 0:8936 0:9742 0:9994 (1 1 1:5) 0:9983 1 1 (1 1 2) 1 1 1 (2 2 1:2) 0:6022 0:7772 0:9402 (2 2 1:4) 0:9857 0:9991 1 (2 2 2) 1 1 1 We have two comments drawn from the results showing in Table 4: 1. In the case that (a b) = (1 1), regression parameters in two growth models are set equal. The power performance is showing satisfactory. This indicates that a test for equality of regression parameters is not su cient for comparison of reference charts..

(22) 11. 2. In cases with (a b) = (2 2), the regression parameters and error variances are both unequal. It also shows satisfactory in power performance. We consider the hypothesis. H0 : y = z y = z :. p. (3.5). By Bonferroni, by simulation, we choose = 1 ; 0:9434 and set the rule as rejecting H0 if Fbeta > f1 (2 (m + k)(n ; 2)) or Fsigma f=2(mn ; 2 kn ; 2) or f1 =2(mn ; 2 kn ; 2): (3.6) ;. ;. The simulated probability of type I error is 0:049, very close to the aquired signi cance level 0:05.. Table 5. Power comparison when regression parameters and error variances are both changed ( = 0:05) parameters power (0y z2) (1:1 1:1) (1:2 1:1) (1:3 1:1) (1:1 1:2) (1:1 1:4). 0:4705 0:9486 0:9997 0:802 0:9992. (0y z2 1z ) (1:1 1:1 1:02) (1:3 1:1 1:01). 0:9987 0:8623. parameters. power. (0y y2 z2) (1:1 1:3 1:1) (1:1 1:4 1:1) (1:1 1:6 1:1). 0:7249 0:9407 0:9998. (1y y2 z2) (1:01 1:2 1:1) (1:01 1:3 1:1) (1:01 1:4 1:1) (1:02 1:1 1:1). 0:8515 0:9257 0:9842 1. We have several comments for the power simulation results: (a) Consider the departures in 0y and z2. The power is low when two parameters are changed slightly as (1:1 1:1). However, it rapidly increse when there departures slightly bigger in either one. (b) The case (0y z2 1z ) = (1:1 1:1 1:02) has power 0:9987, by comparison with case for (0y z2) = (1:1 1:1), showing that a very small shift in regression slope parameter 1z makes power greatly improved..

(23) 12. (c) The settings for (0y y2 z2) x changes in 0y and z2 and then the power performance corresponds to the change in y2 is mild but satisfactory. (d) The settings for (1y y2 z2) also results that power is much more sensitive for a change in slope parameter 1y than changes in intercept or error variances.. Table 6. Power comparison when regression parameters and error variances are both changed ( = 0:1) parameters power (0y z2) (1:1 1:1) (1:2 1:1) (1:3 1:1) (1:1 1:2) (1:1 1:4). 0:6193 0:9771 1 0:8874 0:9994. parameters. power. (0y y2 1z ) (1:1 1:1 1:01) (1:3 1:1 1:01) (1:4 1:1 1:01) (1:1 1:3 1:01) (1:1 1:1 1:02). 0:636 0:9251 0:9988 0:9892 0:9999. (0y y2 z2) (1:1 1:3 1:1) 0:8297 (1:1 1:5 1:1) 0:9979 When the signi cance level is 0:1, the simulated results also show that the power performances for intercept and error variances are not so sensitive as it did for slope parameter. However, these results show that this technique is satisfactory.. 4. Reference Charts Comparison for Random Intercept Model. In this section, we consider one interesting repeated measurements regression model as example to formulate speci cation of hypothesis in (2.5) and develop procedures to test hypothesis of equal reference charts in this model. The random intercept model for one individual is of the form. y(tj ) = 0y + Vy + 1y x1 (tj ) +

(24) y (tj ) j = 1 ::: n 2 ) and

(25) (t ) s are independent where Vy has normal distributions N (0 vy y j 2 normal distributions N (0 y ). Also, variables Vy and

(26) y (tj ) are assumed to 0.

(27) 13. be independent. This is a model of (2.1) with x(t) y = 0y + 1y x1(t) and y (t) = Vy +

(28) y (t). This random intercept model allows each individual to have its own intercept term and then the starting level for this individual is 0y + vy where various subjects may have dierent observations vy s of Vy . This random intercept regression model has the form of (2.6) with designed matrix X and parameter vector y as 0 1 x1 1 1 x2 C B 0y B C X = @ .. (4.1) A y = 1y y = y2In + vy2 J . 1 xn with xj = x(tj ) and J is n n matrix of 1's. The set of reference charts for this random intercept model is 0. 0. q. 2 + 2z : t 2 S g Cy ( ) = f0y + 1y x1 (t) + vy y . (4.2). where z is the th quantile point for the standard normal distribution. Since the covariances are identical, this random intercept model is also called the uniform correlation model. The random intercept model for subject from another population is. z(tj ) = 0z + Vz + 1z x1 (tj ) +

(29) z (tj ) j = 1 ::: n 2 ) and

(30) (t ) s are independent where Vz has normal distributions N (0 vz z j 2 normal distributions N (0 z ). The set of reference charts for this random intercept model is 0. p. 2 + 2 z : t 2 S g: Cz ( ) = f0z + 1z x1 (t) + vz z . The equality of reference charts from these two populations indicates. q. p. 2 + 2 z = + x (t)+ 2 + 2 z for all x (t) and 2 (0 1) 0y +1y x1 (t)+ vy 0z 1z 1 1 y vz z . giving that testing equalities of reference charts is equivalent to test the following hypothesis 2 2 Href : y = z y2 + vy = z2 + vz :. (4.3).

(31) 14. When we test the hypothesis H , acceptance of y = z obviously does not 2 = indicate equality of reference charts since equality of q variance sum y2 +vy p2 indicates 2 may not be guaranteed. Since = and 2 = z2 + vz vy y z vz 2 2 2 2 that y + vy = z + vz is true, so, when we test hypothesis H and the hypothesis is accepted then we are sure that the two reference charts are equal. However, there is a risk that these two reference charts are really 2 = 2 + 2 doesn't indicate that equal when we reject H since y2 + vy z vz y = z and vy = vz are true. When we are allowed to assume that y = z . The hypotsheis is reduced to the following y = z y = z vy = vz (4.4) and then testing hypothesis H is then appropriate. We now restrict on the following no-covariate random intercept model. From (3.1), the random intercept model may be re-written as 2 y = 0y 1mn + y E (y ) = 1m 0n cov(y ) = (y2 + vy )(Im 0y ). where 0y is n n matrix with diagonal elements 1's and o diagonal ele2 ments y2+vyvy 2 . Similarly the random intercept model for z may be re-written as 2 z = 0z 1kn + z E (z ) = 1k 0n cov(z ) = (z2 + vz )(Ik 0z ):. Equality of refernce charts in this no-covariate model requires to test the following hypothesis 2 2 H0 : 0y = 0z y2 + vy = z2 + vz :. (4.5). We rst consider hypothesis. H0 : 0y = 0z : We assume that the response variable for subjects in two groups both have random intercept models that leads us to consider the intercept parameter.

(32) 15. estimators ^0y and ^0z so that their dierence ^0y ; ^0z has (asymptotic) 2 covariance natrix achieving the Cramer-Rao lower bound ma (1n 0y1 1n ) 1 + b2 (1 1 1 ) 1 . Available estimators of and for this kind include 0y 0z k n 0z n mle and generalized least squares estimator. Suppose that we have error variances â2 and ^b2 and covariance matrix estimators ^ 0y and ^ 0z . With normality assumption, the F statistic for this hypothesis is 0. 0. ;. ;. ;. ;. ^0y ; ^0z )2 ( Fbeta = â2 ^b2 + 1 ^ ^ m1n 0y 1n k1n 0z1 1n 0. ;. 0. ;. when H0 is true and the rule for testing H0 is rejecting H0 if Fbeta > 2 (1): First we consider mle method to perform the above test. Let us denote the followings (see the formulas in McCulloch and Searle ( 2001)):. Pn y Pm Pn y Pm y m X n X ij j =1 ij i =1 yi: = n y:j = m y:: = i=1 nmj=1 ij SSTy = (yij ; y::)2 i=1 j =1 m m n X XX 2 2 SSAy =. i=1. n(yi: ; y::) SSEy =. i=1 j =1. (yij ; yi: ). y MSA = SSAy MSE = SSEy MSTy = SST y m;1 y m(n ; 1) mn MSE 1 )MSA MSE if (1 ; y 2 m ^y = MST if (1 ; 1 )MSAy < MSEy y y m (1 m1y)MSAy MSE y if (1 ; 1 )MSA MSE y y n m1 ^v2y = 0 if (1 ; m )MSAy < MSEy m m X X ^0y = (m1n ^ 0y1 1n ) 1 1n ^ 0y1 yj â2 = mn1; 1 (yj ; ^0y 1n ) ^ 0y1 (yj ; ^0y 1n ) j =1 j =1 ;. 0. ;. ;. ;. 0. ;. 0. ;.

(33) 16. and. Pk Pn z Pn z Pk z k X n X ij i =1 j =1 ij j =1 ij i =1 ( z ; z:: )2 SST = zi: = n z:j = k z:: = ij z nm SSAz =. k X i=1. n(zi: ; z:: ) SSEz = 2. k X n X i=1 j =1. i=1 j =1. (zij ; zi: )2. z MSA = SSAz MSE = SSEz MSTz = SST z k;1 z k(n ; 1) kn MSE ; k1 )MSAz MSEz ^z2 = MST z ifif (1 (1 ; k1 )MSAz < MSEz z (1 k1 )MSAz MSEz if (1 ; k1 )MSAz MSEz 2 n ^vz = 0 if (1 ; k1 )MSAz < MSEz k k ^0z = (k1n ^ 0z1 1n ) 1 X 1n ^ 0z1 zj ^b2 = 1 X(zj ; ^0z 1n ) ^ 0z1 (zj ; ^0z 1n ): kn ; 1 j=1 j =1 ;. 0. ;. ;. ;. 0. ;. 0. ;. Based on the above settings, the test for hypothesis H0 : 0y = 0z is 2 rejecting H0 if F = ^â2 f=2(mn ; 1 kn ; 1) or f1 =2(mn ; 1 kn ; 1): b To test if this test is appropriate, we conduct a simulation with m = k, n = 30 and = 0:05. From now on, each simulation performs 10 000 replications and the parameters not been speci ed are given values 1's. The following we list the simulated probabilities of type I error: m = k = 2 10 15 20 30 ;. ;;. ;;. ;;. ;;. ;;. 0.242 0.057 0.058 0.059 0.051 This test seems to be reasonable for its size in testing dierences of location parameters based on the simulated results when n = k 10. We next conduct a simulation to verufy the power of this test. Table 7. Power performance for hypothesis 0y = 0z based mle with 0z = 1 0y Power 0y Power 0y = 1:5 0:192 0y = 3:0 0:983 2:0 0:514 3:2 0:993 2:5 0:859 3:5 1.

(34) 17. The power gradually increases when 0y is moving away from 0z . The results indicates that this test is satisfactory. We now propose a new estimation of these covariance matrices. Model for the j -th subject is. yj = (0y + vy )1n +

(35) y : Instead of estimating regression parameters through whole set of data, we consider combination of estimates from individual regression models. We introduce the idea and process in this estimation in the followings: (a) Suppose that the least squares estimator of this model for j -th subject is b0j . The estimate b0j actually estimate the sum of intercept 0y and variable vy . Since Vy is a variable having zero mean, we may denote ^0y = 1 Pm b as estimate of intercept parameter . 0y m j =1 0j (b) By letting v^yj = b0j ; ^0y j = 1 ::: m, then v^yj 's are appropriate as predictors of intercept error variables Vyj 's. We then have variance estimator P v^2 . 2 ^GLvy = m1 m j =1 yj (c) Since b0j is estimate of vector 0y +vy , the residual yj ;b0j is one predictor of error vector

(36) y . Hence the estimator of variance y2 is appropriate to be P (y ; b 1 ) (y ; b 1 ). The n n matrix estimates 2 ^GLy = m(n1 1) m j =1 j 0j n j 0j n are 0. ;. ^ 2 GLy : diagonal elements 1's, o diagonal elements ^ 2 GLvy 2 GLy + ^GLvy(4.6) and 2 ^GLvz GLz : diagonal elements 1's, o diagonal elements ^ 2 + ^ 2 : GLz GLvz(4.7) The generalized LSEs for regression paramters 0y and 0z and error.

(37) 18. 2 + 2 and 2 = 2 + 2 are, respectively, variances a2 = vy vz z y b. ^GL0y. Pm 1 ^ 1y = j=1 n 0y j ^ ;. 0. Pk 1 ^ 1z = j=1 n 0z j 0. ;. m1n ^ 0y1 1n GL0z k1n ^ 0z1 1n m X 1 1 2 (yj ; ^GL0y 1n ) and ^GLa = mn ; 1 (yj ; ^GL0y 1n ) ^ GLy j =1 0. ;. 0. 0. 1. ^GLb = kn ; 1 2. k X j =1. ;. ;. 1 (zj ; ^GL0z 1n ) ^ GLz (zj ; ^GL0z 1n ): 0. ;. To test if this test is appropriate, we conduct a simulation with m = k, n = 30 and = 0:05. The following we list the simulated probabilities of type I error: m = k = 2 10 15 20 30 ;;. ;;. ;;. ;;. ;;. 0.235 0.052 0.049 0.054 0.053 Besides m = k = 2, this technique has size close to the speci ed signi cance level. We then further study the power performance of this test.. Table 8. Power performance for hypothesis 0y = 0z based generalized. LSE with 0z = 1 0y Power 0y Power 0y = 1:5 0:193 0y = 3:0 0:988 2:0 0:52 3:2 0:996 2:5 0:855 3:5 1 The power performance based on generalized least squares estimator is also acceptable for its monotone power property. We next consider the following hypothesis 2 2 H0 : y2 + vy = z2 + vz :. From normal theory, a F test is de ned as ^ 2 f=2 (mn;1 kn;1) or f1 =2 (mn;1 kn;1): rejecting H0 if Fsigma = ^GLa 2 GLb ;. We now consider the generalized LSE method. We conduct a simulation with m = k, n = 30 and = 0:05. The following we list the simulated.

(38) 19. probabilities of type I error:. m=k=5. 10. ;;. ;;. 0.055. 15. 20. ;;. ;;. 25. ;;. 30. ;;. 0.057 0.062 0.056 0.052 0.059. This is also close to the speci ed signi cance level although they are a bit larger than that.. Table 9. Power performance for hypothesis a2 = b2 based generalized LSE Sample size m=k=2 n = 30 n = 50 n = 100 m = k = 10 n = 30 n = 50 n = 100. 2 = 2 = 1 vz z. 1:5. 2. 3. 0:061 0:059 0:08. 0:367 0:52 0:768. 0:739 0:921 0:997. 0:989 1 1. 0:057 0:069 0:085. 0:923 0:991 1. 1 1 1. 1 1 1. The power is gradually increasing, although uctuated in several cases, in 2 and 2 . The test for equality of reference charts of sample size n and vz z two variance sum is appropriate. It is then considered the combination of the above two tests for hypothesis (3.5). By choosing signi cance level 0:05, we set = 0:05 with the following test: rejecting H0 if F =. (^0y ; ^0z )2 2 (1) ^b2 â2 m1n ^ 0y1 1n + k1n ^ 0z1 1n 0. ;. 0. ;. ^ 2 f=2 (mn ; 1 kn ; 1) or f1 =2 (mn ; 1 kn ; 1): or Fsigma = ^GLa 2 GLb ;. The powers are listed in Table 10.. Table 10. Powers for simultaneous F test.

(39) 20. power. power. 0y vy 0z 0y vy z 2:5 1:1 1:1 0:788 2:5 1:1 1:1 0:837 2:7 1:1 1:1 0:877 2:7 1:1 1:1 0:909 3 1:1 1:1 0:974 3 1:1 1:1 0:98 1:1 2 1:1 0:829 1:1 2:2 1:1 0:853 1:1 2:2 1:1 0:93 1:1 2:4 1:1 0:927 1:1 2:5 1:1 0:99 1:1 2:8 1:1 0:987 We set two cases in departure of null hypothesis H0. The simulated results indicate that the test for equality of reference charts is appropriate by the generalized least squares method. 2. 2. 2. 5. Reference Charts Comparisons for Random Intercept with Covariate Model. In this section, we consider the random intercept with covariate model. From (3.1), this model may be re-written as. y = (1m X )y + y E (y ) = 1m 0n cov(y ) = a2(Im 0y ). (5.1). 2 , is n n matrix with diagonal elements 1's and o where a2 = y2 + vy 0y 2 vy diagonal elements y2 +vy 2 . Similarly the random intercept model for z may be re-written as. z = (1k X )z + z E (z ) = 1k 0n :cov(z ) = b2(Ik 0z ). (5.2). 2 and is n n matrix with diagonal elements 1's and where b2 = z2 + vz 0z 2 o diagonal elements z2 +vzvzy 2 . We rst consider the following hypothesis. H0 : y = z :. (5.3). We assume that the response variable for subjects in two groups both have random intercept models that leads us to consider the regression parameter estimators ^y and ^z so that their dierence ^y ; ^z has (asymptotic) 2 covariance natrix achieving the Cramer-Rao lower bound ma (X 0y1 X ) 1 + b2 (X 1 X ) 1 . Available estimators of and for this kind include y z 0z k mle and generalized least squares estimator. Suppose that we have error 0. 0. ;. ;. ;. ;.

(40) 21. variances â2 and ^b2 and covariance matrix estimators ^ 0y and ^ 0z . With normality assumption, the F statistic for this hypothesis is 2 2 Fbeta = (^y ; ^z ) ( ^ma (X ^ 0y1 X ) 1 + ^kb (X ^ 0z1X ) 1 ) 1(^y ; ^z ) when H0 is true and the rule for testing H0 is 0. 0. ;. ;. 0. ;. rejecting H0 if Fbeta > 2 (2):. ;. ;. (5.4). We now propose a new estimation of these covariance matrices. Model for the j -th subject is. . yj = X 0y+ vy +

(41) y : (5.5) 1 Again, we consider combination of estimates from individual regression models. We introduce the idea and process in this estimation in the followings: (a) Suppose that the least squares estimator of this model for j -th subject is ^j = bb0j . The estimate b0j actually estimate the sum of regression 1j intercept 0y and variable vy . Since Vy is a variable having zero mean, we P b as estimate of regression intercept parameter may denote ^0 = m1 m j =1 0j 0y . (b) By letting v^yj = b0j ; ^0 j = 1 ::: m, then v^yj 's are appropriate as predictions of intercept error variables Vyj 's. We then have variance estimator P v^2 . 2 ^GLvy = m1 m j =1 yj + v 0 y y (c) Since ^j is estimate of vector , the residual yj ; X ^j is 1 one predictor of error vector

(42) y . Hence the estimator of variance y2 is P (y ; X ^ ) (y ; X ^ ). The n n 2 appropriate to be ^GLy = m(n1 2) m j j j j =1 j matrix estimates are 2 ^GLvy GLy : diagonal elements 1's, o diagonal elements ^ 2 + ^ 2 GLy GLvy(5.6) and 2 ^GLvz GLz : diagonal elements 1's, o diagonal elements ^ 2 + ^ 2 : GLz GLvz(5.7) 0. ;.

(43) 22. The generalized LSEs for regression paramters y and z and error vari2 + 2 and 2 = 2 + 2 are, respectively, ances a2 = vy y vz z b. ^GLy = m1 (X ^ 0y1 X ) ;. 0. 2 ^GLa = mn1; 2 2 ^GLb = kn1; 2. 1. ;. m X. j =1 k X j =1. m k X X 1 1 1 1 ^ ^ ^ X 0y yj GLz = k (X 0z X ) X ^ 0z1 zj 0. ;. 0. ;. 0. ;. j =1. ;. j =1. 1 (yj ; X ^GLy ) ^ GLy (yj ; X ^GLy ) and 0. ;. 1 (zj ; X ^GLz ) ^ GLz (zj ; X ^GLz ): 0. ;. To verify if the test of (5.4) based on generalized least squares estimators is appropriate, we conduct a simulation with m = k and = 0:05. The following we list the simulated probabilities of type I error: m = k = 2 10 15 20 30 ;;. ;;. ;;. ;;. ;;. 0:2012 0:0557 0:0544 0:0522 0:0494 Again, m = k = 2 has error probabilities too large. The others are satisfactory. To evaluate the power performance, we design shifts in some regression parameters and we list the powers in the following table.. Table 11. Powers for testing H0 : y = z with generalized LSE (n = 30) Parameter shift k = m = 10 0y = 1:5 2:0 2:5 2:8 3 3:2. Power. Parameter shift. Power. 0:149 0:400 0:759 0:8989 0:952 0:9806. 0y = 1:2 1y = 1:02 0y = 1:2 1y = 1:05 0y = 1:5 1y = 1:02 0y = 1:5 1y = 1:05 0y = 2:0 1y = 1:02 0y = 2:0 1y = 1:05. 0:48 1 0:618 0:998 0:834 1. 1y = 0:96 0:97 0:98 1:02 1:04 1:05. 0:973 0:788 0:407 0:405 0:982 0:999.

(44) 23. The powers resulted from simulation seem to be very strong when there is shift in slope parameter 1y . This ful lls our expectation in regression model. The second hypothesis of interest is 2 2 H0 : vy + y2 = vz + z2:. Suppose that we have regression parameter estimators ^y and ^z and estima2 tors of covariance matrices ^ 0y and ^ 0z . We have estimates of a2 = y2 + vy 2 as and b2 = z2 + vz m X 1 â = mn ; 2 (yj ; X ^y ) ^ 0y1 (yj ; X ^y ) and j =1 2. 1. ^b = kn ; 2 2. ;. 0. k X j =1. (zj ; X ^z ) ^ 0z1 (zj ; X ^z ): 0. ;. From normal theory, a F test is de ned as 2 rejecting H0 if F = ^â2 f=2(mn ; 2 kn ; 2) or f1 =2(mn ; 2 kn ; 2): b ;. To perform this test, it requires estimates of y z and covariance matrices 0y and 0z . Again, we consider mle method and generalized least squares estimation for comparison. We rst consider the mle estimates ^mley ^mlez and ^ mley ^ mlez . Let us denote the followings:. Pn y Pm Pn y Pn x Pm y ij j =1 ij i =1 j =1 ij j =1 j i =1 yi: = n y:j = m y:: = x = nm n n n X X 2 Syx =. (y:j ; y::)(xj ; x) Sxx =. j =1 n X Szx (z:j ; z:: )(xj ; x): j =1. j =1. (xj ; x) . The following maximum likelihood estimates are used to construct the.

(45) 24. estimators of covariance matrices: Pm Pn (y ; (y + ^ (x ; x)))2 S yx 2 ^0y = y:: ; ^1y x ^1y = S ^y = i=1 j=1 ijm(n ;:: 1) 1y j xx P m yi: ; y::)2 ; ^ 2)g 2 ^vy = maxf0 n1 ( n i=1 ( y m Pk Pn (z ; (z + ^ (x ; x)))2 ^0z = z:: ; ^1z x ^1z = SSzx ^z2 = i=1 j=1 ijk(n ;:: 1) 1z j xx P k n ( z ; z:: )2 1 i: 2 i =1 ^vz = maxf0 n ( ; ^z2 )g: k We now state the mle estimates of regression parameters and nn covariance matrices 0y and 0z as ^ ^ 0 y 0 z 2 2 ^mley = ^ ^mlez = ^ â2 = ^y2 + ^vy ^b2 = ^z2 + ^vz 1z 1y (5.8) ^ 2 mley : matrix, diagonal elements 1's, o diagonal elements ^vy2 a (5.9) and 2 mlez : matrix, diagonal elements 1's, o diagonal elements ^^vz2 : b (5.10) The power performance based on mle method is listed in the follwoing table. Table 12. Powers for F test with y2 = vy2 = 1 with mle method 2 = 2 = 1 sample size vz 2 3 5 z k=m=2 n = 30(df = 2 30 ; 28) 0:0517 0:4853 0:8596 0:9947 n = 50(2 50 ; 57) 0:0544 0:6082 0:9398 0:9996 n = 100(2 100 ; 145) 0:0520 0:7136 0:9809 1 k = m = 10 n = 30(10 30 ; 150) 0:05 0:99 1 1 n = 50(10 50 ; 308) 0:051 0:998 1 1 n = 100(10 100 ; 761) 0:05 1 1 1 2 = 2 = 1 are probabilities of type I error. These The powers when vz z values shows that this test is appropriate. Then the powers for departure of null hypothesis H0 are also reasonably good..

(46) 25. The generalized least method apply regression parameter estimates ^GLy and ^GLz and variance estimates ^GLa and ^GLb. To test if this test is appropriate, we conduct a simulation with m = k and = 0:05. The following we list the simulated probabilities of type I error: m = k = 5 10 15 20 25 30 ;;. ;;. ;;. ;;. ;;. ;;. 0:056 0:069 0:049 0:056 0:052 0:06 The type I error probabilities are a bit too large. This needs more studies and modi cations. Table 13. Powers for F test with y2 = vy2 = 1 with generalized least squares method 2 = 2 = 1 sample size vz 1:5 2 3 5 z k=m=2 n = 30 0:062 0:341 0:764 0:986 1 n = 50 0:068 0:514 0:904 0:999 1 n = 100 0:089 0:781 0:994 1 1 k = m = 10 n = 30 0:061 0:944 1 1 1 n = 50 0:056 0:995 1 1 1 n = 100 0:091 1 1 1 1 The power performance as the other cases looks satisfactory. We now consider simultaneous testing for hypothesis 2 2 H0 : y = z and y2 + vy = z2 + vz : The F test for this hypothesis is 2 rejecting H0 if Fbeta > 2 (2) or F = ^â2 f=2 (mn;2 kn;2) or f1 =2(mn;2 kn;2): b Table 14. Powers for simultaneous F test 2 2 2 2 0y vy power 1y vy power vz vz 1:1 2 1:1 0:7 1:03 1:1 1:1 0:795 1:1 2:2 1:1 0:828 1:04 1:1 1:1 0:973 1:1 2:4 1:1 0:928 1:05 1:1 1:1 0:995 1:1 1:1 2 0:683 1:06 1:1 1:1 1 2 2 1:1 1:1 2:3 0:9 0y 1y vy z 1:1 1:1 2:6 0:972 1:1 1:03 1:1 1:1 0:752 1:1 1:04 1:1 1:1 0:966 1:1 1:05 1:1 1:1 0:998 ;.

(47) 26. Table 15. Powers for simultaneous F test 0y y2 z2. 2:5 1:1 1:1 2:7 1:1 1:1 3 1:1 1:1 1:1 2 1:1 1:1 2:2 1:1 1:1 2:5 1:1 The simulated results show power.. power. 0:849 0:938 0:99 0:817 0:926 0:99 that this proposal. References. 0y y2 0z. power. 2:5 1:1 1:1 2:7 1:1 1:1 3 1:1 1:1 1:1 2 1:1 1:1 2:2 1:1. 0:806 0:889 0:978 0:91 0:967. has good performance in. Chi, E. M. and Weerahandi, S. (1998). Comparing treatments under growth curve models: exact tests using generalized p-values. Journal of Statistical Planning and Inference, 71, 179-189. Cole, T. J. (1988). Fitting smoothed centile curves to reference data (with Discussion). Journal of the Royal Statistical Society, Series A, 151, 385-418. Cole, T. J. and Green, P . J. (1992). Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in Medicine, 11, 1305-1319. Gri ths, J. K., Iles, T. C., Koduah, M. and Nix, A. B. J. (2004). Centile charts II: alternative nonparametric approach for establishing timespeci c refrence centiles and assessment of the sample size required. Clinical Chemistry. 50, 907-914. Heckman, N. E. and Zamar, R. H. (2000). Comparing the shapes of regression functions. Biometrika. 87, 135-144. Henry, JJ. (1992). Routine growth monitoring and assessment of growth disorder. Journal of Pediatr Health Care. 6, 291-301. Hoel, P. G. (1964). Methods for comparing growth type curves. Biometrics. Koenker, R. and Bassett, G. J. (1978). Regression quantiles. Econometrica. 46, 33-50..

(48) 27. Laird and Ware (1982). Random-eects models for longitudinal data. Biometrics. 38, 963-974. McCulloch, C. E. and Searle, S. R. (2001). Generalized, Linear, and Mixed Models. Wiley: New York. Pan, H. and Cole, T. J. (2004). A comparison of goodness of t for agerelated reference ranges. Statistics In Medicine. 23, 1749-1765. Richard, R. S., Engeman, R. M., Zerbe, G. O. and Bury, R. B. (1989). A nonparametric comparison of monomolecular growth curves: application to westrn painted turtle data. Growth Development & Aging. 53, 47-56. Scheike, T. H. and Zhang, M.-J. (1998). Cumulative regression function tests for regression models for longitudinal data. The Annals of Statistics. 26, 1328-1355. Scheike, T. H., Zhang, M.-J. and Juul, A. (1999). Comparing reference charts. Biometrical Journal. 6, 679-687. Wright, E. and Royston, P. (1997). Simpli ed estimation of age-speci c reference intervals for skwed data. Statistics in Medicine, 16, 27852803. Zhang, M.-Z. (2010). Hypothesis testing for equality of reference charts. Master thesis: National Chiao Tung University. Zucker, D. M., Zerbe, G. O. and Wu, M. C. (1995). Inference for the association between coe cients in a multivariate growth curve model. Biometrics. 51, 413-424..

(49)