Survival time intervals starting from 0- errors from Surv() and IcenReg

I'm trying to model the risk of diabetic retinopathy as a function of time between screening appointments. The outcome is interval-censored because we only know that the event occurred between the previous screening $$T_{1}$$ and the next one $$T_{2}$$ . Each appointment is accompanied by a number of covariates (that can vary between appointments).

The end goal is to be able to tell the time at which the conditional probability reaches a user-defined threshold of risk. i.e. if we accept that patients should be screened when their risk is 5%, when exactly should this patient be seen next?

The reason I want to model it as time between appointments instead of time from diabetes diagnosis is that retinopathy screening is iterative as it is highly dependent on how well the patient is controlling their biomarkers. In the sample data code, I've made 3 different ways to include time in the model.

1. as time from diagnosis to event
• I want to avoid this
2. time from the previous appointment to the current one, i.e. left (side of the) interval = 0 and right interval = time between the appointments
• ideally, I could use this
3. same as 2), but adding the previous right side of the interval or time between the previous appointments, or lag(r_event).
• I don't like this but is the only way I can think of to make the model "focus" on the current state and also work.

I have two issues that I believe:

1. IcenReg throws errors such as

Lapack routine dgsev: system is exactly singular: U[4,4] = 0

which I assume is an issue with intervals starting from 0?

1. IcenReg makes use of survival's Surv() of type "interval2". in 98% of all appointments the event hasn't occurred. This means that 98% of my intervals will be L: 0, R:Inf.

Is there a better way to approach this?

data <- data.frame()
patients <- 1:100
set.seed(421)

for( i in patients) {
e <- sample(5:15, 1)
sex <- sample(c("F", "M"), 1)
sex <- rep(sex, times = e)
hba1c <- rnorm(e, 8.5, 1.2)
age_diagnosis <- sample(40:80, 1)
age_diagnosis <- rep(age_diagnosis, e)
t_diag_event <- ceiling(rnorm(e, 380, 90)) # when using rnorm(e, 300, 40), there will be issues with plotting survival funciton for "two_event" times
t_event_event <- t_diag_event
t_diag_event <- t_diag_event + rnorm(1, 4*365, 365)
t_diag_event <- cumsum(t_diag_event)
outcome <- rep("No", e-1)
outcome <- c(outcome,
sample(c("Yes", "No"), prob = c(0.98, 0.02), 1))
id <- rep(i, e)
data <- rbind(data,
cbind(id, outcome, t_event_event, t_diag_event,
hba1c, sex, age_diagnosis))
}

data <- data %>%
mutate(across(c(t_event_event, t_diag_event, hba1c, age_diagnosis),
~ as.numeric(.)),
across(c(outcome, sex),
~as.factor(.))) %>%
group_by(id) %>%
mutate(consecutive_appt = row_number() + 1, # 1st appt at t = 0 is not included
l_event = 0,
r_event = t_event_event,
l_diag = t_diag_event- t_event_event,
r_diag = l_diag + t_event_event,
l_two_event = l_event + lag(r_event),
r_two_event = lag(r_event) + t_event_event,
ever_has_outcome= max(as.numeric(outcome)-1),
age_diagnosis = ifelse(ever_has_outcome == 1,
age_diagnosis- rnorm(1, 5, 2),
age_diagnosis),
hba1c = ifelse(ever_has_outcome == 1,
hba1c + rnorm(1, 2, 1),
hba1c)) %>%
ungroup() %>%
# Surv's type interval2
mutate(r_diag = ifelse(outcome == "No",
Inf,
r_diag),
r_event = ifelse(outcome == "No",
Inf,
r_event),
r_two_event = ifelse(outcome == "No",
Inf,
r_two_event)) %>%
dplyr::select(id,
starts_with("t"),
ends_with(c("two_event", "diag", "event")),
everything()) %>%
subset(consecutive_appt >2)

library(icenReg)

# Modelling using time from diagnosis to event
semi_p_diag <- ic_sp(Surv(l_diag, r_diag,
type = "interval2") ~ sex + age_diagnosis + hba1c,
data = data)

diag_baseline(semi_p_diag,
dists = c("exponential", "weibull", "gamma", "loglogistic"))

# Modelling using time from event to event
semi_p_event <- ic_sp(Surv(l_event, r_event,
type = "interval2") ~ sex + age_diagnosis + hba1c,
data = data)

diag_baseline(semi_p_event,
dists = c("exponential", "weibull", "gamma", "loglogistic"))

# Modelling using time from event to event, but where the time is defined as
# t_event_event = time between appointments
# l = 0 + lag(t_event_event)
# r = t_event_event + lag(t_event_event)
semi_p_two_event <- ic_sp(Surv(l_two_event, r_two_event,
type = "interval2") ~ sex + age_diagnosis + hba1c,
data = data)

diag_baseline(semi_p_two_event,
dists = c("exponential", "weibull", "gamma")) # what happened to my semi-parametric survival function? (in case of two_event times)

diag_baseline(semi_p_two_event,
dists = c("exponential", "weibull", "gamma",
"loglogistic")) # will throw an error when the times between events are more uniform

$$$$

• Could you please show the way that you are calling the function and samples of the data as they are formatted? I'm not sure that you are handling the interval2 format correctly, but I can't tell without seeing both a sample of labeled data and the corresponding function call together. It's not clear in particular why you have any intervals starting at 0. In the interval2 format, the first time value should be the right censoring time, and the second value should be Inf to notate that as a right censoring time. There are additional issues here, but let's fix the data format first.
– EdM
Commented Dec 6, 2022 at 21:52
• Also, please provide any such data and code in your edited question via the {} `code tool on the toolbar (equivalent to a block of lines each starting with 4 or more spaces). Unlike with images of data, that makes it possible for text-to-speech software to translate it, and that makes it possible to copy examples directly from the question.
– EdM
Commented Dec 7, 2022 at 14:45
• Thanks for the comments- hope this is what you were looking for @EdM! Commented Dec 7, 2022 at 17:31
• In cancer epidemiology and clinical trials, we handle the interval censoring issue by taking the event timepoint as being on the date of the positive screen. To get an unbiased estimate of incidence density, subjects without events are censored to the last evaluable scan, and subjects missing two consecutive scans prior to an event are censored to the last evaluable scan. Commented Dec 7, 2022 at 18:03