Chapter 5A

(AST405) Lifetime data analysis

5 Inference Procedures for Log-location-scale Distributions

5.1 Inference for location-scale distributions

• Location-scale distributions have survivor function of the form $$\begin{array}{}\text{(5.1)}& S(y;u,b)=S_0(\frac{y-u}{b})-\mathrm{\infty }<y<\mathrm{\infty }\end{array}$$ $$-\mathrm{\infty }<u<\mathrm{\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{array}{rl}{S}_T(t;\alpha ,\beta )& =S_0(\frac{\log t-u}{b})\\ & =S_0^{\star }[(t/\alpha {)}^{\beta }]\end{array}$$

  • $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.

Likelihood based methods

• 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

Likelihood function

• For a censored sample $$\{(t_i,{\delta }_i),i=1,\dots ,n\},$$ the likelihood function $$\begin{array}{r}L(u,b)=\prod _{i=1}^n[\frac{1}{b}{f}_0(\frac{y_i-u}{b}){]}^{{\delta }_i}[S_0(\frac{y_i-u}{b}){]}^{1-{\delta }_i}\end{array}$$

  • $y_i=\log t_i$

  • $f_0(z)=-d{S}_0(z)/dz\to$ pdf

---

• The standardized variable $$z_i=\frac{y_i-u}{b}$$

• The likelihood function $$\begin{array}{r}L(u,b)=\prod _{i=1}^n[\frac{1}{b}{f}_0(z_i){]}^{{\delta }_i}[S_0(z_i){]}^{1-{\delta }_i}\end{array}$$

---

• The corresponding log-likelihood function $$\begin{array}{rl}\ell (u,b)& =-r\log b+\sum _{i=1}^n[{\delta }_i\log f_0(z_i)+(1-{\delta }_i)\log S_0(z_i)]\\ & =-r\log b+\sum _{i=1}^n{\ell }_i(z_i,{\delta }_i)\end{array}$$

  • $r=\sum _{i=1}^n{\delta }_i$

Score functions

$$\begin{array}{rl}\ell (u,b)& =-r\log b+\sum _{i=1}^n{\ell }_i(z_i,{\delta }_i)\\ & =-r\log b+\sum _{i=1}^n[{\delta }_i\log f_0(z_i)+(1-{\delta }_i)\log S_0(z_i)]\end{array}$$

$$\begin{array}{rl}\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}\\ & =\sum _{i=1}^n\frac{\partial {\ell }_i(z_i,{\delta }_i)}{\partial z_i}\times (\frac{-1}{b})\\ & =-\frac{1}{b}\sum _{i=1}^n[{\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}]\end{array}$$

---

$$\ell (u,b)=-r\log b+\sum _{i=1}^n[{\delta }_i\log f_0(z_i)+(1-{\delta }_i)\log S_0(z_i)]$$

$$\begin{array}{rl}\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}\\ & =-\frac{r}{b}+\sum _{i=1}^n\frac{\partial {\ell }_i(z_i,{\delta }_i)}{\partial z_i}\times (\frac{-z_i}{b})\\ & =-\frac{r}{b}-\frac{1}{b}\sum _{i=1}^n{z}_i[{\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}]\end{array}$$

Hessian matrix

$$\frac{\partial \ell (u,b)}{\partial u}=-\frac{1}{b}\sum _{i=1}^n[{\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}]$$

$$\begin{array}{rl}\frac{{\partial }^2\ell (u,b)}{\partial u^2}& =\frac{\partial }{\partial u}[\frac{\partial \ell (u,b)}{\partial u}]\\ & =\frac{1}{b^2}\sum _{i=1}^n[{\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}]\end{array}$$

---

$$\frac{\partial \ell (u,b)}{\partial b}=-\frac{r}{b}-\frac{1}{b}\sum _{i=1}^n{z}_i[{\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}]$$ $$\begin{array}{rl}\frac{{\partial }^2\ell (u,b)}{\partial b^2}& =\frac{r}{b^2}+\frac{2}{b^2}\sum _{i=1}^n{z}_i[{\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}]\\ & +\frac{1}{b^2}\sum _{i=1}^n{z}_i^2[{\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}]\end{array}$$

