R code for fitting a two-failure mode model

A system or component can fail in different ways. For instance, a device may fail due to 1) an electrical surge, or 2) wearout. When a system fails due to two such failure modes, it may be better to model the failure times with a two-failure mode model (which is also known as a competing risks model).

The following R code fits a two-failure mode model. This R code assumes that the failure times of both failure modes follow a Weibull distribution. However, the code can be easily adapted to implement other distributions as well (such as the lognormal distribution). Furthermore, the R code can be extended to include more than two failure modes.

two failure mode surival model

The R code also demonstrates how to construct normal-approximation and likelihood based intervals for the failure probabilities of a two-failure mode model. In addition, the code computes the median and mean life time for a two-failure mode model.

An introduction to multiple-failure mode models is given by Meeker and Escobar in their 1998 book Statistical methods for reliability data.

Procedures for computing normal-approximation and likelihood based intervals for the failure probabilities of a multiple-failure mode model can be found this paper by Hong and Meeker.

Finally, see this blog post for computing likelihood based intervals for the life time quantiles of a two-failure mode 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(Matrix)

##data

#Device G data (see: Meeker & Escobar, Example 15.6, pp. 383-385)
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)



##fit failure-time models to data

##censoring

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

#failure modes S and W combined
#1=failure due to mode S or W; 0=censored
devG$censSW <- ifelse(devG$Cens==1, 0, 1)

##failure-time models

#Weibull analysis of failure mode S only
modS <- survreg(Surv(Kilocycles, censS) ~ 1, data=devG, dist="weibull")
summary(modS)
#95% confidence interval for mu
confint(modS)
#95% confidence interval for sigma
SEparmsS <- sqrt(diag(vcov(modS))) #extract standard errors
exp(log(modS$scale) + c(1,-1)*qnorm(.05/2)*SEparmsS[2])

#Weibull analysis of failure mode W only
modW <- survreg(Surv(Kilocycles, censW) ~ 1, data=devG, dist="weibull")
summary(modW)
#95% confidence interval for mu
confint(modW)
#95% confidence interval for sigma
SEparmsW <- sqrt(diag(vcov(modW))) #extract standard errors
exp(log(modW$scale) + c(1,-1)*qnorm(.05/2)*SEparmsW[2])

#Weibull analysis of failure modes S and W combined
modSW <- survreg(Surv(Kilocycles, censSW) ~ 1, data=devG, dist="weibull")
summary(modSW)
#95% confidence interval for mu
confint(modSW)
#95% confidence interval for sigma
SEparmsSW <- sqrt(diag(vcov(modSW))) #extract standard errors
exp(log(modSW$scale) + c(1,-1)*qnorm(.05/2)*SEparmsSW[2])



##probability plot

#compute plotting positions observations
times <- devG$Kilocycles
event <- devG$censSW
timesSorted <- times[order(times)]
eventSorted <- event[order(times)]
timesi <- unique(timesSorted[which(eventSorted==1)])

cumSurvRaw <- survfit(Surv(times, event)~ 1)
cumSurv <- unique(rep(cumSurvRaw$surv, cumSurvRaw$n.event))
cumFail <- 1-cumSurv
lagFail <- c(0, cumFail[-length(cumFail)])
Prob <- .5*(cumFail+lagFail)
yi <- qsev(Prob)


#function for predicting life time quantiles (Weibull distribution)
Qp <- function (mu, sigma, p){exp(mu+qsev(p)*sigma)}

#function for predicting failure probabilities (Weibull distribution)
Ft <- function (mu, sigma, t){psev( (log(t)-mu) / sigma)}


#predict life time quantiles for failure mode S
timesQuantiles <- seq(0,1,length.out=200)
predTimeS <- sapply(timesQuantiles, Qp, mu=coef(modS), sigma=modS$scale)
piS <- qsev(timesQuantiles)

#predict life time quantiles for failure mode W
timesQuantiles <- seq(0,1,length.out=200)
predTimeW <- sapply(timesQuantiles, Qp, mu=coef(modW), sigma=modW$scale)
piW <- qsev(timesQuantiles)

#predict failure probabilities of 2-failure mode model
#note that if a component/system fails due to two (or more) independent failure modes,
#then the failure probabilities can be modeled with a series-system model
#(see Meeker and Escobar, p. 382)
failureTimes <- seq(1, 500,length.out=10)
FS <- sapply(failureTimes, Ft, mu=coef(modS), sigma=modS$scale)
FW <- sapply(failureTimes, Ft, mu=coef(modW), sigma=modW$scale)
#series-system prediction of failure probability
Fseries <- 1-(1-FS)*(1-FW)


#construct probability plot
plot(timesi, yi, log="x", xlab="Time", ylab="Linearized CDF", xlim=c(1, 2000))
#predictions of failure time model for mode S
lines(predTimeS, piS, col="blue")
#predictions of failure time model for mode W
lines(predTimeW, piW, col="darkgreen")
#predictions of 2-failure mode model
lines(failureTimes, qsev(Fseries), col="red")
legend("bottomright", lty=1, col=c("blue", "darkgreen", "red"),
       legend=c("Mode S", "Mode W", "2-failure mode model"))



