# Local Influence Analysis for Quasi-Likelihood Nonlinear Models with Random Effects

Authors: Tian Xia; Jiancheng Jiang; Xuejun Jiang
Journal: Journal of Probability and Statistics
 (2018)

---

## Abstract

We propose a quasi-likelihood nonlinear model with random effects, which is a hybrid extension of quasi-likelihood nonlinear models and generalized linear mixed models. It includes a wide class of existing models as examples. A novel penalized quasi-likelihood estimation method is introduced. Based on the Laplace approximation and a penalized quasi-likelihood displacement, local influence of minor perturbations on the data set is investigated for the proposed model. Four concrete perturbation schemes are considered in the local influence analysis. The effectiveness of the proposed methodology is illustrated by some numerical examinations on a pharmacokinetics dataset.


---

## Body

## 1. Introduction

In this paper, we propose a quasi-likelihood nonlinear model with random effects (QLNMWRE) and investigate local influence of the model. The QLNMWRE is a hybrid generalization of quasi-likelihood nonlinear models [1, 2] and generalized linear mixed models, and it combines the advantages of both models. Generalized linear mixed models (GLMMs) are extensions of the well-known generalized linear models [3] by adding random effects to the linear predictor. GLMMs are effective and flexible for modeling nonnormal responses, repeated measurements, and other forms of clustered data. Efficient inference for the GLMMs depends on the underlying distribution of the data. Nevertheless, the exact distribution is rarely known in practice. In contrast, the quasi-likelihood method [4] requires only the first and second moments assumptions about the distribution and has been widely applied in the theory and practice of statistics (see, e.g., [5–8]).Detecting influential observations is important in data analysis. The local influence analysis has become a general tool for detecting a group of points with great influence on the fitted model through perturbation schemes [9]. This approach has been successfully applied in many models, such as mixed models [10, 11], generalized linear models [12], generalized linear mixed models [13], exponential family nonlinear models [14], nonlinear reproductive dispersion mixed model [15], nonlinear mixed-effect models [16, 17], and multivariate threshold time series models [1]. However, in these references the local influence method severely depends on the likelihood displacement, which is rarely known in practice. Instead, quasi-likelihood methods do not require the exact likelihood function except the first two moments of the response variables. Hence, we conduct influence analysis of the QLNMWRE using a novel penalized quasi-likelihood estimation method. The proposed methodology is illustrated by analyzing the pharmacokinetics dataset.The remainder of this paper is organized as follows. In Section2, we introduce the QLNMWRE and the corresponding estimation method. A Fisher-scoring iteration algorithm is advanced to calculate the estimators. In Section 3, a penalized quasi-likelihood displacement (PQLD) is proposed, and assessment of local influence under four different perturbation schemes is investigated. In Section 4, the pharmacokinetics dataset is employed to illustrate the effectiveness of the proposed methodology. Finally, we make discussion in Section 5.


## 2. Models and Estimation Method

