Quantiles of the response variable in generalised linear model with binary outcome (survival analysis)

Question

I have a dataset that looks more or less like this (interval-censored outcome, multiple screening events per patient, not synchronised):

ID <- c(1,1,1,1,2,2,2,2,3,3,3)
sex <- c(rep("F", 8), rep("M", 3))
screen_time_since_diagnosis <- c(48, 400, 760, 900, 20, 360, 800, 1400, 0, 400, 600)
outcome <- c(0,0,0,1,0,0,0,0,0,0,1)
diabetes <- as.data.frame(cbind(ID, sex, screen_time_since_diagnosis, outcome), 
                          stringsAsFactors = T)

And am running a generalised linear model with complementary log-log link function like this:

model <- glm(outcome~ sex + screen_time_since_diagnosis,
    data = diabetes, 
    family = binomial(link = "cloglog"))

However, risk calculation is not my end goal. I want to know the time at which the risk of having an event for an individual patient is equal to a certain threshold. In Cox PH it can be done with the quantile.surfit method. There is a package called ciTools that allows one to calculate quantiles of the response variable of a glm object, but it doesn't work for binary variables (when I try to run it I get a message Error in add_quantile.glm(df = newdata, model, p = 0.2) : Prediction intervals for Bernoulli response variables aren't useful.

1. Why are prediction intervals for bernoulli response variables not useful?
2. Is there any other way to calculate this to get similar output to when using quantile() on survfit or Surv objects?

Many thanks!

Answer 1

The first question is straightforward. With a Bernoulli outcome at predicted probability $p$, the variance is $p(1-p)$. The simulation approach of ciTools has nothing to add. You could use the Wald confidence intervals around the linear-predictor value that you specify (based on the covariate values and their coefficients), to get CI around a predicted probability, which is probably what you want. The emmeans package might help with that.

The second question is more complicated. A discrete-time survival model directly evaluates the probability of an event during each time interval, conditional on having been at risk during that interval, as a function of the covariates in place. That's the discrete-time hazard, $h_{t_i}$. The survival curve is the product of the complements of those discrete-time hazards:

$$S(t) =P(T>t)=\prod_{t_i \le t}(1-h_{t_i})$$

I want to know the time at which the risk of having an event for an individual patient is equal to a certain threshold.

For that, you need to use a predict() function on the model, with new data that represent a complete person-period record for the individual. That is, the new data need to include a separate row for the hypothesized individual at each event time in the model. If you ask for the probability prediction based on those new data, you get a vector of hazards over time. Then apply the above formula (or the equivalent estSurv() function in the discSurv package) to get survival over time. Solve for the time at which the survival function equals the desired quantile.

That approach can be extended by similar evaluation of upper and lower confidence limits at each event time and finding when those confidence intervals intersect with the desired quantile. That's how quantile.survfit() calculates CI for a quantile.

Your data and non-synchronized follow-up times

The above holds for a properly constructed discrete-time survival model. The way you formatted your data isn't correct for that, and it's not clear that your type of data can be handled adequately by a discrete-time survival model.

Discrete-time survival data analyzed with a glm() binomial model need to be in person-period format. You need to start with all of the event times in the data, then make sure that each individual (while at risk) has a separate row in the data for all of those event times, with the time, status, and covariate values included in the row.

Your data format doesn't do that. Consider the event at time = 600 for ID3. Model evaluation at time = 600 needs to take into account the fact that both ID1 and ID2 are also at risk at that time, but neither individual has a data row for that value of time. Thus glm() has no way of knowing to include those individuals in the evaluation. Similarly for the event at time = 900 for ID1, glm() has no way of knowing to include the at-risk ID2 in the evaluation. You need to put your data into the correct person-period format.

Even then, with unsynchronized follow-up times, you might end up having to make assumptions about event times that make the analysis questionable. Say that ID4 has an event at time = 100 while ID5 has interval-censored survival with no event at time = 90 but an event before time = 110. When you set up the data row for ID5 at time = 100, will you treat that time for ID5 as no-event or as an event? With the unsynchronized observation times you can't make an informed decision.

Although discrete-time binomial models with complementary log-log links make sense for synchronized observation times, they thus can lose precision with unsynchronized observation times. You would be better off with modeling that takes unsynchronized interval censoring into account, as provided for example by the R icenReg package.