Solution: Maximum likelihood estimators for a shared frailty model (Exam 2022)

The Laplace transform of a stochastic variable \(Z\) is defined as

\[ {\mathscr L}(c) = \mbox{E}[e^{-cZ}].\]

For a \(Z\) that are exponentially distributed with mean \(1/\lambda\) we then get

\[\begin{align*}{\mathscr L}(c) &= \int_0^\infty e^{-cz}\lambda e^{-\lambda z} dz\\ &= \lambda \int_0^\infty e^{-(\lambda + c)z} dz\\ &= \lambda \left[ - \frac{1}{\lambda+c}e^{-(\lambda+c)z}\right]_0^\infty\\ &= \lambda \cdot \left(- \frac{1}{\lambda+c}\cdot 0 - \left(-\frac{1}{\lambda+c}\cdot e^0\right)\right)\\ &= \underline{\underline{\frac{\lambda}{\lambda + c}}} \end{align*} \]

To show the given formula for \({\mathscr L}^{(r)}(c)\) we first check that it is correct for \(r=0\),

\[ {\mathscr L}^{(0)}(c) = (-1)^0 \frac{\lambda \cdot 0!}{(\lambda + c)^{0+1}} = \frac{\lambda}{\lambda + c}, \]

which is equal to \( {\mathscr L}(c) \) as it should. Then we do the induction step, i.e. we assume the given formula to be correct for \( r=s \) and use this to check the formula for \( r=s+1 \). As we now have assumed the formula to be correct for \(r=s\) we have that

\[ \begin{align*} {\mathscr L}^{(s)}(c) &= (-1)^s \frac{\lambda \cdot s!}{(\lambda+c)^{s+1}}\\ &= (-1)^s \lambda (s!) (\lambda + c)^{-(s+1)}. \end{align*} \]

Taking the derivative of this expression with respect to \(c\) we get

\[ \begin{align*} {\mathscr L}^{(s+1)}(c) &= (-1)^s \lambda (s!) (-(s+1)) (\lambda+c)^{-(s+1)-1}\\ &= (-1)^{s+1} \lambda ((s+1)!) (\lambda+c)^{-((s+1)+1)}\\ &= (-1)^{s+1} \frac{\lambda \cdot (s+1)!}{(\lambda+c)^{(s+1)+1}}. \end{align*} \]

We see that this last expression is identical to what we get by inserting \(r=s+1\) in the given expression for \({\mathscr L}^{(r)}(c)\), and we have thereby shown that the given formula is correct for \(r=0,1,2,\ldots\).

The situation given in the problem text is identical to the situation discussed in Section 7.2.2 in ABG, so the log-likelihood is given by (7.3) in ABG,

\[ \ell (\lambda,k) = \sum_{i=1}^m \left[ \sum_{j=1}^{n_i} D_{ij}\ln (\alpha(\widetilde{T}_{ij})) + \ln \left( (-1)^{D_{i\bullet}}{\mathscr L}^{(D_{i\bullet})}(V_i)\right)\right], \]

where

\[ V_i = \sum_{j=1}^{n_i} A(\widetilde{T}_{ij}).\]

From the given formula for \(\alpha(t|Z)\) we have that \(\alpha(t) = t^{k-1}\), which gives

\[ A(t) = \int_0^t \alpha(u) du = \int_0^t u^{k-1} du = \left[ \frac{u^k}{k}\right]_0^t = \frac{t^k}{k}, \]

which in turn gives

\[V_i = \sum_{j=1}^{n_i} \frac{\widetilde{T}^k}{k} = \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k.\]

Inserting the expressions we have for \(\alpha(t)\), \(V_i\) and \({\mathscr L}^{(r)}(c)\) into the expression for \(\ell(\lambda,k)\) we get

\[ \begin{align*} \ell (\lambda,k) &= \sum_{i=1}^m \left[ \sum_{j=1}^{n_i} D_{ij} \ln (\widetilde{T}_{ij}^{k-1}) + \ln \left( (-1)^{D_{i\bullet}}\cdot (-1)^{D_{i\bullet}} \frac{\lambda \cdot D_{i\bullet}!}{(\lambda + \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k)^{D_{i\bullet}+1}}\right)\right]\\ &= \sum_{i=1}^m \sum_{j=1}^{n_i} D_{ij}(k-1)\ln (\widetilde{T}_{ij}) + \sum_{i=1}^m \left[ \ln(\lambda)+\ln(D_{i\bullet}!) - (D_{i\bullet}+1)\ln \left(\lambda + \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k\right) \right]\\ &= (k-1)\sum_{i=1}^m\sum_{j=1}^{n_i} D_{ij}\ln(\widetilde{T}_{ij}) + m\ln(\lambda) + \sum_{i=1}^m\ln(D_{i\bullet}!) - \sum_{i=1}^m (D_{i\bullet}+1)\ln \left(\lambda + \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k\right). \end{align*}\]

To (try to) find expressions for the maximum likelihood estimators we need to find expressions for the partial derivatives of the log-likelihood function,

\[\begin{align*} \frac{\partial \ell}{\partial \lambda} &= \frac{m}{\lambda} - \sum_{i=1}^m \frac{D_{i\bullet}+1}{\lambda + \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k},\\ \frac{\partial \ell}{\partial k} &= \sum_{i=1}^m\sum_{j=1}^{n_i} D_{ij}\ln (\widetilde{T}_{ij}) - \sum_{i=1}^m (D_{i\bullet}+1) \cdot \frac{-\frac{1}{k^2}\sum_{j=1}^{n_i}\widetilde{T}^k_{ij} + \frac{1}{k}\sum_{j=1}^{n_i} \widetilde{T}^k_{ij}\ln(\widetilde{T}_{ij})}{\lambda + \frac{1}{k}\sum_{j=1}^{n_i}\widetilde{T}_{ij}^k}. \end{align*}\]

We see that in both partial derivatives both \(\lambda\) and \(k\) appears in the denominator of the fraction inside the sum over \(i\), so when equating them to zero we are not able to solve any of the two resulting equations with respect to any of the two parameters. To find the maximum likelihood estimates we therefore need to resort to numerical optimisation, for example the Newton-Raphson algorithm.