LetY be a response vector of length n, and let X and Z be n
×
k and n
×
q matrices of explanatory variables associated with fixed and random effects, respectively. Conditional on the q
-dimensional vector of random effects, b, the observations, {
y
i
,
i
=
1
,
…
,
n
}, are independent and satisfy that(1)
E
y
i
∣
b
=
h
x
i
,
β
+
z
i
T
b
≜
μ
i
v
a
r
y
i
∣
b
=
σ
2
v
μ
i
,where β
=
(
β
1
,
β
2
,
…
,
β
p
)
T
(
p
<
n
) is an unknown parameter vector defined in a compact subset B
⊂
R
p, x
i and z
i are defined in a subset X of R
k and a subset Z of R
q, respectively, v
(
·
) is a known variance function, σ
2 is a dispersion parameter that is known or can be estimated separately, h
(
·
,
·
) is a continuously differentiable function such that the derivative matrix D
=
D
β
(
β
)
=
∂
h
(
β
)
/
∂
β
T has rank p for all β, with h
(
β
)
=
(
h
(
x
1
,
β
)
,
…
,
h
(
x
n
,
β
)
)
T, and the random effects b are assumed to be multivariate normally distributed:(2)
b
~
N
0
,
σ
2
Σ
,with Σ being a known nonnegative definite matrix. Following [2, 3, 18, 19], the conditional log quasi-likelihood on b is defined as(3)
Q
1
β
;
y
=
∑
i
=
1
n
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
,where μ
i
=
h
(
x
i
,
β
)
+
z
i
T
b
≜
μ
i
(
β
)
. The model defined by (1)-(3) is the so-called QLNMWRE.Clearly, this QLNMWRE encompasses some important special cases. IfΣ
=
0, then the above model is just the quasi-likelihood nonlinear model discussed by [2]; if μ
i
(
β
)
=
h
(
x
i
T
β
)
+
z
i
T
b, and y
i are independently drawn from a one-parameter exponential family of distributions with density (4)
exp
⁡
θ
i
T
y
i
-
k
θ
i
d
γ
y
i
,
i
=
1
,
…
,
n
,where γ
(
·
) is a measure, then it reduces to generalized linear models with random effects (see [20, 21]). Hence, the QLNMWRE is a hybrid extension of the quasi-likelihood nonlinear models and the generalized linear models with random effects.Letp
(
b
|
σ
2
) be a probability density function of random effect b. Then the joint log quasi-likelihood function of y
=
(
y
1
,
…
,
y
n
)
T and b is(5)
Q
β
,
σ
2
;
y
,
b
=
Q
1
β
;
y
+
log
⁡
p
b
∣
σ
2
=
∑
i
=
1
n
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
-
q
2
log
⁡
2
π
σ
2
-
1
2
log
⁡
Σ
-
1
2
σ
2
b
T
Σ
-
1
b
.Similar to the relationship between the joint log-likelihood function and the marginal log -likelihood function, we have (6)
Q
β
,
σ
2
;
y
,
b
=
Q
β
,
σ
2
;
y
+
Q
β
,
σ
2
;
b
∣
y
,where Q
(
β
,
σ
2
;
y
) is the marginal log quasi-likelihood function of y and Q
(
β
,
σ
2
;
b
∣
y
) is the log quasi-likelihood function of b given y, i.e., (7)
Q
β
,
σ
2
;
y
=
log
⁡
∫
exp
⁡
Q
β
,
σ
2
;
y
,
b
d
b
,
Q
β
,
σ
2
;
b
∣
y
=
log
⁡
exp
⁡
Q
β
,
σ
2
;
y
,
b
∫
exp
⁡
Q
β
,
σ
2
;
y
,
b
d
b
.Following the arguments in [20], the integrated log quasi-likelihood function used to estimate β is defined by(8)
e
x
p
Q
β
,
σ
2
;
y
∝
Σ
-
1
/
2
∫
exp
⁡
-
1
2
σ
2
∑
i
=
1
n
d
i
y
i
,
μ
i
-
1
2
σ
2
b
T
Σ
-
1
b
d
b
,where d
i
(
y
;
μ
)
=
-
2
∫
μ
y
(
y
-
t
/
v
t
)
d
t denotes the deviance measure of fit. If, conditional on b, y
i is a member of the exponential family, then -
d
i
(
y
i
;
μ
i
)
/
(
2
σ
2
) is the conditional log-likelihood of y
i given b, and E
b
[
∑
i
=
1
n
d
i
(
y
i
;
μ
i
)
/
(
2
σ
2
)
] is the log-likelihood function.In general, no analytical expressions are available for the integral in (8) and approximate techniques are needed. The simplest approach is the Laplace approximation [22, 23]. Obviously, the right-hand side of (8) is (9)
c
Σ
-
1
/
2
∫
e
-
m
b
d
b
,where m
(
b
)
=
(
1
/
2
σ
2
)
∑
i
=
1
n
d
i
y
i
,
μ
i
+
(
1
/
2
σ
2
)
b
T
Σ
-
1
b. When the Laplace method is applied to approximate the integrated quasi-likelihood function (8), estimates of β for fixed σ
2 are obtained by maximizing the penalized quasi-likelihood (PQL) (8):(10)
Q
p
β
=
-
1
2
σ
2
∑
i
=
1
n
d
i
y
i
,
μ
~
i
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
=
∑
i
=
1
n
∫
y
i
μ
~
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
~
Σ
-
1
b
~
=
Q
1
β
;
Y
μ
i
=
μ
~
i
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
i
=
f
(
x
i
,
β
)
+
z
i
T
b
~, and b
~
≜
b
~
(
β
) is the root of ∂
m
(
b
)
/
∂
b
=
0 for fixed β. We will use the penalized quasi-likelihood Q
p
(
β
) to estimate β and to conduct local influence analysis. To this end, we need the following assumptions.Assumption A.
(i)E
(
v
(
μ
i
)
)
-
1
(
y
i
-
μ
i
)
μ
i
=
μ
~
i
=
0
,
∀
i
=
1
,
…
,
n;
(ii) there exists some constantM
>
0 and some compact subset B
1
⊂
B such that (11)
sup
i
≥
1
,
β
⊆
B
1
⁡
E
v
μ
i
-
1
y
i
-
μ
i
μ
i
=
μ
~
i
2
≤
M
.It is easily seen that Assumption A holds in generalized linear mixed models and exponential family nonlinear random effects models. Assumption A guarantees existence of the variance-covariance matrix ofe
~, where e
~
=
(
e
~
1
,
⋯
,
e
~
n
)
T with e
~
i
=
v
(
μ
i
)
-
1
(
y
i
-
μ
i
)
μ
i
=
μ
~
i. Let Z
=
(
z
1
,
…
,
z
n
)
T and K
=
∂
e
~
/
∂
μ
T
=
d
i
a
g
(
k
1
,
…
,
k
n
), where (12)
k
i
=
∂
∂
μ
i
y
i
-
μ
i
v
μ
i
μ
i
=
μ
~
i
.Put Ω
1
=
K
-
1
-
Z
Σ
Z
T, W
=
∂
2
h
(
β
)
/
∂
β
∂
β
T, Ω
=
K
0
-
1
-
Z
Σ
Z
T, and K
0
=
E
y
(
K
). Under Assumption A, we have the following result.Theorem 1.
For the model defined by (1)-(3), conditional on b
~, the quasi-score function, the quasi-observed information matrix, and the quasi-Fisher information matrix for β admit the following representations:(13)
S
n
β
≜
Q
˙
p
β
=
σ
-
2
∂
h
β
∂
β
T
T
e
~
=
σ
-
2
D
T
e
~
,
(14)
H
n
β
≜
-
Q
¨
p
β
=
-
∂
2
Q
p
β
∂
β
∂
β
T
=
-
σ
-
2
e
~
T
W
+
D
T
Ω
1
D
,
(15)
F
n
β
≜
E
y
-
Q
¨
p
β
=
-
σ
-
2
D
T
Ω
-
1
D
,where ·
[
·
] indicates the array multiplication.Letβ
^
n denote the maximum quasi-likelihood estimator (MQLE) of β, which is the solution of equation Q
˙
p
(
β
)
=
0. Then the Fisher-scoring iteration method can be used for computing β
^
n by iteratively solving the following equation (see [14, 24]):(16)
β
i
+
1
=
β
i
+
F
n
-
1
β
i
S
n
β
i
=
β
i
-
D
T
Ω
-
1
D
-
1
D
T
e
~
i
=
0,1
,
2
,
⋯where D
,
Ω, and e
~ are all evaluated at β
i and b
i.On the other hand, it follows from (5) that the quasi-score function and the quasi-Fisher information matrix for b can be, respectively, expressed as (17)
S
b
≜
∂
Q
β
,
σ
2
;
y
,
b
∂
b
=
σ
-
2
Z
T
e
-
Σ
-
1
b
,
F
b
≜
E
y
-
∂
2
Q
β
,
σ
2
;
y
,
b
∂
b
∂
b
T
=
σ
-
2
Z
T
V
-
1
Z
+
Σ
-
1
,where e
=
e
(
β
)
=
(
e
1
,
…
,
e
n
)
T with e
i
=
(
y
i
-
μ
i
)
/
v
(
μ
i
), and (18)
V
-
1
=
d
i
a
g
v
-
1
μ
1
,
…
,
v
-
1
μ
n
.Hence, the Fisher-scoring iteration algorithm for computing the predictor of b
j under known β
j is given by(19)
b
j
i
+
1
=
b
j
i
+
Z
T
V
-
1
Z
+
Σ
-
1
-
1
Z
T
e
-
Σ
-
1
b
j
i
,
j
i
=
0,1
,
2
,
…
,where V and e are all evaluated at b
(
j
i
) and β
j. As the iteration scheme (19) converges, b
j
i converges to b
j.In general, the choice of initial valueβ
0 is important for the Fisher-scoring iteration algorithm. We use the algorithm in [2] for quasi-likelihood nonlinear models to find the starting values of parameter β for QLNMWRE with b
0
=
⋯
=
b
q
=
0. Hence, the MQLE β
^ of β can be obtained by solving (16) and (19) until convergence.In order to investigate the statistical diagnostic measures for QLNMWRE, we rewrite (16) (20)
β
i
+
1
=
D
T
Ω
-
1
D
-
1
D
T
Ω
-
1
G
β
i
,where G
=
D
β
i
-
Ω
e
~. When β
i converges to β
^, β
^ can be expressed as(21)
β
^
=
D
T
Ω
-
1
D
-
1
D
-
1
Ω
-
1
G
β
^
,where G
=
D
β
^
-
Ω
e
^, D
,
Ω and e
^ are all evaluated at β
^.


