How to model competing risks with frailties in R

I am trying to figure out if (a) a competing risks model with frailty-terms is the right method for my data and (b) how to estimate such models in R.

My data consists of N ~ 150 individuals with a total of ~ 50.000 (social interaction) events. For each social interaction we know the participant ID (id), how long it took since the last interaction event (time), and what type of social interaction it was (type), e.g., interaction with a friend, interaction with family member. For now, I assume that (a) these interaction types are mutually exclusive and (b) individuals are alone between social interactions (i.e., in the gap time). The aim of my analyses is to model the transition rates (from the alone state) to each type of social interaction. Is a competing risk model with independent frailty terms the correct model for this kind of data?

I have "simulated" such data like this:

set.seed(123)
simdat2 <- frailtySurv::genfrail(N = 30, K = 20,
beta = c(log(2)),
frailty = "gamma", theta = 2,
censor.rate = 0,
lambda_0=function(t, tau=4.6, C=0.01) (tau*(C*t)^tau)/t)

colnames(simdat2) <- c("id","rep","time","event","Covariate1")
simdat2$type <- factor(sample(c("EventType1","EventType2"), nrow(simdat2), replace = T)) # create fake covariate 3 that is associated with eventtype and time simdat2$$Covariate2 <- ifelse(simdat2$$type == "EventType1" & simdat2$time > 200,1,0) +
rnorm(nrow(simdat2),mean = 0.1, sd = 0.1)


So far, this is approach to model these data as competing risks to transition into either EventType1 or EventType2 accounting for the repeated-measures design (i.e., observations nested within individuals) with independent frailty terms:

#For event type 1
coxph(Surv(time, type == "EventType1") ~ Covariate2 + frailty(id), data = simdat2)
#For event type 2
coxph(Surv(time, type == "EventType2") ~ Covariate2 + frailty(id), data = simdat2)



According to Putter et al., 2007, p. 2404, the analysis above should be fine, but I am not sure, as they do not include a frailty term.

And: According to the authors of the survival R-package and their vignette on competing risk models, the following method should also work, but for me it does not. Why?

Error in coxph(Surv(time, type) ~ Covariate2 + frailty(id), id = id, data = simdat2) :
data set has overlapping intervals for one or more subjects


My questions are:

1. Is this the correct statistical approach to model these data?
2. Is this how you would estimate a competing risk model with independent frailty terms?
3. How would you interpret the two coefficients X_type1 = 0.35, X_type2 = -1.75? And why is type 2 larger and significant although I used type 1 to generate the Covariate2?
4. How would I go about in R to allow for the frailty terms to correlate?

I'd be happy for any hints or suggestions on these issues!

• Your simulation might not be modeling your situation well. You don't have competing risks in the usual sense. As the R survival vignette puts it: "The case of multiple event types, but only one event per subject is commonly known as competing risks." You have a multi-state repeated events situation, which is best handled by the counting-process Surv(time1, time2,eventType) data form. And are you modeling the transition back to the "alone" state from the interactions? If not, a point-process model might be better. – EdM Feb 15 at 16:04
• Thanks for your response! I am not modeling the transition back to the "alone" state, because there is no data on the duration of the interaction. If I do this by "faking" data where the time gap between any "interaction" and the "alone" stat is tiny (e.g., 1 sec), the estimates for the transition to the interaction states are similar but not identical between the counting-process and the approach I explained above. Do you know why these estimates are not identical? – teeglaze Feb 17 at 15:13

I would recommend modeling with a robust (cluster) variance estimator if you are having problems with a frailty model. (I have little experience with frailty modeling in survival analysis, accounting in no small part for my recommendation.) That will still provide control for intra-individual correlations. That's the way that the main R survival vignette deals with multi-state/multi-event modeling in Section 3.4. Nothing prevents you from generating your sample data based on frailty and then analyzing the data that way.
The data need to be set up in the counting process (start, stop, event) format, with your multi-level event factor. The id setting in the coxph() call, pointing to the column with the labels for the individuals, accounts for intra-individual correlations and generates robust error estimates that take those correlations into account. I think you can specify that the initial state for each time period is always "alone" by including that value as a column in the data, with that column indicated by the istate setting in the coxph() call. (I've never tried that myself, however.) If I'm correct, then you don't have to model the return to the "alone" state from each type of interaction; all of your modeling would be from the "alone" state to one of the interaction states.
I think that some of your problems might have arisen from using the simple (time, event) format for your data and analysis. That's OK for the usual "competing risks" scenario, in which only 1 event out of multiple possibilities occurs per individual. I'm not sure how well that works in the situation that you are modeling.