Modelling the association in bivariate survival data by using a Bernstein copula

Abstract

Bivariate or multivariate survival data arise when a sample consists of clusters of two or more subjects which are correlated. This paper considers clustered bivariate survival data which is possibly censored. Two approaches are commonly used in modelling such type of correlated data: random effect models and marginal models. A random effect model includes a frailty model and assumes that subjects are independent within a cluster conditionally on a common non-negative random variable, the so-called frailty. In contrast, the marginal approach models the marginal distribution directly and then imposes a dependency structure through copula functions. In this manuscript, Bernstein copulas are used to account for the correlation in modelling bivariate survival data. A two-stage parametric estimation method is developed to estimate in the first stage the parameters in the marginal models and in the second stage the coefficients of the Bernstein polynomials in the association. Hereby we use a penalty parameter to make the fit desirably smooth. In this aspect linear constraints are introduced to ensure uniform univariate margins and we use quadratic programming to fit the model. We perform a Simulation study and illustrate the method on a real data set.

Appendices

Appendix A: Penalty matrix for penalizing second order derivatives

In this appendix we show the penalty matrix based on the second order derivatives.

$$\begin{aligned}&\int \Big (\frac{\partial ^2c_B(u_1,u_2)}{(\partial u_1)^2}\Big )^2 du_1du_2\\&\quad ={\varvec{\beta }^T}\int \Big [\frac{\partial ^2}{(\partial u_1)^2}\Big \{{\varvec{\phi }}_{\mathbf{m}}(u_{1})\otimes {\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big \}\Big ]^T \Big [\frac{\partial ^2}{(\partial u_1)^2}\Big \{{\varvec{\phi }}_{\mathbf{m}}(u_{1})\otimes {\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big \}\Big ]du_1du_2{\varvec{\beta }}\\&\quad = {\varvec{\beta }^T} \int \Big [\Big (\frac{\partial ^2}{(\partial u_1)^2}{\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big )^T\Big (\frac{\partial ^2}{(\partial u_1)^2}{\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big )\Big ]du_1 \otimes \Big [\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big )^T\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big )\Big ]{\varvec{\beta }}\\&\quad ={\varvec{\beta }^T}P_{u_1}{\varvec{\beta }} \end{aligned}$$

where,

$$\begin{aligned} P_{u_1} =\int \Big [\Big (\frac{\partial ^2}{(\partial u_1)^2}{\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big )^T\Big (\frac{\partial ^2}{(\partial u_1)^2}{\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big )\Big ]du_1 \otimes \Big [\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big )^T\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big )\Big ] \end{aligned}$$

The integral of the second order derivatives of Bernstein polynomials are calculated easily. The second order derivative of (3) equals (Doha et al. 2011)

$$\begin{aligned} \frac{\partial ^2}{(\partial u)^2}\phi _{mv}(u)=\frac{(m+1)!}{(m-2)!}\sum _{k = \max (0, v+2-m)}^{min(v,2)}(-1)^{k+2}\left( {\begin{array}{c}2\\ k\end{array}}\right) \phi _{m-2,v-k}(u) \end{aligned}$$

This is rewritten as,

$$\begin{aligned} \frac{\partial ^2}{(\partial u)^2}\phi _{mv}(u)=(\phi _{m-2,v}(u)B)w \end{aligned}$$

with

$$\begin{aligned} B= \begin{bmatrix} 1 &{} -2 &{} 1 &{} 0 &{}\ldots &{} 0 \\ 0 &{} 1 &{} -2 &{} 1 &{}\ddots &{} \vdots \\ \vdots &{} \ddots &{}\ddots &{}\ddots &{} \ddots &{} 0\\ 0 &{} \ldots &{} 0 &{} 1 &{} -2 &{} 1\\ \end{bmatrix} ,\quad B\in {\mathbf {R}}^{(m-1)\times (m+1)} \end{aligned}$$

and $w=\frac{(m+1)!}{(m-2)!}$. Therefore, the matrix $P_{u_1}$ and $P_{u_2}$ are equivalent to

$$\begin{aligned} P_{u_1}&= \Big [wB^T\int \phi _{m-2,v_1}(u_1)\phi _{m-2,v_1}(u_1)du_1Bw\Big ] \otimes \Big [\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big )^T{\varvec{\phi }}_{\mathbf{m}}(u_{2})\Big ]\\ P_{u_2}&= \Big [\Big ({\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big )^T{\varvec{\phi }}_{\mathbf{m}}(u_{1})\Big ]\otimes \Big [wB^T\int \phi _{m-2,v_2}(u_2)\phi _{m-2,v_2}(u_2)du_2Bw\Big ] \end{aligned}$$

So, the penalty can be written as a quadratic form $\lambda \varvec{\beta }^T P_{int} \varvec{\beta }$, where $\lambda$ is the penalty parameter steering the amount of smoothness and $P_{int}=P_{u_1}+P_{u_2}$.

Appendix B: tables and Graphs

See Tables 6 and 7 and Fig. 5

Table 6 Corrected AIC of the estimated model based on data-driven $\lambda$ from the simulated data with 25% censoring

Table 7 KLD and IAE of the estimated model based on data-driven $\lambda$ and penalized integrated second order derivatives

Fig. 5

Density difference plots for the Clayton (first row), Gumbel (second row), Frank (third row) and Gaussian copula (fourth row) with 25% censoring