# A Variance Shift Model for Detection of Outliers in the Linear Measurement Error Model

Authors: Babak Babadi; Abdolrahman Rasekh; Ali Akbar Rasekhi; Karim Zare; Mohammad Reza Zadkarami
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We present a variance shift model for a linear measurement error model using the corrected likelihood of Nakamura (1990). This model assumes that a single outlier arises from an observation with inflated variance. The corrected likelihood ratio and the score test statistics are proposed to determine whether theith observation has an inflated variance. A parametric bootstrap procedure is used to obtain empirical distributions of the test statistics and a simulation study has been used to show the performance of proposed tests. Finally, a real data example is given for illustration.


---

## Body

## 1. Introduction

Outliers are observations that appear inconsistent with the rest of the data set and can have a profound destructive influence on the statistical analysis. To detect these kinds of observations in the linear models, different approaches have been suggested, among those one can refer to the case-deletion and variance shift models. The first approach is based on the assumption that outliers result from a shift in the mean of contaminated observations (see Barnett and Lewis [1] or Weisberg [2]) and the second procedure takes into account the assumption that an outlier arises from an error term with an increased variance (Cook and Weisberg [3]). Cook et al. [4] indicated that the maximum likelihood estimates for the position of the outlier under two methods could be different, unless the largest absolute studentized residual corresponds to the largest absolute residual. Using the residual maximum likelihood (REML), Thompson [5] showed that the residual variance and outlier position are the same under both methods.In the linear mixed models, case deletion method, variance shift outlier model, and related diagnostics are studied widely by different authors. Christensen et al. [6] presented case deletion diagnostics for both fixed effects and variance components models. Banerjee and Frees [7] introduced case deletion diagnostics for both fixed effects and random subject effect in linear longitudinal models. Xuping and Bocheng [8] presented a unified diagnostic method for linear mixed models based upon the joint likelihood given by Robinson [9]. They showed that the estimates of parameters in case deletion method are equivalent to those in mean shift outlier model. Haslett and Dillane [10] proved a “delete = replace” identity in linear models and applied it to deletion diagnostics for estimators of variance components. Zewotir and Galpin [11] provided routine diagnostic tools for fixed effects, random effects, and variance components, which are computationally inexpensive. Li et al. [12] considered subset deletion diagnostics for fixed effects, random effects and one variance component in varying coefficient mixed models. Gumedze et al. [13] extended the variance shift outlier model (VSOM) to the linear mixed model.In linear regression models, independent variables are often susceptible to nonnegligible errors, and then it will be more appropriate to consider the measurement error models (see Fuller [14] and Stefanski [15]). However, in measurement error models ordinary maximum likelihood (ML) estimates lose the consistency. In order to correct the bias on parameter estimation, a method in which the score function itself is corrected for measurement errors, is available for the estimation of parameters. This method is based on the corrected log-likelihood of Nakamura [16] (see also Giménez and Bolfarine [17] for more discussion).On diagnostic methods for measurement error models, some previous works are due to Kelly [18] and Wellman and Gunst [19]. Zhong et al. [20] obtained case deletion and mean shift outlier model for linear measurement error models based upon the corrected likelihood of Nakamura [16]. Rasekh [21] studied multiple outlier detection in multivariate functional measurement error models based on the suitable definition of standardized residuals. Giménez and Galea [22] studied influence measures on corrected score estimators in functional heteroscedastic measurement error models and a local influence study on functional comparative calibration models with replicated data is developed by Giménez and Patat [23].In this paper, we concentrate on the variance shift model of Cook et al. [4], for the linear measurement error model, using the corrected likelihood [16]. In Section 2, we present the basis of the corrected score method and obtain the estimates of parameters of the model. In Section 3, a variance shift model for the linear measurement error model is derived and the joint corrected maximum likelihood estimates are characterized. In Section 4, we develop the likelihood ratio and the score test statistics. Furthermore, a parametric bootstrap procedure is used to generate the empirical distribution of these statistics. In Section 5, to verify the performance of the proposed test statistics, a simulation study is reported and finally, an illustrative example is given in Section 6.


## 2. Corrected Log-Likelihood of Measurement Error Models

