Bivariate Chen distributions with BvCh

Introduction

BvCh provides the probability density, survival and hazard functions, random generation, graphical displays and maximum likelihood estimation for three bivariate Chen type lifetime models:

Model	Construction	Singular part?
BvCh	Marshall-Olkin type, `min` of three Chen variables	yes, on the diagonal
BEC	Marshall-Olkin type with extended Chen variables, adding `lambda`	yes, on the diagonal
IFGM-Chen	Inverse FGM copula with Chen marginals	no

All univariate Chen building blocks are taken from the reliaR package.

library(BvCh)

The BvCh model

With \(U_i \sim \mathrm{Chen}(\alpha_i,\beta)\) independent and \(Z_1=\min(U_1,U_3)\), \(Z_2=\min(U_2,U_3)\), the joint survival function is

\[S(z_1,z_2)=e^{\alpha_1(1-e^{z_1^\beta})}e^{\alpha_2(1-e^{z_2^\beta})}e^{\alpha_3(1-e^{w^\beta})}, \qquad w=\max(z_1,z_2),\]

which is a product of three univariate Chen survival functions. The density is piecewise, with a branch on each side of the diagonal and a singular branch on the diagonal itself.

a1 <- 1.5; a2 <- 2.5; a3 <- 0.5; b <- 1.8
z1 <- c(0.4, 0.9, 0.7, 1.3)
z2 <- c(0.6, 0.4, 0.7, 1.1)
cbind(z1, z2,
      pdf = dbvch(z1, z2, a1, a2, a3, b),
      surv = sbvch(z1, z2, a1, a2, a3, b),
      haz = hbvch(z1, z2, a1, a2, a3, b))
#>       z1  z2          pdf         surv        haz
#> [1,] 0.4 0.6 1.4067219426 1.673615e-01   8.405291
#> [2,] 0.9 0.4 0.8899046491 4.488011e-02  19.828487
#> [3,] 0.7 0.7 0.0507469024 4.431486e-02   1.145144
#> [4,] 1.3 1.1 0.0004204033 1.196383e-06 351.395199

The survival function can be checked against the product form, and the hazard is exactly \(f/S\).

S <- sbvch(z1, z2, a1, a2, a3, b)
Sref <- schen(z1, b, a1) * schen(z2, b, a2) * schen(pmax(z1, z2), b, a3)
c(max_abs_diff_surv = max(abs(S - Sref)),
  max_abs_diff_haz  = max(abs(hbvch(z1, z2, a1, a2, a3, b) -
                                dbvch(z1, z2, a1, a2, a3, b) / S)))
#> max_abs_diff_surv  max_abs_diff_haz 
#>      2.775558e-17      1.776357e-15

The mixture representation

Proposition 1 writes the density as a mixture of an absolutely continuous part and a singular part on the diagonal. Both components are available directly and the mixture identity holds to machine precision.

K <- a1 + a2 + a3
ac <- dbvch_ac(z1, z2, a1, a2, a3, b)
sg <- dbvch_sing(z1, z2, a1, a2, a3, b)
mix <- (a1 + a2) / K * ac + a3 / K * sg
data.frame(ac = ac, sing = sg, mixture = mix, dbvch = dbvch(z1, z2, a1, a2, a3, b))
#>             ac      sing      mixture        dbvch
#> 1 1.5825621855 0.0000000 1.4067219426 1.4067219426
#> 2 1.0011427303 0.0000000 0.8899046491 0.8899046491
#> 3 0.0000000000 0.4567221 0.0507469024 0.0507469024
#> 4 0.0004729538 0.0000000 0.0004204033 0.0004204033

Mass carried by each region

Numerically integrating the density over the three regions recovers \(\alpha_i/K\) exactly, confirming that the density is properly normalised.

m1 <- integrate(function(z) vapply(z, function(v)
  integrate(function(v2) dbvch(v, v2, a1, a2, a3, b), v, Inf,
            rel.tol = 1e-10, subdivisions = 2000L)$value, numeric(1)),
  0, Inf, rel.tol = 1e-8, subdivisions = 2000L)$value
m2 <- integrate(function(z) vapply(z, function(v) if (v > 0)
  integrate(function(v2) dbvch(v, v2, a1, a2, a3, b), 0, v,
            rel.tol = 1e-10, subdivisions = 2000L)$value else 0, numeric(1)),
  0, Inf, rel.tol = 1e-8, subdivisions = 2000L)$value