## 3. Local Influence

The aim of local influence analysis is to investigate the behavior of some influence measureT
(
ω
) when small perturbations are made in the model/data, where ω is an m-dimensional vector of perturbations restricted to some open subset Θ
∈
R
m. For simple statistical models, Cook constructed in [9] the likelihood displacement L
D
(
ω
) and used it to assess the local influence of a minor perturbation. Although this approach is very useful, serious difficulties are encountered when applying it to complicated models, because of the intractable likelihood function. For the sake of coping with those difficulties, some authors have considered alternatives to replace LD(
ω
). For instance, Zhu et al. proposed in [25] the Q-likelihood displacement and established an approach to assess local influence of statistical models with incomplete data, and Jung presented in [26] a quasi-likelihood displacement to obtain local influence analysis in generalized estimating equations. Inspired by [25, 26], we define in this work a new penalized quasi-likelihood displacement and then adapt the local influence approach introduced by [9] to the QLNMWRE.LetQ
p
(
β
) and Q
p
(
β
|
ω
) be the penalized quasi-likelihood for the unperturbed and perturbed models, respectively. We assume that there is an ω
0 such as Q
p
(
β
|
ω
0
)
=
Q
p
(
β
). Let β
^ and β
^
(
ω
) be the MQLE of β under the unperturbed and perturbed models, respectively. Similar to the likelihood displacement [9], we define the penalized quasi-likelihood displacement (PQLD) as (22)
P
Q
L
D
ω
=
2
Q
p
β
^
-
Q
p
β
^
ω
.The influence graph is defined as α
(
ω
)
=
(
ω
T
,
P
Q
L
D
(
ω
)
)
T. Following the approach developed in [9, 25, 26], the normal curvature C
l of α
(
ω
) at ω
0 in the direction of some unit vector l can be used to summarize the local behavior of the penalized quasi-likelihood displacement. As shown in [9], the normal curvature C
l in the unit direction l
(
l
=
1
) at ω
0 is given by(23)
C
l
=
2
l
T
F
¨
l
,where F
¨
=
-
(
∂
2
Q
p
(
β
∣
ω
)
/
∂
ω
∂
ω
T
)
ω
=
ω
0
=
-
△
T
Q
¨
p
-
1
△, in which △
=
∂
2
Q
p
(
β
|
ω
)
/
∂
β
∂
ω
T is a p
×
m matrix evaluated at β
=
β
^ and ω
=
ω
0, Q
¨
p
=
∂
2
Q
p
(
β
)
/
∂
β
∂
β
T is a p
×
p Hessian matrix evaluated at β
=
β
^. The maximum curvature C
m
a
x, which is the largest absolute eigenvalue of 2
F
¨, and the corresponding direction vector l
m
a
x are usually used for identifying locally influential observations. A large value of C
m
a
x is an indication of a serious local problem, and if the i-th element in l
m
a
x is relatively large special attention should be paid to the element being perturbed by ω
i. To apply the local influence method in [9] to the QLNMWRE, we consider the following four perturbation schemes.