Consider the linear measurement error models:(1)
y
=
Z
β
+
ε
,
ε
~
N
0
,
σ
2
I
n
,
X
=
Z
+
U
,
U
~
N
0
,
I
n
⊗
Λ
,
where y
=
(
y
1
,
y
2
,
…
,
y
n
)
′, Z is an n
×
pmatrix of unobservable regressors, β is a p
×
1 vector of unobservable parameters, and σ
2 is the unknown common variance. The matrix X is the observed value of Z with the measurement error U. Furthermore, ε and U are independent, Λ is a p
×
p matrix of known values with nonnegative diagonal elements (Fuller [14]) and I
n is an n
×
n identity matrix. The model (1) is known as functional linear measurement error model.The log-likelihood ofy is given by
(2)
l
β
,
σ
2
,
Z
,
y
=
-
n
2
log
⁡
⁡
2
π
σ
2
-
1
2
σ
2
y
-
Z
β
′
y
-
Z
β
.
If we replace Z by X without considering the measurement errors, then the ML estimates are not consistent in general. To correct for the effects of measurement errors on parameter estimation, we use the corrected score method proposed by Nakamura [16]. This method proposes a corrected log-likelihood l
*
(
β
,
σ
2
,
X
,
y
), which satisfies
(3)
E
*
∂
∂
β
l
*
β
,
σ
2
,
X
,
y
=
∂
∂
β
l
β
,
σ
2
,
Z
,
y
,
E
*
∂
∂
σ
2
l
*
β
,
σ
2
,
X
,
y
=
∂
∂
σ
2
l
β
,
σ
2
,
Z
,
y
,
where E
* denotes the conditional mean with respect to X given y. For model (1), the appropriate corrected log-likelihood is suggested by Nakamura [16] as
(4)
l
*
β
,
σ
2
,
X
,
y
=
-
n
2
log
⁡
⁡
2
π
σ
2
-
1
2
σ
2
y
-
X
β
′
y
-
X
β
-
n
β
′
Λ
β
.
By solving the equations (
∂
/
∂
β
)
l
*
(
β
,
σ
2
,
X
,
y
)
=
0 and (
∂
/
∂
σ
2
)
l
*
(
β
,
σ
2
,
X
,
y
)
=
0, the corrected score estimates of β and σ
2, respectively, are given by (see [16])
(5)
β
^
=
X
′
X
-
n
Λ
-
1
X
′
y
,
σ
^
2
=
1
n
y
-
X
β
^
′
y
-
X
β
^
-
n
β
^
′
Λ
β
^
.


## 3. A Variance Shift Model in the Linear Measurement Error

