Statistical Hypothesis Testing of Failure-Time Data in Time-to-Event (or Survival) Analysis

Introduction

The following article is my attempt at applying the concepts I learnt in the chapter on non-parametric hypothesis testing¹ to failure time data sets in reliability.

The log rank test discussed in this article deals with testing for difference in failure rates between different groups. The intent is to check whether there is a practical and statistical difference between the failure rates of the different groups. One would want to make such comparisons for various reasons. For example, suppose we have failure data of two components from different vendors and we want to find out which component has a lower failure rate. Another reason could be to check for any differences in failure rates before and after a process improvement change has been implemented, or between components operating on different mine sites or different conditions. These methods would be applicable to all such datasets that are right censored (or left truncated and right censored).

Two datasets² will be analysed. Each dataset consists of data from at least two groups. All datasets are available on the student companion website of the aforementioned book. They are also publicly available on DataShare, Iowa State University’s open data repository which can be accessed through this link.

Non-parametric methods do not make any assumptions for the distribution of failure times. These methods are normally tried first before moving on to parametric models, which, if the assumptions are met, could lead to more precise conclusions. This blog article discusses the application of the log rank test. Other non-parametric tests like test for trend, stratified tests etc will be discussed in subsequent blogs.

The survival package in R has been used for a lot of the analysis in this article. Additionally, I have also used a few other packages like survminer which helps in plotting survival curves in a much more informative and visually appealing format. I noticed some issues which lead to incorrect results (p-value) in the survminer package, which I discuss in more detail below. It is best to verify your results, whenever possible, either by calculating some results manually or by cross verifying with other methods. For example, the p-value of a log rank test should be more or less the same as the p-value of the score test from a cox proportional hazards model. I have performed this check for one of the datasets below and was able to conclude that the p-value I was getting from the survminer package wasn’t what I was expecting.

The R-markdown file for this article contains all the code used here and is available on this link.

Dataset I: Snubber

This dataset, first given by Nelson (1981)³ contains information on accelerated life tests of two snubber designs for a pop-up toaster. The new design reduces the manufacturing cost, but does it affect the reliability of the component in comparison to the old design? A logrank test can be performed answer this question. We want to test whether the failure rates are different for the two designs. We set up the null and alternate hypothesis as follows.

\[ \begin{aligned} &H_0: h_{old}(t) = h_{new}(t) \quad \text{for all } t \leq \tau, \quad \text{versus} \\ \\ &H_A: h_{old}(t) \, \text{is different to} \, h_{new}(t) \, \text{for some} \, t \leq \tau. \quad \end{aligned} \] Here, \(h(t)\) is the failure rate and \(\tau\) is the largest time at which both groups have at least one observation at risk.

The hypothesis test is conducted by calculating the weighted difference between the actual failures and expected failures under the null hypothesis of no difference in failure rates⁴. This weighted difference, along with its variance-covariance matrix is used to calculate a chi-squared statistic from which the p-value is calculated. The weights allow us to test for early or late departures between the hazard rates. For the log rank test, the value of the weight is 1, which means equal weightage over the whole range⁵.

In this article, we only use the log rank weight of 1.

Next, we load the snubber dataset, rename columns so that they don’t include spaces and create a new column for status.

snubber <- read_csv("Data/Snubber.csv")

snubber <- snubber %>%
  dplyr::mutate(Status = case_when(`Censoring Indicator` == "Failure" ~ 1, `Censoring Indicator` == "Censored" ~ 0)) %>% dplyr::rename(Cycles = `Toaster Cycles`)

knitr::kable(snubber, caption = "Table 1. Snubber", align = rep('c', 5), table.envir = 'table*') %>%
 kableExtra::kable_styling("striped", full_width = T, position = "left") %>%
kableExtra::scroll_box(height = "200px")

The first column indicates the number of toaster cycles to failure or the number of toaster cycles when right censored. The third column is the frequency counts. The fourth column indicates whether the data point is from the new or old design.

In the next step, a Kaplan-Meir survival curve is fit to the two groups.

fit_snubber <- survival::survfit(
  survival::Surv(time = Cycles, event = Status, type = 'right') ~ Design,
  data = snubber,
  weights = Count,
  conf.type = "log",
  conf.int = 0.95
)

ggsurvplot(
  fit_snubber,
  data = snubber,
  risk.table = T,
  conf.int = T
)

