On a convergent stochastic estimation algorithm for frailty models

Abstract

A maximum likelihood estimation procedure is presented for the frailty model. The procedure is based on a stochastic Expectation Maximization algorithm which converges quickly to the maximum likelihood estimate. The usual expectation step is replaced by a stochastic approximation of the complete log-likelihood using simulated values of unobserved frailties whereas the maximization step follows the same lines as those of the Expectation Maximization algorithm. The procedure allows to obtain at the same time estimations of the marginal likelihood and of the observed Fisher information matrix. Moreover, this stochastic Expectation Maximization algorithm requires less computation time. A wide variety of multivariate frailty models without any assumption on the covariance structure can be studied. To illustrate this procedure, a Gaussian frailty model with two frailty terms is introduced. The numerical results based on simulated data and on real bladder cancer data are more accurate than those obtained by using the Expectation Maximization Laplace algorithm and the Monte-Carlo Expectation Maximization one. Finally, since frailty models are used in many fields such as ecology, biology, economy, …, the proposed algorithm has a wide spectrum of applications.

Appendix: Additional assumptions and sketch of the proof

Roughly speaking to prove Theorem 1, it is sufficient to apply the main convergence result of Kuhn and Lavielle (2004). By this purpose, we have to verify that its assumptions are fulfilled, in particular assumptions (M1) and (M3) (see Kuhn and Lavielle 2004, p. 117).

Let us start by establishing assumption (M1). By assumption (F5), the probability density function f _η can be written as follows:

$$ f_{\eta}(b) = \exp \bigl(-\psi(\eta) + \bigl \langle \hat{S}(b) , \phi(\eta) \bigr \rangle \bigr),$$

(6)

where 〈⋅ ,⋅〉 denotes the usual scalar product, $\hat{S}$ is a Borel function on ℝ taking its values in an open subset, and ψ and ϕ are Borel functions on ℝ^b.

Combining (6) with assumptions (F1–F4), it implies that the considered frailty model belongs also to the exponential family. Consider the two vectors

and the number

$$\varPsi(\theta)=-N \ \psi(\eta)+\sum_{i=1}^{N}\sum_{j=1}^{n_i} \; \; \delta_{ij}\bigl(\log \lambda_0(y_{ij})+x_{ij}^t \beta \bigr).$$

The complete data likelihood function L _N given in Eq. (2) can be written as follows:

$$L_{N}(\mathbf{y},\boldsymbol{\delta},\mathbf{b};\theta) = \exp \bigl(-\varPsi(\theta) +\bigl \langle \tilde {S}(\mathbf{b}) , \varPhi (\theta) \bigr \rangle \bigr),$$

where $\tilde {S}$ is a Borel function on ℝ^N taking its values in an open subset $\mbox {$\mathcal {S}$}$ of $\mbox {$\mathbb {R}^{m}$}$, the value of m depending of model (1) and of the distribution of the frailty.

The assumption (M1) is therefore fulfilled.

Let us state the assumption (H1) which is the same as the assumption (M2) of Kuhn and Lavielle (2004).

H1: Sufficient statistic assumptions.:

The convex hull of $\tilde {S}(\mbox {$\mathbb {R}^{N}$})$ is included in $\mbox {$\mathcal {S}$}$. The function $\bar{s} : \varTheta\rightarrow \mbox {$\mathcal {S}$}$ defined as follows

$$\bar{s}(\theta) = \int_{\mbox {$\mathbb {R}^{N}$}} \ \tilde {S}(\mathbf{b}) \ \pi_\theta(\mathbf{b}|\mathbf{y},\boldsymbol{\delta}) d\mathbf{b}$$

exists and is continuously differentiable on Θ.

Let us turn to assumption (M3). First, let us introduce assumption (H2) which corresponds to its regularity part.

H2: Regularity assumptions.:

The function $\log L_{N}^{obs}: \varTheta\rightarrow \mbox {$\mathbb {R}^{}$}$ is continuously differentiable on Θ and satisfies

$$\partial_{\theta}\ \int\ L_{N}(\mathbf{y},\boldsymbol{\delta},\mathbf{b}; \theta) d\mathbf{b}= \int \partial_{\theta}\ L_{N}(\mathbf{y},\boldsymbol{\delta},\mathbf{b}; \theta) d\mathbf{b}.$$

Denote by $L : \mbox {$\mathcal {S}$}\times\varTheta\to \mbox {$\mathbb {R}^{}$}$ the function defined as follows:

$$L(s; \theta)= - \varPsi(\theta) + \bigl \langle s , \varPhi(\theta) \bigr \rangle .$$

There exists a continuously differentiable function ${\widehat {\theta }}: \mbox {$\mathcal {S}$}\rightarrow\varTheta$, such that:

$$\forall s \in \mbox {$\mathcal {S}$}, \ \forall\theta\in\varTheta,\quad L\bigl(s; {\widehat {\theta }}(s)\bigr)\geq L(s;\theta).$$

To fulfill assumption (M3), it suffices to show that Ψ and Φ are twice continuously differentiable. This is a straight consequence of assumptions (F1–F5).

Finally, let us state the assumptions (SAEM1–SAEM3) which are the same of those of Kuhn and Lavielle (2004).

SAEM1: Step-size sequence assumptions.:

For all k in $\mbox {$\mathbb {N}$}$, γ _k ∈[0,1], such that $\sum_{k=1}^{\infty}\gamma_{k} = \infty$ and $\sum_{k=1}^{\infty}\gamma_{k}^{2} <\infty$.

SAEM2: Additional regularity assumptions.:

The function $\log L_{N}^{obs}: \varTheta \rightarrow \mbox {$\mathbb {R}^{}$}$ and ${\widehat {\theta }}: \mbox {$\mathcal {S}$}\rightarrow\varTheta$ are m times continuously differentiable.

SAEM3: MCMC assumptions.:

• 1. The Markov chain (b _k ) _k≥0 takes its values in a compact subset $\mbox {$\mathcal {E}$}$ of $\mbox {$\mathbb {R}^{N}$}$.

• 2. For any compact subset V of Θ, there exists a real constant C such that for any (θ,θ′) in V ²

  $$\sup_{(a,b)\in\mathcal{E}^2} \bigl|\varPi_\theta(a,b)-\varPi_{\theta '}(a,b)\bigr|\leq C \bigl|\theta-\theta'\bigr|.$$

• 3. The transition probability Π _θ generates a uniformly ergodic chain whose invariant probability is the conditional distribution π _θ (⋅|y,δ):

  

  where ∥⋅∥ _TV denotes the total variation norm. We suppose that:

  $$\sup_\theta K_\theta< +\infty\quad\text{and}\quad \sup_\theta \rho_\theta< 1.$$

• 4. The function $\tilde {S}$ is bounded on $\mbox {$\mathcal {E}$}$.

Remark 7

The first part of the assumption (SAEM3) can be slacked off using truncation on random boundary. We refer to Allassonnière et al. (2010) for more information on this extension.

The proof of Theorem 1 is now complete.