Suppose theith observation is considered with inflated error variance. A variance shift model for this observation takes the form
(6)
y
=
Z
β
+
b
i
d
i
+
ε
,
X
=
Z
+
U
,
where d
i is an n
×
1 vector with value 1 in the ith element and zero elsewhere and b
i is an unknown random coefficient of the form N
(
0
,
α
i
). Model (6) can be considered as a linear mixed measurement error model in which b
i is a random effect with the variance α
i and the covariance matrix of the data are as follows:
(7)
Var
⁡
y
=
α
i
d
i
d
i
′
+
σ
2
I
n
=
σ
2
w
i
-
1
d
i
d
i
′
+
I
n
=
σ
2
H
i
,
where w
i
=
α
i
/
σ
2
+
1  and  Var
⁡
y
i
=
w
i
σ
2, for w
i
≥
1.The log-likelihood and the corrected log-likelihood for this model are given by(8)
l
i
β
,
σ
2
,
w
i
,
Z
,
y
=
-
n
2
log
⁡
⁡
2
π
σ
2
-
1
2
log
⁡
H
i
-
1
2
σ
2
y
-
Z
β
′
H
i
-
1
y
-
Z
β
,
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
-
n
2
log
⁡
⁡
2
π
σ
2
-
1
2
log
⁡
H
i
-
1
2
σ
2
y
-
X
β
′
H
i
-
1
y
-
X
β
-
tr
⁡
H
i
-
1
β
′
Λ
β
,
respectively, which have the following properties:
(9)
E
∂
∂
β
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
∂
∂
β
l
i
β
,
σ
2
,
w
i
,
Z
,
y
,
E
∂
∂
σ
2
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
∂
∂
σ
2
l
i
β
,
σ
2
,
w
i
,
Z
,
y
,
E
∂
∂
w
i
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
∂
∂
w
i
l
i
β
,
σ
2
,
w
i
,
Z
,
y
.
Now, for w
i fixed, the corrected score estimates of β and σ
2 will be obtained with differentiating from the corrected log-likelihood of the variance shift model given in (8) with respect to β and σ
2 [24]. Then we have
(10)
β
^
i
w
i
=
X
′
H
i
-
1
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
H
i
-
1
y
,
(11)
σ
^
i
2
w
i
=
1
n
y
-
X
β
^
i
′
H
i
-
1
y
-
X
β
^
i
-
tr
⁡
H
i
-
1
β
^
i
′
Λ
β
^
i
,
respectively. In the following theorem we derive the asymptotic expressions of the β
^
i
w
i and σ
^
i
2
w
i as functions of β
^ and σ
^
2 given no outliers.Theorem 1.
For model (6), we have
(12)
β
^
i
w
i
=
β
^
-
w
i
-
1
1
+
r
i
i
w
i
-
1
X
′
X
-
n
Λ
-
1
x
i
v
^
i
+
O
p
n
-
1
,
σ
^
i
2
w
i
=
σ
^
v
2
1
-
r
i
i
w
i
-
1
1
+
r
i
i
w
i
-
1
t
i
2
n
-
β
^
′
Λ
β
^
+
O
p
n
-
1
,
where r
i
i is the ith diagonal element of R
=
I
-
X
(
X
′
X
-
n
Λ
)
-
1
X
′, v
^
i
=
y
i
-
x
i
′
β
^ is the ith residual and t
i
=
v
^
i
/
σ
^
ν
r
i
i is the ith studentized residual of the model, in which σ
^
υ
2
=
σ
^
2
+
β
^
′
Λ
β
^ [25].Proof.
From the estimate ofβ given in (10), we have
(13)
X
′
H
i
-
1
X
-
tr
⁡
H
i
-
1
Λ
-
1
=
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
w
i
-
1
w
i
X
′
d
i
d
i
′
X
-
1
=
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
+
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
d
i
·
w
i
-
1
w
i
-
1
-
d
i
′
X
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
d
i
-
1
·
d
i
′
X
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
.
On the other hand, we know that, X
′
X
=
O
p
n, X
′
y
=
O
p
n, d
i
′
X
X
′
d
i
=
O
p
1, X
′
d
i
d
i
′
X
=
O
p
1, (
X
′
X
-
n
Λ
)
-
1
=
O
p
(
n
-
1
), and  [
I
+
O
p
(
n
-
1
)
]
-
1
=
I
+
O
p
(
n
-
1
) (see the appendix for more details). Therefore, the first term in the right hand side of (13) will be
(14)
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
=
X
′
X
-
n
Λ
+
w
i
-
1
w
i
Λ
-
1
=
X
′
X
-
n
Λ
-
1
-
w
i
-
1
w
i
X
′
X
-
n
Λ
-
1
·
I
+
w
i
-
1
w
i
Λ
X
′
X
-
n
Λ
-
1
-
1
Λ
X
′
X
-
n
Λ
-
1
=
X
′
X
-
n
Λ
-
1
+
O
p
n
-
1
I
+
O
p
n
-
1
-
1
O
p
n
-
1
=
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
,
and the second term is
(15)
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
d
i
·
w
i
-
1
w
i
-
1
-
d
i
′
X
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
d
i
-
1
·
d
i
′
X
X
′
X
-
tr
⁡
H
i
-
1
Λ
-
1
=
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
X
′
d
i
·
w
i
-
1
w
i
-
1
h
h
h
h
h
h
-
d
i
′
X
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
X
′
d
i
w
i
-
1
w
i
-
1
-
1
·
d
i
′
X
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
=
w
i
-
1
1
+
r
i
i
w
i
-
1
X
′
X
-
n
Λ
-
1
X
′
d
i
d
i
′
X
X
′
X
-
n
Λ
-
1
+
O
p
n
-
3
.
Combining both terms of (13), we have
(16)
X
′
H
i
-
1
X
-
tr
⁡
H
i
-
1
Λ
-
1
=
X
′
X
-
n
Λ
-
1
+
w
i
-
1
1
+
r
i
i
w
i
-
1
X
′
X
-
n
Λ
-
1
X
′
d
i
d
i
′
X
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
.
Consequently, we have
(17)
β
^
i
w
i
=
X
′
H
i
-
1
X
-
tr
⁡
H
i
-
1
Λ
-
1
X
′
H
i
-
1
y
=
X
′
X
-
n
Λ
-
1
+
w
i
-
1
1
+
r
i
i
w
i
-
1
hhhh
·
X
′
X
-
n
Λ
-
1
X
′
d
i
d
i
′
X
X
′
X
-
n
Λ
-
1
+
O
p
n
-
2
w
i
-
1
1
+
r
i
i
w
i
-
1
·
X
′
I
-
w
i
-
1
w
i
d
i
d
i
′
y
.
Multiplying the above matrix expression out and simplifying them, β
^
i
w
i will be derived. Next, substituting β
^
i
w
i in the estimate of σ
2 given in (11), we can write
(18)
n
σ
^
i
2
w
i
=
y
-
X
β
^
i
′
H
i
-
1
y
-
X
β
^
i
-
tr
⁡
H
i
-
1
β
^
i
′
Λ
β
^
i
=
y
′
H
i
-
1
y
-
β
^
i
′
X
′
H
i
-
1
y
=
y
′
y
-
w
i
-
1
w
i
y
′
d
i
d
i
′
y
-
β
^
-
w
i
-
1
1
+
r
i
i
w
i
-
1
X
′
X
-
n
Λ
-
1
x
i
v
^
i
+
O
p
n
-
1
′
·
X
′
y
-
w
i
-
1
w
i
X
′
d
i
d
i
′
y
=
y
′
y
-
w
i
-
1
w
i
y
′
d
i
d
i
′
y
-
β
^
′
X
′
y
+
w
i
-
1
w
i
β
^
′
X
′
d
i
d
i
′
y
+
w
i
-
1
1
+
r
i
i
w
i
-
1
d
i
′
X
X
′
X
-
n
Λ
-
1
X
′
y
v
^
i
-
w
i
-
1
2
w
i
1
+
r
i
i
w
i
-
1
d
i
′
X
X
′
X
-
n
Λ
-
1
X
′
d
i
d
i
′
y
v
^
i
+
O
p
1
=
n
σ
^
2
-
w
i
-
1
1
+
r
i
i
w
i
-
1
v
^
i
2
+
O
p
1
=
n
σ
^
υ
2
1
-
r
i
i
w
i
-
1
1
+
r
i
i
w
i
-
1
t
i
2
n
-
n
β
^
′
Λ
β
^
+
O
p
1
,
and hence, σ
^
i
2
w
i will be obtained.In the rest of paper, we defineβ
~
i
w
i
=
β
^
-
(
(
w
i
-
1
)
/
(
1
+
r
i
i
(
w
i
-
1
)
)
)
(
X
′
X
-
n
Λ
)
-
1
x
i
v
^
i  and  σ
~
i
2
w
i
=
n
σ
^
v
2
(
1
-
(
r
i
i
(
w
i
-
1
)
/
(
1
+
r
i
i
(
w
i
-
1
)
)
)
(
t
i
2
/
n
)
)
-
n
β
^
′
Λ
β
^. It is obvious that for w
i
=
1 or r
i
i
=
0, β
~
i
w
i
=
β
^, σ
~
i
2
w
i
=
σ
^
2 and for r
i
i
>
0 we have
(19)
lim
⁡
w
i
→
+
∞
β
~
i
w
i
=
β
^
-
r
i
i
-
1
X
′
X
-
n
Λ
-
1
x
i
v
^
i
=
β
^
m
i
,
lim
⁡
w
i
→
+
∞
σ
~
i
2
w
i
=
σ
^
v
2
1
-
t
i
2
n
-
β
^
′
Λ
β
^
=
σ
^
m
i
2
.Remark 2.
β
^
m
i and σ
^
m
i
2 are the corrected score estimates of β and σ
2, respectively, in a mean shift outlier model for ith observation, given by [20, 25].
Nowl
i
*
(
β
,
σ
2
,
w
i
,
X
,
y
) evaluated at (
β
,
σ
2
)
=
(
β
^
i
(
w
i
)
,
σ
^
i
2
(
w
i
)
) is, except for an additive constant, proportional to
(20)
h
w
i
=
-
n
log
⁡
σ
^
i
2
w
i
-
log
⁡
w
i
.
Using Taylor series expansion, an approximate to the h
w
i analogues to one given by Cook et al. [4] will be as
(21)
h
~
w
i
=
-
n
log
⁡
σ
^
2
-
w
i
-
1
1
+
r
i
i
w
i
-
1
v
^
i
2
n
-
log
⁡
w
i
.
The existence of the value of w
i, say w
~
i, over the range of [
1
,
+
∞
) that maximizes h
~
w
i, was proved by Cook et al. [4].


