# Likelihood Inference of Nonlinear Models Based on a Class of Flexible Skewed Distributions

Authors: Xuedong Chen; Qianying Zeng; Qiankun Song
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper deals with the issue of the likelihood inference for nonlinear models with a flexible skew-t-normal (FSTN) distribution, which is proposed within a general framework of flexible skew-symmetric (FSS) distributions by combining with skew-t-normal (STN) distribution. In comparison with the common skewed distributions such as skew normal (SN), and skew-t (ST) as well as scale mixtures of skew normal (SMSN), the FSTN distribution can accommodate more flexibility and robustness in the presence of skewed, heavy-tailed, especially multimodal outcomes. However, for this distribution, a usual approach of maximum likelihood estimates based on EM algorithm becomes unavailable and an alternative way is to return to the original Newton-Raphson type method. In order to improve the estimation as well as the way for confidence estimation and hypothesis test for the parameters of interest, a modified Newton-Raphson iterative algorithm is presented in this paper, based on profile likelihood for nonlinear regression models with FSTN distribution, and, then, the confidence interval and hypothesis test are also developed. Furthermore, a real example and simulation are conducted to demonstrate the usefulness and the superiority of our approach.


---

## Body

## 1. Introduction

The common assumption of distribution for random error is normal in statistical modeling. This assumption may lack the robustness against departures from normality and/or outliers and may result in misleading inferential results [1, 2]. For the past few years, there is an increasing interest in developing more flexible parametric families capable of adopting as closely as possible real data, which exhibit quite substantial nonnormal characteristics such as skewness and heavy tails. In a variety of applications, one popular option is to modify a symmetric probability density function of a variable, thereby introducing skewness. An important advantage of this sort of approach compared with other approaches to robustness is an explicit statement of the probabilistic setting, leading to a clear interpretation of the results [3]. Following this idea, the skw-normal (SN) distribution was firstly introduced by [4], and, then, the skew-t (ST) distribution was introduced by [5]; the skew-t-normal (STN) was introduced by [6]; moreover, some extensions to these multivariate cases were studied by [7, 8] and so on. Since then, several authors have tried to extend these results to more general forms of skew-symmetric distributions, of which here we would like to mention [9], in this paper; they proposed a general framework of distributions which is called flexible skew-symmetric (FSS) distribution. As pointed out by [10] that this distribution family enjoys a sufficient flexibility in that with different choice of submodel settings, the FSS distribution includes several known distributions such as the SN and ST as its special cases.However, in many practical applications, it is not rare at all to encounter a multimodality, sometimes with an even irregular shape, and, for this case, all the distributions mentioned above appear to be unsufficient to describe the multimodal feature of the data. A solution to this problem is to use finite mixture models. In [11], the authors worked with a mixture model with component densities belonging to the STN distribution and a computationally feasible EM-type algorithm was developed for calculating the maximum likelihood (ML) estimates of parameters. Unfortunately, although the proposed methodology is useful for analyzing multimodal asymmetric data, it suffers from the problem of “model identification” as the number of the parameters to be estimated is usually large. As a result, in this paper, we deal with a new extension of the class of FSS distributions, which is referred to as flexible skew-t-normal (FSTN) distribution. This new distribution is proposed within the general framework of the FSS distributions in combination with the definition of STN distribution. In practical applications, it is able to regulate the density in a more flexible way to offer robustness and it can be treated as an appealing option for accommodating data with skewness and heavy tails as well as multimodality jointly.On the other hand, nonlinear regression models are widely applied in the fields of economics, engineering, biomedical research, and so forth, where the nonlinear function of unknown parameters is used to explain or investigate the nonlinear relationship of random phenomena under study. More recently, several authors have used a class of skewed distributions in the context of nonlinear regression models, and some valuable results were obtained. For example, [12] developed the robust estimation and the local influence analysis for regression model with SMSN distribution. From Bayesian point of view, [13] considered the Bayesian estimation and the case influence diagnostics for nonlinear regression models with SMSN distributions. More related literature could be found in [14–17]. Generally speaking, for model fitting of the nonlinear regression with skewed distributions, a popular approach is to consider the hierarchical representation of variables with a specific distribution, in which the postulated distribution is expressed as several conditional distributions of simpler forms such as normal and Student’s t and Gamma. Based on that, EM algorithm or Bayesian hierarchical approach then can be implemented effectively for conducting model estimation and statistical inference.In this paper, our aim is to develop an approach to likelihood inference of nonlinear regression models with FSTN assumption. As there is no stochastic representation for FSTN distribution, all the methods cited above become unavailable for our considered problem and an alternative way is to return to the original Newton-Raphson iterative procedure for model estimation. Under the nonlinear regression paradigm, the accuracy of estimates is affected by the strength of nonlinearity and the corresponding confidence interval and hypothesis test require the assumption of normality of the estimators or distribution, which is too restrictive. Besides, considering that in many practical applications, rather than the total parameters, we are usually interested in a proper subset of them. By taking all these factors into account, in this paper we focus on the parameters of interest and propose a modified Newton-Raphson iterative algorithm for calculating the ML estimates based on profile likelihood. Furthermore, the confidence interval and hypothesis test for the parameters of interest are also considered. We conduct an application and a simulation study to compare the algorithm effectiveness and distribution robustness for nonlinear regression model in terms of fitting performance and model selection. The results from the numerical examples illustrate the usefulness and the superiority of our methodology.The remainder of this paper is organized as follows. In Section2, we briefly discuss the FSS distribution and FSTN distributions. In Section 3, we present the likelihood inference including the quantities of the first- and the second-order derivatives as well as the standard Newton-Raphson iterative formula. In Section 4, we give an introduction of profile inference for our proposed model, where the confidence estimation and hypothesis test are presented too. Section 5 gives numerical examples using both simulated and real data to illustrate the performance of the proposed methodology. Finally, some concluding remarks are given in Section 6.


