WEIBULL EXTENSIONS OF BIVARIATE EXPONENTIAL REGRESSION MODEL FOR A ...
WEIBULL EXTENSIONS OF BIVARIATE EXPONENTIAL
REGRESSION MODEL FOR A SURVIVAL DATA
DAVID D. HANAGAL
Department of Statistics, University of Pune, Pune-411007, India.
Email: ddh@stats.unipune.ernet.in, david−hanagal@yahoo.co.in
Abstract
In this paper, we propose two bivariate Weibull regression models for the survival data
based on bivariate exponential distribution of Gumbel (1960). There are some biometrical
applications which motivated to study these particular models. We obtain the estimation of
regression parameters and derive the test procedure for testing the significance of regression
parameters based on censored samples.
Key words: Bivariate Weibull model, Parametric regression, Survival times.
1
Introduction
Gumbel (1960) proposed bivariate exponential (BVE) distribution which has interesting prop-
erties like marginal exponentials and absolute continuity. Lee (1979) extended this BVE model
to bivariate Weibull (BVW) distribution by introducing shape parameters as done by Lu (1989)
and Hanagal (1996) for BVE of Marshall-Olkin (1967) model. Hanagal (2004, 2005) proposed
bivariate regression analysis based on Weibull model with identical covariates for both com-
ponents. There are some situations where we find some non-identical covariates in addition
to identical covariates for a paired components in a system. For example, failure times of a
pair of dental implants in a jaw of a patient. Here the identical covariates are age and sex of
a patient and non-identical covariates may be 1) different materials (ceramic, metal) of dental
implant, 2) different shapes (screw, anchor, pillar, hollow cylinder) of dental implant, 3) dental
implants in different locations (front, premolar, molar) of a jaw, 4) dental implants in different
jaws (lower, upper).
Dental implant study in medical science is one branch called Implantology. Implantology
studies the scientific technique of installation of dental implant through surgical procedure after
loss of natural tooth. After installation of dental implant, one usually waits three to six months
1
for the placement of an artificial tooth on this implant. With some procedures, placement of
an artificial tooth after the installation of the implant is immediate. Implants within a patient
are dependent and their life times are correlated. See Haas et al (1996) and Ivanoff et al (1999).
We propose BVW regression model based on identical and non-identical covariates on both
components. Here the censoring of life times of two implants is due to withdrawals or death
of a patient or termination of the study. The censoring time is independent of the life times
of the paired implants. Here we can use homogeneous (univariate) censoring scheme given by
Hanagal (1992a, 1992b) because withdrawals or death of a patient or termination of the study
will censor both life times of dental implants. We try to collect data on this but it is incomplete
due to very few failures of dental implants, very high censoring and not an adequate sample
size. Unfortunately we are not in a position to give real data for our proposed model.
The objective of this paper is to focus on the BVW regression model for a paired components
which are randomly censored, estimate the regression parameters and test the significance of
regression parameters. In Section 2, we introduce BVW regression model and in Section 3, we
obtain estimation of regression parameters in the proposed model. We present test procedures
for testing the significance of regression parameters. In Section 4, we present another Weibull
extension of Gumbel (1990) with regression covariates.
2
Bivariate Weibull Regression Model
The Weibull extension of BVE-type I of Gumbel (1960) with survival function is given by
F (t1, t2) = P [T1 > t1, T2 > t2]
= exp{−λ1tc1 − λ
− δλ
tc2},
t
1
2tc2
2
1λ2tc1
1
2
1, t2 > 0
(1)
where λ1, λ2 > 0, 0 ≤ δ ≤ 1
Here the marginal distribution of T1 and T2 are distributed as Weibull with scale parameters
λ1 and λ2 and shape parameters c1 and c2 respectively. The parameter δ corresponds to the
dependence parameter in BVW model.
As we have seen in the univariate Weibull regression, the scale parameter of the univariate
Weibull distribution can be expressed in terms of regression coefficients.
If λ is the scale
parameter of the exponential distribution, then λ = e−β z or λ = eβ z where β is the vector of
2
regression parameters and z is the vector of regressors or covariates. In the similar manner, the
scale parameters λ1 and λ2 can be expressed in terms of regression parameters in the following
way.
λ1 = exp{−(β z
z
0 0 + β1 1)}
λ2 = exp{−(β z
z
0 0 + β2 2)}
(2)
where
β
= (β
0
01, ...., β0p),
−∞ < β0 < ∞
β
= (β
1
11, ...., β1q ),
−∞ < β1 < ∞
β
= (β
2
21, ...., β2q ),
−∞ < β2 < ∞
z
= (z
0
01, ...., z0p)
z
= (z
1
11, ...., z1q )
z
= (z
2
21, ...., z2q )
The exponent terms in the above expressions, we can take either positive or negative but
in either case λ1, λ2 > 0. β z
0 0 corresponds to the term containing identical covariates for both
components. β z
z
1 1 corresponds to the term containing covariates for first component and β2 2
corresponds to the term containing covariates for second component.
Now the survival function of BVW in terms of covariates is given by
F (t
z0+β z
z
z
z
z
z
1 1)
0+β2 2)
0+β1 1+β2 2)
1, t2)
= exp{−tc1e−(β0
− tc2e−(β0
− δtc1tc2e−(2β0
}
(3)
1
2
1
2
3
Estimation of the Parameters
For the bivariate life time distribution, we use univariate censoring scheme given by Hanagal
(1992a, 1992b) because the individuals do not enter at the same time and withdrawal or death
of an individual will censor both life times of the components. Here the censoring time is
independent of the life times of both components. This is the standard univariate right censoring
for both failure times T1 and T2.
3
Suppose that there are n independent pairs of components under study and i-th pair of the
components have life times (t1i, t2i) and a censoring time (yi). The life times associated with
i-th pair of the components are given by
(T1i, T2i) = (t1i, t2i),
max(t1i, t2i) < yi
= (t1i, yi),
t1i < yi < t2i
= (yi, t2i),
t2i < yi < t1i
= (yi, yi),
yi < min(t1i, t2i).
(4)
Discarding factors which do not contain any of the parameters, we want to estimate the pa-
rameters in the proposed model. Now the likelihood of the sample of size n is given by
n1
n2
n3
n4
L = (
f1)(
f2)(
f3)(
F )
(5)
i=1
i=1
i=1
i=1
where
f
z0i+β z
z
z
1 1i)
0i+β2 2i)
1
= c1c2tc1−1tc2−1e−(β0
e−(β0
1i
2i
[(1 + δtc2e−(β z
z
z
z
0 0i+β2 2i))(1 + δtc1 e−(β0 0i+β1 1i)) − δ]F (t
2i
1i
1i, t2i),
0 < t1i, t2i < yi
f
z0i+β z
z
z
1 1i)
0i+β2 2i)
2
= c1tc1−1e−(β0
(1 + δyc2e−(β0
)F (t
1i
i
1i, yi),
0 < t1i < yi < t2i
f
z0i+β z
z
z
2 2i)
0i+β1 1i)
3
= c2tc2−1e−(β0
(1 + δyc1e−(β0
)F (y
2i
i
i, t2i),
0 < t2i < yi < t1i
F
= P [T1i > yi, T2i > yi],
0 < yi < min(t1i, t2i)
n1, n2, n3 and n4 are the random number of observations observed to fail in the range space
corresponding to f1, f2, f3 and F respectively. f1 is the pdf with respect to Lebesque measure
in R2 and f2 and f3 are the pdf with respect to Lebesque measure in R1.
The likelihood equations can be obtained by taking first order partial derivatives of the
loglikelihood and equating to zero. The likelihood equations are not easy to solve. One may
obtain maximum likelihood estimators (MLEs) by Newton-Raphson procedure. But we came
to know from the simulation study in Section 3.2, the likelihood equations do not converge at
all in the Newton-Raphson procedure and the method of maximum likelihood (ML) fails to
estimate all the parameters simultaneously. One can obtain estimates of the parameters by two
stage MLE method or conditional MLE method. In the first stage, estimate the parameters
c1, c2, δ by ML method under the base line model by conditioning β0 = β1 = β2 = 0 and then
in the second stage estimate the parameters β0, β1, β2 by ML method after substituting MLEs
4
of c1, c2, δ obtained from the first stage. These are conditional MLEs. When (i) there are no
consistent estimators in closed form (ii) iterative procedures fail to converge, one can adopt two
stage MLE procedure for estimating the parameters. We are mainly interested in estimating
and testing the regression parameters and the other parameters (c1, c2, δ) are involved in the
base line model which are here nuisance parameters. The observed fisher information matrix
I1 with appropriate second order partial derivatives based on first stage is
∂2logL
∂2logL
∂2logL
∂c2
∂c
∂c
1
1∂c2
1∂δ
I
∂2logL
∂2logL
∂2logL
1
= −
.
(6)
∂c2∂c1
∂c2
∂c2∂δ
2
∂2logL
∂2logL
∂2logL
∂δ∂c1
∂δ∂c2
∂δ2
The observed fisher information matrix I2 which is of the order (p+2q)×(p+2q) with appropriate
second order partial derivatives based on the second stage is
∂2logL
∂2logL
∂2logL
∂β0i∂β0j
∂β0i∂β1j
∂β0i∂β2j
I
∂2logL
∂2logL
∂2logL
2
= −
.
(7)
∂β1i∂β0j
∂β1i∂β1j
∂β1i∂β2j
∂2logL
∂2logL
∂2logL
∂β2i∂β0j
∂β2i∂β1j
∂β2i∂β2j
The inverse of the above observed information matrix (I2) is the observed variance-covariance
matrix ( ˆ
Σ11 = I−1) of the MLEs ˆ
β = ( ˆ
β01, ..., ˆ
β0p, ˆ
β11, ..., ˆ
β1q, ˆ
β21, ..., ˆ
β2q) of the parameters
β = (β01, ..., β0p, β11, ..., β1q, β21, ..., β2q) .
√
Thus
n( ˆ
β − β) has asymptotic multivariate normal distribution with mean vector zero and
variance-covariance matrix Σ11, where Σ11 is (p + 2q) × (p + 2q) variance covariance matrix of
ˆ
β = ( ˆ
β01, ..., ˆ
β0p, ˆ
β11, ..., ˆ
β1q, ˆ
β21, ..., ˆ
β2q)
3.1
Test for Regression Coefficients
The hypotheses about β can be frequently put in the form Ho : β11 = 0, with β partitioned as
β = (β11, β22) where θ11 is k × 1, (k < p + 2q). To test H0 against the alternative that β11 = 0
one can use
Λ
ˆ
ˆ
1
=
ˆ
β Σ−1β
11
22
11
(8)
where ˆ
Σ22 is k × k asymptotic observed variance-covariance matrix of ˆ
β11. Under H0, Λ1 is
asymptotically chi-square with k d.f..
5
3.2
Simulation study
We generate 1000 samples of sizes n=40, 60 and 80 from BVW model and obtain conditional
MLEs of the parameters based on first stage. We observed from the simulation study as in
Table 1 that MLEs are very close to the known values of the parameters in BVW model. In
Table 2, we obtain regression parameters and also obtain the power of test statistics based on
chisquare test at the levels of significance (α) = .01 and .05. The following are three hypotheses
of the tests discussed in Section 3.1.
1) H0: β01 = 0 Vs H1: β01 = .5
2) H0: β11 = 0 Vs H1: β11 = .5
3) H0: β21 = 0 Vs H1: β21 = .5
The first test is for testing common regression parameter corresponding to both components
equal to zero, the second test is for testing the regression parameter corresponding to first
component equal to zero and the third test is for testing the regression parameter corresponding
to second component equal to zero. It is observed from Table 2 that the tests are very powerful
as sample size approaches to 80. The distribution of the censoring time is taken as exponential
with failure rate .03.
6
Table 1: MLEs of the Parameters in First BVW Model.
Parameters
c1
c2
δ
values
1.5
2.0
0.5
n = 40
Est.
1.5018
2.0129
0.5016
s.e.
0.0195
0.0210
0.0315
n = 60
Est.
1.5010
2.0125
0.5012
s.e.
0.0185
0.0198
0.0257
n = 80
Est.
1.5009
2.0116
0.5008
s.e.
0.0174
0.0181
0.0216
7
Table 2: Conditional MLEs of the Regression Parameters and
Power of the Tests in First BVW Regression Model.
Parameters
β01
β11
β21
values
0.5
0.5
0.5
n = 40
Est.
0.5018
0.5013
0.5017
s.e.
0.0197
0.0198
0.0342
Power, α = .01
0.697
0.712
0.654
Power, α = .05
0.771
0.797
0.753
n = 60
Est.
0.5016
0.5012
0.5020
s.e.
0.0172
0.0171
0.0288
Power, α = .01
0.736
0.756
0.687
Power, α = .05
0.854
0.868
0.814
n = 80
Est.
0.5010
0.5005
0.5009
s.e.
0.0153
0.0152
0.0257
Power, α = .01
0.826
0.848
0.789
Power, α = .05
0.914
0.947
0.898
4
Another Weibull Extension of Gumbel
The another Weibull extension of BVE-type III of Gumbel (1960) is given by
F (t1, t2) = P [T1 > t1, T2 > t2]
= exp{−(λ1tc1/δ + λ
)δ},
t
1
2tc2/δ
2
1, t2 > 0
(9)
where λ1, λ2 > 0, 0 ≤ δ ≤ 1
Here the marginal distribution of T1 and T2 are distributed as Weibull with scale parameters
λ1 and λ2 and shape parameters c1 and c2 respectively. The parameter δ corresponds to the
dependence parameter in BVW model.
The survival function in terms of covariates is given by
8
F (t
z0+β z
z
z
1 1)
0+β2 2)
1, t2)
= exp{−(tc1/δe−(β0
+ tc2/δe−(β0
)δ}
(10)
1
2
Now we want to estimate the parameters in the proposed model. The likelihood of the
sample of size n is written in the same way as in Eqn (5) with
f
z0i+β z
z
z
1 1i)
0i+β2 2i)
1
= c1c2tc1/δ−1tc2/δ−1e−(β0
e−(β0
1i
2i
[(tc1/δ−1e−(β z
z
z
z
0 0i+β1 1i) + tc2/δ−1e−(β0 0i+β2 2i))δ + 1 − 1/δ]
1i
2i
(tc1/δ−1e−(β z
z
z
z
0 0i+β1 1i) + tc2/δ−1e−(β0 0i+β2 2i))δ−2F (t
1i
2i
1i, t2i),
0 < t1i, t2i < yi
f
z0i+β z
z
z
1 1i)
0i+β2 2i)
2
= (tc1/δ−1e−(β0
+ yc2/δ−1e−(β0
)δ−1
1i
i
c
z0i+β z
1 1i)
1tc1/δ−1e−(β0
F (t
1i
1i, yi),
0 < t1i < yi < t2i
f
z0i+β z
z
z
1 1i)
0i+β2 2i)
3
= (tc1/δ−1e−(β0
+ yc1/δ−1e−(β0
)δ−1
1i
i
c
z0i+β z
2 2i)
2tc2/δ−1e−(β0
F (y
2i
i, t2i),
0 < t2i < yi < t1i
F
= P [T1i > yi, T2i > yi],
0 < yi < min(t1i, t2i)
Here also the likelihood equations are not easy to solve. We observed from the simulation
study that Newton-Raphson procedure for the likelihood equations do not converge and hence
we proceed to two stage MLE method as discussed in Section 3. We obtain test for regression
parameters in the similar way as in Section 3.1 and do the simulation study based on the sample
sizes n = 40, 60 and 80. Table 3 and 4 give the parameter estimates and power of the tests of
hypothesis discussed in Section 3.2.
9
Table 3: MLEs of the Parameters in Second BVW Model.
Parameters
c1
c2
δ
values
1.5
2.0
0.5
n = 40
Est.
1.5027
2.0029
0.5032
s.e.
0.0243
0.0256
0.0287
n = 60
Est.
1.5019
2.0022
0.5026
s.e.
0.0172
0.0171
0.0288
n = 80
Est.
1.5011
2.0015
0.5017
s.e.
0.0153
0.0162
0.0257
10
Table 4: Conditional MLEs of the Regression Parameters and
Power of the Tests in Second BVW Regression Model.
Parameters
β01
β11
β21
values
0.5
0.5
0.5
n = 40
Est.
0.5031
0.5033
0.5027
s.e.
0.0157
0.0148
0.0182
Power, α = .01
0.7432
0.7654
0.7832
Power, α = .05
0.8321
0.8452
0.8577
n = 60
Est.
0.5024
0.5027
0.5018
s.e.
0.0142
0.0131
0.0125
Power, α = .01
0.8342
0.8424
0.8528
Power, α = .05
0.8924
0.9107
0.9216
n = 80
Est.
0.5010
0.5014
0.5006
s.e.
0.0133
0.0122
0.0117
Power, α = .01
0.9122
0.9218
0.9287
Power, α = .05
0.9532
0.9611
0.9637
5
DISCUSSIONS:
I have simulated 1000 samples each of size n = 40, 60 and 80. If I take smaller sample sizes
for the simulation, there is a problem of convergence of estimates of the parameters in N-
R procedure. In the survival data, one should remember that the number of failures should
be less than the sample size. In the simulation process, the percentage of censoring changes
from sample to sample for fixed sample size. So the effective sample size for the parametric
model is the number of failures. In our case, we have a bivariate Weibull model with three
parameters under the base line model. The sample sizes 20 and 30 are very small for bivariate
Weibull model with three parameters and censoring scheme. The efficiency and convergence of
estimators depend on three things as follows:
1) sample size,
11
2) percentage of censoring,
3) number of parameters in the model.
When the sample size is very small and it is highly censored and there are more number
of parameters in the model, the probability of convergence of the estimates of the parameters
is very less. If we take into account the above things, the power of the tests based on these
estimates will perform better.
ACKNOWLEDGEMENTS
I thank the referee for the useful suggestions and comments which resulted in an improved
version of an earlier paper.
References
References
[1] Gumbel, E.J. (1960). Bivariate exponential distribution. Journal of the Amer. Statist.
Assoc., 55, 698-707. Haas, R.,
[2] Mensdorff-Pouilly, N., Mailath, G. and Watzek, G.(1996). Survival of 1,920 IMZ implants
followed for up to 100 months. Int. J. Oral Maxillofac Implants 11, 581-88.
[3] Hanagal, D.D. (1992a). Some inference results in bivariate exponential distributions based
on censored samples. Comm.Statist., Theory and Methods, 21, 1273-95.
[4] Hanagal, D.D. (1992b). Some inference results in modified Freund’s bivariate exponential
distribution. Biometrical Journal, 34(6), 745-56.
[5] Hanagal, D.D. (1996). A multivariate Weibull distribution. Economic Quality Control, 11,
193-200.
[6] Hanagal, D.D. (2004). Parametric bivariate regression analysis based on censored samples:
A Weibull model. Economic Quality Control, 19, 83-90.
[7] Hanagal, D.D. (2005). Bivariate Weibull regression model based on censored samples.
Statistical Papers, (to appear).
12
[8] Ivanoff, C.T., Grondahl, K., Sennerby, L., Bergstorm, C. and Lekholm, W.(1999). Influence
of variations in implant diameters: a 3 to 5 year retrospective clinical report. Int. J. Oral
Maxillofac Implants, 14, 173-80.
[9] Lee, L. (1979). Multivariate distributions having Weibull properties. Journal of Multivari-
ate Analysis, 9, 267-77.
[10] Lu, J.C. (1989). Weibull extensions of the Freund and Marshall-Olkin bivariate exponential
models. IEEE Transactions on Reliability, 38, 615-19.
[11] Marshall, A.W. and Olkin, I.(1967). A multivariate exponential distribution. Journal of
Amer. Statist. Assoc., 62, 30-44.
13