---

$$\begin{array}{rl}\frac{\partial \ell (u,b)}{\partial u}& =-\frac{1}{b}\sum _{i=1}^n[{\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}]\\ \frac{{\partial }^2\ell (u,b)}{\partial u\partial b}& =\frac{1}{b^2}\sum _{i=1}^n[{\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}]\\ & +\frac{1}{b^2}\sum _{i=1}^n{z}_i[{\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}]\end{array}$$

Score function and information matrix

$$\begin{array}{rl}U(u,b)& =\left[\begin{array}{c}\frac{\partial \ell (u,b)}{\partial u}\\ \frac{\partial \ell (u,b)}{\partial b}\end{array}\right]\\ I(u,b)& =-H(u,b)=-\left[\begin{array}{cc}\frac{{\partial }^2\ell (u,b)}{\partial u^2}& \frac{{\partial }^2\ell (u,b)}{\partial u\partial b}\\ \frac{{\partial }^2\ell (u,b)}{\partial b\partial u}& \frac{{\partial }^2\ell (u,b)}{\partial b^2}\end{array}\right]\end{array}$$

Statistical inference

• MLE

$$(\hat{u},\hat{b}{)}^{\prime }={\arg max}_{(u,b{)}^{\prime }\in \mathrm{\Theta }}\ell (u,b)$$

• Variance-covariance matrix

$$\text{var}(\hat{u},\hat{b})=[I(\hat{u},\hat{b}){]}^{-1}=\hat{V}$$

• Sampling distribution

$$\left(\begin{array}{c}\hat{u}\\ \hat{b}\end{array}\right)\sim {\mathcal{N}}_2\left(\begin{array}{c}\left[\begin{array}{c}u\\ b\end{array}\right],[I(\hat{u},\hat{b}){]}^{-1}\end{array}\right)$$

Wald type CIs

• For a large $n$, $(\hat{u},\hat{b}{)}^{\prime }$ follows a bivariate normal distribution with mean $(u,b{)}^{\prime }$ 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{array}{r}{Z}_1=\frac{\hat{u}-u}{se(\hat{u})},Z_2=\frac{\hat{b}-b}{se(\hat{b})},Z_2^{\prime }=\frac{\log \hat{b}-\log b}{se(\log \hat{b})}\end{array}$$

  • $se(\log \hat{b})=se(\hat{b})/\hat{b}$

• $(1-p)100\mathrm{\%}$ confidence intervals $$\begin{array}{rl}\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{array}$$

Quantiles

• $p$th quantile of log lifetime $Y$ $$\begin{array}{rl}P(Y\le y_p)=p\Rightarrow S_0(\frac{y_p-u}{b})& =1-p\\ \frac{y_p-u}{b}& =S_0^{-1}(1-p)\\ y_p& =u+b{w}_p\end{array}$$

  • $w_p=S_0^{-1}(1-p)=F_0^{-1}(p)\to$ $p$th quantile of $S_0(z)$, the standardize distribution

---

• Estimate of $p$th quantile and the corresponding standard error $$\begin{array}{rl}{\hat{y}}_p& =\hat{u}+\hat{b}{w}_p\\ se({\hat{y}}_p)& =\sqrt{{\hat{V}}_{11}+w_p^2{\hat{V}}_{22}+2w_p{\hat{V}}_{12}}\end{array}$$

• Pivotal quantity $$Z_p=\frac{{\hat{y}}_p-y_p}{se({\hat{y}}_p)}\sim \mathcal{N}(0,1)$$

• $(1-q)100\mathrm{\%}$ confidence interval for $y_p$ $${\hat{y}}_p\pm z_{1-q/2}se({\hat{y}}_p)$$

Likelihood ratio procedures