## 4. Analogue of Likelihood Ratio Test and Score Test Statistics

In this section we derive analogues of likelihood ratio test and score test statistics for testing the null hypothesis that the observation is not unusual (w
i
=
1) against the alternative that it has an inflated variance, w
i
>
1.


### 4.1. Corrected Likelihood Ratio Test

Letl
i
*
β
^
i
w
~
i
,
σ
^
i
2
w
~
i
,
w
~
i
,
X
,
y and l
*
(
β
^
,
σ
^
2
,
X
,
y
) be the corrected log-likelihood evaluated at (
β
^
i
(
w
~
i
)
,
σ
^
i
2
(
w
~
i
)
,
w
~
i
) and under the null hypothesis, respectively, for testing. Consider
(22)
H
0
:
w
i
=
1
versus
H
1
:
w
i
>
1
.
The corrected likelihood ratio test is defined as
(23)
CLR
T
i
=
-
2
l
*
β
^
,
σ
^
2
,
X
,
y
-
l
i
*
β
^
i
w
~
i
,
σ
^
i
2
w
~
i
,
w
~
i
,
X
,
y
.


### 4.2. Score Test Statistic

In order to derive the score test statistic, we only require the estimate of parameters under the null hypothesis. In the following theorem based on the observed information matrix we obtain this test statistic under the hypothesis ofw
i
=
1.Theorem 3.
The score test statistic for theith observation (
S
C
i
), based on the observed information matrix for testing H
0
:
w
i
=
1, is given by
(24)
S
C
i
w
i
=
1
=
n
A
2
r
i
i
t
i
2
-
1
2
2
n
-
1
2
A
r
i
i
t
i
2
-
1
+
1
-
A
2
r
i
i
t
i
2
-
1
2
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhh
i
f
r
i
i
t
i
2
>
1
,
where A
=
σ
^
v
2
/
σ
^
2.Proof.
Let the corrected observed information matrix ofy for σ
2 and w
i be J
(
σ
2
,
w
i
), then the score test statistic (see [26]) for testing H
0
:
w
i
=
1 against H
1
:
w
i
>
1 is
(25)
SC
i
=
∂
∂
w
i
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
2
J
w
w
,
where J
w
w is the lower right corner of J
-
1
(
σ
2
,
w
i
). Substituting w
i
=
1, β
=
β
^, and σ
2
=
σ
^
2 into the elements of (25), we have
(26)
∂
∂
w
i
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
-
1
2
+
1
2
σ
^
2
v
^
i
2
-
β
^
′
Λ
β
^
=
1
2
A
r
i
i
t
i
2
-
1
,
J
σ
2
,
w
i
=
n
2
σ
^
4
1
2
σ
^
4
v
^
i
2
-
β
^
′
Λ
β
^
1
2
σ
^
4
v
^
i
2
-
β
^
′
Λ
β
^
-
1
2
+
1
σ
^
2
v
^
i
2
-
β
^
′
Λ
β
^
=
n
2
σ
^
4
1
2
σ
^
2
A
r
i
i
t
i
2
-
1
+
1
1
2
σ
^
2
A
r
i
i
t
i
2
-
1
+
1
A
r
i
i
t
i
2
-
1
+
1
2
,
where J
w
w
=
2
n
/
{
n
-
1
[
2
A
(
r
i
i
t
i
2
-
1
)
+
1
]
-
A
2
(
r
i
i
t
i
2
-
1
)
2
}, and then substituting in (25) the result is achieved.Because the null hypothesis is on the boundary of the parameter space, the standard asymptotic theory does not apply in this case [27]. Therefore, a parametric bootstrap procedure can be used to approximate the distributions of the likelihood ratio and the score test statistic (see Section 4.3 and [13]).


