Parametric simultaneous robust inferences for regression coefficient under generalized linear models

(1)

On: 24 December 2014, At: 18:45 Publisher: Taylor & Francis

Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Click for updates

Journal of Statistical Computation and

Simulation

Publication details, including instructions for authors and subscription information:

http://www.tandfonline.com/loi/gscs20

Parametric simultaneous robust

inferences for regression coefficient

under generalized linear models

Li-Chu Chiena & Tsung-Shan Tsoubcd a

Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan, Republic of China

b

Institute of Statistics, National Central University, Jhongli, Taiwan, Republic of China

c

Institute of Systems Biology and Bioinformatics, Center for Biotechnology and Biomedical Engineering, National Central University, Jhongli, Taiwan, Republic of China

d

Cathay Medical Research Institute, Cathay General Hospital, Taipei, Taiwan, Republic of China

Published online: 15 Oct 2012.

To cite this article: Li-Chu Chien & Tsung-Shan Tsou (2014) Parametric simultaneous robust

inferences for regression coefficient under generalized linear models, Journal of Statistical Computation and Simulation, 84:4, 850-867, DOI: 10.1080/00949655.2012.731409

To link to this article: http://dx.doi.org/10.1080/00949655.2012.731409

PLEASE SCROLL DOWN FOR ARTICLE

Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

(2)

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions

(3)

Vol. 84, No. 4, 850–867, http://dx.doi.org/10.1080/00949655.2012.731409

Parametric simultaneous robust inferences for regression

coefficient under generalized linear models

Li-Chu Chiena* and Tsung-Shan Tsoub,c,d

a_{Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan, Republic of China;}b_{Institute of}

Statistics, National Central University, Jhongli, Taiwan, Republic of China;cInstitute of Systems Biology and Bioinformatics, Center for Biotechnology and Biomedical Engineering, National Central University, Jhongli, Taiwan, Republic of China;d_{Cathay Medical Research Institute, Cathay General Hospital, Taipei,}

Taiwan, Republic of China

(Received 4 October 2011; final version received 14 September 2012)

In this article, the parametric robust regression approaches are proposed for making inferences about regression parameters in the setting of generalized linear models (GLMs). The proposed methods are able to test hypotheses on the regression coefficients in the misspecified GLMs. More specifically, it is demonstrated that with large samples, the normal and gamma regression models can be properly adjusted to become asymptotically valid for inferences about regression parameters under model misspecification. These adjusted regression models can provide the correct type I and II error probabilities and the correct coverage probability for continuous data, as long as the true underlying distributions have finite second moments.

Keywords: generalized linear models; robust normal regression; robust gamma regression

1. Introduction

Generalized linear models (GLMs) were introduced by Nelder and Wedderburn [1] as a unifying family of models for non-standard cross-sectional regression analysis with non-normal responses. The statistical analysis of such models is based on the asymptotic properties of the maximum likelihood estimator (MLE). Fahrmeir and Kaufmann [2] presented mild general conditions, which, respectively, assure weak or strong consistency or asymptotic normality of the MLE. More on this study can be found in [3]. More generally, Fahrmeir [4] dealt with the asymptotic behaviour of the quasi-MLE in misspecified GLMs.

Cantoni and Ronchetti [5] proposed a natural class of robust estimation techniques for GLMs. Their method is more reliable than the classical estimation procedures in providing the accurate statistical inference when the data include outlying points. Adimari and Ventura [6] also studied robust inference for GLMs. They derived a robust quasi-profile log-likelihood function that was obtained from an estimating function that defines the class of Mallows-type robust estimators considered by Cantoni and Ronchetti [5]. Li and Hsiao [7] suggested a method for consistently estimating GLMs with measurement errors without making any prior distributional assumption on

*Corresponding author. Email: lcchien@nctu.edu.tw

(4)

Journal of Statistical Computation and Simulation 851

the measurement error or the latent variables. However, the robustness of their proposed method requires the knowledge of the probability distribution of latent variables. Sinha [8] developed a robust method for analysing GLMs with non-ignorable missing covariates. Recently, Bianco et

al. [9] introduced a resistant procedure to test hypotheses on the regression parameter in GLMs

with missing responses.

On the other hand, Heagerty and Kurland [10] evaluated the impact of model violations on the estimate of a regression coefficient in generalized linear mixed models (GLMMs). Jiang and Zhang [11] proposed robust methods to estimate parameters of interest in settings of GLMMs, in which only the conditional means of the responses given the random effects are specified. Yau and Kuk [12] proposed robust estimation procedures for GLMMs based on the notion of maximum quasi-likelihood and residual maximum quasi-quasi-likelihood. Sinha [13] developed a robust method for identifying and downweighting the outliers when estimating the parameters in the GLMMs. Sinha [14] further described a robust quasi-likelihood method for fitting the GLMMs to longitudinal data.

In addition, robust restricted maximum likelihood (robust REML) in mixed linear models are introduced by Richardson and Welsh [15] who made classical REML robust by bounding the influence of outlying observations on the estimate. Yun and Lee [16] discussed the robust estimation in mixed linear models with non-monotone missingness. Jacqmin-Gadda et al. [17] investigated the robustness of the MLE of fixed effects from a linear mixed model when the error terms are either correlated or non-Gaussian or of non-constant variance.

Royall and Tsou [18] advocated the robust likelihood function concept. They developed a technique that adjusts a working likelihood function, making it robust. The resulting adjusted robust likelihood function remains valid evidential representation of the parameter, even when the working model is incorrect. Motivated by the above results, Tsou [19] proposed a parametric robust way for comparing two population means and two population variances in misspecified models. Tsou and Cheng [20] applied the robust likelihood techniques to analyse contaminated data in regression settings. Tsou [21] further extended the robust likelihood concept to analyse count data. In this article, the robust likelihood techniques are used to make inferences about regression parameters in the GLM setting.

This article is organized as follows. Section 2 contains a brief review of the idea of robust likelihood functions introduced by Royall and Tsou [18]. The robust normal regression (RNR) and robust gamma regression (RGR) are briefly introduced in Section 3. Section 4 presents a simulation study which shows the advantage of the RNR and RGR models with respect to (w.r.t.) the ordinary normal and gamma regression models. Section 5 concludes with a brief discussion. Some technical background material from the previous sections is deferred into the appendix.

2. Robust likelihood functions

Suppose that Y1, Y2, . . . , YNis a sequence of independent random variables. On the basis of a

pri-ori knowledge or convenience, we postulate a working model for the probability distributions of Yi’s,{fi= fi(•; ψ) = f (•; ηi(ψ)), i= 1, 2, . . . , N, ψ ∈ }, where ψ is a fixed-dimensional vector of unknown parameters. For example, under normal regression settings, ηi(ψ)= (xtiγ, σ2), ψ=