## 2. Models and Notation

The class of skewed distributions such as SN, ST, and STN perform to be plausible for modeling skewness or (and) heavy tails underlying the observations. The actual situation is that it is not rare at all to encounter multimodality, sometimes with an even more irregular shape, and, for this case, the aforementioned distributions become unsufficient. In this paper, with the adoption of a sufficiently flexible class of distributions, we consider one of these extensions, referred to as the family of flexible skew-symmetric (FSS) distributions which is introduced by [9] with the following density function of type:
(1)
f
(
y
;
μ
,
ω
,
α
)
=
2
ω
-
1
f
0
(
z
)
G
{
P
K
(
z
)
}
,
where f
0 and G are symmetric univariate density and distribution function, respectively, that is, f
0
(
x
)
=
f
0
(
-
x
), G
(
-
x
)
=
1
-
G
(
x
), and P
K
(
x
)
=
α
1
x
+
α
3
x
3
+
⋯
α
2
K
-
1
x
2
K
-
1 is an odd polynomial of degree K (i.e., a polynomial including only terms of odd degree), z
=
ω
-
1
(
x
-
μ
), μ
∈
R is the location parameter, ω
>
0 is the scale parameter, and α
=
(
α
1
,
α
3
,
…
,
α
2
K
-
1
) and α
1
,
α
3
,
…
,
α
2
K
-
1
∈
R are shape parameters.In general, the density function of STN distribution can be represented as2
ω
-
1
t
(
z
;
v
)
Φ
(
α
z
), where t and Φ, respectively, denote the univariate standard Student’s t density function and the univariate standard normal distribution function and v is the degrees of freedom. The skewness is regulated by the shape parameter α and the tail thickness of the distribution is controlled by v. Commonly, in comparison with the SN distribution, the STN distribution exhibits obvious feature of heavy tails when v
≤
10.In this paper, we work with one version of (1) and the specific definition can be presented as follows. Let f
0
(
z
)
=
t
(
z
;
v
) and G
(
x
)
=
Φ
(
x
), where Φ
(
x
) is defined as before and
(2)
t
(
z
;
v
)
=
Γ
(
v
+
1
)
/
2
ω
v
π
1
/
2
Γ
(
v
/
2
)
1
+
z
2
v
-
v
+
1
/
d
,
that is, the density function of univariate Student’s t distribution with 0 location, ω scale, and v degrees of freedom. The above extension of (1) is referred to as flexible skew-t-normal (FSTN) distribution, denoted by FSTN
(
μ
,
ω
,
α
,
v
).It is noted that the FSTN distribution is proposed within the general framework of FSS distribution by combining with the definition of STN distribution and, as a consequence, it shares analogous feature with these two distributions. For all of that, the FSTN distribution presents some interesting and peculiar features and is able to regulate the density in a more flexible way. To be particular, except for modeling skewness and tail thickness, the FSTN distribution allows for multimodality, depending on the specific setting ofα. For the purpose of comparison and illustration, we assume f
0
(
z
)
=
ϕ
(
z
) and G
(
z
)
=
Φ
(
z
) in FSS distribution, which is denoted by FSN, and, then, we set K
=
2 in P
K
(
z
), that is, P
K
(
x
)
=
α
1
x
+
α
3
x
3 in STN distribution, and, for this case, different selections of α
1 and α
3 determine whether the density is unimodal or bimodal. Moreover, the same assumption for G and P in FSTN distribution is made.Figure1 displays the density functions of FSN and STN as well as FSTN distributions with four different situations considered, namely, α
1
=
1, α
3
=
0, and v
=
10; α
1
=
α
3
=
1 and v
=
6; α
1
=
1, α
3
=
-
1, and v
=
4; α
1
=
-
1, α
3
=
1, and v
=
4, respectively, with μ
=
1, ω
=
1.5 for all cases. By examination of Figure 1, we can detect how these three densities change with different combinations of α and v. For instance, in Figure 1(a), FSN, STN, and FSTN appear to be very close, while STN and FSTN are heavy tailed to a little extent. In Figure 1(b), both FSN and FSTN are unimodal when α
1 and α
3 keep the same sign, and the ranking for the degree of skewness is STN, FSTN, and FSN in turn; that is, the STN and FSTN distributions have thicker tails compared to FSN distribution. With opposite sign of α
1 and α
3, both FSN and FSTN distributions are bimodal and highly skewed in Figures 1(c) and 1(d); moreover, in the same direction, FSTN has thicker tails than FSN distribution. Our proposed FSTN distribution can be treated as a proper compromise between the FSS distribution and the STN distribution. It allows for a wider range of tail behavior compared to FSS distribution whereas it is able to accommodate multimodality which cannot be described by STN distribution. From the applied viewpoint, the FSTN distribution is an appealing option which can be expected to yield robust inferential results in the presence of outlying observations.((a)–(d)) Density functions of FSN, STN, and FSTN distributions with four different situations.
(a)
(b)
(c)
(d)


## 3. Likelihood Inference

