Generalized inverse-Gaussian frailty models with application to TARGET neuroblastoma data

Abstract

A new class of survival frailty models based on the generalized inverse-Gaussian (GIG) distributions is proposed. We show that the GIG frailty models are flexible and mathematically convenient like the popular gamma frailty model. A piecewise-exponential baseline hazard function is employed, yielding flexibility for the proposed class. Although a closed-form observed log-likelihood function is available, simulation studies show that employing an EM-algorithm is advantageous concerning the direct maximization of this function. Further simulated results address the comparison of different methods for obtaining standard errors of the estimates and confidence intervals for the parameters. Additionally, the finite-sample behavior of the EM-estimators is investigated and the performance of the GIG models under misspecification assessed. We apply our methodology to a TARGET (Therapeutically Applicable Research to Generate Effective Treatments) data about the survival time of patients with neuroblastoma cancer and show some advantages of the GIG frailties over existing models in the literature.

Appendix

1.1 A.1 Observed information matrix

Define $c_i(\theta )=\displaystyle \sum _{j=1}^{n_i}H_0(t_{ij})e^{x_{ij}^\top \beta }$ and $a_{irs}(\theta )=\displaystyle \sum _{j=1}^{n_i}H_0(t_{ij})e^{x_{ij}^\top \beta }x_{ijr}x_{ijs}$, for $r,s=1,\ldots ,p$, $\varDelta _i\equiv \varDelta _i(\theta )=\bigg \{\alpha ^{-1}\bigg (\alpha ^{-1}+2\displaystyle \sum _{j=1}^{n_i}H_0(t_{ij})e^{x_{ij}^\top \beta }\bigg )\bigg \}^{1/2}$, and $\lambda _i^*=\lambda +\sum _{j=1}^{n_i}\delta _{ij}$, for $i=1,\ldots ,m$. We have that the elements of the observed information matrix $J_n(\theta _*)=-\partial ^2\ell (\theta )/\partial \theta _*\partial \theta _*^\top$ are given by

$$\begin{aligned} \dfrac{\partial ^2\ell (\theta )}{\partial \beta _r\beta _s}= & {} \dfrac{1}{\alpha ^2}\sum _{i=1}^m\dfrac{K''_{\lambda _i^*}(\varDelta _i)K_{\lambda _i^*}(\varDelta _i)\varDelta _i-K'_{\lambda _i^*}(\varDelta _i)[K'_{\lambda _i^*}(\varDelta _i)\varDelta _i+K_{\lambda _i^*}(\varDelta _i)]}{K_{\lambda _i^*}(\varDelta _i)^2\varDelta _i^3}a_{ir}(\theta )a_{is}(\theta )\\&+\dfrac{1}{\alpha }\sum _{i=1}^m\dfrac{K'_{\lambda _i^*}(\varDelta _i)}{K_{\lambda _i^*}(\varDelta _i)\varDelta _i}a_{irs}(\theta ) -\sum _{i=1}^m\dfrac{{\lambda _i^*}a_{irs}(\theta )}{\alpha ^{-1}+2c_i(\theta )}+2\sum _{i=1}^m\dfrac{{\lambda _i^*}a_{ir}(\theta )a_{is}(\theta )}{(\alpha ^{-1}+2c_i(\theta ))^2}, \\ \dfrac{\partial ^2\ell (\theta )}{\partial \beta _r\partial \alpha }= & {} -\dfrac{1}{\alpha ^2}\sum _{i=1}^m\dfrac{K'_{\lambda _i^*}(\varDelta _i)}{K_{\lambda _i^*}(\varDelta _i)\varDelta _i}a_{ir}(\theta )-\dfrac{1}{\alpha ^2}\sum _{i=1}^m\frac{{\lambda _i^*}a_{ir}(\theta )}{(\alpha ^{-1}+2c_i(\theta ))^{2}}\\&-\dfrac{1}{\alpha ^3}\sum _{i=1}^m\dfrac{K''_{\lambda _i^*}(\varDelta _i)K_{\lambda _i^*}(\varDelta _i)\varDelta _i-K'_{\lambda _i^*}(\varDelta _i)[K'_{\lambda _i^*}(\varDelta _i)\varDelta _i+K_{\lambda _i^*}(\varDelta _i)]}{K_{\lambda _i^*}(\varDelta _i)^2\varDelta ^3_i}\\&\quad (\alpha ^{-1}+c_i(\theta ))a_{ir}(\theta ), \end{aligned}$$