Looking at the curves in the plot, there doesn’t seem to be much difference between the two groups. A hypothesis test for difference in hazard rates will formalize our comparison of the two groups.

The survdiff function from the survival package can be used to conduct the log-rank tests. The syntax is as follows:

survival::survdiff(
  survival::Surv(time = Cycles, event = Status, type = 'right') ~ Design,
  data = snubber,
  rho = 0
) ##rho = 0 gives log rank test

## Call:
## survival::survdiff(formula = survival::Surv(time = Cycles, event = Status, 
##     type = "right") ~ Design, data = snubber, rho = 0)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## Design=New 26       13     13.9    0.0547     0.124
## Design=Old 25       13     12.1    0.0626     0.124
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7

The p-value is 0.7, which means we do not have sufficient evidence at the 5% significance level to suggest that there is a difference in the hazard rates. But, there is a problem in the result of survdiff shown above! If you notice the syntax for survdiff, it does not have an argument for “Counts”, which is the frequency of the number of events (or number of censored observations). survdiff treats every row as one count, whereas we want frequency of observations to be factored in just like it was done in the survfit function using the “weights = Count” argument. This can also be verified from the numbers under the “Observed” column in the results of survdiff shown above. The column shows that there are 13 failures in the old design as well as the new design whereas the number of failures in the old design is 18 and in the new design is 16, as shown below.

(sum(snubber$Count[which(snubber$Status ==1 & snubber$Design == "Old")]))

## [1] 18

(sum(snubber$Count[which(snubber$Status ==1 & snubber$Design == "New")]))

## [1] 16

Hence, if we want to use survdiff, we need to modify our dataset in such a way that every row counts as only one observation. We do that in the next chunk of code.

count_gt1.index <- which(snubber$Count > 1) #identify rows with count greater than 1.
count_adj <- snubber$Count[count_gt1.index] - 1

snubber.split <- snubber[rep(row.names(snubber)[count_gt1.index],count_adj),]

snubber.new <- rbind(snubber,snubber.split) %>% dplyr::select(!c("Count")) %>% dplyr::arrange(desc(Design))

sum(snubber$Count) == dim(snubber.new)[1] #check whether sum of counts equals number of rows in the new dataset

## [1] TRUE

knitr::kable(snubber.new, caption = "Table 2. Snubber (Modified)", align = rep('c', 4), table.envir = 'table*') %>%
 kableExtra::kable_styling("striped", full_width = T, position = "left") %>%
kableExtra::scroll_box(height = "200px")

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

The modified table now has 106 rows, with one row for each observation. We can now apply the survdiff function from the survival package.

survival::survdiff(
  survival::Surv(time = Cycles, event = Status, type = 'right') ~ Design,
  data = snubber.new,
  rho = 0
) ##rho = 0 gives log rank test

## Call:
## survival::survdiff(formula = survival::Surv(time = Cycles, event = Status, 
##     type = "right") ~ Design, data = snubber.new, rho = 0)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## Design=New 54       16     17.1    0.0697     0.152
## Design=Old 52       18     16.9    0.0704     0.152
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.7

The “Observed” column above shows the right number of failures. The p-value is still 0.7, which again indicates that we do not have sufficient evidence at the 5% significance level to support the alternate hypothesis of difference in the failure rates of the new and old design.

I also calculated the log rank test statistic and the p-value manually and was able to get a similar result⁶.

We now plot the Kaplan Meir curves along with the p-value on the plot.

survminer::ggsurvplot(
  fit_snubber,
  pval = 0.7,
  risk.table = T,
  conf.int = T
)

There is another note of caution that is warranted here. The ggsurvplot function from the survminer package which we have used to generate the plot above, also has the ability to conduct tests of 2 or more groups with different weights. The calculated p-value is plotted on the survival curve through additional arguments as shown below.

survminer::ggsurvplot(
  fit_snubber,
  pval = T, #default is weight = 1 for log rank test
  pval.method = T,
  risk.table = T,
  conf.int = T
)

Since the Kaplan-Meir curve model “fit_snubber” is an input to ggsurvplot, one could easily mistake the p-value reported on the plot to be the correct one that takes into account the frequency counts from the “fit_snubber” model object. Unfortunately, this is not the case. ggsurvplot calls the survdiff function by default on the unmodified dataframe in the “fit_snubber” model object, which leads to the incorrect p-value. Hence, if you are planning to use ggsurvplot to plot Kaplan Meir curves, it is recommended that you hardcode the p-value like we did initially after the correct calculations through survdiff, instead of relying on ggsurvplot calculations.