### 4.3. Empirical Distribution

Based on the Gumedze et al. [13] the following parametric bootstrap procedures, for test statistics CLR
T
i and SC
i, can be used to derive the empirical distributions of these statistics under the hypothesis of no outliers exist in the observations:Step  1. Fit model (1) to the data and calculate estimates Z
^ (see Zare and Rasekh [28]), β
^ and σ
^
2, where
(27)
Z
^
′
=
X
′
+
σ
^
υ
-
2
Λ
β
^
ν
^
′
.Step  2a. Generate a new data vector from
(28)
y
*
=
Z
^
β
^
+
ε
*
,
X
*
=
Z
^
+
U
,
where ε
* is randomly generated as N
(
0
,
σ
^
2
I
n
) and U is randomly generated as N
0
,
I
n
⊗
Λ. Fit model (1) to y
*.Step  2b. Compute the test statistic (CLRT
i or SC
i) for i
=
1
,
2
,
3
,
…
,
n, by fitting a variance shift model to simulated data y
* for each observation in turn and save the order statistics of the set {
CLRT
i or SC
i
:
i
=
1
,
2
,
3
,
…
,
n
}.Step  3. Repeat Steps  2a and 2b R times, for R acceptably large, for example, R
=
10000. Therefore, an empirical distribution of size R is generated for each order statistic.Step  4. Calculate the 100
(
1
-
α
) percentile for each order statistic for the level of size α.The percentile of thekth order statistic can be considered as a threshold for the kth largest value of the test statistic from the original data and if the k largest values of the test statistic from the original data all exceed their respective thresholds, then it is concluded that these are all outliers.


## 4.1. Corrected Likelihood Ratio Test

Letl
i
*
β
^
i
w
~
i
,
σ
^
i
2
w
~
i
,
w
~
i
,
X
,
y and l
*
(
β
^
,
σ
^
2
,
X
,
y
) be the corrected log-likelihood evaluated at (
β
^
i
(
w
~
i
)
,
σ
^
i
2
(
w
~
i
)
,
w
~
i
) and under the null hypothesis, respectively, for testing. Consider
(22)
H
0
:
w
i
=
1
versus
H
1
:
w
i
>
1
.
The corrected likelihood ratio test is defined as
(23)
CLR
T
i
=
-
2
l
*
β
^
,
σ
^
2
,
X
,
y
-
l
i
*
β
^
i
w
~
i
,
σ
^
i
2
w
~
i
,
w
~
i
,
X
,
y
.


## 4.2. Score Test Statistic