##normal-approximation confidence intervals for failure probabilities of 2-failure mode model

#function for computing normal-approximation confidence intervals
#notes:
#-the 2-failure mode model is identical to a series-system model (see Meeker and Escobar, p. 382)
#-the current function assumes the failure times of both failure
# modes follow a Weibull distribution
#-these normal-approximation confidence intervals may display "bend-back" behavior in the tails
ciFseries <- function(t,mu1,sigma1,vcov1,mu2,sigma2,vcov2,alpha=.05) {
  #compute cumulative failure probability for series system
  Fseries <- 1-(1-Ft(t=t, mu=mu1, sigma=sigma1))*(1-Ft(t=t, mu=mu2, sigma=sigma2))
  
  #apply the delta method for "undoing" the log transformation of sigma1
  vcov1[1,2] <- vcov1[2,1] <- t(c(0,exp(log(sigma1))))%*%vcov1%*%c(1,0)
  vcov1[2,2] <- t(c(0,exp(log(sigma1))))%*%vcov1%*%c(0,exp(log(sigma1)))
  #apply the delta method for "undoing" the log transformation of sigma2
  vcov2[1,2] <- vcov2[2,1] <- t(c(0,exp(log(sigma2))))%*%vcov2%*%c(1,0)
  vcov2[2,2] <- t(c(0,exp(log(sigma2))))%*%vcov2%*%c(0,exp(log(sigma2)))
  #construct variance-covariance matrix for series-system
  vcovSeries <- as.matrix(bdiag(vcov1,vcov2))
  
  #take partial derivaties of the function Fseries=1-(1-F1)*(1-F2)
  #[where F1=psev((log(t)-mu1)/sigma1) (=cumulative failure probability mode 1),
  #and F1=psev((log(t)-mu2)/sigma2) (=cumulative failure probability mode 2)]
  #with respect to mu1
  dF.dmu1 <- -dsev((log(t)-mu1)/sigma1)*1/sigma1*(1-psev((log(t)-mu2)/sigma2))
  #with respect to sigma1
  dF.dsigma1 <- -dsev((log(t)-mu1)/sigma1)*((log(t)-mu1)/sigma1)*1/sigma1*(1-psev((log(t)-mu2)/sigma2))
  #with respect to mu2
  dF.dmu2 <- -dsev((log(t)-mu2)/sigma2)*1/sigma2*(1-psev((log(t)-mu1)/sigma1))
  #with respect to sigma2
  dF.dsigma2 <- -dsev((log(t)-mu2)/sigma2)*((log(t)-mu2)/sigma2)*1/sigma2*(1-psev((log(t)-mu1)/sigma1))
  #combine partial derivatives in vector
  dF.dtheta <- c(dF.dmu1, dF.dsigma1, dF.dmu2, dF.dsigma2)
  #apply delta method for computing the standard error
  SEFseries <- sqrt(t(dF.dtheta)%*%vcovSeries%*%dF.dtheta)
  
  #compute limits of confidence interval (using a logit transformation)
  w <- exp((qnorm(1-alpha/2)*SEFseries)/(Fseries*(1-Fseries)))
  lcl <- Fseries/(Fseries+(1-Fseries)*w)
  ucl <- Fseries/(Fseries+(1-Fseries)/w)
  
  #return values
  c(time=t, Fhat=as.vector(Fseries), Std.Err=SEFseries, lcl=lcl, ucl=ucl)}

#2-failure mode model's predictions of failure probabilities (including confidence intervals)
#notice the unrealistic "bend-back" behavior of the confidence interval in the upper tail
ts <- seq(1, 500, by=1)
Fs <- data.frame(t(sapply(ts, ciFseries,
                          mu1=coef(modS), sigma1=modS$scale, vcov1=vcov(modS),
                          mu2=coef(modW), sigma2=modW$scale, vcov2=vcov(modW))))
plot(timesi, Prob, log="x", ylim=c(0,1), xlim=c(1,500),
     xlab="Time", ylab="Cumulative failure probability")
lines(Fs$t, Fs$Fhat, col="blue")
lines(Fs$t, Fs$ucl, lty=2, col="red")
lines(Fs$t, Fs$lcl, lty=2, col="red")

#inspect failure probabilities at specified times (including confidence intervals)
ts <- c(2, 5, 10, 20, 50, 100, 200, 500)
Fs <- data.frame(t(sapply(ts, ciFseries,
                          mu1=coef(modS), sigma1=modS$scale, vcov1=vcov(modS),
                          mu2=coef(modW), sigma2=modW$scale, vcov2=vcov(modW))))
#display estimates
Fs



##likelihood based confidence intervals for failure probabilities of 2-failure mode model
#these confidence interval do not display "bend-back" behavior in the tails

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