• 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 $${\mathrm{\Lambda }}_1(u_0)=2\ell (\hat{u},\hat{b})-2\ell (u_0,\tilde{b}(u_0))$$

• MLEs $$\begin{array}{rl}(\hat{u},\hat{b}{)}^{\prime }& ={\arg max}_{(u,b{)}^{\prime }\in \mathrm{\Theta }}\ell (u,b)\text{unrestricted}\\ \tilde{b}(u_0)& ={\arg max}_{b\in {\mathrm{\Theta }}_1}\ell (u_0,b)\text{under }{H}_0\end{array}$$

---

• Under $H_0:u=u_0$, asymptotically ${\mathrm{\Lambda }}_1(u_0)\sim {\chi }_{(1)}^2$

  • Approximate two-sided $(1-p)100\mathrm{\%}$ confidence interval for $u$ can be obtained as the set of values of $u_0$ for which $${\mathrm{\Lambda }}_1(u_0)\le {\chi }_{(1),1-p}^2$$

---

Homework - 1

• Obtain the expression of likelihood ratio test statistic based confidence interval for the scale parameter $b$

---

• The $p$th quantile of location-scale distribution can be expressed as $$y_p=u+w_{p}b,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 $$\begin{array}{}\text{(5.2)}& \mathrm{\Lambda }(y_{p_0})=2\ell (\hat{u},\hat{b})-2\ell (\tilde{u},\tilde{b})\end{array}$$

---

Steps

• The estimates $\hat{u}$ and $\hat{b}$ are MLEs under $H_1$ $$(\hat{u},\hat{b}{)}^{\prime }={\arg max}_{(u,b{)}^{\prime }\in \mathrm{\Theta }}\ell (u,b)$$

• Steps to obtain MLEs $\tilde{u}$ and $\tilde{b}$, under $H_0:y_p=y_{p_0}$

  1. Under $H_0$, $y_{p_0}=u+w_{p}b\Rightarrow u=y_{p_0}-b{w}_p$

  2. $\tilde{b}={\arg max}_{b\in \mathrm{\Theta }}\ell (y_{p_0}-w_{p}b,b)$

  3. $\tilde{u}=y_{p_0}-w_p\tilde{b}$

• $(1-q)100\mathrm{\%}$ Confidence interval for $y_p$ can be obtained from the set of $y_{p_0}$ values such that $$\mathrm{\Lambda }(y_{p_0})\le {\chi }_{(1),1-q}^2$$

---

• 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 $$\mathrm{\Lambda }(s_0)=2\ell (\hat{u},\hat{b})-2\ell (\tilde{u},\tilde{b})$$

---

Steps

• Steps for obtaining MLEs $\tilde{u}$ and $\tilde{b}$ under $H_0$

  1. Under $H_0:S(y_0)=s_0$ $$S(y_0)=S_0(\frac{y_0-u}{b})=s_0\Rightarrow u=y_0-S_0^{-1}(s_0)b$$

  2. $\tilde{b}={\arg max}_{b\in \mathrm{\Theta }}\ell (y_0-S_0^{-1}(s_0)b,b)$

  3. $\tilde{u}=y_0-S_0^{-1}(s_0)\tilde{b}$

• The $(1-p)100\mathrm{\%}$ confidence interval for $S(y_0)$ can be defined as the set of $s_0$ values such that $$\begin{array}{r}\mathrm{\Lambda }(s_0)\le {\chi }_{(1),1-p}^2\end{array}$$

---

• 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

5.2 Weibull and extreme-value distributions

• The pdf of Weibull distribution $$\begin{array}{r}f(t;\alpha ,\beta )=\frac{\beta }{\alpha }(\frac{t}{\alpha }{)}^{\beta -1}\exp [-(t/\alpha {)}^{\beta }]\end{array}$$

  • $\alpha >0$ and $\beta >0$ are scale and shape parameters, respectively

---

