Chapter 6

(AST405) Lifetime data analysis

6 Parametric Regression Models

6.1 Log-location-scale (Accelerated Failure Time) Regression Models

Linear regression model

• Distributional assumption for the response $$(Y|x)=Y(x)\sim \mathcal{N}(\mu (x),{\sigma }^2)$$
• Regression model for the parameters $$\begin{array}{rl}\mu (x)& ={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p={\mathbf{x}}^{\prime }\mathit{\beta }\\ \mathrm{var}(Y|x)& ={\sigma }^2\end{array}$$

---

• Instead of the parameters, linear regression model can be defined in terms of other functions, such as survivor function $$\begin{array}{}\text{(6.1)}& \begin{array}{rl}{S}_Y(y)& =Pr(Y>y)\\ & =1-\mathrm{\Phi }(\frac{y-\mu (x)}{\sigma })\end{array}\end{array}$$

Regression models for lifetimes

• Similar to continuous and binary responses, regression analysis of lifetimes involves specifications for the distribution of a lifetime $(T)$ given a vector of $p$-dimensional (say) covariate $\mathbf{x}$ $$(T|\mathbf{x})=T(\mathbf{x})$$

---

• For parametric regression models for lifetimes $T$, parameters (e.g. scale and shape parameters) need to be defined as a function of measured covariates (linear predictors)

• It requires selecting a link function (e.g. identity, log, logit, etc.) for relating model parameters with linear predictors

• Similar to linear and logistic regression models, maximum likelihood method of estimation is used to estimate parameters of the model

Log-location-scale AFT model

• For a lifetime that follows a distribution of the log-location-scale family of distributions, the survivor function of lifetime $T$ for a given covariate vector $\mathbf{x}$ is defined as $$\begin{array}{}\text{(6.2)}& S(t|\mathbf{x})=S_0^{\star }([t/\alpha (\mathbf{x}){]}^{\delta })\end{array}$$

  • Scale parameter $\alpha (\mathbf{x})$ is defined as a function of covariate vector $\mathbf{x}$

  • Shape parameter $\delta$ does not depend on $\mathbf{x}$

  • Survivor function of the corresponding standardized distribution $S_0^{\star }(x)=S_0(\log x)$ is defined earlier

---

• For a log-lifetime that follows a distribution of the location-scale family of distribution, the survivor function of log-lifetime $Y$ for a given covariate vector $\mathbf{x}$ is defined as $$\begin{array}{}\text{(6.3)}& S(y|\mathbf{x})=S_0(\frac{y-u(\mathbf{x})}{b})\end{array}$$

  • Location parameter $u(\mathbf{x})$ is defined as a function covariate vector $\mathbf{x}$

  • Scale parameter $b$ does not depend on $\mathbf{x}$

• The model (Equation 6.3) for log-lifetime is similar to the linear regression model (Equation 6.1) with $$\mu (\mathbf{x})=u(\mathbf{x}),\sigma =b,\text{and}\mathrm{\Phi }(x)=1-S_0(x)$$

---

• The model for lifetime (Equation 6.2) or log-lifetime (Equation 6.3) is known as accelerated failure time (AFT) model $$\begin{array}{rl}S(t|\mathbf{x})& =S_0^{\star }([t/\alpha (\mathbf{x}){]}^{\delta })\\ S(y|\mathbf{x})& =S_0(\frac{y-u(\mathbf{x})}{b})\end{array}$$

• Models for the parameters $\alpha (\mathbf{x})$ and $u(\mathbf{x})$ are defined so that associated parametric restrictions are satisfied, $\alpha (\mathbf{x})>0$ and $-\mathrm{\infty }<u(\mathbf{x})<\mathrm{\infty }$, e.g. $$\begin{array}{rl}u(\mathbf{x})& ={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p={\mathbf{x}}^{\prime }\mathit{\beta }\\ \alpha (\mathbf{x})& =\exp ({\mathbf{x}}^{\prime }\mathit{\beta })\end{array}$$

---

• AFT model can also be expressed as $$\begin{array}{}\text{(6.4)}& \frac{Y-u(\mathbf{x})}{b}=Z\Rightarrow Y=u(\mathbf{x})+bZ\end{array}$$

  • $Z\sim S_0(z)$, i.e. $Z$ follows a standardized log-location-scale distribution, e.g. standard normal or extreme-value distributions with location 0 and scale 1, etc.

• Linear regression model (Equation 6.1) can also be expressed as Equation 6.4: $$Y=\mu (\mathbf{x})+\sigma Z,Z\sim \mathcal{N}(0,1)$$

---

• In AFT model defined in terms of the distribution of lifetime $T$, covariates alter the time scale

  • If $\alpha (\mathbf{x})=\exp ({\mathbf{x}}^{\prime }\mathit{\beta })>1$, the effect of covariate vector is to increase time (decelerate time)

  • If $\alpha (\mathbf{x})=\exp ({\mathbf{x}}^{\prime }\mathit{\beta })<1$, the effect of covariate vector is to shorten time (accelerate time)

---

• The accelerated failure time model is a general model for survival data, in which explanatory variables measured on an individual are assumed to act multiplicatively on the time-scale

• Log-location-scale AFT models are a special case of AFT models where the log of survival time follows a location-scale distribution.

• AFT models assume that covariates accelerate or decelerate the time to event.

---

• The following example is described in Collett (2015)

• Suppose patients are randomized to receive one of the two treatments $A$ (standard) and $B$ (new)

• Under an accelerated failure time model, the survival time of an individual on the new treatment is taken to be a multiple of the survival time for an individual on the standard treatment.

• Thus, the effect of the new treatment is to “speed up” or “slow down” the passage of time

---

• For a specific time $t$ $$S(t|\text{trt}=B)=S(t\alpha |\text{trt}=A)$$

• One interpretation of this model is that the lifetime of an individual on the new treatment ($B$) is $\alpha$ times the lifetime that the individual would have experienced under the standard treatment $(A)$

• When the end-point of concern is the death of a patient

  • $\alpha >1$ new treatment is promoting longevity
  • $\alpha <1$ new treatment is worse (accelerating death)

• The quantity $\alpha$ is therefore termed the acceleration factor

---

• The acceleration factor can also be interpreted in terms of the median survival times of patients on the new and standard treatments, $t_A(50)$ and $t_B(50)$ $$S_B\{t_B(50)\}=S_A\{t_A(50)\}=0.50$$
• Under AFT model $$S_B\{t_B(50)\}=S_A\{t_B(50)/\alpha \}\Rightarrow t_B(50)=\alpha t_A(50)$$
• Under the AFT model, the median survival time of a patient on the new treatment is $\alpha$ times that of a patient on the standard treatment

---

• Under AFT model, the survivor functions with covariate vectors ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ can be compared as $$S(t|{\mathbf{x}}_1)=S(ct|{\mathbf{x}}_2)$$

  • If $c>1$, subjects with covariate ${\mathbf{x}}_2$ survives longer compared to subjects with covariate vector ${\mathbf{x}}_1$

  • If $c<1$ subjects with covariate ${\mathbf{x}}_2$ survives shorter compared to subjects with covariate vector ${\mathbf{x}}_1$

---

• Under AFT model, $S_1(t)=S_2(ct)$ for $c>0$, we can express the mean survival time ${\mu }_2$ of Population 2 can be expressed in terms of ${\mu }_1$, mean survival time of Population 1 as $$\begin{array}{rl}{\mu }_2& ={\int }_0^{\mathrm{\infty }}{S}_2(t)dt\\ & =c{\int }_0^{\mathrm{\infty }}{S}_2(cu)du\\ & =c{\int }_0^{\mathrm{\infty }}{S}_1(u)du\\ & =c{\mu }_1\end{array}$$

---

• In general, let $\phi$ is a population quantity such that $S(\phi )=\theta$ for some $\theta \in (0,1)$ and $$S_2({\phi }_2)=\theta =S_1({\phi }_1)=S_2(c{\phi }_1)$$
• Then ${\phi }_2=c{\phi }_1$, i.e., under the AFT model, the expected survival time, median survival time of population 2 all are $c$ times as much as those of population 1

---

Comparison between two log-location density functions with covariate vectors $x_1$ and $x_2$, where $u(x_2)>u(x_1)$

---

Comparison between two log-location survival functions with covariate vectors $x_1$ and $x_2$, where $u(x_2)>u(x_1)$

---

Comparison between two log-location survival functions with the same scale parameters, but different location parameters

---

Comparison between two log-location survival functions with the same location parameters, but different scale parameters

Proportional hazards model

• There are two approaches to regression modeling for lifetimes

  1. AFT model, where the effects of covariates are assessed by comparing corresponsing time scales

  2. Hazards model, where effects of covariates on the hazard function are studied

---

• The most common hazards model is the proportional hazards model (Cox 1972), where hazard function for lifetime $T$ given $\mathbf{x}$ is defined as $$\begin{array}{rl}h(t|\mathbf{x})& =h_0(t)r(\mathbf{x})\end{array}$$

  • $r(\mathbf{x})\to$ a positive-valued function of linear predictor, e.g. $r(\mathbf{x})=\exp ({\mathbf{x}}^{\prime }\mathit{\beta })$, which does not include the intercept term

  • $h_0(t)\to$ a positive-valued function, which is known as baseline hazards function, i.e. $h(t|\mathbf{x}=\mathbf{0})=h_0(t)$

  • $h_0(t)$ could be either fully parametric or unspecified

---

If you take two individuals with covariates $x_1$ and $x_2$:

$$\frac{h(t|x_1)}{h(t|x_2)}=\frac{h_0(t)e^{\beta x_1}}{h_0(t)e^{\beta x_2}}=e^{\beta (x_1-x_2)}$$

This ratio does not depend on time (t), this is exactly the proportional hazards property.

---

$$h(t|x)=h_0(t)e^{x\beta }$$

• For a binary predictor $x$ (1=male, 0=female), the hazard ratio can be defined as $$\begin{array}{rl}\frac{h(t|x=1)}{h(t|x=0)}& =\frac{h_0(t)e^{\beta }}{h_0(t)}=e^{\beta }\\ h(t|x=1)& =h(t|x=0)e^{\beta }\end{array}$$

• $\beta >0\Rightarrow$ Hazard of the event is higher for male compared to female

---

• Under proportional hazards model, the cumulative hazard function is defined as $$\begin{array}{rl}H(t|\mathbf{x})& ={\int }_0^{t}h(u|\mathbf{x})du\\ & =r(\mathbf{x}){\int }_0^t{h}_0(u)du\\ & =r(\mathbf{x})H_0(t)\end{array}$$

---

• Under proportional hazards model, the survivor function is defined as $$\begin{array}{rl}S(t|\mathbf{x})& =e^{-H(t|\mathbf{x})}=e^{-r(\mathbf{x})H_0(t)}=[S_0(t){]}^{r(\mathbf{x})}\end{array}$$

  • $S_0(t)\to$ baseline survivor function and $r(\mathbf{x})>0$

  • Interpret the survival probabilities for the following cases $$(a)r(\mathbf{x})>1\text{and}(b)r(\mathbf{x})<1$$

---

Under proporitonal hazards model, comparison between baseline survivor function $S_0(t)$ and $S_1(t|x)=[S_0(t){]}^{0.5}$

---

Under proporitonal hazards model, comparison between baseline survivor function $S_0(t)$ and $S_1(t|x)=[S_0(t){]}^{1.5}$

---

Under proporitonal hazards model, comparison between hazard functions $H(t|x_1)$ and $H(t|x_2)=1.5H(t|x_1)$

Parametric proportional hazards model

• Depending on whether the baseline hazard function $h_0(t)$ is fully parametric or not, a PH model $h(t|\mathbf{x})=h_0(t)r(\mathbf{x})$ could be either parametric or semi-parametric

  • PH model is parametric if $h_0(t)=h_1(\mathit{\alpha },t)$ for some parameter vector $\mathit{\alpha }$

  • PH model is semi-parametric if $h_0(t)$ is unspecified

---

• Weibull model can be defined as both AFT and PH model

Weibull regression model

• Weibull as an AFT model $$\begin{array}{rl}S(t|\mathbf{x})& =\exp (-[t/\alpha (\mathbf{x}){]}^{\delta })\end{array}$$ where $$\begin{array}{}\text{(6.5)}& \begin{array}{rl}\alpha (\mathbf{x})& =\exp ({\mathbf{x}}^{\prime }{\mathit{\beta }}_{\mathsf{A}\mathsf{F}\mathsf{T}})\end{array}\end{array}$$

---

• Weibull as a PH model $$\begin{array}{rl}h(t|\mathbf{x})& =\frac{\delta }{\alpha (\mathbf{x})}[\frac{t}{\alpha (\mathbf{x})}{]}^{\delta -1}\\ & =(\delta t^{\delta -1})[\alpha (\mathbf{x}){]}^{-\delta }\\ & =h_1(\delta ,t)r(\mathbf{x})\end{array}$$

  • Assume $$\begin{array}{}\text{(6.6)}& \begin{array}{rl}r(\mathbf{x})=\exp ({\mathbf{x}}^{\prime }{\mathit{\beta }}_{\mathsf{P}\mathsf{H}})& =[\alpha (\mathbf{x}){]}^{-\delta }\\ \Rightarrow \exp (-{\mathbf{x}}^{\prime }{\mathit{\beta }}_{\mathsf{P}\mathsf{H}}/\delta )& =\alpha (\mathbf{x})\end{array}\end{array}$$

---

• Equating the expression of $\alpha (\mathbf{x})$ from the AFT (Equation 6.5) and PH (Equation 6.6) Weibull model, we can show $$\begin{array}{rl}\exp (-{\mathbf{x}}^{\prime }{\mathit{\beta }}_{\mathsf{P}\mathsf{H}}/\delta )& =\alpha (\mathbf{x})=\exp ({\mathbf{x}}^{\prime }{\mathit{\beta }}_{\mathsf{A}\mathsf{F}\mathsf{T}})\\ \Rightarrow {\mathit{\beta }}_{\mathsf{P}\mathsf{H}}& =-\delta {\mathit{\beta }}_{\mathsf{A}\mathsf{F}\mathsf{T}}=-\frac{1}{b}{\mathit{\beta }}_{\mathsf{A}\mathsf{F}\mathsf{T}}\end{array}$$

AFT and PH model

• Survivor function for some constants $c>0$ and $r(\mathbf{x})>0$ $$\begin{array}{rl}S(t|{\mathbf{x}}_2)& =S(ct|{\mathbf{x}}_1)\\ S(t|{\mathbf{x}}_2)& =[S(t|{\mathbf{x}}_1){]}^{r({\mathbf{x}}_1)/r({\mathbf{x}}_2)}\\ H(t|{\mathbf{x}}_2)& =[r({\mathbf{x}}_2)/r({\mathbf{x}}_1)]H(t|{\mathbf{x}}_1)\end{array}$$

6.2 Inference for Log-location-scale AFT Models

Likelihood methods

• Data $$\{(y_i,{\delta }_i,{\mathbf{x}}_i),i=1,\dots ,n\}$$

  • Log-lifetime or log-censoring $y_i=\log t_i$

  • Censoring indicator ${\delta }_i=I(\text{ith observation is a failure})$

  • ${\mathbf{x}}_i=(1,x_{i1},\dots ,x_{ip}{)}^{\prime }$ is a vecor of covariates

---

• Assume $Y_i$ follows a location-scale distribution with location parameter $u({\mathbf{x}}_i;\mathit{\beta })$ and scale parameter $b$

• Regression model $$\begin{array}{rl}u({\mathbf{x}}_i;\mathit{\beta })& ={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_p{x}_{ip}={\mathbf{x}}_i^{\prime }\mathit{\beta }\end{array}$$

  • Vector of regression parameters $$\mathit{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_p{)}^{\prime }$$

  • Covariate vector ${\mathbf{x}}_i$ contains both categorical and quantitative variables, and for accurate computation, quantitative variables are centered

---

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

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

• $z_i=\frac{y_i-u({\mathbf{x}}_i;\mathit{\beta })}{b}$

• $u({\mathbf{x}}_i;\mathit{\beta })={\mathbf{x}}_i^{\prime }\mathit{\beta }$

---

Score functions

Elements of $(p+2)$-dimensional vector of score function $$\begin{array}{rl}{U}_j(\mathit{\beta },b)& =\frac{\partial \ell (\mathit{\beta },b)}{\partial {\beta }_j},j=0,1,\dots ,p\\ U_b(\mathit{\beta },b)& =\frac{\partial \ell (\mathit{\beta },b)}{\partial b}\end{array}$$

• Homework: Obtain the expressions of score function (Eq. 6.3.3 and 6.3.4 of textbook)

---

Information matrix

Elements of observed information matrix $$I(\mathit{\beta },b)=-\left[\begin{array}{cc}\frac{{\partial }^2\ell }{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\prime }}& \frac{{\partial }^2\ell }{\partial \mathit{\beta }\partial b}\\ \frac{{\partial }^2\ell }{\partial b\partial {\mathit{\beta }}^{\prime }}& \frac{{\partial }^2\ell }{\partial b^2}\end{array}\right]$$

• Homework: Obtain the expressions of information matrix (Eq. 6.3.5, 6.3.6 and 6.3.7 of textbook)

---

MLEs

$$\begin{array}{r}({\hat{\mathit{\beta }}}^{\prime },\hat{b}{)}^{\prime }={\arg max}_{({\mathit{\beta }}^{\prime },b{)}^{\prime }\in \mathrm{\Theta }}\ell (\mathit{\beta },b)\end{array}$$

• Iterative procedures (e.g. Newton-Raphson method) is used obtain MLE for $\mathit{\beta }$ and $b$

• MLEs $({\hat{\mathit{\beta }}}^{\prime },\hat{b}{)}^{\prime }$ follow a $(p+2)$-variate normal distribution with mean $({\mathit{\beta }}^{\prime },b{)}^{\prime }$ and variance matrix $$\begin{array}{r}\hat{V}=[I(\hat{\mathit{\beta }},\hat{b}){]}^{-1}\end{array}$$

• Large sample based tests and confidence intervals can be obtained using the sampling distribution of $({\hat{\mathit{\beta }}}^{\prime },\hat{b}{)}^{\prime }$

Test of hypothesis

Let ${\mathit{\beta }}^{\prime }=({\mathit{\beta }}_1^{\prime },{\mathit{\beta }}_2^{\prime })$, where ${\mathit{\beta }}_1$ is a $k$-dimensional vector of regression parameters, where $k<p$ $$\begin{array}{r}{H}_0:{\mathit{\beta }}_1={\mathit{\beta }}_1^0\end{array}$$

• Likelihood ratio test statistic $$\begin{array}{r}{\mathrm{\Lambda }}_1=2\ell ({\hat{\mathit{\beta }}}_1,{\hat{\mathit{\beta }}}_2,\hat{b})-2\ell ({\mathit{\beta }}_1^0,{\tilde{\beta }}_2,\tilde{b})\end{array}$$

  • $({\tilde{\mathit{\beta }}}^{\prime },\tilde{b}{)}^{\prime }={\arg max}_{({\mathit{\beta }}_2^{\prime },b{)}^{\prime }\in \mathrm{\Theta }}\ell ({\mathit{\beta }}_1^0,{\mathit{\beta }}_2,b)$

  • Under $H_0$, ${\mathrm{\Lambda }}_1\sim {\chi }_{(k)}^2$

---

Let ${\mathit{\beta }}^{\prime }=({\mathit{\beta }}_1^{\prime },{\mathit{\beta }}_2^{\prime })$, where ${\mathit{\beta }}_1$ is a $k$-dimensional vector of regression parameters, where $k<p$ $$\begin{array}{r}{H}_0:{\mathit{\beta }}_1={\mathit{\beta }}_1^0\end{array}$$

• Wald statistic $$\begin{array}{r}{\mathrm{\Lambda }}_2=({\mathit{\beta }}_1-{\mathit{\beta }}_1^0{)}^{\prime }{V}_{11}^{-1}({\mathit{\beta }}_1-{\mathit{\beta }}_1^0)\end{array}$$

  • $V_{11}=var({\hat{\mathit{\beta }}}_1)$ is a $k\times k$ matrix and $$\begin{array}{r}\hat{V}=[I(\hat{\mathit{\beta }},\hat{b}){]}^{-1}=\left[\begin{array}{cc}{V}_{11}& V_{12}\\ V_{12}^{\prime }& V_{22}\end{array}\right]\end{array}$$

    • Under $H_0$, ${\mathrm{\Lambda }}_1\sim {\chi }_{(k)}^2$

---

Null hypothesis $$\begin{array}{r}{H}_0:{\beta }_j=0\end{array}$$

• Test statistic $$\begin{array}{r}{Z}_j=\frac{{\hat{\beta }}_j}{se({\hat{\beta }}_j)}\end{array}$$

• $100(1-\alpha )\mathrm{\%}$ confidence interval for ${\beta }_j$ $${\hat{\beta }}_j\pm z_{1-\alpha /2}se({\hat{\beta }}_j)$$

• For a small sample, LRT statistic can be used to test the hypothesis and to obtain confidence interval

Quantiles

• The $pth$ quantile of $Y$ given $\mathbf{x}$ $$\begin{array}{r}{y}_p(\mathbf{x})={\mathbf{x}}^{\prime }\mathit{\beta }+b{w}_p\end{array}$$

• Estimate and corresponding SEs of $pth$ quantile
  $$\begin{array}{rl}{\hat{y}}_p(\mathbf{x})& ={\mathbf{x}}^{\prime }\hat{\mathit{\beta }}+\hat{b}{w}_p\\ \text{and using delta method},se({\hat{y}}_p(\mathbf{x}))& ={\mathbf{a}}^{\prime }V\mathbf{a}\end{array}$$

  • $\mathbf{a}=({\mathbf{x}}^{\prime },w_p)$ and $w_p=S_0^{-1}(1-p)$

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

Survival probability

• We are interest to obtain confidence interval for $S(y_0)$, which can be expressed in terms of the parameters of location-scale distribution as $$\begin{array}{rl}S(y_0)& =S_0(\frac{y_0-{\mathbf{x}}^{\prime }\mathit{\beta }}{b})\\ S_0^{-1}(S(y_0))& =\frac{y_0-{\mathbf{x}}^{\prime }\mathit{\beta }}{b}=\psi (\mathbf{x})\end{array}$$

---

• Estimate and the corresponding SE of $\psi (\mathbf{x})$ $$\begin{array}{r}\hat{\psi }(\mathbf{x})=\frac{y_0-{\mathbf{x}}^{\prime }\hat{\mathit{\beta }}}{\hat{b}}\text{and using delta method,}se(\hat{\psi }(\mathbf{x}))=[{\mathbf{a}}^{\prime }V\mathbf{a}{]}^{1/2}\end{array}$$

  • ${\mathbf{a}}^{\prime }=(-1/\hat{b})({\mathbf{x}}^{\prime },\hat{\psi }(\mathbf{x}))$

• $(1-\alpha )100\mathrm{\%}$ confidence interval for $\psi (\mathbf{x})$ $$\hat{\psi }(\mathbf{x})\pm z_{1-\alpha /2}se(\hat{\psi }(\mathbf{x}))$$

---

• Wald-type $(1-\alpha )100\mathrm{\%}$ confidence interval for $S(y_0)$ $$\begin{array}{rl}\hat{\psi }(\mathbf{x})-z_{1-\alpha /2}se(\hat{\psi }(\mathbf{x}))& <\psi (\mathbf{x})<\hat{\psi }(\mathbf{x})+z_{1-\alpha /2}se(\hat{\psi }(\mathbf{x}))\\ L& <\psi (\mathbf{x})<U\\ L& <S_0^{-1}(S(y_0))<U\\ S_0(L)& <S(y_0)<S_0(U)\end{array}$$

6.3 Weibull AFT

• Distributional assumption $$\begin{array}{rl}T(\mathbf{x})=(T|\mathbf{x})& \sim \text{Weib}(\alpha (\mathbf{x}),\delta )\\ Y(\mathbf{x})=(Y|\mathbf{x})=(\log T|\mathbf{x})& \sim \text{EV}(u(\mathbf{x}),b)\end{array}$$

• Regression model for the parameters $$\begin{array}{rl}u(\mathbf{x})& ={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p={\mathbf{x}}^{\prime }\mathit{\beta }\\ \alpha (\mathbf{x})& =\exp ({\mathbf{x}}^{\prime }\mathit{\beta })\end{array}$$

  • $\mathbf{x}=(1,x_1,\dots ,x_p{)}^{\prime }$

  • $\mathit{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_p{)}^{\prime }$

---

• Regression model for the response $$\begin{array}{rl}Y(\mathbf{x})& ={\mathbf{x}}^{\prime }\mathit{\beta }+bZ\end{array}$$

  • $Z\sim \text{EV}(0,1)$

  • $f_0(z)=\exp (z-e^z)$

  • $S_0(z)=\exp (-e^z)$

---

• Log-likelihood function $$\begin{array}{rl}\ell (\mathit{\beta },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({\delta }_i{z}_i-e^{z_i})\end{array}$$
  • $z_i=(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })/b$
• We can now obtain score functions, information matrix, and MLE’s for $\beta$ and $b$ (according to Section 6.2.)

---

• We’ve already seen that the Weibull model implies a proportional hazard model

• It is the only parametric model that is both an AFT model and a Proportional Hazards (PH) model at the same time

Leukemia survival times

• Data on survival times for 33 leukemia patients are available, where survival times are in weeks from diagnosis

• Data on two covariates are also available

  • White blood cell count (WBC) at diagnosis

  • Binary variable AG indicates a positive (AG=1) or negative (AG=0) test related to white blood cell characteristics

---

---

tab6_1

# A tibble: 33 × 5
    time   wbc    AG status   lwbc
   <dbl> <dbl> <int>  <dbl>  <dbl>
 1    65  2.3      1      1  0.833
 2   140  0.75     1      0 -0.288
 3   100  4.3      1      1  1.46 
 4   134  2.6      1      1  0.956
 5    16  6        1      1  1.79 
 6   106 10.5      1      0  2.35 
 7   121 10        1      1  2.30 
 8     4 17        1      1  2.83 
 9    39  5.4      1      1  1.69 
10   121  7        1      0  1.95 
# ℹ 23 more rows

---

• Consider Weibull AFT model with covariates $x_1=AG$ and $x_2=\log (wbc)$ $$\begin{array}{}\text{(6.8)}& \begin{array}{r}Y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+bZ\end{array}\end{array}$$

  • $Z\sim \text{EV}(0,1)$

---

Fit Weibull regression model Equation 6.8 using R

mod62 <- survreg(Surv(time, status) ~ AG + lwbc, 
                 data = tab6_1, dist = "weibull")

mod62E <- survreg(Surv(log(time), status) ~ AG + lwbc, 
                 data = tab6_1, dist = "extreme")

---

MLEs of model parameters

tidy(mod62, conf.int = T) |> 
  mutate(p.value = scales::pvalue(p.value))

# A tibble: 4 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl> <chr>      <dbl>     <dbl>
1 (Intercept)    3.84      0.534     7.19  <0.001     2.79     4.89  
2 AG             1.18      0.427     2.76  0.006      0.340    2.01  
3 lwbc          -0.366     0.150    -2.45  0.014     -0.660   -0.0731
4 Log(scale)     0.112     0.147     0.765 0.444     NA       NA     

---

Fitted model with $x_1=AG$ and $x_2=\log (wbc)$

$$\hat{Y}=3.841+1.177x_1-0.366x_2+\exp (1.119)Z$$

---

• Variance matrix of the estimated parameters

vcov(mod62) %>% round(3)

            (Intercept)     AG   lwbc Log(scale)
(Intercept)       0.286 -0.130 -0.067      0.003
AG               -0.130  0.182  0.016      0.005
lwbc             -0.067  0.016  0.022     -0.005
Log(scale)        0.003  0.005 -0.005      0.021

---

Summary of Weibull AFT model fit

term	estimate	std.error	statistic	p.value	conf.low	conf.high
${\beta }_0$	3.841	0.534	7.188	<0.001	2.794	4.889
${\beta }_1$	1.177	0.427	2.757	0.006	0.340	2.014
${\beta }_2$	-0.366	0.150	-2.449	0.014	-0.660	-0.073
$\log b$	0.112	0.147	0.765	0.444	NA	NA

• AG and WBC have significant effects on leukemia survival times. Positive AG and low WBC count are associated with more prolonged survival

• Since $\log b$ is not significant, i.e. there is not enough evidence to reject $H_0:\log b=1$, exponential AFT model would be appropriate for analyzing this data

---

Interpretations

$$\begin{array}{r}\exp ({\hat{\beta }}_1)=\exp (1.177)=3.246\end{array}$$

• A specific quantile (say median) lifetime of a patient with a positive AG value (i.e. $x_1=1$) is 3.2 times that of a patient with a negative AG (i.e. $x_1=0$) value provided WBC value remains constant

• Note this interpretation is true for any quantile (Why?)

---

$$\begin{array}{r}\exp ({\hat{\beta }}_2)=\exp (-0.366)=0.693\end{array}$$

• A specific quantile (say median) lifetime of a patient decreases 30.7 percent with one unit increase of log(WBC) [or 2718 unit increase of true WBC count] provided AG value remains constant

---

Fitted values

augment(mod62, type.predict = "response") |> 
  select(1:4) |> 
  slice(1:3)

# A tibble: 3 × 4
  `Surv(time, status)`    AG   lwbc .fitted
                <Surv> <int>  <dbl>   <dbl>
1                  65      1  0.833   111. 
2                 140+     1 -0.288   168. 
3                 100      1  1.46     88.6

---

augment(mod62, type.predict = "link") |> 
  mutate(.fittedE = exp(.fitted)) |> 
  select(2:4, .fittedE) |> 
  slice(1:3)

# A tibble: 3 × 4
     AG   lwbc .fitted .fittedE
  <int>  <dbl>   <dbl>    <dbl>
1     1  0.833    4.71    111. 
2     1 -0.288    5.12    168. 
3     1  1.46     4.48     88.6

• Estimate for a subject with $AG=1$ and $\log (wbc)=.833$

$$\hat{u}={\hat{\beta }}_0+{\hat{\beta }}_1(1)+{\hat{\beta }}_2(.833)=(3.841)+(1.177)(1)+(-0.366)(.833)=4.713$$

---

#predict(object = mod62, newdata = tibble(AG = 1, lwbc = .833), 
#        predict = "response")
augment(x = mod62, newdata = tibble(AG = 1, lwbc = .833), 
        type.predict = "response")

# A tibble: 1 × 4
     AG  lwbc .fitted .se.fit
  <dbl> <dbl>   <dbl>   <dbl>
1     1 0.833    111.    41.3

• Estimate for a subject with $AG=1$ and $\log (wbc)=.833$

$$\begin{array}{rl}\hat{\alpha }=\exp ({\hat{\beta }}_0+{\hat{\beta }}_1(1)+{\hat{\beta }}_2(.833))& =\exp ((3.841)+(1.177)(1)+(-0.366)(.833))\\ & =\exp (4.713)=111.399\end{array}$$

---

LRT

• Likelihood ratio tests for $H_0:{\beta }_1=0$ $$\begin{array}{r}{\mathrm{\Lambda }}_1(0)=2\ell ({\hat{\beta }}_0,{\hat{\beta }}_1,{\hat{\beta }}_2,\log \hat{b})-2\ell ({\tilde{\beta }}_0,0,{\tilde{\beta }}_2,\log \tilde{b})\end{array}$$

  • ${\mathrm{\Lambda }}_1(0)\sim {\chi }_{(1)}^2$

• The corresponding $Z$ statistic $$Z=\text{sign}({\hat{\beta }}_1){\mathrm{\Lambda }}_1^{1/2}\sim \mathcal{N}(0,1)$$

---

Estimate of model parameters under $H_0:{\beta }_1=0$

# mod62a <- update(mod62, formula = . ~ . - AG)
mod62a <- survreg(Surv(time, status) ~  lwbc, 
                 data = tab6_1, dist = "weibull")
tidy(mod62a)

# A tibble: 3 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    4.85      0.500      9.71 2.67e-22
2 lwbc          -0.500     0.165     -3.03 2.41e- 3
3 Log(scale)     0.222     0.146      1.52 1.28e- 1

---

LRTa <- anova(mod62a, mod62)

Terms	Resid. Df	-2*LL	Df	Deviance	Pr(>Chi)
lwbc	30	271.931	NA	NA	NA
AG + lwbc	29	265.013	1	6.918	0.009

• ${\mathrm{\Lambda }}_1(0)=6.918\Rightarrow Z=2.63$

---

Comparison between Wald- and LRT-type Z statistics

term	estimate	Wald	LRT
${\beta }_0$	3.841	7.188	NA
${\beta }_1$	1.177	2.757	2.63
${\beta }_2$	-0.366	-2.449	-2.46
$\log b$	0.112	0.765	NA

---

Quantiles

$${\hat{y}}_p={\mathbf{x}}^{\prime }\hat{\mathit{\beta }}+\log (-\log (1-p))\hat{b}$$

• Consider a subject with covariate values $x_1=1$ and $x_2=\log (10)$, the linear predictor ${\mathbf{x}}^{\prime }\hat{\mathit{\beta }}$ $$\begin{array}{rl}\hat{u}=x{x}^{\prime }\hat{\mathit{\beta }}& ={\hat{\beta }}_0+{\hat{\beta }}_1+\log (10){\hat{\beta }}_2=4.175\end{array}$$

---

• Median survival time of the patient with covariate values $x_1=1$ and $x_2=\log (10)$ $$\begin{array}{rl}{\hat{y}}_{.50}& =4.175+(-0.367)(1.119)=3.765\\ {\hat{t}}_{.50}& =\exp (3.765)=43.163\text{weeks}\end{array}$$

• Homework: Obtain a 95% confidence interval of the median survival time of a patient with covariate values $x_1=1$ and $x_2=\log (10)$

---

Survival probability

$$S(y_0)=\exp \{-\exp [(y_0-{\mathbf{x}}^{\prime }\hat{\mathit{\beta }})/\hat{b}]\}$$ For a patient with covariate values $x_1=1$ and $x_2=\log (10)$, obtain $S(\log 10)$ $$S(\log 10)=\exp \{-\exp [(\log 10-4.175)/1.119]\}=0.816$$

• Homework: Obtain the 95% CI for $S(\log 10)$

6.4 Log-normal AFT

• Distributional assumption $$\begin{array}{rl}T(\mathbf{x})=(T|\mathbf{x})& \sim \text{log-Norm}(\mu (\mathbf{x}),{\sigma }^2)\\ Y(\mathbf{x})=(Y|\mathbf{x})=(\log T|\mathbf{x})& \sim \mathcal{N}(\mu (\mathbf{x}),{\sigma }^2)\end{array}$$

---

• Regression model for the parameters $$\begin{array}{rl}\mu (\mathbf{x})& ={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p={\mathbf{x}}^{\prime }\mathit{\beta }\end{array}$$

• Regression model for the response $$\begin{array}{rl}Y(\mathbf{x})& ={\mathbf{x}}^{\prime }\mathit{\beta }+\sigma Z\end{array}$$

  • $Z\sim \mathcal{N}(0,1)$

  • $f_0(z)=\varphi (z)$

  • $S_0(z)=1-\mathrm{\Phi }(z)$

Times to pulmonary exacerbation

• Patients with cystic fibrosis are susceptible to an accumulation of mucus in the lungs, which leads to pulmonary exacerbation and deterioration of lung function

• A clinical trial was conducted to investigate the efficacy of the new drug DNase-1

  • Subjects are randomly assigned to a new treatment or a placebo

• Time of interest is the time to first exacerbation after randomization, and data on fev (forced expiratory volume at the time of randomization) are also measured

---

---

# A tibble: 761 × 13
      id   trt  time   fev  inst entry.dt   end.dt     ivstart ivstop time0
   <int> <int> <dbl> <dbl> <int> <date>     <date>       <dbl>  <dbl> <dbl>
 1     1     1   168  28.8     1 1992-03-20 1992-09-04      NA     NA   168
 2     2     1   169  64       1 1992-03-24 1992-09-09      NA     NA   169
 3     3     0    65  67.2     1 1992-03-24 1992-09-08      65     75   168
 4     4     1   168  57.6     1 1992-03-26 1992-09-10      NA     NA   168
 5     5     0   171  57.6     1 1992-03-24 1992-09-11      NA     NA   171
 6     6     1   166  25.6     1 1992-03-27 1992-09-09      NA     NA   166
 7     7     0   168  86.4     1 1992-03-27 1992-09-11      NA     NA   168
 8     8     0    90  32       1 1992-03-28 1992-09-10      90    104   166
 9     9     1   169  86.4     2 1992-02-27 1992-08-14      NA     NA   169
10    10     0     8  28.8     2 1992-03-06 1992-08-22       8     22   169
# ℹ 751 more rows
# ℹ 3 more variables: status <dbl>, fevm <dbl>, visit <int>

---

• Assume survival time $T(\mathbf{x})$ follows a log-normal distribution with scale parameter $\alpha (\mathbf{x})$ and shape parameter $\delta$

• Consider following AFT model for log survival time $$Y(\mathbf{x})={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\sigma Z$$

  • $Z\sim \mathcal{N}(0,1)$

  • $x_1=I(\text{trt}=1)$

  • $x_2=\text{fev}-\text{mean(fev)}$

---

• R codes for fitting the AFT model

mod63a <- survreg(Surv(log(time), status) ~ trt + fevm, 
                  dist = "gaussian",
                  data = tab1_4) 

---

tidy(mod63a)

# A tibble: 4 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   5.09     0.0684      74.4  0       
2 trt           0.336    0.0951       3.53 4.19e- 4
3 fevm          0.0159   0.00197      8.09 5.91e-16
4 Log(scale)    0.137    0.0408       3.36 7.84e- 4

---

• AFT model $$\begin{array}{rl}Y& ={\mathbf{x}}^{\prime }\mathit{\beta }+bZ\Rightarrow T=\exp ({\mathbf{x}}^{\prime }\mathit{\beta })\exp (bZ)\end{array}$$

  • ${\mathbf{x}}^{\prime }\mathit{\beta }={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p$

• For a binary predictor $x_j$ $$T=\exp ({\mathbf{x}}^{\prime }\mathit{\beta })\exp (bZ)=\{\begin{array}{ll}\exp (bZ)& \text{for control}\\ \exp ({\beta }_j)\exp (bZ)& \text{for treatment}\end{array}$$

• It can be shown that $$T_{trt}=\exp (\beta )T_{control}$$

---

• ${\beta }_{trt}=0.336\Rightarrow \exp ({\beta }_{trt})=1.399$

  • Treatment increases the time to first pulmonary exacerbation by about 40% compared to the control when fev is fixed

• ${\beta }_{fev}=0.016\Rightarrow \exp ({\beta }_{fev})=1.016$

  • One-unit increase in fev results about 2% increase in lifetime provided treatment is constant

---

6.5 Log-logistic AFT

• Distributional assumptions $$\begin{array}{rl}T(\mathbf{x})=(T|\mathbf{x})& \sim \text{log-logistic}(\alpha (\mathbf{x}),\beta )\\ Y(\mathbf{x})=(Y|\mathbf{x})=(\log T|\mathbf{x})& \sim \text{logistic}(u(\mathbf{x}),b)\end{array}$$

---

• Regression model for the parameters $$\begin{array}{rl}u(\mathbf{x})& ={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_p{x}_p={\mathbf{x}}^{\prime }\mathit{\beta }\\ \alpha (\mathbf{x})& =e^{u(\mathbf{x})}\end{array}$$

• Regression model for the response $$\begin{array}{rl}Y(\mathbf{x})& ={\mathbf{x}}^{\prime }\mathit{\beta }+bZ\end{array}$$

  • $Z\sim \text{Logistic}(0,1)$

  • $f_0(z)=e^z[1+e^z{]}^{-2}$

  • $S_0(z)=[1+e^z{]}^{-1}$

---

• Lifetime distribution $$T(xx)\sim \text{Log-Logistic}(\alpha (xx),\delta )$$

• The survivor function $$S(t|xx)=\frac{1}{1+(t/\alpha (xx){)}^{\delta }}\Rightarrow \frac{1-S(t|xx)}{S(t|xx)}=(t/\alpha (xx){)}^{\delta }$$

  • $(t/\alpha (xx){)}^{\delta }\to$ the odds of failure at time $t$ for a subject with covariate vector $xx$

---

• For two subjects with covariate vectors $x{x}_1$ and $x{x}_2$ $$\frac{[1-S(t|x{x}_2)]/S(t|x{x}_2)}{[1-S(t|x{x}_1)]/S(t|x{x}_1)}=[\frac{\alpha (x{x}_1)}{\alpha (x{x}_2)}{]}^{\delta },\text{independent of }t$$

• A model of the form $$\frac{1-S(t|xx)}{S(t|xx)}=(t/\alpha (xx){)}^{\delta }\Rightarrow \log \frac{1-S(t|xx)}{S(t|xx)}=\delta \log (t)-\delta \log \alpha (xx)$$ is known as the proportional odds model

---

• Consider a model $\log \alpha (x)={\beta }_0+{\beta }_{1}x$ $$\frac{[1-S(t|x=1)]/S(t|x=1)}{[1-S(t|x=0)]/S(t|x=0)}=e^{-\delta {\beta }_1}=e^{-{\beta }_1^{\star }}$$

  • The odds of failure at time $t$ for a subject with $x=1$ is $\exp (-{\beta }^{\star })$ times that of the odds of failure for a subject with $x=0$

Times to pulmonary exacerbation

• R codes for fitting AFT model

mod63b <- survreg(Surv(log(time), status) ~ trt + fevm, 
                  dist = "logistic",
                  data = tab1_4) 

---

tidy(mod63b)

# A tibble: 4 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   5.08     0.0600      84.6  0       
2 trt           0.293    0.0861       3.41 6.55e- 4
3 fevm          0.0145   0.00181      8.00 1.20e-15
4 Log(scale)   -0.489    0.0466     -10.5  8.08e-26

---

• ${\beta }_{trt}=0.293\Rightarrow \exp ({\beta }_{trt})=1.341$

  • Treatment increases the time to first pulmonary exacerbation by about 34% compared to the control when fev is fixed

• ${\beta }_{fev}=0.014\Rightarrow \exp ({\beta }_{fev})=1.015$

  • One-unit increase in fev results in a 1.5% increase in lifetime provided treatment is constant

---

• Interpret the treatment effect in terms of odds of failure $$\frac{[1-S(t|trt=1,fev=x)]/S(t|trt=1,fev=x)}{[1-S(t|trt=0,fev=x)]/S(t|trt=0,fev=x)}=\exp (-\hat{\delta }{\hat{\beta }}_1)=0.62$$

  • $\hat{\delta }=\exp (-\log \hat{b})=\exp (0.489)=1.631$

  • The odds of failure is 38% lower in the treatment group compared to the control group provided fev value is fixed

---

---

Comparison between normal and logistic regression models in analysing time to pulmonary exacerbation data

term	est	se	est	se
(Intercept)	5.093	0.068	5.078	0.060
trt	0.336	0.095	0.293	0.086
fevm	0.016	0.002	0.014	0.002
Log(scale)	0.137	0.041	-0.489	0.047

Other regression models

• Additive hazards model $$h(t|xx)=h_0(t;\alpha )+r(xx;\beta \beta )$$

6.6 Graphical methods and model assessment

• Graphical methods are helpful in summarizing information and suggesting possible models

• These methods also provide ways to check assumptions concerning the form of a lifetime distribution and its relationship to covariates

• Exploratory analysis of a lifetime distribution given covariates would helpful to select the appropiate Model for the analysis

---

• For a single quantitative covariate, a plot of lifetime or log-lifetime against the covariate or a function of it could indicate the nature of the relationship between lifetime and the covariate

• If the proportion of censoring is small, such a plot would be helpful, different symbols can be used in those plots for censored and failure times

• When more than one quantitative covariate and light censoring, one can consider grouping individuals so that within a group, individuals will have similar values of important covariates

• Let there are $J$ such groups and ${\hat{S}}_j$ is the Kaplan-Meier estimate for the group $j=1,\dots ,J$

---

AFT model $$S(t|\mathbf{x})=S_0[\frac{\log t-u(\mathbf{x})}{b}]$$

• If $u(\mathbf{x})$ is approximately constant for individuals within each group $j=1,\dots ,J$, and if an AFT model is appropriate, the plots of $$\log [-\log S(t|\mathbf{x})]\text{ vs}\log t$$ should be roughly parallel in horizontal direction ($\log t$)

---

Proportional hazards model $$S(t|\mathbf{x})=[S_0(t){]}^{r(\mathbf{x})}$$

• If $r(\mathbf{x})$ is approximately constant for individuals within each group $j=1,\dots ,J$, and if a proportional hazards model is appropriate, the plots of $$\log [-\log S(t|\mathbf{x})]\text{vs}\log t$$ should be roughly parallel in vertical direction

---

• If the plots of $\log [-\log S(t|\mathbf{x})]$ vs $\log t$ is roughly linear then Weibull models are suggested

• In addition to linear, if the plots are parallel, then Weibull models with a constant shape parameter are suggested, in that case, both AFT and PH models can be considered

• Statistical analysis of data is an iterative process involving exploration, model fitting, and model assessment

References

Collett, David. 2015. Modelling Survival Data in Medical Research. Chapman; Hall/CRC. https://doi.org/10.1201/b18041.

Cox, David R. 1972. “Regression Models and Life-Tables.” Journal of the Royal Statistical Society: Series B (Methodological) 34 (2): 187–202.