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