(AST405) Lifetime data analysis
Location-scale distributions have survivor function of the form \[ S(y; u, b) = S_0\bigg(\frac{y - u}{b}\bigg) \quad -\infty < y < \infty \tag{5.1}\] \[-\infty < u < \infty\;\; \text{and} \;\; \;b >0\]
Log-lifetime \(Y=\log T\) has a location-scale distribution with survivor function of the form Equation 5.1
Lifetime variable \(T\) has a log-location-scale distribution with the survivor function \[\begin{align} S_T(t;\alpha, \beta) &= S_0\bigg(\frac{\log t - u}{b}\bigg) \notag \\[.25em] & = S_0^\star\Big[(t/\alpha)^\beta\Big] \end{align}\]
\(S_0^\star(w) = S_0(\log w)\) for \(w>0\)
\(u = \log \alpha\)
\(b = (1/\beta)\)
Lifetime and log-lifetime distributions
exponential, Weibull, log-logistic, log-normal, etc.
extreme-value, logistic, normal, etc.
Goal is to estimate the parameters \((u, b)\) or \((\alpha, \beta)\)
Some advantages of estimating \((u, b)\)
Log-likelihood function for \((u, b)\) is more closer to quadratic than that for \((\alpha, \beta)\)
Large sample normal approximations for \((\hat u, \hat b)\) tend to be more accurate that those for \((\hat \alpha, \hat \beta)\)
A better choice of parameters for obtaining MLEs and implementing normal approximations is \((u, \log b)\), which is used by most statistical software
For a censored sample \[\big\{(t_i, \delta_i), \;i=1, \ldots, n\big\},\] the likelihood function \[\begin{align} L(u, b) = \prod_{i=1}^n \Bigg[\frac{1}{b}\,f_0\bigg(\frac{y_i - u}{b}\bigg)\Bigg]^{\delta_i}\;\Bigg[S_0\bigg(\frac{y_i - u}{b}\bigg)\Bigg]^{1-\delta_i} \end{align}\]
\(y_i = \log t_i\)
\(f_0(z) = -d S_0(z)/dz\;\rightarrow\) pdf
The standardized variable \[z_i=\frac{y_i - u}{b}\]
The likelihood function \[\begin{align} L(u, b) = \prod_{i=1}^n \Bigg[\frac{1}{b}\,f_0\big(z_i)\Bigg]^{\delta_i}\;\Big[S_0\big(z_i)\Big]^{1-\delta_i} \end{align}\]
The corresponding log-likelihood function \[\begin{align} \ell(u, b) & = -r\log b +\sum_{i=1}^n\Big[\delta_i\log f_0(z_i) + (1-\delta_i)\log S_0(z_i) \Big] \notag \\[.25em] & = -r\log b + \sum_{i=1}^n\ell_i(z_i, \delta_i) \end{align}\]
\[\begin{aligned} {\color{purple} \ell(u, b)}& = -r\log b + \sum_{i=1}^n\ell_i(z_i, \delta_i) \\& = {\color{purple} -r\log b +\sum_{i=1}^n\Big[\delta_i\log f_0(z_i) + (1-\delta_i)\log S_0(z_i)\Big]}\end{aligned}\]
\[\begin{align} \frac{\partial \ell(u, b)}{\partial u} &= \sum_{i=1}^n\frac{\partial \ell_i(z_i, \delta_i)}{\partial z_i}\times \frac{\partial z_i}{\partial u}\notag \\ &= \sum_{i=1}^n\frac{\partial \ell_i(z_i, \delta_i)}{\partial z_i}\times \bigg(\frac{-1}{b}\bigg)\notag\\[.5em] & = -\frac{1}{b} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\notag \end{align}\]
\[ {\color{purple} \ell(u, b) = -r\log b +\sum_{i=1}^n\Big[\delta_i\log f_0(z_i) + (1-\delta_i)\log S_0(z_i) \Big]} \]
\[\begin{align} \frac{\partial \ell(u, b)}{\partial b} &= -\frac{r}{b} + \sum_{i=1}^n\frac{\partial \ell_i(z_i, \delta_i)}{\partial z_i}\times \frac{\partial z_i}{\partial b}\notag \\ &= -\frac{r}{b} + \sum_{i=1}^n \frac{\partial \ell_i(z_i, \delta_i)}{\partial z_i}\times \Bigg(\frac{-z_i}{b}\Bigg)\notag\\[.5em] & = -\frac{r}{b} -\frac{1}{b} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\notag \end{align}\]
\[ {\color{purple} \frac{\partial \ell(u, b)}{\partial u} = -\frac{1}{b} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]} \]
\[\begin{align*} \frac{\partial^2 \ell(u, b)}{\partial u^2} & = \frac{\partial}{\partial u} \bigg[\frac{\partial \ell(u, b)}{\partial u}\bigg]\\[.35em] %& = \frac{\partial}{\partial %z_i}\bigg[\frac{\partial \ell(u, b)}{\partial u}\bigg]\; %\frac{\partial z_i}{\partial u}\notag \\[.35em] & = \frac{1}{b^2} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg] \end{align*}\]
\[ {\color{purple}\frac{\partial \ell(u, b)}{\partial b} = -\frac{r}{b} -\frac{1}{b} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]} \] \[\begin{align*} \frac{\partial^2 \ell(u, b)}{\partial b^2} & = \frac{r}{b^2} +\frac{2}{b^2} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\\[.25em] & \;\;\;+\frac{1}{b^2}\sum_{i=1}^nz_i^2\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg] \end{align*}\]
\[\begin{align*} {\color{purple} \frac{\partial \ell(u, b)}{\partial u}} &= {\color{purple}-\frac{1}{b} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]} \\[.5em] \frac{\partial^2 \ell(u, b)}{\partial u\,\partial b} & = \frac{1}{b^2} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\\[.25em] &\;\; \;\;+ \frac{1}{b^2} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg] \end{align*}\]
\[ \begin{aligned} U(u, b) &= \begin{bmatrix} \frac{\partial \ell(u, b)}{\partial u} \\[.5em] \frac{\partial \ell(u, b)}{\partial b} \end{bmatrix} \\[1em] I(u, b) & = -H(u, b) = -\begin{bmatrix} \frac{\partial^2 \ell(u, b)}{\partial u^2} & \frac{\partial^2 \ell(u, b)}{\partial u\,\partial b} \\[.25em] \frac{\partial^2 \ell(u, b)}{\partial b\,\partial u} & \frac{\partial^2 \ell(u, b)}{\partial b^2} \end{bmatrix} \end{aligned} \]
\[ (\hat u, \hat b)' = {\arg\,\max}_{(u, b)' \in \Theta} \ell(u, b) \]
\[ \text{var}(\hat u, \hat b) = \Big[I(\hat u, \hat b)\Big]^{-1} = \hat{V} \]
\[ \begin{pmatrix}\hat u \\ \hat b\end{pmatrix} \sim \mathcal{N}_2\begin{pmatrix} \begin{bmatrix} u \\ b\end{bmatrix}, \Big[I(\hat u, \hat b)\Big]^{-1} \end{pmatrix} \]
For a large \(n\), \((\hat u, \hat b)'\) follows a bivariate normal distribution with mean \((u, b)'\) and variance matrix \(\hat V\)
Standard error of \(\hat u\) and \(\hat b\) can be obtained from the diagonal elements of \(\hat V\) \[ se(\hat u) = \hat{V}_{11}^{1/2}\;\;\text{and}\;\;se(\hat b) = \hat{V}_{22}^{1/2} \]
Following pivotal quantities follow standard normal distributions \[\begin{align*} Z_1 = \frac{\hat u - u}{se(\hat u)},\;\; \;\;\; Z_2 = \frac{\hat b - b}{se(\hat b)}, \;\;\;\;\; Z_2' = \frac{\log \hat b - \log b}{se(\log \hat b)} \end{align*}\]
\((1-p)100\%\) confidence intervals \[\begin{align*} \hat u &\pm z_{1-p/2}\,se(\hat u) \\ \hat b &\pm z_{1-p/2}\,se(\hat b) \\ \hat b &\exp\{\pm z_{1-p/2}\,se(\log \hat b)\} \end{align*}\]
\(p\)th quantile of log lifetime \(Y\) \[ \begin{aligned} P(Y\leq y_p) = p \;\Rightarrow\;S_0\Big(\frac{y_p - u}{b}\Big) &= 1- p\\ \frac{y_p - u}{b} &= S_0^{-1}(1-p)\\ y_p &= u + bw_p \end{aligned} \]
Estimate of \(p\)th quantile and the corresponding standard error \[\begin{aligned} \hat y_p &= \hat u + \hat b w_p\\[.3em] se(\hat y_p) &= \sqrt{\hat V_{11} + w_p^2 \hat V_{22} +2 w_p\hat V_{12}} \end{aligned}\]
Pivotal quantity \[ Z_p = \frac{\hat y_p - y_p}{se(\hat y_p)} \sim \mathcal{N}(0, 1) \]
\((1-q)100\%\) confidence interval for \(y_p\) \[ \hat y_p \pm z_{1-q/2}\, se(\hat y_p) \]
Normal approximation based confidence intervals could be inaccurate for small samples
An alternative to normal approximation, bootstrap simulations can be used to estimate the distributions of pivots
All these methods can perform poorly in small samples with heavy censoring
Implementation of likelihood ratio based confidence intervals is relatively complicated, but LRT based CI often performs better in small and medium-size samples
To test the hypothesis \(H_0: u = u_0\), the following likelihood ratio test statistic can be used \[ \Lambda_1(u_0) = 2\ell(\hat u, \hat b) - 2\ell(u_0, \tilde b(u_0)) \]
MLEs \[ \begin{aligned} (\hat u, \hat b)' & = {\arg\,\max}_{(u, b)' \in \Theta}\, \ell(u, b) \;\;\text{unrestricted} \\[.35em] \tilde b(u_0) & = {\arg\,\max}_{b \in \Theta_1} \,\ell(u_0, b)\;\;\text{under $H_0$} \end{aligned} \]
Under \(H_0: u = u_0\), asymptotically \(\Lambda_1(u_0)\sim\chi^2_{(1)}\)
The \(p\)th quantile of location-scale distribution can be expressed as \[ y_p = u + w_pb, \;\;where\;\;w_p = S_0^{-1}(1-p) \]
To obtain confidence intervals for a quantile, consider the null hypothesis \[H_0: y_p = y_{p_0}\]
The corresponding likelihood ratio test statistic \[ \Lambda(y_{p_0}) = 2\ell(\hat u, \hat b) - 2\ell(\tilde u, \tilde b) \tag{5.2}\]
The estimates \(\hat u\) and \(\hat b\) are MLEs under \(H_1\) \[(\hat u, \hat b)' = {\arg\,\max}_{(u, b)'\in\Theta} \ell(u, b)\]
Steps to obtain MLEs \(\tilde u\) and \(\tilde b\), under \(H_0: y_p = y_{p_0}\)
Under \(H_0\), \(y_{p_0} = u + w_pb\;\Rightarrow\;{\color{purple}u = y_{p_0} - bw_p}\)
\(\tilde b = {\arg\,\max}_{b\in\Theta} \ell(y_{p_0} - w_pb, b)\)
\(\tilde u = y_{p_0} - w_p \tilde b\)
\((1-q)100\%\) Confidence interval for \(y_p\) can be obtained from the set of \(y_{p_0}\) values such that \[ \Lambda(y_{p_0}) \leq \chi^2_{(1), 1- q} \]
To obtain confidence interval for \(S(y_0)\), consider the null hypothesis \[H_0: S(y_0) = s_0\]
The same likelihood ratio statistic Equation 5.2 can be used to test the hypothesis \[ \Lambda({s_0}) = 2\ell(\hat u, \hat b) - 2\ell(\tilde u, \tilde b) \]
Steps for obtaining MLEs \(\tilde u\) and \(\tilde b\) under \(H_0\)
Under \(H_0: S(y_0)=s_0\) \[ S(y_0) = S_0\Big(\frac{y_0 - u}{b}\Big) = s_0 \;\;\Rightarrow\;\;{\color{purple}u = y_0 - S_0^{-1}(s_0)b} \]
\(\tilde b ={\arg\,\max}_{b \in \Theta} \;\ell\big(y_0 - S_0^{-1}(s_0)b, b\big)\)
\(\tilde u = y_0 - S_0^{-1}(s_0)\tilde b\)
The \((1-p)100\%\) confidence interval for \(S(y_0)\) can be defined as the set of \(s_0\) values such that \[\begin{align} \Lambda(s_0) \leq \chi^2_{(1), 1-p} \end{align}\]
The likelihood ratio procedure can provide quite accurate confidence intervals when the number of failures is about 20 or more
Two-sided intervals perform better than one-sided intervals as the former giving more closer to nominal coverage than the other
The pdf of Weibull distribution \[\begin{align} f(t; \alpha, \beta) = \frac{\beta}{\alpha}\;\bigg(\frac{t}{\alpha}\bigg)^{\beta-1} \exp\Big[-\big(t/\alpha\big)^\beta\Big] \end{align}\]
The pdf of extreme-value distribution \[\begin{align} f(y; u, b) &= \frac{1}{b}\;\exp\big[(y-u)/b\big]\;\exp\Big[-e^{(y-u)/b}\Big] \\ & = \frac{1}{b}\,f_0\Big(\frac{y-u}{b}\Big) \end{align}\]
\(u=\log\alpha\)
\(b=(1/\beta)\)
Extreme-value distribution is used to make inferences about Weibull distribution
Censored sample \[ \big\{(t_i, \delta_i),\;i=1, \ldots, n\big\} \]
Define \[y_i =\log t_i \;\;\text{and}\;\; z_i = (y_i - u)/b\]
General expression of likelihood function \[ {\color{purple} \ell(u, b) = -r\log b +\sum_{i=1}^n\Big[\delta_i\log f_0(z_i) + (1-\delta_i)\log S_0(z_i) \Big]} \]
For extreme-value distribution \[\begin{align*} S_0(z) &= \exp(-e^z)\\[.25em] f_0(z) &= -\frac{d}{dz}S_0(z)=\exp(z-e^z) \end{align*}\]
\[ {\color{purple} \ell(u, b) = -r\log b +\sum_{i=1}^n\Big[\delta_i\log f_0(z_i) + (1-\delta_i)\log S_0(z_i) \Big]} \]
The general expression for location-scale family can also help us find the expressions for the first (and also second) derivatives of \(\ell(u, b)\).
General expressions \[\begin{align} \frac{\partial \ell(u, b)}{\partial u} & = -\frac{1}{b} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\notag \\[.25em] \frac{\partial \ell(u, b)}{\partial b} & = -\frac{r}{b} -\frac{1}{b} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\notag \end{align}\]
For extreme-value distribution \[\begin{align*} \frac{\partial \log f_0(z)}{\partial z} & = \frac{\partial}{\partial z} \log\big\{\exp(z-e^z)\big\} =1 - e^z \\[.25em] \frac{\partial \log S_0(z)}{\partial z} & = \frac{\partial}{\partial z} \log\big\{\exp(-e^z)\big\} =- e^z \end{align*}\]
These gives straightforward expressions for the first and second derivatives of \(\ell(u, b)\), which can be used to find MLEs, \({\hat u}, {\hat b}\)
\[\begin{align*} \frac{\partial^2 \ell(u, b)}{\partial u^2} & = \frac{1}{b^2} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg] =-\frac{1}{b^2} \sum_{i=1}^n e^{z_i} \end{align*}\] where \[\begin{align*} \frac{\partial^2 \log f_0(z)}{\partial z^2} & = \frac{\partial}{\partial z} \big\{1 - e^z\big\} = -e^z= \frac{\partial^2 \log S_0(z)}{\partial z^2} \end{align*}\]
\[\begin{align} {\color{purple} \frac{\partial^2 \ell(u, b)}{\partial b^2} \bigg\vert_{u=\hat u, \,b=\hat b}=-\frac{1}{\hat b^2} \sum_{i=1}^n e^{\hat z_i}=-\frac{r}{\hat b^2}} \end{align}\]
From Section 5.1.4, it can be shown that the Hessian matrix will be
\[\begin{align*} \frac{\partial^2 \ell(u, b)}{\partial b^2} & = \frac{r}{b^2} +\frac{2}{b^2} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\\[.25em] & \;\;\;+\frac{1}{b^2}\sum_{i=1}^nz_i^2\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg]\\[.25em] & = \frac{r}{b^2} +\frac{2}{b^2} \sum_{i=1}^nz_i \big[\delta_i(1-e^{z_i}) -(1-\delta_i)e^{z_i}\big] -\frac{1}{b^2}\sum_{i=1}^nz_i^2 e^{z_i} \\[.25em] &=\frac{r}{b^2} +\frac{2}{b^2} \sum_{i=1}^nz_i\big[\delta_i-e^{z_i}\big] -\frac{1}{b^2}\sum_{i=1}^nz_i^2 e^{z_i} \end{align*}\]
\[\begin{align} {\color{purple} \frac{\partial^2 \ell(u, b)}{\partial b^2} \bigg\vert_{u=\hat u,\,b=\hat b}= -\frac{r}{\hat b^2}-\frac{1}{\hat b^2}\sum_{i=1}^n\hat z_i^2 e^{\hat z_i}} \end{align}\]
\[\begin{align*} \frac{\partial^2 \ell(u, b)}{\partial u\,\partial b} & = \frac{1}{b^2} \sum_{i=1}^n\bigg[\delta_i\,\frac{\partial \log f_0(z_i)}{\partial z_i} + (1-\delta_i)\, \frac{\partial \log S_0(z_i)}{\partial z_i}\bigg]\\[.25em] &\;\; \;\;+ \frac{1}{b^2} \sum_{i=1}^nz_i\bigg[\delta_i\,\frac{\partial^2 \log f_0(z_i)}{\partial z_i^2} + (1-\delta_i)\, \frac{\partial^2 \log S_0(z_i)}{\partial z_i^2}\bigg] \\[.25em] & = \frac{1}{b^2} \sum_{i=1}^n\Big[\delta_i(1-e^{z_i}) - (1-\delta_i)e^{z_i}-z_ie^{z_i}\Big]\\[.25em] & = \frac{1}{b^2} \sum_{i=1}^n\Big[\delta_i - e^{z_i}-z_ie^{z_i}\Big] = \frac{1}{b^2} \Big[r -\sum_{i=1}^n e^{z_i}-\sum_{i=1}^n z_ie^{z_i} \Big] \end{align*}\]
\[\begin{align} {\color{purple} \frac{\partial^2 \ell(u, b)}{\partial u\,\partial b} \bigg\vert_{u=\hat u,\,b=\hat b}= -\frac{1}{\hat b^2}\sum_{i=1}^n\hat z_i e^{\hat z_i}} \end{align}\]
Observed information matrix at MLEs \(\hat u\) and \(\hat b\) \[\begin{align*} I(\hat u, \hat b) &= -H(\hat u, \hat b) \\[.25em] & = \frac{1}{{\hat{b}}^2}\begin{bmatrix} r & \sum_{i=1}^n \hat z_i e^{\hat z_i}\\[.5em] \sum_{i=1}^n \hat z_i e^{\hat z_i} & r + \sum_{i=1}^n \hat z_i^2 e^{\hat z_i} \end{bmatrix} \end{align*}\]
Covariance matrix of \((\hat u, \hat b)'\) \[\begin{align} \hat V & = \Big[I(\hat u, \hat b)\Big]^{-1} \end{align}\]
MLEs of \(\alpha\) and \(\beta\) (Weibull model parameters) \[ \hat\alpha = e^{\hat u}\;\;\text{and}\;\;\hat\beta=1/\hat b \]
Covariance matrix of \((\hat\alpha, \hat\beta)'\) (using multivariate delta method) \[ {\color{purple} \operatorname{var}(\hat\alpha, \hat\beta) = \boldsymbol{G}\,\hat V\,\boldsymbol{G}'} \] where \[ G = \begin{bmatrix} \frac{\partial g_1(u, b)}{du} & \frac{\partial g_1(u, b)}{\partial b} \\[.25em] \frac{\partial g_2(u, b)}{\partial u} & \frac{\partial g_2(u, b)}{\partial b} \end{bmatrix} = \begin{bmatrix}e^{\hat u} & 0 \\[.25em] 0 & -\frac{1}{\hat b^2} \end{bmatrix} \]
\(\alpha = g_1(u, b) = e^u\)
\(\beta = g_2(u, b)=(1/b)\)
LRT statistics \[ \Lambda(b_0) = 2\ell(\hat u, \hat b) -2\ell(\tilde u(b_0), b_0) \]
\(100(1-p)\%\) CI for \(b\) is defined by the set of \(b_0\) values such that \[ \Lambda_1(b_0)\leq \chi^2_{(1), 1-p} \]
Similarly, confidence interval for \(u\) can be obtained using the corresponding LRT statistics (Homework)
The estimate of \(p\)th quantile \[ \hat{y}_p = \hat{u} + w_p \hat{b} \]
Standard error of \(\hat y_p\) (using the multivariate delta method) \[ \operatorname{var}(\hat y_p) = \begin{bmatrix}1 & w_p\end{bmatrix} \hat V \begin{bmatrix} 1 \\ w_p\end{bmatrix} = \hat V_{11} + \hat V_{22}w_p^2 + 2\hat V_{12}w_p \]
Large sample based \(100(1-q)\%\) confidence interval for \(y_p\) \[ \hat y_p\pm z_{1-q/2} \;se(\hat y_p) \]
Find the \(100(1-q)\%\) confidence interval for \(t_p\)
To obtain LRT statistic based confidence interval for the quantile \(y_p\), consider the following null hypothesis \[ H_0: y_p = y_{p_0} \]
The corresponding LRT statistic \[ \Lambda(y_{p_0}) = 2\ell(\hat u, \hat b) - 2\ell(\tilde u, \tilde b) %\tag{\ref{eq-lrt-quantile}} \]
LRT statistic based \((1-q)100\%\) confidence interval for \(y_p\) can be obtained from the set of \(y_{p_0}\) values such that \[ \Lambda(y_{p_0}) \leq \chi^2_{(1), 1- q} \]
To obtain confidence interval for survival probability \[ S(y_0) = S_0\Big(\frac{y_0-u}{b}\Big)=\exp\Big[-\exp\Big(\frac{y_0-u}{b}\Big)\Big] \]
We can defined \[ \psi = S_0^{-1}\Big(S(y_0)\Big)=\log\big[-\log\big(S(y_0)\big)\big]=\frac{y_0 - u}{b} \]
MLE and SE \[ \begin{aligned} \hat \psi &= \frac{y_0 - \hat u}{\hat b}\\ \operatorname{var}(\hat \psi) & = \boldsymbol{a}'\hat V\boldsymbol{a} = \begin{bmatrix}-1/\hat b & - \hat\psi/b\end{bmatrix} \begin{bmatrix} \hat V_{11} & \hat V_{12} \\ \hat V_{21} & \hat V_{22}\end{bmatrix} \begin{bmatrix} -1/\hat b \\ -\hat\psi/b\end{bmatrix} \end{aligned} \]
\((1-p)100\%\) CI for \(\psi\) \[ \begin{aligned} \hat\psi - se(\hat\psi)\,z_{1-p2} &< \psi \leq \hat\psi + se(\hat\psi)\,z_{1-p/2} \\ L & < \psi \leq U \end{aligned} \]
Confidence interval for \(S(y_0)\) \[ \begin{aligned} L &< \log\big[-\log\big(S(y_0)\big)\big] \leq U \\ \exp\big[-\exp(U)\big] &<S(y_0) \leq \exp\big[-\exp(L)\big] \end{aligned} \]
Consider the null hypothesis \(H_0: S(y_0)=s_0\), where \[ S(y_0) = \exp\Big[-\exp\Big(\frac{y_0-u}{b}\Big)\Big] \]
The \((1-p)100\%\) confidence interval for \(S(y_0)\) can be defined as the set of \(s_0\) values such that \(\Lambda(s_0) \leq \chi^2_{(1), 1-p}\), where \[ \Lambda(s_0) = 2\ell(\hat u, \hat b) - 2\ell(\tilde u, \tilde b) %\tag{\ref{lrt-s0}} \]
Leukemia remission time data were given in Example 1.1.7 and used as an example for the non-parametric methods (e.g. Kaplan-Meier method) described in Chapter 3
Two groups of patients (6MP and placebo), each group has 21 patients, were followed up to observed either remission or censoring times (in weeks)
glimpse(gehan65)Rows: 42
Columns: 3
$ time <dbl> 6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 3…
$ status <dbl> 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, …
$ drug <chr> "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP",…| drug | status | n |
|---|---|---|
| 6-MP | 0 | 12 |
| 6-MP | 1 | 9 |
| placebo | 1 | 21 |
Two separate Weibull distributions are assumed for the failure times of two treatment groups, e.g.
6MP group: \[T\sim \text{Weibull}(\alpha_1, \beta_1),\; Y=\log T\sim EV(u_1, b_1)\]
Placebo group: \[T\sim \text{Weibull}(\alpha_2, \beta_2),\; Y=\log T\sim EV(u_2, b_2)\]
Objectives: Drawing inference about the parameters
Observed data \[ \big\{(t_i, \delta_i), i=1, \ldots, n\big\} \]
Log-likelihood function \[ \ell(\alpha, \beta) = \sum_{i=1}^n\Big[\delta_i\log f(t_i; \alpha, \beta) + (1-\delta_i)\log S(t_i; \alpha, \beta)\Big] \]
MLEs \[(\hat\alpha, \hat\beta)' = {\arg\,\max}_{(\alpha, \beta)'\in \Theta}\,\ell(\alpha, \beta)\]
Define \(y = \log t\) and corresponding probability density and survivor function \[\begin{align} f(y; u, b) &= \frac{1}{b}\,\exp\Big[(y-u)/b - e^{(y-u)/b}\Big] \\ S(y; u, b) &= \exp\Big[ - e^{(y-u)/b}\Big] \end{align}\]
Log-likelihood function \[\begin{align} \ell_{ev}(u, b) = \log \prod_{i=1}^n\Big[f(y_i; u, b)\Big]^{\delta_i}\big[S(y_i; u, b)\Big]^{1-\delta_i} \end{align}\]
MLEs \[(\hat u, \hat b)' = {\arg\,\max}_{(u, b)'\in \Theta}\,\ell_{ev}(u, b)\]
survreg functionsurvreg() can also be used to fit distributions of log-location-scale family, its syntax is similar to the syntax of survfit()
survreg(formula, data, dist)In formula, response is a Surv object, e.g. to model the variables time and status \[\texttt{formula = Surv(time, status)}\sim \texttt{1}\]
Lifetime or log-lifetime distributions can be passed to survreg by the argument dist
Available lifetime or log-lifetime distributions include “weibull”, “exponential”, “gaussian”, “logistic”, “lognormal”, “loglogistic”, “extreme”
The time argument of Surv function is either a lifetime or a log-lifetime depending on whether the mentioned dist is a lifetime (e.g. “weibull”) or a log-lifetime (e.g. “extreme”) \[\begin{align*}
\text{weibull}\rightarrow\;\; & \texttt{formula = Surv(time, status)}\sim \texttt{1}\\
\text{extreme}\rightarrow\;\; & \texttt{formula = Surv(log(time), status)}\sim \texttt{1}
\end{align*}\]
d6mp <- gehan65 |>
filter(drug == "6-MP")
glimpse(d6mp)Rows: 21
Columns: 3
$ time <dbl> 6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 3…
$ status <dbl> 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0
$ drug <chr> "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP", "6-MP",…survreg functionw_sreg_6mp <- survreg(Surv(time, status) ~ 1,
data = d6mp, dist = "weibull")ev_sreg_6mp <- survreg(Surv(log(time), status) ~ 1,
data = d6mp, dist = "extreme")broom::tidy(ev_sreg_6mp)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.52 0.273 12.9 6.28e-38
2 Log(scale) -0.303 0.278 -1.09 2.77e- 1vcov(ev_sreg_6mp) (Intercept) Log(scale)
(Intercept) 0.07473057 0.03305811
Log(scale) 0.03305811 0.07750538survreg functionThe survreg() function returns estimates of \((u, \log b)'\) and corresponding variance matrix
For making inference about Weibull distribution, followings are required
estimate of \((u, b)'\) and the corresponding variance matrix
estimate of \((\alpha, \beta)'\) and the corresponding variance matrix
It is important to understand the methods to obtain estimates and the corresponding variance of \((u, b)'\) and \((\alpha, \beta)'\) from the estimates and the corresponding variance of \((u, \log b)'\)
95% CI using the sampling distribution of \((\hat u, \log \hat b)\) \[\begin{align*} \hat u \pm z_{.975} \,se(\hat u) &= 2.984 \;\text{to}\; 4.055\\ \hat \alpha \pm z_{.975} \,se(\hat\alpha) &= \exp(2.984) \;\text{to}\; \exp(4.055) \\ &= 19.76 \;\text{to}\; 57.698 \end{align*}\]
(Obtain the variance matrix of \((\hat u, \hat b)\) using the sampling distribution of \((\hat u, \log\hat b)'\))
\[\begin{align*} \hat u \pm z_{.975} \,se(\hat u) &= 2.984 \;\text{to}\; 4.055\\ \hat \alpha \pm z_{.975} \,se(\hat\alpha) &= \exp(2.984) \;\text{to}\; \exp(4.055) \\ &= 19.76 \;\text{to}\; 57.698 \end{align*}\]
| parameter | method | 6-MP |
|---|---|---|
| \(\alpha\) | Wald \((\hat \alpha)\) | (15.674, 51.856) |
| NA | Wald \((\hat u)\) | (19.76, 57.698) |
| NA | LRT | (21.933, 76.708) |
| \(\beta\) | Wald \((\hat \beta)\) | (0.615, 2.092) |
| NA | Wald \((\log \hat b)\) | (0.784, 2.336) |
| NA | Wald \((\hat b)\) | (0.876, 2.979) |
| NA | LRT | (0.726, 2.203) |
dplacebo <- gehan65 %>%
filter(drug == "placebo")Model fit with the data of placebo group
w_sreg_p <- survreg(Surv(time, status) ~ 1,
data = dplacebo, dist = "weibull")Estimates of model parameters
broom::tidy(w_sreg_p)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.25 0.168 13.4 5.72e-41
2 Log(scale) -0.315 0.174 -1.82 6.94e- 2| parameter | method | 6-MP | Placebo |
|---|---|---|---|
| \(\alpha\) | Wald \((\hat \alpha)\) | (15.674, 51.856) | (6.363, 12.601) |
| NA | Wald \((\hat u)\) | (19.76, 57.698) | (6.824, 13.175) |
| NA | LRT | (21.933, 76.708) | (6.659, 13.25) |
| \(\beta\) | Wald \((\hat \beta)\) | (0.615, 2.092) | (0.904, 1.837) |
| NA | Wald \((\log \hat b)\) | (0.784, 2.336) | (0.975, 1.926) |
| NA | Wald \((\hat b)\) | (0.876, 2.979) | (1.023, 2.077) |
| NA | LRT | (0.726, 2.203) | (0.951, 1.868) |
Estimate of \(p\)th quantile \[ \hat y_p = \hat u + \hat b w_p \]
\(w_p = \log[-\log(1-p)]\)
\(\hat u= 3.519\) and \(\hat b = 0.739\) (for treatment group)
Wald-type CI (see Section 5.2.6.1 for detail)
\[
\hat y_p \pm se(\hat y_p) z_{1-q/2}
\]
Note the estimate of \(\hat y_p\) depends on the estimate of \(\hat u\) and \(\hat b\), and the corresponding variance matrix
survreg() returns estimate and variance matrix for \(\hat u\) and \(\log \hat b\)
| \(p\) | \(w_p\) | \(\hat y_p\) | \(se(\hat y_p)\) | lower | upper |
|---|---|---|---|---|---|
| 0.25 | -1.246 | 2.599 | 0.655 | 3.726 | 48.559 |
| 0.50 | -0.367 | 3.249 | 0.264 | 15.357 | 43.225 |
| 0.75 | 0.327 | 3.761 | 0.395 | 19.822 | 93.241 |
| \(p\) | lower | upper | lower | upper |
|---|---|---|---|---|
| 0.25 | 3.726 | 48.559 | 6.586 | 23.058 |
| 0.50 | 15.357 | 43.225 | 16.289 | 51.342 |
| 0.75 | 19.822 | 93.241 | 27.522 | 112.730 |
| \(p\) | lower | upper | lower | upper |
|---|---|---|---|---|
| 0.25 | 1.362 | 10.707 | 2.031 | 5.927 |
| 0.50 | 4.499 | 11.708 | 5.755 | 9.488 |
| 0.75 | 8.592 | 16.863 | 8.873 | 16.996 |
For Weibull distribution, the expression of survivor function \[ S(t; \alpha, \beta) = \exp\big(-(t/\alpha)^{\beta}) \]
Estimated survivor function \[ S(t; \hat\alpha, \hat\beta) = \exp\big(-(t/\hat\alpha)^{\hat\beta}) \]
| par | 6-MP | placebo |
|---|---|---|
| \(\alpha\) | 33.765 | 9.482 |
| \(\beta\) | 1.354 | 1.370 |
