Robust Variance Estimate with Cox model in Complex Survey Setting

Question

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!

Answer 1

The answers are correct (though they may not be what you want).

When the first stage of sampling is treated as independent (ie, no fpc= argument, what survey statisticians misleadingly call "with replacement"), additional stages have no effect on the variance estimate.

Answer 2

In your last block of code you have the correct specification of a design, but you use the non-existent function design rather than svydesign and more importantly you don't pass the design object to your function so it doesn't do anything.

You want

design3 <- svydesign(id=~id_cluster,strat=~id_stratum,data=dataset) model3 <- svycoxph(Surv(time,censor)~exposure+confounder, design=design3)