u <- seq(1e-8, 6, length.out = 200001L)
m0 <- a3 / K * sum(reliaR::dchen(u, b, K)) * (u[2] - u[1])
rbind(f1 = m1, f2 = m2, diagonal = m0, total = m1 + m2 + m0,
      expected = c(a1/K, a2/K, a3/K, 1))
#>               [,1]      [,2]      [,3]      [,4]
#> f1       0.3333333 0.3333333 0.3333333 0.3333333
#> f2       0.5555556 0.5555556 0.5555556 0.5555556
#> diagonal 0.1111111 0.1111111 0.1111111 0.1111111
#> total    1.0000000 1.0000000 1.0000000 1.0000000
#> expected 0.3333333 0.5555556 0.1111111 1.0000000

Simulation

Random generation draws three independent Chen variables and applies the Marshall-Olkin construction.

sim <- rbvch(100000, a1, a2, a3, b)
q <- seq(0.4, 1.6, length.out = 6)
data.frame(q = q,
           theo = 1 - schen(q, b, a1) * schen(q, b, a3),
           emp  = sapply(q, function(v) mean(sim[, 1] <= v)))
#>      q      theo     emp
#> 1 0.40 0.3454302 0.34552
#> 2 0.64 0.6769215 0.67878
#> 3 0.88 0.9116514 0.91249
#> 4 1.12 0.9919098 0.99189
#> 5 1.36 0.9999162 0.99996
#> 6 1.60 1.0000000 1.00000

The frequency of ties also matches \(\alpha_3/K\): ## The BEC model

The Bivariate Extended Chen model replaces the Chen building block with the extended Chen distribution, whose survival function is \([1+\lambda(e^{x^\beta}-1)]^{-\alpha}\).

x1 <- c(0.4, 0.9, 0.7, 1.3); x2 <- c(0.6, 0.4, 0.7, 1.1)
cbind(x1, x2,
      pdf  = dbec(x1, x2, a1, a2, a3, b, lambda = 2),
      surv = sbec(x1, x2, a1, a2, a3, b, lambda = 2),
      haz  = hbec(x1, x2, a1, a2, a3, b, lambda = 2))
#>       x1  x2         pdf         surv        haz
#> [1,] 0.4 0.6 0.904656287 0.0758479551 11.9272337
#> [2,] 0.9 0.4 0.504457568 0.0323650572 15.5864878
#> [3,] 0.7 0.7 0.019214169 0.0200094337  0.9602555
#> [4,] 1.3 1.1 0.004864604 0.0001719257 28.2947963

Off the diagonal the joint hazard factorises into the componentwise hazard gradients:

off <- x1 != x2
h12 <- hbec1(x1, x2, a1, a2, a3, b, 2) * hbec2(x1, x2, a1, a2, a3, b, 2)
data.frame(hbec = hbec(x1, x2, a1, a2, a3, b, 2),
           hbec1_times_hbec2 = h12,
           off_diagonal = off)
#>         hbec hbec1_times_hbec2 off_diagonal
#> 1 11.9272337          11.92723         TRUE
#> 2 15.5864878          15.58649         TRUE
#> 3  0.9602555          74.68934        FALSE
#> 4 28.2947963          28.29480         TRUE

The IFGM-Chen model

This model pairs the Inverse FGM copula with Chen marginals and is absolutely continuous.

t1 <- 1.2; b1 <- 1.3; t2 <- 1.5; b2 <- 1.1; g <- 0.6; om <- 0.4
cbind(x = 0.5, y = 0.9,
      reliability = relifgmchen(0.5, 0.9, t1, b1, t2, b2, g, om),
      pdf  = difgmchen(0.5, 0.9, t1, b1, t2, b2, g, om),
      haz  = hifgmchen(0.5, 0.9, t1, b1, t2, b2, g, om))
#>        x   y reliability       pdf      haz
#> [1,] 0.5 0.9  0.05424035 0.5088101 9.380656

data.frame(C = c(cifgm(0.3, 0.7, g, om), cifgm(0.5, 0.5, g, om)),
           tau = tau_ifgm(g, om), rho = rho_ifgm(g, om))
#>           C         tau        rho
#> 1 0.1963698 -0.06989602 -0.1048831
#> 2 0.2306328 -0.06989602 -0.1048831

data.frame(gamma = c(0, 0.3, 0.6, 0.9, 1),
           tau = sapply(c(0, 0.3, 0.6, 0.9, 1), tau_ifgm, omega = 0.4),
           rho = sapply(c(0, 0.3, 0.6, 0.9, 1), rho_ifgm, omega = 0.4))