for $r,s=1,\ldots ,p$, and

$$\begin{aligned} \dfrac{\partial ^2\ell (\theta )}{\partial \alpha ^2}= & {} \dfrac{1}{2\alpha ^2}\sum _{i=1}^m{\lambda _i^*}-\dfrac{2m}{\alpha ^3}\dfrac{K'_\lambda (\alpha ^{-1})}{K_\lambda (\alpha ^{-1})}-\dfrac{m}{\alpha ^4}\dfrac{K''_\lambda (\alpha ^{-1})K_\lambda (\alpha ^{-1})-K'_\lambda (\alpha ^{-1})^2}{K_\lambda (\alpha ^{-1})^2}\\&+ \dfrac{2}{\alpha ^3}\sum _{i=1}^m\dfrac{K'_{\lambda _i^*}(\varDelta _i)}{K_{\lambda _i^*}(\varDelta _i)\varDelta _i}(\alpha ^{-1}+c_i(\theta ))+ \dfrac{1}{\alpha ^4}\sum _{i=1}^m\dfrac{K'_{\lambda _i^*}(\varDelta _i)}{K_{\lambda _i^*}(\varDelta _i)\varDelta _i}\\&-\dfrac{1}{\alpha ^3}\sum _{i=1}^m\dfrac{{\lambda _i^*}}{\alpha ^{-1}+2c_i(\theta )}+\dfrac{1}{2\alpha ^4}\sum _{i=1}^m\dfrac{{\lambda _i^*}}{(\alpha ^{-1}+2c_i(\theta ))^2}\\&+\dfrac{1}{\alpha ^4}\sum _{i=1}^m\dfrac{K''_{\lambda _i^*}(\varDelta _i)K_{\lambda _i^*}(\varDelta _i)\varDelta _i-K'_{\lambda _i^*}(\varDelta _i)[K'_{\lambda _i^*}(\varDelta _i)\varDelta _i+K_{\lambda _i^*}(\varDelta _i)]}{K_{\lambda _i^*}(\varDelta _i)^2\varDelta ^3_i}(\alpha ^{-1}+c_i(\theta ))^2. \end{aligned}$$

1.2 A.2 Louis information matrix

From Louis (1982), we have that the information matrix obtained from the EM-algorithm, say $\mathbf{I}_n(\theta _*)$, is given by

$$\begin{aligned} I_n(\theta _*)=E\left( -\dfrac{\partial ^2\ell _c(\theta )}{\partial \theta _*\partial {\theta _*^\top }}\Big | Y^{obs}\right) -E\left( \dfrac{\partial \ell _c(\theta )}{\partial \theta _*}\dfrac{\partial \ell _c(\theta )}{\partial \theta _*}^\top \Big |Y^{obs}\right) , \end{aligned}$$

(10)

where we have defined $Y^{obs}=\{(t_{ij},\delta _{ij}),\, j=1,\ldots ,n_i,\, i=1,\ldots ,m\}$.

Let $\tau _i(\theta )=E(Z_i^2|Y^{obs})$ and $\nu _i(\theta )=E(Z_i^{-2}|Y^{obs})$ for $i=1,\ldots ,m$, where explicit expressions are directly available by using (3) and (9). The elements of the information matrix (10) are given by