In order to derive the score test statistic, we only require the estimate of parameters under the null hypothesis. In the following theorem based on the observed information matrix we obtain this test statistic under the hypothesis ofw
i
=
1.Theorem 3.
The score test statistic for theith observation (
S
C
i
), based on the observed information matrix for testing H
0
:
w
i
=
1, is given by
(24)
S
C
i
w
i
=
1
=
n
A
2
r
i
i
t
i
2
-
1
2
2
n
-
1
2
A
r
i
i
t
i
2
-
1
+
1
-
A
2
r
i
i
t
i
2
-
1
2
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhh
i
f
r
i
i
t
i
2
>
1
,
where A
=
σ
^
v
2
/
σ
^
2.Proof.
Let the corrected observed information matrix ofy for σ
2 and w
i be J
(
σ
2
,
w
i
), then the score test statistic (see [26]) for testing H
0
:
w
i
=
1 against H
1
:
w
i
>
1 is
(25)
SC
i
=
∂
∂
w
i
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
2
J
w
w
,
where J
w
w is the lower right corner of J
-
1
(
σ
2
,
w
i
). Substituting w
i
=
1, β
=
β
^, and σ
2
=
σ
^
2 into the elements of (25), we have
(26)
∂
∂
w
i
l
i
*
β
,
σ
2
,
w
i
,
X
,
y
=
-
1
2
+
1
2
σ
^
2
v
^
i
2
-
β
^
′
Λ
β
^
=
1
2
A
r
i
i
t
i
2
-
1
,
J
σ
2
,
w
i
=
n
2
σ
^
4
1
2
σ
^
4
v
^
i
2
-
β
^
′
Λ
β
^
1
2
σ
^
4
v
^
i
2
-
β
^
′
Λ
β
^
-
1
2
+
1
σ
^
2
v
^
i
2
-
β
^
′
Λ
β
^
=
n
2
σ
^
4
1
2
σ
^
2
A
r
i
i
t
i
2
-
1
+
1
1
2
σ
^
2
A
r
i
i
t
i
2
-
1
+
1
A
r
i
i
t
i
2
-
1
+
1
2
,
where J
w
w
=
2
n
/
{
n
-
1
[
2
A
(
r
i
i
t
i
2
-
1
)
+
1
]
-
A
2
(
r
i
i
t
i
2
-
1
)
2
}, and then substituting in (25) the result is achieved.Because the null hypothesis is on the boundary of the parameter space, the standard asymptotic theory does not apply in this case [27]. Therefore, a parametric bootstrap procedure can be used to approximate the distributions of the likelihood ratio and the score test statistic (see Section 4.3 and [13]).


## 4.3. Empirical Distribution

Based on the Gumedze et al. [13] the following parametric bootstrap procedures, for test statistics CLR
T
i and SC
i, can be used to derive the empirical distributions of these statistics under the hypothesis of no outliers exist in the observations:Step  1. Fit model (1) to the data and calculate estimates Z
^ (see Zare and Rasekh [28]), β
^ and σ
^
2, where
(27)
Z
^
′
=
X
′
+
σ
^
υ
-
2
Λ
β
^
ν
^
′
.Step  2a. Generate a new data vector from
(28)
y
*
=
Z
^
β
^
+
ε
*
,
X
*
=
Z
^
+
U
,
where ε
* is randomly generated as N
(
0
,
σ
^
2
I
n
) and U is randomly generated as N
0
,
I
n
⊗
Λ. Fit model (1) to y
*.Step  2b. Compute the test statistic (CLRT
i or SC
i) for i
=
1
,
2
,
3
,
…
,
n, by fitting a variance shift model to simulated data y
* for each observation in turn and save the order statistics of the set {
CLRT
i or SC
i
:
i
=
1
,
2
,
3
,
…
,
n
}.Step  3. Repeat Steps  2a and 2b R times, for R acceptably large, for example, R
=
10000. Therefore, an empirical distribution of size R is generated for each order statistic.Step  4. Calculate the 100
(
1
-
α
) percentile for each order statistic for the level of size α.The percentile of thekth order statistic can be considered as a threshold for the kth largest value of the test statistic from the original data and if the k largest values of the test statistic from the original data all exceed their respective thresholds, then it is concluded that these are all outliers.


## 5. Simulation Study

