Chapter 5 Statistical Methods in Survival Analysis

5.1 Likelihood Theory

There are always merits in obtaining raw data (that is, exact individual failure times), as opposed to grouped data (counts of the number of failures over prescribed time intervals). Given raw data, we can always construct grouped data, but the converse is typically not true; therefore, we limit discussion in this chapter to the raw data case.

The random variable T has denoted a random lifetime in previous chapter. So it is natural to use [latex]T_1, \, T_2, \, \ldots, \, T_n[/latex] to denote a random sample of n such lifetimes, where n is the number of items on test. When specific values are given for realizations of such lifetimes, which is typically the case from this point forward, they are denoted by [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex]. In other words, [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] are the experimental values of the mutually independent and identically distributed random variables [latex]T_1, \, T_2, \, \ldots, \, T_n[/latex]. The associated ordered observations, or order statistics, are denoted by [latex]t_{(1)}, \, t_{(2)}, \, \ldots, \, t_{(n)}[/latex].

The Greek letter θ is often used to denote a generic unknown parameter. We will refer to [latex]\hat \theta[/latex] in the abstract as a point estimator; when [latex]\hat \theta[/latex] assumes a specific numeric value, it will be referred to as a point estimate. The probability distribution of a statistic is referred to as a sampling distribution.

Assume that there is a single unknown parameter θ in the probability model for T. Assume further that the data values [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] are mutually independent and identically distributed random variables. The joint probability density function of the data values is the product of the marginal probability density functions of the individual observations:

$$\begin{array}{l}L(t_1,t_2,\dots ,t_n,\theta )=\prod _{i=1}^{n}f(t_i;\theta ).\end{array}$$

This function is the likelihood function. In order to simplify the notation, the likelihood function is often written as simply

$$\begin{array}{l}L(\theta )=\prod _{i=1}^{n}f(t_i).\end{array}$$

The maximum likelihood estimator of θ, which is denoted by [latex]\hat \theta[/latex], is the value of θ that maximizes [latex]L(\theta)[/latex].

The next example reviews the associated notions of the log likelihood function, score vector, maximum likelihood estimator, Fisher information matrix, and observed information matrix for a two-parameter lifetime model. We assume for now that there are no censored observations in the data set; all of the failure times are observed.

In some cases, it is possible to find the exact distribution of a pivotal quantity which results in exact statistical inference (that is, constructing exact confidence intervals and performing exact hypothesis tests). It is more often the case that exact statistical inference is not possible, and asymptotic properties associated with the likelihood function must be relied on for approximate inference. The next section reviews some asymptotic properties that arise in likelihood theory. When a large data set of lifetimes is available, these properties often lead to approximate statistical methods of inference.

5.2 Asymptotic Properties

When the number of items on test n is large, there are some asymptotic results concerning the likelihood function that are useful for constructing confidence intervals and performing hypothesis tests associated with a vector of p unknown parameters [latex]\boldsymbol{\theta} = \left( \theta_1 , \, \theta_2 , \, \ldots , \, \theta_p \right) ^ \prime[/latex]. As indicated in the example in the last section, the [latex]p \times 1[/latex] score vector [latex]{\boldsymbol U} (\boldsymbol{\theta})[/latex] has elements

$$\begin{array}{l}{U}_i(\mathit{\theta })=\frac{\partial \ln L(\mathit{t},\mathit{\theta })}{\partial {\theta }_i}=\frac{\partial }{\partial {\theta }_i}\sum _{j=1}^n\ln f(t_j,\mathit{\theta })\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, p[/latex]. Therefore, each element of the score vector is a sum of mutually independent random variables, and, when n is large, the elements of [latex]\boldsymbol{U}( \boldsymbol{\theta})[/latex] are asymptotically normally distributed by the central limit theorem. More specifically, the score vector [latex]\boldsymbol{U} (\boldsymbol{\theta})[/latex] is asymptotically normal with population mean [latex]{\bf 0}[/latex] and variance–covariance matrix [latex]I(\boldsymbol{\theta})[/latex], where [latex]I(\boldsymbol{\theta})[/latex] is the Fisher information matrix. This means that when the true value for the parameter vector is [latex]\boldsymbol{\theta}_0[/latex] then

$$\begin{array}{l}{\mathit{U}}^{\mathrm{\prime <!--\; \prime \; -->}}({\mathit{\theta }}_0)I({\mathit{\theta }}_0{)}^{-1}\mathit{U}({\mathit{\theta }}_0)\end{array}$$

is asymptotically chi-square with p degrees of freedom. This can be used to determine confidence intervals and perform hypothesis tests with respect to [latex]\boldsymbol{\theta}[/latex].

The maximum likelihood estimator for the parameter vector [latex]\hat{\boldsymbol{\theta}}[/latex] can also be used for confidence intervals and hypothesis testing. Since [latex]\hat{\boldsymbol{\theta}}[/latex] is asymptotically normal with population mean [latex]\boldsymbol{\theta}[/latex] and variance–covariance matrix [latex]I^{-1} (\boldsymbol{\theta})[/latex], when [latex]\boldsymbol{\theta} = \boldsymbol{\theta}_0[/latex],

