# Cox models in R: cluster, strata, robust, and frailty ... modeling grouped/clustered survival data

Consider a situation where we have individual patient survival data from a series of clinical trials of patients treated in a similar way. We might have a dataset that looks like this, with n total subjects from k total trials, experiencing events (event = 1) or censored (event = 0) at time t.

subject trial time event
1 1 3.6 1
2 1 5.4 1
3 1 4.6 0
4 2 5.1 0
... ... ... ...
n k - -

Note that there are no covariates (e.g. treatment, age, gender). We are interested in the survival of a group of people, whose information is collected from various trials, like might occur in a meta-analysis for example. Following Glidden DV et al. [Stat Med 2004] terminology, we could have a number of possible models for this data. Using R packages survival and coxme, we could fit them:

1. Marginal model, normal variance

f1 <- coxph(Surv(time, event) ~ 1)


2. Marginal model, robust variance

f2a <- coxph(Surv(time, event) ~ cluster(trial))

f2b <- coxph(Surv(time, event) ~ 1, robust = T)


3. Conditional model, fixed effects

f3 <- coxph(Surv(time, event) ~ trial)


4. Conditional model, stratified

f4 <- coxph(Surv(time, event) ~ strata(trial))


5. Conditional model, random effects (frailty)

f5a <- coxph(Surv(time, event) ~ frailty(trial))

f5b <- coxme(Surv(time, event) ~ (1 | trial))


I have a few questions about these models.

1. In the absence of covariates, are models f2a and f2b really different from f1?
2. Why does fitting model f2a result in an error while f2b does not?
3. How do I examine the baseline hazards of the various models (e.g. plot them) to see how these models differ in their treatment of the hazards?

I appreciate any insight that can be provided

• As suggested by John, muhaz has kernel density estimates of the baseline hazards which are very useful. I have also discovered the package bshazard which has smoothed baseline hazards and excellent plot functionality. Feb 2 at 14:01

In the absence of covariates, f2b and f1 are effectively equivalent as single models because robust=TRUE only affects the standard errors for the regression coefficients and there aren't any. They aren't exactly interchangeable if you consider them as base models that some other model is nested in (because that model will have regression coefficients).

f2a seems to have a bug that might not be worth fixing, because (a) the use of cluster in that way is discouraged and (b) there isn't any real point to asking for robust standard errors for coefficients when you aren't estimating any coefficients.

The predict.coxph function will give you predicted hazards (it says how to convert these into survival) for all the coxph objects.

• Thank you for your response and for confirming that f2a and f2b should give a similar result but don't in this case. Feb 2 at 13:57

I think f2b should be the same as f1. Do you get the same result?

It seems like f2a should be different, but can you explain what error you get with f2a? Can you also try this instead

f2a <- coxph(Surv(time, event) ~ 1,cluster=trial)
`

I would use the package muhaz to estimate the hazard functions within each trial and graph them. I don't know of any way to estimate the hazard function from these Cox regression models.

You can plot the estimated survival functions by using
predict(f3, newdata=...,type="survival")
newdata should have the time points where you want to estimate the survival.

• Thank you for your comments. Indeed f2b and f1 do appear identical, which is why I was curious about the application of robust standard errors when there are no covariates in the model. Feb 2 at 13:56