SEs rms cph vs survey svycoxph

Here are the model outputs.

Nathan

Survey package

ca.ATE.design <- svydesign(ids = ~ id, weights = ~ get.weights(ca.ATE.ps,
stop.method = 'ks.mean'), data = ca.dt)
Independent Sampling design (with replacement)

svydesign(ids = ~id, weights = ~get.weights(ca.ATE.ps, stop.method =
"ks.mean"), 
    data = ca.dt)

ca.ATE.dexmg.svy

Call:
svycoxph(formula = Surv(daysfromsurgerytodeath, as.logical(deceased)) ~
    dexamethasonemg + paincontrol + histgrade + adjuvant + stage +
        anesthetictransfusionunits, design = ca.ATE.design)


                              coef exp(coef) se(coef)      z       p
dexamethasonemg            -0.0863     0.917   0.0339 -2.550 1.1e-02
paincontrolNot Epidural     0.6027     1.827   0.2370  2.543 1.1e-02
histgradeg2                 0.9340     2.545   0.4307  2.168 3.0e-02
histgradeg3                 1.2749     3.578   0.4453  2.863 4.2e-03
adjuvantyes                -0.5810     0.559   0.2529 -2.298 2.2e-02
stageib                    -0.4394     0.644   0.6056 -0.726 4.7e-01
stageiia                    1.6565     5.241   0.5193  3.190 1.4e-03
stageiib                    1.6928     5.435   0.4902  3.453 5.5e-04
stageiii                    1.8211     6.179   0.5130  3.550 3.9e-04
stageiv                     2.3251    10.227   0.6940  3.350 8.1e-04
anesthetictransfusionunits  0.1963     1.217   0.0400  4.908 9.2e-07

Likelihood ratio test=  on 11 df, p=  n= 144, number of events= 102

rms package

ca.ATE.dexmg.rms2

Cox Proportional Hazards Model

cph(formula = Surv(daysfromsurgerytodeath, as.logical(deceased)) ~
    dexamethasonemg + paincontrol + histgrade + adjuvant + stage +
        anesthetictransfusionunits + cluster(id), data = ca.dt,
    weights = get.weights(ca.ATE.ps, stop.method = "ks.mean"),
    robust = T, x = T, y = T, se.fit = T, surv = T, time.inc = 30)

                    Model Tests       Discrimination
                                         Indexes
Obs       144    LR chi2    117.80    R2       0.559
Events    102    d.f.           11    Dxy     -0.459
Center 2.4016    Pr(> chi2) 0.0000    g        1.083
                 Score chi2 122.57    gr       2.953
                 Pr(> chi2) 0.0000

                           Coef    S.E.   Wald Z Pr(>|Z|)
dexamethasonemg            -0.0863 0.0192 -4.49  <0.0001
paincontrol=Not Epidural    0.6027 0.1203  5.01  <0.0001
histgrade=g2                0.9340 0.2209  4.23  <0.0001
histgrade=g3                1.2749 0.2612  4.88  <0.0001
adjuvant=yes               -0.5810 0.1741 -3.34  0.0008
stage=ib                   -0.4394 0.1899 -2.31  0.0207
stage=iia                   1.6565 0.2097  7.90  <0.0001
stage=iib                   1.6928 0.1979  8.55  <0.0001
stage=iii                   1.8211 0.2411  7.55  <0.0001
stage=iv                    2.3251 0.1886 12.33  <0.0001
anesthetictransfusionunits  0.1964 0.0214  9.17  <0.0001





From:  Thomas Lumley <tlumley at uw.edu>
Date:  Tuesday, February 25, 2014 at 3:09 PM
To:  "Nathan Leon Pace, MD, MStat" <n.l.pace at utah.edu>
Cc:  r help list <r-help at r-project.org>
Subject:  Re: [R] SEs rms cph vs survey svycoxph


On Tue, Feb 25, 2014 at 2:51 PM, Nathan Pace

I?ve used twang to get ATE propensity scores.

I?ve done multivariable, case weighted Cox PH models in survey using
svycoxph and in rms using cph with id(cluster) set to get robust estimates.

The model language is identical.

The point estimates are identical, but the CIs are considerably wider with
svycoxph estimates.

There is a note in the svycoxph help page stating the SEs should agree
closely unless the model fits poorly.




The actual note on the svycoxph help page says
"The standard errors agree closely with survfit.coxph for independent
sampling when the model fits well, but are larger when the model fits
poorly. "
That is, the note is for the survival curve rather than the coefficients.

It's still surprising that there's a big difference, but I think we need
more information.  

   -thomas

Thomas Lumley
Professor of Biostatistics
University of Auckland