Contents lists available atScienceDirect
Computational Statistics and Data Analysis
journal homepage:www.elsevier.com/locate/csda
A goodness-of-fit test for Archimedean copula models in the presence of
right censoring
Takeshi Emura, Chien-Wei Lin, Weijing Wang
∗Institute of Statistics, National Chiao-Tung University, Hsin-Chu, Taiwan, ROC
a r t i c l e i n f o
Article history:
Received 28 November 2008
Received in revised form 15 March 2010 Accepted 16 March 2010
Available online 2 April 2010 Keywords:
Archimedean copulas Bivariate failure times Clayton model Concordance estimator Goodness of fit test U-statistics
a b s t r a c t
A goodness-of-fit testing procedure for Archimedean copula (AC) models is developed based on right-censored data. The proposed approach extends an existing method, which is suitable for the Clayton model, to general AC models. Asymptotic properties of the proposed test statistics under the true model assumption are derived. Simulation analysis shows that the proposed test has reasonable performance. Finally, two real data examples are analyzed for illustrative purposes.
© 2010 Elsevier B.V. All rights reserved.
1. Introduction
Copula models are popular choices for describing multivariate data. They have the major advantage that the dependence structure can be studied separately from marginal distributions. Archimedean copula (AC) models form a useful family of copula models in which the dependence structure can be further characterized by a univariate function (Nelsen, 2006, Section 4). See Frees and Valdez (1998) for insurance applications and Wang and Wells(2000) for biomedical applications. In this article, we consider the bivariate case and propose a testing procedure for checking the goodness-of-fit for the imposed AC model assumption.
If the failure time variables
(
X,
Y)
follow a copula model, their joint survival distribution S(
x,
y) =
Pr(
X>
x,
Y>
y)
can be written asS
(
x,
y) =
C{
SX(
x),
SY(
y)},
where SX
(
x) =
S(
x,
0)
and SY(
y) =
S(
0,
y)
are the marginal survival functions, and C(
u, v) : [
0,
1]
2→ [
0,
1]
is the copula function. For AC models, the bivariate copula function can be further simplified asCα
(
u, v) = φ
−α1[
φ
α(
u) + φ
α(v)],
(1)where
φ
α(·) : [
0,
1] → [
0, ∞]
is a univariate function which has two continuous derivatives satisfyingφ
α(
1) =
0, φ
0α
(
t) =
∂φα(t)
∂t
<
0 andφ
00α(
t) =
∂ 2φα(t)
∂t2
>
0. The parameterα
measures the degree of association related to Kendall’sτ
, which is a rank-invariant correlation measure. Specifically, let(
Xi,
Yi)
and(
Xj,
Yj)
be two random replications of(
X,
Y)
. Kendal’s tau is defined asτ =
Pr(
∆ij=
1) −
Pr(
∆ij=
0),
∗Corresponding author. Tel.: +886 3 5728746; fax: +886 2 5728745.
E-mail address:wjwang@stat.nctu.edu.tw(W. Wang).
0167-9473/$ – see front matter©2010 Elsevier B.V. All rights reserved.
where∆ij
=
I{
(
Xi−
Xj)(
Yi−
Yj) >
0}
. It holds thatτ =
4Z
∞ 0Z
∞ 0 S(
x,
y)
dS(
x,
y) −
1=
4Z
1 0Z
1 0 C(
u, v)
dC(
u, v) −
1=
4Z
1 0φ
α(
t)
φ
0 α(
t)
dt+
1,
where the last identity specifies the relationship between
τ
andα
for an AC model. The following cross ratio function (Oakes, 1989) is useful for describing the local association:θ
∗(
x
,
y) =
S(
x,
y) · ∂
2S
(
x,
y)/∂
x∂
y∂
S(
x,
y)/∂
x·
∂
S(
x,
y)/∂
y.
For AC models,Oakes(1989) showed that
θ
∗(
x,
y)
can be written asθ
∗(
x
,
y) = θ
α{
S(
x,
y)},
(2)where
θ
α(v) = −vφ
00(v)/φ
0(v)
is a univariate function.Semiparametric inference of the association parameter
α
for a chosen AC model has been a popular research topic in the literature. The form ofφ
α(·)
is assumed but the marginal distributions are not specified. Therefore how to choose an appropriateφ
α(·)
based on the data at hand is an important issue for practical applications.Oakes(1989) andGenest and Rivest(1993) derived analytic properties of AC models that are useful for model selection. Specifically,Oakes(1989) showed thatE
[
∆ij| ˜
Xij=
x, ˜
Yij=
y] =
θ
α{
S(
x,
y)}
θ
α{
S(
x,
y)} +
1,
(3)whereX
˜
ij=
Xi∧
XjandY˜
ij=
Yi∧
Yj, where a∧
b=
min(
a,
b)
.Genest and Rivest(1993) derived that the variable V=
S(
X,
Y)
has the distribution functionKα
(v) =
Pr(
V≤
v) = v − λ
α(v),
(4)where
λ
α(v) = φ
α(v)/φ
0α
(v)
for 0< v ≤
1.Eq.(4)has been adopted for model selection byGenest and Rivest(1993) for complete data and byWang and Wells (2000) for bivariate censored data. The idea is to measure the goodness-of-fit based on the distance between nonparametric and model-based estimators of K
(v) =
Pr(
V≤
v)
which characterizes the model. However, the distribution of the distance statistics (i.e., the supremum norm or L2 form of a function) is quite complicated. Alternatively, some authors suggested approximating the distribution of a test statistics under the null hypothesis based on re-sampling techniques such as parametric bootstrap procedures (Dobric and Schmid, 2007;Nikoloulopoulos and Karlis, 2008;Genest et al.,2009). However, the computational cost of re-sampling methods may restrict their applicability.A different strategy, originally proposed byGill and Schumacher(1987), is to compare the discrepancy between two point estimators which are derived under the same class of estimating functions but with different weight functions. If the imposed model assumption is correctly specified, the two estimators will be close to each other. On the other hand, when the model assumption is wrong, the values of the two estimators will fall apart if the weights are properly chosen. Formal goodness-of-fit tests can be constructed since asymptotic properties are easier to handle.Shih(1998) adopted this idea to check the assumption of the Clayton model with
φ
α(
t) = (
t−α−
1)/α
. In this article we extend this approach to general AC models based on Eq.(3). The proposed goodness-of-fit statistics are introduced in Section2. Section3contains simulation results and Section4contains real data applications. Concluding remarks are given in Section5.2. The proposed test We consider testing
H0
:
C(
u, v) = φ
α−1[
φ
α(
u) + φ
α(v)]
for someα ∈
R.
The alternative hypothesis can be any copula other than the one specified as H0. We will examine the power performance when the alternative hypothesis follows a different AC model or a non-AC model. We will first ignore censoring to illustrate the main idea. In Section2.5, the proposed method will be modified to account for the censoring effect.
2.1. The main idea
In the absence of censoring, the observed data can be written as
{
(
Xi,
Yi); (
i=
1, . . . ,
n)}
, and the concordance indicator ∆ijis observable. A class of estimating equations based on Eq.(3)can be constructed asUk
(α) =
X
i<j Wk( ˜
Xij, ˜
Yij, α)
"
∆ij−
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1#
,
(5)where Wk
(·, ·, ·)
is a weight function andSˆ
(
x,
y)
is an estimator of S(
x,
y)
. For complete data, we can use the empirical sur-vival functionSˆ
(
x,
y) =
1nP
iI
(
Xi≥
x,
Yi≥
y)
. To construct a goodness-of-fit test, we need to choose two weight functions, and then obtainα
ˆ
kwhich solves Uk(α) =
0 for k=
1,
2. Under H0and some regularity conditions, n1/2{ ˆ
α
1− ˆ
α
2}
converges to a mean-zero normal distribution. The power of the test depends on the choice of two weights. In principle, one weight func-tion should utilize more model informafunc-tion while the other should be less model-dependent. It is expected that n1/2{ ˆ
α
1
− ˆ
α
2}
will show a large discrepancy when H0does not hold. We propose to choose one weight function using some likelihood in-formation, which will be discussed in Section2.2. For the other weight function, we suggest the unweighted version since it does not contain any model information. Sometimes, for better normal approximation, one may use n1/2{ ˆ
γ
1
− ˆ
γ
2}
, whereγ
is a monotone function ofα
.2.2. Conditional likelihood and the weight function
We generalize the likelihood approach proposed byClayton(1978) and show that the resulting estimating function can also be written as the form of Uk
(α)
in(5). Define the set of grid points,ψ =
(
(
x,
y) |
nX
i=1 I(
Xi=
x,
Yi≥
y) =
1,
nX
i=1 I(
Xi≥
x,
Yi=
y) =
1)
.
Define D
(
x,
y) = P
ni=1I(
Xi=
x,
Yi=
y)
, which measures the number of observed failures at(
x,
y) ∈ ψ
, and R(
x,
y) =
P
ni=1I
(
Xi≥
x,
Yi≥
y)
, which measures the number at risk at(
x,
y) ∈ ψ
. Assuming that the data contains no ties and conditional on R(
x,
y)
for(
x,
y) ∈ ψ
, D(
x,
y)
is a Bernoulli random variable with the success probability Pr{
D(
x,
y) =
1|
R
(
x,
y) =
r}
.Clayton(1978) derived this probability under the Clayton model withφ
α(v) = (v
−α−
1)/α
. For an AC model, we show thatPr
{
D(
x,
y) =
1|
R(
x,
y) =
r} =
θ
α{
S(
x,
y)}
r
−
1+
θ
α{
S(
x,
y)}
.
Combining all the points in
ψ
under the working assumption of independence among different grids, the likelihood function can be written as L(α) =
Y
(x,y)∈ψθ
α{
S(
x,
y)}
R(
x,
y) −
1+
θ
α{
S(
x,
y)}
D(x,y) R(
x,
y) −
1 R(
x,
y) −
1+
θ
α{
S(
x,
y)}
1−D(x,y).
(6)Motivated by the paper ofOakes(1986), we apply some algebraic operations (given inAppendix A) to show that the resulting log-likelihood can be written as
l
(α) =
X
i logθ
α{
S(
Xi,
Yi)}
Rii−
1+
θ
α{
S(
Xi,
Yi)}
+
X
i<j(
1−
∆ij)
log"
Rij−
1 Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
#
,
(7)where Rij
=
R( ˜
Xij, ˜
Yij) =
nSˆ
( ˜
Xij, ˜
Yij)
. The resulting score function becomesU1
(α) =
X
i<j˙
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1]
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[
Rij−
1+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}]
"
∆ij−
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1#
=
X
i<j W1( ˜
Xij, ˜
Yij)
"
∆ij−
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1#
,
(8a)where
θ
˙
α(v) = ∂θ
α(v)/∂α
. The second estimating function is given byU2
(α) =
X
i<j"
∆ij−
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1#
.
(8b)By solving Uk
(α) =
0, we haveα
ˆ
kfor k=
1,
2.Though we have assumed the absence of ties in the derivation, the resulting Eqs.(8a)and(8b)can still be calculated for tied data without any tie breaking method. Note that the validity of(8a)and(8b)is due to the moment equation(3)that is derived under the continuity assumption of S
(
x,
y)
. If the underlying distribution is continuous but ties occur due to random round-off of the observation, the resulting bias in(8a)and(8b)should be modest. On the other hand, if ties occur due to the discrete mass in the distribution, the proposed method may not be appropriate.For the Clayton model with
φ
α(
t) = (
t−α−
1)/α
,θ
α{
S(
x,
y)} = α +
1, which is a constant. In this special case, Eqs.(8a)and(8b)reduce to the results inShih(1998). For the Gumbel model withφ
α(v) = {−
log(v)}
α+1, we haveθ
α{
S(
x,
y)} =
1−
α/{
log S(
x,
y)}
and W1( ˜
Xij, ˜
Yij) =
2 logSˆ
( ˜
Xij, ˜
Yij) − α
{
logSˆ
( ˜
Xij, ˜
Yij) − α}{α −
RijlogSˆ
( ˜
Xij, ˜
Yij)}
.
Therefore, U1(α) =
X
i<j 2 logSˆ
( ˜
Xij, ˜
Yij) − α
{
logSˆ
( ˜
Xij, ˜
Yij) − α}{α −
RijlogSˆ
( ˜
Xij, ˜
Yij)}
"
∆ij−
logSˆ
( ˜
Xij, ˜
Yij) − α
2 logSˆ
( ˜
Xij, ˜
Yij) − α
#
.
(9)We will use the Gumbel model for illustration in subsequent discussions.
2.3. Asymptotic distributional theory
In this section, we derive the asymptotic distribution of n1/2
( ˆα
1− ˆ
α
2)
and also n1/2( ˆγ
1− ˆ
γ
2)
if the latter provides a better approximation to the normal distribution. The theoretical challenge comes from the fact that U1(α)
and U2(α)
involve the estimatorSˆ
(
x,
y)
(see Eq.(9)for the Gumbel case).Lemma 1states that U1(α)
and U2(α)
can be approximated by the following two U-statistics:˜
U1(α) =
X
i<j˙
θ
α{
S( ˜
Xij, ˜
Yij)}[θ
α{
S( ˜
Xij, ˜
Yij)} +
1]
θ
α{
S( ˜
Xij, ˜
Yij)}
S( ˜
Xij, ˜
Yij)
"
∆ij−
θ
α{
S( ˜
Xij, ˜
Yij)}
θ
α{
S( ˜
Xij, ˜
Yij)} +
1#
,
(10a) and˜
U2(α) =
X
i<j"
∆ij−
θ
α{
S( ˜
Xij, ˜
Yij)}
θ
α{
S( ˜
Xij, ˜
Yij)} +
1#
.
(10b)Lemma 1. Under the correct model and regularity conditions inAppendix B.1,
n
n 2 −1 U1(α) =
n 2 −1˜
U1(α) +
oP(
1),
n 2 −1 U2(α) =
n 2 −1˜
U2(α) +
oP(
1),
where oP(
1)
is uniform inα
.If the parameter
α
has a positive value as in the Gumbel model, the natural logarithm transformation can improve the normal approximation. The following two results are derived by takingγ =
logα
andγ
ˆ
k=
logα
ˆ
kfor k=
1,
2.Lemma 2. Under the correct model and the regularity conditions inAppendix B.1,
n1/2
( ˆγ
1− ˆ
γ
2) =
n1/2 n 2 −1X
i<j h{
(
Xi,
Yi), (
Xj,
Yj)} +
oP(
1),
where the function h is symmetric in its arguments and is defined as
h
{
(
Xi,
Yi), (
Xj,
Yj)} ≡
1α
˙θ
α{
S( ˜
Xij, ˜
Yij)}[θ
α{
S( ˜
Xij, ˜
Yij)} +
1]
ALθ
α{
S( ˜
Xij, ˜
Yij)}
S( ˜
Xij, ˜
Yij)
−
1 A! "
∆ij−
θ
α{
S( ˜
Xij, ˜
Yij)}
θ
α{
S( ˜
Xij, ˜
Yij)} +
1#
,
A≡
E˙
θ
α{
S( ˜
X12, ˜
Y12)}
[
θ
α{
S( ˜
X12, ˜
Y12)} +
1]
2!
and AL≡
E[ ˙
θ
α{
S( ˜
X12, ˜
Y12)}]
2θ
α{
S( ˜
X12, ˜
Y12)}[θ
α{
S( ˜
X12, ˜
Y12)} +
1]
!
.
Based onLemma 2and applying the central limit theorem for U-statistics (p. 162, Theorem 12.3 ofVan Der Vaart, 1998), we can obtain the following theorem.
Theorem 1. Under the correct model and the regularity conditions inAppendix B.1, n1/2
( ˆγ
1
− ˆ
γ
2)
converges in distribution to anormal distribution with mean zero and variance
σ
2=
4E[
h{
(
X2.4. The proposed testing procedure
For testing H0, we propose to reject H0at the significance level
κ
if| ˆ
γ
1− ˆ
γ
2|
/ ˆσ
is greater than the upperκ/
2 percentile of the standard normal distribution, whereσ
ˆ
is an estimator ofσ/
√
n.The asymptotic variance
σ
2can be estimated based on its analytical expression by averaging over all possible observation triples. However this analytic formula becomes complicated for right-censored data. Alternatively, we can use the jackknife estimator defined asˆ
σ
2 Jack=
n−
1 n nX
i=1{ ˆ
γ
1(−i)− ˆ
γ
2(−i)−
( ˆγ
1(·)− ˆ
γ
2(·))}
2,
where
γ
ˆ
k(−i)is the statisticγ
ˆ
kignoring the i-th observation andγ
ˆ
1(·)− ˆ
γ
2(·)=
1n
P
n i=1( ˆγ
(−i) 1− ˆ
γ
(−i) 2)
.2.5. Modification for right censoring
When
(
Xi,
Yi)
are subject to right censoring by(
Ai,
Bi)
such that one can only observe(
^ Xi,
^ Yi, δ
iX, δ
iY) (
i=
1, . . . ,
n)
, whereX^i=
Xi∧
Ai,
^Yi
=
Yi∧
Bi, δ
iX=
I(
Xi≤
Ai)
andδ
iY=
I(
Yi≤
Bi)
. We assume that(
Ai,
Bi)
is independent of(
Xi,
Yi)
. The order of Xiand Xjis known if and only ifX˜
ij≤ ˜
Aij, whereX˜
ij=
Xi∧
XjandA˜
ij=
Ai∧
Aj. Similarly, the order of Yiand Yj can be known if and only ifY˜
ij≤ ˜
Bij. Define Zij=
I( ˜
Xij≤ ˜
Aij, ˜
Yij≤ ˜
Bij)
, which indicates whether the ordering relationship is certain or not. We can modify Uk(α)
in(5)by selecting only ‘‘orderable’’ pairs with Zij=
1 such thatUk
(α) =
X
i<j ZijWk( ˜
Xij, ˜
Yij, α, ˆ
S)
"
∆ij−
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} +
1#
(
k=
1,
2),
(11a) where∆ij=
I{
(
^ Xi−
^ Xj)(
^ Yi−
^Yj
) >
0}
andSˆ
(
x,
y)
is an estimator of S(
x,
y)
suitable for right-censored data. Note that, to simplify the presentation, we may use the same notations but change their definitions for the censoring case. The first suggested weight function has the formW1
( ˜
Xij, ˜
Yij, α,
S) =
˙
θ
α{
S( ˜
Xij, ˜
Yij)}[θ
α{
S( ˜
Xij, ˜
Yij)} +
1]
θ
α{
S( ˜
Xij, ˜
Yij)}[
Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}]
,
where Rij=
P
n l=1I(
^ Xl≥ ˜
Xij,
^Yl
≥ ˜
Yij)
and the second one is W2( ˜
Xij, ˜
Yij, α,
S) =
1.A possible candidate ofS
ˆ
(
x,
y)
is the nonparametric estimator proposed byDabrowska(1988). Since the estimator is somewhat complicated, one may use the following model-based estimator:˜
Sα
(
x,
y) = φ
α−1{
φ
α[ˆ
SX(
x)] + φ
α[ˆ
SY(
y)]},
whereS
ˆ
X(
x)
andSˆ
Y(
y)
are marginal Kaplan–Meier estimators of SX(
x)
and SY(
y)
, respectively. To obtainα
ˆ
k, one can solve the equationX
i<j ZijWk( ˜
Xij, ˜
Yij, α, ˜
S(l−1))
"
∆ij−
θ
α{˜
S(l−1)( ˜
Xij, ˜
Yij)}
θ
α{˜
S(l−1)( ˜
X ij, ˜
Yij)} +
1#
=
0 (11b)successively for l
=
1,
2, . . .
, until some convergence criterion is met, whereS˜
(l)is defined asS˜
α(
x,
y)
, withα
being the estimated value in the l-step. The initial valueα
0may be obtained by inverting a naïve estimate of Kendall’s tau based on{
(
X^i,
^Yi
) : (
i=
1, . . . ,
n)}
. Our simulation analysis shows that the two methods for handling the survival function yield close results.The asymptotic normality results can be established using similar arguments as in the uncensored case. The key step is to approximate Uk
(α)
in(11a)by a U-statistic. For instance,˜
U2(α) =
X
i<j Zij"
∆ij−
θ
α{
S( ˜
Xij, ˜
Yij)}
θ
α{
S( ˜
Xij, ˜
Yij)} +
1#
,
is a U-statistic that can approximate U2
(α)
. Under regularity conditions such as those listed in Proposition 4.1 ofDabrowska (1988),Sˆ
(
x,
y)
(orS˜
α(
x,
y)
) converges in probability to S(
x,
y)
uniformly in(
x,
y) ∈ [
0, τ
1) × [
0, τ
2)
for some(τ
1, τ
2)
. Then, it follows from the same arguments as the proof ofLemma 1thatn
n 2 −1 U2(α) =
n 2 −1˜
U2(α) +
oP(
1).
Applying a Taylor series expansion and the central limit theorem for U-statistics, as in the proof ofLemma 1, n1/2
(
logα
ˆ
1−
Table 1A
Empirical probabilities of not rejecting the Gumbel model at the 5% significance level.
τ =0.3 τ =0.4 τ =0.5 τ =0.6 τ =0.7 n=100 CEN%=(0, 0) 0.99 0.97 0.96 0.96 0.95 CEN%=(20, 20) 0.98 0.93 0.95 0.97 0.96 CEN%=(50, 50) 0.97 0.96 0.95 0.96 0.99 n=200 CEN%=(0, 0) 0.97 0.95 0.92 0.93 0.93 CEN%=(20, 20) 0.95 0.94 0.94 0.93 0.96 CEN%=(50, 50) 0.95 0.97 0.98 0.97 0.96
Note: CEN% denotes the two marginal censoring rates, 100×Pr(A<X)and 100×Pr(B<Y). The probabilities are calculated based on 100 replications.
Table 1B
Means and standard deviations (in parentheses) of the test statistic| ˆγ1− ˆγ2|/ ˆσJackunder H0.
τ =0.3 τ =0.4 τ =0.5 τ =0.6 τ =0.7 n=100 CEN%=(0, 0) 0.06 (0.89) 0.18 (0.90) 0.28 (1.00) 0.12 (0.89) 0.32 (0.91) CEN%=(20, 20) −0.01 (0.94) 0.30 (1.00) 0.13 (1.04) 0.19 (1.03) 0.14 (0.89) CEN%=(50, 50) 0.15 (0.99) −0.13 (0.94) 0.14 (0.98) 0.17 (0.81) 0.06 (0.81) n=200 CEN%=(0, 0) 0.07 (0.94) 0.13 (1.01) 0.22 (0.94) 0.05 (0.94) 0.06 (0.94) CEN%=(20, 20) 0.12 (0.99) 0.19 (1.02) 0.10 (1.02) 0.05 (0.89) 0.07 (0.97) CEN%=(50, 50) −0.01 (0.83) 0.03 (0.92) 0.12 (0.84) 0.13 (0.97) 0.04 (0.90) Note: CEN% denotes the two marginal censoring rates, 100×Pr(A<X)and 100×Pr(B<Y). The means and standard deviations are calculated based on 100 replications.
3. Simulation results
The simulation study contains three components. We first examine the performance of the proposed test when H0is correctly specified. Then, we study the robustness of the proposed test under dependent censoring, a condition which violates the assumption of independent censorship. Finally, we study the power of the test when H0does not hold. The level of significance is set to be 5%.
Here we use the Gumbel model for illustration. The null hypothesis is given by
H0
:
C(
u, v) =
exp−
(−
log u)
α+1+
(−
logv)
α+1 1 α+1 for someα >
0.
The two marginal distributions both follow an exponential distribution with the hazard rate equal to 1. Specifically, we let
X
= −
U1/(α+1)log(
V),
Y= −
(
1−
U)
1/(α+1)log(
V),
where U
∼
Uniform(
0,
1)
and V has the distribution function Kα(v) = v − v
log(v)/(α +
1)
. Five values ofα
whose correspondingτ
is equal to 0.3, 0.4, 0.5, 0.6 and 0.7 are chosen. The bivariate censoring variables(
A,
B)
are generated independently from exponential distributions such that the censoring proportion is 0, 0.2 and 0.5 respectively in each coordinate. Based on the data, we obtainα
ˆ
kwhich solves Eq.(11b), and thenγ
ˆ
k=
logα
ˆ
kfor k=
1,
2. The jackknife estimatorσ
ˆ
2Jack, illustrated in Section2.4, is used for variance estimation. The null hypothesis is rejected if
| ˆ
γ
1− ˆ
γ
2|
/ ˆσ
Jack exceeds 1.96.Table 1Areports the empirical probabilities of not rejecting H0based on 100 replications. When n
=
100, the type-Ierror rates are slightly lower than the nominal 5% level in some cases but the results improve as the sample size increases to n
=
200.Table 1Breports the means and standard deviations of| ˆ
γ
1− ˆ
γ
2|
/ ˆσ
Jackwhich can be used to evaluate the validity of normal approximation. We see that, when n increases, the approximation also improves. In general, the means are close to zero and the standard deviations are slightly less than one. This may be due to the fact that the jackknife algorithm tends to overestimate the variance, which results in slightly lower type-I error rate.We also evaluate the robustness of the proposed test under dependent censoring in which
(
A,
B)
and(
X,
Y)
are no longer independent. Let A=
C1I(
X< −
log(
0.
5)) +
C2I(
X≥ −
log(
0.
5))
, where C1and C2are independent exponential distributions with hazard ratesθ/
2 and 1/
2, respectively. It is easy to see that X and A are positively correlated when 1< θ
and negatively correlated when 0< θ <
1. We define B in the same way.Table 2Asummarizes the empirical probabilities of not rejecting H0based on 100 replications under the dependent
censoring setting. Surprisingly, the empirical type-I error rates are still close to the nominal 5% level. In fact,Table 2Bshows that
| ˆ
γ
1− ˆ
γ
2|
/ ˆσ
Jackcan still be approximated by the standard normal distribution. However, we find that bothγ
ˆ
1andγ
ˆ
2Table 2A
Empirical probabilities of not rejecting the Gumbel model at the 5% significance level under dependent censoring.
τ =0.3 τ =0.4 τ =0.5 τ =0.6 τ =0.7 n=100 θ =1/4 CEN%=(28, 28) 0.98 0.96 0.96 0.95 0.96 θ =1/2 CEN%=(30, 30) 0.98 0.97 0.95 0.96 0.95 θ =2 CEN%=(39, 39) 0.97 0.98 0.97 0.99 0.98 θ =4 CEN%=(47, 47) 0.99 0.96 0.94 0.99 0.97 n=200 θ =1/4 CEN%=(28, 28) 0.93 0.95 0.93 0.93 0.94 θ =1/2 CEN%=(30, 30) 0.95 0.95 0.96 0.95 0.97 θ =2 CEN%=(39, 39) 0.91 0.94 0.93 0.97 0.97 θ =4 CEN%=(47, 47) 0.97 0.97 0.97 0.97 0.96
Note:θis related to the degree of dependence between failure times and censoring times. CEN% denotes the two marginal censoring rates, 100×Pr(A<X) and 100×Pr(B<Y). The probabilities are calculated based on 100 replications.
Table 2B
Means and standard deviations (in parentheses) of the test statistics| ˆγ1− ˆγ2|/ ˆσJackunder H0and dependent censoring.
τ =0.3 τ =0.4 τ =0.5 τ =0.6 τ =0.7 n=100 θ =1/4 CEN%=(28, 28) 0.11 (0.84) 0.01 (0.96) 0.16 (0.91) 0.17 (0.95) 0.17 (0.90) θ =1/2 CEN%=(30, 30) 0.21 (0.80) −0.01 (0.98) 0.23 (0.90) 0.12 (0.85) 0.19 (0.99) θ =2 CEN%=(39, 39) 0.08 (0.99) 0.13 (0.77) 0.06 (0.84) −0.10 (0.82) −0.01 (0.82) θ =4 CEN%=(47, 47) 0.08 (0.85) 0.01 (1.03) −0.31 (0.98) −0.14 (0.79) −0.06 (0.86) n=200 θ =1/4 CEN%=(28, 28) 0.06 (1.15) 0.12 (0.99) 0.26 (1.05) 0.28 (1.01) 0.27 (0.95) θ =1/2 CEN%=(30, 30) 0.06 (0.98) 0.15 (0.94) 0.20 (0.97) −0.05 (0.91) 0.17 (0.81) θ =2 CEN%=(39, 39) 0.13 (1.01) 0.12 (1.00) 0.14 (1.06) −0.04 (0.92) −0.10 (0.95) θ =4 CEN%=(47, 47) 0.03 (0.90) 0.01 (0.94) −0.01 (0.93) −0.28 (0.85) −0.24 (0.87) Note:θis related to the degree of dependence between failure times and censoring times.τdenotes the Kendall’s tau between X and Y . CEN% denotes the two marginal censoring rates, 100×Pr(A<X)and 100×Pr(B<Y). The means and standard deviations are calculated based on 100 replications.
censoring affects the biasness of the two estimators in the same direction and with similar magnitude. We also investigated different setups of dependent censoring, and the results were similar. More thorough analysis for the robustness under dependent censoring deserves future study.
To assess the power of the test, data is used from the Clayton and Frank models, respectively. The power is the probability of correctly rejecting Gumbel’s assumption. The empirical power functions are plotted inFigs. 1Afor n
=
100 andFig. 1B for n=
200, respectively. In each figure, we see that the power deceases asτ
decreases. This is reasonable since these three models coincide with each other as Kendall’sτ
approaches zero, and it becomes more difficult to distinguish between sim-ilar models. As the level of censoring increases, the power performance becomes worse. As expected, the powers become larger when the sample size increases.We further examine the power when the alternative model is similar to that specified in the null hypothesis. Here we consider the Galambos copula model, which is not an AC model but is similar to the Gumbel AC model in some key aspects. For example, both copulas have upper tail dependences. In terms of the Jeffreys divergence measure, the distance between the two models is only 0.001 even when
τ ≥
0.
5 (Nikoloulopoulos and Karlis, 2008). We report the results underτ =
0.
5. When the censoring proportion is 50%, the power is only 4% for n=
100 and 11% for n=
200. Nevertheless, the power in-creases to 23% for n=
100 and 40% for n=
200 when the censoring proportion reduces to 20%. In the absence of censoring, the power further increases to 48% for n=
100 and 79% for n=
200.It is worth noting that the method in Shih(1998) was developed to verify the Clayton model assumption. In her simulations, the Clayton assumption was set as the null hypothesis and the power under the Gumbel alternative was evaluated. Comparing Shih’s analysis with the left panels ofFigs. 1Aand1B, the two results are close, which implies that reversing the roles of the null and alternative hypotheses does not matter.
4. Real data applications
The first example is from the Australian Twin Study (Duffy et al., 1990), in which
(
X,
Y)
represents the ages at appendicectomy measured for each twin pair. As inPrentice and Hsu(1997), a subset of the original data containing 748 dizygotic pairs is chosen. In the sample, 82 observations are uncensored, 117 are censored for X , 105 are censored for Y and 444 are censored for both X and Y . The proposed method is applied to test the goodness-of-fit for four AC model candidates individually, and the results are summarized inTable 3A. The Gumbel model provides the best fit to the data while the Clayton and Log-copula models are rejected at the 5% significance level.Kendall’s tau Clayton Model Kendall’s tau Frank Model 0% censored 20% censored 50% censored 0.3 0.4 0.5 0.6 0.7 Power (%) 100 80 60 40 20 0 Power (%) 100 80 60 40 20 0 0% censored 20% censored 50% censored 0.3 0.4 0.5 0.6 0.7
Fig. 1A. Empirical powers with n = 100 under H0: Gumbel versus Ha: Not Gumbel. Powers are the rates of rejecting H0with 5% significance during 100 replications. Kendall’s tau Clayton Model 0% censored 20% censored 50% censored Frank Model 0% censored 20% censored 50% censored 0.3 0.4 0.5 0.6 0.7 Kendall’s tau 0.3 0.4 0.5 0.6 0.7 Power (%) 100 80 60 40 20 0 Power (%) 100 80 60 40 20 0
Fig. 1B. Empirical powers with n =200 under H0: Gumbel versus Ha: Not Gumbel. Powers are the rates of rejecting H0with 5% significance during 100 replications.
Table 3A
The goodness-of-fit test results for four AC models based on the Australian Twin Study (Duffy et al., 1990).
ˆ α1(ˆτ1) αˆ2(ˆτ2) ( ˆγ1− ˆγ2)/ ˆσJack p-value Clayton 1.446 (0.182) 1.717 (0.264) −1.867 0.000 Frank 1.308 (0.143) 1.496 (0.163) −1.090 0.117 Gumbel 0.115 (0.103) 0.114 (0.102) 0.084 0.497 Log-copula 1.447 (0.295) 1.147 (0.248) 1.351 0.034
In the second example, X and Y represent the time to the first and second recurrence of infection in kidney patients on dialysis (McGilchrist and Aisbett, 1991). Among 38 patients, 6 observations are censored for X , 12 observations are censored for Y and 3 observations are censored for both X and Y . This data has been analyzed byWang and Wells(2000), in which Eq.(4)was used to construct a distance measure between nonparametric and model-based curves. For the kidney data, their approach selected the Gumbel model which gives the smallest distance among several alternatives. Note that their method is not a formal statistical test. Applying the proposed test to four model candidates, the p-values are 0.452 (Gumbel), 0.189 (Clayton), 0.365 (Frank) and 0.423 (Log-copula). Due to the small sample size, no model is rejected. However, the Gumbel model still gives the largest p-value, which supports the conclusion ofWang and Wells(2000). Under the Gumbel model, we obtain
α
ˆ
1=
0.
282(ˆτ
1=
0.
220), ˆα
2=
0.
262(ˆτ
2=
0.
208)
.5. Conclusion
Copula models are often defined without specifying the marginal distributions. Likelihood-based strategies for model selection, such as the Akaike Information Criterion (AIC), cannot be directly applied to such a semiparametric structure. The proposed test extends the approach used byShih(1998) from the Clayton model to a class of AC models. The property
of the conditional moment in(3)is utilized to construct a class of estimating functions for the association parameter. One suggested weight function is derived based on the conditional likelihood function, while the other is the unweighted version. The contrast between the two weights provides a way to detect any violation of the model assumption. Compared with existing methods for model selection for copulas, the proposal yields a formal goodness-of-fit test.
In the presence of right censoring, we propose to delete non-orderable pairs from the estimating equation since the remaining observations still maintain the unbiased property. This strategy has been used byOakes (1982,1986)andShih (1998). Despite its simplicity, data deletion may result in efficiency loss especially when censoring is heavy. The techniques of handling missing data, including imputation and weighting, may be applied to handle the situation that∆ijis unknown due to censoring (Hsieh, 2010). This direction deserves further investigation. However, if alternative approaches involve extra estimation for nuisance parameters, efficiency gain is not guaranteed.
Acknowledgements
The authors are grateful to the editor and the two referees for their helpful comments on the first version of this paper. We also thank Ms. Sophia Shing for proofreading the article. The research was financially supported by the research Grant 098-2811-M-009-025 funded by National Science Council of Taiwan, ROC.
Appendix A. Derivation of the log-likelihood function Taking the log of Eq.(6), we obtain the log-likelihood function
l
(α) =
X
(x,y)∈ϕ D(
x,
y)
logθ
α{
S(
x,
y)}
R(
x,
y) −
1+
θ
α{
S(
x,
y)}
+ {
1−
D(
x,
y)}
log R(
x,
y) −
1 R(
x,
y) −
1+
θ
α{
S(
x,
y)}
=
X
iX
j:Xj≥Xi,Yj≤Yi D(
Xi,
Yj)
logθ
α{
S(
Xi,
Yj)}
R(
Xi,
Yj) −
1+
θ
α{
S(
Xi,
Yj)}
+ {
1−
D(
Xi,
Yj)}
log R(
Xi,
Yj) −
1 R(
Xi,
Yj) −
1+
θ
α{
S(
Xi,
Yj)}
=
X
i logθ
α{
S(
Xi,
Yi)}
R(
Xi,
Yi) −
1+
θ
α{
S(
Xi,
Yi)}
+
X
iX
j:Xj>Xi,Yj<Yi log R(
Xi,
Yj) −
1 R(
Xi,
Yj) −
1+
θ
α{
S(
Xi,
Yj)}
.
As long as Xj
>
Xiand Yj<
Yihold, Xi= ˜
Xijand Yj= ˜
Yij, and hence S( ˜
Xij, ˜
Yij) =
S(
Xi,
Yj)
and R(
Xi,
Yj) =
Rij. Thus, the second term in the previous equation becomesL
≡
X
iX
j:Xj>Xi,Yj<Yi log"
Rij−
1 Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
#
=
X
iX
j:Xj>Xi I(
Yj<
Yi)
log"
Rij−
1 Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
#
.
It follows that L=
X
iX
j:Xj>Xi(
1−
∆ij)
log"
Rij−
1 Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
#
=
X
i<j(
1−
∆ij)
log"
Rij−
1 Rij−
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
#
.
Appendix B B.1. Regularity conditions(i) A parameter spaceΘfor
α
is open and the true value lies insideΘ.(ii)
θ
α(v)
is twice differentiable with respect toα
and differentiable with respect tov
. (iii) 1/θ
α(v)
andθ
0α
(v) =
∂v∂θ
α(v)
are bounded with respect to parameters(α, v)
. (iv) A and ALexist and are strictly positive.B.2. Proof ofLemma 1
First, we show that
n 2 −1{
U2(α) − ˜
U2(α)}
converges in probability to zero uniformly inα
. It is straightforward to see that|
U2(α) − ˜
U2(α)| =
X
i<jθ
α{ˆ
S( ˜
Xij, ˜
Yij)}
1+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}
−
θ
α{
S( ˜
Xij, ˜
Yij)}
1+
θ
α{
S( ˜
Xij, ˜
Yij)}
≤
X
i<jθ
α{ˆ
S( ˜
Xij, ˜
Yij)} − θ
α{
S( ˜
Xij, ˜
Yij)}
.
Let M=
supα,v|
θ
0α
(v)| < ∞
based on the regularity condition (iii) ofAppendix B.1. Applying the Taylor expansion to the last term, it follows that n 2 −1|
U2(α) − ˜
U2(α)| ≤
M sup x,yˆ
S(
x,
y) −
S(
x,
y)
.
Based on the Glivenko–Cantelli theorem, the last term converges almost surely to zero uniformly in
α
. Next, we show that n 2 −1{
nU1(α) − ˜
U1(α)}
converges in probability to zero uniformly inα
. Similarly,|
nU1(α) − ˜
U1(α)| ≤
X
i<j˙
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[
1+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}]
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[ˆ
S( ˜
Xij, ˜
Yij) −
1/
n+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}/
n]
−
˙
θ
α{
S( ˜
Xij, ˜
Yij)}[
1+
θ
α{
S( ˜
Xij, ˜
Yij)}]
θ
α{
S( ˜
Xij, ˜
Yij)}
S( ˜
Xij, ˜
Yij)
+
X
i<j˙
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[
1+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}]
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}[ˆ
S( ˜
Xij, ˜
Yij) −
1/
n+
θ
α{ˆ
S( ˜
Xij, ˜
Yij)}/
n]
|
θ
α{ˆ
S( ˜
Xij, ˜
Yij)} − θ
α{
S( ˜
Xij, ˜
Yij)}|.
Now we show that the last two terms, multiplied by n 2 −1, converge in probability to zero uniformly in
α
. According to the Chebyshev inequality, it is sufficient to show that the quantitiesI1
=
E˙
θ
α{ˆ
S( ˜
X12, ˜
Y12)}[
1+
θ
α{ˆ
S( ˜
X12, ˜
Y12)}]
θ
α{ˆ
S( ˜
X12, ˜
Y12j)}[ˆ
S( ˜
X12, ˜
Y12) −
1/
n+
θ
α{ˆ
S( ˜
X12, ˜
Y12)}/
n]
−
˙
θ
α{
S( ˜
X12, ˜
Y12)}[
1+
θ
α{
S( ˜
X12, ˜
Y12)}]
θ
α{
S( ˜
X12, ˜
Y12)}
S( ˜
X12, ˜
Y12)
and I2
=
E˙
θ
α{ˆ
S( ˜
X12, ˜
Y12)}[
1+
θ
α{ˆ
S( ˜
X12, ˜
Y12)}]
θ
α{ˆ
S( ˜
X12, ˜
Y12j)}[ˆ
S( ˜
X12, ˜
Y12) −
1/
n+
θ
α{ˆ
S( ˜
X12, ˜
Y12)}/
n]
|
θ
α{ˆ
S( ˜
X12, ˜
Y12)} − θ
α{
S( ˜
X12, ˜
Y12)}|
converge to zero. Showing this requires substantial technical efforts due to the instability of the denominator terms. The key property in the proof is
E
1θ
α{
S( ˜
X12, ˜
Y12)}
S( ˜
X12, ˜
Y12)
< ∞,
which can be shown from Corollary 1 ofShih(1998), under the regularity conditions listed inAppendix B.1. After tedious evaluations similar to the previous arguments and applying the Glivenko–Cantelli theorem toS, we can show that I
ˆ
1and I2 converge to zero.B.3. Proof ofLemma 2
Lemma 1and the Taylor expansion for U2
(
eγ)
allow us to use the asymptotic expansionsn1/2
(γ
2−
γ ) =
(
−
n 2 −1α
U˙˜
2(α
∗)
)
−1(
n1/2 n 2 −1˜
U2(α) +
oP(
1)
)
,
whereU
˙˜
2(α) = ∂ ˜
U2(α)/∂α
andα
∗is betweenα
andα
ˆ
2. Based on the strong law of large number for U-statistics, the term−
n2
−1˙˜
U2
(α
∗)
converges almost surely to a constant A. Applying the same argument to U1(
eγ)
, we obtain the asymptotic expansion n1/2( ˆγ
1−
γ ) =
n1/2 n 2 −1˜
U1(α)
α
AL+
oP(
1),
n1/2( ˆγ
2−
γ ) =
n1/2 n 2 −1˜
U2(α)
α
A+
oP(
1).
Combining these formulas, we obtain the U-statistics expression n1/2
( ˆγ
1− ˆ
γ
2) =
n1/2 n 2 −1X
i<j 1α
˙θ
α{
S( ˜
Xij, ˜
Yij)}[θ
α{
S( ˜
Xij, ˜
Yij)} +
1]
ALθ
α{
S( ˜
Xij, ˜
Yij)}
S( ˜
Xij, ˜
Yij)
−
1 A! "
∆ij−
θ
α{
S( ˜
Xij, ˜
Yij)}
θ
α{
S( ˜
Xij, ˜
Yij)} +
1#
+
oP(
1).
Appendix C. Two examples of AC models
Example 1 (Clayton Model). The generating function can be written as
φ
α(v) = (v
−α−
1
)/α
forα ∈ (
0, ∞)
. The joint survival function can be written asPr
(
X>
x,
Y>
y) =
SX(
x)
−α+
SY(
y)
−α−
1 −1/α.
It follows that
τ = α/(α +
2)
andKα
(v) = v − φ
α(v)/φ
α0(v) = v + v(
1−
v
α)/α.
A special property of the Clayton model reflects in its local odds ratio, which can be expressed as
θ
∗(
x,
y) = α +
1. The odds ratio does not depend on(
x,
y)
.Example 2 (Gumbel Model). The generating function can be written as
φ
α(v) = {−
log(v)}
α+1, α ∈ [
0, ∞)
. The joint survival function can be written asPr
(
X>
x,
Y>
y) =
exp−
{−
log SX(
x)}
α+1+ {−
log SY(
y)}
α+1 α+11.
It follows that
τ = α/(α +
1)
andKα
(v) = v − φ
α(v)/φ
α0(v) = v − v
log(v)/(α +
1).
The local odds ratio can be expressed asθ
∗(
x,
y) =
1−
αlog S(x,y). Compared with the Clayton model, however,
θ
∗
(
x,
y)
depends on
(
x,
y)
. ReferencesClayton, D.G., 1978. A model for association in bivariate life tables and its application to epidemiological studies of familial tendency in chronic disease incidence. Biometrika 65, 141–151.
Dabrowska, D., 1988. Kaplan–Meier estimate on the plane. The Annals of Statistics 16, 1475–1489.
Dobric, J., Schmid, F., 2007. A goodness of fit test for copulas based on Rosenblatt’s transformation. Computational Statistics and Data Analysis 51, 4633–4642.
Duffy, D.L., Martin, N.G., Matthews, J.D., 1990. Appendectomy in Australian twins. American Journal of Human Genetics 47, 590–592. Frees, E.W., Valdez, E., 1998. Understanding the relationships using copulas. North American Actuarial Journal 2, 1–25.
Genest, C., Remillard, B., Beaudoin, D., 2009. Goodness-of-fit tests for copulas: a review and a power study. Insurance: Mathematics and Economics 44, 199–213.
Genest, C., Rivest, L.-P., 1993. Statistical inference procedures for bivariate Archimedean copulas. Journal of the American Statistical Association 88, 1034–1043.
Gill, R., Schumacher, M., 1987. A simple test of the proportional hazards assumption. Biometrika 74, 289–300.
Hsieh, , 2010. Estimation of Kendall’s tau from censored data. Computational Statistics and Data Analysis 54, 1613–1621. McGilchrist, C.A., Aisbett, C.W., 1991. Regression with frailty in survival analysis. Biometrics 47, 461–466.
Nelsen, R.B., 2006. An Introduction to Copulas, 2nd ed. In: Springer Series in Statistics, Springer-Verlag, New York.
Nikoloulopoulos, A.K., Karlis, D., 2008. Copula model evaluation based on parametric bootstrap. Computational Statistics and Data Analysis 52, 3342–3353. Oakes, D., 1982. A model for association in bivariate survival data. Journal of the Royal Statistical Society, Series B 44, 414–422.
Oakes, D., 1986. Semiparametric inference in a model for association in bivariate survival data. Biometrika 73, 353–361. Oakes, D., 1989. Bivariate survival models induced by frailties. Journal of the American Statistical Association 84, 487–493.
Prentice, R.L., Hsu, Li., 1997. Regression on hazard ratios and cross ratios in multivariate failure time analysis. Biometrika 84, 349–363. Shih, J.H., 1998. A goodness-of-fit test for association in a bivariate survival model. Biometrika 85, 189–200.
Van Der Vaart, A.W., 1998. Asymptotic Statistics. In: Cambridge Series in Statistics and Probabilistic Mathematics, Cambridge University Press, Cambridge. Wang, W., Wells, M., 2000. Model selection and semiparametric inference for bivariate failure-time data. Journal of the American Statistical Association