### 3.1. Case-Weights Perturbation

Letω be an n
×
1 perturbation vector such that ω
0
=
(
1
,
⋯
,
1
)
T. The joint log quasi-likelihood function for the perturbed model is given by (24)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
i
=
1
n
ω
i
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C
=
-
q
/
2
log
⁡
2
π
σ
2
-
1
/
2
Σ. Then the penalized quasi-likelihood function can be expressed as(25)
Q
p
β
∣
ω
=
∑
i
=
1
n
ω
i
∫
y
i
μ
~
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
i
=
h
(
x
i
,
β
)
+
z
i
T
b
~, Z
=
(
z
1
,
z
2
,
…
,
z
n
)
T, and b
~ satisfies (26)
∑
i
=
1
n
w
i
y
i
-
μ
~
i
v
μ
~
i
z
i
-
Σ
-
1
b
~
=
0
.Hence, b
~
=
Σ
Z
T
W
¯
e
~, where W
¯
=
d
i
a
g
(
ω
1
,
…
,
ω
n
). Then (27)
∂
2
Q
p
β
∣
w
∂
β
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
δ
i
e
^
i
+
σ
-
2
D
T
K
·
K
-
1
Ω
-
1
K
-
1
-
I
δ
i
e
^
i
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
δ
i
e
^
i
,and thus (28)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
E
,where E
=
d
i
a
g
(
e
^
1
,
…
,
e
^
n
).


### 3.2. Response Variable Perturbation

