I solved this problem but cannot post the answer due to my low (<10) reputation. Instead, I'm posting the answer here.
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:
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")
g <- survreg(a ~ x, dist="exponential")
print( c(exp( coef(g)[1] ), coef(g)[2] ) ) # should be about 50 and 1.