Assessing Proportional Hazards Assumption of a Cox Model with Caseweights

Question

I have data that I have used the R package twang with multinomial propensity scoring, to assign caseweights based on an ATT model. Once these caseweights have been assigned, they are typically used with a survey procedure. In R, there is the survey package which contains svycoxph, but this can also be done using the survival package with + cluster(id) in the model, to compute robust standard errors (gives same result for my caseweights).

I am have satisfactorily fit the non-propensity scored data as such, after reviewing cox.zph test and the Schoenfeld residual plots:

> coxph1 <- coxph(Surv(time,status) ~ x1 + x2 + age + tt(age), 
data=data, tt=function(x,t,...) x * log(t))

> cox.zph(coxph1)
            rho    chisq     p
x1      -0.0306     1.54 0.215
x2      -0.0367     2.21 0.137
age     -0.5752 37284.88 0.000
tt(age)  0.5789 48667.57 0.000
GLOBAL       NA 48676.43 0.000

x1 and x2 are categorical (1=yes, 0=no) predictors.

I have removed some variables for simplicity of explanation, disregard the age variable.

When I fit a model using the weights and robust standard error measurement, I am concerned that the use of the weights are obscuring the cox.zph and Schoenfeld plots:

> coxph2 <- coxph(Surv(time,status)~ x1 + x2 + age + tt(age) + cluster(id), weights=w, data=data, tt=function(x,t,...) x * log(t))
> cox.zph(coxph2)
           rho   chisq                   p
x1       0.182    60.1 0.00000000000000888
x2      -0.128    42.5 0.00000000006934375
age     -0.581 54338.6 0.00000000000000000
tt(age)  0.586 70624.8 0.00000000000000000
GLOBAL      NA 70734.7 0.00000000000000000

The Schoenfeld plots look accordingly bad.

This gets even worse as more predictors are added (in terms of how it looks in comparison to the unweighted version).

Side question: x1 has a maximum follow up time of 60, whereas x2 and the reference category have follow up to 100 months. Is there any chance this is having an impact on the caseweights/their relation to the residual?

This issue persists if using the survey package and svycox.

I have two questions: 1) Is there in fact an issue with using caseweights and or robust SE estimation in how the Schoenfeld residual is constructed, and should these be avoided? 2) If 1 is true, is there another preferred method for testing the PH assumption? So far, I have been adding a time transformation tt() for each and seeing if that is significant, but that is far from an ideal approach. The only other alternative I see is to treat those that violate the assumption in the unweighted model the same way in the weighted model.

Answer 1

Update:

I think I have found a solution to the problem. I edited the original cox.zph function to included weighted Schoenfeld residuals:

Here was the original:

cox.zph <- function(fit, transform='km', global=TRUE) {
call <- match.call()
if (!inherits(fit, 'coxph')) stop ("Argument must be the result of coxph")
if (inherits(fit, 'coxph.null'))
stop("The are no score residuals for a Null model")

sresid <- resid(fit, 'schoenfeld')
varnames <- names(fit$coefficients)
nvar <- length(varnames)
ndead<- length(sresid)/nvar
if (nvar==1) times <- as.numeric(names(sresid))
else         times <- as.numeric(dimnames(sresid)[[1]])

And the change:

cox.zph <- function(fit, transform='km', global=TRUE) {
call <- match.call()
if (!inherits(fit, 'coxph')) stop ("Argument must be the result of coxph")
if (inherits(fit, 'coxph.null'))
stop("The are no score residuals for a Null model")

sresid <- resid(fit, 'schoenfeld', **weighted=TRUE**)
varnames <- names(fit$coefficients)
nvar <- length(varnames)
ndead<- length(sresid)/nvar
if (nvar==1) times <- as.numeric(names(sresid))
else         times <- as.numeric(dimnames(sresid)[[1]])

Does anyone see potential problems with solution? It appears to have made both the zph test and the scaled schoenfeld residual plots much more reasonable.