A parametric bootstrap simulation study is carried out to demonstrate the empirical performance of the various proposed test statistics in terms of the probability of a typeI error and power on a single unit.The response variabley
j is generated from the model y
j
=
Z
β
+
ε
j; j
=
1,2
,
…
,
1100, where y
j
=
(
y
1
j
,
y
2
j
,
…
,
y
n
j
)
′, Z
=
1
n
,
z, 1
n is a n
×
1 vector all of whose elements are 1’s, z
′
=
z
1
,
z
2
,
…
,
z
n and ε
j is rewritten in accordance with y
j. We consider the following combinations for simulation: n
=
40 or 100 where β
=
(
β
0
,
β
1
)
=
10,0.3 or (
β
0
,
β
1
)
=
20,1, z
i
~
U
100,900, ε
i
j
~
N
(
0
,
σ
2
), i
=
1,2
,
…
,
n, σ
2
=
2 or 8 and Λ
=
diag
⁡
(
λ
1
,
λ
2
)
=
diag
⁡
0,0.25 or Λ
=
diag
⁡
0,1. The simulation study was conducted using the R software and the codes are available from the second author upon request.For each simulated data set, the CLRT and the score test statistics were calculated for the first observation. The choice of the first observation was arbitrary. To generate an empirical distribution of the test statistic under the null hypothesis, data sets fork
=
1
,
2
,
…
,
2500 were simulated as y
j
k
*
=
Z
^
j
β
^
j
+
ε
j
k
*, where ε
j
k
*
~
N
(
0
,
σ
^
j
2
I
n
), β
^
j, Z
^
j, and σ
^
j
2 are the corrected estimate of β, Z and σ
2 from y
j. The probability of a type I error estimate for a given test statistic and α
=
0.05 was calculated as the number of data sets for which the test statistic exceeded the 95th percentile of the empirical distribution, divided by 1100 [13].The CLRT and the score test statistics were performed for a variance shift model for the first observation of each simulated data and 95th percentiles from the empirical distribution of each test statistics were used as threshold values for the test statistics observed on the original data sety
j. The empirical probability of type I errors for thresholds derived from the empirical distribution under the null hypothesis are calculated for the corrected likelihood ratio and score test statistics for α
=
0.05(Tables 1, 2, and 3). A glance at the results of these tables indicate that in general the probability of a type I errors of both CLRT and score test statistics are close to the nominal value of 0.05.Table 1
The probability of a typeΙ error (α
=
0.05) and power of CLRT and score test statistics for a variance shift model in a linear measurement error with combination different of parameters β, σ
2, Λ, and b
=
1.
n
β
σ
2
(λ
1, λ
2)
CLRT
Score test
Sig level
Power
Sig level
Power
40
(20, 1)
2
(0, 1)
0.0457
0.0620
0.0459
0.0621
(0, 0.25)
0.0478
0.0870
0.0477
0.0873
8
(0, 1)
0.0478
0.0511
0.0477
0.0513
(0, 0.25)
0.0411
0.0600
0.0414
0.0603
(10, 0.3)
2
(0, 1)
0.0443
0.0909
0.0441
0.0909
(0, 0.25)
0.0388
0.0944
0.0387
0.0945
8
(0, 1)
0.0388
0.0600
0.0387
0.0603
(0, 0.25)
0.0354
0.0619
0.0360
0.0621
100
(20, 1)
2
(0, 1)
0.0450
0.0801
0.0450
0.0801
(0, 0.25)
0.0520
0.0950
0.0522
0.0954
8
(0, 1)
0.0520
0.0649
0.0522
0.0648
(0, 0.25)
0.0529
0.0668
0.0531
0.0666
(10, 0.3)
2
(0, 1)
0.0520
0.0950
0.0522
0.0954
(0, 0.25)
0.0520
0.0989
0.0522
0.0990
8
(0, 1)
0.0529
0.0631
0.0531
0.0630
(0, 0.25)
0.0520
0.0636
0.0522
0.0639Table 2
The probability of a typeΙ error (α
=
0.05) and power of CLRT and score test statistics for a variance shift model in a linear measurement error with combination different of parameters β, σ
2, Λ, and b
=
3.
n
β
σ
2
(λ
1, λ
2)
CLRT
Score test
Sig level
Power
Sig level
Power
40
(20, 1)
2
(0, 1)
0.0455
0.2840
0.0453
0.2844
(0, 0.25)
0.0476
0.4311
0.0478
0.4311
8
(0, 1)
0.0478
0.1385
0.0475
0.1386
(0, 0.25)
0.0419
0.1558
0.0419
0.1557
(10, 0.3)
2
(0, 1)
0.0413
0.4768
0.0417
0.4770
(0, 0.25)
0.0364
0.4832
0.0364
0.4833
8
(0, 1)
0.0367
0.1592
0.0363
0.1593
(0, 0.25)
0.0368
0.1608
0.0366
0.1611
100
(20, 1)
2
(0, 1)
0.0450
0.3660
0.0453
0.3663
(0, 0.25)
0.0525
0.4823
0.0529
0.4824
8
(0, 1)
0.0520
0.1584
0.0521
0.1584
(0, 0.25)
0.0537
0.1629
0.0534
0.1629
(10, 0.3)
2
(0, 1)
0.0524
0.5239
0.0522
0.5239
(0, 0.25)
0.0531
0.5394
0.0528
0.5391
8
(0, 1)
0.0529
0.1653
0.0534
0.1656
(0, 0.25)
0.0525
0.1672
0.0525
0.1674Table 3
The probability of a typeΙ error (α
=
0.05) and power of CLRT and score test statistics for a variance shift model in a linear measurement error with combination different of parameters β, σ
2, Λ, and b
=
5.
n
β
σ
2
(λ
1, λ
2)
CLRT
Score test
Sig level
Power
Sig level
Power
40
(20, 1)
2
(0, 1)
0.0457
0.7172
0.0459
0.7174
(0, 0.25)
0.0478
0.8755
0.0475
0.8758
8
(0, 1)
0.0476
0.3169
0.0474
0.3168
(0, 0.25)
0.0411
0.3536
0.0414
0.3537
(10, 0.3)
2
(0, 1)
0.0449
0.9002
0.0446
0.9001
(0, 0.25)
0.0389
0.9082
0.0387
0.9082
8
(0, 1)
0.0388
0.3600
0.0387
0.3600
(0, 0.25)
0.0366
0.3650
0.0364
0.3654
100
(20, 1)
2
(0, 1)
0.0457
0.7893
0.0453
0.7894
(0, 0.25)
0.0523
0.9035
0.0521
0.9037
8
(0, 1)
0.0525
0.3547
0.0528
0.3546
(0, 0.25)
0.0535
0.3897
0.0532
0.3897
(10, 0.3)
2
(0, 1)
0.0525
0.9229
0.0522
0.9226
(0, 0.25)
0.0521
0.9272
0.0521
0.9271
8
(0, 1)
0.0536
0.4024
0.0535
0.4023
(0, 0.25)
0.0520
0.4059
0.0522
0.4059In order to access the relative sensitivity of the CLRT and score test statistics, we introduce the shift values 1, 3, and 5 for the first observation and again for each combination of parameters, 1100 data sets are generated from the following model:(29)
y
j
=
Z
β
+
b
d
1
+
ε
j
,
for b
=
1, 3, or 5, where d
1 is an n
×
1 vector with value 1 in the first element and zero elsewhere. The CLRT, the score test statistic, and their empirical distribution were calculated as for the probability of a type I error given above, consequently the power of the test statistics are also derived. Results of Tables 1–3 show that with increase of the displacement, b, the power of the CLRT and score test statistic, increases in general. Moreover, we can see that power of the test statistics also increase as sample size increases. These tables also show the result of the CLRT and score test statistics are nearly identical in the empirical probability of a type I errors and power.


