I have a single-stage cluster sample, and I am trying to estimate the hazard ratio of a given exposure after controlling for a confounder. My dataset contains 10 strata, with 20 clusters within each stratum. Furthermore, my dataset contains repeated observations, so I want to obtain a robust variance estimate that accounts for the clusters due to these repeated observations. The strata are identified by $id\_stratum$, the clusters are identified by $id\_cluster$, and the repeated observations are identified by $id\_unique$.
If we ignore the complex design features, we can easily obtain robust variance estimates that account for the repeated observations by specifying $cluster(id\_unique)$, i.e.
model <- coxph(Surv(time,censor)~exposure+confounder+cluster(id_unique),data=dataset)
coef exp(coef) se(coef) robust se z Pr(>|z|)
exposure 0.66509 1.94466 0.07115 0.10234 6.499 8.08e-11 ***
confounder 0.16797 1.18290 0.06465 0.01614 10.408 < 2e-16 ***
However, I don't know how to obtain robust variance estimates while accounting for the complex design features. I have tried the following:
1) Adding $cluster(id\_unique)$ to the svycoxph function after specifying the design using svydesign argument.
design1 <- design(id=~id_cluster,strat=~id_stratum,data=dataset)
model1 <- coxph(Surv(time,censor)~exposure+confounder+cluster(id_unique),data=dataset)
coef exp(coef) se(coef) robust se z p
exposure 0.6651 1.9447 0.0712 0.0992 6.71 2.0e-11
confounder 0.1680 1.1829 0.0646 0.0173 9.68 < 2e-16
2) Specifying another level of clustering for $id\_unique$ in the svydesign argument.
design2 <- design(id=~id_cluster+id_cluster_unique,strat=~id_stratum,data=dataset)
model2 <- coxph(Surv(time,censor)~exposure+confounder,data=dataset)
coef exp(coef) se(coef) robust se z p
exposure 0.6651 1.9447 0.0712 0.0992 6.71 2.0e-11
confounder 0.1680 1.1829 0.0646 0.0173 9.68 < 2e-16
While the above results coincide, they are also equal to the variance estimates that I get when I don't consider $id\_unique$ at all, which suggests that I am not obtaining the right answers.
design3 <- design(id=~id_cluster,strat=~id_stratum,data=dataset)
model3 <- coxph(Surv(time,censor)~exposure+confounder,data=dataset)
coef exp(coef) se(coef) robust se z p
exposure 0.6651 1.9447 0.0712 0.0992 6.71 2.0e-11
confounder 0.1680 1.1829 0.0646 0.0173 9.68 < 2e-16
I have read through the documentation for the survey package, and I have also read Thomas Lumley's book, but I can't find anything that addresses the issue of robust variance estimation in a complex survey setting. Any suggestions would be greatly appreciated. Thanks!