#log-likelihood function of series-system (=2-failure mode model)
llikweibullSeries <- function (theta, y, d1, d2) {
  mu1 <- theta[1]
  sigma1 <- theta[2]
  mu2 <- theta[3]
  sigma2 <- 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 of series-system
loglikw <- llikweibullSeries(theta=c(mu1, sigma1, mu2, sigma2),
                             y=devG$Kilocycles, d1=devG$censS, d2=devG$censW)
loglikw


#log-likelihood function (needed for profiling the failure probabilities)
llikweibullSeriesF <- function (theta, y, d1, d2, te, Fs) {
  sigma1 <- theta[1]
  mu2 <- theta[2]
  sigma2 <- theta[3]
  mu1 <- log(te) - sigma1*qsev(1-((1-Fs)/(1-psev( (log(te)-mu2)/sigma2 ))) )
  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
profile.Fseries <- function(y, d1, d2, te, Fss) {
  plkmu <- rep(0, length(Fss))
  
  for (i in 1:length(Fss)) {
    temp <- optim(c(sigma1,mu2,sigma2), llikweibullSeriesF,
                  method='BFGS', control=list(fnscale=-1),
                  y=y, d1=d1, d2=d2, te=te, Fs=Fss[i] )
    plkmu[i] <- temp$value
  }
  plkmu
}

#specify time at which to compute the failure probability confidence interval
t <- 100

#MLE for failure probability of 2-failure mode model at a specified time
(FsMLE <- 1-(1-Ft(t=t, mu=mu1, sigma=sigma1))*(1-Ft(t=t, mu=mu2, sigma=sigma2)))
#fixed values of failure probabilities
Fss <- seq(.1, .6, .01)

#likelihood profile
ll <- profile.Fseries(Fss, y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)

#plot profile likelihood
plot(Fss, ll, type='l', xlab=paste("F(",t,")", sep=""), ylab="log-likelihood")
alpha <- .95 #specify confidence level
limit <- loglikw - .5*qchisq(alpha, df=1)
abline(h=limit,col="blue",lty=2)
abline(v=FsMLE,col="red",lty=2)

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

lowerci <- uniroot(limitlkh, c(.1,FsMLE), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)
upperci <- uniroot(limitlkh, c(FsMLE,.5), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)
c(lowerci$root, upperci$root)

#compare with normal-approximation confidence interval
round(ciFseries(t=100, mu1=coef(modS), sigma1=modS$scale, vcov1=vcov(modS),
                mu2=coef(modW), sigma2=modW$scale, vcov2=vcov(modW)), 3)



##compute confidence interval in upper tail

#specify time at which to compute the failure probability confidence interval
t <- 500

#MLE for failure probability of 2-failure mode model at a specified time
(FsMLE <- 1-(1-Ft(t=t, mu=mu1, sigma=sigma1))*(1-Ft(t=t, mu=mu2, sigma=sigma2)))
#fixed values of failure probabilities
Fss <- seq(.75, 1, length.out=100)

#likelihood profile
ll <- profile.Fseries(Fss, y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)

#plot profile likelihood
plot(Fss, ll, type='l', xlab=paste("F(",t,")", sep=""), ylab="log-likelihood")
alpha <- .95 #specify confidence level
limit <- loglikw - .5*qchisq(alpha, df=1)
abline(h=limit,col="blue",lty=2)
abline(v=FsMLE,col="red",lty=2)

#compute limits of confidence interval by hand
lowerci <- uniroot(limitlkh, c(.75,FsMLE), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)
upperci <- uniroot(limitlkh, c(FsMLE,1), y=devG$Kilocycles, d1=devG$censS, d2=devG$censW, te=t)
c(lowerci$root, upperci$root)

#compare with normal-approximation confidence interval
#note: this interval is unrealistic due to the "bend-back" behavior
round(ciFseries(t=500, mu1=coef(modS), sigma1=modS$scale, vcov1=vcov(modS),
                mu2=coef(modW), sigma2=modW$scale, vcov2=vcov(modW)),3)



##compute median and mean life time for 2-failure mode model

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

#1) median life time

#plot cumulative failure probability as a function of time for 2-failure mode model
t <- seq(1, 500, 1)
ps <- 1- (1-psev((log(t)-mu1)/sigma1)) *(1-psev((log(t)-mu2)/sigma2))
plot(t, ps, type="l", col="blue", log="x", ylab="Cumulative failure probability",
     xlab="Time")

#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 median life time (at p=.5)
medLT <- uniroot(Qss, c(0,500), mu1=mu1, sigma1=sigma1, mu2=mu2, sigma2=sigma2,
                 p=.5)$root
medLT #median life time


#2) mean life time (or MTTF)

#set up integrand for numerical integration
integrand <- function(x, mu1, mu2, s1, s2) {
  #survival function for mode S
  f1 <- psev((log(x)-mu1)/s1)
  #survival function for mode W
  f2 <- psev((log(x)-mu2)/s2)
  1-(1-(1-f1)*(1-f2))
}

#estimate of MTTF (using numerical integration)
integrate(integrand, 0, Inf, mu1=coef(modS), mu2=coef(modW),
          s1=modS$scale, s2=modW$scale)$value