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
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.