(γt, σ2₎t_{and f}

i= fi(yi; ψ)= exp{−(yi− xtiγ)2/2σ2}/ √

2π σ . Here xirepresents the p character-istics that are specific to yi, and γ represents the p regression coefficients that describe how xi affects the expected value of Yi. Note that this model is a collection of probability distributions, each of which is identified by a unique value of ψ.

Now partition ψ as ψt= (θt, ϕt₎_{, where θ is the p-vector of parameters of interest and ϕ is the} remaining fixed-dimensional nuisance parameters. Let θ0and ϕ0 denote the limiting values of

(5)

the MLEs, ˆθ and ˆϕ, based on the working model f = (f1, f2, . . . , fN), when the Yi’s are actually generated from the family{hi= h(•; τi(θ, λ)), i= 1, 2, . . . , N}, where λ is the nuisance parameter vector under h= (h1, h2, . . . , hN). Now suppose that the parameters of inference under the working model f , namely θ, remain the parameters of interest under h, so that θ0has the same interpretation of the true values of the parameters of interest. This result is what Royall and Tsou [18] referred to as the first condition of robustness (FCR). This condition is crucial for the working model to be adjustable for valid inferences. Note that White [22] showed that, more often than not, the FCR is not satisfied once f = h.

Write lθ and lϕ for the first derivatives of the log-likelihood function l(θ, ϕ) w.r.t. θ and ϕ,

respectively, whose derivatives w.r.t. ϕ are correspondingly denoted by lθϕand lϕϕ. Now, let Ihθϕ and Ihϕϕ be the limiting values of−lθϕ/N and−lϕϕ/N, respectively, under h and the limiting

values of−lθθ/Nand−lϕθ/N, under h, are denoted by Ihθθand Ihϕθ, respectively. Note that these limiting values are all evaluated at θ0and ϕ0.

Now define the following two p× p matrices:

A= Ihθθ− IhθϕI−1_hϕϕIhϕθ (1) and

B= Vhθθ− IhθϕI−1hϕϕVhϕθ− VhθϕI−1hϕϕIhϕθ+ IhθϕI−1hϕϕVhϕϕI−1hϕϕIhϕθ. (2) Here Vhθθ= limN→∞Eh[lθ(θ0, ϕ0)ltθ(θ0, ϕ0)/N],Vhθϕ= limN→∞Eh[lθ(θ0, ϕ0)lϕt(θ0, ϕ0)/N] and Vhϕϕ= limN→∞Eh[lϕ(θ0, ϕ0)lϕt(θ0, ϕ0)/N], where Eh stands for the expectation evaluated under h.

Let ˆθ be the MLE of θ and ˆA and ˆB be the empirical versions of A and B. A direct application of Taylor’s expansion shows that the adjusted Wald statistic N(ˆθ− θ0)tˆA ˆB

−1_{ˆA(ˆθ − θ}

0)has an asymp-totic χ2

p distribution for general hi, i= 1, 2, . . . , N, that have finite second moments. Here χp2is denoted as a chi-squared distribution with p degrees of freedom.Another asymptotically equivalent counterpart, the adjusted score statistic N−1{lt

θ(θ0,ˆϕ(θ0))} ˆB −1

(θ0,ˆϕ(θ0)){lθ(θ0,ˆϕ(θ0))}, where

ˆϕ(θ0)and ˆB(θ0, ϕ(θ0))are the constrained MLEs of ϕ and B given θ0, respectively, has the same limiting χ2

p distribution even if the working model assumptions fail.

3. Robust regression models

Consider a set of observations y1, y2, . . . , yN corresponding to N independent not identically distributed random variables Y1, Y2, . . . , YN. Under GLMs, the mean response, μi, depends on the

p covariates (xi,0, xi,1, . . . , xi,p−1)= xti, by μi= g(ηi), where ηi= xtiγ= γ0xi,0+ γ1xi,1+ · · · +

γp−1xi,p−1is a linear predicator with the p regression coefficients (γ0, γ1, . . . , γp−1)= γt, and g(•) is a monotonic and differentiable response function.

3.1. Robust normal regression

Under a normal working model, the log-likelihood function for the ith observation yiis

li= − 1 2log σ 2₋1 2log 2π− (yi− μi)2 2σ2 .

(6)

The log-likelihood equation for γj−1is 1 σ2 N i=1 μ_i(yi− μi)xi,j−1= 0, j = 1, 2, . . . , p, (3)

where μiis the first derivative of μiw.r.t. ηi. The solutions of Equation (3) are the maximum quasi-likelihood (MQL) estimators [23–25] or M-estimators [26], when the observations y1, y2, . . . , yN are not necessarily from normal distributions. McCullagh [27] showed that, under mild regularity conditions, the consistency of the MQL estimates under model misspecification depends only on the correct specification of the regression. In other words, the normal working model provides the consistent estimates of regression parameters under incorrectly specified models. Thus, the FCR is fulfilled, so that the normal working model can be properly adjusted to become asymptotically legitimate for regression parameters of interest under model misspecification.

Without loss of generality, let γp−w, γp−w+1, . . . , γp−2, γp−1 be the w parameters of interest and let (γp−1, γp−2, . . . , γp−w+1, γp−w)be denoted by (β1, β2, . . . , βw−1, βw)= βt for notational convenience. Let μi,0 and μi,0 be, respectively, the true values of μi and μi. Let Varh(Yi), i= 1, 2, . . . , N, be the true variances of Yi, i= 1, 2, . . . , N. Let Z = (z0, z1, . . . , zp₋₁)be the N× p design matrix, so that Zt= (x1, x2, . . . , xN).

After lengthy derivations, it shows that the (u, v), u, v= 1, 2, . . . , w, elements of the w × w adjusting matrices An and Bn of AnB−1n Anthat make the normal regression model robust can be written in the forms (for details, see the appendix):

An(uv)= lim N_→∞ 1 Nσ₀2 N i=1 (μ_i,0)2 ⎛ ⎝xi,p−u− p_−w j=1 |Δnj(u)| |Δn| xi,j−1 ⎞ ⎠ ⎛ ⎝xi,p−v− p_−w j=1 |Δnj(v)| |Δn| xi,j−1 ⎞ ⎠ and Bn(uv)= lim N_→∞ 1 Nσ₀4 N i₌₁ Varh(Yi)(μ_i,0)2 ⎛ ⎝xi,p−u− p−w j₌₁ |Δnj(u)| |Δn| xi,j−1 ⎞ ⎠ × ⎛ ⎝xi,p−v− p−w j₌₁ |Δnj(v)| |Δn| xi,j−1 ⎞ ⎠ .

Here |Δn| represents the determinant of the matrix Δn, where Δn= WtV−1n W with W=

(z0, . . . , zj−2, zj−1, zj, . . . , zp−w−1) and Vn= diag((1/μ_1,0)2, (1/μ_2,0)2, . . . , (1/μ_N,0)2) being a diagonal matrix of order N. On the other hand, Δnj(u)= WtV−1n Wj(u)and Δnj(v)= WtV−1n Wj(v)with Wj(u)= (z0, . . . , zj−2, zp−u, zj, . . . , zp−w−1)and Wj(v)= (z0, . . . , zj−2, zp−v, zj, . . . , zp−w−1)derived by the jth column zj−1 of W replaced by zp−u and zp−v, respectively. Here σ₀2is the limit of the MLE of σ2_, _ˆσ2_{, that has the same interpretation of the limit of}N

i=1Varh(Yi)/N. Note that the interpretation of σ2

0 depends on h and is, therefore, unknown.

In the special case with all the regression coefficients of interest, let (γp₋₁, γp₋₂, . . . , γ1, γ0)=

(β1, β2, . . . , βp₋₁, βp)= βt. Then, the adjusting matrices Anand Bncan be simplified as follows:

An(uv)= lim N_→∞ 1 Nσ2 0 N i=1

(μ_i,0)2(xi,p−u)(xi,p−v) and Bn(uv)= lim N→∞ 1 Nσ₀4 N i₌₁

Varh(Yi)(μ_i,0)2(xi,p−u)(xi,p−v).

(7)

In applications, consistent estimates ˆAn and ˆBn of An and Bn can be obtained by Varh(Yi) replaced by (yi− ˆμi)2with ˆμ being the MLE of μ and other unknown quantities replaced by their respective empirical versions.

Let β0 be the true value of β and consider the null hypothesis H0: β= β0. Let ˆβ be the MLE of β based on the normal working model and let ˆϕ(β0) and ˆBn(β0,ˆϕ(β0)) be the restricted MLEs of ϕ and Bn given β0. Here the vector of the nuisance parameters, ϕ, con-tains the scale parameter σ2 _{and some regression coefficients that are not to be tested. Under}

H0, the adjusted Wald statistic N( ˆβ− β0) tˆA

nˆB −1

n ˆAn( ˆβ− β0) and the adjusted score statis-tic N−1{lt

β(β0, ϕ(β0))} ˆB −1

n (β0,ˆϕ(β0)){lβ(β0,ˆϕ(β0))} are asymptotically equivalent and have an asymptotic χw2 distribution as long as the second moments of the true underlying distributions exist. Note that ˆAnˆB

−1

n ˆAn is free of σ2. Thus, with large samples, the effect of σ2 is actually removed. Hence, σ2_{can be treated known, a priori, as any arbitrary positive value.}

3.2. Robust gamma regression

Under a gamma working model, the log-likelihood function for the ith observation yiis

li= r log r − r log μi+ (r − 1) log yi− rμ−1i yi− log (r). The score functions

r N i=1 μ_i μi yi− μi μi xi,j−1, j= 1, 2, . . . , p,

have zero expectation as long as μi, i= 1, 2, . . . , N, are correctly specified. Hence, the regression parameters of interest can be consistently estimated by the gamma working model, whatever h is. Thus, the FCR is satisfied.

Calculations parallel to An and Bnshow that the (u, v), u, v= 1, 2, . . . , w, components of the adjusting matrices Agand Bgunder the gamma working model are of the forms (for details, see the appendix): Ag(uv)= lim N_→∞ r0 N N i=1 _μ i,0 μi,0 2⎛ ⎝xi,p−u− p_−w j=1 |Δgj(u)| |Δg| xi,j−1 ⎞ ⎠ ⎛ ⎝xi,p−v− p_−w j=1 |Δgj(v)| |Δg| xi,j−1 ⎞ ⎠ and Bg(uv)= lim N→∞ r₀2 N N i₌₁ Varh(Yi) μ2_i,0 _μ i,0 μi,0 2⎛ ⎝x_i,p_−u₋p−w j₌₁ |Δgj(u)| |Δg| xi,j−1 ⎞ ⎠ × ⎛ ⎝xi,p−v− p−w j=1 |Δgj(v)| |Δg| xi,j−1 ⎞ ⎠ ,

where Δg=WtV−1g W, Δgj(u)=WtV−1g Wj(u)and Δgj(v)=WtV−1g Wj(v)with Vg= diag((μ1,0/μ1,0)2,

(μ2,0/μ2,0)2, . . . , (μN,0/μN,0)2). Here r0is the limit of the MLE of r, whose interpretation depends on h and is, therefore, unknown.

(8)

In the special case with all the regression parameters of interest, Agand Bgreduce to

Ag(uv)= lim N_→∞ r0 N N i₌₁ _μ i,0 μ_i,0 2

(xi,p−u)(xi,p−v)

and Bg(uv)= lim N_→∞ r2 0 N N i₌₁ Varh(Yi) μ2_i,0 _μ i,0 μ_i,0 2

(xi,p−u)(xi,p−v).

In application, consistent estimates ˆAgand ˆBgof Agand Bgcan be derived by replacing the unknown components in Ag and Bg by their respective empirical analogues, just as we dealt with ˆAn and ˆBn. Note that ˆAgˆB

−1

g ˆAg is free of r, so that with large samples, the effect of r is completely eliminated. Therefore, r can be treated known, a priori, in the beginning as any pos-itive value. The resulting adjusted Wald statistic N( ˆβ− β₀)tˆA

gˆB −1

g ˆAg( ˆβ− β₀)and the resulting adjusted score statistic N−1{ltβ(β0,ˆϕ(β0))} ˆB

−1

g (β0,ˆϕ(β0)){lβ(β0,ˆϕ(β0))}, under H0: β= β0, are asymptotically distributed as χ2

wfor general h with the finite second moments. Note that here all MLEs are derived under the gamma working model.

4. Simulation studies

To investigate the performance of the RNR and RGR models in the finite sample situation, sim-ulation studies are conducted using N = 450, 900 and 1350 replicated samples, respectively, generated from the three regression models

Model 1: μi= exp(ηi), Model 2: μi= (2.5ηi+ 2/3)3, Model 3: μi= η2i

with the linear predicator ηigiven by

ηi= xi,0+ γ1xi,1+ γ2xi,2 for i= 1, 2, . . . , N,

where the values of xi,0, i= 1, 2, . . . , N, are set by 1 and the values of xi,j, i= 1, 2, . . . , N, j = 1, 2, are independently generated from a uniform distribution between 0 and 1. Here regression coefficients γ1and γ2are considered as the parameters of interest. For simplicity, let γt= (γ1, γ2) and γt₀= (1.0, 1.0). We test the null hypothesis H0 : γ= γ0 and the two alternative hypotheses

HA: γt= (0.4, 1.0) and γt = (0.7, 1.3), respectively.

Simulated data sets are generated from three sources including the Weibull, inverse Gaussian and chi-squared distributions, respectively. A Weibull distribution with the shape parameter λ and the scale parameter k, W (k, λ), has a simple relationship between the second central moment and the first moment, that is, Var(Y )= aμ2_{, where a > 0 is a function of the shape parameter λ. For} example, when λ= 1 and k = μ, Var(Y) = μ2_{. Similarly, an inverse Gaussian distribution with} the mean μ and the shape parameter λ, IG(μ, λ), has a variance proportional to the cubic of its mean value, that is, Var(Y )= μ3_/λ_{. On the other hand, a non-central chi-squared distribution} with ν degrees of freedom and a non-centrality parameter μ− ν > 0, χ2

v(μ− v), has a mean value of μ and a variance of 2(2μ− ν), so that χ2

v(μ− v) has a variance roughly proportional to its mean.

(9)

To demonstrate the robustness characters of the adjusted Wald and score statistics under the normal and gamma working models, in our simulations, the observations, yi, i= 1, 2, . . . , N, are sampled in the following way. First, the first 0.3N observations, yi, i= 1, 2, . . . , 0.3N, are independently generated from the Weibull distributions, W (μi, 1), with the shape parameter of 1 and the scale parameter of μi, i= 1, 2, . . . , 0.3N, respectively. Then, the next 0.3N obser-vations, yi, i= 0.3N + 1, 0.3N + 2, . . . , 0.6N, are independently generated from the inverse Gaussian distributions, IG(μi, 100), with the shape parameter of 100 and the mean value of

μi, i= 0.3N + 1, 0.3N + 2, . . . , 0.6N, respectively. Finally, the rest of the 0.3N observations,

yi, i= 0.6N + 1, 0.6N + 2, . . . , N, are independently generated from the non-central chi-squared distributions, χ₁2(μi− 1), with one degree of freedom and the non-centrality parameter of

μi− 1, i = 0.6N + 1, 0.6N + 2, . . . , N, respectively.

Three additional test statistics are also included for contrast. They are the maximum likelihood ratio test statistic, the Wald test statistic and the score test statistic, respectively. The maximum likelihood ratio test statistic for testing the null hypothesis H0: γ = γ0is defined by

QL = 2{l(ˆγ, ˆϕ) − l(γ0,ˆϕ(γ0))}.

Then, its two asymptotically equivalent test statistics, the Wald test statistic and the score test statistic, are defined by

QW= N(ˆγ − γ0)tˆA(ˆγ − γ0) and

QS= N−1{lγt(γ0,ˆϕ(γ0))} ˆA −1

(γ₀,ˆϕ(γ₀)){lγ(γ₀,ˆϕ(γ₀))},

respectively. Here all notation definitions are given as in the previous sections. Each of the test statistics QL, QWand QS, under the null hypothesis H0 : γ= γ0, has an asymptotic chi-squared distribution with degrees of freedom equal to the dimension of γ. Thus, in our simulations, each of the maximum likelihood ratio test statistic, the Wald test statistic and the score test statistic rejects H0, when each of the test statistics QL, QW and QSexceeds the critical value of χ2,0.952 , where χ2

2,0.95 represents the 95th quantile of the chi-squared distribution χ22. More discussions

about the test statistics QL, QWand QScan be found in [28, Section 9.3].

The simulation performance are carried out for 3000 simulation runs with the xi,1’s and xi,2’s being regenerated after every 50 simulation runs. The empirical type I error probabilities based on the adjusted Wald statistic, the adjusted score statistic, the maximum likelihood ratio test statistic, the Wald test statistic and the score test statistic are labelled as AWα, ASα, Lα, Wα and Sα, respectively. On the other hand, AWcp, AScp, Lcp, Wcp and Scp symbolize the coverage probabilities of the nominal 95% confidence interval constructed using the adjusted Wald statistic, the adjusted score statistic, the maximum likelihood ratio test statistic, the Wald test statistic and the score test statistic, respectively. The empirical type I error probability is computed as the proportion of rejections of the null hypothesis H0: γ= γ0at the nominal 5% significance level, when the data are actually generated from H0. On the other hand, when the data are sampled from the alternative hypothesis HA, the empirical type I error probability exhibits the power of the test. Results from the adjusted Wald statistic, the adjusted score statistic, the maximum likelihood ratio test statistic, the Wald test statistic and the score test statistic based on the normal and gamma working models are tabulated in the tables below. The average of the 3000 ˆγ values and their sample covariance matrix are termed as mean(ˆγ) and S2_{, respectively. In order to contrast} the differences between the covariance matrix estimates based on the adjusted and unadjusted test statistics, the average of the unadjusted covariance matrix estimate ofˆγ, namely ˆA−1/N, denoted

(10)

Journal of Statistical Computation and Simulation 857

Table 1. Model 1: μi= exp(ηi), i= 1, 2, . . . , N. Working

model mean(ˆγ) S2 _Var

A(ˆγ) Var(ˆγ) AWα AWcp ASα AScp Lα Lcp Wα Wcp Sα Scp

N_{= 450} H0: γ= 1.0 1.0 Normal 1.0030 1.0024 0.0223 0.0054 0.0242 0.0219 0.0040 0.0216 0.0161 0.0001 0.0161 0.0643 0.9357 0.0587 0.9413 0.1190 0.8810 0.1223 0.8777 0.1157 0.8843 Gamma 0.9992 0.9998 0.0143 0.0005 0.0150 0.0145 0.0004 0.0143 0.0165 0.0001 0.0165 0.0597 0.9403 0.0407 0.9530 0.0307 0.9693 0.0290 0.9710 0.0320 0.9680 HA: γ= 0.4 1.0 Normal 0.3980 1.0038 0.0198 0.0018 0.0203 0.0193 0.0014 0.0194 0.0152 0.0001 0.0177 0.9753 0.9377 0.9673 0.9413 0.9860 0.9133 0.9883 0.9137 0.9827 0.9187 Gamma 0.3980 1.0016 0.0154 0.0003 0.0158 0.0154 0.0001 0.0154 0.0184 0.0001 0.0184 0.9950 0.9350 0.9947 0.9473 0.9923 0.9690 0.9927 0.9690 0.9920 0.9687 HA: γ= 0.7 1.3 Normal 0.7027 1.3032 0.0227 0.0049 0.0248 0.0219 0.0036 0.0227 0.0147 0.0001 0.0181 0.8280 0.9393 0.8150 0.9443 0.8573 0.8803 0.8587 0.8730 0.8577 0.8850 Gamma 0.7005 1.2987 0.0145 0.0005 0.0151 0.0145 0.0003 0.0145 0.0166 0.0001 0.0166 0.8970 0.9453 0.8860 0.9530 0.8603 0.9647 0.8637 0.9653 0.8573 0.9667 N= 900 H0: γ= 1.0 1.0 Normal 1.0033 1.0003 0.0120 0.0026 0.0112 0.0113 0.0024 0.0109 0.0081 0.0000 0.0080 0.0603 0.9397 0.0593 0.9407 0.1187 0.8813 0.1203 0.8797 0.1183 0.8817 Gamma 1.0012 0.9995 0.0074 0.0001 0.0072 0.0073 0.0002 0.0072 0.0083 0.0000 0.0082 0.0593 0.9407 0.0517 0.9483 0.0357 0.9643 0.0357 0.9643 0.0357 0.9643 HA: γ= 0.4 1.0 Normal 0.3992 1.0009 0.0102 0.0010 0.0098 0.0099 0.0008 0.0096 0.0076 0.0000 0.0088 0.9997 0.9403 0.9993 0.9437 0.9997 0.9120 1.0000 0.9100 0.9997 0.9137 Gamma 0.3984 0.9997 0.0082 0.0002 0.0076 0.0078 0.0001 0.0077 0.0092 0.0000 0.0091 1.0000 0.9407 1.0000 0.9477 1.0000 0.9647 1.0000 0.9633 1.0000 0.9657 HA: γ= 0.7 1.3 Normal 0.7019 1.3011 0.0118 0.0023 0.0115 0.0113 0.0022 0.0113 0.0073 0.0000 0.0090 0.9853 0.9400 0.9847 0.9437 0.9900 0.8793 0.9893 0.8763 0.9890 0.8790 Gamma 0.7008 1.3005 0.0073 0.0003 0.0070 0.0073 0.0001 0.0072 0.0083 0.0000 0.0082 0.9953 0.9453 0.9953 0.9510 0.9937 0.9647 0.9937 0.9643 0.9930 0.9963

(11)

L.-C. Chien and T .-S. Tsou Table 1. Continued Working

N= 1350 H0: γ= 1.0 1.0 Normal 1.0022 1.0015 0.0076 0.0014 0.0073 0.0073 0.0014 0.0073 0.0054 0.0000 0.0054 0.0560 0.9440 0.0543 0.9457 0.1123 0.8877 0.1130 0.8870 0.1100 0.8900 Gamma 1.0013 1.0014 0.0049 0.0000 0.0047 0.0048 0.0000 0.0048 0.0055 0.0000 0.0055 0.0527 0.9473 0.0450 0.9550 0.0320 0.9680 0.0327 0.9673 0.0323 0.9677 HA: γ= 0.4 1.0 Normal 0.4003 1.0033 0.0068 0.0004 0.0068 0.0064 0.0005 0.0065 0.0050 0.0000 0.0059 1.0000 0.9367 1.0000 0.9393 1.0000 0.9093 1.0000 0.9083 1.0000 0.9093 Gamma 0.4003 1.0023 0.0054 −0.0001 0.0054 0.0052 0.0000 0.0052 0.0061 0.0000 0.0061 1.0000 0.9410 1.0000 0.9460 1.0000 0.9680 1.0000 0.9670 1.0000 0.9693 HA: γ= 0.7 1.3 Normal 0.7002 1.3016 0.0075 0.0013 0.0078 0.0073 0.0012 0.0076 0.0049 0.0000 0.0060 0.9983 0.9373 0.9983 0.9430 0.9983 0.8793 0.9983 0.8790 0.9983 0.8813 Gamma 0.6993 1.3005 0.0049 0.0002 0.0048 0.0048 0.0001 0.0049 0.0055 0.0000 0.0055 0.9993 0.9427 0.9993 0.9463 0.9993 0.9663 0.9993 0.9623 0.9993 0.9637

(12)

Journal of Statistical Computation and Simulation 859 Table 2. Model 2: μi= (2.5ηi+ 2/3)3, i= 1, 2, . . . , N. Working

model mean(ˆγ) S2 _Var

N= 450 H0: γ= 1.0 1.0 Normal 1.0067 1.0107 0.0681 0.0342 0.0674 0.0566 0.0261 0.0531 0.0302 0.0017 0.0299 0.1123 0.8877 0.1090 0.8910 0.1723 0.8277 0.1763 0.8273 0.1650 0.8350 Gamma 0.9903 0.9946 0.0144 −0.0001 0.0141 0.0133 0.0002 0.0132 0.0099 −0.0003 0.0099 0.0877 0.9123 0.0640 0.9360 0.1190 0.8810 0.1200 0.8800 0.1167 0.8833 HA: γ= 0.4 1.0 Normal 0.3974 1.0059 0.0221 0.0061 0.0247 0.0209 0.0058 0.0225 0.0132 0.0005 0.0160 0.9753 0.9027 0.9690 0.9103 0.9767 0.8677 0.9783 0.8650 0.9737 0.8707 Gamma 0.3953 0.9963 0.0081 −0.0001 0.0090 0.0076 0.0001 0.0085 0.0066 −0.0001 0.0068 1.0000 0.9257 1.0000 0.9457 1.0000 0.9093 1.0000 0.9080 1.0000 0.9083 HA: γ= 0.7 1.3 Normal 0.7022 1.3157 0.0637 0.0317 0.0785 0.0553 0.0246 0.0605 0.0283 0.0017 0.0350 0.6840 0.8790 0.6817 0.8853 0.7183 0.8220 0.7190 0.8197 0.7127 0.8287 Gamma 0.6915 1.2931 0.0131 −0.0001 0.0149 0.0123 0.0002 0.0140 0.0096 −0.0002 0.0100 0.9270 0.9137 0.9107 0.9337 0.9530 0.8760 0.9557 0.8750 0.9497 0.8760 N= 900 H0: γ= 1.0 1.0 Normal 1.0094 1.0093 0.0289 0.0139 0.0271 0.0269 0.0131 0.0268 0.0152 0.0009 0.0152 0.0910 0.9090 0.0897 0.9103 0.1533 0.8467 0.1560 0.8440 0.1507 0.8493 Gamma 0.9997 0.9994 0.0068 0.0000 0.0065 0.0068 0.0001 0.0067 0.0049 −0.0002 0.0049 0.0657 0.9343 0.0510 0.9490 0.1047 0.8953 0.1067 0.8933 0.1037 0.8963 HA: γ= 0.4 1.0 Normal 0.4030 1.0058 0.0114 0.0032 0.0111 0.0107 0.0032 0.0112 0.0067 0.0002 0.0081 0.9983 0.9230 0.9967 0.9247 0.9980 0.8730 0.9987 0.8723 0.9980 0.8760 Gamma 0.3998 0.9998 0.0038 0.0000 0.0041 0.0038 0.0032 0.0043 0.0033 −0.0001 0.0034 1.0000 0.9443 1.0000 0.9500 1.0000 0.9217 1.0000 0.9223 1.0000 0.9223 HA: γ= 0.7 1.3 Normal 0.7081 1.3132 0.0289 0.0129 0.0302 0.0269 0.0125 0.0301 0.0146 0.0009 0.0177 0.9177 0.9047 0.9190 0.9060 0.9187 0.8457 0.9150 0.8447 0.9183 0.8500 Gamma 0.6998 1.2995 0.0063 0.0000 0.0071 0.0063 0.0001 0.0072 0.0048 −0.0002 0.0049 0.9983 0.9313 0.9973 0.9497 0.9987 0.8933 0.9990 0.8917 0.9987 0.8910

(13)

N_{= 1350} H0: γ= 1.0 1.0 Normal 1.0010 1.0048 0.0170 0.0083 0.0173 0.0169 0.0081 0.0174 0.0099 0.0005 0.0100 0.0903 0.9097 0.0897 0.9103 0.1560 0.8440 0.1563 0.8437 0.1513 0.8487 Gamma 0.9982 1.0005 0.0045 0.0002 0.0044 0.0045 0.0001 0.0045 0.0033 −0.0001 0.0033 0.0610 0.9390 0.0500 0.9500 0.1183 0.8817 0.1167 0.8833 0.1157 0.8843 HA: γ= 0.4 1.0 Normal 0.3990 1.0034 0.0070 0.0020 0.0073 0.0070 0.0020 0.0074 0.0044 0.0001 0.0053 1.0000 0.9313 1.0000 0.9340 1.0000 0.8737 1.0000 0.8727 1.0000 0.8743 Gamma 0.3994 1.0006 0.0025 0.0001 0.0028 0.0026 0.0001 0.0029 0.0022 0.0000 0.0023 1.0000 0.9460 1.0000 0.9530 1.0000 0.9107 1.0000 0.9113 1.0000 0.9103 HA: γ= 0.7 1.3 Normal 0.7001 1.3077 0.0172 0.0078 0.0193 0.0172 0.0077 0.0193 0.0095 0.0005 0.0116 0.9850 0.9123 0.9870 0.9160 0.9880 0.8430 0.9863 0.8400 0.9877 0.8450 Gamma 0.6986 1.3008 0.0042 0.0002 0.0048 0.0042 0.0001 0.0049 0.0032 −0.0001 0.0033 0.9997 0.9417 0.9997 0.9497 1.0000 0.8907 1.0000 0.8917 1.0000 0.8897

(14)

Journal of Statistical Computation and Simulation 861 Table 3. Model 3: μi= η2i, i= 1, 2, . . . , N. Working

model mean(ˆγ) S2 _Var

N= 450 H0: γ= 1.0 1.0 Normal 1.0023 0.9991 0.0242 0.0040 0.0232 0.0237 0.0038 0.0236 0.0203 0.0009 0.0203 0.0573 0.9427 0.0527 0.9473 0.0793 0.9207 0.0810 0.9190 0.0777 0.9223 Gamma 1.0020 0.9989 0.0162 −0.0001 0.0163 0.0163 −0.0003 0.0164 0.0191 −0.0008 0.0192 0.0487 0.9513 0.0377 0.9623 0.0270 0.9730 0.0263 0.9737 0.0293 0.9707 HA: γ= 0.4 1.0 Normal 0.3966 0.9971 0.0180 0.0015 0.0183 0.0181 0.0014 0.0178 0.0154 0.0004 0.0166 0.9833 0.9460 0.9813 0.9490 0.9877 0.9310 0.9880 0.9307 0.9873 0.9343 Gamma 0.3977 0.9968 0.0129 −0.0004 0.0137 0.0133 −0.0002 0.0134 0.0155 −0.0003 0.0162 0.9980 0.9370 0.9960 0.9550 0.9967 0.9697 0.9967 0.9690 0.9967 0.9680 HA: γ= 0.7 1.3 Normal 0.7001 1.3000 0.0244 0.0041 0.0256 0.0243 0.0034 0.0242 0.0199 0.0008 0.0215 0.7823 0.9383 0.7750 0.9413 0.7927 0.9187 0.7927 0.9163 0.7887 0.9217 Gamma 0.7001 1.2998 0.0158 −0.0004 0.0170 0.0161 −0.0003 0.0164 0.0187 −0.0007 0.0195 0.8360 0.9413 0.8287 0.9593 0.7793 0.9757 0.7850 0.9740 0.7723 0.9730 N= 900 H0: γ= 1.0 1.0 Normal 1.0004 1.0010 0.0128 0.0019 0.0117 0.0121 0.0020 0.0117 0.0102 0.0004 0.0101 0.0567 0.9433 0.0553 0.9447 0.0743 0.9257 0.0757 0.9243 0.0723 0.9277 Gamma 1.0003 1.0023 0.0087 −0.0003 0.0082 0.0083 −0.0002 0.0082 0.0096 −0.0005 0.0095 0.0627 0.9373 0.0503 0.9497 0.0367 0.9633 0.0380 0.9620 0.0350 0.9650 HA: γ= 0.4 1.0 Normal 0.4011 1.0015 0.0095 0.0008 0.0089 0.0092 0.0008 0.0088 0.0077 0.0002 0.0082 1.0000 0.9477 1.0000 0.9480 1.0000 0.9313 1.0000 0.9317 1.0000 0.9337 Gamma 0.4029 1.0022 0.0069 −0.0002 0.0069 0.0067 −0.0001 0.0067 0.0077 −0.0002 0.0080 1.0000 0.9383 1.0000 0.9513 1.0000 0.9717 1.0000 0.9703 1.0000 0.9697 HA: γ= 0.7 1.3 Normal 0.7019 1.3034 0.0130 0.0017 0.0115 0.0125 0.0019 0.0120 0.0100 0.0004 0.0107 0.9743 0.9450 0.9747 0.9487 0.9783 0.9207 0.9780 0.9227 0.9780 0.9247 Gamma 0.7015 1.3033 0.0081 0.0000 0.0082 0.0082 −0.0002 0.0082 0.0093 −0.0005 0.0096 0.9897 0.9433 0.9910 0.9503 0.9863 0.9673 0.9870 0.9647 0.9853 0.9683

(15)

N= 1350 H0: γ= 1.0 1.0 Normal 1.0025 1.0020 0.0082 0.0013 0.0080 0.0080 0.0013 0.0080 0.0068 0.0002 0.0068 0.0627 0.9373 0.0603 0.9397 0.0860 0.9140 0.0860 0.9140 0.0857 0.9143 Gamma 1.0007 0.9999 0.0058 −0.0001 0.0054 0.0055 −0.0001 0.0055 0.0064 −0.0003 0.0064 0.0580 0.9420 0.0507 0.9493 0.0307 0.9693 0.0340 0.9660 0.0303 0.9697 HA: γ= 0.4 1.0 Normal 0.4021 1.0004 0.0061 0.0005 0.0063 0.0061 0.0005 0.0060 0.0051 0.0001 0.0055 1.0000 0.9423 1.0000 0.9427 1.0000 0.9243 1.0000 0.9243 1.0000 0.9253 Gamma 0.4028 1.0003 0.0046 0.0000 0.0047 0.0045 −0.0001 0.0045 0.0051 −0.0001 0.0054 1.0000 0.9380 1.0000 0.9447 1.0000 0.9660 1.0000 0.9650 1.0000 0.9643 HA: γ= 0.7 1.3 Normal 0.7009 1.3015 0.0085 0.0009 0.0081 0.0082 0.0012 0.0082 0.0066 0.0002 0.0072 0.9960 0.9453 0.9957 0.9473 0.9950 0.9157 0.9950 0.9150 0.9947 0.9187 Gamma 0.7007 1.3014 0.0056 −0.0002 0.0055 0.0055 −0.0001 0.0055 0.0062 −0.0003 0.0064 0.9990 0.9480 0.9993 0.9513 0.9983 0.9707 0.9987 0.9687 0.9987 0.9710

(16)

by Var(ˆγ) and the average of the adjusted covariance matrix estimate of ˆγ, namely ˆA−1ˆB ˆA−1/N, denoted by VarA(ˆγ) are also included. Note that because with large samples the adjusted Wald and score statistics under the normal and gamma working models are free of σ2and r, the non-regression parameters σ2and r in the RNR and RGR models are treated known, a priori, as the same arbitrarily chosen value of 1, respectively.

From Tables 1–3, it is evident that the adjusting matrices successfully correct the normal and gamma working models and make them robust. As can be seen from Tables 1–3, the averages of the adjusted covariance matrix estimates, VarA(ˆγ), are nearly equivalent to the sample covariance matrix of ˆγ, S2_{, whereas the averages of the unadjusted covariance matrix estimates, Var(}_{ˆγ), are} different from S2_.

It is also observed that when the simulated data sets are generated under the null hypothesis

H0, the adjusted Wald and score statistics are more effective than the test statistics QL, QWand QS in providing the correct type I error probabilities. As can be seen from Tables 1–3, when the data are generated from H0, the values of AWα and ASα are more close to the nominal significance level 0.05, in contrast with the values of Lα, Wα and Sα.

On the other hand, it is noted that when the simulated data sets are generated under the alternative hypothesis HA, the adjusted Wald and score statistics not only rightly reject the null hypothesis H0 but also provide the right confidence region. As can be seen from Tables 1–3, when the data are generated from HA, the values of AWα and ASα gradually approach the value of 1 and the values of AWcp and AScp inchmeal approximate to the nominal confidence level 0.95, as the sample size N increases. On the contrary, the test statistics QL, QWand QS, under HA, only succeed in rejecting H0, but they do not provide the exactly correct confidence region. For example, in the case of Model 2 with the sample size N = 1350 and the alternative hypothesis HA: γt= (0.7, 1.3), respectively, the values of Lcp, Wcp and Scp under the gamma working model, 0.8907, 0.8917, and 0.8897, are far from the nominal confidence level 0.95, in comparison with the values of AWcp and AScp under the gamma working model, 0.9417 and 0.9497.

Obviously, from the results of Tables 1–3, it is enough to verify that the adjusted Wald and score statistics based on the normal and gamma working models furnish a foundation for valid inferences for the regression parameters of interest, even though the true underlying distributions are not from these two working models. Despite the fact that the RNR and RGR models remain the robustness property in misspecified models, some finite sample differences are revealed in the numerical performances.

The results in Tables 1–3 apparently display that the adjusted covariance matrix estimates, Var(ˆγ), under the gamma working model are smaller than that under the normal working model. This explicitly means that in terms of the hypothesis testing, the RGR model is more powerful than the RNR model for non-negative continuous data.

5. Concluding remarks

We propose the parametric robust regression methods in the GLM setting. The proposed meth-ods can provide the valid inferences about the regression parameters of interest under model misspecification.

The adjusting matrices for the normal and gamma working models are submitted here. They successfully adjust these two working models into robust models, whatever the true underlying distributions are, as long as their second moments exist. The two adjusted models, namely the RNR and RGR models, warrant the asymptotically legitimate inferences under model misspecification. Simulation studies illustrate that the RGR model is more efficient for more general non-negative continuous random variables.

(17)

One of the many attractive features of our proposed methods is that with large samples, the effect of σ2 _{of the normal model and the effect of r of the gamma model are entirely purged} by their respective adjustments. Hence, although the non-regression parameters σ2 _{and r are} artificially given positive values, the asymptotic validity of the RNR and RGR models are always obtained.

Finally, we noted that the above discussion was centred on the case of all continuous random variables. For robust inferences for count data, we refer the interested readers to [21].

Acknowledgements

The authors are grateful for helpful comments from the associate editor and the reviewer, which improved the quality of this work.

References

[1] J.A. Nelder and R.W.M. Wedderburn, Generalized linear models, J. R. Statist. Soc. Ser. A 135 (1972), pp. 370–384. [2] L. Fahrmeir and H. Kaufmann, Consistency and asymptotic normality of the maximum likelihood estimator in

generalized linear models, Ann. Stat. 13 (1985), pp. 342–368.

[3] L. Fahrmeir, Asymptotic testing theory for generalized linear models, Statistics 18 (1987), pp. 65–76.

[4] L. Fahrmeir, Maximum likelihood estimation in misspecified generalized linear models, Statistics 21 (1990), pp. 487–502.

[5] E. Cantoni and E. Ronchetti, Robust inference for generalized linear models, J. Am. Statist. Assoc. 96 (2001), pp. 1022–1030.

[6] C. Adimari and L. Ventura, Robust inference for generalized linear models with application to logistic regression, Statist. Probab. Lett. 55 (2001), pp. 413–419.

[7] T. Li and C. Hsiao, Robust estimation of generalized linear models with measurement errors, J. Econom. 118 (2004), pp. 51–65.

[8] S.K. Sinha, Robust methods for generalized linear models with nonignorable missing covariates, Canad. J. Stat. 36 (2008), pp. 277–299.

[9] A.M. Bianco, G. Boente, and I.M. Rodrigues, Robust tests in generalized linear models with missing responses, Comput. Stat. Data Anal. (in press), doi: 10.1016.csda.2012.05.008.

[10] P.J. Heagerty and B.F. Kurland, Misspecified maximum likelihood estimates and generalised linear mixed models, Biometrika 88 (2001), pp. 973–985.

[11] J. Jiang and W. Zhang, Robust estimation in generalised linear mixed models, Biometrika 88 (2001), pp. 753–765.

[12] K.K.W. Yau, and A.Y.C. Kuk, Robust estimation in generalized linear mixed models, J. R. Statist. Soc. Ser. B 64 (2002), pp. 101–117.

[13] S.K. Sinha, Robust analysis of generalized linear mixed models, J. Am. Statist. Assoc. 99 (2004), pp. 451–460. [14] S.K. Sinha, Robust inference in generalized linear models for longitudinal data, Canad. J. Stat. 34 (2006),

pp. 261–278.

[15] A.M. Richardson and A.H. Welsh, Robust restricted maximum likelihood in mixed linear models, Biometrics 51 (1995), pp. 1429–1439.

[16] S. Yun and Y. Lee, Robust estimation in mixed linear models with non-monotone missingness, Stat. Med. 25 (2006), pp. 3877–3892.

[17] H. Jacqmin-Gadda, S. Sibillot, C. Proust, J.-M. Molina, and R. Thiebaut, Robustness of the linear mixed model to

misspecified error distribution, Comput. Stat. Data Anal. 51 (2007), pp. 5142–5154.

[18] R.M. Royall and T.-S. Tsou, Interpreting statistical evidence using imperfect models: Robust adjusted likelihood

functions, J. R. Statist. Soc. Ser. B 65 (2003), pp. 391–404.

[19] T.-S. Tsou, Comparing two population means and variances – a parametric robust way, Comm. Stat. Theory Methods 32 (2003), pp. 2013–2019.

[20] T.-S. Tsou and K.-F. Cheng, Parametric robust regression analysis of contaminated data, Comm. Stat. Theory Methods 33 (2004), pp. 1887–1898.

[21] T.-S. Tsou, Robust Poisson regression, J. Statist. Plann. Inference 136 (2006), pp. 3173–3186. [22] H. White, Maximum likelihood estimation of misspecified models, Econometrica 50 (1982), pp. 1–25.

[23] R.W.M. Wedderburn, Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method, Biometrika 61 (1974), pp. 439–447.

[24] P. McCullagh and J.A. Nelder, Generalized Linear Models, Chapman and Hall, New York, 1989. [25] C.C. Heyde, Quasi-likelihood and Its Application, Springer-Verlag, New York, 1997.

[26] P.J. Huber, Robust Statistics, John Wiley, New York, 1981.

[27] P. McCullagh, Quasi-likelihood functions, Ann. Stat. 11 (1983), pp. 59–67.

[28] D.R. Cox and D.V. Hinkley, Theoretical Statistics, Chapman and Hall, New York, 1996.

(18)

Appendix

Here some details regarding the quantities required to calculate the adjusting matrices Anand Bnof AnB−1n Anthat correct

the normal regression model are provided.

To facilitate calculation of the adjusting matrices, let βt= (β1, β2, . . . , βw)= (γp−1, γp−2, . . . , γp−w)be the w-vector of

parameters of interest and let ϕ be the (p− w + 1)-dimensional nuisance parameters with the non-regression parameter

σ2 and the (p− w) regression coefficients (γ0, γ1, . . . , γp−w−1). Under the normal working model, Ihββand Ihβϕ are

approximately equal to the w× w matrix

1 Nσ2 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N i=1 (μi,0)2xi,p2−1 . . . N i=1

(μi,0)2xi,p−1xi,p−w

. . . . . . . . . N i=1(μ

i,0)2xi,p−wxi,p−1 . . . N i=1(μ i,0)2x 2 i,p−w ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

and the w× (p − w + 1) matrix

1 Nσ2 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 N i=1(μ

i,0)2xi,p−1xi,0 . . . N

i=1(μ

i,0)2xi,p−1xi,p−w−1

0

N

i=1(μ

i,0)2xi,p−2xi,0 . . . N

i=1(μ

i,0)2xi,p−2xi,p−w−1

. . . . . . . . . . . . 0 N i=1(μ

i,0)2xi,p−wxi,0 . . . N

i=1(μ

i,0)2xi,p−wxi,p−w−1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,

respectively. Ihϕϕis asymptotically expressed by

1 Nσ02 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ iϕϕ1 0 . . . 0 0 N i=1(μ i,0)2x2i,0 . . . N i=1(μ

i,0)2xi,0xi,p−w−1

. . . . . . . . . . . . 0 N i=1(μ

i,0)2xi,p−w−1xi,0 . . . N i=1(μ i,0)2x 2 i,p−w−1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , where iϕϕ1= −σ02lσ2_σ2.

For simplicity of notation, let Δn be the (p− w) × (p − w) matrix with the jth row as (Ni=1(μi,0)2xi,j−1xi,0, . . . ,

N

i=1(μi,0)2xi,j−1xi,p−w−1)for j= 1, 2, . . . , p − w. Then, Ihϕϕis approximately written in the form

Ihϕϕ≈ 1 Nσ02 iϕϕ1 0 0 Δn ,

where 0 is a (p− w)-vector and consists of only zeros. Its inverse I−1_hϕϕis approximately given by

I−1_hϕϕ≈Nσ 2 0 |Δn| ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ |Δn| iϕϕ1 0 . . . 0 0 R1,1 . . . Rp−w,1 . . . . . . . . . . . . 0 R1,p−w . . . Rp−w,p−w ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,

where Rl,m= (−1)l+m|M(lm)(Δn)| is the (l,m)th cofactor of Δn.

(19)

Similarly, Vhββ, Vhβϕand Vhϕϕare approximately equal to 1 Nσ4 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N i=1

Varh(Yi)(μi,0)2x2i,p−1 . . .

N

i=1

Varh(Yi)(μi,0)2xi,p−1xi,p−w

. . . . . . . . . N i=1Varh(Yi)(μ

i,0)2xi,p−wxi,p−1 . . .

N i=1Varh(Yi)(μ i,0)2x2i,p−w ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , 1 Nσ4 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 N i=1Varh(Yi)(μ

i,0)2xi,p−1xi,0 . . . N

i=1Varh(Yi)(μ

i,0)2xi,p−1xi,p−w−1

0

N

i=1Varh(Yi)(μ

i,0)2xi,p−2xi,0 . . . N

i=1Varh(Yi)(μ

i,0)2xi,p−2xi,p−w−1

. . . . . . . . . . . . 0 N i=1Varh(Yi)(μ

i,0)2xi,p−wxi,0 . . . N

i=1Varh(Yi)(μ

i,0)2xi,p−wxi,p−w−1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and 1 Nσ₀4 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ vϕϕ1 vσ2γ0 . . . vσ2γp−w−1 vγ0σ2 N i=1Varh(Yi)(μ i,0)2x2i,0 . . . N i=1Varh(Yi)(μ

i,0)2xi,0xi,p−w−1

. . . . . . . . . . . . v_γ_p_−w−1_σ2 N i=1Varh(Yi)(μ

i,0)2xi,p−w−1xi,0 . . . N i=1Varh(Yi)(μ i,0)2x2i,p−w−1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,

respectively. Here vϕϕ1= σ04Eh[lσ2(β₀, ϕ₀)l_σ2(β₀, ϕ₀)] and v_σ2_γ_j = v_γ_j_σ2= σ₀4Eh[lσ2(β₀, ϕ₀)lγj(β0, ϕ0)], j = 0, 1, . . . ,

p− w − 1.

According to Equation (1), the (u, v) entries of Anare derived as follows:

An(uv)= Ihβuβv− IhβuϕI−1hϕϕIhϕβv− IhβvϕI−1hϕϕIhϕβu+ IhβvϕI−1hϕϕIhϕϕI−1hϕϕIhϕβu = lim N→∞ 1 Nσ02 ⎧ ⎨ ⎩ N i=1 (μi,0) 2_x

i,p−uxi,p−v−_|Δ1

n| N i=1 (μi,0) 2 ⎛ ⎝xi,p−v p−w j=1 xi,j−1|Δnj(u)| ⎞ ⎠ − 1 |Δn| N i=1 (μi,0) 2 ⎛ ⎝xi,p−u p−w j=1 xi,j−1|Δnj(v)| ⎞ ⎠ + 1 |Δn|2 N i=1 (μi,0) 2 ⎛ ⎝p−w j=1 xi,j−1|Δnj(v)| ⎞ ⎠ ⎛ ⎝p−w j=1 xi,j−1|Δnj(u)| ⎞ ⎠ ⎫ ⎬ ⎭ = lim N→∞ 1 Nσ02 N i=1 (μi,0) 2 ⎛ ⎝xi,p−u− p−w j=1 |Δnj(u)| |Δn| xi,j−1 ⎞ ⎠ ⎛ ⎝xi,p−v− p−w j=1 |Δnj(v)| |Δn| xi,j−1 ⎞ ⎠ ,

where|Δnj(u)| =pm−w=1(Rj,mNi=1(μi,0)2xi,m−1xi,p−u).

(20)

Journal of Statistical Computation and Simulation 867 According to Equation (2), the (u, v) entries of Bnare derived as follows:

Bn(uv)= Vhβuβv− IhβuϕI−1hϕϕVhϕβv− VhβuϕI−1hϕϕIhϕβv+ IhβuϕI−1hϕϕVhϕϕI−1hϕϕIhϕβv = lim N→∞ 1 Nσ4 0 ⎧ ⎨ ⎩ N i=1

Varh(Yi)(μi,0)2xi,p−uxi,p−v−

1 |Δn| N i=1 Varh(Yi)(μi,0)2 ⎛ ⎝xi,p−v p−w j=1 xi,j−1|Δnj(u)| ⎞ ⎠ −_|Δ1 n| N i=1 Varh(Yi)(μi,0)2 ⎛ ⎝xi,p−u p−w j=1 xi,j−1|Δnj(v)| ⎞ ⎠ + _|Δ1 n|2 N i=1 Varh(Yi)(μi,0)2 ⎛ ⎝p−w j=1 xi,j−1|Δnj(u)| ⎞ ⎠ ⎛ ⎝p−w j=1 xi,j−1|Δnj(v)| ⎞ ⎠ ⎫ ⎬ ⎭ = lim N→∞ 1 Nσ4 0 N i=1 Varh(Yi)(μi,0)2 ⎛ ⎝xi,p−u− p−w j=1 |Δnj(u)| |Δn| xi,j−1 ⎞ ⎠ ⎛ ⎝xi,p−v− p−w j=1 |Δnj(v)| |Δn| xi,j−1 ⎞ ⎠ . The adjusting matrices Agand Bgof the gamma working model are derived in the analogous way.