Incorporating interval-censored event times into standard survival models

Question

I am trying to model the time until some event occurs for individuals observed over a 24 month period.

For about 75% of people, no event occurs. For 15% of people, we know exact time of the event. For the other 10%, we only know a time window in which it occurs. So some sample "survival times" might be:

>24 
2.5 
>24 
5.0 
0 to 6
6 to 12 
>24 
18 to 24 

and so on. Can this type of data be accommodated within "standard" survival models? If so, some hint at how to construct the Surv() object in R and how to pass it to survreg would be appreciated.

Answer 1

I solved this problem. Yes, you can do this within classical models and it is simple to do this using the survival package in R.

In constructing the Surv() object, you need to use the code 3 for the event argument to Surv(), in conjunction with type='interval'. Here's a simple example where the first 50 observations have no "event" within 24 months, the second 50 have exactly known event times and the last 50 only have 6-month windows:

x <- rnorm(150)
time <- c(rep(24,50), ceiling(24*runif(50)), ceiling(18*runif(50)))
time2 <- c(rep(24,50), time[51:100], time[101:150]+6)
status <- c( rep(0,50), rep(1,50), rep(3,50) )  
a <- Surv(time, time2, status, type="interval")
model <- survreg(a ~ x, dist="exponential")

I did some simulations to verify that this does recover the population structure correctly, even if all of the event times are interval censored (p=1 below).

n <- 1000 # sample size
p <- .5 # proportion of event times (when event time is known) that are interval censored 

# simulate baseline waiting time of 50 and covariate effect of 1. 
t <- rexp(n, rate=1/50)
x <- rnorm(n) 
t <- t*exp(x) 

# censor everything after 24. 
w <- which(t>24) 
t[w] <- 24 
status <- rep(1,n)
status[w] <- 0

# randomly pick some ones to interval censor and generate a random window around them
w <- which(status==1) 
s <- sample(w,floor(p*length(w)))
status[s] <- 3 
ot <- t
t[s] <- floor( t[s]-6*runif(length(s)) ) 
t1 <- t 
t2 <- t
t2[s] <- t2[s] + 6 
t1[ which(t1<=0) ] <- 1e-10 
t2[ which(t2<=0) ] <- t1[which(t2<=0)] + 6

# fit the model. 
a <- Surv(t1, t2, status, type="interval")
model <- survreg(a ~ x, dist="exponential")
print( c(exp( coef(model)[1] ), coef(model)[2] ) ) # should be about 50 and 1.