• The pdf of extreme-value distribution $$\begin{array}{rl}f(y;u,b)& =\frac{1}{b}\exp [(y-u)/b]\exp [-e^{(y-u)/b}]\\ & =\frac{1}{b}{f}_0(\frac{y-u}{b})\end{array}$$

  • $u=\log \alpha$

  • $b=(1/\beta )$

• Extreme-value distribution is used to make inferences about Weibull distribution

Likelihood based inference procedures

• Censored sample $$\{(t_i,{\delta }_i),i=1,\dots ,n\}$$

• Define $$y_i=\log t_i\text{and}{z}_i=(y_i-u)/b$$

---

• General expression of likelihood function $$\ell (u,b)=-r\log b+\sum _{i=1}^n[{\delta }_i\log f_0(z_i)+(1-{\delta }_i)\log S_0(z_i)]$$

• For extreme-value distribution $$\begin{array}{rl}{S}_0(z)& =\exp (-e^z)\\ f_0(z)& =-\frac{d}{dz}{S}_0(z)=\exp (z-e^z)\end{array}$$

---

$$\ell (u,b)=-r\log b+\sum _{i=1}^n[{\delta }_i\log f_0(z_i)+(1-{\delta }_i)\log S_0(z_i)]$$

• Log-likelihood function for EV distribution $$\begin{array}{}\text{(5.3)}& \ell (u,b)=-r\log b+\sum _{i=1}^n({\delta }_i{z}_i-e^{z_i})\end{array}$$
  • $r=\sum _i{\delta }_i$
• This log-likelihood function $\ell (u,b)$ is easily maximized to give $\hat{u},\hat{b}$ (using software)

Score functions

• 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{array}{rl}\frac{\partial \ell (u,b)}{\partial u}& =-\frac{1}{b}\sum _{i=1}^n[{\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}]\\ \frac{\partial \ell (u,b)}{\partial b}& =-\frac{r}{b}-\frac{1}{b}\sum _{i=1}^n{z}_i[{\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}]\end{array}$$

---

• For extreme-value distribution $$\begin{array}{rl}\frac{\partial \log f_0(z)}{\partial z}& =\frac{\partial }{\partial z}\log \{\exp (z-e^z)\}=1-e^z\\ \frac{\partial \log S_0(z)}{\partial z}& =\frac{\partial }{\partial z}\log \{\exp (-e^z)\}=-e^z\end{array}$$

• 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}$

Hessian matrix

Hessian matrix

• Hessian matrix at MLEs $\hat{u}$ and $\hat{b}$ $$H(\hat{u},\hat{b})=-\frac{1}{{\hat{b}}^2}\left[\begin{array}{cc}r& \sum _{i=1}^n{\hat{z}}_i{e}^{{\hat{z}}_i}\\ \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{array}\right]$$

Proof

• For extreme-value distribution

$$\begin{array}{rl}\frac{{\partial }^2\ell (u,b)}{\partial u^2}& =\frac{1}{b^2}\sum _{i=1}^n[{\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}]=-\frac{1}{b^2}\sum _{i=1}^n{e}^{z_i}\end{array}$$ where $$\begin{array}{rl}\frac{{\partial }^2\log f_0(z)}{\partial z^2}& =\frac{\partial }{\partial z}\{1-e^z\}=-e^z=\frac{{\partial }^2\log S_0(z)}{\partial z^2}\end{array}$$