A perturbation of the response variables(
y
1
,
⋯
,
y
n
)
T is introduced by replacing y
i by y
i
ω
=
y
i
+
ω
i, where ω
=
(
ω
1
,
⋯
,
ω
n
)
T, and ω
0
=
(
0
,
…
,
0
)
T represents the situation with no perturbation. In this case, the joint log quasi-likelihood function for the perturbed model is given by (29)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
i
=
1
n
∫
y
i
+
ω
i
μ
i
y
i
+
ω
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C is a constant. It follows from Section 2 that the penalized quasi-likelihood function is (30)
Q
p
β
∣
ω
=
∑
i
=
1
n
∫
y
i
+
ω
i
μ
~
i
y
i
+
ω
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
i
=
h
(
x
i
,
β
)
+
z
i
T
b
~, and b
~ satisfies (31)
∑
i
=
1
n
y
i
+
w
i
-
μ
~
i
v
μ
~
i
z
i
-
Σ
-
1
b
~
=
0
.It follows that b
~
=
Σ
Z
T
(
e
~
+
W
v
), where W
v
=
(
w
v
1
,
…
,
w
v
n
)
T with w
v
i
=
(
w
i
/
v
μ
i
)
μ
i
=
μ
~
i. Then (32)
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
∂
ω
i
+
∂
W
v
∂
ω
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
∂
μ
~
T
T
∂
μ
~
∂
ω
i
+
δ
i
T
1
v
μ
~
i
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
-
Ω
1
δ
i
T
1
v
μ
~
i
+
δ
i
T
1
v
μ
~
i
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
δ
i
T
1
v
μ
~
i
,and (33)
△
=
∂
2
Q
p
β
∣
w
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
E
∗
,where E
∗
=
d
i
a
g
(
e
1
∗
,
…
,
e
n
∗
) with e
i
∗
=
1
/
v
(
μ
~
i
) and μ
^
i
=
h
(
x
i
,
β
^
)
+
z
i
T
b
^.


### 3.3. Explanatory Variables Perturbation

In this case, we focus on the perturbation of a specific explanatory variable. Under this condition we have the perturbed explanatory matrixX
ω
=
(
x
1
,
⋯
,
x
t
ω
,
⋯
,
x
n
)
T with x
t
ω
=
x
t
+
ω, where x
t is a single explanatory variable of matrix X
ω corresponding to y
t and ω
0
=
(
0
,
…
,
0
)
T denotes no perturbation. Then the joint log quasi-likelihood function for the perturbed model is (34)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
j
=
1
n
∫
y
i
μ
j
w
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
Σ
-
1
b
,where C is a constant, μ
j
ω
=
h
(
x
j
,
β
)
+
z
j
T
b
(
j
≠
i
), and μ
t
ω
=
h
(
x
t
ω
β
)
+
z
t
T
b. It follows from Section 2 that (35)
Q
p
β
∣
ω
=
∑
j
=
1
n
∫
y
i
μ
~
j
w
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
Σ
-
1
b
~
,where μ
~
j
ω
=
h
(
x
j
,
β
)
+
z
j
T
b
~
(
j
≠
t
), μ
~
t
ω
=
h
(
x
t
ω
,
β
)
+
z
t
T
b
~, and b
~ satisfies (36)
∑
j
=
1
n
y
i
-
μ
~
j
w
v
μ
~
j
w
z
i
-
Σ
-
1
b
~
=
0
.Therefore, b
~
=
Σ
Z
T
e
~
w and e
~
w
=
(
(
y
1
-
μ
~
1
)
/
v
(
μ
~
1
)
,
⋯
,
(
y
n
-
μ
~
n
)
/
v
(
μ
~
n
)
)
T. Let h
t
T
=
∂
h
(
x
t
,
β
)
/
∂
x
t
T
β
=
β
^ and H
b
t
=
∂
2
h
(
β
)
/
∂
β
∂
x
t
T
β
=
β
^. Then (37)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
e
^
ω
T
∂
∂
ω
T
∂
h
β
∂
β
ω
0
,
β
^
+
σ
-
2
D
T
∂
e
~
w
∂
μ
~
∂
μ
~
∂
ω
T
ω
0
,
β
^
=
σ
-
2
e
^
T
H
b
t
+
σ
-
2
D
T
Ω
1
-
1
δ
t
h
t
T
,where ·
[
·
] indicates the array multiplication.


### 3.4. Perturbation of Covariates in Random Effects

