R code for constructing likelihood based confidence intervals for the life time quantiles of a competing risks model

Competing risks model are employed when, for instance, a device has two different causes of failure (also referred to as failure modes). The R code below shows how to model the failure data of such a device with a competing risks model.
Furthermore, the R code demonstrates how to compute likelihood based confidence intervals for the life time quantiles of the competing risks model. In computing these confidence intervals the R code assumes that the observed failure times follow a Weibull distribution. However, it will also be demonstrated how to adapt the code in case the failure times follow a lognormal distribution.
It should be noted that the competing risks model in the R code focuses on a situation in which a device has two failure modes. Nevertheless, the code can easily be extended to include more than two failure modes.

competing risks model ci

An introduction to competing risks models is given by Meeker and Escobar in their 1998 book Statistical methods for reliability data (pp. 382-385).

A procedure for computing likelihood based intervals for the life time quantiles of a competing risks model can be found this paper by Hong and Meeker.

Finally, see this blog post for computing likelihood based intervals for the failure probabilities of a competing risks model.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s R code blog posts.

R code

library(survival)
library(SPREDA)
library(pracma)

##data

#Device G data (see: Meeker & Escobar, Example 15.6, pp. 383-385)
#note that Device G fails due to either:
#1) electrical surge (=failure mode S)
#2) product wear (=failure mode W)
devG <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/DeviceG.txt", header=T)
devG$Cens <- as.numeric(devG$Cause)-1
head(devG, n=10)

#specify censoring for competing risks model

#failures due to electrical surge (=failure mode S)
#1=failure due to mode S; 0=censored
devG$censS <- ifelse(devG$Cens==0, 1, 0)

#failure due to product wear (=failure mode W)
#1=failures due to mode W; 0=censored
devG$censW <- ifelse(devG$Cens==2, 1, 0)



##fit failure-time models to data
#assuming that the failure times for both failure modes follow a Weibull distribution

#analysis of failure mode S only
modS <- survreg(Surv(Kilocycles, censS) ~ 1, data=devG, dist="weibull")
summary(modS)

#analysis of failure mode W only
modW <- survreg(Surv(Kilocycles, censW) ~ 1, data=devG, dist="weibull")
summary(modW)

#extract MLEs for model parameters
mu1 <- coef(modS)
sigma1 <- modS$scale
mu2 <- coef(modW)
sigma2 <- modW$scale



##likelihood based confidence intervals for life time quantiles tp
#tp is the time (t) at which a specified proportion (p) of the units fails

#specify proportion
p <- .01

#function for computing abs(p-f(mu1,sigma1,mu2,sigma2)),
#where p is the specified proportion, and f(.) some function
#of mu1, sigma1, mu2, and sigma2
mu1Min <- function (mu1Par, sigma1, mu2, sigma2, tpHat){
  z1 <- (log(tpHat)-mu1Par)/sigma1
  z2 <- (log(tpHat)-mu2)/sigma2
  abs(p - (1- (1-psev(z1)) * (1-psev(z2)) ) )}

#log-likelihood function for competing risks model
llikweibullcr <- function (theta, y, d1, d2) {
  mu1 <- theta[1]
  sigma1 <- exp(theta[2])
  mu2 <- theta[3]
  sigma2 <- exp(theta[4])
  
  z1 <- (log(y)-mu1)/sigma1
  z2 <- (log(y)-mu2)/sigma2
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log( ((( 1/sigma1*y) * dsev(z1))^d1)*( (1-psev(z1)) ^(1-d1) ) )) +
    sum(log( ((( 1/sigma2*y) * dsev(z2))^d2)*( (1-psev(z2)) ^(1-d2) ) ))}

#compute log-likelihood for competing risks model
(loglikw <- llikweibullcr(theta=c(mu1, log(sigma1), mu2, log(sigma2)),
                          y=devG$Kilocycles, d1=devG$censS, d2=devG$censW))

#log-likelihood function (needed for profiling the life time quantiles)
mu1SE <- sqrt(diag(vcov(modS)))[1] #extract standard error for mu1

