Solved – Competing risk survival analysis with time-dependent covariates

rsurvival

Can anyone recommend an R package that handles left-truncation, right censoring, AND time-dependent covariates?

I have a data set that consists of distance-to-event data rather than time-to-event as seen in traditional survival analysis. I am studying the disturbance effects of cruise ship travel in a Bay in Alaska to a small seabird that floats on the water when resting or feeding. In reaction to an approaching ship a bird can either 1.Flush from the water 2.Dive underwater 3. allow the ship to pass (treated as censored).

Since birds are selected at a certain distance from the ship and this distance decreases as the ship advances I have transformed the distances measured so that at the ship progresses the distance increases rather than decreases (Subtracted each distance from the furthest distance observed for all birds). The distances at which the bird ultimately reacts will then be back transformed when presenting the results. This allows the data to be passed to typical Survival Analysis software as if it was time-to-event (increasing each repeated observation).

I am exploring the effects of multiple different fixed covariates (weather, distance to shore, ship speed, etc.) but have one time (or distance rather) dependent covariate, the birds bearing off the bow in relation to the ships course. Unless the bird is directly in the line of the ship its relative bearing off the ships course will change as the ship approaches closer. Eventually if the bird remains on the water it will be at 90 degrees tot he bow and the ship will pass by the bird. Therefore I have repeated distance and bearing measurements on each bird allowing this variable to change over the course of the observation. I have considered this covariate to be binary (0-44 degrees and 45-90 degrees).

I am mainly interested in how the covariates effect the competing events (Flush or Dive) and whether these two avoidance strategies change based on distance to the ship (preliminary analysis appears that they do).

Best Answer

The frailtypack package in its last version allows you to consider time-varying covariates effects in Cox, shared and joint models. Maybe It could help.

Related Solutions

R – Visualize Survival Analysis with Time-Dependent Covariates

I haven't been able to gain access to the Simon and Makuch article mentioned above, but having researched the topic I found:

Steven M Snapinn, Qi Jiang & Boris Iglewicz (2005) Illustrating the Impact of a Time-Varying Covariate With an Extended Kaplan-Meier Estimator, The American Statistician, 59:4, 301-307.

That article proposes a time-dependent Kaplan-Meier plot (KM) by simply updating the cohorts at all event times. It also cites the Simon and Makuch article for proposing a similar idea. The regular KM does not allow this, it only allows a fixed division into groups. The proposed method actually splits survival time according to covariate status - just like one could do when estimating a Cox model with piecewise constant covariates. For the Cox model this is a viable idea, and a standard one. It is, however, more intricate when doing a KM plot. Let me illustrate it with a simulation example.

Let's assume that we have no censoring, but some event (e.g., giving birth) that might or might not occur before time of death. Let's also assume constant hazards for the sake of simplicity. We will also assume that giving birth does not alter the hazard of dying. We will now follow the procedure prescribed in the above article. The article clearly states how this is done in R, simply split your subjects on time of giving birth such that they're constant in your grouping variable. Then use the counting process formulation in the Surv function. In code

library(survival)
library(ggplot2)
n          <- 10000
data       <- data.frame(id = seq(n), 
                         preg = rexp(n, 1), 
                         death = rexp(n, .5), 
                         enter = 0, 
                         per = NA, 
                         event = 1)
data$exit  <- data$death

data0      <- data
data0$exit <- with(data, pmin(preg, death))
data0$per  <- 0
data0$event[with(data0, preg < death)] <- 0

data1       <- subset(data, preg < death)
data1$enter <- data1$preg
data1$per   <- 1

data <- rbind(data0, data1)
data <- data[order(data$id), ]

Sfit <- survfit(Surv(time = enter, time2 = exit, event = event) ~ per, data = data)
autoplot(Sfit, censSize = 0)$plot

I'm more or less splitting it "by hand". We could use survSplit as well. The procedure actually gives me a very nice estimate.

No bias

We get almost identical estimates for the two groups as we should. But actually, my simulation was perhaps a bit unrealistic. Let's say a woman cannot give birth in the first two time units for some reason. This is at least reasonable in your example: there'll be some time between two pregnancies corresponding to the same woman. Making a small addition to the code

data <- data.frame(id = seq(n), 
                   preg = rexp(n, 1) + 2, 
                   death = rexp(n, .5), 
                   enter = 0, 
                   preg = NA, 
                   event = 1)

we get the following plot:

Bias

The same thing would be happening with your data. You will not see any 3rd pregnancies for at least some initial period of time, meaning that your estimate will be 1 for that group and that period of time. This is in my opinion a misrepresentation of your data. Consider my simulation. The hazards are identical, but for every time point the per1 estimate is greater than the per0 estimate.

You could consider different remedies to this problem. You propose pasting them together at some point (let the per1-curve start from a certain point on the per0-curve). I like this idea. If I do it on the simulation data, we get:

enter image description here

In our specific case, I think this represents data way better, but I don't know of any published results that support this approach. Heuristically, one can use the argument I presented in another answer:

KM plot with time-varying coefficient

Logistic – How to Weight Observations in Binary Logistic Regression

Weighting observations differently is appropriate when some measurements are made with greater precision than others.

In this case the problem is much more likely with the functional form of your model: assuming a linear relationship between the logit of flight probability & ship distance. Use a model that's suggested by theory, or one that's flexible enough to incorporate curvature: polynomials & restricted cubic splines are the standard choices; mine would probably be splines for the situation you described

Best Answer

Related Solutions

R – Visualize Survival Analysis with Time-Dependent Covariates

Logistic – How to Weight Observations in Binary Logistic Regression

Related Question