Consider perturbing the data for thekth explanatory variable of Z, by modifying the data matrix Z as Z
w
=
Z
+
ω
d
k
T, where d
k is a q
-vector with 1 at kth position and zeros elsewhere. Under this situation, the perturbed joint log quasi-likelihood can be expressed as (38)
Q
β
,
σ
2
;
y
,
b
∣
ω
=
C
+
∑
i
=
1
n
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C is a quantity that does not depend on β and ω, and μ
=
h
(
β
)
+
Z
w
b. When ω
0
=
0, it indicates no perturbation. It follows from Section 2 that (39)
Q
p
β
∣
ω
=
∑
i
=
1
n
∫
y
i
μ
~
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
=
h
(
β
)
+
Z
w
b
~, and b
~ satisfies (40)
Z
w
T
e
~
w
-
Σ
-
1
b
~
=
0
,and therefore, b
~
=
Σ
Z
w
T
e
~
w with e
~
w
=
(
(
y
1
-
μ
~
1
)
/
v
(
μ
~
1
)
,
⋯
,
(
y
n
-
μ
~
n
)
/
v
(
μ
~
n
)
)
T. Then (41)
∂
2
Q
p
β
∣
ω
∂
β
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
w
∂
μ
~
T
T
∂
μ
~
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
K
I
-
Z
Σ
Z
T
K
-
1
δ
i
b
^
k
+
Z
Σ
d
k
e
^
i
=
σ
-
2
D
T
Ω
1
-
1
δ
i
b
^
k
+
Z
Σ
d
k
e
^
iHence, (42)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
δ
1
b
^
k
δ
2
b
^
k
⋯
δ
n
b
^
k
+
Z
Σ
d
k
e
^
1
d
k
e
^
2
⋯
d
k
e
^
n
=
σ
-
2
D
T
Ω
1
-
1
I
n
b
^
k
+
Z
Σ
d
k
e
^
T
.


## 3.1. Case-Weights Perturbation

Letω be an n
×
1 perturbation vector such that ω
0
=
(
1
,
⋯
,
1
)
T. The joint log quasi-likelihood function for the perturbed model is given by (24)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
i
=
1
n
ω
i
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C
=
-
q
/
2
log
⁡
2
π
σ
2
-
1
/
2
Σ. Then the penalized quasi-likelihood function can be expressed as(25)
Q
p
β
∣
ω
=
∑
i
=
1
n
ω
i
∫
y
i
μ
~
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
i
=
h
(
x
i
,
β
)
+
z
i
T
b
~, Z
=
(
z
1
,
z
2
,
…
,
z
n
)
T, and b
~ satisfies (26)
∑
i
=
1
n
w
i
y
i
-
μ
~
i
v
μ
~
i
z
i
-
Σ
-
1
b
~
=
0
.Hence, b
~
=
Σ
Z
T
W
¯
e
~, where W
¯
=
d
i
a
g
(
ω
1
,
…
,
ω
n
). Then (27)
∂
2
Q
p
β
∣
w
∂
β
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
δ
i
e
^
i
+
σ
-
2
D
T
K
·
K
-
1
Ω
-
1
K
-
1
-
I
δ
i
e
^
i
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
δ
i
e
^
i
,and thus (28)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
E
,where E
=
d
i
a
g
(
e
^
1
,
…
,
e
^
n
).


## 3.2. Response Variable Perturbation

A perturbation of the response variables(
y
1
,
⋯
,
y
n
)
T is introduced by replacing y
i by y
i
ω
=
y
i
+
ω
i, where ω
=
(
ω
1
,
⋯
,
ω
n
)
T, and ω
0
=
(
0
,
…
,
0
)
T represents the situation with no perturbation. In this case, the joint log quasi-likelihood function for the perturbed model is given by (29)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
i
=
1
n
∫
y
i
+
ω
i
μ
i
y
i
+
ω
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C is a constant. It follows from Section 2 that the penalized quasi-likelihood function is (30)
Q
p
β
∣
ω
=
∑
i
=
1
n
∫
y
i
+
ω
i
μ
~
i
y
i
+
ω
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
i
=
h
(
x
i
,
β
)
+
z
i
T
b
~, and b
~ satisfies (31)
∑
i
=
1
n
y
i
+
w
i
-
μ
~
i
v
μ
~
i
z
i
-
Σ
-
1
b
~
=
0
.It follows that b
~
=
Σ
Z
T
(
e
~
+
W
v
), where W
v
=
(
w
v
1
,
…
,
w
v
n
)
T with w
v
i
=
(
w
i
/
v
μ
i
)
μ
i
=
μ
~
i. Then (32)
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
∂
ω
i
+
∂
W
v
∂
ω
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
∂
μ
~
T
T
∂
μ
~
∂
ω
i
+
δ
i
T
1
v
μ
~
i
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
-
Ω
1
δ
i
T
1
v
μ
~
i
+
δ
i
T
1
v
μ
~
i
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
δ
i
T
1
v
μ
~
i
,and (33)
△
=
∂
2
Q
p
β
∣
w
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
K
-
1
E
∗
,where E
∗
=
d
i
a
g
(
e
1
∗
,
…
,
e
n
∗
) with e
i
∗
=
1
/
v
(
μ
~
i
) and μ
^
i
=
h
(
x
i
,
β
^
)
+
z
i
T
b
^.