## 6. Example: Concrete Compressive Strength Data

These data were given by Wellman and Gunst [19] and contain comprehensive strength measurements of 41 sample of concrete. It was desired to use a linear regression model to predict comprehensive strength of concrete 28 days after pouring from the strength measurements taken two days after pouring. Zhong et al. [20] analyzed this data set using the linear measurement error model with Λ
=
diag
⁡
0
,
1. The zero diagonal element in Λ corresponds to the constant term, a predictor variable measured without error. They indicated that the sample 21 exhibits a strong influence on the fitted model. Here we consider a variance shift measurement error model for this data set.Figures1, 2, and 3 show plots of the square of studentized residual t
i
2 of the data under model (1), estimates of the variance shift parameter w
~
i, and estimated variance σ
^
i
2 under model (6), versus case numbers, respectively. From these figures it is obvious that case 21 stands out as a possible outlier with relatively large values of t
i
2 and w
~
i and a small estimated σ
^
i
2.Figure 1
Plot of square of the studentized residual(
t
i
2
) versus case number.Figure 2
Plot of variance shift estimates(
w
~
i
) versus case number.Figure 3
Plot of error variance estimates (σ
^
i
2) versus case number.Next, the corrected likelihood ratio and score test statistics were calculated for each observation under model (6), and then 10000 simulated data sets were generated from the fitted model under the null hypothesis (
H
0
:
w
i
=
1
). In each simulation, a variance shift model was fitted for each observation and the test statistics were sorted and used to generate the empirical distribution of the order statistics for each test [13]. Figures 4 and 5 give plots of the test statistics from the real data and 95th percentile from the empirical distribution of the first, second and third largest values for each test statistic. These figures show that the statistics for observation 21 is larger than the 95th percentile of the distribution of the corresponding order statistics.Figure 4
Corrected likelihood ratio test statistic for each case, with 95th ercentile of the empirical distribution underH
0 shown for the first k order statistics for corrected likelihood ratio test: k
=
1 (solid line), k
=
2 (dashed line), and k
=
3 (dotted line).Figure 5
Score test statistic for each case, with 95th percentile of the empirical distribution underH
0 shown for the first k order statistics for score test: k
=
1 (solid line), k
=
2 (dashed line), and k
=
3 (dotted line).Finally, the fitted regression line measurement error model and the variance shift measurement error model for case 21 are shown in Figure6.Figure 6
Concrete compressive strengths, in pounds per square inch, at 2 and 28 days, measurement error model (dashed line); a variance shift model (solid line).


## 7. Conclusions

We extended the variance shift model to the linear measurement error models based on the corrected likelihood of Nakamura [16]. We derived the approximate estimate of parameters of the proposed model under the variance shift and indicated that if the variance shift parameter tends to infinity, these estimates will be the same as those obtained from a mean shift outlier model. Also, we proposed a corrected likelihood ratio test and derived the score test statistic for testing that an observation stands out as possible outlier and it is shown that the score test statistic is function of studentized residuals of model. The performance of both the corrected likelihood ratio and the score test statistics is studied using a parametric bootstrap simulation, and it was found out that with the increase of b or sample size, the power of both test statistics increases.



---

*Source: 396875-2014-09-14.xml*