R code for constructing likelihood based confidence intervals for stress-strength models

In my previous blog post I demonstrated how to compute bootstrap confidence intervals for stress-strength models.
The following R-code may be used for computing likelihood based confidence intervals for stress-strength models.

Definition of the unreliability in case of a stress-strength model
A component may fail when the stress (or load) exceeds the strength. Accordingly, the unreliability of the component is defined as:

U_{hat} =\int^{+\infty}_{-\infty} f_l(l)F_s(l) \ dl

where Uhat is the unreliability of the component, fl(·) the probability density function of the load, and Fs(·) the cumulative density function of the strength.

Note that both the load and strength random variables never take negative values. Since both these random variables are non-negative, the lower limit of the above integral changes. That is, the lower limit changes from -∞ to 0, since fl(·) and Fs(·) will be 0 for negative values (i.e., for l<0). As a consequence, the above integral is given by:

U_{hat} =\int^{+\infty}_{0} f_l(l)F_s(l) \ dl

The reliability of the component is Rhat=1-Uhat.

Computing likelihood based intervals for Uhat
For computing the likelihood, a location-scale parameterization for the load and strength distributions will be employed. The location of the distribution for the load random variable will be denoted by μl, and its scale by σl. In addition, μs and σs are, respectively, the location and scale of strength random variable.
A consequence of using a location-scale based distribution is that fl(·) in the above integral becomes a function of μl and σl. Similarly, Fs(·) will be some function of μs and σs.
Furthermore, the likelihood function for the load (denoted by Ll) can now be defined in terms of μl and σl, and the likelihood function for the strength (denoted by Ls) in terms of μs and σs.
Finally, under the assumption that the load and strength random variables are independent, the likelihood function for stress-strength model is given by Lss=Ll(μl, σl)*Ls(μs, σs).

The likelihood function Lss will be used in computing the likelihood based intervals for Uhat. However, for computing these Uhat confidence intervals we need to reparameterize Lss. For this reparameterization to work, remember that Uhat was defined as some function of the parameters μl, σl, μs, and σs , denoted by Uhat=f(μl, σl, μs, σs). As a consequence, it is also possible to express μl as some function of Uhat, σl, μs, and σs , denoted by μl=g(Uhat, σl, μs, σs).

The reparameterization of Lss requires that μl is substituted by g(Uhat, σl, μs, σs) in the original likelihood function for Lss. In that way, Lss becomes a function of Uhat, σl, μs, σs. This latter parameterization of Lss makes it possible to construct likelihood based confidence intervals for Uhat.

For our specific reparameterization of Lss, it is thus needed to express μl as a function of the parameters Uhat, σl, μs, and σs. In other words, we need to compute μl given any values for Uhat, σl, μs, and σs.
The definition of Uhat showed that Uhat=f(μl, σl, μs, σs), and therefore that Uhat f(μl, σl, μs, σs) = 0. This latter equation may yield our sought for value for μl, given any values for the remaining parameters in the equation.
One way of obtaining this value for μl is by relying on numerical optimization. More specifically, μl can be found for any given values Uhat, σl, μs, and σs by minimizing abs(Uhat – f(μl, σl, μsσs)) over all possible values of μl.stress strength model likelihood ci

The following R code focuses at first on a situation in which the stress (or load) and strength both follow a normal distribution. However, any other combination of distributions may be used in the R code (e.g, having the load follow a lognormal distribution, while strength a Weibull distribution). In fact, it will be demonstrated how to change the code such that the stress and strength both follow a Weibull distribution.

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(StressStrength)


##stress and strength both follow a normal distribution

##data

#generate artificial data

#stress data
stress <- rnorm(20, 5000, 190)

#strength data
strength <- rnorm(20, 6000, 600)


#fit stress and strength distribution (note: both follow a normal distribution)
sts <- survreg(Surv(stress)~1, dist="gaussian")
str <- survreg(Surv(strength)~1, dist="gaussian")

#obtain MLE for mu and sigma
#Lmu
stsMu <- coef(sts)[1]
names(stsMu) <- "stsMu"
#Lsigma
stsSigma <- sts$scale
names(stsSigma) <- "stsSigma"

#Smu
strMu <- coef(str)[1]
names(strMu) <- "strMu"
#Ssigma
strSigma <- str$scale
names(strSigma) <- "strSigma"


##compute likelihood based interval

#specify integrand for numerical integration
integrand <- function(x, mu1, mu2, s1, s2) {
  #pdf for stress
  f1 <- dnorm((x-mu1)/s1) / s1
  #cdf for strength
  f2 <- pnorm((x-mu2)/s2)
  f1*f2
}