## 3.3. Explanatory Variables Perturbation

In this case, we focus on the perturbation of a specific explanatory variable. Under this condition we have the perturbed explanatory matrixX
ω
=
(
x
1
,
⋯
,
x
t
ω
,
⋯
,
x
n
)
T with x
t
ω
=
x
t
+
ω, where x
t is a single explanatory variable of matrix X
ω corresponding to y
t and ω
0
=
(
0
,
…
,
0
)
T denotes no perturbation. Then the joint log quasi-likelihood function for the perturbed model is (34)
Q
β
,
σ
2
;
y
,
b
,
ω
=
C
+
∑
j
=
1
n
∫
y
i
μ
j
w
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
Σ
-
1
b
,where C is a constant, μ
j
ω
=
h
(
x
j
,
β
)
+
z
j
T
b
(
j
≠
i
), and μ
t
ω
=
h
(
x
t
ω
β
)
+
z
t
T
b. It follows from Section 2 that (35)
Q
p
β
∣
ω
=
∑
j
=
1
n
∫
y
i
μ
~
j
w
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
Σ
-
1
b
~
,where μ
~
j
ω
=
h
(
x
j
,
β
)
+
z
j
T
b
~
(
j
≠
t
), μ
~
t
ω
=
h
(
x
t
ω
,
β
)
+
z
t
T
b
~, and b
~ satisfies (36)
∑
j
=
1
n
y
i
-
μ
~
j
w
v
μ
~
j
w
z
i
-
Σ
-
1
b
~
=
0
.Therefore, b
~
=
Σ
Z
T
e
~
w and e
~
w
=
(
(
y
1
-
μ
~
1
)
/
v
(
μ
~
1
)
,
⋯
,
(
y
n
-
μ
~
n
)
/
v
(
μ
~
n
)
)
T. Let h
t
T
=
∂
h
(
x
t
,
β
)
/
∂
x
t
T
β
=
β
^ and H
b
t
=
∂
2
h
(
β
)
/
∂
β
∂
x
t
T
β
=
β
^. Then (37)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
e
^
ω
T
∂
∂
ω
T
∂
h
β
∂
β
ω
0
,
β
^
+
σ
-
2
D
T
∂
e
~
w
∂
μ
~
∂
μ
~
∂
ω
T
ω
0
,
β
^
=
σ
-
2
e
^
T
H
b
t
+
σ
-
2
D
T
Ω
1
-
1
δ
t
h
t
T
,where ·
[
·
] indicates the array multiplication.


## 3.4. Perturbation of Covariates in Random Effects

Consider perturbing the data for thekth explanatory variable of Z, by modifying the data matrix Z as Z
w
=
Z
+
ω
d
k
T, where d
k is a q
-vector with 1 at kth position and zeros elsewhere. Under this situation, the perturbed joint log quasi-likelihood can be expressed as (38)
Q
β
,
σ
2
;
y
,
b
∣
ω
=
C
+
∑
i
=
1
n
∫
y
i
μ
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
T
Σ
-
1
b
,where C is a quantity that does not depend on β and ω, and μ
=
h
(
β
)
+
Z
w
b. When ω
0
=
0, it indicates no perturbation. It follows from Section 2 that (39)
Q
p
β
∣
ω
=
∑
i
=
1
n
∫
y
i
μ
~
i
y
i
-
t
σ
2
v
t
d
t
-
1
2
σ
2
b
~
T
Σ
-
1
b
~
,where μ
~
=
h
(
β
)
+
Z
w
b
~, and b
~ satisfies (40)
Z
w
T
e
~
w
-
Σ
-
1
b
~
=
0
,and therefore, b
~
=
Σ
Z
w
T
e
~
w with e
~
w
=
(
(
y
1
-
μ
~
1
)
/
v
(
μ
~
1
)
,
⋯
,
(
y
n
-
μ
~
n
)
/
v
(
μ
~
n
)
)
T. Then (41)
∂
2
Q
p
β
∣
ω
∂
β
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
∂
e
~
w
∂
μ
~
T
T
∂
μ
~
∂
w
i
ω
0
,
β
^
=
σ
-
2
D
T
K
I
-
Z
Σ
Z
T
K
-
1
δ
i
b
^
k
+
Z
Σ
d
k
e
^
i
=
σ
-
2
D
T
Ω
1
-
1
δ
i
b
^
k
+
Z
Σ
d
k
e
^
iHence, (42)
△
=
∂
2
Q
p
β
∣
ω
∂
β
∂
ω
T
ω
0
,
β
^
=
σ
-
2
D
T
Ω
1
-
1
δ
1
b
^
k
δ
2
b
^
k
⋯
δ
n
b
^
k
+
Z
Σ
d
k
e
^
1
d
k
e
^
2
⋯
d
k
e
^
n
=
σ
-
2
D
T
Ω
1
-
1
I
n
b
^
k
+
Z
Σ
d
k
e
^
T
.