Considern independent observations satisfying a nonlinear regression model as [
y
i
]
~
FSTN
(
μ
i
,
ω
,
α
,
v
) with μ
i
=
η
(
x
i
,
β
) for i
=
1,2
,
…
,
n. Here, x
i is a p-dimensional vector and β is a p
×
1 vector of parameters. Also, let X
=
(
x
1
,
…
,
x
n
)
T be the n
×
p design matrix; η
(
·
,
·
) is a known twice differentiable function. Then, the corresponding log-likelihood for parameter θ
=
(
β
,
α
,
ω
,
v
) is given by
(3)
l
(
θ
)
=
constant
-
n
log
⁡
ω
+
∑
i
=
1
n
‍
log
⁡
t
(
z
i
)
+
∑
i
=
1
n
‍
log
⁡
Φ
{
P
K
(
z
i
)
}
.
And we have
(4)
∂
l
(
θ
)
∂
β
=
-
1
ω
∑
i
=
1
n
‍
t
′
(
z
i
)
t
(
z
i
)
+
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
,
∂
l
(
θ
)
∂
ω
=
-
1
ω
∑
i
=
1
n
‍
t
′
(
z
i
)
t
(
z
i
)
z
i
+
z
i
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
+
n
,
∂
l
(
θ
)
∂
v
=
∑
i
=
1
n
‍
z
i
2
t
(
z
i
)
+
n
2
D
G
v
+
1
2
-
D
G
v
2
-
1
2
,
∂
l
(
θ
)
∂
α
2
j
-
1
=
∑
i
=
1
n
‍
z
i
2
j
-
1
ϕ
P
K
z
i
Φ
P
K
z
i
P
K
′
z
i
,
j
=
1,2
,
…
,
K
,
where η
i
=
η
(
x
i
,
β
) and D
G
(
x
)
=
d
[
log
⁡
Γ
(
x
)
]
/
d
x. The corresponding second-order derivatives of (3) can be shown as
(5)
∂
2
l
θ
∂
β
∂
β
T
=
∑
i
=
1
n
‍
∂
η
i
∂
β
T
1
ω
2
h
(
z
i
)
+
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
-
1
ω
t
′
(
z
i
)
t
(
z
i
)
+
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
2
η
i
∂
β
∂
β
T
,
∂
2
l
θ
∂
β
∂
ω
=
∑
i
=
1
n
‍
1
ω
2
t
′
z
i
t
z
i
+
P
K
′
z
i
ϕ
P
K
z
i
Φ
P
K
z
i
+
z
i
ω
2
h
(
z
i
)
+
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
,
∂
2
l
(
θ
)
∂
β
∂
v
=
-
1
ω
∑
i
=
1
n
‍
(
1
+
v
-
1
)
z
i
3
-
t
(
z
i
)
z
i
2
v
2
t
2
(
z
i
)
∂
η
i
∂
β
,
∂
2
l
θ
∂
β
∂
α
2
j
-
1
=
∑
i
=
1
n
‍
z
i
2
j
-
1
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
×
∂
η
i
∂
β
,
j
=
1,2
,
…
,
K
,
∂
2
l
(
θ
)
∂
ω
2
=
1
ω
2
∑
i
=
1
n
‍
t
′
(
z
i
)
t
(
z
i
)
+
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
3
z
i
+
z
i
2
h
(
z
i
)
+
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
+
P
K
′′
z
i
ϕ
P
K
z
i
Φ
P
K
z
i
,
∂
2
l
(
θ
)
∂
ω
∂
v
=
-
1
ω
∑
i
=
1
n
‍
(
1
+
v
-
1
)
z
i
3
-
t
(
z
i
)
z
i
2
v
2
t
2
(
z
i
)
,
∂
2
l
(
θ
)
∂
ω
∂
α
2
j
-
1
=
-
∑
i
=
1
n
‍
1
ω
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
z
i
×
z
i
2
j
-
1
,
j
=
1,2
,
…
,
K
,
∂
2
l
(
θ
)
∂
v
2
=
n
2
T
G
v
+
1
2
-
T
G
v
2
+
∑
i
=
1
n
‍
log
⁡
t
(
z
i
)
,
∂
2
l
(
θ
)
∂
α
2
j
-
1
∂
α
2
k
-
1
=
∑
i
=
1
n
‍
z
i
j
+
k
Q
{
P
K
(
z
i
)
}
,
j
,
k
=
1,2
,
…
,
K
,
∂
2
l
(
θ
)
∂
v
∂
α
2
j
-
1
=
0
,
j
=
1,2
,
…
,
K
,
where T
G
(
x
)
=
d
2
[
log
⁡
Γ
(
x
)
]
/
d
x
2, h
(
x
)
=
t
′′
(
x
)
t
(
x
)
-
[
t
′
(
x
)
]
2
/
t
2
(
x
), and Q
(
x
)
=
ϕ
(
x
)
/
Φ
(
x
)
[
x
+
ϕ
(
x
)
/
Φ
(
x
)
].AssumeU
(
θ
)
=
(
U
T
(
β
)
,
U
(
ω
)
,
U
(
v
)
,
U
T
(
α
)
)
T is the gradient or score vector and H
(
θ
) is the Hessian matrix composed of the above second-order derivatives. To obtain the ML estimate of θ, the Newton-Raphson iteration algorithm is defined by
(6)
θ
^
(
k
+
1
)
=
θ
^
(
k
)
-
[
H
(
θ
)
θ
=
θ
^
(
k
)
-
1
U
(
θ
)
θ
=
θ
^
(
k
)
.
It is noted that the above iterative procedure is an unpartitioned algorithm; that is, all the parameters including nonlinear regression coefficients β, scale parameter ω, and shape parameter α as well as tail thickness parameter v are estimated simultaneously. For our considered problem, there are at least two difficulties that may be encountered for (6); the first one is that once the number of the parameters to be estimated becomes large, the corresponding computational burden turns to be heavy with an unacceptable estimation error, and the second one is as follows: when the strength of nonlinearity of the link function η
(
·
,
·
) changes, the iterative process may become unstable or even nonconvergent, leading to the poor estimation results. Considering the needs of practical problems, rather than the total parameter set, we are usually interested in a proper subset of it. To improve the efficiency of the algorithm and to facilitate statistical inference of the nonlinear models with FSTN distribution, we put forward the following profile likelihood method based on (3) and (6).


## 4. Profile Likelihood Inference

### 4.1. Profile Estimation Algorithm

Letθ
=
(
β
T
,
γ
T
) be a partition of θ, where β is a parameter vector of nonlinear regression, an interest parameter, and γ is a q-dimension (
q
=
K
+
2
) nuisance parameter with γ
=
(
α
,
ω
,
v
)
T. Similarly, the partition of H and U is given as H
(
θ
)
=
H
11
H
12
H
21
H
22, where H
11
=
∂
2
l
(
θ
)
/
∂
β
∂
β
T, H
12
=
∂
2
l
(
θ
)
/
∂
β
∂
γ
T
=
(
∂
2
l
(
θ
)
/
∂
β
∂
ω
,
∂
2
l
(
θ
)
/
∂
β
∂
v
,
∂
2
l
(
θ
)
/
∂
β
∂
α
1
,
…
,
∂
2
l
(
θ
)
/
∂
β
∂
α
2
k
-
1
), H
12
=
H
21
T, and H
22
=
∂
2
l
(
θ
)
/
∂
γ
∂
γ
T, and the diagonal elements of H
22 are given by ∂
2
l
(
θ
)
/
∂
ω
2
,
∂
2
l
(
θ
)
/
∂
v
2
,
∂
2
l
(
θ
)
/
∂
α
1
2
,
…
,
∂
2
l
(
θ
)
/
∂
α
2
k
-
1
2, respectively, and the off-diagonal elements in H
22 can be obtained similarly. Let U
(
θ
)
=
(
U
1
,
U
2
), corresponding to the partition of θ.In the subsequent context, we focus on the estimation ofθ based on profile likelihood method [18]. Firstly, suppose β is known and we rewrite the original likelihood function (3) as
(7)
l
(
θ
)
=
l
(
β
,
γ
)
=
l
β
(
γ
)
,
where notation l
β
(
γ
) denotes that β is fixed but γ varies. For each β, to estimate γ we can obtain
(8)
γ
^
β
=
arg
max
⁡
γ
⁡
l
β
(
γ
)
.
Alternatively, to estimate β, we evaluate the maximum value of l
β
(
γ
) over γ
^
β and have
(9)
β
^
p
=
arg
max
⁡
β
⁡
l
β
(
γ
^
β
)
=
arg
max
⁡
β
⁡
l
(
β
,
γ
^
β
)
,
where l
(
β
,
γ
^
β
) and β
^
p are referred to as the profile likelihood function and the profile ML estimation, respectively.Following [19], we define the profile Newton-Raphson iteration formula as follows:
(10)
β
^
p
(
k
+
1
)
=
β
^
p
(
k
)
-
H
β
∣
γ
-
1
U
β
∣
γ
,
γ
^
p
(
k
+
1
)
=
γ
^
p
(
k
)
-
H
22
-
1
U
1
+
H
22
-
1
H
21
H
β
∣
γ
-
1
U
β
∣
γ
,
where U
β
∣
γ
=
U
1
-
H
12
(
H
22
)
-
1
U
2, H
β
∣
γ
=
H
11
-
H
12
(
H
22
)
-
1
H
21, and all the matrices and the vectors on the right hand side of (10) are evaluated at β
^
p
(
k
) and γ
^
p
(
k
).Note that, both in (6) and in (10), the strength of nonlinearity of link function η
(
·
,
·
) is reflected by ∂
2
η
i
/
∂
β
∂
β
T to large extents. Therefore, by examination of the expression of ∂
2
l
(
θ
)
/
∂
β
∂
β
T, we find that when the element of
(11)
ω
t
′
(
z
i
)
t
(
z
i
)
+
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
2
η
i
∂
β
∂
β
T
is much less than the corresponding element of
(12)
∂
η
i
T
∂
β
h
z
i
+
P
K
′
z
i
Q
P
K
z
i
+
P
K
′′
z
i
ϕ
P
K
z
i
Φ
P
K
z
i
×
ϕ
P
K
z
i
Φ
P
K
z
i
∂
η
i
∂
β
,
the following approximating result can be obtained as
(13)
H
11
*
=
∑
i
=
1
n
∂
η
i
∂
β
T
1
ω
2
h
(
z
i
)
+
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
T
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
.
And then, the iteration formulas of U
β
∣
γ
*, [
H
β
∣
γ
*
]
-
1, β
^
p, and β
^
p
* can also be obtained just as before. The above estimation procedure is referred to as profile modified Newton-Raphson iteration algorithm.The stopping rule for the above algorithm can be presented as that iteration proceeds until some distance involving two successive evaluations of the profile likelihoodl
(
β
,
γ
^
β
), such as l
(
β
^
(
k
+
1
)
,
γ
^
β
^
(
k
+
1
)
)
-
l
(
β
^
(
k
)
,
γ
^
β
^
(
k
)
) or l
(
β
^
(
k
+
1
)
,
γ
^
β
^
(
k
+
1
)
)
/
l
(
β
^
(
k
)
,
γ
^
β
^
(
k
)
)
-
1, is small enough; for example, ɛ
<
1
0
-
4 is adopted in this paper.The asymptotic covariance matrix of the ML estimates for profile likelihood can be evaluated by inverting the expected information matrix; however, it does not have a closed form expression; the observed information matrixJ
(
β
)
=
-
∂
2
l
(
β
,
γ
^
β
)
/
∂
β
∂
β
T can be used as a replacement which is estimated by H
11
-
H
12
(
H
22
)
-
1
H
21, where H
11, H
12, and H
22 can be obtained as similar as above.The choice of the initial values plays an important role in nonlinear regression fitting; in this paper, the specific steps for choosing the starting values are implemented as follows:(i)
compute the initial valueβ
^
(
0
) based on the nonlinear regression model with standard normal assumption;
(ii)
withβ
^
(
0
) fixed, compute the initial values ω
(
0
) and α
(
0
) for the SN and finite mixture SN assumptions, respectively.In order to simplify the estimation of parameterv, we have fixed integral values for v from 2 to 40 by one; choose the value of v that maximizes the profile likelihood as v
(
0
), and, then, the initial values of β
^
(
0
), ω
(
0
), α
(
0
), v
(
0
) that are required in the estimation procedure are all obtained.


### 4.2. Profile Confidence Estimation and Hypothesis Test

Confidence interval and hypothesis test play an important role in statistical inference and, in the subsequent content, we will consider the profile confidence estimation and hypothesis test for the parameters of interest in nonlinear regression models. SupposeJ
(
θ
)
=
E
[
-
H
(
θ
)
]; the following regular conditions for likelihood inference are assumed:(R-i)
∫
‍
∂
log
⁡
L
(
x
;
θ
)
/
∂
θ
·
L
(
x
;
θ
)
d
x
=
0;
(R-ii)
∂
/
∂
θ
∫
‍
L
(
x
;
θ
)
=
∫
‍
∂
L
(
x
;
θ
)
/
∂
θ
d
x
=
0;
(R-iii)
∂
/
∂
θ
∫
‍
g
(
x
)
L
(
x
;
θ
)
d
x
=
∫
‍
g
(
x
)
∂
L
(
x
;
θ
)
/
∂
θ
d
x, where g
(
x
) is arbitrary measurable function;
(R-iv)
E
∂
log
⁡
L
(
x
;
θ
)
/
∂
θ
2
>
0.Apart from the above assumptions, in this paper, some additional conditions are assumed to hold for that:(A-i)
η
(
x
i
,
β
) is twice continuously differentiable with respect to β
∈
R
p;
(A-ii)
D
=
(
∂
η
1
/
∂
β
,
…
,
∂
η
n
/
∂
β
)
T is full rank for all β;
(A-iii)
assumeθ
0 is the true value of θ; then, there exist a neighbour region N
0 of θ
0 and a constant c such that λ
min
⁡
(
J
(
θ
)
)
>
c for any θ
∈
N
0.Under the above assumption, following [20], we have
(14)
2
{
l
(
β
^
p
,
γ
^
β
^
p
)
-
l
(
β
0
,
γ
^
β
0
)
}
⟶
χ
2
(
p
)
,
where β
0 is the true value of β and the convergence is under the meaning of convergence in probability. Based on the profile likelihood theory, the confidence region of β with the confidence level 100
(
1
-
α
)
% is given by
(15)
C
(
Y
)
=
{
β
:
2
[
l
(
β
^
p
,
γ
^
β
^
p
)
-
l
(
β
,
γ
^
β
)
]
≤
χ
2
(
p
,
α
)
}
,
where χ
2
(
p
,
α
) denotes the upper α percentile of chi-square distribution with p degrees of freedom. Furthermore, the hypothesis test
(16)
H
0
:
β
=
β
0
against
H
1
:
β
≠
β
0
can also be considered and the corresponding test statistics are presented by 2
[
l
(
β
^
,
γ
^
β
^
)
-
l
(
β
0
,
γ
^
β
0
)
]. Unlike standard likelihood method, profile likelihood confidence intervals and hypothesis test do not need assumption of normality of the estimator or distribution which is too restrictive; they are based on an asymptotic chi-square distribution of the log profile likelihood ratio test statistics and these properties bring a lot of convenience and feasibility in practical computation.


## 4.1. Profile Estimation Algorithm

Letθ
=
(
β
T
,
γ
T
) be a partition of θ, where β is a parameter vector of nonlinear regression, an interest parameter, and γ is a q-dimension (
q
=
K
+
2
) nuisance parameter with γ
=
(
α
,
ω
,
v
)
T. Similarly, the partition of H and U is given as H
(
θ
)
=
H
11
H
12
H
21
H
22, where H
11
=
∂
2
l
(
θ
)
/
∂
β
∂
β
T, H
12
=
∂
2
l
(
θ
)
/
∂
β
∂
γ
T
=
(
∂
2
l
(
θ
)
/
∂
β
∂
ω
,
∂
2
l
(
θ
)
/
∂
β
∂
v
,
∂
2
l
(
θ
)
/
∂
β
∂
α
1
,
…
,
∂
2
l
(
θ
)
/
∂
β
∂
α
2
k
-
1
), H
12
=
H
21
T, and H
22
=
∂
2
l
(
θ
)
/
∂
γ
∂
γ
T, and the diagonal elements of H
22 are given by ∂
2
l
(
θ
)
/
∂
ω
2
,
∂
2
l
(
θ
)
/
∂
v
2
,
∂
2
l
(
θ
)
/
∂
α
1
2
,
…
,
∂
2
l
(
θ
)
/
∂
α
2
k
-
1
2, respectively, and the off-diagonal elements in H
22 can be obtained similarly. Let U
(
θ
)
=
(
U
1
,
U
2
), corresponding to the partition of θ.In the subsequent context, we focus on the estimation ofθ based on profile likelihood method [18]. Firstly, suppose β is known and we rewrite the original likelihood function (3) as
(7)
l
(
θ
)
=
l
(
β
,
γ
)
=
l
β
(
γ
)
,
where notation l
β
(
γ
) denotes that β is fixed but γ varies. For each β, to estimate γ we can obtain
(8)
γ
^
β
=
arg
max
⁡
γ
⁡
l
β
(
γ
)
.
Alternatively, to estimate β, we evaluate the maximum value of l
β
(
γ
) over γ
^
β and have
(9)
β
^
p
=
arg
max
⁡
β
⁡
l
β
(
γ
^
β
)
=
arg
max
⁡
β
⁡
l
(
β
,
γ
^
β
)
,
where l
(
β
,
γ
^
β
) and β
^
p are referred to as the profile likelihood function and the profile ML estimation, respectively.Following [19], we define the profile Newton-Raphson iteration formula as follows:
(10)
β
^
p
(
k
+
1
)
=
β
^
p
(
k
)
-
H
β
∣
γ
-
1
U
β
∣
γ
,
γ
^
p
(
k
+
1
)
=
γ
^
p
(
k
)
-
H
22
-
1
U
1
+
H
22
-
1
H
21
H
β
∣
γ
-
1
U
β
∣
γ
,
where U
β
∣
γ
=
U
1
-
H
12
(
H
22
)
-
1
U
2, H
β
∣
γ
=
H
11
-
H
12
(
H
22
)
-
1
H
21, and all the matrices and the vectors on the right hand side of (10) are evaluated at β
^
p
(
k
) and γ
^
p
(
k
).Note that, both in (6) and in (10), the strength of nonlinearity of link function η
(
·
,
·
) is reflected by ∂
2
η
i
/
∂
β
∂
β
T to large extents. Therefore, by examination of the expression of ∂
2
l
(
θ
)
/
∂
β
∂
β
T, we find that when the element of
(11)
ω
t
′
(
z
i
)
t
(
z
i
)
+
P
K
′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
2
η
i
∂
β
∂
β
T
is much less than the corresponding element of
(12)
∂
η
i
T
∂
β
h
z
i
+
P
K
′
z
i
Q
P
K
z
i
+
P
K
′′
z
i
ϕ
P
K
z
i
Φ
P
K
z
i
×
ϕ
P
K
z
i
Φ
P
K
z
i
∂
η
i
∂
β
,
the following approximating result can be obtained as
(13)
H
11
*
=
∑
i
=
1
n
∂
η
i
∂
β
T
1
ω
2
h
(
z
i
)
+
P
K
′
(
z
i
)
Q
{
P
K
(
z
i
)
}
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
T
+
P
K
′′
(
z
i
)
ϕ
{
P
K
(
z
i
)
}
Φ
{
P
K
(
z
i
)
}
∂
η
i
∂
β
.
And then, the iteration formulas of U
β
∣
γ
*, [
H
β
∣
γ
*
]
-
1, β
^
p, and β
^
p
* can also be obtained just as before. The above estimation procedure is referred to as profile modified Newton-Raphson iteration algorithm.The stopping rule for the above algorithm can be presented as that iteration proceeds until some distance involving two successive evaluations of the profile likelihoodl
(
β
,
γ
^
β
), such as l
(
β
^
(
k
+
1
)
,
γ
^
β
^
(
k
+
1
)
)
-
l
(
β
^
(
k
)
,
γ
^
β
^
(
k
)
) or l
(
β
^
(
k
+
1
)
,
γ
^
β
^
(
k
+
1
)
)
/
l
(
β
^
(
k
)
,
γ
^
β
^
(
k
)
)
-
1, is small enough; for example, ɛ
<
1
0
-
4 is adopted in this paper.The asymptotic covariance matrix of the ML estimates for profile likelihood can be evaluated by inverting the expected information matrix; however, it does not have a closed form expression; the observed information matrixJ
(
β
)
=
-
∂
2
l
(
β
,
γ
^
β
)
/
∂
β
∂
β
T can be used as a replacement which is estimated by H
11
-
H
12
(
H
22
)
-
1
H
21, where H
11, H
12, and H
22 can be obtained as similar as above.The choice of the initial values plays an important role in nonlinear regression fitting; in this paper, the specific steps for choosing the starting values are implemented as follows:(i)
compute the initial valueβ
^
(
0
) based on the nonlinear regression model with standard normal assumption;
(ii)
withβ
^
(
0
) fixed, compute the initial values ω
(
0
) and α
(
0
) for the SN and finite mixture SN assumptions, respectively.In order to simplify the estimation of parameterv, we have fixed integral values for v from 2 to 40 by one; choose the value of v that maximizes the profile likelihood as v
(
0
), and, then, the initial values of β
^
(
0
), ω
(
0
), α
(
0
), v
(
0
) that are required in the estimation procedure are all obtained.


## 4.2. Profile Confidence Estimation and Hypothesis Test

Confidence interval and hypothesis test play an important role in statistical inference and, in the subsequent content, we will consider the profile confidence estimation and hypothesis test for the parameters of interest in nonlinear regression models. SupposeJ
(
θ
)
=
E
[
-
H
(
θ
)
]; the following regular conditions for likelihood inference are assumed:(R-i)
∫
‍
∂
log
⁡
L
(
x
;
θ
)
/
∂
θ
·
L
(
x
;
θ
)
d
x
=
0;
(R-ii)
∂
/
∂
θ
∫
‍
L
(
x
;
θ
)
=
∫
‍
∂
L
(
x
;
θ
)
/
∂
θ
d
x
=
0;
(R-iii)
∂
/
∂
θ
∫
‍
g
(
x
)
L
(
x
;
θ
)
d
x
=
∫
‍
g
(
x
)
∂
L
(
x
;
θ
)
/
∂
θ
d
x, where g
(
x
) is arbitrary measurable function;
(R-iv)
E
∂
log
⁡
L
(
x
;
θ
)
/
∂
θ
2
>
0.Apart from the above assumptions, in this paper, some additional conditions are assumed to hold for that:(A-i)
η
(
x
i
,
β
) is twice continuously differentiable with respect to β
∈
R
p;
(A-ii)
D
=
(
∂
η
1
/
∂
β
,
…
,
∂
η
n
/
∂
β
)
T is full rank for all β;
(A-iii)
assumeθ
0 is the true value of θ; then, there exist a neighbour region N
0 of θ
0 and a constant c such that λ
min
⁡
(
J
(
θ
)
)
>
c for any θ
∈
N
0.Under the above assumption, following [20], we have
(14)
2
{
l
(
β
^
p
,
γ
^
β
^
p
)
-
l
(
β
0
,
γ
^
β
0
)
}
⟶
χ
2
(
p
)
,
where β
0 is the true value of β and the convergence is under the meaning of convergence in probability. Based on the profile likelihood theory, the confidence region of β with the confidence level 100
(
1
-
α
)
% is given by
(15)
C
(
Y
)
=
{
β
:
2
[
l
(
β
^
p
,
γ
^
β
^
p
)
-
l
(
β
,
γ
^
β
)
]
≤
χ
2
(
p
,
α
)
}
,
where χ
2
(
p
,
α
) denotes the upper α percentile of chi-square distribution with p degrees of freedom. Furthermore, the hypothesis test
(16)
H
0
:
β
=
β
0
against
H
1
:
β
≠
β
0
can also be considered and the corresponding test statistics are presented by 2
[
l
(
β
^
,
γ
^
β
^
)
-
l
(
β
0
,
γ
^
β
0
)
]. Unlike standard likelihood method, profile likelihood confidence intervals and hypothesis test do not need assumption of normality of the estimator or distribution which is too restrictive; they are based on an asymptotic chi-square distribution of the log profile likelihood ratio test statistics and these properties bring a lot of convenience and feasibility in practical computation.


## 5. Simulation Study

To investigate the experimental performance of our methodology, we undertake a simulation study to compare the fitting performance of misspecified distribution as well as understanding the large sample properties of the ML estimates. To realize this purpose, we generate the artificial data from the following two nonlinear models:Model(
1
): drug responsiveness model:
(17)
y
i
=
β
0
-
β
0
1
+
x
i
/
β
1
β
2
+
ε
i
,Model(
2
): nonlinear growth-curve model:
(18)
y
i
=
β
1
1
+
β
2
exp
⁡
(
-
β
3
x
)
+
ε
i
,
with two kinds of skewed distributions for random error as follows: Case (I): ɛ
~
STN
(
μ
,
ω
,
α
,
v
) and Case (II): ɛ
~
FSTN
(
μ
,
ω
,
α
,
v
).The true parameter population is chosen as follows:β
=
(
2,5
,
-
0.8
) is assumed to be the same for the above two models. In Model (
1
), let μ
=
δ
1, ω
=
1.5, α
=
1.5, and v
=
4 and in Model (
2
) let μ
=
δ
2, ω
=
1.5, α
=
(
1
,
-
0.8
), and v
=
4. Note that, for the values of the degrees of freedom v, a relative value (v
=
4) can yield a heavy-tailed distribution as we need.In this simulation, there are totally four simulated data sets corresponding to two nonlinear models, Model(
1
) and Model (
2
), along with two distributions, Case (I) and Case (II), for random error. Similar to previous analysis, each simulated data set is fitted under STN and FSN as well as FSTN scenarios using three different estimation algorithms. To shed light on the experimental performance of our methodology, an interesting comparison can be made by examining how often we can recognize the true model.Table1 shows the absolute value of the average bias between the true and the estimated parameters and the percentages of each model being ranked as the best model based on AIC criterion out of 500 replications are also presented in Table 1.Table 1
Results of the simulation study.
Method
Model
Error
Fittingmodel
Bias and perc.
β
0
β
1
β
2
(%)
NR
1
Case I
STN
0.048
0.044
0.034
38.4
FSN
0.052
0.043
0.031
30.2
FSTN
0.051
0.048
0.035
31.4
Case II
STN
0.941
0.445
0.374
10.5
FSN
0.035
0.014
0.011
43.8
FSTN
0.033
0.010
0.010
45.7
2
Case I
STN
0.067
0.061
0.053
36.2
FSN
0.070
0.057
0.055
30.2
FSTN
0.073
0.072
0.048
33.6
Case II
STN
0.963
0.469
0.391
8.5
FSN
0.047
0.023
0.017
44.6
FSTN
0.043
0.022
0.021
46.9
P-NR
1
Case I
STN
0.049
0.040
0.032
32.6
FSN
0.050
0.047
0.038
40.1
FSTN
0.052
0.049
0.036
27.3
Case II
STN
1.029
0.420
0.248
11.4
FSN
0.039
0.018
0.010
39.6
FSTN
0.038
0.017
0.009
49.0
2
Case I
STN
0.070
0.050
0.047
30.3
FSN
0.067
0.057
0.056
42.3
FSTN
0.073
0.071
0.050
27.4
Case II
STN
1.052
0.446
0.265
10.2
FSN
0.053
0.031
0.020
39.6
FSTN
0.052
0.022
0.017
50.2
MP-NR
1
Case I
STN
0.045
0.049
0.033
32.3
FSN
0.051
0.047
0.036
38.3
FSTN
0.045
0.044
0.031
29.4
Case II
STN
0.773
0.373
0.266
14.5
FSN
0.031
0.018
0.011
38.4
FSTN
0.028
0.010
0.008
47.1
2
Case I
STN
0.068
0.072
0.062
28.1
FSN
0.076
0.074
0.062
41.2
FSTN
0.073
0.060
0.054
30.7
Case II
STN
0.794
0.394
0.290
11.3
FSN
0.037
0.026
0.017
39.4
FSTN
0.036
0.021
0.013
49.3By examination of Table1, we can find the following: (i) the difference for the estimation results between drug responsiveness model and growth-curve model is significant, indicating that the nonlinearity of model imposes an impact on parameter estimation; (ii) when the true distribution for random error is STN, the three fitting distributions have roughly the same behavior; it is hard to tell which approach for parameter estimation is good and which approach is bad and the similar results can be obtained for model selection based on AIC; (iii) when the true distribution for random error is FSTN, the FSN and the FSTN distributions outperform the STN distribution by producing estimates with lower bias and higher AIC proportion; furthermore, the PNR and MPNR methods perform better than the traditional NR method in general.To study the consistence properties of ML estimate, we focus on the situation that the true distribution for random error is Case (II) whereas the fitting distribution is FSTN too. For this case, two estimation algorithms, PNR and MPNR, are adopted for parameter estimation and samples of different sizes (n
=
50, 100, and 200) are generated from Models (
1
) and (
2
). We compute the 95% confidence interval for the parameters of interest and the mean square error (MSE) for different model, where MSE
(
β
)
=
1
/
500
∑
i
=
1
500
β
(
i
)
-
β
2. The length of the 95% confidence interval and the MSE results are summarized in Table 2. From Table 2 we can see that the length of the confidence interval and the MSE tend to decrease as the sample size increases as expected.Table 2
The table shows the length of 95% confidence interval of parameters of interest and MSE of specified model.
Model
Sample size
Profile-NR
Modified profile-NR
β
0
β
1
β
2
MSE
β
0
β
1
β
2
MSE
1
50
1.015
0.675
0.374
0.195
0.959
0.589
0.341
0.167
100
0.558
0.365
0.152
0.051
0.519
0.254
0.132
0.073
200
0.172
0.134
0.076
0.004
0.120
0.141
0.081
0.005
2
50
1.258
0.828
0.449
0.240
1.152
0.692
0.390
0.154
100
0.569
0.448
0.208
0.099
0.516
0.427
0.211
0.022
200
0.232
0.181
0.122
0.009
0.214
0.145
0.120
0.008Tables1 and 2 show that in general FSTN distribution enjoys more robustness and flexibility in modeling data with skewness and heavy tails as well as multimodality in comparison with other skewed alternatives and the implementation of MPNR method brings more accuracy and improvement for model estimation in the context of nonlinear regression with this new distribution.


## 6. Conclusion

We have proposed a new skewed distribution based on the general FSS distribution framework, called the FSTN distribution, which is allowed to accommodate multimodality, asymmetry, and heavy tails jointly to offer greater flexibility than SN and STN counterparts. Moreover, we have developed a modified profile of Newton-Raphson iterative algorithm for estimating the parameters of interest of nonlinear model with FSTN distribution and the interval estimation and hypothesis test in a profile likelihood paradigm are also considered.Numerical studies reveal that, in the context of nonlinear regression analysis, if the true distribution is STN or FSN whereas the fitting distribution is FSTN, the estimation results are not influenced by this misspecification of distribution assumption. However, once the true distribution is FSTN while the fitting distribution is STN, the estimation results appear to be somewhat disappointing, which shows the robustness of FSTN distribution. In general, the combination of FSTN distribution with MPNR method brings more accuracy and improvement on the estimation of nonlinear regression.So far the present methodology is limited to the complete data analysis; the extensions of this paper include missing data version as well as Bayesian analysis of this model, which will be reported in another paper.



---

*Source: 542985-2014-12-03.xml*