R code for modeling with left truncated and right/interval censored data

Burn-in testing is used to screen out units or systems with short lifetimes. Units or systems that survived a burn-in test may give rise to left truncated data that is either right or interval censored.

Left truncated and right censored data
Tobias and Trindade reported in their 2012 book Applied Reliability on field failure times of units that survived a burn-in test of 5000 hours (Example 4.6, p. 109). These field failure times represent an example of left truncation in combination with right censoring.

Left truncated and interval censored data
Meeker and Escobar described in their 1998 book Statistical Methods for Reliability Data a field-tracking study of units that survived a 1000 hours burn-in test (Example 11.11, pp. 269-270). The data in Meeker and Escobar’s study is an example of left truncation in combination with interval censoring.

Model fitting using maximum likelihood optimization

The R code fits a Weibull (or lognormal) model to left truncated data that is either right or interval censored. The fitting of these models is done by log-likelihood optimization (using the optim function in R).

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

Weibull/lognormal model fitted to left truncated and right censored data

library(survival)
library(SPREDA)

#data reported by Tobias & Trindade in Example 4.6 (p. 109)
failtimes <- c(454,620,974,1063,1265,1936,2351,2364,2976,3095,3381,3475,4055,
               5592,5641,5765,5822,7913,8660,8682,8999,9002,9712,9960,10338,
               11337,11368,11904,12210,12776,12859,14036,15000)
#truncation time
TruncTime <- 5000
#weights
count <- c(rep(1,32), 69)
#failure=1, censored=0
status <- c(rep(1,32), 0)

burnT <- data.frame(TruncTime, failtimesfield=failtimes, status, count)
burnT$failtimes <- burnT$failtimesfield+TruncTime


#obtain start values for mu and sigma by fitting a model that
#ignores the units came from a truncated population
init <- survreg(Surv(failtimes, status)~1,
                data=burnT, weights=count, dist="weibull")

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

#Weibull model fitted to the left truncated and right censored data
#for a definition of the used likelihood, see: Meeker & Escobar, p. 268, eq. (11.4)

#log-likelihood function
lliktruncr <- function (theta, y, d, wt=rep(1,length(y)), trunctime, dist) {
  mu <- theta[1]
  sigma <- theta[2]
  
  z <- (log(y)-mu)/sigma
  ztrunc <- (log(trunctime)-mu)/sigma
  
  if (dist=="lognormal") {
    sum(wt*log(( (( 1/(sigma*y) * dnorm(z))/(1-pnorm(ztrunc)) )^d)*
                 (( (1-pnorm(z))/(1-pnorm(ztrunc)))^(1-d) )))}
  else if (dist=="weibull") {
    sum(wt*log(( (( 1/(sigma*y) * dsev(z))/(1-psev(ztrunc)) )^d)*
                 (( (1-psev(z))/(1-psev(ztrunc)))^(1-d) )))}
}

#optimization of log-likelihood
weibull.mle <- optim(c(mu, sigma), lliktruncr,
                     method='BFGS', control=list(fnscale=-1),
                     y=burnT$failtimes, d=burnT$status,
                     wt=burnT$count, trunctime=burnT$TruncTime,
                     dist="weibull",
                     hessian=TRUE)

#MLE for mu and sigma
weibull.mle$par
#variance covariance matrix for mu and sigma
solve(-1*weibull.mle$hessian)
#Tobias & Trindade obtained a MLE for alpha of 40432, and for beta of 0.745
exp(weibull.mle$par[1])
1/weibull.mle$par[2]

Weibull/lognormal model fitted to left truncated and interval censored data

library(survival)
library(SPREDA)


##data
#data reported by Meeker & Escobar in Example 11.11 (pp. 269-270)
#lower limit of time interval
HoursL <- c(1000,2000,5000,6000,7000,8000,9000,10000,11000)
#upper limit of time interval
HoursU <- c(2000,5000,6000,7000,8000,9000,10000,11000,NA)
#failure=1, censored=0
Status <- c(1,1,1,1,1,1,1,1,0)
#weights
Count <- c(2,5,6,11,7,14,10,14,4924)
#truncation time
TruncTime <- 1000

burnT <- data.frame(HoursL,HoursU,Status,Count,TruncTime)



##Weibull model
#obtain start values for mu and sigma by fitting a model that
#ignores the units came from a truncated population
mod1 <- survreg(Surv(HoursL, HoursU, type="interval2") ~ 1,
                weights=Count, data=burnT,
                dist="weibull")

mu <- coef(mod1)[1]
names(mu) <- "mu"
sigma <- mod1$scale
names(sigma) <- "sigma"



#Weibull model fitted to the left truncated and interval censored data
#for a definition of the used likelihood, see: Meeker & Escobar, p. 268

#log-likelihood function
lliktrunci <- function (theta, yl, yu, d, wt=rep(1,length(yl)), trunctime, dist) {
  mu <- theta[1]
  sigma <- exp(theta[2])
  
  zl <- (log(yl)-mu)/sigma
  zu <- (log(yu)-mu)/sigma
  ztrunc <- (log(trunctime)-mu)/sigma
  
  if (dist=="lognormal") {
    sum(wt*log(( ( (pnorm(zu)-pnorm(zl))/(1-pnorm(ztrunc)) )^d)*
                 (( (1-pnorm(zl))/(1-pnorm(ztrunc)))^(1-d) )))}
  else if (dist=="weibull") {
    sum(wt*log(( ( (psev(zu)-psev(zl))/(1-psev(ztrunc)) )^d)*
                 (( (1-psev(zl))/(1-psev(ztrunc)))^(1-d) )))}
}

#optimization of log-likelihood
#log transformation of sigma ensures that sigma takes a positive value
weibull.mle <- optim(c(mu, log(sigma)), lliktrunci,
                     method='BFGS', control=list(fnscale=-1),
                     yl=burnT$HoursL, yu=burnT$HoursU, d=burnT$Status,
                     wt=burnT$Count, trunctime=burnT$TruncTime,
                     dist="weibull")

#MLE of mu
(muMLE <- weibull.mle$par[1])
#MLE of sigma
(sigmaMLE <- exp(weibull.mle$par[2]))


#F(t), compare with Meeker & Escobar's estimate for F(1000) of 3.8e-5 (p.269)
round(pweibull(1000, scale=exp(muMLE), shape=1/sigmaMLE), 6)





##lognormal model
#obtain start values for mu and sigma by fitting a model that
#ignores the units came from a truncated population
mod2 <- survreg(Surv(HoursL, HoursU, type="interval2") ~ 1,
                weights=Count, data=burnT, dist="lognormal")
summary(mod2)

mu <- coef(mod2)[1]
names(mu) <- "mu"
sigma <- mod2$scale
names(sigma) <- "sigma"



#lognormal model fitted to the left truncated and interval censored data
lognormal.mle <- optim(c(mu, log(sigma)), lliktrunci,
                       method='BFGS', control=list(fnscale=-1),
                       yl=burnT$HoursL, yu=burnT$HoursU, d=burnT$Status,
                       wt=burnT$Count, trunctime=burnT$TruncTime,
                       dist="lognormal")

#MLE of mu
(muMLE <- lognormal.mle$par[1])
#MLE of sigma
(sigmaMLE <-exp(lognormal.mle$par[2]))


#F(t), compare with Meeker & Escobar's estimate for F(1000) of 9.5e-6 (p.269)
round(plnorm(1000, meanlog=muMLE, sdlog=sigmaMLE), 7)