## 4. Numerical Results

To illustrate how to use the proposed methodology, we consider the data set reported by [27]. The data came from a study of the pharmacokinetics of indomethacin following bolus intravenous injection of the same dose in six human volunteers. For each subject, the plasma concentrations of indomethacin were measured at 11 time points from 15 min to 8 hours postinjection. Davidian et al. used nonlinear repeated model to analyze the dataset in [28]; we model it using the following QLNMWRE:(43)
μ
i
j
=
f
x
i
j
,
β
+
b
i
i
=
1
,
⋯
,
6
;
j
=
1
,
⋯
,
11
,where response variables y
i
j
|
b
i belong to the Gumbel distribution (cf. [29]) with the density function (44)
p
y
i
j
∣
b
i
=
exp
⁡
y
i
j
-
θ
-
exp
⁡
y
i
j
-
θ
,
-
∞
<
y
i
j
<
∞
,
-
∞
<
θ
<
∞
,
b
i
~
N
(
0
,
a
), and f
(
x
,
β
)
=
e
β
1
exp
⁡
(
-
e
β
2
x
)
+
e
β
3
exp
⁡
(
-
e
β
4
x
). By [29], we have E
(
y
)
=
θ
-
γ and V
a
r
(
y
)
=
π
2
/
6
=
σ
2
v
(
μ
), where γ
=
0.5772 is called the Euler constant, σ
2
=
π
2
/
6 and v
(
·
)
=
1. It is easily shown that Assumption A holds for our proposed model. Therefore, we can apply our proposed methodology to estimate the parameters in model (43). Using the algorithm in Section 2, we obtain the MQLE of β, the predictive values of b
1
,
⋯
,
b
6 as follows: β
^
=
(
0.8317,0.0446
,
-
13.2203
,
-
1.2535
)
T and (45)
b
^
=
-
0.1133,0.0825,0.1600,0.0500
,
-
0.0625,0.1133
T
.Now we present local influence analysis for the above fitting results. Under case-weight perturbation, cases 23, 45, and 56 are most influential, as depicted as in Figure1(a). Cases 1, 12, 23, 45, and 56 are identified as influential points, and case 23 is the most influential, as shown in Figure 1(b). The index plots of -
F
¨ and l
m
a
x for perturbation on explanatory variables are given in Figures 2(a) and 2(b), respectively. From Figure 2(a) we can see that cases 12, 23, 45, and 56 are identified as influential points. Figure 2(b) shows that cases 1, 12, 23, 34, 45, and 56 are influential. Figure 3 displays the index plots of l
m
a
x for the perturbation of random effects. For these types of perturbation, case 23 is identified as being the most influential. Note that case 23 exerts great influence in each perturbation scheme, which indicates that the results obtained through different perturbation schemes are quite consistent. Special attention should be paid to those influential cases, which may be worthwhile to consider a more formal test to check whether they are outliers.Figure 1
Index plots of -F
i
i and l
m
a
x for case-weights perturbation.
(a)
(b)Figure 2
Index plots of perturbation of explanatory variables.
(a)
(b)Figure 3
Index plots ofl
m
a
x for perturbation of random effects design matrix.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusion

In this work, we have assessed the local influence of minor perturbations of our proposed models. The key idea of the previous approach is to study the behavior of the likelihood displacement obtained from a relevant perturbation. However, it is difficult to apply it directly to the proposed model due to the fact that the marginal quasi-likelihood function of the QLNMWRE involves the intractable integral. To solve this problem, we have employed Laplace’s method to approximate the marginal quasi-likelihood function of the QLNMWRE, which results in the penalized quasi-likelihood (PQL). Based on the PQL and the penalized quasi-likelihood displacement, the estimates of unknown parameters have been proposed, and local influence analysis has been investigated. Our numerical example has demonstrated that our proposed local influence technique is rather useful in the detection of influential points. Although the focus of this article is on the assessment of influential points in the QLNMWRE, the local influence approach can be extended to other complicated models.



---

*Source: 4878925-2018-08-08.xml*