$$\begin{array}{r}\frac{{\partial }^2\ell (u,b)}{\partial b^2}{|}_{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{array}$$

From Section 5.1.4, it can be shown that the Hessian matrix will be

$$\begin{array}{rl}\frac{{\partial }^2\ell (u,b)}{\partial b^2}& =\frac{r}{b^2}+\frac{2}{b^2}\sum _{i=1}^n{z}_i[{\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}]\\ & +\frac{1}{b^2}\sum _{i=1}^n{z}_i^2[{\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}]\\ & =\frac{r}{b^2}+\frac{2}{b^2}\sum _{i=1}^n{z}_i[{\delta }_i(1-e^{z_i})-(1-{\delta }_i)e^{z_i}]-\frac{1}{b^2}\sum _{i=1}^n{z}_i^2{e}^{z_i}\\ & =\frac{r}{b^2}+\frac{2}{b^2}\sum _{i=1}^n{z}_i[{\delta }_i-e^{z_i}]-\frac{1}{b^2}\sum _{i=1}^n{z}_i^2{e}^{z_i}\end{array}$$

$$\begin{array}{r}\frac{{\partial }^2\ell (u,b)}{\partial b^2}{|}_{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{array}$$

$$\begin{array}{rl}\frac{{\partial }^2\ell (u,b)}{\partial u\partial b}& =\frac{1}{b^2}\sum _{i=1}^n[{\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}]\\ & +\frac{1}{b^2}\sum _{i=1}^n{z}_i[{\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}]\\ & =\frac{1}{b^2}\sum _{i=1}^n[{\delta }_i(1-e^{z_i})-(1-{\delta }_i)e^{z_i}-z_i{e}^{z_i}]\\ & =\frac{1}{b^2}\sum _{i=1}^n[{\delta }_i-e^{z_i}-z_i{e}^{z_i}]=\frac{1}{b^2}[r-\sum _{i=1}^n{e}^{z_i}-\sum _{i=1}^n{z}_i{e}^{z_i}]\end{array}$$

$$\begin{array}{r}\frac{{\partial }^2\ell (u,b)}{\partial u\partial b}{|}_{u=\hat{u},b=\hat{b}}=-\frac{1}{{\hat{b}}^2}\sum _{i=1}^n{\hat{z}}_i{e}^{{\hat{z}}_i}\end{array}$$

Covariance matrix

• Observed information matrix at MLEs $\hat{u}$ and $\hat{b}$ $$\begin{array}{rl}I(\hat{u},\hat{b})& =-H(\hat{u},\hat{b})\\ & =\frac{1}{{\hat{b}}^2}\left[\begin{array}{cc}r& \sum _{i=1}^n{\hat{z}}_i{e}^{{\hat{z}}_i}\\ \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{array}\right]\end{array}$$

• Covariance matrix of $(\hat{u},\hat{b}{)}^{\prime }$ $$\begin{array}{rl}\hat{V}& =[I(\hat{u},\hat{b}){]}^{-1}\end{array}$$

---

• 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 }{)}^{\prime }$ (using multivariate delta method) $$\mathrm{var}(\hat{\alpha },\hat{\beta })=\mathit{G}\hat{V}{\mathit{G}}^{\prime }$$ where $$G=\left[\begin{array}{cc}\frac{\partial g_1(u,b)}{du}& \frac{\partial g_1(u,b)}{\partial b}\\ \frac{\partial g_2(u,b)}{\partial u}& \frac{\partial g_2(u,b)}{\partial b}\end{array}\right]=\left[\begin{array}{cc}{e}^{\hat{u}}& 0\\ 0& -\frac{1}{{\hat{b}}^2}\end{array}\right]$$

  • $\alpha =g_1(u,b)=e^u$

  • $\beta =g_2(u,b)=(1/b)$

---

• Wald-type statistics based $100(1-p)\mathrm{\%}$ CI for $u$ and $b$ $$\begin{array}{rl}\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{array}$$

CI for $(u,b)$ (LRT based)

• Log-likelihood function corresponding to $H_0:b=b_0$ is (from Equation 5.3) $$\ell (u,b_0)=-r\log b_0+\sum _{i=1}^n[{\delta }_i(\frac{y_i-u}{b_0})-e^{(y_i-u)/b_0}]$$
• MLE of $u$ under $H_0:b=b_0$ $$\begin{array}{rl}\frac{\partial \ell (u,b_0)}{\partial u}{|}_{u=\tilde{u}}=0\Rightarrow & -\frac{1}{b_0}[r-\sum _{i=1}^n{e}^{(y_i-\tilde{u})/b_0}]=0\\ \Rightarrow & \tilde{u}(b_0)=b_0\log [\frac{1}{r}\sum _{i=1}^n{e}^{y_i/b_0}]\end{array}$$

