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.
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.datall-analyse.nl/blog_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