R code for a Fatigue-Limit Model

The following R code implements a fatigue-limit model. Details of the implemented fatigue-limit model can be found in this 1997 paper by Pascual and Meeker, and in Meeker and Escobar’s 1998 book Statistical Methods for Reliability Data (pp. 593-597).

The implemented model assumes that the fatigue-limit is fixed (i.e., identical fatigue-limit for all specimens). But in some practical cases it may be more realistic to assume that the fatigue-limits of individual specimens are considerably different. Under such circumstances a random fatigue-limit model should be used, which allows for a different fatigue-limit for each specimen. A random fatigue-limit model is described by Pascual and Meeker in this 1998 paper. However, the implemented fixed fatigue-limit model can still be applied when the differences in individual fatigue limits between specimens are small.

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)

#Superalloy data (see: Meeker & Escobar, pp. 593-597)
alloy <- read.table("http://www.public.iastate.edu/~wqmeeker/anonymous/Stat533_data/Splida_text_data/superalloy.txt", header=T)
#failure=1, runout=0
alloy$status <- abs(as.numeric(alloy$Status)-2)
alloy$Cycles <- 1000*alloy$kCycles
head(alloy)

#log-log S-N curve
plot(x=alloy$kCycles, y=alloy$PseudoSress, log="xy", pch=alloy$status,
     xlab="Thousands of Cycles", ylab="Stress")
legend("topright", pch=c(1,0), legend=c("failure","runout"))



#obtain start values for Beta0mu and Beta0sigma of the fatigue limit model
#by fitting a model that ignores that mu and sigma are related to stress
init <- survreg(Surv(Cycles, status)~1,
                data=alloy, dist="lognormal")

mu <- coef(init)
names(mu) <- "mu"
sigma <- init$scale
names(sigma) <- "sigma"

#obtain start value for gamma (=fatigue limit) from log-log S-N curve
gamma <- 70

#log-likelihood function of fatigue limit model
llikfm <- function (beta, t, x, d) {
  mu <- beta[1] + beta[2]*log(x-beta[5])
  sigma <- exp(beta[3] + beta[4]*log(x))
  sum( d*(log(dnorm( ((log(t)-mu)/sigma))) - log(sigma*t) ) +
         (1-d)*log(1-pnorm(((log(t)-mu)/sigma))))
}

#optimization of log-likelihood
theta <- c(mu, 0, sigma, 0, gamma)
names(theta) <- c("Beta0mu","Beta1mu","Beta0sigma","Beta1sigma","Gamma")

fm.mle <- optim(c(theta), llikfm,
                method='BFGS', control=list(fnscale=-1),
                t=alloy$Cycles, d=alloy$status,
                x=alloy$PseudoSress)

#MLEs for parameters of fatigue limit model
#these MLEs are similar to those reported by Pascual & Meeker (1997) in Table 2
#(or see: Meeker & Escobar, Table 22.2, p. 596)
fm.mle$par
#log-likelihood model
(loglikModel <- fm.mle$value)



#log-log S-N curve including quantiles
plot(x=alloy$kCycles, y=alloy$PseudoSress, log="xy", pch=alloy$status,
     ylim=c(70,150), xlim=c(5,500),
     xlab="Thousands of Cycles", ylab="Stress")
legend("topright", pch=c(1,0), legend=c("failure","runout"))
#fatigue limit
abline(h=fm.mle$par[5], col="red", lty=2)
#percentiles
Stress <- seq(76, 150, by=1)
MuCycles <- fm.mle$par[1] + fm.mle$par[2]*log(Stress-fm.mle$par[5])
VarCycles <- exp(fm.mle$par[3] + fm.mle$par[4]*log(Stress))
#50th percentile (=median)
lines(x=exp(MuCycles + VarCycles*qnorm(.5))/1000, y=Stress, type="l",
      col="blue", lty=2)
#5th percentile
lines(x=exp(MuCycles + VarCycles*qnorm(.05))/1000, y=Stress, type="l",
      col="blue", lty=2)
#95th percentile
lines(x=exp(MuCycles + VarCycles*qnorm(.95))/1000, y=Stress, type="l",
      col="blue", lty=2)



#likelihood based confidence intervals for the parameters
#of the fatigue limit model

#calculate confidence interval for the gamma-parameter
#(the confidence intervals for the other 4 parameters of the model are
#are calculated in exactly the same way)

#log-likelihood for profiling gamma
llikp <- function (beta, t, x, d, fixed) {
  mu <- beta[1] + beta[2]*log(x-fixed)
  sigma <- exp(beta[3] + beta[4]*log(x))
  sum( d*(log(dnorm( ((log(t)-mu)/sigma))) - log(sigma*t) ) +
         (1-d)*log(1-pnorm(((log(t)-mu)/sigma))))
}

#profile function for gamma

#start values for theta:
#theta[1]=Beta0mu, theta[2]=Beta1mu, theta[3]=Beta0sigma, theta[4]=Beta1sigma
theta <- c(mu, 0, sigma, 0)

profile.gamma <- function(t, x, d, gammas) {
  plkmu <- rep(0, length(gammas))
  for (i in 1:length(gammas)) {
    temp <- optim(theta, llikp,
                  method='BFGS', control=list(fnscale=-1),
                  t=t, x=x, d=d, fixed=gammas[i])
    plkmu[i] <- temp$value
  }
  plkmu
}

#MLE of gamma
fm.mle$par[5]
#fixed values for gamma
gammas <- seq(40, 80, .5)

#likelihood profile
ll <- profile.gamma(gammas, t=alloy$Cycles, x=alloy$PseudoSress,
                    d=alloy$status)

#plot profile likelihood for gamma including a 95% confidence interval
plot(gammas,ll, type='l', xlab=expression(gamma), ylab="log-likelihood")
#include lower limit for log-likelihood values
limit <- loglikModel - .5*qchisq(.95, df=1)
abline(h=limit,col="blue",lty=2)
#red lines represent the boundaries of the 95% likelihood based interval
#these boundaries were reported by Pascual & Meeker (1997) in Table 2
abline(v=c(49.98, 79.79), col="red", lty=2)

#calculate the 95% confidence interval by hand
limitlkh <- function(t, x, d, gammas) {
  profile.gamma(t, x, d, gammas) - limit
}

lowerci <- uniroot(limitlkh,c(45,55), t=alloy$Cycles, x=alloy$PseudoSress,
                   d=alloy$status)
upperci <- uniroot(limitlkh,c(72,80), t=alloy$Cycles, x=alloy$PseudoSress,
                   d=alloy$status)
c(lowerci$root, upperci$root)