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.

*

*In the absence of covariates, are models f2a and f2b really different from f1?

*Why does fitting model f2a result in an error while f2b does not?

*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
 A: 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.
A: 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.