Next, we will look at an example where there is a big change in the p-value before and after we modify the dataframe to make it suitable for applying the survdiff function.

Dataset II: Circuit Packs

This dataset contains life test data from a test conducted to compare failure time distribution of components from two different vendors. It was first reported by Hooper and Amster (1998)⁷

Just like the previous example, we load the dataset and clean up column names.

circuitPack <- read_csv("Data/CircuitPack05.csv")

## Rows: 41 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Censoring Indicator, Vendor Number
## dbl (2): Days, Count
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

circuitPack <- circuitPack %>%
  dplyr::mutate(Status = case_when(`Censoring Indicator` == "Failure" ~ 1, `Censoring Indicator` == "Censored" ~ 0)) %>% 
  dplyr::rename(Vendor = `Vendor Number`)

knitr::kable(circuitPack, caption = "Table 3. Circuit Pack", align = rep('c', 5), table.envir = 'table*') %>%
 kableExtra::kable_styling("striped", full_width = T, position = "left") %>%
kableExtra::scroll_box(height = "200px")

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

Notice that the number of censored components in the “Count” column at 365 days are much higher than the number of failures.

We set up the null and alternate hypothesis as follows.

\[ \begin{aligned} &H_0: h_{Vendor 1}(t) = h_{Vendor 2}(t) \quad \text{for all } t \leq \tau, \quad \text{versus} \\ \\ &H_A: h_{Vendor 1}(t) \, \text{is different to} \, h_{Vendor 2}(t) \, \text{for some} \, t \leq \tau. \end{aligned} \]

Here is the Kaplan Meir curve for this scenario.

fit <- survival::survfit(
  survival::Surv(time = Days, event = Status, type = 'right') ~ Vendor,
  data = circuitPack,
  weights = Count,
  conf.type = "log",
  conf.int = 0.95
)

ggsurvplot(fit, data = circuitPack, risk.table = T, ylim = c(0.95,1),  conf.int = T)

Again, the survival probability plot doesn’t suggest much difference between the two curves.

The log rank test using survdiff directly on the dataset is shown below.

survival::survdiff(
  survival::Surv(time = Days, event = Status, type = 'right') ~ Vendor,
  data = circuitPack,
  rho = 0
) ##rho = 0 gives log rank test

## Call:
## survival::survdiff(formula = survival::Surv(time = Days, event = Status, 
##     type = "right") ~ Vendor, data = circuitPack, rho = 0)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## Vendor=Vendor1 22       21     20.2    0.0347    0.0724
## Vendor=Vendor2 19       18     18.8    0.0372    0.0724
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.8

The column “N” in the output of survdiff shows 22 for Vendor 1 and 19 for Vendor 2. But the total number of observations are much higher than these values, as shown in the next chunk of code. Also, notice that the p-value is 0.8.

(sum(circuitPack$Count[which(circuitPack$Vendor == "Vendor1")]))

## [1] 1041

(sum(circuitPack$Count[which(circuitPack$Vendor == "Vendor2")]))

## [1] 1245

There are 1041 observations from Vendor 1 and 1245 observations from Vendor. This discrepancy will lead to a big change in the p-value once we modify the data frame like we did in the first example.

count_gt1.index <- which(circuitPack$Count > 1) #identify rows with count greater than 1.
count_adj <- circuitPack$Count[count_gt1.index] - 1

circuitPack.split <- circuitPack[rep(row.names(circuitPack)[count_gt1.index],count_adj),]

circuitPack.new <- rbind(circuitPack,circuitPack.split) %>% dplyr::select(!c("Count")) %>% dplyr::arrange((Vendor))

sum(circuitPack$Count) == dim(circuitPack.new)[1] #check whether sum of counts equals number of rows in the new dataset

## [1] TRUE

knitr::kable(circuitPack.new, caption = "Table 4. Circuit Pack (Modified)", align = rep('c', 4), table.envir = 'table*') %>%
 kableExtra::kable_styling("striped", full_width = T, position = "left") %>%
kableExtra::scroll_box(height = "200px")

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")