$$\begin{array}{l}{(\hat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{\prime <!--\; \prime \; -->}}I({\mathit{\theta }}_0)(\hat{\mathit{\theta }}-{\mathit{\theta }}_0)\end{array}$$

is also asymptotically chi-square with p degrees of freedom. Two statistics that are asymptotically equivalent to this statistic that can be used to estimate the value of the chi-square random variable are

$$\begin{array}{l}{(\hat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{\prime <!--\; \prime \; -->}}I\left(\hat{\mathit{\theta }}\right)(\hat{\mathit{\theta }}-{\mathit{\theta }}_0)\end{array}$$

and

$$\begin{array}{l}{(\hat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{\prime <!--\; \prime \; -->}}O\left(\hat{\mathit{\theta }}\right)(\hat{\mathit{\theta }}-{\mathit{\theta }}_0).\end{array}$$

A third asymptotic result involves the likelihood ratio statistic

$$\begin{array}{l}-2[\ln L(\mathit{\theta })-\ln L(\hat{\mathit{\theta }})]=-2\ln \left[\frac{L(\mathit{\theta })}{L(\hat{\mathit{\theta }})}\right],\end{array}$$

which is asymptotically chi-square with p degrees of freedom. The conditions necessary for these asymptotic properties to apply are cited at the end of the chapter.

These three asymptotic results are summarized in the result below, where the a above the ∼ is shorthand for “asymptotically distributed.”

All of the statistical methods developed thus far have assumed that we are able to observe all n of the items on test fail. The associated lifetimes are denoted by [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex]. Although this is ideal and might be the case in some settings, a short testing time or items with long lifetimes might result in some items that survive the test. The lifetimes of the items which do not fail during the test are known as right-censored observations. The lifetimes of these items are not observed, but are known to exceed the time at which the test is concluded. If a decision concerning the acceptability of the items must be made with some of the items still operating at the end of the test, then a statistical model must be formulated to account for the unobserved lifetimes of these items. The next section introduces the important topic of censoring, which is pervasive in survival analysis.

5.3 Censoring

Censoring occurs in lifetime data sets when only an upper or lower bound on the lifetime is known. Censoring occurs frequently in lifetime data sets because it is often impossible or impractical to observe the lifetimes of all the items on test. A data set for which all failure times are known is called a complete data set. Figure 5.1 illustrates a complete data set of [latex]n = 5[/latex] items placed on test simultaneously at time [latex]t = 0[/latex], where the [latex]\times[/latex]‘s denote failure times. Consider the two endpoints of each of the horizontal segments. It is critical to provide an unambiguous definition of the time origin (for example, the time a product is purchased or the time a cancer is diagnosed). Likewise, failure must be defined in an unambiguous fashion. This is easier to define for a light bulb or a fuse than for a ball bearing or a sock. Outside of a reliability setting, a data set of lifetimes is often generically referred to as time-to-event data, corresponding to the time between the time origin and the event of interest. A censored observation occurs when only a bound is known on the time of failure. If a data set contains one or more censored observations, it is called a censored data set.

The most common type of censoring is known as right censoring. In a right-censored data set, one or more items have only a lower bound known on their lifetime. The term sample size is now vague. From this point forward, we use n to denote the number of items on test and use r to denote the number of observed failures. In an industrial life testing situation, for example, [latex]n = 12[/latex] cell phones are put on a continuous, rigorous life test on January 1, and [latex]r = 3[/latex] of the cell phones have failed by December 31. These failed cell phones are discarded upon failure. The remaining [latex]n - r = 12 - 3 = 9[/latex] cell phones that are still operating on December 31 have lifetimes that exceed 365 days, and are therefore right-censored observations. Right censoring is not limited to just reliability applications. In a medical study in which T is the survival time after the diagnosis of a particular type of cancer, for example, a patient can either (a) still be alive at the end of a study, (b) die of a cause other than the particular type of cancer, constituting a right-censored observation, or (c) lose contact with the study (for example, if they leave town), constituting a right-censored observation.

Three special cases of right censoring are common in survival analysis. The first is Type II or order statistic censoring. As shown in Figure 5.2, this corresponds to terminating a study upon one of the ordered failures. The diagram corresponds to a set of [latex]n = 5[/latex] items placed on a test simultaneously at time [latex]t = 0[/latex]. The test is terminated when [latex]r = 3[/latex] failures are observed. Time advances from left to right in Figure 5.2 and the failure of the first item (corresponding to the third ordered observed failure) terminates the test. The lifetimes of the third and fourth items are right censored. Observed failure times are indicated by an [latex]{\large{\times}}[/latex] and right-censoring times are indicated by a [latex]\, \Large{\circ}[/latex]. In Type II censoring, the time to complete the test is random.

The second special case is Type I or time censoring. As shown in Figure 5.3, this corresponds to terminating the study at a particular time. The diagram shows a set of [latex]n = 5[/latex] items placed on a test simultaneously at [latex]t = 0[/latex] that is terminated at the time indicated by the dotted vertical line. For the realization illustrated in Figure 5.3, there are [latex]r = 4[/latex] observed failures. In Type I censoring, the number of failures r is random.

Finally, random censoring occurs when individual items are withdrawn from the test at any time during the study. Figure 5.4 illustrates a realization of a randomly right-censored life test with [latex]n = 5[/latex] items on test and [latex]r = 2[/latex] observed failures. It is usually assumed that the failure times and the censoring times are mutually independent random variables and that the probability distribution of the censoring times does not involve any unknown parameters from the failure time distribution. In other words, in a randomly censored data set, items cannot be more or less likely to be censored because they are at unusually high or low risk of failure.

Although other types of censoring exist, such as left censoring and interval censoring, the focus of this chapter will be on right censoring because it is the most common type of censoring. In the case of right censoring, the ratio [latex]r / n[/latex] is the fraction of items which are observed to fail. When [latex]r / n[/latex] is close to one, the data set is referred to as a lightly censored data set; when [latex]r / n[/latex] is close to zero, the data set is referred to as a heavily censored data set. In the reliability setting, many data sets are heavily censored because the items have long lifetimes. In the biomedical setting, certain cancers have long remission times, resulting in heavily censored data sets.

Of the following three approaches to handling the problem of censoring, only one is both valid and practical. The first approach is to ignore all the censored values and to perform analysis only on those items that were observed to fail. Although this simplifies the mathematics involved, it is not a valid approach. If, for example, this approach is used on a right-censored data set, the analyst is discarding the right-censored values, and these are typically the items that have survived the longest. In this case, the analyst arrives at an overly pessimistic result concerning the lifetime distribution because the best items (that is, the right-censored observations) have been excluded from the analysis. A second approach is to wait for all the right-censored observations to fail. Although this approach is valid statistically, it is not practical. In an industrial setting, waiting for the last light bulb to burn out or the last machine to fail may take so long that the product being tested will not get to market in time. In a medical setting, waiting for the last patient to die from a particular disease may take decades. For these reasons, the proper approach is to handle censored observations probabilistically, including the censored values in the likelihood function.

The likelihood function for a censored data set can be written in several different equivalent forms. Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be mutually independent observations denoting lifetimes sampled randomly from a population. The corresponding right-censoring times are denoted by [latex]c_1, \, c_2, \, \ldots, \, c_n[/latex]. The t_i and c_i values are assumed to be independent, for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. In the case of Type I right censoring, [latex]c_1 = c_2 = \cdots = c_n = c[/latex]. The set U contains the indexes of the items that are observed to fail during the test (that is, the uncensored observations):

$$\begin{array}{l}U=\{i|t_i\le c_i\}.\end{array}$$

The set C contains the indexes of the items whose failure time exceeds the corresponding censoring time (that is, those that are right censored):

$$\begin{array}{l}C=\{i|t_i>c_i\}.\end{array}$$

This notation, along with an important notion known as alignment, are illustrated in the next example.

The usual form for right-censored lifetime data is given by the pairs [latex](x_i, \, \delta_i)[/latex], where [latex]x_i = \min \{t_i, \, c_i\}[/latex] and δ_i is a censoring indicator variable:

$$\begin{array}{l}{\delta }_i=\{\begin{array}{lll}0& & t_i>c_i\\ 1& & t_i\le c_i\end{array}\end{array}$$

for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. The [latex](x_i, \, \delta_i)[/latex] pairs can be reconstructed from the [latex](t_i, \, c_i)[/latex] pairs and vice versa. Hence, δ_i is 1 if the failure of item i is observed and 0 if the failure of item i is right censored, and x_i is the failure time (when [latex]\delta_i = 1[/latex]) or the censoring time (when [latex]\delta_i = 0[/latex]). For the vector of unknown parameters [latex]\boldsymbol{\theta} = (\theta_1, \, \theta_2, \, \ldots, \, \theta_p) ^ \prime[/latex], ignoring a constant factor, the likelihood function is

$$\begin{array}{l}L(\mathit{x},\mathit{\theta })=\prod _{i=1}^{n}f(x_i,\mathit{\theta }{)}^{{\delta }_i}S(x_i,\mathit{\theta }{)}^{1-{\delta }_i}=\prod _{i\in U}f(t_i,\mathit{\theta })\prod _{i\in C}S(c_i,\mathit{\theta })\end{array}$$

where [latex]S (c_i, \, \boldsymbol{\theta})[/latex] is the survivor function of the population distribution with parameters [latex]\boldsymbol{\theta}[/latex] evaluated at censoring time c_i, [latex]i \, \in \, C[/latex]. The reason that the survivor function is the appropriate term in the likelihood function for a right-censored observation is that [latex]S(c_i, \, \boldsymbol{\theta})[/latex] is the probability that item i survives to c_i. The log likelihood function is

$$\begin{array}{l}\ln L(\mathit{x},\mathit{\theta })=\sum _{i\in U}\ln f(t_i,\mathit{\theta })+\sum _{i\in C}\ln S(c_i,\mathit{\theta }),\end{array}$$

or

$$\begin{array}{l}\ln L(\mathit{x},\mathit{\theta })=\sum _{i\in U}\ln f(x_i,\mathit{\theta })+\sum _{i\in C}\ln S(x_i,\mathit{\theta }).\end{array}$$

Since the probability density function is the product of the hazard function and the survivor function, the log likelihood function can be simplified to

$$\begin{array}{l}\ln L(\mathit{x},\mathit{\theta })=\sum _{i\in U}\ln h(x_i,\mathit{\theta })+\sum _{i\in U}\ln S(x_i,\mathit{\theta })+\sum _{i\in C}\ln S(x_i,\mathit{\theta })\end{array}$$

or

$$\begin{array}{l}\ln L(\mathit{x},\mathit{\theta })=\sum _{i\in U}\ln h(x_i,\mathit{\theta })+\sum _{i=1}^n\ln S(x_i,\mathit{\theta }),\end{array}$$

where the second summation now includes all n items on test. Finally, to write the log likelihood in terms of the hazard and cumulative hazard functions only,

$$\begin{array}{l}\ln L(\mathit{x},\mathit{\theta })=\sum _{i\in U}\ln h(x_i,\mathit{\theta })-\sum _{i=1}^{n}H(x_i,\mathit{\theta }),\end{array}$$

since [latex]H(t) = -\ln \, S(t)[/latex]. The choice of which of these three expressions for the log likelihood to use for a particular distribution depends on the particular forms of [latex]S(t)[/latex], [latex]f(t)[/latex], [latex]h(t)[/latex], and [latex]H(t)[/latex]. In other words, one of the distribution representations may possess a mathematical form that is advantageous over the others.

The next example will use the last version of the log likelihood function to find a maximum likelihood estimator and an asymptotically exact confidence interval for an unknown parameter.

To summarize the material introduced so far in this chapter, point estimators are statistics calculated from a data set to estimate an unknown parameter. Confidence intervals reflect the precision of a point estimator. The most common technique for determining a point estimator for an unknown parameter is maximum likelihood estimation, which involves finding the parameter value(s) that make the observed data values the most likely. The maximum likelihood estimators are usually found by using calculus to maximize the log likelihood function. Most population lifetime distributions do not have exact confidence intervals for unknown parameters, so the asymptotic properties of the likelihood function can be used to generate approximate confidence intervals for unknown parameters. Finally, many data sets in reliability are censored, which means that only a bound is known on the lifetime for one or more of the data values. The most common censoring mechanism is known as right censoring, where only a lower bound on the lifetime is known. The number of items on test is denoted by n and the number of observed failures is denoted by r.

The next section applies the techniques developed so far in this chapter to the exponential distribution.

5.4 Exponential Distribution

The exponential distribution is popular due to its tractability for parameter estimation and inference. The exponential distribution can be parameterized by either its population rate λ or its population mean [latex]\mu = 1 / \lambda[/latex]. Using the rate to parameterize the distribution, the survivor, density, hazard, and cumulative hazard functions are

$$\begin{array}{l}S(t,\lambda )=e^{-\lambda t}f(t,\lambda )=\lambda e^{-\lambda t}h(t,\lambda )=\lambda H(t,\lambda )=\lambda t\end{array}$$

for [latex]t \ge 0[/latex]. Note that the unknown parameter λ has been added as an argument in these lifetime distribution representations because it is now also an argument in the likelihood function and is estimated from data.

All the analysis in this and subsequent sections assumes that a random sample of n items from a population has been placed on a test and subjected to typical environmental conditions. Equivalently, [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] are independent and identically distributed random lifetimes from a particular population distribution (exponential in this section). As with all statistical inference, care must be taken to ensure that a random sample of lifetimes is collected. Consequently, random numbers should be used to determine which n items to place on test. In a reliability setting, laboratory conditions should adequately mimic field conditions. Only representative items should be placed on test because items manufactured using a previous design may have a different failure pattern than those with the current design. This is more difficult in a biomedical setting because of inherent differences between patients.

Four classes of data sets (complete, Type II right censored, Type I right censored, and randomly right censored) are considered in separate subsections. In all cases, n is the number of items placed on test and r is the number of observed failures.

5.4.3 Type I Censored Data Sets

The analysis for Type I censored data sets is similar to that for the Type II censoring case. The test is terminated at time c. The censoring times for each item on test are the same: [latex]c_1 = c_2 = \cdots= c_n= c[/latex]. The number of observed failures, r, is a random variable. The total time on test in this case is

$$\begin{array}{l}\sum _{i=1}^n{x}_i=\sum _{i\in U}{t}_i+\sum _{i\in C}{c}_i=\sum _{i=1}^r{t}_{(i)}+(n-r)c.\end{array}$$

As before, the log likelihood function is

$$\begin{array}{l}\ln L(\lambda )=\sum _{i\in U}\ln h(x_i,\lambda )-\sum _{i=1}^{n}H(x_i,\lambda )=r\ln \lambda -\lambda \sum _{i=1}^n{x}_i,\end{array}$$

and the score statistic is

$$\begin{array}{l}U(\lambda )=\frac{r}{\lambda }-\sum _{i=1}^n{x}_i.\end{array}$$

The maximum likelihood estimator for [latex]r > 0[/latex] is

$$\begin{array}{l}\hat{\lambda }=\frac{r}{\sum _{i=1}^n{x}_i},\end{array}$$

the information matrix is

$$\begin{array}{l}I(\lambda )=\frac{r}{{\lambda }^2},\end{array}$$

and the observed information matrix is

$$\begin{array}{l}O(\hat{\lambda })=\frac{r}{{\hat{\lambda }}^2}.\end{array}$$

The functional form of the maximum likelihood estimator is identical to the Type II censoring case. For identical values of r, Type I censoring has a larger total time on test [latex]\sum_{\,i\,=\,1}^{n} x_i[/latex] than the corresponding Type II censoring case because a Type I test ends between failures r and [latex]r + 1[/latex]. Thus the expected value of [latex]\hat{\lambda}[/latex] is smaller for Type I censoring than for Type II censoring. One problem that arises with Type I censoring is that the sampling distribution of [latex]\sum_{\,i\,=\,1}^{n} x_i[/latex] is no longer tractable, so an exact confidence interval for λ has not been established. Although many more complicated methods exist, one of the best approximation methods is to assume that [latex]2 \lambda \sum_{\,i\,=\,1}^{n} x_i[/latex] has the chi-square distribution with [latex]2r+1[/latex] degrees of freedom. This approximation, illustrated in Figure 5.13, is based on the fact that if [latex]c = t_{(r)}[/latex], then [latex]2 \lambda \sum_{\,i\,=\,1}^{n} x_i \sim \chi^2 (2r)[/latex], and if [latex]c = t_{{(r}+1)}[/latex], then [latex]2 \lambda \sum_{{\,i\,=\,1}}^{n} x_i \sim \chi^2 (2r+2)[/latex]. Since c is between [latex]t_{(r)}[/latex] and [latex]t_{(r + 1)}[/latex], [latex]2 \lambda \sum_{{\,i\,=\,1}}^{n} x_i[/latex] will be approximately chi-square with [latex]2r+1[/latex] degrees of freedom. This constitutes a proof, after a little algebra, of the following result.

5.5 Weibull Distribution

The Weibull distribution is typically more appropriate for modeling the lifetimes of items with a strictly increasing or decreasing hazard function, such as mechanical items. Rather than looking at each censoring mechanism (for example, no censoring, Type II censoring, Type I censoring) individually, we proceed directly to the general case of random censoring.

Maximum likelihood estimators. As before, let [latex]t_1 , \, t_2 , \, \ldots , \, t_n[/latex]be the failure times, [latex]c_1 , \, c_2 , \, \ldots , \, c_n[/latex] be the associated censoring times, and [latex]{x_i = \min \{t_i, \, c_i \}}[/latex] for [latex]i= 1, \, 2, \, \ldots, \, n[/latex]. The Weibull distribution has hazard and cumulative hazard functions

$$\begin{array}{l}h(t,\lambda ,\kappa )=\kappa \lambda (\lambda t{)}^{\kappa -1}t\ge 0\end{array}$$

and

$$\begin{array}{l}H(t,\lambda ,\kappa )=(\lambda t{)}^{\kappa }t\ge 0.\end{array}$$

When there are r observed failures, the log likelihood function is

$$\begin{array}{ll}\ln L(\lambda ,\kappa )& =\sum _{i\in U}\ln h(x_i,\lambda ,\kappa )-\sum _{i=1}^{n}H(x_i,\lambda ,\kappa )\\ & =\sum _{i\in U}(\ln \kappa +\kappa \ln \lambda +(\kappa -1)\ln x_i)-\sum _{i=1}^n(\lambda x_i{)}^{\kappa }\\ & =r\ln \kappa +\kappa r\ln \lambda +(\kappa -1)\sum _{i\in U}\ln x_i-{\lambda }^{\kappa }\sum _{i=1}^n{x}_i^{\kappa },\end{array}$$

and the [latex]2 \times 1[/latex] score vector has elements

$$\begin{array}{l}{U}_1(\lambda ,\kappa )=\frac{\partial \ln L(\lambda ,\kappa )}{\partial \lambda }=\frac{\kappa r}{\lambda }-\kappa {\lambda }^{\kappa -1}\sum _{i=1}^n{x}_i^{\kappa }\end{array}$$

and

$$\begin{array}{l}{U}_2(\lambda ,\kappa )=\frac{\partial \ln L(\lambda ,\kappa )}{\partial \kappa }=\frac{r}{\kappa }+r\ln \lambda +\sum _{i\in U}\ln x_i-\sum _{i=1}^n(\lambda x_i{)}^{\kappa }\ln (\lambda x_i).\end{array}$$

When these equations are set equal to zero, the simultaneous equations have no closed-form solution for [latex]\hat{\lambda}[/latex] and [latex]\hat \kappa[/latex]:

$$\begin{array}{l}\frac{\kappa r}{\lambda }-\kappa {\lambda }^{\kappa -1}\sum _{i=1}^n{x}_i^{\kappa }=0,\end{array}$$

$$\begin{array}{l}\frac{r}{\kappa }+r\ln \lambda +\sum _{i\in U}\ln x_i-\sum _{i=1}^n(\lambda x_i{)}^{\kappa }\ln (\lambda x_i)=0.\end{array}$$

One piece of good fortune, however, to avoid solving a [latex]2 \times 2[/latex] set of nonlinear equations, is that the first equation can be solved for λ in terms of [latex]\kappa[/latex] as

$$\begin{array}{l}\lambda ={\left(\frac{r}{\sum _{i=1}^n{x}_i^{\kappa }}\right)}^{1/\kappa }.\end{array}$$

Notice that λ reduces to the maximum likelihood estimator for the exponential distribution when [latex]\kappa = 1[/latex]. Using this expression for λ in terms of [latex]\kappa[/latex] in the second element of the score vector yields a single, albeit more complicated, expression with [latex]\kappa[/latex] as the only unknown. After applying some algebra, this equation reduces to

$$\begin{array}{l}g(\kappa )=\frac{r}{\kappa }+\sum _{i\in U}\ln x_i-\frac{r\sum _{i=1}^n{x}_i^{\kappa }\ln x_i}{\sum _{i=1}^n{x}_i^{\kappa }}=0,\end{array}$$

which must be solved iteratively. One technique that can be used to solve this equation is the Newton–Raphson procedure, which uses

$$\begin{array}{l}{\kappa }_{j+1}={\kappa }_j-\frac{g({\kappa }_j)}{g^{\mathrm{\prime <!--\; \prime \; -->}}({\kappa }_j)},\end{array}$$

where [latex]\kappa_0[/latex] is an initial estimator. The iterative procedure can be repeated until the desired accuracy for [latex]\kappa[/latex] is achieved; that is, [latex]| \kappa_{j+1} - \kappa_j | < \epsilon[/latex], for some small positive real number ϵ. When the accuracy is achieved, the maximum likelihood estimator [latex]\hat{ \kappa}[/latex] is used to calculate [latex]\hat{\lambda} = \big(r / \sum_{\,i\,=\,1}^{n} x_i^{\hat{\kappa}}\big)^{1 / {\hat \kappa}}[/latex]. The derivative of [latex]g(\kappa)[/latex] reduces to

$$\begin{array}{l}{g}^{\mathrm{\prime <!--\; \prime \; -->}}(\kappa )=-\frac{r}{{\kappa }^2}-\frac{r}{{\left(\sum _{i=1}^n{x}_i^{\kappa }\right)}^2}[\left(\sum _{i=1}^n{x}_i^{\kappa }\right)(\sum _{i=1}^n(\ln x_i{)}^2{x}_i^{\kappa })-{\left(\sum _{i=1}^n{x}_i^{\kappa }\ln x_i\right)}^2].\end{array}$$

Determining an initial estimator [latex]\kappa_0[/latex] is not trivial. When there are no censored observations, Menon’s initial estimator for [latex]\kappa_0[/latex] is

$$\begin{array}{l}{\kappa }_0={\left\{\frac{6}{(n-1){\pi }^2}[\sum _{i=1}^n(\ln t_i{)}^2-\frac{{\left(\sum _{i=1}^n\ln t_i\right)}^2}{n}]\right\}}^{-1/2}.\end{array}$$

Least squares estimation can be used in the case of a right-censored data set. The Newton–Raphson procedure can fail to converge to the maximum likelihood estimators. A bisection algorithm or fixed point algorithm often provides more reliable convergence.

Fisher and observed information matrices. The [latex]2 \times 2[/latex] Fisher and observed information matrices are based on the following partial derivatives:

$$\begin{array}{ll}\frac{-{\partial }^2\ln L(\lambda ,\kappa )}{\partial {\lambda }^2}& =\frac{\kappa r}{{\lambda }^2}+\kappa (\kappa -1){\lambda }^{\kappa -2}\sum _{i=1}^n{x}_i^{\kappa },\\ \frac{-{\partial }^2\ln L(\lambda ,\kappa )}{\partial \lambda \partial \kappa }& =-\frac{r}{\lambda }+[\left(\kappa {\lambda }^{\kappa -1}\right)\left(\sum _{i=1}^n{x}_i^{\kappa }\ln x_i\right)+\left(\sum _{i=1}^n{x}_i^{\kappa }\right)(\kappa {\lambda }^{\kappa -1}\ln \lambda +{\lambda }^{\kappa -1})]\\ & =-\frac{r}{\lambda }+{\lambda }^{\kappa -1}[\kappa \sum _{i=1}^n{x}_i^{\kappa }\ln x_i+(1+\kappa \ln \lambda )\sum _{i=1}^n{x}_i^{\kappa }],\\ \frac{-{\partial }^2\ln L(\lambda ,\kappa )}{\partial {\kappa }^2}& =\frac{r}{{\kappa }^2}+\sum _{i=1}^n(\lambda x_i{)}^{\kappa }(\ln \lambda x_i{)}^2.\end{array}$$

The expected values of these quantities are not tractable, so the Fisher information matrix does not have closed-form elements. The observed information matrix, however, can be determined by using [latex]\hat{\lambda}[/latex] and [latex]\hat{\kappa}[/latex] as arguments in these expressions.

5.6 Proportional Hazards Model

Parameter estimation for the proportional hazards model, which was introduced in Section 4.6, is considered in this section. Since there is now a vector of covariates in addition to a failure or censoring time for each item on test, special notation must be established to accommodate the covariates. The proportional hazards model has the unique feature that the baseline distribution need not be defined in order to estimate the regression coefficients associated with the covariates.

A lifetime model that incorporates a vector of covariates [latex]{\boldsymbol z} = ( z_1, \, z_2, \, \ldots, \, z_q) ^ \prime[/latex] models the impact of the covariates on survival. The reason for including this vector may be to determine which covariates significantly affect survival, to determine the distribution of the lifetime for a particular setting of the covariates, or to fit a more complicated distribution from a small data set, as opposed to fitting separate distributions for each level of the covariates.

The proportional hazards model was defined in Section 4.6 by

$$\begin{array}{l}h(t,\mathit{z})=\psi (\mathit{z})h_0(t),\end{array}$$

for [latex]t \ge 0[/latex], where [latex]h_0 (t)[/latex] is a baseline hazard function. The covariates increase the hazard function when [latex]\psi ({\boldsymbol z}) > 1[/latex] or decrease the hazard function when [latex]\psi ({\boldsymbol z}) < 1[/latex]. The goal of this section is to develop techniques for estimating the [latex]q \times 1[/latex] vector of regression coefficients [latex]\boldsymbol \beta[/latex] from a data set consisting of n items on test and r observed failure times.

The notation used to describe a data set in a lifetime model involving covariates will borrow some notation established earlier in this chapter, but also establish some new notation. As before, n is the number of items on test and r is the number of observed failures. The failure time of the ith item on test, t_i, is either observed or right censored at time c_i, for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. As before, let [latex]x_i = \min \{t_i, \, c_i\}[/latex] and δ_i be a censoring indicator variable (1 for an observed failure and 0 for a right-censored value), for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. In addition, a [latex]q \times 1[/latex] vector of covariates [latex]{\boldsymbol z}_i = (z_{i1}, \, z_{i2}, \, \ldots, \, z_{iq}) ^ \prime[/latex] is collected for each item on test, for [latex]i = 1, \, 2, \, \ldots, \, n[/latex]. Thus, [latex]z_{ij}[/latex] is the value of covariate j for item i, for [latex]i = 1, \, 2, \, \ldots, \, n[/latex] and [latex]j = 1, \, 2, \, \ldots, \, q[/latex]. This formulation of the problem can be stated in matrix form as

$$\begin{array}{l}\mathit{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]\mathit{\delta }=\left[\begin{array}{c}{\delta }_1\\ {\delta }_2\\ ⋮\\ {\delta }_n\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}\mathit{Z}=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \dots & z_{1q}\\ z_{21}& z_{22}& \dots & z_{2q}\\ ⋮& ⋮& \ddots & ⋮\\ z_{n1}& z_{n2}& \dots & z_{nq}\end{array}\right].\end{array}$$

Each row in the [latex]{\boldsymbol Z}[/latex] matrix consists of the values of the q covariates collected on a particular item. The matrix approach is useful because complicated systems of equations can be expressed compactly and operations on data sets can be performed efficiently by a computer. For parameter estimation, the survivor, density, hazard, and cumulative hazard functions now have the extra arguments [latex]{\boldsymbol z}[/latex] and [latex]\boldsymbol \beta[/latex] associated with them:

$$
S(t, \, {\boldsymbol z} , \, \boldsymbol \theta , \, \boldsymbol \beta) \qquad \qquad
f(t, \, {\boldsymbol z} , \, \boldsymbol \theta , \, \boldsymbol \beta) \qquad \qquad
h(t, \, {\boldsymbol z} , \, \boldsymbol \theta , \, \boldsymbol \beta) \qquad \qquad
H(t, \, {\boldsymbol z} , \, \boldsymbol \theta , \, \boldsymbol \beta) ,
$$

for [latex]t \ge 0[/latex], where the vector [latex]\boldsymbol \theta = (\theta_1, \, \theta_2, \, \ldots, \, \theta_p) ^ \prime[/latex] consists of the p unknown parameters associated with the baseline distribution, which must be estimated along with the regression coefficients [latex]\boldsymbol \beta[/latex].

Parameter estimation for the proportional hazards model can be divided into two cases. The first case is when the baseline distribution is known. This case applies when previous test results have indicated that a particular functional form of the baseline distribution is appropriate. The second case is when the baseline distribution is unknown. This is almost certainly the case when looking at a data set of lifetimes and covariates for the first time without any guidance with respect to an appropriate baseline distribution.

5.7 Exercises

• 5.1 Consider a large batch of light bulbs whose lifetimes are known to have exponential(1) lifetimes. Gina knows that the population distribution is exponential, but she does not know the value of the population mean. She estimates the population mean lifetime of the light bulbs by averaging n observed lifetimes from bulbs chosen at random from the batch. Find the smallest value of n that assures, with probability of at least 0.95, that the sample mean is within 0.2 of the population mean

  1. exactly,
  2. approximately, using the central limit theorem.

• 5.2 Libby is a statistician for a light bulb company. She knows that the lifetimes of the 60-watt bulbs that her company manufactures are exponentially distributed with population mean 1500 hours. She conducts a life test in which 39 of their 60-watt bulbs are placed on life test until they fail and the average of the failure times is recorded. Find the probability that the sample mean exceeds 1600 hours using

  1. the central limit theorem,
  2. the exact distribution of the sample mean.

• 5.3 Let [latex]t_1, \, t_2, \, \ldots, \, t_{n}[/latex] be a random sample from an exponential[latex](\lambda)[/latex] population, where λ is a positive unknown failure rate parameter. Find an unbiasing constant c_n so that [latex]c_n t_{(1)}[/latex] is an unbiased estimator of [latex]1 / \lambda[/latex], where [latex]t_{(1)} = \min \left\{ t_1, \, t_2, \, \ldots, \, t_n \right\}[/latex] is the first order statistic. Hint: the unbiasing constant c_n is a function of the number of items on test n.

• 5.4 Debbie purchases a laptop computer with a random lifetime T whose probability distribution is a special case of the log logistic distribution with survivor function

  $$\begin{array}{l}S(t)=\frac{1}{1+\lambda t}t>0,\end{array}$$

  where λ is a positive unknown scale parameter. From just a single observation of the lifetime of her laptop computer, find an exact two-sided 90% confidence interval for λ.

• 5.5 Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from an exponential population with mean θ, where θ is a positive unknown parameter. An exact two-sided 90% confidence interval for θ is

  $$\begin{array}{l}27<\theta <55.\end{array}$$

  Carol is not concerned about large values of θ. Only small values of θ are of concern. What is an exact one-sided 95% confidence interval of the form [latex]\theta > k[/latex], for some constant k?

• 5.6 If [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] are n mutually independent observations from a log normal distribution with probability density function

  $$\begin{array}{l}f(t)=\frac{1}{\sqrt{2\pi }\sigma t}{e}^{-\frac{1}{2}{\left(\frac{\ln t-\mu }{\sigma }\right)}^2}t\ge 0\end{array}$$

  for [latex]\sigma > 0[/latex] and [latex]-\infty < \mu < \infty[/latex], find the maximum likelihood estimators of μ and σ and exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence intervals for μ and σ in terms of [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex].

• 5.7 Let [latex]t_1, \, t_2, \, \ldots , \, t_7[/latex] be a random sample of the lifetimes of [latex]n = 7[/latex] items on test drawn from an exponential population with positive unknown mean θ.

  1. Find an exact two-sided 90% confidence interval for the median by finding a pivotal quantity based on the sample median [latex]t_{(4)}[/latex].
  2. Give an exact two-sided 90% confidence interval for the median for the [latex]n = 7[/latex] rat survival times in the treatment group from Efron and Tibshirani (1993, page 11):

     $$\begin{array}{l}16,23,38,94,99,141,197.\end{array}$$

  3. Conduct a Monte Carlo simulation experiment to provide convincing numerical evidence that the exact two-sided 90% confidence interval for the median is indeed an exact two-sided 90% confidence interval for an exponential population when θ is arbitrarily set to 1.

• 5.8 This chapter has emphasized confidence intervals. Another type of statistical interval is known as a prediction interval, which contains a future value of an observation with a prescribed probability. Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from an exponential population with a positive unknown mean θ. Conduct a Monte Carlo simulation experiment that provides convincing numerical evidence that the [latex]{100(1 - \alpha)}\%[/latex] prediction interval for [latex]t_{n+1}[/latex]

  $$\begin{array}{l}\frac{\overline{t}}{F_{2n,2,\alpha /2}}<t_{n+1}<\frac{\overline{t}}{F_{2n,2,1-\alpha /2}}\end{array}$$

  is an exact prediction interval for the arbitrary parameter settings [latex]n = 11[/latex], [latex]\alpha = 0.05[/latex], and [latex]\theta = 19[/latex].

• 5.9 Let [latex]T_1, \, T_2, \, T_3[/latex] be mutually independent random variables such that T_i is exponentially distributed with mean [latex]i \kern 0.02em \theta[/latex], for [latex]i = 1, \, 2, \, 3[/latex], where θ is a positive unknown parameter.

  1. Find the maximum likelihood estimator [latex]\hat \theta[/latex].
  2. Find the probability density function of the maximum likelihood estimator [latex]\hat \theta[/latex].
  3. Is [latex]\hat \theta[/latex] an unbiased estimator of θ?
  4. Find an exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for θ.
  5. Perform a Monte Carlo simulation experiment to evaluate the coverage of the confidence interval for [latex]\theta = 10[/latex] and [latex]\alpha = 0.1[/latex].

• 5.10 Let [latex]t_1, \, t_2, \, \ldots , \, t_n[/latex] be a random sample from a population with probability density function

  $$\begin{array}{l}f(t)=\frac{\theta }{t^{\theta +1}}t\ge 1,\end{array}$$

  where θ is a positive unknown parameter.

  1. Find the maximum likelihood estimator of θ.
  2. Use the invariance property to find the maximum likelihood estimator of the median of the distribution.

• 5.11 Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from a population with probability density function

  $$\begin{array}{l}f(t)=\sqrt{\frac{\lambda }{2\pi t^3}}{e}^{-\lambda (t-1{)}^2/(2t)}t>0,\end{array}$$

  where λ is a positive unknown parameter. This distribution is known as the standard Wald distribution which is a special case of the inverse Gaussian distribution Find the maximum likelihood estimator of λ.

• 5.12 Let [latex]T_1, \, T_2, \, \ldots , \, T_n[/latex] be mutually independent and identically distributed random variables from a population having probability density function

  $$\begin{array}{l}f(t)=7e^{-7(t-\theta )}t\ge \theta .\end{array}$$

  Find the limiting distribution of [latex]n \left( T_{(1)} - \theta \right)[/latex]. Support this limiting distribution by conducting a Monte Carlo simulation experiment.

• 5.13 Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from a population with probability density function

  $$\begin{array}{l}f(t)=\frac{\theta }{(1+t{)}^{\theta +1}}t\ge 0,\end{array}$$

  where θ is a positive unknown parameter. Calculate an asymptotically exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for θ based on the asymptotic normality of the maximum likelihood estimator.

• 5.14 Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from a population with probability density function

  $$\begin{array}{l}f(t)=\sqrt{\frac{1}{2\pi t^3}}{e}^{-(t-\theta {)}^2/(2t{\theta }^2)}t>0,\end{array}$$

  where θ is a positive unknown parameter. This population distribution is a special case of the inverse Gaussian distribution. Calculate an asymptotically exact two-sided [latex]{100(1 - \alpha)}\%[/latex] confidence interval for θ based on the asymptotic normality of the maximum likelihood estimator. Hint: the expected value of T is [latex]E[T] = \theta[/latex].

• 5.15 Let [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] be a random sample from a population with probability density function

  $$\begin{array}{l}f(t)=\frac{\theta }{(1+\theta t{)}^2}t\ge 0,\end{array}$$

  where θ is a positive unknown parameter. This is a special case of the log logistic distribution.

  1. Find the maximum likelihood estimator of θ. Hint: The maximum likelihood estimator cannot be expressed in closed form.
  2. Find the maximum likelihood estimate of θ for the [latex]n = 7[/latex] rat survival times (in days) of the treatment group from Efron and Tibshirani (1993, page 11):

     $$\begin{array}{l}16,23,38,94,99,141,197.\end{array}$$

  3. Find an asymptotically exact two-sided 95% confidence interval for θ based on the likelihood ratio statistic for the rat survival times from part (b).

• 5.16 If n items from an exponential population with failure rate λ are placed on a life test that is terminated after r failures have occurred, show that

  $$\begin{array}{l}V\left[\sum _{i=1}^n{x}_i\right]=\frac{r}{{\lambda }^2},\end{array}$$

  where [latex]x_i = \min \{t_i, \, c_i \}[/latex], t_i is the time to failure of the ith item, and c_i is the right-censoring time for the ith item, [latex]i = 1, \, 2, \, \ldots, \, n[/latex].

• 5.17 If n items from an exponential population with failure rate λ are placed on a life test that is terminated after r failures have occurred, find the expected time to complete the test if

  1. failed items are not replaced,
  2. failed items are immediately replaced with new items.

• 5.18 Find the score, maximum likelihood estimator, and Fisher information matrix for a Type II censored random sample from a population with

  $$\begin{array}{l}f(t)=\frac{\theta }{t^{\theta +1}}t\ge 1,\end{array}$$

  where θ is a positive unknown parameter.

• 5.19 The lifetimes of studio light bulbs, measured in days, is exponentially distributed with an unknown failure rate λ. James places n studio light bulbs on test at noon on one day and subsequently checks for failed bulbs at noon on subsequent days until all bulbs have failed. Let [latex]r_1, \, r_2, \, \ldots, \, r_k[/latex] be the number of observed bulb failures, some of which may be zero, on the k days that the bulbs are inspected. Find the maximum likelihood estimator for λ. Also, give the maximum likelihood estimate for the data values [latex]r_1 = 8[/latex], [latex]r_2 = 5[/latex], [latex]r_3 = 2[/latex], [latex]r_5 = 1[/latex], and all other r_i values equal zero.

• 5.20 James’s friend Alexandra decides to simplify matters from the previous question by assuming that all failures that occur during any interval occur at midnight. What is Alexandra’s maximum likelihood estimator for λ as a function of n and [latex]r_1, \, r_2, \, \ldots, \, r_k[/latex]?

• 5.21 Dre conducts a life test on n items from an exponential population with mean θ. He observes only the value of a single order statistic [latex]t_{(k)}[/latex], where k is known. So [latex]k - 1[/latex] lifetimes are left censored at [latex]t_{(k)}[/latex], one lifetime is observed at [latex]t_{(k)}[/latex], and [latex]n - k[/latex] lifetimes are right censored at [latex]t_{(k)}[/latex].

  1. What is the score statistic for estimating θ?
  2. What is the maximum likelihood estimator for θ when [latex]n = 30[/latex], [latex]k = 11[/latex], and [latex]t_{(11)} = 15.5[/latex]?

• 5.22 Consider a Type II right-censored life test with n items on test and [latex]r = 1[/latex] failure is observed at time [latex]t_{(1)}[/latex]. Assume that the items placed on the life test have lifetimes that are well described by a Rayleigh(λ) population.

  1. What is the maximum likelihood estimator for λ?
  2. What is an exact confidence interval for λ?
  3. What is the expected width of the confidence interval from part (b)?
  4. Verify the coverage and expected width of the exact confidence interval for [latex]\lambda = 2[/latex], [latex]n = 7[/latex], and [latex]\alpha = 0.05[/latex] via Monte Carlo simulation.

• 5.23 A randomly right-censored data set is collected from a population with hazard function

  $$\begin{array}{l}h(t)=\theta (1+t)t\ge 0,\end{array}$$

  where θ is a positive parameter.

  1. Find the maximum likelihood estimator [latex]\hat \theta[/latex].
  2. Give an expression for the observed information matrix.
  3. Give an asymptotically exact confidence interval for θ based on the observed information matrix.

• 5.24 Candice conducts a life test in which n items are simultaneously placed on test at time 0. The test is concluded at time [latex]c > 0[/latex]. Assuming that the lifetimes of the items are from an exponential population with mean θ, find the distribution of the number of failures that occur by time c.

• 5.25 Show that when a random sample is drawn an exponential(λ) population with Type II right censoring

  $$\begin{array}{l}\frac{2r\lambda }{\hat{\lambda }}\sim {\chi }^2(2r),\end{array}$$

  where [latex]\chi^2 (2r)[/latex] is the chi-square distribution with 2r degrees of freedom.

• 5.26 Consider a Type II right censored sample of n items on test and r observed failures drawn from an exponential population with mean θ. Show that the maximum likelihood estimate [latex]\hat{\theta}[/latex] is unbiased.

• 5.27 Assume that a life test without replacement is conducted on n items from an exponential population with failure rate λ. The exact failure times are not known, but the test is terminated upon the rth ordered failure at time [latex]t_{(r)}[/latex]. Find a point estimator for λ.

• 5.28 Consider a population of items with exponential(λ) lifetimes. A life test with replacement is terminated when r failures occur or at time c, whichever occurs first. This is a combination of Type I and Type II right censoring. Find the expected number of items that fail during the test as a function of λ.

• 5.29 For a life test of n items with exponential(λ) lifetimes (items are not replaced upon failure) which is continued until all items fail, show that

  $$\begin{array}{l}E[\hat{\lambda }]=\frac{n}{n-1}\lambda ,\end{array}$$

  where λ is the population failure rate and [latex]\hat{\lambda}[/latex] is the maximum likelihood estimator for λ. Thus, an unbiasing constant for [latex]\hat{\lambda}[/latex] is [latex]u_n = {(n - 1)} / n[/latex]. Equivalently,

  $$\begin{array}{l}E\left[\frac{n-1}{n}\hat{\lambda }\right]=E\left[u_n\hat{\lambda }\right]=\lambda .\end{array}$$

  Find an unbiasing constant for the case of Type II right censoring.

• 5.30 Give a point and 95% interval estimator for the median lifetime of the 6–MP treatment group assuming that the data have been drawn from an exponential(λ) population.

• 5.31 Consider the following Type II right censored data set for the lifetime of a product ([latex]n = 5[/latex] and [latex]r = 3[/latex]) drawn from an exponential population with failure rate λ:

  $$\begin{array}{l}3.63.98.5.\end{array}$$

  1. Find the maximum likelihood estimator for the mean of the population.
  2. Find the maximum likelihood estimator for [latex]S(5)[/latex].
  3. Find an exact two-sided 80% confidence interval for [latex]E\left[T ^ 3 \right][/latex].
  4. Find an exact one-sided 95% lower confidence interval for [latex]S(5)[/latex].
  5. Find the p-value for the test [latex]H_0: \lambda = 0.04[/latex] versus [latex]H_1: \lambda > 0.04[/latex].
  6. Find the value of the log likelihood function at the maximum likelihood estimate.
  7. Find the value of the observed information matrix.
  8. Assume the data values

     $$\begin{array}{l}3.84.66.09.6\end{array}$$

     constitute a complete data set for a different product. Find an exact two-sided 90% confidence interval for the ratio of the failure rates of the two products if both are assumed to come from exponential populations.

• 5.32 Sara observes a single observed lifetime T from an exponential(λ) population, where λ is a positive unknown rate parameter. Find an exact two-sided 95% confidence interval for λ.

• 5.33 Justin places a single item is placed on test ([latex]n = 1[/latex]). The only information that is available is that the item failed between times a and b, where [latex]a < b[/latex]. In other words, the single item’s lifetime is interval censored. Assuming that the population time to failure is exponential(λ), what is the maximum likelihood estimator of λ?

• 5.34 Natalie conducts a life test with [latex]n = 19[/latex] items on test and random right censoring. Let [latex]t_1, \, t_2, \, \ldots , \, t_{19}[/latex] be the independent exponential(2) times to failure. Let [latex]c_1, \, c_2, \, \ldots , \, c_{19}[/latex] be the independent exponential(1) censoring times, which are independent of the times to failure. Use Monte Carlo simulation to estimate the actual coverage of the following approximate confidence interval procedures for the population failure rate λ at for [latex]\alpha = 0.05[/latex].

  1. The confidence interval consisting of all λ satisfying

     $$\begin{array}{l}\frac{\hat{\lambda }{\chi }_{2r,1-\alpha /2}^2}{2r}<\lambda <\frac{\hat{\lambda }{\chi }_{2r,\alpha /2}^2}{2r}.\end{array}$$

  2. The confidence interval consisting of all λ satisfying

     $$\begin{array}{l}\hat{\lambda }-z_{\alpha /2}O(\hat{\lambda }{)}^{-1/2}<\lambda <\hat{\lambda }+z_{\alpha /2}O(\hat{\lambda }{)}^{-1/2}.\end{array}$$

  3. The confidence interval consisting of all λ satisfying

     $$\begin{array}{l}2[\ln L(\hat{\lambda })-\ln L(\lambda )]<{\chi }_{1,\alpha }^2.\end{array}$$

  Replicate the experiment so as to estimate the actual coverages to three digits of accuracy.

• 5.35 Sixty-watt light bulb lifetimes are known to be exponentially distributed with unknown positive population mean θ from previous test results. The company that produces these light bulbs would like to estimate θ by testing n bulbs to failure at one facility and m bulbs to failure at a second facility. Let [latex]X_1, \, X_2, \, \ldots , \, X_n[/latex] be the independent lifetimes of the bulbs tested at the first facility; let [latex]Y_1, \, Y_2, \, \ldots , \, Y_m[/latex] be the independent lifetimes of the bulbs tested at the second facility. An unbiased estimate of θ is the convex combination

  $$\begin{array}{l}\hat{\theta }=p{\hat{\theta }}_X+(1-p){\hat{\theta }}_Y,\end{array}$$

  where [latex]0 < p < 1[/latex], [latex]\hat \theta_X = \bar X[/latex] is the maximum likelihood estimator of θ for the data from the first facility, and [latex]\hat \theta_Y = \bar Y[/latex] is the maximum likelihood estimator of θ for the data from the second facility. Find the value of p that minimizes [latex]V\big[ \kern 0.01em \hat \theta \kern 0.04em \big][/latex].

• 5.36 Ash would like to test the hypothesis

  $$\begin{array}{l}{H}_0:\lambda =17\end{array}$$

  versus

  $$\begin{array}{l}{H}_1:\lambda >17\end{array}$$

  using a single value T from an exponential(λ) population, where λ is a positive unknown population failure rate. The null hypothesis is rejected if [latex]T < 0.01[/latex]. Find the significance level α for the test.

• 5.37 Let T be an observation from an exponential population with positive unknown population mean θ. This observation is used to test

  $$\begin{array}{l}{H}_0:\theta =6\end{array}$$

  versus

  $$\begin{array}{l}{H}_1:\theta =2.\end{array}$$

  1. Find the critical value for the test for a fixed significance level α.
  2. Find β for a fixed significance level α.

• 5.38 Paul collects a random sample [latex]t_1, \, t_2, \, \ldots , \, t_n[/latex] from an exponential population with positive unknown mean θ. Show that the sample mean, [latex]\bar{t}[/latex], and n times the first order statistic, [latex]n t_{(1)}[/latex], are both unbiased estimators of θ.

• 5.39 Jessica and Mary collect a random sample [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] of light bulb lifetimes drawn from an exponential(λ) population, where λ is a positive unknown failure rate. The bulbs are stamped with “1000 hour MTTF,” indicating that the mean time to failure equals 1000 hours. They would like to determine whether there is statistical evidence in the sample that indicates the bulbs last longer than 1000 hours.

  1. State the appropriate H₀ and H₁.
  2. Jessica uses the test statistic [latex]\bar{t}[/latex] and Mary uses the test statistic [latex]nt_{(1)}[/latex] to test the hypothesis. Find the critical values for their test statistics when [latex]\alpha = 0.05[/latex] and [latex]n = 10[/latex].
  3. Draw the power curves associated with each of the test statistics from part (b) on the same set of axes using a computer. Again assume that [latex]\alpha = 0.05[/latex] and [latex]n = 10[/latex].

• 5.40 Camille observes a single lifetime T from an exponential population with a positive unknown population mean θ. She would like to test

  $$\begin{array}{l}{H}_0:\theta =1\end{array}$$

  versus

  $$\begin{array}{l}{H}_1:\theta >1\end{array}$$

  at [latex]\alpha = 0.07[/latex] using T as a test statistic.

  1. Find the critical value c for this test.
  2. Plot the power function for this test.

• 5.41 Ellen collects a random sample [latex]t_1, \, t_2, \, \ldots, \, t_{10}[/latex] of light bulb lifetimes from an exponential(λ) population, where λ is a positive unknown failure rate. Ellen is a reliability engineer. She is confident from previous test results that the time to failure for these light bulbs is exponentially distributed. She is interested in testing whether a manufacturer’s claim that the population mean time to failure for the bulbs is 1000 hours. So she would like to test

  $$\begin{array}{l}{H}_0:\lambda =0.001\end{array}$$

  versus

  $$\begin{array}{l}{H}_1:\lambda >0.001.\end{array}$$

  She is in a hurry. She places ten bulbs on test and only observes the first bulb fail at [latex]t_{(1)} = 14[/latex] hours, and would like to draw a conclusion at 14 hours. Give the p-value for the test based on the value of this single order statistic.

• 5.42 Liz collects a random sample of lifetimes [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] from an exponential(λ) distribution, where λ is a positive unknown failure rate parameter. She conducts a significance test of

  $$\begin{array}{l}{H}_0:\lambda =1\end{array}$$

  versus

  $$\begin{array}{l}{H}_0:\lambda \ne 1,\end{array}$$

  which achieves a p-value of [latex]p = 0.07[/latex] for a particular data set. If she then computes an exact two-sided 95% confidence interval for λ for this particular data set, will the confidence interval contain 1?

• 5.43 Karen fits the ball bearing data set to the Weibull distribution parameterized as

  $$\begin{array}{l}S(t)=e^{-(\lambda t{)}^{\kappa }}t\ge 0,\end{array}$$

  yielding maximum likelihood estimates [latex]\hat{\lambda} = 0.0122[/latex] and [latex]\hat{\kappa} = 2.10[/latex]. Ute also wants to fit the same data set to the Weibull distribution, but she uses the parameterization

  $$\begin{array}{l}S(t)=e^{-\rho t^{\beta }}t\ge 0.\end{array}$$

  What will be the maximum likelihood estimates [latex]\hat{\rho}[/latex] and [latex]\hat{\beta}[/latex] that Ute obtains for the ball bearing data set?

• 5.44 Jay conducts a life test with [latex]n = 5[/latex] items on test which is terminated when [latex]r = 3[/latex] items have failed. Failed items are not replaced in this traditional Type II right-censored data set. Assuming that the time to failure of an item in the population has a Weibull([latex]\lambda, \, \kappa[/latex]) distribution with known, positive parameters λ and [latex]\kappa[/latex], what is the probability density function of the time to complete the life test?

• 5.45 Jennie collects a random sample [latex]t_1, \, t_2, \, \ldots, \, t_7[/latex] from a Rayleigh population with probability density function

  $$\begin{array}{l}f(t)=2{\theta }^{-2}t{e}^{-(t/\theta {)}^2}t>0,\end{array}$$

  where θ is a positive unknown parameter. She would like to test

  $$\begin{array}{l}{H}_0:\theta =10\end{array}$$

  versus

  $$\begin{array}{l}{H}_1:\theta >10\end{array}$$

  using the test statistic [latex]t_{(1)} = \min \left\{ t_1, \, t_2, \, \ldots, \, t_7 \right\}[/latex], which assumes the value [latex]t_{(1)} = 6[/latex]. Find the p-value for her test.

• 5.46 Mildred collects a random sample [latex]t_1, \, t_2, \, \ldots, \, t_n[/latex] from a Rayleigh[latex](\lambda)[/latex] population with survivor function

  $$\begin{array}{l}S(t)=e^{-(\lambda t{)}^2}t>0,\end{array}$$

  where λ is a positive unknown parameter.

  1. Find the maximum likelihood estimator of λ.
  2. Show that the log likelihood function is maximized at the maximum likelihood estimator [latex]\hat{\lambda}[/latex].
  3. Given that the expected value of T is [latex]E[T] = {\sqrt{\pi}} / (2 \lambda)[/latex], find the method of moments estimator of λ.

• 5.47 Find the elements of the score vector for the log logistic distribution for a randomly right-censored data set.

• 5.48 Bryan places n items on test and observes r failures. Assuming that the failure times of the items follow the log logistic distribution and censoring is random, set up an expression for the boundary of a 95% confidence region for the shape parameter [latex]\kappa[/latex] and scale parameter λ of the log logistic distribution based on the likelihood ratio statistic. Assume that the survivor function for the log logistic distribution is

  $$\begin{array}{l}S(t)=\frac{1}{1+(\lambda t{)}^{\kappa }}t\ge 0,\end{array}$$

  for [latex]\lambda > 0[/latex] and [latex]\kappa > 0[/latex]. It is not necessary to solve for the maximum likelihood estimators.

• 5.49 Consider a proportional hazards model with [latex]n = 3[/latex] items on test and distinct failure times [latex]t_1, \, t_2, \, t_3[/latex]. Compute the joint probability mass function values for the [latex]3! = 6[/latex] possible rank vectors, and show that they sum to 1.

• 5.50 Give the equations that must be solved in order to find the maximum likelihood estimators [latex]\hat{\lambda}[/latex], [latex]\hat{\kappa}[/latex], and [latex]\hat {\boldsymbol beta}[/latex] for a proportional hazards model with log logistic baseline distribution and log linear link function. A random right-censoring scheme is used.

• 5.51 Joyce fits the Cox proportional hazards model with unknown baseline distribution given in Examples 5.14, 5.15, and 5.16 to the [latex]n = 3[/latex] light bulb failure times. The purpose of the study was to determine the effect of wattage on survival for 60-watt and 100-watt light bulbs.

  1. What is the value of the regression coefficient for wattage if it were coded as [latex]z = 60[/latex] and [latex]z = 100[/latex] rather than as a binary covariate?
  2. Write a short paragraph indicating whether or not these two approaches are fundamentally equivalent ways of coding the covariate. If they differ, is one method of coding the covariate superior to the other for the purpose of the study?

• 5.52 Survival times (in weeks) for two groups of leukemia patients (AG positive and AG negative blood types), along with an additional covariate, white blood cell count are given in Feigl, P. and Zelen, M., “Estimation of Exponential Survival Probabilities with Concomitant Information,” Biometrics, Vol. 21, No. 4, pp. 826–838, 1965, and are displayed below.

  AG positive group		AG negative group
  Survival time	White blood count	Survival time	White blood count
  65	2300	56	4400
  156	750	65	3000
  100	4300	17	4000
  134	2600	7	1500
  16	6000	16	9000
  108	10500	22	5300
  121	10000	3	10000
  4	17000	4	19000
  39	5400	2	27000
  143	7000	3	28000
  56	9400	8	31000
  26	32000	4	26000
  22	35000	3	21000
  1	100000	30	79000
  1	100000	4	100000
  5	52000	43	100000
  65	100000

  1. Fit the Cox proportional hazards model to the survival times. Code the blood type as the indicator variable z₁, using 1 for AG positive and 0 for AG negative. The second covariate z₂ is the natural logarithm of the white blood cell counts minus the sample mean of the natural logarithms of the white blood cell counts. Include the interaction term [latex](z_1 - \bar z_1) z_2[/latex] in the model. Use the Breslow method for handling tied survival times.
  2. Write a sentence interpreting the sign of [latex]\hat \beta_1[/latex], [latex]\hat \beta_2[/latex], and [latex]\hat \beta_3[/latex] in terms of risk to the patient.
  3. Give a 95% confidence interval for β₁.
  4. If covariates associated with p-values that are less than 0.10 are considered statistically significant, what is the fitted hazard function for a leukemia patient with baseline hazard function [latex]h_0(t)[/latex], white blood cell count 9000 who has AG positive blood type? Hint: The sample mean of the natural logarithms of the white blood cell types is 9.52 and the mean of the blood types coded as an indicator variable is [latex]17 / 33 = 0.515[/latex].

• 5.53 Consider the Cox proportional hazards model with a single ([latex]q = 1[/latex]) binary covariate z₁, an exponential(λ) baseline distribution, and a log linear link function. The baseline distribution can be absorbed into the link function by creating an artificial covariate [latex]z_0 = 1[/latex] and setting [latex]\lambda = e ^ {\beta_0 z_0}[/latex].

  1. For a randomly right-censored data set, find the score vector.
  2. For a randomly right-censored data set, find closed-form expressions for the maximum likelihood estimators [latex]\hat \beta_0[/latex] and [latex]\hat \beta_1[/latex].
  3. For the [latex]n = 3[/latex] observations given in vector form below, calculate the maximum likelihood estimates [latex]\hat \beta_0[/latex] and [latex]\hat \beta_1[/latex].

     $$\begin{array}{l}\mathit{x}=\left[\begin{array}{c}80\\ 20\\ 50\end{array}\right]\mathit{\delta }=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]\mathit{Z}=\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right].\end{array}$$

  4. What is the hazard function of the fitted model for the data from part (c)?
  5. Use the observed information matrix to give approximate two-sided 95% confidence intervals for β₀ and β₁ for the data from part (c).
  6. Give the p-values for testing the hypotheses

     $$\begin{array}{ll}{H}_0:& {\beta }_i=0\\ H_1:& {\beta }_i\ne 0\end{array}$$

     for [latex]i = 0, \, 1[/latex], for the data from part (c).

• 5.54 The wattage of the [latex]n = 3[/latex] light bulbs in Example 5.16 was coded as the covariate [latex]z_1 = 0[/latex] for a 60-watt bulb and [latex]z_1 = 1[/latex] for a 100-watt bulb. When the Cox proportional hazards model with an unspecified baseline hazard function was fit to the data set, the point estimate for the regression parameter was [latex]\hat \beta_1 = -0.347[/latex]. Without doing the derivation from scratch, what is the point estimate for the regression parameter if the wattage (that is, 60 watts or 100 watts) of the bulb were used as the covariate.

• 5.55 Mark fits a Cox proportional hazards model with unknown baseline distribution to a data set of drill bit failure times (measured in number of items drilled) with [latex]q = 2[/latex], for which the covariates denote the turning speed (revolutions per minute, rpm) and the hardness of the material (Brinell hardness number, BHN) being drilled. The turning speeds range from 2400 to 4800 rpm and the hardness of the materials ranges from 250 to 440 BHN. Interactions are not considered and the variables are not centered. The fitted model has estimated regression vector [latex]\hat {\boldsymbol \beta} = (0.014, \, 0.45) ^ \prime[/latex], and the inverse of the observed information matrix is

  $$\begin{array}{l}{O}^{-1}(\hat{\mathbf{\beta }})=\left[\begin{array}{cc}0.000081& 0.000016\\ 0.000016& 0.010000\end{array}\right].\end{array}$$

  Write a paragraph interpreting these results.