Bias correction in a multivariate normal regression model with general parameterization

Bias correction in a m ultiv ariate normal regression mo del with general parameterization Alexandr e G. P atri ota, Artur J. Lemonte Dep artamento de Estat ´ ıstic a, Universidade de S˜ ao Paulo, R ua do Mat˜ ao, 1010 , S˜ ao Paulo/SP, 055 08-090, Br azil Abstract This pap er dev elops a bias correction sc heme for a multiv ariate normal mo del under a general parameterization. In the mo del, the mean v ector and the cov ariance matrix s h a re the same parameters. It includes many imp or- tan t regression m odels av ailable in the literature as sp ecial cases, suc h as (non)linear regression, errors-in-v ariables mo dels, and so forth. Moreo ver, heteroscedastic situations ma y also b e studied within our framework. W e deriv e a general expression for the second-order biases of maxim um like li- ho od estimates of the mo del p a rameters and show that it is alwa ys p ossible to obtain the second order b ia s by m e ans of ord inary weigh ted lest-squares regressions. W e en lighten s uc h general expr e ssion w i th an errors-in-v ariables mo del and also conduct some simulatio ns in order to v erify the p erformance of the corrected estimates. The s imulation r esults sho w that the bias correc- tion sc heme yields nearly unbiase d estimators. W e also present an emp iric al ilustration. Key W o r ds: Bias c orr e ction, err ors-in-variables mo del, maximum likeliho o d estimation, multivariate r e gr ession. 1 In tro duction Applications of m ultiv ariate normal mo dels are commonly found in the literature. There are simple mo dels that do no t require asymptotic appro ximations for t h e maxim um lik eliho od estimators. Nev ertheless, in the ma j orit y of problems, the es- timation pro cedure in suc h m ultiv ariate normal mo dels embrace on the asymptotic theory . F or instance, it is hard to compute the exact distribution of maxim um lik e- liho o d estimators (MLEs) for nonlinear m ultiv ariate regressions, errors-in-v ariables mo dels and man y others when the sample size is finite (practical situatio ns ). Then, in the practical applications, the asymptotic distribution of the MLE is us ed as an appro ximation to its exact distribution. It considerable simplifies t he inferen- tial pro cess. In general, under some regular it y conditions, the MLEs are consisten t and efficien t, i.e., asymptotically , their biases con v erge to zero and their v ariance- co v ariance matrices approac h the in v erse of the Fisher information. Moreov er, under 1 suc h regular ity conditions, the MLEs a r e asymptotically normally distributed. Al- though the MLEs ha v e these imp ortan t features, they may b e stro ng ly biased for small o r ev en mo derate samples sizes when more complex mo de ls a re considered. Th us, a bias correction can play an imp ortant role to improv e the estimation of the mo del par a me ters. An imp ortan t a re a o f researc h in statistics is the study of the finite-sample b e- ha vior of MLEs. It is w ell k on wn that MLEs are oftentimes bia s ed, t hus displayin g systematic error. This is not a s erious pro blem for relativ ely large sample sizes, since the bias is t ypically of order O ( n − 1 ), whereas the asymptotic standard errors are of order O ( n − 1 / 2 ). Ho wev er, for small or ev en mo derate v alues of the sample size n , bias can constitute a problem. Th us, a v ailability of form ulae for its approximate computation is imp ortan t for a c curate estimation of man y mo dels that are used in a num b er of applications. Bias correction of MLEs is particularly imp ortant when the sample size, or the to tal information, is small (V asconcellos and Cribari–Neto, 2005). Bias adjustment has b een extensiv ely studied in the statistical literature. Bo x (1971) giv es a general expression for the n − 1 bias in multiv ariate nonlinear mo dels where cov aria nc e matrices are known. F or nonlinear regression mo dels, Co ok et al. (1986) relate bias to the p osition of the explanatory v ariables in the sample space. Cordeiro and McCullagh (1991) giv e general matrix form ulae for bias correction in generalized linear m o dels. Cordeiro and V asconcellos (1997) obtained general matrix form ulae for bia s correction in multiv ar iate nonlinear regression mo dels with normal errors, while V asconcellos and Cordeiro (199 7 ) obtained general formulae for bias in m ultiv ariate nonlinear heteroscedastic regression. Also, Cordeiro and V a sconcellos (1999) obtained second order biases of the maximum lik eliho od estimators in von Mises regression mo dels, while Cordeiro et al. (2000) obtained bias correction for symmetric nonlinear regression mo dels. V a s concellos and Cordeiro ( 2 000) obtained bias correction for m ultiv ariate nonlinear Studen t t r egr ession mo dels. More recen tly , V asconcellos and Cribari–Neto (2005) obtained bias correction in a new class of b eta regressions. Cordeiro and Dem ´ etrio (2008) deriv e form ulae fo r the second-order biases of the quasi-maxim um likelihoo d estimators, while Cordeiro and T o y ama (2008) deriv e g eneral for mulae for the second-order biases in generalized nonlinear mo dels with disp e rsion co v ariates. In this pap er w e study a multiv aria te normal mo del with general parameteriza- tion and deriv e the second-order biases of the maxim um lik eliho od estimates. Here, the general parameterization means a sort of unification of sev eral imp ortan t mo dels whic h can b e constructed using the m ultiv ariate normal mo del. F or instance, the m ultiv ariate nonlinear regressions studied by Cordeiro and V asconcellos (1 997) and their heteroscedastic v ersion ( V asconcellos and Cordeiro, 1997) are just particular cases of our pro p osal. In this pap er w e prop ose a mo del in whic h the mean µ and the v ariance Σ of the observ ed v ariables are indexed by the same ve ctor o f par ame ters, sa y θ . The existing w orks on bias correction assume that the mean and v ar ianc e do not share any parameters, how ev er, in errors-in-v ariables mo dels, for example, this assumption is not realistic. Indeed , that assumption mak es the computation of the bias formulae less complicated, but it restricts t h e applicabilit y of the approa c h to a sp ecial class of mo dels. In view of that, the main goal of this article is to extend 2 the bias correction to a wide class of multiv ar ia te mo dels whic h has not y et b een considered in the statistical litera t u re. The o u tline o f the pap er is as follo ws. Section 2 presen ts the main mo del and computes the first three deriv at ives of the log-lik eliho od function and their exp ec- tations. In Section 3, w e presen t matrix formulae fo r the second-order biases of the MLEs for the general mo del. In Section 4, w e presen t some useful examples of the prop osed formu lation. Monte Carlo sim ulation results a r e presen ted and discussed in Section 5. The num erical results show that the bias correction w e deriv e is ef- fectiv e in small samples; it deliv ers estimators that are nearly un biased and displa y sup e rior finite-sample b eha vior. Section 6 contains a n empirical ilustration. Finally , Section 7 concludes the pa p er. 2 Mo del sp e c ification W e consider the situation in whic h n indep enden t m ultiv ariate random v a r ia ble s Y 1 , . . . , Y n are observ ed a nd the n um ber of resp onses measured in eac h observ atio n is q . W e also admit that auxiliary co v ariates can b e observ ed, sa y X 1 , . . . , X n . The m ultiv ariate regression mo del can then b e represen ted as Y i = µ i ( θ ) + u i , i = 1 , 2 , . . . , n, (1) where Y i is a q × 1 v ector o f dep enden t v ariables, µ i ( θ ) ≡ µ i ( θ , X i ) is a mean func- tion (the shap e is assumed known) whic h is three times con tin uously differen tiable with resp ect to each elemen t o f θ and X i is a n m × 1 v ector of kno wn explanatory v ariables asso ciated with the i th observ ed resp onse Y i . Also, θ is a p × 1 vec tor of unkno wn para me ters of in terest. Additionally , as the foundation for estimation b y maxim um lik eliho od and h ypo t h esis testing, w e assume that t he indep enden t random errors u i ’s follow a m ultiv ariate norma l N q ( 0 , Σ i ( θ )) distribution, where Σ i ( θ ) is a q × q nonsingular co v ariance matrix and the en tries of Σ i ( θ ) ar e assumed three times contin uously differen tiable in eac h elemen t of θ . W e are assuming, in addition, that the functions µ i ( θ ) and Σ i ( θ ) are defined in a con v enien t wa y since θ mu st b e identifiable in mo del (1) . The class of mo dels presen ted ab o v e is very rich and includes man y imp ortan t regression mo dels. F or example, in an errors-in-v aria b les model w e observ e tw o v ariables, namely Y i and X i whose relationship is g iv en b y Y i = α + β x i + e i and X i = x i + u i , (2) where x i ∼ N ( µ x , σ 2 x ), e i ∼ N (0 , σ 2 ) and u i ∼ N (0 , σ 2 u ), with σ 2 u kno wn and, additionally , x i , e i and u i are m utually uncorrelated. Then, denoting Y i = ( Y i , X i ) ⊤ and θ = ( α, β , µ x , σ 2 x , σ 2 ) ⊤ w e ha v e that Y i ∼ N 2 ( µ ( θ ) , Σ ( θ )), where µ ( θ ) =  α + β µ x µ x  and Σ ( θ ) =  β 2 σ 2 x + σ 2 β σ 2 x β σ 2 x σ 2 x + σ 2 u  . This is a simple linear regression in whic h the cov aria t e is sub ject to measure- men t erro r s . This is a go o d example where the usual approach (assuming that Σ 3 and µ do not share any parameter) is not applicable. Measuremen t error mo de ls ha v e b een largely used in epidemiology (Kulathinal et al., 2002; de Castro et al., 2008), astrophy sics (Akritas and Bershady , 1996; Kelly, 2007; Kelly et al., 2008) and analytical c hemistry (Cheng a nd Riu, 2006) to av oid inconsisten t estimators (see F uller, 198 7, for f urt h er details). Other sp e cial cases of mo del (1 ) are: m ulti- v ariate heteroscedastic nonlinear errors-in-v ariables mo dels , m ultiv ariate no n linear heteroscedastic mo dels, univ ariate nonlinear mo dels, factor analysis, mixed mo dels and so on. As can b e seen, mo del (1) can encompass a wide class o f mo dels. T o simplify t h e no tation, define Y = v ec( Y 1 , Y 2 , . . . , Y n ), µ = ve c ( µ 1 ( θ ) , . . . , µ n ( θ )), Σ = diag { Σ 1 ( θ ) , . . . , Σ n ( θ ) } and u = Y − µ , where ve c( · ) is the ve c op erator, whic h transforms a matrix into a v ector by stac king the columns of the matrix one un- derneath the other. T hen, the log-likelihoo d asso c iated with (1), apart fro m an unimp ortan t constant, is ℓ ( θ ) ∝ − 1 2 log | Σ | − 1 2 tr { Σ − 1 uu ⊤ } . (3) W e mak e some assumptions (Co x and Hinkley, 1974, Ch. 9) on t he behav ior of ℓ ( θ ) as the sample size n approac hes infinity , suc h as the regularit y of the first three deriv ative s of ℓ ( θ ) with resp ect to θ and the uniqueness of the MLE of θ , b θ . W e now in tro duce the f o llo wing tot al log-lik eliho od deriv atives , in whic h the indices r , s and t range fro m 1 to p . Le t U r = ∂ ℓ ( θ ) /∂ θ r , U sr = ∂ 2 ℓ ( θ ) /∂ θ s ∂ θ r and U tsr = ∂ 3 ℓ ( θ ) /∂ θ t ∂ θ s ∂ θ r b e the first three deriv atives o f ℓ ( θ ). Th e standard notation for the momen ts of t ho s e log- lik eliho o d deriv ativ es is used (Lawley , 195 6), namely: κ sr = E ( U sr ), κ s,r = E ( U s U r ), κ tsr = E ( U tsr ) and so o n. F urthermore, w e define the deriv ativ e of κ sr with resp ect to θ t as κ ( t ) sr = ∂ κ sr /∂ θ t . Not all κ ’s are functionally indep enden t ; e.g., κ s,r = − κ sr , whic h is the t ypical elemen t of t he information matrix K θ , assumed to b e nonsingular. All κ ’s refer to a total ov er the sample and are, in general, of order n . Finally , let κ s,r denote the corresp onding elemen t o f K − 1 θ . Define the f o llo wing quan tities: a r = ∂ µ ∂ θ r , a sr = ∂ 2 µ ∂ θ s ∂ θ r , a tsr = ∂ 3 µ ∂ θ t ∂ θ s ∂ θ r , C r = ∂ Σ ∂ θ r , C sr = ∂ 2 Σ ∂ θ s ∂ θ r , C tsr = ∂ 3 Σ ∂ θ t ∂ θ s ∂ θ r , A r = ∂ Σ − 1 ∂ θ r = − Σ − 1 C r Σ − 1 and A sr = ∂ A r ∂ θ s , where r , s, t = 1 , 2 , . . . , p . W e assume that these deriv ativ es do exist. T o compute the deriv ative s of ℓ ( θ ) we mak e use of metho ds in matrix differen tial calculus, as describe in Magnus a n d Neudec k er (1988). Th us, the first deriv ativ e of (3) with resp e ct to the r th elemen t of θ is U r = 1 2 tr { A r ( Σ − uu ⊤ ) } + tr { Σ − 1 a r u ⊤ } . By using some simple matrix pro perties, the score function for θ can b e written in matrix form a s U θ ≡ U θ ( θ ) = e D ⊤ Σ − 1 u − 1 2 e V ⊤ e Σ − 1 v ec( Σ − uu ⊤ ) , 4 where e D = ( a 1 , . . . , a p ), e V = (ve c( C 1 ) , . . . , vec( C p )), e Σ = Σ ⊗ Σ and ⊗ is the Kronec k er pro duct. Let e F = e D e V ! , f H =  Σ 0 0 2 e Σ  − 1 and e u =  u − v ec( Σ − uu ⊤ )  . (4) Then, note that the score function can b e written as U θ = e F ⊤ f H e u . The second and third deriv ativ es are giv en, resp ectiv ely , b y U sr = 1 2 tr { ( A s Σ A r + A r Σ A s − Σ − 1 C sr Σ − 1 )( Σ − uu ⊤ ) } + 1 2 tr { A r ( C s + a s u ⊤ + u a ⊤ s ) } + tr { ( A s a r + Σ − 1 a sr ) u ⊤ − Σ − 1 a r a ⊤ s } and U tsr = U (1) tsr + U (2) tsr + U (3) tsr + U (4) tsr , where U (1) tsr = − 1 2 tr { ( Σ − 1 C tsr Σ − 1 + Σ − 1 C r s A t + A t C sr Σ − 1 )( Σ − uu ⊤ ) + Σ − 1 C sr Σ − 1 C t + Σ − 1 C sr Σ − 1 ( a t u ⊤ + u a ⊤ t ) } , U (2) tsr = 1 2 tr { ( A st Σ A r + A s ( C t A r + Σ A r t ))( Σ − uu ⊤ ) + ( A r t Σ A s + A r ( C t A s + Σ A st ))( Σ − uu ⊤ ) + ( A s Σ A r + A r Σ A s )( C t + a t u ⊤ + u a ⊤ t ) } , U (3) tsr = 1 2 tr { A r t ( C s + a s u ⊤ + u a ⊤ s ) + A r ( C ts + a ts u ⊤ − a s a ⊤ t + a t a ⊤ s + u a ⊤ ts ) } and U (4) tsr = tr { ( A st a r + A s a r t + A t a sr + Σ − 1 a sr t ) u ⊤ − ( A s a r + Σ − 1 a sr ) a ⊤ t } − tr { Σ − 1 a r a ⊤ ts + A t a r a ⊤ s + Σ − 1 a s a ⊤ tr } . Note that E ( u ) = 0 a nd E ( uu ⊤ ) = Σ . Knowing these prop erties, the exp ectation of U r , U sr and U tsr are easily obtained. The quan tities κ sr , κ tsr and κ ( r ) ts ( r , s, t = 1 , 2 , . . . , p ) are given, resp ec tiv ely , by κ sr = 1 2 tr { A r C s } − a ⊤ s Σ − 1 a r , (5) κ tsr = tr { ( A r Σ A s + A s Σ A r ) C t } + 1 2 tr { A s C tr + A r C ts + A t C sr } − ( a ⊤ t A s a r + a ⊤ s A t a r + a ⊤ s A r a t + a ⊤ t Σ − 1 a sr + a ⊤ ts Σ − 1 a r + a ⊤ s Σ − 1 a tr ) (6) 5 and κ ( r ) ts = 1 2 tr { ( A r Σ A s + A s Σ A r ) C t + A t C r s + A s C r t } − ( a ⊤ r t Σ − 1 a s + a ⊤ t A r a s + a ⊤ t Σ − 1 a r s ) . (7) Again, b y using some matrix prop erties on expression (5 ) , w e can written the exp ected Fisher inf o rmation as K θ ≡ K θ ( θ ) = e D ⊤ Σ − 1 e D + 1 2 e V ⊤ e Σ − 1 e V . Using e F and f H given in (4), w e can write K θ in the form K θ = e F ⊤ f H e F . (8) The MLE b θ satisfy the equation U θ = 0 . Af t er some matrix manipulat io ns , the Fisher scoring metho d can b e used to estimate θ b y iterativ ely solving the equation ( e F ( m ) ⊤ f H ( m ) e F ( m ) ) θ ( m +1) = e F ( m ) ⊤ f H ( m ) e u ∗ ( m ) , m = 0 , 1 , 2 , . . . , (9) where e u ∗ ( m ) = e F ( m ) θ ( m ) + e u ( m ) . Eac h lo op, through t he iterativ e sc heme (9), consists of an iterative re-w eigh ted least squares algorithm to optimize the log -lik eliho o d (3). Using equation (9 ) and any soft w are ( MAPLE , MATLAB , Ox , R , SAS ) with a w eigh ted lin- ear regression routine one can compute t he MLE b θ iteratively . It is also noteworth y that the MLE in even m uc h complex mo dels, such as multiv aria te heteroscedastic nonlinear errors-in- v aria bles mo dels, may b e atta in ed via iterativ e fo rm ula (9). It is w ell kno wn t ha t MLEs are consisten t, asymptotically efficien t and asymp- totically normal distributed. W e can write b θ a ∼ N p ( θ , K − 1 θ ), when n is lar ge, a ∼ denoting approxim ately distributed. Hence, h ypo t h eses testing can b e carried out using this asymptotic distribution. 3 Biases o f e stimates of θ In this section w e compute the biases of ML estimates o f θ fo r mo dels defined b y (1). Let B ( b θ a ) b e the n − 1 bias of b θ a , a = 1 , 2 , . . . , p . It follows fro m the general expression for the multiparameter n − 1 biases of MLEs giv en by Co x and Snell (1968) t h at B ( b θ a ) = X ′ t,s,r κ a,t κ s,r  1 2 κ tsr − κ ts,r  , where P ′ denotes t he summation o v er all com binations of the parameters θ 1 , . . . , θ p . F ollowing Cordeiro and Klein (1994), we write t he ab o v e equation in matrix nota tion to obtain n − 1 bias vector B ( b θ ) of b θ in the form B ( b θ ) = K − 1 θ W v ec ( K − 1 θ ) , 6 where W = ( W (1) , . . . , W ( p ) ) is a p × p 2 partitioned matr ix, each W ( r ) referring to the r th comp onen t of θ b eing a p × p matrix with typical ( t, s )th elemen t give n by w ( r ) ts = 1 2 κ tsr + κ ts,r = κ ( r ) ts − 1 2 κ tsr . Notice that fro m (6) and (7) w e ha v e that w ( r ) ts = 1 4 tr { A t C sr + A s C tr − A r C ts } − 1 2 ( a ⊤ t Σ − 1 a sr + a ⊤ s Σ − 1 a tr − a ⊤ r Σ − 1 a ts ) + 1 2 ( a ⊤ s A t a r + a ⊤ t A s a r − a ⊤ t A r a s ) . (10) Since K θ is a symmetric matrix and w e are interested in the m ultiplication r esult of W v ec( K − 1 θ ), many terms of (10 ) cancel. Indeed, note that the t th elemen t of W v ec( K − 1 θ ) is give n b y w (1) t 1 κ 1 , 1 + ( w (1) t 2 + w (2) t 1 ) κ 1 , 2 + · · · + ( w ( s ) tr + w ( r ) ts ) κ s,r + · · · + ( w ( p − 1) tp + w ( p ) t ( p − 1) ) κ p − 1 ,p + w ( p ) tp κ p,p and w ( s ) tr + w ( r ) ts = 1 2 tr( A t C sr ) − a ⊤ t Σ − 1 a sr + a ⊤ s A t a r . Therefore, we can replace the elemen t w ( r ) ts b y 1 4 tr( A t C sr ) − 1 2 a ⊤ t Σ − 1 a sr + 1 2 a ⊤ s A t a r and W ( r ) ma y b e written in an equiv alen t w a y as W ( r ) = e F ⊤ f H Φ r , r = 1 , . . . , p , where Φ r = − 1 2 ( G r + J r ) with G r =  a 1 r · · · a pr v ec( C 1 r ) · · · v ec( C pr )  and J r =  0 2( I nq ⊗ a r ) e D  , where I m denotes an m × m iden tit y matrix. That is, the matrix W can b e written as W = e F ⊤ f H ( Φ 1 , . . . , Φ p ). Then, w e arriv e at the fo llo wing theorem. Theorem . The n − 1 bias ve ctor B ( b θ ) of b θ is given by B ( b θ ) = ( e F ⊤ f H e F ) − 1 e F ⊤ f H e ξ , (11) wher e e ξ = ( Φ 1 , . . . , Φ p ) ve c (( e F ⊤ f H e F ) − 1 ) . In order to in terpret formulae (11) it is helpful to exploit the relationship b et wee n the n − 1 bias of b θ and a linear regr ession. The bia s v ector B ( b θ ) is simply the set co effic ien ts from the ordinary w eigh ted lest-squares regression of e ξ on the columns of e F , using w eigh ts in f H . As expression (11) mak es clear, for a n y particular mo del of the class of mo dels presen ted in Section 2, it is alwa ys p ossible to express the bias of b θ a s the solution of an o rdinary w eigh ted lest-squares regression. Equation (11) is easily handled alg e braically f o r any type of nonlinear mo del, since it only in v olv es simple op erations on matr ices and ve ctors. This equation, in conjunction with a computer algebra system suc h as MAPLE (Ab ell and Braselton, 1994 ), will yield B ( b θ ) algebraically with minimal effort. Also, w e can compute the bias B ( b θ ) nume rically via soft w are with nume rical linear alg e bra facilities suc h as Ox (D o ornik, 200 6 ) a nd R ( R D evelopme n t Core T eam, 2006). [Note that we hav e describ ed a pro ce dure to attain a corrected estimator in a general form ulation that cov ers a wide class of mo dels. In the next section w e shall presen t some sp ecial cases to shed light on the applicabilit y of our general f o rm ulation.] Therefore, we are able to compute the n − 1 biases of the MLEs for the general mo del (1) using form ula (11 ) . On the righ t-hand side of express ion (1 1), whic h is 7 of order n − 1 , consisten t estimates of the parameter θ can b e inserted to define the corrected MLE e θ = b θ − b B ( b θ ), where b B ( · ) denotes the v alue o f B ( · ) a t b θ . The bias-corrected estimate e θ is exp ected to ha v e b etter sampling prop erties than t he uncorrected estimator, b θ . In fact, w e presen t some sim ulations in Section 5 that indicate that e θ has smaller bia s than its corresp onding MLE, th us suggesting that the bias corrections ab o v e ha v e the effect of shrinking the mo dified estimates tow ard to the true parameter v alues. F ollowing Cordeiro (2 0 08), it is w orth emphasizing that there are o the r metho ds to bias-correcting MLEs. In regular para metric problems, Firth (1993) dev elop ed the so-called “prev en tive” metho d, whic h also a llo ws fo r the remov al of the second- order biases. His metho d consists of mo difying the original score function to remo v e the first- order term from the asymptotic bias of these estimates. In exp onen tial families with canonical parameterization his correction sc heme consists in p enalizing the lik eliho o d by the Jeffreys in v ariant prior. This is a prev en tive approac h to bias adjustmen t whic h has its merits, but the connections b et w een our results and his w ork are not pursued in this pap er since they will b e dev elop e d in future researc h. Additionally , we sh ould also stress that it is possible to av oid cum b e rsome and tedious a lgeb ra on cum ulant calculations b y using Efron’s bo otstrap (Efron and Tibshirani, 19 93). W e use t he analytical approach here since it leads to a closed-form solution and we do not need to run extensiv e n umerical resamples. Moreov er, the application of the analytical bias seems to generally b e the most feasible pro cedure to use and it contin ues to receiv e at te n tion in the literature. 4 Sp ecial mo del s It is useful to consider some examples to illustrate the applicabilit y o f the results in the previous section and clarif y the no t ation used. Other imp ortan t sp ecial mo dels could also b e easily handled since fo rm ula ( 1 1) only requires simple op erations on matrices and vec tors. First, consider a univ ar ia te no n linear mo del ( q = 1) in whic h Σ = σ 2 I n . Note that t his mo de l is a pa r tic ular case of model (1) with θ = ( β ⊤ , σ 2 ) ⊤ and µ = ( µ 1 ( β ) , . . . , µ n ( β )) ⊤ , where the comp onen ts of µ and Σ are unrelated and v ary indep e nden tly . Let p − 1 b e the dimension o f β . H ere, e D = ( a 1 , . . . , a p − 1 , 0 ) and e V = ( 0 , vec( C p )). Also, e F = diag { e D (1) , e V (2) } , where e D (1) = ( a 1 , . . . , a p − 1 ) and e V (2) = v ec( C p ). Then, from (8), the exp ected Fisher informat ion for θ can b e written as K θ = e F ⊤ f H e F = diag { K β , K σ 2 } , where K β = e D (1) ⊤ e D (1) /σ 2 is Fisher’s information for β and K σ 2 = n/ 2 σ 4 is the information relativ e to σ 2 . Since K θ is blo c k-diag onal, β and σ are globally o rthogonal (Co x a nd R eid, 1987). F rom (11) , note that ( e F ⊤ f H e F ) − 1 e F ⊤ f H = " ( e D (1) ⊤ e D (1) ) − 1 e D (1) ⊤ 0 0 1 n e V (2) ⊤ ( I n ⊗ I n ) # . Also, e ξ = e ξ 1 e ξ 2 ! = " − σ 2 2 ¨ G v ec { ( e D (1) ⊤ e D (1) ) − 1 } − P p − 1 k =1 ( I n ⊗ a k ) e D (1) K − 1 β k # , 8 where ¨ G = ( a β 1 , . . . , a β ( p − 1) ) with a β k = ( a 1 k , . . . , a ( p − 1) k ) and K − 1 β k is the k th column of K − 1 β . Then, B ( b θ ) =  B ( b β ) B ( b σ 2 )  = " ( e D (1) ⊤ e D (1) ) − 1 e D (1) ⊤ e ξ 1 1 n e V (2) ⊤ ( I n ⊗ I n ) e ξ 2 # . Note that B ( b β ) = ( e D (1) ⊤ e D (1) ) − 1 e D (1) ⊤ e ξ 1 agrees with the result due to Co ok et al. (1 9 86, Equation (3)). Additionally , w e obtain the following simple form originally first giv en by Beale (1960 ): B ( b σ 2 ) = − ( p − 1) σ 2 /n ; note that e V (2) ⊤ ( I n ⊗ I n ) p − 1 X k =1 ( I n ⊗ a k ) e D (1) K − 1 β k = p − 1 X k =1 v ec( C p ) ⊤ ( I n ⊗ a k ) e D (1) K − 1 β k = p − 1 X k =1 tr { a k K − 1 β k e D (1) ⊤ } = tr { e D (1) K − 1 β e D (1) ⊤ } = ( p − 1) σ 2 . As a second application, consider the m ultiv ariate nonlinear heteroscedastic re- gression mo del studied b y V asconcellos a nd Cordeiro (1997) . Note that this mo del is a particular case of mo del (1), with θ = ( β ⊤ , σ ⊤ ) ⊤ , µ = v ec( µ 1 ( β ) , . . . , µ n ( β )) and Σ = diag { Σ 1 ( σ ) , . . . , Σ n ( σ ) } . Therefore, the comp onen ts of µ and Σ are un- related and v ary indep ende n tly . Let p 1 and p 2 = p − p 1 b e t he dimensions of β and σ , respective ly . Here, e D = ( a 1 , . . . , a p 1 , 0 ) a n d e V = ( 0 , ve c( C p 1 +1 ) , . . . , v ec( C p )). Let e D (1) = ( a 1 , a 2 , . . . , a p 1 ) and e V (2) = (v ec( C p 1 +1 ) , v ec( C p 1 +2 ) . . . , v ec( C p )), then e F = diag { e D (1) , e V (2) } . F rom (8), the exp ected Fisher info rmation fo r θ can b e written as K θ = e F ⊤ f H e F = diag { K β , K σ } , where K β = e D (1) ⊤ Σ − 1 e D (1) is Fisher’s information for β and K σ = 1 2 e V (2) ⊤ e Σ − 1 e V (2) is the informatio n relativ e to σ . Since K θ is blo c k-diagonal, β and σ are globally orthogonal. F rom (11 ), it follows that ( e F ⊤ f H e F ) − 1 e F ⊤ f H = " ( e D (1) ⊤ Σ − 1 e D (1) ) − 1 e D (1) ⊤ Σ − 1 0 0 ( e V (2) ⊤ e Σ − 1 e V (2) ) − 1 e V (2) ⊤ e Σ − 1 # . Also, e ξ = e ξ 1 e ξ 2 ! = " − 1 2 ¨ G v ec { ( e D (1) ⊤ Σ − 1 e D (1) ) − 1 } −  ¨ W vec { ( e V (2) ⊤ e Σ − 1 e V (2) ) − 1 } + P p 1 k =1 ( I nq ⊗ a k ) e D (1) K − 1 β k  # , where ¨ G = ( a β 1 , . . . , a β p 1 ) with a β k = ( a 1 k , . . . , a p 1 k ) and ¨ W = ( v σ ( p 1 +1) , . . . , v σ p ) with v σ k = (v ec( C ( p 1 +1) k ) , . . . , vec( C pk )) and K − 1 β k is the k th column of K − 1 β . There- fore, B ( b θ ) =  B ( b β ) B ( b σ )  = " ( e D (1) ⊤ Σ − 1 e D (1) ) − 1 e D (1) ⊤ Σ − 1 e ξ 1 ( e V (2) ⊤ e Σ − 1 e V (2) ) − 1 e V (2) ⊤ e Σ − 1 e ξ 2 # . Note that B ( b β ) = ( e D (1) ⊤ Σ − 1 e D (1) ) − 1 e D (1) ⊤ Σ − 1 e ξ 1 agrees with the result due to V asconcellos and Cordeiro (1 997, Equation (3.2)). Additionally , note that B ( b σ ) = 9 ( e V (2) ⊤ e Σ − 1 e V (2) ) − 1 e V (2) ⊤ e Σ − 1 e ξ 2 also reduces to V a s concellos and Cordeiro’s (1997) Eq. (3.8), since e V (2) ⊤ e Σ − 1 p 1 X k =1 ( I nq ⊗ a k ) e D (1) K − 1 β k = e V (2) ⊤ e Σ − 1 v ec( ∆ ∗ ) , where ∆ ∗ is as defined b y V a s concellos and Cordeiro ( 1 997, p. 148). Next, unlik e the tw o mo dels discussed previously , we consider a mo del where the elemen ts of µ and Σ are related and do not v ary indep enden tly . Consider the nonlinear heteroscedastic errors-in-v ariables mo del Y i = α + β x i + exp ( γ z i ) + e i and X i = x i + u i , where x i ∼ N ( µ x , σ 2 x ) and u i ∼ N (0 , σ 2 u ) are the measuremen t errors with σ 2 u kno wn and e i ∼ N (0 , σ 2 exp { η z i } ). The co v ariate z i is known. In this example, the v ector of parameters is θ = ( α, β , γ , µ x , σ 2 x , σ 2 , η ) ⊤ and the mean and v ariance functions f or the i th observ ation ( Y i , X i ) are given b y µ i =  α + β µ x + exp ( γ z i ) µ x  and Σ i =  β 2 σ 2 x + σ 2 exp( η z i ) β σ 2 x β σ 2 x σ 2 x + σ 2 u  . Then, a 1 = 1 n ⊗  1 0  , a 2 = 1 n ⊗  µ x 0  , a 3 = v ec  z 1 exp( γ z 1 ) 0  · · ·  z n exp( γ z n ) 0  , a 4 = 1 n ⊗  β 1  and a 5 = a 6 = a 7 = 0 , where 1 n denotes an n × 1 v ector of ones. Also, a r s = 0 for a ll r, s except f o r a 24 = a 42 = 1 n ⊗  1 0  and a 33 = v ec  z 2 1 exp( γ z 1 ) 0  · · ·  z 2 n exp( γ z n ) 0  . Also, C r = 0 for a ll r except for C 2 = I n ⊗  2 β σ 2 x σ 2 x σ 2 x 0  , C 5 = I n ⊗  β 2 β β 1  , C 6 = n ⊕ i =1  exp( η z i ) 0 0 0  and C 7 = n ⊕ i =1  z i σ 2 exp( η z i ) 0 0 0  , where ⊕ is the direct sum of matrices. Additiona lly , C r s = 0 f or all r, s except for C 22 = I n ⊗  2 σ 2 x 0 0 0  , C 25 = C 52 = I n ⊗  2 β 1 1 0  , C 67 = C 76 = n ⊕ i =1  z i exp( η z i ) 0 0 0  and C 77 = n ⊕ i =1  z 2 i σ 2 exp( η z i ) 0 0 0  . 10 Th us, e D = ( a 1 , a 2 , a 3 , a 4 , 0 , 0 , 0 ), e V = ( 0 , v ec( C 2 ) , 0 , 0 , v ec( C 5 ) , v ec( C 6 ) , v ec( C 7 )) and the matrix formula (11) can b e used to compute the second-order bias for this mo del. Notice that, as v ec( C 2 ) is not equal to zero, t he deriv atio n of algebraic expression using matrix form ula (11) b ecomes v ery tedious, since the structure of K θ is not blo c k-dia g onal unlike the t w o previous examples. How ev er, using MAPLE , for example, the deriv ation can b e easily done. Also, the n − 1 bias v ector B ( b θ ) can b e obtained n umerically via softw are with n umerical linear algebra facilities with minimal effort suc h as R and Ox . 5 Sim ulation s tudy W e recall that, f or lar ge samples the biases of the MLEs are neglible. Ho w ev er, for small and mo derate sample sizes the second-order biases ma y b e large and can b e used to impro v e the estimation. W e shall use Monte Carlo sim ulation to ev aluate the finite sample p erformance of the original MLEs and their corrected v ersions. All sim ulations w ere p erformed using R (R Deve lopmen t Core T eam, 2006). The sample sizes considered w ere n = 15 , 25 , 35 , 5 0 and 100, the num b er of Mon te Carlo replications w as 5 ,000. W e consider the simple errors-in-v ariables mo del as describ ed in F uller (1987): Y i = α + β x i + e i and X i = x i + u i , where x i ∼ N ( µ x , σ 2 x ) and u i ∼ N (0 , σ 2 u ) are the measuremen t errors with σ 2 u kno wn and e i ∼ N (0 , σ 2 ), with i = 1 , 2 , . . . , n . Here, θ = ( α, β , µ x , σ 2 x , σ 2 ) ⊤ and µ = 1 n ⊗  α + β µ x µ x  and Σ = I n ⊗  β 2 σ 2 x + σ 2 β σ 2 x β σ 2 x σ 2 x + σ 2 u  . F rom the previous expressions, w e hav e immediately that a 1 = 1 n ⊗  1 0  , a 2 = 1 n ⊗  µ x 0  , a 3 = 1 n ⊗  β 1  , a 4 = a 5 = 0 and a r s = 0 for a ll r , s except for a 23 = a 32 = 1 n ⊗  1 0  . Also, C r = 0 for a ll r except for C 2 = I n ⊗  2 β σ 2 x σ 2 x σ 2 x 0  , C 4 = I n ⊗  β 2 β β 1  and C 5 = I n ⊗  1 0 0 0  . Additionally , C r s = 0 for a ll r, s except f o r C 22 = I n ⊗  2 σ 2 x 0 0 0  and C 24 = C 42 = I n ⊗  2 β 1 1 0  . Th us, e D = ( a 1 , a 2 , a 3 , 0 , 0 ) and e V = ( 0 , vec( C 2 ) , 0 , v ec( C 4 ) , v ec( C 5 )). Therefore, all the quantities necessary to calculate B ( b θ ) using expression (11) ar e g iv en. 11 In order to analyze the p oin t estimation results, w e computed, f o r each sample size and for each estimator: the relativ e bias (the relativ e bia s of an estimator b θ is defined a s { E ( b θ ) − θ } /θ , its estimate b eing obt a ine d b y estimating E ( b θ ) b y Monte Carlo) and t he ro ot mean square error , i.e., √ MSE, where MSE is the mean squared error estimated fro m the 5,000 Mon te Carlo replications. Without loss of generality , the true v alues of the regression par a me ters were set a t α = 67, β = 0 . 42, µ x = 70, σ 2 x = 247 and σ 2 = 43. The para m eter setting w ere choosen in order to represen t the dataset (yields of corn on Marshall soil in Io w a) presen ted in F uller (198 7, p. 18). The kno wn measuremen t error v ariance is σ 2 u = 57 (whic h w as attained through a previous exp erimen t). T able 1 give s the relativ e biases and √ MSE of b oth uncorrected and corrected estimates. The figures in this table confirm that the bias-corrected estimates ar e generally closer to the true pa r ame ter v alues than the unadjusted estimates. W e observ e that, in absolute v alue, t he estimated relativ e bia ses o f the bias-corrected estimator were smaller than t ha t of the original MLE fo r all sample sizes considered, th us sho wing the effectiv eness o f the bias correction sc hemes used in the definition of these estimators. F or instance, when n = 15, the estimated relativ e bias of the estimators of α , β , µ x , σ 2 x and σ 2 a v erage − 0 . 051 8 whereas the biases of the fiv e corresp onding bias- adjusted estimators a v erage − 0 . 0056; that is, the av erage bias (in v a lue absolute) of the MLEs is almost ten times la r g e r than that of the bias-corrected estimators. This suggests that the second-order bias of MLEs should no t b e ignor ed in samples of small to mo derate sizes since they can b e nonnegligible. W e can readily see tha t the MLEs of σ 2 x and σ 2 are on av erag e far from the true parameter v alue, t hus displa ying large bias, for the differen t sample sizes considered, ev en when n = 100. This stress es the imp ortance of using a bias correction. F or instance, when n = 50, the relativ e biases of b σ 2 x and b σ 2 (MLEs) w ere − 0 . 0226 and − 0 . 0563, resp ec tiv ely , while the relativ e biases of e σ 2 x and e σ 2 (BCEs) w ere 0 . 001 6 (sixteen times lesser) and − 0 . 0011 (fifty times lesser), resp ectiv ely . Observ e that the MLEs a re alwa ys underestimating the true v alues of σ 2 x and σ 2 , since their biases are alw a ys negative s. Note also that r o ot mean square erro r decrease with n , as exp e cted. Additionally , w e note that all estimators hav e similar ro ot mean squared errors. 6 An empirical i l ustration Next, as an empirical ilustration, consider a small data set give n by F uller (1 987, p. 18) . The data s et is presen ted in T able 2. The data are yields of corn and determinations of av aila b le soil nitrogen collected at 11 sites on Marshall soil in Io w a. F ollo wing F uller (198 7, p. 18 ), we a s sume that the estimates of soil nitrog en con tain measuremen t erros arise from t w o sources. F ir st, only a small sample o f soil is selected from eac h plot and, as a result, there is the sampling error asso c iated with the use of sample to r epresen t the whole. Second, there is a measuremen t erro r asso ciated with the c hemical analysis used to determined the lev el of nitrogen in the soil sample. The v ar ianc e arising from these tw o sources is σ 2 u = 57. According to 12 T able 1: Relativ e biases a nd √ MSE of uncorrected and cor r ected estimates for an errors-in-v aria bles mo del. MLE BCE n θ Rel. bias √ MSE Rel. bias √ MSE 15 α − 0 . 0240 12.46 0 . 0232 11.29 β 0 . 0547 0.17 − 0 . 052 6 0.16 µ x 0 . 0014 4.48 0 . 0014 4.48 σ 2 x − 0 . 0796 108.49 − 0 . 0029 113.8 1 σ 2 − 0 . 1807 19.52 0 . 0031 20.38 25 α − 0 . 0198 9.0 5 0 . 0009 8.14 β 0 . 0440 0.13 − 0 . 002 9 0.11 µ x 0 . 0004 3.43 0 . 0004 3.43 σ 2 x − 0 . 0553 85.73 − 0 . 0082 88.05 σ 2 − 0 . 1198 15.48 − 0 . 0104 15.73 35 α − 0 . 0117 7.0 5 0 . 0010 6.68 β 0 . 0267 0.10 − 0 . 002 3 0.09 µ x − 0 . 0001 2.9 6 − 0 . 0 0 01 2.96 σ 2 x − 0 . 0424 71.36 − 0 . 0084 72.64 σ 2 − 0 . 0799 12.83 − 0 . 0014 13.04 50 α − 0 . 0080 5.6 9 0 . 0002 5.50 β 0 . 0190 0.08 0 . 0005 0.08 µ x − 0 . 0007 2.4 5 − 0 . 0 0 07 2.45 σ 2 x − 0 . 0226 60.76 0 . 0016 61.71 σ 2 − 0 . 0563 10.75 − 0 . 0011 10.89 100 α − 0 . 0025 3.83 0 . 0013 3.78 β 0 . 0057 0.05 − 0 . 002 9 0.05 µ x 0 . 0002 1.72 0 . 0002 1.72 σ 2 x − 0 . 0131 42.24 − 0 . 0009 42.54 σ 2 − 0 . 0298 7.6 3 − 0 . 0 0 21 7.67 BCE: bias-corrected estimator. 13 F uller (198 7 , p. 18 ), mo del (2) is a v alid represen ta tion to these data. T able 2: Yields of corn on Marshall soil in Iow a. Soil Soil Yield Nitrogen Yield Nitrogen Site ( Y ) ( X ) Site ( Y ) ( X ) 1 8 6 70 7 99 50 2 115 97 8 96 70 3 9 0 53 9 99 94 4 8 6 64 10 104 69 5 110 95 11 96 51 6 9 1 64 The MLEs, t h e large-sample estimates of the corresp onding standard error s, the biases and the bias-corrected estimates are giv en T able 3. Thes e estimates w ere obtained using R . The figures in this table show that the biases of the estimates of α and β are muc h less than standard errors of the corresp onding estimates. In cases of margina l statistical significance, biases of this maginitude could hav e a small effect on the conclusions. Ho w ev er, note t ha t the MLEs of σ 2 x and σ 2 are strongly biased, as evidenced by our sim ulations studies, i.e., they underestimate the mo del v ariances. Therefore, the bias-corrected estimates may b e used instead of the MLEs to mak e p oin t inferences. T able 3: MLEs and bias-corrected estimativ es. P arameter MLEs S.E. Bias BCEs α 66.8606 11.7272 − 2 . 5334 69.3939 β 0.4331 0.1633 0 . 0359 0.3973 µ x 70.6364 5.0194 0 . 0000 70.6364 σ 2 x 220.1405 118.173 1 − 25 . 1946 245.3351 σ 2 38.4058 20.9357 − 10 . 3344 48.740 2 BCE: bias-corrected estimate. Figure 1 presen ts the scatterplot of the data together with the fitted lines ob- tained using the MLEs and BCEs. Notice that the line pro duced by the bias cor- rection sc heme the inclination sligh tly atten uated and in tercept increased relativ e to the non- c orrected one. 7 Conclus ions This pap er prop osed a bias corr e ction for a m ultiv ariate no r ma l mo del with a quite general pa rame terization. The main result cen ters on mo dels where t he mean and the v ariance share the same ve ctor of pa r a me ters. Man y mo dels are particular cases of the prop osed mo del such as (non)linear regr essions, errors-in-v ariables mo de ls, 14 50 60 70 80 90 85 90 95 100 105 110 115 Soil Nitrogen (X) Yield (Y) MLE BCE Figure 1: Scatterplot f o r the data set together with the fitted lines. mixed mo dels, factor analysis a nd so forth. W e ha v e sho wn that it is alw a ys p ossible to express the second order bias v ector of t h e maxim um lik eliho o d estimates as an ordinary w eigh ted least-squares regression. Moreov er, we deriv ed a bia s -adjustmen t sc heme that nearly eliminates the bias of the maxim um lik eliho o d estimator in small and mo derate samples. Our sim ulation results suggest that the bias correction w e ha v e deriv ed is very effectiv e, ev en when the sample size is small. Indeed, the bias correction mec hanism prop osed in this pap er yields mo dified maxim um lik eliho o d estimators that are nearly un biased. W e also presen t e d an empirical ilustration that illustrates the prop osed bias-adjustmen t sch eme. Ac kno w ledgmen ts W e gratefully ackno wledge g ran ts from F APE SP . W e wish to thank F rancisco Cribari- Neto for v alua b le commen ts on this manus cript. References Ab ell, M.L., Braselton, J.P . (1994). The M a ple V Handb o ok . AP Professional, New Y ork. Akritas, M.G., Bershady , M.A. (1996). Linear regression for astronomical data w it h mea- surement errors and int rinsic scatter. The Astr ophysic al Journal . 470 :706– 714. Beale, E . M.L. (1960). Confid e nce regions in non-linear estimation (with discussion). Jour- nal of the R oyal Statistic al So ciet B . 22 :41–88. Bo x, M.J. (1971 ). Bias in nonlin ear estimation (with discussion). Journal of the R oyal Statistic al So ciet B . 33 :171–20 1. 15 Cheng, C.L., Riu J. (2006 ). On estimating lin ear relationships when b oth v ariables are sub j e ct to heteroscedastic measurement err ors . T e chnometrics . 48 :511–519 . Co ok, R.D., Ts a i, C., W ei, B.(1986). Bias in nonlinear regression. Biometrika . 73 :615 –623. Cordeiro, G.M. (2008). Corrected maximum lik elihoo d estimators in linear heteroske dastic regression mo dels. Br azilian R eview of Ec onometrics . 28 :53–67 . Cordeiro, G.M., Dem ´ etrio, C.G.B. (2008). Corrected estimators in extended quasi- lik elihoo d m odels. Communic ations in Statistics, The ory and Me t ho ds . 37 :873– 880. Cordeiro, G.M., F err a ri, S .L.P ., Urib e-Opazo, M.A., V asconcellos, K.L.P . (2000). C or- rected maximum-lik eliho od estimation in a class of symmetric nonlinear regression mo d- els. Statistics and Pr ob ability L etters . 46 :317– 328. Cordeiro, G.M., Klein, R. (1994 ). Bias correction in ARMA mo dels. Statistics and Pr ob- ability L e tters . 19 :169–176 . Cordeiro, G.M., McCullagh, P . (1991). Bias correction in generalized lin e ar m o dels. Jour- nal of the R oyal Statistic al So ciety B . 53 :629–643 . Cordeiro, G.M., V asconcellos, K.L.P . (1997 ). Bias correction for a class of m ultiv ariate nonlinear regression mo dels. Statistics and P r ob ability L e tters . 35 :155–1 64. Cordeiro, G.M., V asconcellos, K .L.P . (1999) . Second-order b ia ses of the maxim um lik eli- ho od estimates in von Mises regression mo dels. Austr alian and N ew Ze ala nd Journal of Statistics . 41 :901–910 . Cordeiro, G.M., T o y ama, M.C. (200 8). Bias corr ection in generalized nonlinear mod- els with disp ersion co v ariates. Communic ations in Statistics, The ory and Metho ds . 37 :221 9–2225 . Co x, D.R., Hinkley , D.V. (1974). The or etic al Statistics . London: C hapman and Hall. Co x, D.R., Reid, N. (1987). Parameter orth ogonalit y and approximate cond i tional infer- ence (with discussion). J o urnal of the R oya l Statistic al So cie ty B . 40 :1–39 . Co x, D.R., S n el l, E. (1968 ). A general defin it ion of residuals (with discussion). J o urnal of the R oyal Statistic al So cie t y . 30 :248–2 75. de Castro, M., Galea, M., Bolfarine, H. (2008 ). Hyp othesis testing in an err o rs-in-v ariables mo del with h et eroscedastic measurement errors. Statistics in M e dicine . 27 :521 7–5234 . Do o rnik, J.A. (2006) . An Obj e c t -Oriente d M a trix L anguage – O x 4 . Tim b erlak e Consu l- tan ts Pr e ss, London. 5th ed. Efron, B., Tibsh ir a ni, R.J. (1993 ). A n intr o duction to the b o otstr ap . Chapman and Hall, New Y ork. Firth, D. (1993). Bias r e duction of maximum lik elihoo d estimates. Biometrika . 80 :27–38 . F uller, W. (1987). M e asur ement Err or Mo dels . Wiley: Chichester. 16 Kelly , B.C. (2007). Some asp ects of measuremen t error in linear regression of astronomical data. The Astr oph ysic al Journal . 665 :1489–15 06. Kelly , B.C., Bec h told, J., T rump, J.R., V ertergaard, M., Siemigino wsk a, A. (2008). O b- serv ational constrain ts on the dep endence of ratio-quiet quasar X-ra y emission on blac k hole mass and accretion rate. Astr oph ysic al Journal Supplement Series . 176 (2): 355– 373. Kulathinal, S.B., Kuu la smaa, K., Gasbarra, D. (2002). E s t imation of an errors-in-v ariables regression mo del w hen the v ariances of the measurement err or v ary b et w een the obser- v ations. Statistics in Me dici ne . 21 :1089 –1101 . La wley , D.N. (1956) . A general metho d for approxi mating to the distr i bution of lik eliho od ratio criteria. Biometrika . 43 :295–30 3. Magn us, J.R., Neudec k er, H. (1988). Matrix differ ential c alculus . Wiley: New Y ork. R Deve lopmen t Core T eam (2006). R: A Language and E nvironmen t for Statistical Com- puting. Vienna, Au stria. IS BN 3-900051-0 7-0. http ://www.R -project .org. V asconcellos, K.L.P ., Cord eiro, G.M. (1997). Ap p ro ximate b ias for m ultiv ariate nonlinear heteroscedastic r egressions. Br azilian Journal of P r ob ability and Statistics . 11 :141– 159. V asconcellos, K.L.P ., Cordeiro, G.M. (2000). Bias corrected estimates in multiv ariate Stu- den t t r egression mo dels. Communic ations in Statistics, The ory and Metho ds . 29 :797– 822. V asconcellos, K.L.P ., C r ibari–Neto, F. (2005). Improv ed maximum like liho o d estimation in a new class of b eta regression mo dels. Br azilian Journal of P r ob ability and Statistics . 19 :13– 31. 17