- Email: [email protected]
Robust inference in an heteroscedastic measurement error model
Journal of the Korean Statistical Society 39 (2010) 439–447
Contents lists available at ScienceDirect
Journal of the Korean Statistical Society journal homepage: www.elsevier.com/locate/jkss
Robust inference in an heteroscedastic measurement error model Mário de Castro a,∗ , Manuel Galea b a
Universidade de São Paulo, Instituto de Ciências Matemáticas e de Computação, Caixa Postal 668, 13560-970, São Carlos-SP, Brazil
b
Universidad de Valparaíso, Av. Gran Bretaña 1091, Playa Ancha, Valparaíso, Chile
article
info
Article history: Received 25 April 2009 Accepted 27 September 2009 Available online 14 October 2009 AMS 2000 subject classifications: 62J99 62F99 Keywords: Errors-in-variables models Robust inference Student t distribution ECM algorithm
abstract In this paper we deal with robust inference in heteroscedastic measurement error models. Rather than the normal distribution, we postulate a Student t distribution for the observed variables. Maximum likelihood estimates are computed numerically. Consistent estimation of the asymptotic covariance matrices of the maximum likelihood and generalized least squares estimators is also discussed. Three test statistics are proposed for testing hypotheses of interest with the asymptotic chi-square distribution which guarantees correct asymptotic significance levels. Results of simulations and an application to a real data set are also reported. © 2009 The Korean Statistical Society. Published by Elsevier B.V. All rights reserved.
1. Introduction In this paper we deal with robust inference in heteroscedastic measurement error models. Rather than the normal distribution, we postulate a Student t distribution for the observed variables. The Student t distribution and its extensions have been successfully applied in data modeling (see, e. g., Jones, 2008, for a recent account). As typically considered in the literature (Lange, Little, & Taylor, 1989), the relevance of using the Student t distribution is related to its capability of downweighting influential observations. Heteroscedastic measurement error models with known error variances are a common setup in areas such as Epidemiology (Kulathinal, Kuulasmaa, & Gasbarra, 2002), Analytical Chemistry (Cheng & Riu, 2006), and Botany (Veenendaal et al., 2008), to name just a few. Estimates of these variances can be obtained by replications and are then treated as known. We embrace a functional model with known diagonal scale matrices varying between the observations. In a related work, Bolfarine and Arellano-Valle (1994) studied robust inference in an homoscedastic functional model. By combining the robustness provided by the Student t distribution and heteroscedasticity in the errors, our proposal is a robust counterpart to the normal heteroscedastic model in Ripley and Thompson (1987), Walter (1997), and Galea-Rojas, de Castilho, Bolfarine,and de Castro (2003), among others. The remaining of the paper is organized as follows. Section 2 comprises the methodological results. Maximum likelihood (ML) estimates are computed numerically. A simple graphical procedure for model checking is implemented. Consistent estimation of the asymptotic covariance matrix of the ML and generalized least squares (GLS) estimators is also discussed. Three test statistics are proposed for testing hypotheses of interest with the limiting chi-square distribution which guarantees correct asymptotic significance levels. Results of simulations and an application to a real data set are also reported in Sections 3 and 4, respectively. Concluding remarks are given in Section 5.
∗
Corresponding author. Tel.: +55 16 33739567; fax: +55 16 33739175. E-mail address: [email protected] (M. de Castro).
1226-3192/$ – see front matter © 2009 The Korean Statistical Society. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jkss.2009.09.003
440
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
2. Model and parameter estimation Let n be the sample size, Xi the observed value of the covariate in unit i, Yi the observed response in unit i, and xi the unobserved (true) covariate value for unit i. These variables are related by the model Yi = α + β xi + ei
and Xi = xi + ui ,
(1)
where ei and ui are contaminating errors due to diverse causes, i = 1, . . . , n. The unknown true values xi , i = 1, . . . , n, are modeled as constants. Since their number increases with the sample size, they are called incidental parameters. In this case, the above model is called a functional model (Cheng & Van Ness, 1999; Fuller, 1987). Letting Zi = (Yi , Xi )⊤ , instead of the normal distribution (Galea-Rojas et al., 2003; Ripley & Thompson, 1987), we assume that indep.
Zi ∼ t2 (µi , Λi , ν),
(2)
where t2 (µi , Λi , ν) denotes the bivariate Student t distribution with location vector µi = (α +β xi , xi ) , known scale matrix Λi = diag(λei , λui ), and ν degrees of freedom, i = 1, . . . , n. Homoscedastic versions of the model in (2) with unknown Λ were studied for instance by Bolfarine and Arellano-Valle (1994). The log-likelihood function corresponding to the model defined by (1) with assumption (2) can be written as ⊤
n
ℓ(θ, ν, x; Z ) = constant + (ν + 2) log ν −
n ν+2 −
2
2
log(ν + δi2 ),
(3)
i=1
with θ = (α, β)⊤ , x = (x1 , . . . , xn )⊤ , Z = (Z1 , . . . , Zn ), and
δi2 =
(Xi − xi )2 (Yi − α − β xi )2 + , λei λui
(4)
i = 1, . . . , n. After differentiating (3) with respect to xi and solving the resulting likelihood equation, it follows that the ML estimator of xi is given by
xi = xi (θ) =
λei Xi + βλui (Yi − α) . λei + β 2 λui
(5)
Taking into account that Yi − α − β xi =
βλui (Yi − α − β Xi ) λei (Yi − α − β Xi ) and Xi − xi = − , λei + β 2 λui λei + β 2 λui
by replacing xi for xi in (3) we arrive at the profile log-likelihood function n
ℓp (θ, ν; Z ) = constant + (ν + 2) log ν − 2
n ν+2 −
2
log ν +
i=1
(Yi − α − β Xi )2 λei + β 2 λui
.
(6)
Suppose for a while that ν is known. Direct maximization of the function in (6) is not a simple task. Alternatively, by using indep.
i.i.d.
1 −1 2 the stochastic representation Zi |qi ∼ N2 (µi , q− χν , with Z1 , . . . , Zn independent from q1 , . . . , qn , it i Λi ) and qi ∼ ν is possible to formulate an ECM algorithm (McLachlan & Krishnan, 1997; Meng & Rubin, 1993) as an extension of the EM algorithm in Bolfarine and Arellano-Valle (1994). Let Q = (q1 , . . . , qn )⊤ . It follows from the above stochastic representation that the complete data log-likelihood function is given by
ℓc (θ, x; Z , Q ) = constant − (m)
Let θ
n 1−
2 i =1
δ +
qi i2
ν 2
n − i =1
log qi −
n −
qi
.
(7)
i=1
and x(m) be the current estimates of θ and x, where m is the iteration counter. In the M step of the ECM algorithm we (m)
, x(m) ) = E{ℓc (θ, x; Z , Q ) | Z } evaluated . In the M step the estimates of θ and x are updated by maximizing the Q function. Since ν is
compute the expected complete data log-likelihood function given by Q(θ, x; θ (m)
at θ = θ
(m)
and x = x
indep.
2 known, the E step requires only the computation of E{qi | Zi }. Using the fact that qi | Zi ∼ (ν + δi2 )−1 χν+ 2 , we obtain 2 −1 E{qi | Zi } = (ν + δi ) (ν + 2). By replacing qi for the updated E{qi | Zi } in (7), i = 1, . . . , n, we arrive at the Q function. We can show that the Q function is maximized taking x as in (5). Next, estimates of θ = (α, β)⊤ are updated through a iterative weighted least squares scheme. The sketch of the ECM algorithm is as follows:
1. Let α (0) , β (0) , and x(0) be the starting values of α, β , and x, respectively. Initialize the iteration counter (m = 0); 2. Increment m by 1;
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
441
3. E step. Compute (m−1) 2 (m−1) 2 (m) {Yi − α (m−1) − β (m−1) xi } {Xi − xi } δi2 = + λei λui ν+2 (m) = ; qi (m) ν + δ2
and
i
(m)
(m−1)
α (m) = Y w − qi (m) /λei , i = 1, . . . , n, and 4. CM step 1. Compute the weights w i (m) = β (m−1) xw n ∑
(m) Yw
=
i =1 n
n ∑
w i (m) Yi and
∑ i=1
w i (m)
(m−1) xw
=
i=1
, where
(m−1) w i (m) xi n ∑ i =1
; w i (m)
5. CM step 2. Compute n ∑
(m) β
=
i =1
(m−1) {Yi − α (m) } w i (m) xi n ∑ i =1
; w i
(m) (m−1)2
xi
(m)
6. CM step 3. Compute xi using (5), i = 1, . . . , n; 7. Repeat steps 2 to 6 until convergence. Our stopping rule is based on the maximum relative difference between estimates in two successive iterations. Starting values α (0) , β (0) , and x(0) can be taken as the estimates of α , β , and x from the normal model (Galea-Rojas et al., 2003). According to Lange et al. (1989), ν can be interpreted as a robustness tuning parameter. Choosing a suitable set of values for ν , estimation of θ for each element in this set ( θ ν , say) can be achieved as described above. Then, the degrees of freedom, considered fixed and known, are such that ℓp (θ, ν; Z ) in (6) attains its maximum. ∑n Estimation of θ can also be accomplished by a GLS approach, which consists in minimizing i=1 (Yi − α − β Xi )2 /(λei + β 2 λui ). Since this method coincides with the ML method applied to the normal functional model (see, e. g., Cheng & Riu, 2006), estimation of θ is done as in Galea-Rojas et al. (2003) and it will be not repeated here. 2.1. Model checking Goodness of fit received much less attention than inference in the measurement error models literature. We describe the i.i.d.
implementation of a simple graphical device for model checking. Under the model in (2), we have that vi = δi2 /2 ∼ F2,ν , with δi2 as in (4). After applying the Wilson–Hilferty transformation (Johnson, Kotz, & Balakrishnan, 1994) 1/3
1 − 92ν vi ri =
2 9ν
2/3
vi
+
1 9
− 89 1/2 ,
(8)
i.i.d.
we obtain ri ∼ N (0, 1), i = 1, . . . , n, approximately. Such distributional result enables us to assess the fitting capability of the model in practice employing the simulated envelope proposed by Atkinson (1985). In order to implement this graphical tool for model checking, first we simulate J samples from (2) using the ML estimates θ and xi ( θ) in (5). For the jth simulated sample we compute the ML estimate θ and the transformed values rj1 , . . . , rjn from (8), which are ordered as rj(1) ≤ · · · ≤ rj(n) . The ordered pairs (Φ −1 (i − 3/8)/(n + 1/4) , r(i) ) are represented in a graph, where Φ −1 (·) denotes the quantile function of the standard normal distribution and ri is computed from (8) with the ML estimates θ . The limits of the ∑ J J J envelope, given by minj=1 rj(i) and maxj=1 rj(i) , and the line connecting the points (Φ −1 (i − 3/8)/(n + 1/4) , j=1 rj(i) /J ), i = 1, . . . , n, are also drawn in the graph. This plot form the basis to guide us on assessing departures from the postulated model. The technique will be illustrated in Section 4. 2.2. Covariance matrix estimators In this section we discuss the estimation of the covariance matrices of the ML and GLS estimators of θ . Maximization of the ∑n profile log-likelihood function (6) can be seen as the solution of the estimating equation expressed by i=1 Ψ i (Zi , θ) = 0,
442
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
where
∂ (ν + 2)(Yi − α − β Xi ) ℓpi ν(λei + β 2 λui ) + (Yi − α − β Xi )2 (9) Ψ i (Zi , θ) = ∂α = , ∂ ∂ ℓpi ℓpi xi (θ) ∂β ∂α with xi (θ) as in (5), i = 1, . . . , n. Under suitable regularity conditions (Mak, 1982), this estimating equation leads to a consistent and asymptotically normal estimator with mean θ and covariance matrix Ω(θ) = n−1 A(θ)−1 B(θ){A(θ)−1 }⊤ , where A(θ) = n−1
n − ∂ E
i=1
∂θ
Ψ (Z , θ) ⊤ i i
and B(θ) = n−1
n −
cov{Ψ i (Zi , θ)},
(10)
i=1
which can be estimated by A( θ) = n−1
n − ∂ i=1 n
B( θ) = n−1
−
∂θ ⊤
Ψ i (Zi , θ) and (11)
Ψ i (Zi , θ)Ψ i (Zi , θ) , ⊤
i=1
evaluated at θ = θ . For the sake of simplicity we omit the index n from these matrices. The elements of A(θ), represented by a11 , a12 = a21 , and a22 , are computed differentiating the components of (9), that conduces to
ν+2 δi2α δi2α 2 −δ , + iαα 2(ν + δi2 ) ν + δi2 i =1 n − δi2α δi2β ν+2 2 = −δiαβ + , 2(ν + δi2 ) ν + δi2 i =1 n − δi2β δi2β ν+2 2 = −δiββ + , 2(ν + δi2 ) ν + δi2 i =1
a11 =
a12
a22
n −
and
with δi2 given in (4), δi2α = −2(Yi − α − β Xi )/(λei + β 2 λui ), δi2αα = 2/(λei + β 2 λui ), δi2β = Xi δi2α − βλui δi2 δi2αα , δi2αβ = δi2αα (Xi − βλui δi2α ), and δi2ββ = Xi δi2αβ − βλui δi2β δi2αα − λui δi2 (δi2αα )2 (λei − β 2 λui )/2, i = 1, . . . , n. On the other side, we realize that the expectations and covariances in (10) can be readily computed as a consequence of the results presented by Vilca-Labra, Arellano-Valle, and Bolfarine (1998) (see also Vilca-Labra, 1996). We define ci =
1
λui
+
β2 , λei
x2i = xi 2 − 1
n − λei + β 2 λui D(θ) = i=1 ·
ν
,
ci (ν − 2) xi
and
λei + β 2 λui , x2 i
λei + β 2 λui
noticing that E(x2i ) = x2i , i = 1, . . . , n. After some algebraic manipulations the matrices in (10) turn out to be A(θ) =
B(θ) =
ν + 2 −1 n D(θ) and ν+3 2 (ν + 2) (ν + 1)(ν + 3)
n−1
0
D(θ) +
·
n − ν+3 1 . ν − 1 i=1 ci (λei + β 2 λui ) 0
(12)
In this way, considering (11) and (12) we have two estimators of the covariance matrix Ω(θ), denoted by Ωobs ( θ) and Ωexp ( θ), respectively. Furthermore, the covariance matrix of the GLS estimator of θ , denoted by ΩGLS ( θ), can be found
in Cheng and Riu (2006) (see also de Castro, Galea-Rojas, & Bolfarine, 2008). Approximate elliptical confidence regions for θ can be constructed employing these estimators. Hypotheses testing is the subject matter of the next section. The aforementioned regularity conditions aim to ensure, for instance, the application of the Ljapunoff’s central limit theorem. In our model we take ν > 4 and suppose that xi , λei , and λui in (2), i = 1, . . . , n, are bounded sequences.
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
443
Table 1 Mean, standard error (SE), estimated standard error (est. SE) and standard deviation of the estimated SE (sd est. SE) of the parameter estimates from 10 000 simulated samples according to different sample sizes and methods. True values are α = 1.5 and β = 2.0. Sample
Method
size (n)
α
β
Mean
SE
est. SE
sd est. SE
Mean
SE
est. SE
sd est. SE
40
GLS MLobs MLexp
1.503 1.503 1.503
0.1756 0.1596 0.1596
0.1504 0.1527 0.1568
0.04651 0.05853 0.009251
2.003 2.003 2.003
0.08596 0.07839 0.07839
0.07097 0.07445 0.07642
0.02123 0.03082 0.005053
80
GLS MLobs MLexp
1.502 1.501 1.501
0.1433 0.1287 0.1287
0.1330 0.1274 0.1277
0.02834 0.02378 0.005329
2.001 2.000 2.000
0.05981 0.05357 0.05357
0.05389 0.05275 0.05292
0.01282 0.01227 0.002393
160
GLS MLobs MLexp
1.501 1.500 1.500
0.07636 0.06901 0.06901
0.07062 0.06744 0.06803
0.01404 0.01534 0.001718
2.001 2.000 2.000
0.03680 0.03310 0.03310
0.03358 0.03233 0.03249
0.006671 0.006593 0.0009307
2.3. Hypotheses testing Now we deal with the problem of testing hypotheses like H: θ = θ 0 or H: β = β0 , where θ 0 and β0 are known. As a motivating example for the first hypothesis, in the methods comparison problem (Galea-Rojas et al., 2003; Riu & Rius, 1996) the main interest is to test if one method is unbiased with respect to the other, meaning that θ 0 = (0, 1)⊤ . Walter (1997) provides examples involving the hypothesis H: β = β0 , with β0 = 1. Such hypotheses can be tested using the Wald statistics W = ( θ − θ 0 )⊤ {Ω( θ)}−1 ( θ − θ 0 ) and W = ( β − β0 )2 /ω22 , where Ω( θ) indicates an estimator of Ω(θ) and ω22 is the (2,2) entry of Ω( θ). So, for each hypothesis three test statistics arise from the estimators in Section 2.2. They will be denoted by WGLS , Wobs , and Wexp . Under the hypotheses H: θ = θ 0 and H: β = β0 , the three statistics follow a chi-square distribution with two degrees and one degree of freedom, respectively, as n → ∞.
3. Simulation study Since in Section 2 we rely on asymptotic theory, the performance of the estimators and test statistics will be evaluated in a simulation study. Simulated scenarios mimic some characteristics of the data set in our example (Section 4). Sample sizes are 40, 80, and 160, whilst the distribution in (2) has ν = 5 degrees of freedom. The negative of the true unobservable xi , i = 1, . . . , n, are sampled from the gamma distribution with shape and rate parameters equal to 2.8 and 1.2, respectively. Then, for each generated xi we search for the index i∗ corresponding to the nearest observed value of the covariate X (see Fig. 3). The scale parameters for the ith observation are λei∗ and λui∗ taken from the data set in Section 4. The choice of the true θ will be clarified in the sequel. For each sample size, once xi , λei , and λui have been generated, we hold these values fixed throughout a batch of 10 000 simulated samples. The observations Yi and Xi follow the distribution in (2). Computations were performed with a specific purpose R program (R Development Core Team, 2009), which is available from the first author upon request. In Table 1 we present mean, standard error (SE), estimated standard error (est. SE), and standard deviation of the estimated SE (sd est. SE) of the parameter estimates from the simulated samples. Estimated standard errors are the mean of the square root of the diagonal elements of the ΩGLS ( θ), Ωobs ( θ), and Ωexp ( θ) matrices in Section 2.2. MLobs and MLexp methods mean that the standard errors were estimated using Ωobs (θ) and Ωexp ( θ), respectively, with θ replaced by the ML estimator ( θ). True vales are α = 1.5 and β = 2. Regardless the estimation method and the sample size, the means of the estimates are close to the true values. Not surprisingly, SE, est. SE, and sd est. SE decrease with the sample size. The GLS method yields estimates with greater SE than the ML method. Overall, the est. SE underestimates the SE, but they get closer to each other when the sample size increases, as expected. Moreover, SE and est. SE are in better agreement in the MLexp method than in the remaining methods. Estimates of the SE from the MLexp are more precise, as we can conclude by comparing the sd est. cells in Table 1. Completing our simulations we study the performance of the test statistics proposed in Section 2.3. For the hypothesis H: (α, β)⊤ = (0, 1)⊤ , Fig. 1 displays the chi-square (with two degrees of freedom) quantile–quantile (QQ) plots of the test statistics in samples of size 160. The vertical segments identify the critical values 5.99 and 9.21 associated to the significance levels 5% and 1%. From these plots we conclude that the Wexp statistic outperforms the WGLS and Wobs statistics. This assertion is further strengthened by the rejection rates in Table 2. When testing H: β = 1, the true value of α is −0.2. For both hypotheses in Table 2, whichever the sample size the Wexp statistic furnishes rejection rates with the best approximation to the nominal levels. In the simulations with n = 40 the rejection rates are distant from the nominal levels, even for the Wexp statistic, but there is a clear improvement of Wexp over WGLS and Wobs . Bearing in mind the results in Tables 1 and 2, as well as the QQ plots in Fig. 1, our simulations reveal that the ML estimation method and the Wexp test statistic seem to be the recommended tools for inference (at least within the scope of our study).
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
0
5
10
W
15
20
a
25
444
0
5
10
15
20
25
15 W 10 5 0
0
5
10
W
15
20
25
c
20
b
25
Theoretical χ22 quantiles
0
5
10
15
20
25
0
Theoretical χ22 quantiles
5
10
15
20
25
Theoretical χ22 quantiles
Fig. 1. Chi-square QQ plots of 10 000 simulated values of test statistics for the hypothesis H: (α, β)⊤ = (0, 1)⊤ in samples of size 160. (a) WGLS , (b) Wobs , and (c) Wexp . Table 2 Rejection rates (in %) of the hypotheses (a) H: (α, β)⊤ = (0, 1)⊤ and (b) H: β = 1 from 10 000 simulated samples according to different sample sizes, test statistics, and nominal significance levels. Sample size (n)
Test statistic
(a)
(b)
Significance level 1%
5%
1%
5%
40
WGLS Wobs Wexp
8.59 12.1 2.18
18.5 21.5 6.77
5.13 4.44 1.60
12.7 10.5 6.00
80
WGLS Wobs Wexp
3.01 3.17 1.40
9.53 9.28 5.85
2.88 2.03 1.34
9.31 7.42 5.62
160
WGLS Wobs Wexp
2.41 2.65 1.20
8.47 8.22 5.71
1.81 1.72 1.14
7.40 6.25 5.45
4. Application This section is dedicated to an example illustrating the proposed techniques. Table III in Walter (1997) shows metaanalysis data gathered from n = 35 studies of interventions to lower serum cholesterol and reduce total mortality. Response variable (L1 ) and covariate (L2 ) are the logit of the mortality rates in the experimental and control groups, respectively. In our example the scale parameters in (2) are the reciprocals of the study weights. Details about the number of treated participants in the studies and the computation of the weights are given by Walter (1997). Before any inferential analysis, we examine the Student t model as an alternative to the normal model fitted by Walter (1997). The profile log-likelihood function (6) is maximized with ν = 5 degrees of freedom, which seems to suggest that is worth to consider the Student t model. The
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
b
Sample values and simulated envelope –1 0 1 2
2
–2
0 –2
Sample values and simulated envelope
4
3
a
445
–2
–1
0
1
2
–2
Theoretical N(0,1) quantiles
–1
0
1
2
Theoretical N(0,1) quantiles
Fig. 2. Normal QQ plots of the transformed values in (8) and simulated envelopes. (a) Normal model and (b) Student t model with 5 degrees of freedom. Table 3 Parameter estimates and standard errors from different methods. Parameter
α β
Estimate
Standard error
GLS
ML
GLS
MLobs
MLexp
Bootstrap
−0.2040
−0.07459
0.09762 0.07044
0.07803 0.03793
0.07423 0.03706
0.1557 0.07730
0.8678
0.9660
envelopes in Fig. 2 were done with J = 100 samples. Since many transformed values from (8) lie outside the envelope in Fig. 2(a), whereas the sample points in Fig. 2(b) lie within the strip, the Student t model provides a more adequate fit to the data than the normal model, so that the latter will be excluded from the ongoing discussion. Notice that for the normal model, the transformed values are calculated with ν → ∞ in (8). Point estimates and its standard errors compose Table 3. Fig. 3 displays the sample points, square root of the scale parameters (crosses), and regression lines. Compared to the ML estimate, the GLS estimate of β induces a weaker relation between L1 and L2 . Table 3 also contains standard errors obtained via the bootstrap resampling method (Efron & Tibshirani, 1993) based on 10 000 data sets resampled from {(Yi , λei , Xi , λui ), i = 1, . . . , n}. Bootstrap standard errors are greater than the ones computed using asymptotic theory. Bernsen, Tasche, and Nagelkerke (1999) found similar results with the normal model and advocated that resampling methods should be used. According to Walter (1997), β − 1 measures the linear dependence between the variables. Applying the test statistics of Section 2.3 to the estimates in Table 2, the hypothesis H: β = 1 should not be rejected at the significance level of 1%. So, the linear dependence between L1 and L2 is not significant. Finally, we construct in Fig. 4 approximate 95% confidence regions for (α, β). GLS region is markedly different from the ML regions (see also Table 3). The hypothesis H: (α, β)⊤ = (0, 1)⊤ , interpreted as equivalence in the total mortality rates between the experimental and control groups, should not be rejected at a significance level of 5%. 5. Conclusion We have presented inferential procedures for a model with a wealth of applications in many areas. In particular, a simple graphical device for checking the model was implemented. In our example this device was useful in revealing the inaptness of the normal model. The simulations showed that the estimator of covariance matrix Ω(θ) using the matrices in (12) leads to better results when compared to estimator with the matrices in (11). The proposed test statistics can be easily extended to hypotheses other than the ones in Section 2.3, linear or nonlinear in θ . Since it is not possible to enunciate a general theoretical result regarding the merits of the estimation methods and test statistics, as a guidance to the practitioner having to select a statistical test, we suggest a choice based on simulations under
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
L1
–6
–5
–4
–3
–2
–1
0
446
–6
–5
–4
–3 L2
–2
–1
0
0.7
0.8
β
0.9
1.0
Fig. 3. Scatter plot of the logit of the mortality rates in the experimental (L1 ) and control (L2 ) groups, square root of the scale parameters (crosses), and regression lines: ML from the Student t model with 5 degrees of freedom (solid line) and GLS (dashed line).
–0.4
–0.3
–0.2
α
–0.1
0.0
0.1
Fig. 4. Approximate 95% confidence regions for (α, β) and the (0,1) point (): Student t models with 5 degrees of freedom (Ωexp ( θ): solid line, Ωobs ( θ): dashed line) and GLS method (dotted line).
conditions that resemble the problem at hand. Outputs like the ones in Fig. 1 and Tables 1 and 2 might be valuable in this respect. We assume that the scale matrix in (2) is known. Notwithstanding this limitation, such models have proven to be useful in many real-world applications of the normal model. Additional information (replicated observations, for instance), if available, enables us to propose a model in which the elements of the scale matrix are estimated. A Student t model might be envisioned as an alternative to the normal model in Dolby, Cormack, and Sinclair (1987). Such an extension constitutes an area for future research. Acknowledgements The authors would like to thank two anonymous referees for the comments and suggestions. This work is partially supported by FONDECYT (Fondo Nacional de Desarrollo Científico y Tecnológico, Chile, Proyecto #1070919).
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
447
References Atkinson, A. C. (1985). Plots, transformations, and regression. Oxford: Clarendon. Bernsen, R. M. D., Tasche, M. J. A., & Nagelkerke, N. J. D. (1999). Some notes on baseline risk and heterogeneity in meta-analysis. Statistics in Medicine, 18, 233–237. Bolfarine, H., & Arellano-Valle, R. B. (1994). Robust modelling in measurement error models using the t distribution. Brazilian Journal of Probability and Statistics, 8, 67–84. Cheng, C.-L., & Riu, J. (2006). On estimating linear relationships when both variables are subject to heteroscedastic measurement errors. Technometrics, 48, 511–519. Cheng, C.-L., & VanNess, J. W. (1999). Statistical regression with measurement error. London: Arnold. de Castro, M., Galea-Rojas, M., & Bolfarine, H. (2008). Letter to the editor. Technometrics, 50, 86–86. Dolby, G. R., Cormack, R. M., & Sinclair, D. F. (1987). On fitting bivariate functional relationships to unpaired and unequally replicated data. Biometrika, 74, 393–399. Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York: Chapman and Hall. Fuller, W. A. (1987). Measurement error models. New York: Wiley. Galea-Rojas, M., de Castilho, M. V., Bolfarine, H., & de Castro, M. (2003). Detection of analytical bias. The Analyst, 128, 1073–1081. Johnson, N. L., Kotz, S., & Balakrishnan, N. (1994). Continuous univariate distributions. Vol. 1 (2nd ed.). New York: Wiley. Jones, M. C. (2008). The t family and their close and distant relations. Journal of the Korean Statistical Society, 37, 293–302. Kulathinal, S. B., Kuulasmaa, K., & Gasbarra, D. (2002). Estimation of an errors-in-variables regression model when the variances of the measurement errors vary between the observations. Statistics in Medicine, 21, 1089–1101. Lange, K. L., Little, R. J. A., & Taylor, J. M. G. (1989). Robust statistical modeling using the t distribution. Journal of the American Statistical Association, 84, 881–896. Mak, T. K. (1982). Estimation in the presence of incidental parameters. Canadian Journal of Satistics, 10, 121–132. McLachlan, G. J., & Krishnan, T. (1997). The EM algorithm and extensions. New York: Wiley. Meng, X.-L., & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80, 267–278. R Development Core Team. (2009). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Ripley, B. D., & Thompson, M. (1987). Regression techniques for the detection of analytical bias. The Analyst, 112, 377–383. Riu, J., & Rius, F. X. (1996). Assessing the accuracy of analytical methods using linear regression with errors in both axes. Analytical Chemistry, 68, 1851–1857. Veenendaal, E. M., Mantlana, K. B., Pammenter, N. W., Weber, P., Huntsman-Mapila, P., & Lloyd, J. (2008). Growth form and seasonal variation in leaf gas exchange of Colophospermum mopane savanna trees in northwest Botswana. Tree Physiology, 28, 417–424. Vilca-Labra, F. E. (1996). Modelo de regressão elíptico com erros nas variáveis. Ph.D. thesis. Universidade de São Paulo – Instituto de Matemática e Estatística, São Paulo, Brasil (In Portuguese). Vilca-Labra, F., Arellano-Valle, R. B., & Bolfarine, H. (1998). Elliptical functional models. Journal of Multivariate Analysis, 65, 36–57. Walter, S. D. (1997). Variation in baseline risk as an explanation of heterogeneity in meta-analysis. Statistics in Medicine, 16, 2883–2900.
Contents lists available at ScienceDirect
Journal of the Korean Statistical Society journal homepage: www.elsevier.com/locate/jkss
Robust inference in an heteroscedastic measurement error model Mário de Castro a,∗ , Manuel Galea b a
Universidade de São Paulo, Instituto de Ciências Matemáticas e de Computação, Caixa Postal 668, 13560-970, São Carlos-SP, Brazil
b
Universidad de Valparaíso, Av. Gran Bretaña 1091, Playa Ancha, Valparaíso, Chile
article
info
Article history: Received 25 April 2009 Accepted 27 September 2009 Available online 14 October 2009 AMS 2000 subject classifications: 62J99 62F99 Keywords: Errors-in-variables models Robust inference Student t distribution ECM algorithm
abstract In this paper we deal with robust inference in heteroscedastic measurement error models. Rather than the normal distribution, we postulate a Student t distribution for the observed variables. Maximum likelihood estimates are computed numerically. Consistent estimation of the asymptotic covariance matrices of the maximum likelihood and generalized least squares estimators is also discussed. Three test statistics are proposed for testing hypotheses of interest with the asymptotic chi-square distribution which guarantees correct asymptotic significance levels. Results of simulations and an application to a real data set are also reported. © 2009 The Korean Statistical Society. Published by Elsevier B.V. All rights reserved.
1. Introduction In this paper we deal with robust inference in heteroscedastic measurement error models. Rather than the normal distribution, we postulate a Student t distribution for the observed variables. The Student t distribution and its extensions have been successfully applied in data modeling (see, e. g., Jones, 2008, for a recent account). As typically considered in the literature (Lange, Little, & Taylor, 1989), the relevance of using the Student t distribution is related to its capability of downweighting influential observations. Heteroscedastic measurement error models with known error variances are a common setup in areas such as Epidemiology (Kulathinal, Kuulasmaa, & Gasbarra, 2002), Analytical Chemistry (Cheng & Riu, 2006), and Botany (Veenendaal et al., 2008), to name just a few. Estimates of these variances can be obtained by replications and are then treated as known. We embrace a functional model with known diagonal scale matrices varying between the observations. In a related work, Bolfarine and Arellano-Valle (1994) studied robust inference in an homoscedastic functional model. By combining the robustness provided by the Student t distribution and heteroscedasticity in the errors, our proposal is a robust counterpart to the normal heteroscedastic model in Ripley and Thompson (1987), Walter (1997), and Galea-Rojas, de Castilho, Bolfarine,and de Castro (2003), among others. The remaining of the paper is organized as follows. Section 2 comprises the methodological results. Maximum likelihood (ML) estimates are computed numerically. A simple graphical procedure for model checking is implemented. Consistent estimation of the asymptotic covariance matrix of the ML and generalized least squares (GLS) estimators is also discussed. Three test statistics are proposed for testing hypotheses of interest with the limiting chi-square distribution which guarantees correct asymptotic significance levels. Results of simulations and an application to a real data set are also reported in Sections 3 and 4, respectively. Concluding remarks are given in Section 5.
∗
Corresponding author. Tel.: +55 16 33739567; fax: +55 16 33739175. E-mail address: [email protected] (M. de Castro).
1226-3192/$ – see front matter © 2009 The Korean Statistical Society. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jkss.2009.09.003
440
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
2. Model and parameter estimation Let n be the sample size, Xi the observed value of the covariate in unit i, Yi the observed response in unit i, and xi the unobserved (true) covariate value for unit i. These variables are related by the model Yi = α + β xi + ei
and Xi = xi + ui ,
(1)
where ei and ui are contaminating errors due to diverse causes, i = 1, . . . , n. The unknown true values xi , i = 1, . . . , n, are modeled as constants. Since their number increases with the sample size, they are called incidental parameters. In this case, the above model is called a functional model (Cheng & Van Ness, 1999; Fuller, 1987). Letting Zi = (Yi , Xi )⊤ , instead of the normal distribution (Galea-Rojas et al., 2003; Ripley & Thompson, 1987), we assume that indep.
Zi ∼ t2 (µi , Λi , ν),
(2)
where t2 (µi , Λi , ν) denotes the bivariate Student t distribution with location vector µi = (α +β xi , xi ) , known scale matrix Λi = diag(λei , λui ), and ν degrees of freedom, i = 1, . . . , n. Homoscedastic versions of the model in (2) with unknown Λ were studied for instance by Bolfarine and Arellano-Valle (1994). The log-likelihood function corresponding to the model defined by (1) with assumption (2) can be written as ⊤
n
ℓ(θ, ν, x; Z ) = constant + (ν + 2) log ν −
n ν+2 −
2
2
log(ν + δi2 ),
(3)
i=1
with θ = (α, β)⊤ , x = (x1 , . . . , xn )⊤ , Z = (Z1 , . . . , Zn ), and
δi2 =
(Xi − xi )2 (Yi − α − β xi )2 + , λei λui
(4)
i = 1, . . . , n. After differentiating (3) with respect to xi and solving the resulting likelihood equation, it follows that the ML estimator of xi is given by
xi = xi (θ) =
λei Xi + βλui (Yi − α) . λei + β 2 λui
(5)
Taking into account that Yi − α − β xi =
βλui (Yi − α − β Xi ) λei (Yi − α − β Xi ) and Xi − xi = − , λei + β 2 λui λei + β 2 λui
by replacing xi for xi in (3) we arrive at the profile log-likelihood function n
ℓp (θ, ν; Z ) = constant + (ν + 2) log ν − 2
n ν+2 −
2
log ν +
i=1
(Yi − α − β Xi )2 λei + β 2 λui
.
(6)
Suppose for a while that ν is known. Direct maximization of the function in (6) is not a simple task. Alternatively, by using indep.
i.i.d.
1 −1 2 the stochastic representation Zi |qi ∼ N2 (µi , q− χν , with Z1 , . . . , Zn independent from q1 , . . . , qn , it i Λi ) and qi ∼ ν is possible to formulate an ECM algorithm (McLachlan & Krishnan, 1997; Meng & Rubin, 1993) as an extension of the EM algorithm in Bolfarine and Arellano-Valle (1994). Let Q = (q1 , . . . , qn )⊤ . It follows from the above stochastic representation that the complete data log-likelihood function is given by
ℓc (θ, x; Z , Q ) = constant − (m)
Let θ
n 1−
2 i =1
δ +
qi i2
ν 2
n − i =1
log qi −
n −
qi
.
(7)
i=1
and x(m) be the current estimates of θ and x, where m is the iteration counter. In the M step of the ECM algorithm we (m)
, x(m) ) = E{ℓc (θ, x; Z , Q ) | Z } evaluated . In the M step the estimates of θ and x are updated by maximizing the Q function. Since ν is
compute the expected complete data log-likelihood function given by Q(θ, x; θ (m)
at θ = θ
(m)
and x = x
indep.
2 known, the E step requires only the computation of E{qi | Zi }. Using the fact that qi | Zi ∼ (ν + δi2 )−1 χν+ 2 , we obtain 2 −1 E{qi | Zi } = (ν + δi ) (ν + 2). By replacing qi for the updated E{qi | Zi } in (7), i = 1, . . . , n, we arrive at the Q function. We can show that the Q function is maximized taking x as in (5). Next, estimates of θ = (α, β)⊤ are updated through a iterative weighted least squares scheme. The sketch of the ECM algorithm is as follows:
1. Let α (0) , β (0) , and x(0) be the starting values of α, β , and x, respectively. Initialize the iteration counter (m = 0); 2. Increment m by 1;
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
441
3. E step. Compute (m−1) 2 (m−1) 2 (m) {Yi − α (m−1) − β (m−1) xi } {Xi − xi } δi2 = + λei λui ν+2 (m) = ; qi (m) ν + δ2
and
i
(m)
(m−1)
α (m) = Y w − qi (m) /λei , i = 1, . . . , n, and 4. CM step 1. Compute the weights w i (m) = β (m−1) xw n ∑
(m) Yw
=
i =1 n
n ∑
w i (m) Yi and
∑ i=1
w i (m)
(m−1) xw
=
i=1
, where
(m−1) w i (m) xi n ∑ i =1
; w i (m)
5. CM step 2. Compute n ∑
(m) β
=
i =1
(m−1) {Yi − α (m) } w i (m) xi n ∑ i =1
; w i
(m) (m−1)2
xi
(m)
6. CM step 3. Compute xi using (5), i = 1, . . . , n; 7. Repeat steps 2 to 6 until convergence. Our stopping rule is based on the maximum relative difference between estimates in two successive iterations. Starting values α (0) , β (0) , and x(0) can be taken as the estimates of α , β , and x from the normal model (Galea-Rojas et al., 2003). According to Lange et al. (1989), ν can be interpreted as a robustness tuning parameter. Choosing a suitable set of values for ν , estimation of θ for each element in this set ( θ ν , say) can be achieved as described above. Then, the degrees of freedom, considered fixed and known, are such that ℓp (θ, ν; Z ) in (6) attains its maximum. ∑n Estimation of θ can also be accomplished by a GLS approach, which consists in minimizing i=1 (Yi − α − β Xi )2 /(λei + β 2 λui ). Since this method coincides with the ML method applied to the normal functional model (see, e. g., Cheng & Riu, 2006), estimation of θ is done as in Galea-Rojas et al. (2003) and it will be not repeated here. 2.1. Model checking Goodness of fit received much less attention than inference in the measurement error models literature. We describe the i.i.d.
implementation of a simple graphical device for model checking. Under the model in (2), we have that vi = δi2 /2 ∼ F2,ν , with δi2 as in (4). After applying the Wilson–Hilferty transformation (Johnson, Kotz, & Balakrishnan, 1994) 1/3
1 − 92ν vi ri =
2 9ν
2/3
vi
+
1 9
− 89 1/2 ,
(8)
i.i.d.
we obtain ri ∼ N (0, 1), i = 1, . . . , n, approximately. Such distributional result enables us to assess the fitting capability of the model in practice employing the simulated envelope proposed by Atkinson (1985). In order to implement this graphical tool for model checking, first we simulate J samples from (2) using the ML estimates θ and xi ( θ) in (5). For the jth simulated sample we compute the ML estimate θ and the transformed values rj1 , . . . , rjn from (8), which are ordered as rj(1) ≤ · · · ≤ rj(n) . The ordered pairs (Φ −1 (i − 3/8)/(n + 1/4) , r(i) ) are represented in a graph, where Φ −1 (·) denotes the quantile function of the standard normal distribution and ri is computed from (8) with the ML estimates θ . The limits of the ∑ J J J envelope, given by minj=1 rj(i) and maxj=1 rj(i) , and the line connecting the points (Φ −1 (i − 3/8)/(n + 1/4) , j=1 rj(i) /J ), i = 1, . . . , n, are also drawn in the graph. This plot form the basis to guide us on assessing departures from the postulated model. The technique will be illustrated in Section 4. 2.2. Covariance matrix estimators In this section we discuss the estimation of the covariance matrices of the ML and GLS estimators of θ . Maximization of the ∑n profile log-likelihood function (6) can be seen as the solution of the estimating equation expressed by i=1 Ψ i (Zi , θ) = 0,
442
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
where
∂ (ν + 2)(Yi − α − β Xi ) ℓpi ν(λei + β 2 λui ) + (Yi − α − β Xi )2 (9) Ψ i (Zi , θ) = ∂α = , ∂ ∂ ℓpi ℓpi xi (θ) ∂β ∂α with xi (θ) as in (5), i = 1, . . . , n. Under suitable regularity conditions (Mak, 1982), this estimating equation leads to a consistent and asymptotically normal estimator with mean θ and covariance matrix Ω(θ) = n−1 A(θ)−1 B(θ){A(θ)−1 }⊤ , where A(θ) = n−1
n − ∂ E
i=1
∂θ
Ψ (Z , θ) ⊤ i i
and B(θ) = n−1
n −
cov{Ψ i (Zi , θ)},
(10)
i=1
which can be estimated by A( θ) = n−1
n − ∂ i=1 n
B( θ) = n−1
−
∂θ ⊤
Ψ i (Zi , θ) and (11)
Ψ i (Zi , θ)Ψ i (Zi , θ) , ⊤
i=1
evaluated at θ = θ . For the sake of simplicity we omit the index n from these matrices. The elements of A(θ), represented by a11 , a12 = a21 , and a22 , are computed differentiating the components of (9), that conduces to
ν+2 δi2α δi2α 2 −δ , + iαα 2(ν + δi2 ) ν + δi2 i =1 n − δi2α δi2β ν+2 2 = −δiαβ + , 2(ν + δi2 ) ν + δi2 i =1 n − δi2β δi2β ν+2 2 = −δiββ + , 2(ν + δi2 ) ν + δi2 i =1
a11 =
a12
a22
n −
and
with δi2 given in (4), δi2α = −2(Yi − α − β Xi )/(λei + β 2 λui ), δi2αα = 2/(λei + β 2 λui ), δi2β = Xi δi2α − βλui δi2 δi2αα , δi2αβ = δi2αα (Xi − βλui δi2α ), and δi2ββ = Xi δi2αβ − βλui δi2β δi2αα − λui δi2 (δi2αα )2 (λei − β 2 λui )/2, i = 1, . . . , n. On the other side, we realize that the expectations and covariances in (10) can be readily computed as a consequence of the results presented by Vilca-Labra, Arellano-Valle, and Bolfarine (1998) (see also Vilca-Labra, 1996). We define ci =
1
λui
+
β2 , λei
x2i = xi 2 − 1
n − λei + β 2 λui D(θ) = i=1 ·
ν
,
ci (ν − 2) xi
and
λei + β 2 λui , x2 i
λei + β 2 λui
noticing that E(x2i ) = x2i , i = 1, . . . , n. After some algebraic manipulations the matrices in (10) turn out to be A(θ) =
B(θ) =
ν + 2 −1 n D(θ) and ν+3 2 (ν + 2) (ν + 1)(ν + 3)
n−1
0
D(θ) +
·
n − ν+3 1 . ν − 1 i=1 ci (λei + β 2 λui ) 0
(12)
In this way, considering (11) and (12) we have two estimators of the covariance matrix Ω(θ), denoted by Ωobs ( θ) and Ωexp ( θ), respectively. Furthermore, the covariance matrix of the GLS estimator of θ , denoted by ΩGLS ( θ), can be found
in Cheng and Riu (2006) (see also de Castro, Galea-Rojas, & Bolfarine, 2008). Approximate elliptical confidence regions for θ can be constructed employing these estimators. Hypotheses testing is the subject matter of the next section. The aforementioned regularity conditions aim to ensure, for instance, the application of the Ljapunoff’s central limit theorem. In our model we take ν > 4 and suppose that xi , λei , and λui in (2), i = 1, . . . , n, are bounded sequences.
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
443
Table 1 Mean, standard error (SE), estimated standard error (est. SE) and standard deviation of the estimated SE (sd est. SE) of the parameter estimates from 10 000 simulated samples according to different sample sizes and methods. True values are α = 1.5 and β = 2.0. Sample
Method
size (n)
α
β
Mean
SE
est. SE
sd est. SE
Mean
SE
est. SE
sd est. SE
40
GLS MLobs MLexp
1.503 1.503 1.503
0.1756 0.1596 0.1596
0.1504 0.1527 0.1568
0.04651 0.05853 0.009251
2.003 2.003 2.003
0.08596 0.07839 0.07839
0.07097 0.07445 0.07642
0.02123 0.03082 0.005053
80
GLS MLobs MLexp
1.502 1.501 1.501
0.1433 0.1287 0.1287
0.1330 0.1274 0.1277
0.02834 0.02378 0.005329
2.001 2.000 2.000
0.05981 0.05357 0.05357
0.05389 0.05275 0.05292
0.01282 0.01227 0.002393
160
GLS MLobs MLexp
1.501 1.500 1.500
0.07636 0.06901 0.06901
0.07062 0.06744 0.06803
0.01404 0.01534 0.001718
2.001 2.000 2.000
0.03680 0.03310 0.03310
0.03358 0.03233 0.03249
0.006671 0.006593 0.0009307
2.3. Hypotheses testing Now we deal with the problem of testing hypotheses like H: θ = θ 0 or H: β = β0 , where θ 0 and β0 are known. As a motivating example for the first hypothesis, in the methods comparison problem (Galea-Rojas et al., 2003; Riu & Rius, 1996) the main interest is to test if one method is unbiased with respect to the other, meaning that θ 0 = (0, 1)⊤ . Walter (1997) provides examples involving the hypothesis H: β = β0 , with β0 = 1. Such hypotheses can be tested using the Wald statistics W = ( θ − θ 0 )⊤ {Ω( θ)}−1 ( θ − θ 0 ) and W = ( β − β0 )2 /ω22 , where Ω( θ) indicates an estimator of Ω(θ) and ω22 is the (2,2) entry of Ω( θ). So, for each hypothesis three test statistics arise from the estimators in Section 2.2. They will be denoted by WGLS , Wobs , and Wexp . Under the hypotheses H: θ = θ 0 and H: β = β0 , the three statistics follow a chi-square distribution with two degrees and one degree of freedom, respectively, as n → ∞.
3. Simulation study Since in Section 2 we rely on asymptotic theory, the performance of the estimators and test statistics will be evaluated in a simulation study. Simulated scenarios mimic some characteristics of the data set in our example (Section 4). Sample sizes are 40, 80, and 160, whilst the distribution in (2) has ν = 5 degrees of freedom. The negative of the true unobservable xi , i = 1, . . . , n, are sampled from the gamma distribution with shape and rate parameters equal to 2.8 and 1.2, respectively. Then, for each generated xi we search for the index i∗ corresponding to the nearest observed value of the covariate X (see Fig. 3). The scale parameters for the ith observation are λei∗ and λui∗ taken from the data set in Section 4. The choice of the true θ will be clarified in the sequel. For each sample size, once xi , λei , and λui have been generated, we hold these values fixed throughout a batch of 10 000 simulated samples. The observations Yi and Xi follow the distribution in (2). Computations were performed with a specific purpose R program (R Development Core Team, 2009), which is available from the first author upon request. In Table 1 we present mean, standard error (SE), estimated standard error (est. SE), and standard deviation of the estimated SE (sd est. SE) of the parameter estimates from the simulated samples. Estimated standard errors are the mean of the square root of the diagonal elements of the ΩGLS ( θ), Ωobs ( θ), and Ωexp ( θ) matrices in Section 2.2. MLobs and MLexp methods mean that the standard errors were estimated using Ωobs (θ) and Ωexp ( θ), respectively, with θ replaced by the ML estimator ( θ). True vales are α = 1.5 and β = 2. Regardless the estimation method and the sample size, the means of the estimates are close to the true values. Not surprisingly, SE, est. SE, and sd est. SE decrease with the sample size. The GLS method yields estimates with greater SE than the ML method. Overall, the est. SE underestimates the SE, but they get closer to each other when the sample size increases, as expected. Moreover, SE and est. SE are in better agreement in the MLexp method than in the remaining methods. Estimates of the SE from the MLexp are more precise, as we can conclude by comparing the sd est. cells in Table 1. Completing our simulations we study the performance of the test statistics proposed in Section 2.3. For the hypothesis H: (α, β)⊤ = (0, 1)⊤ , Fig. 1 displays the chi-square (with two degrees of freedom) quantile–quantile (QQ) plots of the test statistics in samples of size 160. The vertical segments identify the critical values 5.99 and 9.21 associated to the significance levels 5% and 1%. From these plots we conclude that the Wexp statistic outperforms the WGLS and Wobs statistics. This assertion is further strengthened by the rejection rates in Table 2. When testing H: β = 1, the true value of α is −0.2. For both hypotheses in Table 2, whichever the sample size the Wexp statistic furnishes rejection rates with the best approximation to the nominal levels. In the simulations with n = 40 the rejection rates are distant from the nominal levels, even for the Wexp statistic, but there is a clear improvement of Wexp over WGLS and Wobs . Bearing in mind the results in Tables 1 and 2, as well as the QQ plots in Fig. 1, our simulations reveal that the ML estimation method and the Wexp test statistic seem to be the recommended tools for inference (at least within the scope of our study).
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
0
5
10
W
15
20
a
25
444
0
5
10
15
20
25
15 W 10 5 0
0
5
10
W
15
20
25
c
20
b
25
Theoretical χ22 quantiles
0
5
10
15
20
25
0
Theoretical χ22 quantiles
5
10
15
20
25
Theoretical χ22 quantiles
Fig. 1. Chi-square QQ plots of 10 000 simulated values of test statistics for the hypothesis H: (α, β)⊤ = (0, 1)⊤ in samples of size 160. (a) WGLS , (b) Wobs , and (c) Wexp . Table 2 Rejection rates (in %) of the hypotheses (a) H: (α, β)⊤ = (0, 1)⊤ and (b) H: β = 1 from 10 000 simulated samples according to different sample sizes, test statistics, and nominal significance levels. Sample size (n)
Test statistic
(a)
(b)
Significance level 1%
5%
1%
5%
40
WGLS Wobs Wexp
8.59 12.1 2.18
18.5 21.5 6.77
5.13 4.44 1.60
12.7 10.5 6.00
80
WGLS Wobs Wexp
3.01 3.17 1.40
9.53 9.28 5.85
2.88 2.03 1.34
9.31 7.42 5.62
160
WGLS Wobs Wexp
2.41 2.65 1.20
8.47 8.22 5.71
1.81 1.72 1.14
7.40 6.25 5.45
4. Application This section is dedicated to an example illustrating the proposed techniques. Table III in Walter (1997) shows metaanalysis data gathered from n = 35 studies of interventions to lower serum cholesterol and reduce total mortality. Response variable (L1 ) and covariate (L2 ) are the logit of the mortality rates in the experimental and control groups, respectively. In our example the scale parameters in (2) are the reciprocals of the study weights. Details about the number of treated participants in the studies and the computation of the weights are given by Walter (1997). Before any inferential analysis, we examine the Student t model as an alternative to the normal model fitted by Walter (1997). The profile log-likelihood function (6) is maximized with ν = 5 degrees of freedom, which seems to suggest that is worth to consider the Student t model. The
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
b
Sample values and simulated envelope –1 0 1 2
2
–2
0 –2
Sample values and simulated envelope
4
3
a
445
–2
–1
0
1
2
–2
Theoretical N(0,1) quantiles
–1
0
1
2
Theoretical N(0,1) quantiles
Fig. 2. Normal QQ plots of the transformed values in (8) and simulated envelopes. (a) Normal model and (b) Student t model with 5 degrees of freedom. Table 3 Parameter estimates and standard errors from different methods. Parameter
α β
Estimate
Standard error
GLS
ML
GLS
MLobs
MLexp
Bootstrap
−0.2040
−0.07459
0.09762 0.07044
0.07803 0.03793
0.07423 0.03706
0.1557 0.07730
0.8678
0.9660
envelopes in Fig. 2 were done with J = 100 samples. Since many transformed values from (8) lie outside the envelope in Fig. 2(a), whereas the sample points in Fig. 2(b) lie within the strip, the Student t model provides a more adequate fit to the data than the normal model, so that the latter will be excluded from the ongoing discussion. Notice that for the normal model, the transformed values are calculated with ν → ∞ in (8). Point estimates and its standard errors compose Table 3. Fig. 3 displays the sample points, square root of the scale parameters (crosses), and regression lines. Compared to the ML estimate, the GLS estimate of β induces a weaker relation between L1 and L2 . Table 3 also contains standard errors obtained via the bootstrap resampling method (Efron & Tibshirani, 1993) based on 10 000 data sets resampled from {(Yi , λei , Xi , λui ), i = 1, . . . , n}. Bootstrap standard errors are greater than the ones computed using asymptotic theory. Bernsen, Tasche, and Nagelkerke (1999) found similar results with the normal model and advocated that resampling methods should be used. According to Walter (1997), β − 1 measures the linear dependence between the variables. Applying the test statistics of Section 2.3 to the estimates in Table 2, the hypothesis H: β = 1 should not be rejected at the significance level of 1%. So, the linear dependence between L1 and L2 is not significant. Finally, we construct in Fig. 4 approximate 95% confidence regions for (α, β). GLS region is markedly different from the ML regions (see also Table 3). The hypothesis H: (α, β)⊤ = (0, 1)⊤ , interpreted as equivalence in the total mortality rates between the experimental and control groups, should not be rejected at a significance level of 5%. 5. Conclusion We have presented inferential procedures for a model with a wealth of applications in many areas. In particular, a simple graphical device for checking the model was implemented. In our example this device was useful in revealing the inaptness of the normal model. The simulations showed that the estimator of covariance matrix Ω(θ) using the matrices in (12) leads to better results when compared to estimator with the matrices in (11). The proposed test statistics can be easily extended to hypotheses other than the ones in Section 2.3, linear or nonlinear in θ . Since it is not possible to enunciate a general theoretical result regarding the merits of the estimation methods and test statistics, as a guidance to the practitioner having to select a statistical test, we suggest a choice based on simulations under
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
L1
–6
–5
–4
–3
–2
–1
0
446
–6
–5
–4
–3 L2
–2
–1
0
0.7
0.8
β
0.9
1.0
Fig. 3. Scatter plot of the logit of the mortality rates in the experimental (L1 ) and control (L2 ) groups, square root of the scale parameters (crosses), and regression lines: ML from the Student t model with 5 degrees of freedom (solid line) and GLS (dashed line).
–0.4
–0.3
–0.2
α
–0.1
0.0
0.1
Fig. 4. Approximate 95% confidence regions for (α, β) and the (0,1) point (): Student t models with 5 degrees of freedom (Ωexp ( θ): solid line, Ωobs ( θ): dashed line) and GLS method (dotted line).
conditions that resemble the problem at hand. Outputs like the ones in Fig. 1 and Tables 1 and 2 might be valuable in this respect. We assume that the scale matrix in (2) is known. Notwithstanding this limitation, such models have proven to be useful in many real-world applications of the normal model. Additional information (replicated observations, for instance), if available, enables us to propose a model in which the elements of the scale matrix are estimated. A Student t model might be envisioned as an alternative to the normal model in Dolby, Cormack, and Sinclair (1987). Such an extension constitutes an area for future research. Acknowledgements The authors would like to thank two anonymous referees for the comments and suggestions. This work is partially supported by FONDECYT (Fondo Nacional de Desarrollo Científico y Tecnológico, Chile, Proyecto #1070919).
M. de Castro, M. Galea / Journal of the Korean Statistical Society 39 (2010) 439–447
447
References Atkinson, A. C. (1985). Plots, transformations, and regression. Oxford: Clarendon. Bernsen, R. M. D., Tasche, M. J. A., & Nagelkerke, N. J. D. (1999). Some notes on baseline risk and heterogeneity in meta-analysis. Statistics in Medicine, 18, 233–237. Bolfarine, H., & Arellano-Valle, R. B. (1994). Robust modelling in measurement error models using the t distribution. Brazilian Journal of Probability and Statistics, 8, 67–84. Cheng, C.-L., & Riu, J. (2006). On estimating linear relationships when both variables are subject to heteroscedastic measurement errors. Technometrics, 48, 511–519. Cheng, C.-L., & VanNess, J. W. (1999). Statistical regression with measurement error. London: Arnold. de Castro, M., Galea-Rojas, M., & Bolfarine, H. (2008). Letter to the editor. Technometrics, 50, 86–86. Dolby, G. R., Cormack, R. M., & Sinclair, D. F. (1987). On fitting bivariate functional relationships to unpaired and unequally replicated data. Biometrika, 74, 393–399. Efron, B., & Tibshirani, R. J. (1993). An introduction to the bootstrap. New York: Chapman and Hall. Fuller, W. A. (1987). Measurement error models. New York: Wiley. Galea-Rojas, M., de Castilho, M. V., Bolfarine, H., & de Castro, M. (2003). Detection of analytical bias. The Analyst, 128, 1073–1081. Johnson, N. L., Kotz, S., & Balakrishnan, N. (1994). Continuous univariate distributions. Vol. 1 (2nd ed.). New York: Wiley. Jones, M. C. (2008). The t family and their close and distant relations. Journal of the Korean Statistical Society, 37, 293–302. Kulathinal, S. B., Kuulasmaa, K., & Gasbarra, D. (2002). Estimation of an errors-in-variables regression model when the variances of the measurement errors vary between the observations. Statistics in Medicine, 21, 1089–1101. Lange, K. L., Little, R. J. A., & Taylor, J. M. G. (1989). Robust statistical modeling using the t distribution. Journal of the American Statistical Association, 84, 881–896. Mak, T. K. (1982). Estimation in the presence of incidental parameters. Canadian Journal of Satistics, 10, 121–132. McLachlan, G. J., & Krishnan, T. (1997). The EM algorithm and extensions. New York: Wiley. Meng, X.-L., & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80, 267–278. R Development Core Team. (2009). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Ripley, B. D., & Thompson, M. (1987). Regression techniques for the detection of analytical bias. The Analyst, 112, 377–383. Riu, J., & Rius, F. X. (1996). Assessing the accuracy of analytical methods using linear regression with errors in both axes. Analytical Chemistry, 68, 1851–1857. Veenendaal, E. M., Mantlana, K. B., Pammenter, N. W., Weber, P., Huntsman-Mapila, P., & Lloyd, J. (2008). Growth form and seasonal variation in leaf gas exchange of Colophospermum mopane savanna trees in northwest Botswana. Tree Physiology, 28, 417–424. Vilca-Labra, F. E. (1996). Modelo de regressão elíptico com erros nas variáveis. Ph.D. thesis. Universidade de São Paulo – Instituto de Matemática e Estatística, São Paulo, Brasil (In Portuguese). Vilca-Labra, F., Arellano-Valle, R. B., & Bolfarine, H. (1998). Elliptical functional models. Journal of Multivariate Analysis, 65, 36–57. Walter, S. D. (1997). Variation in baseline risk as an explanation of heterogeneity in meta-analysis. Statistics in Medicine, 16, 2883–2900.