llikweibullcrtp <- function (theta, y, d1, d2, tpHat) {
  sigma1 <- exp(theta[1])
  mu2 <- theta[2]
  sigma2 <- exp(theta[3])
  mu1 <- optim(par=mu1, mu1Min, method="Brent",
               lower=mu1-25*mu1SE, upper=mu1+25*mu1SE,
               tpHat=tpHat, sigma1=sigma1, mu2=mu2, sigma2=sigma2)$par
  
  z1 <- (log(y)-mu1)/sigma1
  z2 <- (log(y)-mu2)/sigma2
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log( ((( 1/sigma1*y) * dsev(z1))^d1)*( (1-psev(z1)) ^(1-d1) ) )) +
    sum(log( ((( 1/sigma2*y) * dsev(z2))^d2)*( (1-psev(z2)) ^(1-d2) ) ))}

#profile function (using the Nelder-Mead optimization method)
profile.crtp <- function(y, d1, d2, tpHats) {
  plktphat <- rep(0, length(tpHats))
  
  for (i in 1:length(tpHats)) {
    temp <- fminsearch(llikweibullcrtp, c(log(sigma1),mu2,log(sigma2)),
                       minimize=FALSE, dfree=TRUE,
                       y=y, d1=d1, d2=d2, tpHat=tpHats[i] )
    plktphat[i] <- temp$fval
  }
  plktphat
}


#MLE for tp

#function for computing life time quantile at a specified proportion (=p)
Qss <- function(t, mu1, sigma1, mu2, sigma2, p) {
  (1- (1-psev((log(t)-mu1)/sigma1)) * (1-psev((log(t)-mu2)/sigma2))) - p}

#compute MLE for life time quantile (at specified p)
(tpHat <- uniroot(Qss, c(0,1e+9), mu1=mu1, sigma1=sigma1, mu2=mu2, sigma2=sigma2,
                  p=p)$root)

#fixed values of life time quantile
tpHats <- seq(5e-3, 5, length.out=50)

#likelihood profile
ll <- profile.crtp(tpHats=tpHats, y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)

#plot profile likelihood
plot(tpHats, ll, type='l', xlab=bquote(t[.(p)]), ylab="log-likelihood")
alpha <- .95 #specify confidence level
limit <- loglikw - .5*qchisq(alpha, df=1)
abline(h=limit,col="blue",lty=2)
abline(v=tpHat,col="red",lty=2)

#compute limits of confidence interval by hand
limitlkh <- function(y, d1, d2, tps) {
  profile.crtp(y, d1, d2, tps) - limit
}

lowerci <- uniroot(limitlkh, c(5e-3,tpHat), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)
upperci <- uniroot(limitlkh, c(tpHat,5), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)
c(lowerci$root, upperci$root)


#include confidence bounds in plot
abline(v=lowerci$root, col="darkgreen", lty=2)
abline(v=upperci$root, col="darkgreen", lty=2)





##failure time-models
#assuming that the failure times for both failure modes follow a lognormal distribution

#analysis of failure mode S only
modS <- survreg(Surv(Kilocycles, censS) ~ 1, data=devG, dist="lognormal")
summary(modS)

#analysis of failure mode W only
modW <- survreg(Surv(Kilocycles, censW) ~ 1, data=devG, dist="lognormal")
summary(modW)

#extract MLEs for model parameters
mu1 <- coef(modS)
sigma1 <- modS$scale
mu2 <- coef(modW)
sigma2 <- modW$scale



##likelihood based confidence intervals for life time quantiles tp
#tp is the time (t) at which a specified proportion (p) of the units fails

#specify proportion
p <- .01

#function for computing abs(p-f(mu1,sigma1,mu2,sigma2)),
#where p is the specified proportion, and f(.) some function
#of mu1, sigma1, mu2, and sigma2
mu1Min <- function (mu1Par, sigma1, mu2, sigma2, tpHat){
  z1 <- (log(tpHat)-mu1Par)/sigma1
  z2 <- (log(tpHat)-mu2)/sigma2
  
  abs(p - (1- (1-pnorm(z1)) * (1-pnorm(z2)) ) )}