---

• LRT statistics $$\mathrm{\Lambda }(b_0)=2\ell (\hat{u},\hat{b})-2\ell (\tilde{u}(b_0),b_0)$$

• $100(1-p)\mathrm{\%}$ CI for $b$ is defined by the set of $b_0$ values such that $${\mathrm{\Lambda }}_1(b_0)\le {\chi }_{(1),1-p}^2$$

• Similarly, confidence interval for $u$ can be obtained using the corresponding LRT statistics (Homework)

CI for quantiles

• The $p$th quantile of $Y\sim EV(u,b)$ $$\begin{array}{rl}S(y_p)& =S_0(\frac{y_0-u}{b})=(1-p)\\ & \exp [-\exp (\frac{y_p-u}{b})]=(1-p)\\ & \frac{y_p-u}{b}=\log [-\log (1-p)]=S_0^{-1}(1-p)=w_p\\ & y_p=u+w_{p}b\end{array}$$

---

CI for quantiles (Wald)

• 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) $$\mathrm{var}({\hat{y}}_p)=\left[\begin{array}{cc}1& w_p\end{array}\right]\hat{V}\left[\begin{array}{c}1\\ w_p\end{array}\right]={\hat{V}}_{11}+{\hat{V}}_{22}{w}_p^2+2{\hat{V}}_{12}{w}_p$$

• Large sample based $100(1-q)\mathrm{\%}$ confidence interval for $y_p$ $${\hat{y}}_p\pm z_{1-q/2}se({\hat{y}}_p)$$

• Find the $100(1-q)\mathrm{\%}$ confidence interval for $t_p$

---

CI for quantiles (LRT)

• 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 $$\mathrm{\Lambda }(y_{p_0})=2\ell (\hat{u},\hat{b})-2\ell (\tilde{u},\tilde{b})$$

  • The procedure of obtaining parameter estimates $\tilde{u}$ and $\tilde{b}$ (under $H_0$) is explained in Section 5.1.9.1)

• LRT statistic based $(1-q)100\mathrm{\%}$ confidence interval for $y_p$ can be obtained from the set of $y_{p_0}$ values such that $$\mathrm{\Lambda }(y_{p_0})\le {\chi }_{(1),1-q}^2$$

CI for $S(\cdot )$ (Wald)

• To obtain confidence interval for survival probability $$S(y_0)=S_0(\frac{y_0-u}{b})=\exp [-\exp (\frac{y_0-u}{b})]$$

• We can defined $$\psi =S_0^{-1}(S(y_0))=\log [-\log (S(y_0))]=\frac{y_0-u}{b}$$

• MLE and SE $$\begin{array}{rl}\hat{\psi }& =\frac{y_0-\hat{u}}{\hat{b}}\\ \mathrm{var}(\hat{\psi })& ={\mathit{a}}^{\prime }\hat{V}\mathit{a}=\left[\begin{array}{cc}-1/\hat{b}& -\hat{\psi }/b\end{array}\right]\left[\begin{array}{cc}{\hat{V}}_{11}& {\hat{V}}_{12}\\ {\hat{V}}_{21}& {\hat{V}}_{22}\end{array}\right]\left[\begin{array}{c}-1/\hat{b}\\ -\hat{\psi }/b\end{array}\right]\end{array}$$

---

• $(1-p)100\mathrm{\%}$ CI for $\psi$ $$\begin{array}{rl}\hat{\psi }-se(\hat{\psi })z_{1-p2}& <\psi \le \hat{\psi }+se(\hat{\psi })z_{1-p/2}\\ L& <\psi \le U\end{array}$$