$$\begin{aligned} E\left( -\dfrac{\partial ^2\ell _c(\theta )}{\partial \beta _r\partial \beta _s}\bigg |Y^{obs}\right)= & {} \sum _{i=1}^m\sum _{j=1}^{n_i}\omega _i(\theta )H_0(t_{ij})e^{x_{ij}^\top \beta }x_{ijr}x_{ijs}, \quad r,s=1,\ldots , p, \\ E\left( -\dfrac{\partial ^2\ell _c(\theta )}{\partial \alpha ^2}\bigg |Y^{obs}\right)= & {} \dfrac{m}{\alpha ^4}\dfrac{K_\lambda ''(\alpha ^{-1})K_\lambda (\alpha ^{-1})-K_\lambda '(\alpha ^{-1})^2}{K_\lambda (\alpha ^{-1})^2}\\&+\dfrac{2m}{\alpha ^3}\dfrac{K_\lambda '(\alpha ^{-1})}{K_\lambda (\alpha ^{-1})}+\dfrac{1}{\alpha ^3}\sum _{i=1}^m\left( \omega _i(\theta )+\kappa _i(\theta )\right) , \\ E\left( \dfrac{\partial \ell _c(\theta )}{\partial \beta _r}\dfrac{\partial \ell _c(\theta )}{\partial \beta _s}\bigg |Y^{obs}\right)= & {} \left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijr}\right) \left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijs}\right) \\&- \left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijs}\right) \left( \sum _{i=1}^m\omega _i(\theta )a_{ir}(\theta )\right) \\&-\left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijr}\right) \left( \sum _{i=1}^m\omega _i(\theta )a_{is}(\theta )\right) + \sum _{i=1}^m\tau _i(\theta )a_{ir}(\theta )a_{is}(\theta )\\&+ \sum _{i\ne i'}^m\omega _i(\theta )\omega _{i'}(\theta )a_{ir}(\theta )a_{i's}(\theta ),\quad r,s=1,\ldots , p, \\ E\left( \left( \dfrac{\partial \ell _c(\theta )}{\partial \alpha }\right) ^2\bigg |Y^{obs}\right)= & {} \dfrac{m^2}{\alpha ^4}\left( \dfrac{K'_\lambda (\alpha ^{-1})}{K_\lambda (\alpha ^{-1})}\right) ^2+ \dfrac{m}{\alpha ^4}\dfrac{K'_\lambda (\alpha ^{-1})}{K_\lambda (\alpha ^{-1})}\sum _{i=1}^m\left( \omega _i(\theta )+\kappa _i(\theta )\right) \\&+\dfrac{1}{4\alpha ^4}\sum _{i=1}^m\left( \tau _i(\theta )+2+\nu _i(\theta )\right) \\&+\dfrac{1}{4\alpha ^4}\sum _{i\ne i'}^m\left( \omega _i(\theta )+\kappa _i(\theta )\right) \left( \omega _{i'}(\theta )+\kappa _{i'}(\theta )\right) , \\ E\left( \dfrac{\partial \ell _c(\theta )}{\partial \beta _r}\dfrac{\partial \ell _c(\theta )}{\partial \alpha }\bigg |Y^{obs}\right)= & {} \dfrac{m}{\alpha ^2}\dfrac{K'_\lambda (\alpha ^{-1})}{K_\lambda (\alpha ^{-1})} \left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijr}-\sum _{i=1}^m\omega _i(\theta )a_{ir}(\theta )\right) \\&+\dfrac{1}{2\alpha ^2}\left( \sum _{i=1}^m\sum _{j=1}^{n_i}\delta _{ij}x_{ijr}\right) \\&\times \left( \sum _{i=1}^m\left( \omega _i(\theta )+\kappa _i(\theta )\right) \right) - \dfrac{1}{2\alpha ^2}\sum _{i=1}^m\left( \tau _i(\theta )+1\right) a_{ir}(\theta )\\&-\dfrac{1}{2\alpha ^2}\sum _{i\ne i'}^m\left( \omega _i(\theta )+\kappa _i(\theta )\right) \omega _{i'}(\theta )a_{i'r}(\theta ), \, r=1,\ldots ,p, \end{aligned}$$

and $E\left( -\dfrac{\partial ^2\ell _c(\theta )}{\partial \beta _r\partial \alpha }\bigg |Y^{obs}\right) =0$, for $r=1,\ldots ,p$, where we have defined $a_{ir}(\theta )=\displaystyle \sum _{j=1}^{n_i}H_0(t_{ij})e^{x_{ij}^\top \beta }x_{ijr}$, $b_{is}(\theta )=\displaystyle \sum _{j=1}^{n_i}H^{(s)}_0(t_{ij})e^{x_{ij}^\top \beta }$, for $i=1,\ldots ,m$, $r=1,\dots ,p$, and $s=1,\ldots ,k+1$.