#log-likelihood function of competing risks model
llikweibullcr <- function (theta, y, d1, d2) {
  mu1 <- theta[1]
  sigma1 <- exp(theta[2])
  mu2 <- theta[3]
  sigma2 <- exp(theta[4])
  
  z1 <- (log(y)-mu1)/sigma1
  z2 <- (log(y)-mu2)/sigma2
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log( ((( 1/sigma1*y) * dnorm(z1))^d1)*( (1-pnorm(z1)) ^(1-d1) ) )) +
    sum(log( ((( 1/sigma2*y) * dnorm(z2))^d2)*( (1-pnorm(z2)) ^(1-d2) ) ))}

#compute log-likelihood of competing risks model
(loglikw <- llikweibullcr(theta=c(mu1, log(sigma1), mu2, log(sigma2)),
                          y=devG$Kilocycles, d1=devG$censS, d2=devG$censW))

#log-likelihood function (needed for profiling the life time quantiles)
mu1SE <- sqrt(diag(vcov(modS)))[1] #extract standard error for mu1

llikweibullcrtp <- function (theta, y, d1, d2, tpHat) {
  sigma1 <- exp(theta[1])
  mu2 <- theta[2]
  sigma2 <- exp(theta[3])
  mu1 <- optim(par=mu1, mu1Min, method="Brent",
               lower=mu1-25*mu1SE, upper=mu1+25*mu1SE,
               tpHat=tpHat, sigma1=sigma1, mu2=mu2, sigma2=sigma2)$par
  
  z1 <- (log(y)-mu1)/sigma1
  z2 <- (log(y)-mu2)/sigma2
  
  #log-likelihood for exact (not censored) and right-censored observations
  sum(log( ((( 1/sigma1*y) * dnorm(z1))^d1)*( (1-pnorm(z1)) ^(1-d1) ) )) +
    sum(log( ((( 1/sigma2*y) * dnorm(z2))^d2)*( (1-pnorm(z2)) ^(1-d2) ) ))}

#profile function (using the Nelder-Mead optimization method)
profile.crtp <- function(y, d1, d2, tpHats) {
  plktphat <- rep(0, length(tpHats))
  
  for (i in 1:length(tpHats)) {
    temp <- fminsearch(llikweibullcrtp, c(log(sigma1),mu2,log(sigma2)),
                       minimize=FALSE, dfree=TRUE,
                       y=y, d1=d1, d2=d2, tpHat=tpHats[i] )
    plktphat[i] <- temp$fval
  }
  plktphat
}

#MLE for tp

#function for computing life time quantile at a specified proportion (=p)
Qss <- function(t, mu1, sigma1, mu2, sigma2, p) {
  (1- (1-pnorm((log(t)-mu1)/sigma1)) * (1-pnorm((log(t)-mu2)/sigma2))) - p}

#compute MLE for life time quantile (at specified p)
(tpHat <- uniroot(Qss, c(0,1e+9), mu1=mu1, sigma1=sigma1, mu2=mu2, sigma2=sigma2,
                  p=p)$root)

#fixed values of life time quantile
tpHats <- seq(5e-3, 10, length.out=50)

#likelihood profile
ll <- profile.crtp(tpHats=tpHats, y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)

#plot profile likelihood
plot(tpHats, ll, type='l', xlab=bquote(t[.(p)]), ylab="log-likelihood")
alpha <- .95 #specify confidence level
limit <- loglikw - .5*qchisq(alpha, df=1)
abline(h=limit,col="blue",lty=2)
abline(v=tpHat,col="red",lty=2)

#compute limits of confidence interval by hand
limitlkh <- function(y, d1, d2, tps) {
  profile.crtp(y, d1, d2, tps) - limit
}

lowerci <- uniroot(limitlkh, c(5e-3,tpHat), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)
upperci <- uniroot(limitlkh, c(tpHat,10), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)
c(lowerci$root, upperci$root)


#include confidence bounds in plot
abline(v=lowerci$root, col="darkgreen", lty=2)
abline(v=upperci$root, col="darkgreen", lty=2)