#estimate of unreliability using numerical integration
uplimit <- qnorm(.99999999)*sts$scale+coef(sts) #upper limit of the integral
Uhat <- integrate(integrand, 0, uplimit, mu1=coef(sts), mu2=coef(str),
                  s1=sts$scale, s2=str$scale)$value
#estimate of Reliability (=Rhat)
(Rhat <- 1-Uhat)

#compare the numerical solution of the integral with the exact solution
SSR(c(coef(str), str$scale), c(coef(sts), sts$scale))


#function for computing abs(Uhat-f(stsMu,stsSigma,strMu,stsSigma))
stsMuMin <- function (stsMuPar, stsSigma, strMu, strSigma, Uhat){
  uplimit <- qnorm(.99999999)*stsSigma+stsMuPar
  abs(Uhat-integrate(integrand, 0, uplimit, mu1=stsMuPar, mu2=strMu,
                     s1=stsSigma, s2=strSigma)$value)}

#log-likelihood function
stsMuSE <- sqrt(diag(vcov(sts)))[1] #extract standard error for stsMu

llik.ss <- function (modpar, stsObs, strObs, Uhat) {
  stssigma <- exp(modpar[1])
  strsigma <- exp(modpar[2])
  strmu <- modpar[3]
  stsmu <- optim(par=stsMu, stsMuMin, method="Brent",
                 lower=stsMu-20*stsMuSE, upper=stsMu+20*stsMuSE,
                 Uhat=Uhat, stsSigma=stssigma, strMu=strmu, strSigma=strsigma)$par
  
  z1 <- (stsObs-stsmu)/stssigma
  z2 <- (strObs-strmu)/strsigma
  
  #assuming that the stress and strength random variables are independent,
  #the log-likelihood function is given by:
  sum(log(dnorm(z1)/stssigma)) +
    sum(log(dnorm(z2)/strsigma))
}

#profile function
profile.ss <- function(stsObs, strObs, Uhats) {
  plkuhat <- rep(0, length(Uhats))
  
  for (i in 1:length(Uhats)) {
    temp <- optim(c(log(stsSigma), log(strSigma), strMu), llik.ss,
                  method='Nelder-Mead', control=list(fnscale=-1),
                  stsObs=stsObs, strObs=strObs, Uhat=Uhats[i] )
    plkuhat[i] <- temp$value
  }
  plkuhat
}

#MLE for Rhat
Rhat

#generate fixed values for Rhat
Rhats <- seq(.75, .999, length.out=25)
Uhats <- 1-Rhats

#likelihood profile (=computationally intensive)
#note: in case of warning messages, or irregular profile plot,
#decrease the width of the interval for possible values of Xmu
#in the llik.ss function
#that is, change the bounds to, for instance, stsMu+/-10*stsMuSE
ll <- profile.ss(Uhats=Uhats, stsObs=stress, strObs=strength)

#plot profile likelihood for Rhat
plot(Rhats,ll, type='l', xlab="Reliability", ylab="log-likelihood")

#include lower limit for log-likelihood values
loglikw <- llik.ss(modpar=c(log(stsSigma), log(strSigma), strMu),
                   stsObs=stress, strObs=strength, Uhat=Uhat)
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for Rhat
abline(v=Rhat, col="red", lty=2)


#compute limits of confidence interval by hand (=computationally intensive)
limitlkh <- function(stsObs, strObs, Uhats) {
  profile.ss(stsObs, strObs, Uhats) - limit
}

lowerci <- uniroot(limitlkh, c(Uhat, 1-.75), stsObs=stress, strObs=strength)
upperci <- uniroot(limitlkh, c(1-1, Uhat), stsObs=stress, strObs=strength)
round(c(1-lowerci$root, 1-upperci$root), 8)

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


#compare this confidence interval with that computed by the Reiser-Guttman method
estSSR(strength, stress, type="RG", alpha=.05)$CI
#the Reiser-Gutmann method yields on average good results
#that is, the coverage probability of these intervals is often very close to the nominal level


#lower bound of one sided interval
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-2*alpha, df=1)

lowerci <- uniroot(limitlkh, c(Uhat, 1-.75), stsObs=stress, strObs=strength)
1-lowerci$root

#compare confidence interval with that computed by the Guo-Krishnamoorthy method
estSSR(strength, stress, type="GK", alpha=.05, twoside=FALSE)$CI
#this Guo-Krishnamoorthy method yields on average excellent results
#that is, the coverage probability of these intervals is often very close to the nominal level





##stress and strength both follow a Weibull distribution

##data

#generate artificial data

#stress data
stress <- exp(7 + .3*rsev(40))

#strength data
strength <- exp(8 + .6*rsev(30))


#fit stress and strength distribution (note: both follow a Weibull distribution)
sts <- survreg(Surv(stress)~1, dist="weibull")
str <- survreg(Surv(strength)~1, dist="weibull")

