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

Question

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?

2. 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 

```