#>   gamma         tau         rho
#> 1   0.0 -0.01097992 -0.01646968
#> 2   0.3 -0.04183077 -0.06280926
#> 3   0.6 -0.06989602 -0.10488312
#> 4   0.9 -0.09572557 -0.14350074
#> 5   1.0 -0.10390898 -0.15571040

With the copula formula exactly as given in the source document, \(\tau\) and \(\rho\) are non-positive for \(\gamma,\omega \geq 0\). The source quotes \(\tau \approx 0.156>0\) for the same parameters, which is inconsistent with the copula formula it also quotes. See ?cifgm for details.

Estimation

fitbvch, fitbec and fitifgmchen estimate the models by maximum likelihood. For the Marshall-Olkin type models the observed-data log-likelihood is exact, because the contribution of an observation is fully determined by which side of the diagonal it falls on, so direct numerical maximisation reaches the same point as the EM algorithm.

d <- rbvch(1000, a1, a2, a3, b)
fit <- fitbvch(d)
fit
#> 
#> BvCh model fit by maximum likelihood
#> 
#> Call: fitbvch(x1 = d)
#> 
#> Coefficients:
#>        Estimate Std. Error
#> alpha1   1.4773     0.0671
#> alpha2   2.5024     0.1069
#> alpha3   0.5197     0.0443
#> beta     1.8739     0.0388
#> 
#> n = 1000   (z1 < z2: 315, z1 > z2: 572, z1 = z2: 113)
#> -log L = 193.2   AIC = 394.5   BIC = 414.1

Standard errors come from the numerical Hessian and are reported on the natural scale.

coef(fit)
#>   alpha1   alpha2   alpha3     beta 
#> 1.477284 2.502355 0.519653 1.873876
fit$se
#> [1] 0.06711179 0.10690942 0.04432055 0.03882814

Goodness of fit

Because \(Z_1=\min(U_1,U_3)\), the marginal of \(Z_1\) is \(\mathrm{Chen}(\alpha_1+\alpha_3,\beta)\); likewise for \(Z_2\), and \(\min(Z_1,Z_2)=\min(U_1,U_2,U_3)\) is \(\mathrm{Chen}(\alpha_1+\alpha_2+\alpha_3,\beta)\).

ksfit_bvch(fit, d)
#>       Series    alpha     beta          D   p.value      mean        sd
#> 1         Z1 1.996937 1.873876 0.01767393 0.9136015 0.5331568 0.2452298
#> 2         Z2 3.022008 1.873876 0.02063283 0.7881814 0.4439171 0.2192070
#> 3 min(Z1,Z2) 4.499292 1.873876 0.01843417 0.8860342 0.3708482 0.1871929

Model comparison

f1 <- fitbvch(rbvch(500, a1, a2, a3, b))
f2 <- fitbec(rbec(500, a1, a2, a3, b, 2))
f3 <- fitifgmchen(rifgmchen(500, t1, b1, t2, b2, g, om))
compare_models(f1, f2, f3)
#>       Model nobs  neglogL      AIC      BIC
#> 1      BvCh  500 105.6994 219.3988 236.2573
#> 2       BEC  500 110.8234 231.6467 252.7198
#> 3 IFGM-Chen  500 190.0897 392.1793 417.4670

The IFGM-Chen copula parameters gamma and omega are only weakly identified: the marginal parameters are estimated precisely but the likelihood is nearly flat in the copula direction. See ?fitifgmchen.

Graphics

The absolutely continuous part of the BvCh density is plotted for a range of \(\beta\) values; as \(\beta\) increases the surface becomes more peaked near the origin.

plot_surface_bvch(alpha1 = a1, alpha2 = a2, alpha3 = a3,
                  beta.vec = c(0.8, 1.8, 5, 10))

Contour plot of the absolutely continuous part, with the singular support marked.

plot_contour_bvch(alpha1 = a1, alpha2 = a2, alpha3 = a3, beta = b)

Reliability surface of the IFGM-Chen model.

plot_reliability_ifgmchen(theta1 = t1, beta1 = b1, theta2 = t2, beta2 = b2,
                          gamma = g, omega = om)

References

Chen, Z. (2000). A new two-parameter lifetime distribution with bathtub shape or increasing failure rate function. Statistics & Probability Letters, 49, 155-161.
Nelsen, R. B. An Introduction to Copulas, 2nd edition, Springer, 2006.
Klein, M. and Schumacher, C. (2002). Bivariate survival and three dimensional Marshall-Olkin type copulas. Journal of the Royal Statistical Society B, 64, 721-734.

c(empirical = mean(sim[, 1] == sim[, 2]), expected = a3 / K)
#> empirical  expected 
#> 0.1120200 0.1111111