• Confidence interval for $S(y_0)$ $$\begin{array}{rl}L& <\log [-\log (S(y_0))]\le U\\ \exp [-\exp (U)]& <S(y_0)\le \exp [-\exp (L)]\end{array}$$

CI for $S(\cdot )$ (LRT)

• Consider the null hypothesis $H_0:S(y_0)=s_0$, where $$S(y_0)=\exp [-\exp (\frac{y_0-u}{b})]$$

• The $(1-p)100\mathrm{\%}$ confidence interval for $S(y_0)$ can be defined as the set of $s_0$ values such that $\mathrm{\Lambda }(s_0)\le {\chi }_{(1),1-p}^2$, where $$\mathrm{\Lambda }(s_0)=2\ell (\hat{u},\hat{b})-2\ell (\tilde{u},\tilde{b})$$

  • The procedure of obtaining parameter estimates $\tilde{u}$ and $\tilde{b}$ (under $H_0$) is explained in Section 5.1.9.2)

Example 5.2.1:

• 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)

---

Remission time data

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

• Each group has 21 subjects, and all subjects of the placebo group were failed (end of remission) and drug group (6MP) has 12 censored times

---

• 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 $$\{(t_i,{\delta }_i),i=1,\dots ,n\}$$

• Log-likelihood function $$\ell (\alpha ,\beta )=\sum _{i=1}^n[{\delta }_i\log f(t_i;\alpha ,\beta )+(1-{\delta }_i)\log S(t_i;\alpha ,\beta )]$$

• MLEs $$(\hat{\alpha },\hat{\beta }{)}^{\prime }={\arg max}_{(\alpha ,\beta {)}^{\prime }\in \mathrm{\Theta }}\ell (\alpha ,\beta )$$

---

Analysis of remission time data (Extreme-value distribution)

• Define $y=\log t$ and corresponding probability density and survivor function $$\begin{array}{rl}f(y;u,b)& =\frac{1}{b}\exp [(y-u)/b-e^{(y-u)/b}]\\ S(y;u,b)& =\exp [-e^{(y-u)/b}]\end{array}$$

• Log-likelihood function $$\begin{array}{r}{\ell }_{ev}(u,b)=\log \prod _{i=1}^n[f(y_i;u,b){]}^{{\delta }_i}[S(y_i;u,b){]}^{1-{\delta }_i}\end{array}$$

• MLEs $$(\hat{u},\hat{b}{)}^{\prime }={\arg max}_{(u,b{)}^{\prime }\in \mathrm{\Theta }}{\ell }_{ev}(u,b)$$

survreg function