#obtain MLE for mu and sigma
#Lmu
stsMu <- coef(sts)[1]
names(stsMu) <- "stsMu"
#Lsigma
stsSigma <- sts$scale
names(stsSigma) <- "stsSigma"

#Smu
strMu <- coef(str)[1]
names(strMu) <- "strMu"
#Ssigma
strSigma <- str$scale
names(strSigma) <- "strSigma"



##compute likelihood based interval

#specify integrand for numerical integration
integrand <- function(x, mu1, mu2, s1, s2) {
  #pdf for stress
  f1 <- dsev((log(x)-mu1)/s1) / (s1*x)
  #cdf for strength
  f2 <- psev((log(x)-mu2)/s2)
  f1*f2
}

#estimate of unreliability using numerical integration
uplimit <- exp(qsev(.99999999)*sts$scale+coef(sts)) #upper limit of the integral
Uhat <- integrate(integrand, 0, uplimit, mu1=coef(sts), mu2=coef(str),
                  s1=sts$scale, s2=str$scale)$value
#estimate of Reliability (=Rhat)
(Rhat <- 1-Uhat)


#function for computing abs(Uhat-f(stsMu,stsSigma,strMu,stsSigma))
stsMuMin <- function (stsMuPar, stsSigma, strMu, strSigma, Uhat){
  uplimit <- exp(qsev(.99999999)*stsSigma+stsMuPar)
  abs(Uhat-integrate(integrand, 0, uplimit, mu1=stsMuPar, mu2=strMu,
                     s1=stsSigma, s2=strSigma)$value)}

#log-likelihood function
stsMuSE <- sqrt(diag(vcov(sts)))[1] #extract standard error for stsMu

llik.ss <- function (modpar, stsObs, strObs, Uhat) {
  stssigma <- exp(modpar[1])
  strsigma <- exp(modpar[2])
  strmu <- modpar[3]
  stsmu <- optim(par=stsMu, stsMuMin, method="Brent",
                 lower=stsMu-20*stsMuSE, upper=stsMu+20*stsMuSE,
                 Uhat=Uhat, stsSigma=stssigma, strMu=strmu, strSigma=strsigma)$par
  
  z1 <- (log(stsObs)-stsmu)/stssigma
  z2 <- (log(strObs)-strmu)/strsigma
  #assuming that the stress and strength random variables are independent,
  #the log-likelihood function is given by:
  sum(log(dsev(z1)/(stsObs*stssigma))) +
    sum(log(dsev(z2)/(strObs*strsigma)))
}

#profile function
profile.ss <- function(stsObs, strObs, Uhats) {
  plkuhat <- rep(0, length(Uhats))
  
  for (i in 1:length(Uhats)) {
    temp <- optim(c(log(stsSigma), log(strSigma), strMu), llik.ss,
                  method='Nelder-Mead', control=list(fnscale=-1),
                  stsObs=stsObs, strObs=strObs, Uhat=Uhats[i] )
    plkuhat[i] <- temp$value
  }
  plkuhat
}

#MLE for Rhat
Rhat

#generate fixed values for Rhat
Rhats <- seq(.60, .99, length.out=25)
Uhats <- 1-Rhats

#likelihood profile (=computationally intensive)
#note: in case of warning messages, or irregular profile plot,
#decrease the width of the interval for possible values of Xmu
#in the llik.ss function
#that is, change the bounds to, for instance, stsMu+/-10*stsMuSE
ll <- profile.ss(Uhats=Uhats, stsObs=stress, strObs=strength)

#plot profile likelihood for Rhat
plot(Rhats,ll, type='l', xlab="Reliability", ylab="log-likelihood")

#include lower limit for log-likelihood values
loglikw <- llik.ss(modpar=c(log(stsSigma), log(strSigma), strMu),
                   stsObs=stress, strObs=strength, Uhat=Uhat)
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-alpha, df=1)
abline(h=limit,col="blue",lty=2)
#include MLE for Rhat
abline(v=Rhat, col="red", lty=2)


#compute limits of confidence interval by hand (=computationally intensive)
limitlkh <- function(stsObs, strObs, Uhats) {
  profile.ss(stsObs, strObs, Uhats) - limit
}

lowerci <- uniroot(limitlkh, c(Uhat, 1-.60), stsObs=stress, strObs=strength)
upperci <- uniroot(limitlkh, c(1-.99, Uhat), stsObs=stress, strObs=strength)
round(c(1-lowerci$root, 1-upperci$root), 8)

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


#lower bound of one sided interval
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-2*alpha, df=1)

lowerci <- uniroot(limitlkh, c(Uhat, 1-.65), stsObs=stress, strObs=strength)
1-lowerci$root