• R function survreg() 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 $$\mathtt{\text{formula = Surv(time, status)}}\sim \mathtt{\text{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{array}{rl}\text{weibull}\to & \mathtt{\text{formula = Surv(time, status)}}\sim \mathtt{\text{1}}\\ \text{extreme}\to & \mathtt{\text{formula = Surv(log(time), status)}}\sim \mathtt{\text{1}}\end{array}$$

---

Data for the treatment (6MP) group

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",…

---

Analysis of the treatment group using survreg function

w_sreg_6mp <- survreg(Surv(time, status) ~ 1, 
                  data = d6mp, dist = "weibull")

ev_sreg_6mp <- survreg(Surv(log(time), status) ~ 1, 
                  data = d6mp, dist = "extreme")

---

Extreme-value distribution

• Estimates of model parameters $u$ and $\log b$

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- 1

---

• Variance-covariance matrix of $(\hat{u},\log \hat{b})$

vcov(ev_sreg_6mp)

            (Intercept) Log(scale)
(Intercept)  0.07473057 0.03305811
Log(scale)   0.03305811 0.07750538

---

Analysis of remission time data using survreg function

• The survreg() function returns estimates of $(u,\log b{)}^{\prime }$ and corresponding variance matrix

• For making inference about Weibull distribution, followings are required

  1. estimate of $(u,b{)}^{\prime }$ and the corresponding variance matrix

  2. estimate of $(\alpha ,\beta {)}^{\prime }$ and the corresponding variance matrix

• It is important to understand the methods to obtain estimates and the corresponding variance of $(u,b{)}^{\prime }$ and $(\alpha ,\beta {)}^{\prime }$ from the estimates and the corresponding variance of $(u,\log b{)}^{\prime }$

---

Homework

• Obtain the variance-covarince matrix of $(\hat{\alpha },\hat{\beta }{)}^{\prime }$ and $(\hat{u},\hat{b}{)}^{\prime }$

CIs of $(\alpha ,\beta )$ (6MP group)

• 95% CI using the sampling distribution of $(\hat{\alpha },\hat{\beta })$ $$\begin{array}{rl}\hat{\alpha }\pm z_{.975}se(\hat{\alpha })& =33.765\pm (1.96)(9.23)\\ & =15.674\text{to}51.856\\ \hat{\beta }\pm z_{.975}se(\hat{\beta })& =1.354\pm (1.96)(0.377)\\ & =0.615\text{to}2.092\end{array}$$

---

• 95% CI using the sampling distribution of $(\hat{u},\log \hat{b})$ $$\begin{array}{rl}\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{array}$$

  • Similarly $$\begin{array}{rl}\log \hat{b}\pm z_{.975}se(\log \hat{b})& =-0.849\text{to}0.243\\ \hat{\beta }\pm z_{.975}se(\hat{\beta })& =1/\exp (0.243)\text{to}1/\exp (-0.849)\\ & =0.784\text{to}2.336\end{array}$$

---

(Obtain the variance matrix of $(\hat{u},\hat{b})$ using the sampling distribution of $(\hat{u},\log \hat{b}{)}^{\prime }$)

• 95% CI using the sampling distribution of $(\hat{u},\hat{b})$

$$\begin{array}{rl}\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{array}$$

• Similarly $$\begin{array}{rl}\hat{b}\pm z_{.975}se(\hat{b})& =0.336\text{to}1.142\\ \hat{\beta }\pm z_{.975}se(\hat{\beta })& =1/1.142\text{to}1/0.336\\ & =0.876\text{to}2.979\end{array}$$

---

• Using the method described in Section 5.2.5, we obtain the LRT-based CIs for $u$ and $b$

Plot of LRT statistic against different null values $u_0$ and 95% confidence interval for $u$ and $\alpha$

---

Plot of LRT statistic against different null values $b_0$ and 95% confidence interval for $\log b$ and $\beta$

---

95% confidence intervals for $\alpha$ and $\beta$ by different methods

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)

Analyses for Placebo group

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

---

95% confidence intervals for $\alpha$ and $\beta$ by different methods

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)

Quantiles and their CIs

• 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}$

---

95% confidence intervals for different quantiles of treatment group (6-MP)

$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

---

Plot of LRT statistic against different null values $y_{p_0}$ and 95% confidence interval for $y_{.25}$ and $t_{.25}$ (6-MP group)

---

Plot of LRT statistic against different null values $y_{p_0}$ and 95% confidence interval for $y_{.5}$ and $t_{.5}$ (6-MP group)

---

95% confidence intervals of different quantiles using Wald and LRT method (6-MP group)

$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

---

95% confidence intervals for different quantiles using Wald and LRT method (placebo group)

$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

Survivor function

• For Weibull distribution, the expression of survivor function $$S(t;\alpha ,\beta )=\exp (-(t/\alpha {)}^{\beta })$$

• Estimated survivor function $$S(t;\hat{\alpha },\hat{\beta })=\exp (-(t/\hat{\alpha }{)}^{\hat{\beta }})$$

par	6-MP	placebo
$\alpha$	33.765	9.482
$\beta$	1.354	1.370

---

Comparison survival probabilities of between two treatment groups

---

Figure 5.1: Comparison of parametric (Weibull) and non-parametric (step-function) estimates of survivor function using remission time data

---

Homework

• Obtain Wald and LRT statistics based confidence interval for the survival probability $S(10)$